MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_trace_field.t
Go to the documentation of this file.
3  use mod_physics
4  implicit none
5 
6 contains
7 
8  subroutine trace_field_multi(xfm,wPm,wLm,dL,numL,numP,nwP,nwL,forwardm,ftype,tcondi)
9  ! trace multiple field lines
10  ! xfm: locations of points at the field lines. User should provide xfm(1:numL,1,1:ndim)
11  ! as seed points, then field line wills be traced from the seed points. xfm(1:numL,1,:)
12  ! are user given and other points are given by the subroutine as feedback.
13  ! numL: number of field lines user wants to trace. user given
14  ! numP: maximum number of points at the field line. User defined. Note that not every
15  ! point of numP is valid, as the tracing will stop when the field line leave the
16  ! simulation box. The number of valid point is given in wLm(1)
17  ! wPm: point variables, which have values at each point of xfm. The way to get wP is
18  ! user defined, with help of IO given by subroutine set_field_w in mod_usr_methods.
19  ! User can calculate the density/temperature at the point and then store the values
20  ! into wPm, or do something else.
21  ! nwP: number of point variables. user given
22  ! wLm: line variables, variables for the lines rather than for each point of the lines.
23  ! For example, wLm(1) stores the number of valid points at the field lines. The way to
24  ! get wLm is also user defined with set_field_w, the same as wPm. User can calculate
25  ! the maximum density/temperature at the field lines and stores it in wLm
26  ! nwL: number of line variables. user given
27  ! dL: length step for field line tracing. user given
28  ! forwardm: true--trace field forward; false--trace field line backward. user given
29  ! ftype: type of field user wants to trace. User can trace velocity field by setting
30  ! ftype='Vfield' or trace magnetic field by setting ftype='Bfield'. It is possible
31  ! to trace other fields, e.g. electric field, where user can define the field with
32  ! IO given by subroutine set_field in mod_usr_methods. user given
33  ! tcondi: user given
35 
36  integer, intent(in) :: numL,numP,nwP,nwL
37  double precision, intent(inout) :: xfm(numL,numP,ndim),wPm(numL,numP,nwP),wLm(numL,1+nwL)
38  double precision, intent(in) :: dL
39  logical, intent(in) :: forwardm(numL)
40  character(len=std_len), intent(in) :: ftype,tcondi
41 
42  integer :: indomain,ipe_now,igrid_now,igrid,j,iL
43  integer :: ipoint_in,ipoint_out,numSend,nRT,nRTmax
44  double precision :: x3d(3),statusF(4+ndim),statusL(numL,4+ndim+nwL),statusS(numL,4+ndim+nwL)
45  logical :: continueL(numL),myL(numL)
46  logical :: stopT,forward
47  integer :: ipointm(numL),igridm(numL)
48  double precision :: xf(numP,ndim),wP(numP,nwP),wL(1+nwL)
49  double precision, allocatable :: data_send(:,:,:),data_recv(:,:,:)
50 
51  if (tcondi/='TRAC') then
52  wpm=zero
53  else
54  wpm=-1
55  endif
56  wlm=zero
57  xfm(1:numl,2:nump,:)=zero
58  stopt=.true.
59  myl=.false.
60  xf=zero
61  wp=zero
62 
63  ! find the pe and igrid for the first point
64  do il=1,numl
65  indomain=0
66  wlm(il,1)=0
67  {if (xfm(il,1,^db)>=xprobmin^db .and. xfm(il,1,^db)<xprobmax^db) indomain=indomain+1\}
68  if (indomain==ndim) then
69  if (tcondi/='TRAC') wlm(il,1)=1
70  continuel(il)=.true.
71  ! find pe and igrid
72  x3d=0.d0
73  do j=1,ndim
74  x3d(j)=xfm(il,1,j)
75  enddo
76  call find_particle_ipe(x3d,igrid_now,ipe_now)
77  stopt=.false.
78  ipointm(il)=1
79  igridm(il)=igrid_now
80  if (mype==ipe_now) then
81  myl(il)=.true.
82  else
83  xfm(il,1,:)=zero
84  endif
85  else
86  continuel(il)=.false.
87  wlm(il,1)=zero
88  endif
89  enddo
90 
91  do while(stopt .eqv. .false.)
92  ! tracing multiple field lines inside pe
93  statuss=zero
94  do il=1,numl
95  if (myl(il) .and. continuel(il)) then
96  igrid=igridm(il)
97  ipoint_in=ipointm(il)
98  xf(ipoint_in,:)=xfm(il,ipoint_in,:)
99  wl(:)=wlm(il,:)
100  forward=forwardm(il)
101  statusf=zero
102  call find_points_in_pe(igrid,ipoint_in,xf,wp,wl,dl,nump,nwp,nwl,forward,ftype,tcondi,statusf)
103  ipoint_out=int(statusf(1))
104  xfm(il,ipoint_in:ipoint_out-1,:)=xf(ipoint_in:ipoint_out-1,:)
105  wpm(il,ipoint_in:ipoint_out-1,:)=wp(ipoint_in:ipoint_out-1,:)
106  ! status for each field line
107  ! 1: index of next point
108  ! 2: ipe of next point
109  ! 3: igrid of next point
110  ! 4: 1-> continue tracing; 0-> stop tracing
111  ! 5:4+ndim: coordinate of next point
112  ! 4+ndim+1:4+ndim+nwL: wL(2:1+nwL)
113  ! for TRAC nwL=2 -> wL(2): current Tcoff; wL(3): Tmax
114  ! for TRAC nwP=2 -> wP(:,1):ipe; wP(:,2):igrid
115  statuss(il,1:4+ndim)=statusf(1:4+ndim)
116  statuss(il,4+ndim+1:4+ndim+nwl)=wl(2:1+nwl)
117  if (tcondi=='TRAC') wlm(il,1)=ipoint_out-1
118  endif
119  enddo
120 
121  ! comunicating tracing results
122  numsend=numl*(4+ndim+nwl)
123  call mpi_allreduce(statuss,statusl,numsend,mpi_double_precision,&
124  mpi_sum,icomm,ierrmpi)
125 
126  ! for next step
127  stopt=.true.
128  myl=.false.
129  do il=1,numl
130  if (continuel(il)) then
131  ipointm(il)=int(statusl(il,1))
132  if (mype==int(statusl(il,2))) myl(il)=.true.
133  igridm(il)=int(statusl(il,3))
134  if (int(statusl(il,4))==0) then
135  stopt=.false.
136  else
137  continuel(il)=.false.
138  endif
139  if (myl(il)) xfm(il,ipointm(il),1:ndim)=statusl(il,4+1:4+ndim)
140  if (tcondi/='TRAC') wlm(il,1)=ipointm(il)-1
141  wlm(il,2:1+nwl)=statusl(il,4+ndim+1:4+ndim+nwl)
142  endif
143  enddo
144  enddo
145 
146  ! communication after tracing
147  if (tcondi/='TRAC') then
148  nrtmax=0
149  do il=1,numl
150  if (nrtmax<int(wlm(il,1))) nrtmax=int(wlm(il,1))
151  enddo
152  numsend=numl*nrtmax*(ndim+nwp)
153 
154  allocate(data_send(numl,nrtmax,ndim+nwp),data_recv(numl,nrtmax,ndim+nwp))
155  data_send(:,:,:)=zero
156  do il=1,numl
157  nrt=int(wlm(il,1))
158  data_send(il,1:nrt,1:ndim)=xfm(il,1:nrt,1:ndim)
159  if (nwp>0) data_send(il,1:nrt,1+ndim:ndim+nwp)=wpm(il,1:nrt,1:nwp)
160  enddo
161  call mpi_allreduce(data_send,data_recv,numsend,mpi_double_precision,&
162  mpi_sum,icomm,ierrmpi)
163  do il=1,numl
164  nrt=int(wlm(il,1))
165  xfm(il,1:nrt,1:ndim)=data_recv(il,1:nrt,1:ndim)
166  if (nwp>0) wpm(il,1:nrt,1:nwp)=data_recv(il,1:nrt,1+ndim:ndim+nwp)
167  enddo
168  deallocate(data_send,data_recv)
169  endif
170 
171  end subroutine trace_field_multi
172 
173  subroutine trace_field_single(xf,wP,wL,dL,numP,nwP,nwL,forward,ftype,tcondi)
174  ! trace a field line
175  ! xf: locations of points at the field line. User should provide xf(1,1:ndim)
176  ! as seed point, then field line will be traced from the seed point. xf(1,:) is
177  ! user given and xf(2:wL(1),:) are given by the subroutine as feedback.
178  ! numP: maximum number of points at the field line. User defined. Note that not every
179  ! point of numP is valid, as the tracing will stop when the field line leave the
180  ! simulation box. The number of valid point is given in wL(1)
181  ! wP: point variables, which have values at each point of xf. The way to get wP is
182  ! user defined, with help of IO given by subroutine set_field_w in mod_usr_methods.
183  ! User can calculate the density/temperature at the point and then store the values
184  ! into wP, or do something else.
185  ! nwP: number of point variables. user given
186  ! wL: line variables, variables for the line rather than for each point of the line.
187  ! For example, wL(1) stores the number of valid points at the field line. The way to
188  ! get wL is also user defined with set_field_w, the same as wP. User can calculate
189  ! the maximum density/temperature at the field line and stores it in wL
190  ! nwL: number of line variables. user given
191  ! dL: length step for field line tracing. user given
192  ! forward: true--trace field forward; false--trace field line backward. user given
193  ! ftype: type of field user wants to trace. User can trace velocity field by setting
194  ! ftype='Vfield' or trace magnetic field by setting ftype='Bfield'. It is possible
195  ! to trace other fields, e.g. electric field, where user can define the field with
196  ! IO given by subroutine set_field in mod_usr_methods. user given
197  ! tcondi: user given
198  use mod_usr_methods
200 
201  integer, intent(in) :: numP,nwP,nwL
202  double precision, intent(inout) :: xf(numP,ndim),wP(numP,nwP),wL(1+nwL)
203  double precision, intent(in) :: dL
204  logical, intent(in) :: forward
205  character(len=std_len), intent(in) :: ftype,tcondi
206 
207  integer :: indomain,ipoint_in,ipe_now,igrid_now,igrid,j
208  integer :: ipoint_out,ipe_next,igrid_next,numRT
209  double precision :: x3d(3),statusF(4+ndim),status_bcast(4+ndim+nwL)
210  logical :: stopT
211  double precision, allocatable :: data_send(:,:),data_recv(:,:)
212 
213  wp=zero
214  wl=zero
215  xf(2:nump,:)=zero
216 
217  ! check whether or the first point is inside simulation box. if yes, find
218  ! the pe and igrid for the point
219  indomain=0
220  wl(1)=0
221  {if (xf(1,^db)>=xprobmin^db .and. xf(1,^db)<xprobmax^db) indomain=indomain+1\}
222  if (indomain==ndim) then
223  wl(1)=1
224 
225  ! find pe and igrid
226  x3d=0.d0
227  do j=1,ndim
228  x3d(j)=xf(1,j)
229  enddo
230  call find_particle_ipe(x3d,igrid_now,ipe_now)
231  stopt=.false.
232  ipoint_in=1
233  if (mype/=ipe_now) xf(1,:)=zero
234  else
235  if (mype==0) then
236  call mpistop('Field tracing error: given point is not in simulation box!')
237  endif
238  endif
239 
240 
241  ! other points in field line
242  do while(stopt .eqv. .false.)
243 
244  if (mype==ipe_now) then
245  igrid=igrid_now
246  ! looking for points in one pe
247  call find_points_in_pe(igrid,ipoint_in,xf,wp,wl,dl,nump,nwp,nwl,forward,ftype,tcondi,statusf)
248  status_bcast(1:4+ndim)=statusf(1:4+ndim)
249  status_bcast(4+ndim+1:4+ndim+nwl)=wl(2:1+nwl)
250  endif
251  ! comunication
252  call mpi_bcast(status_bcast,4+ndim+nwl,mpi_double_precision,ipe_now,icomm,ierrmpi)
253  statusf(1:4+ndim)=status_bcast(1:4+ndim)
254  wl(2:1+nwl)=status_bcast(4+ndim+1:4+ndim+nwl)
255 
256  ! prepare for next step
257  ipoint_out=int(statusf(1))
258  ipe_next=int(statusf(2))
259  igrid_next=int(statusf(3))
260  if (int(statusf(4))==1) then
261  stopt=.true.
262  wl(1)=ipoint_out-1
263  endif
264  if (mype==ipe_next) then
265  do j=1,ndim
266  xf(ipoint_out,j)=statusf(4+j)
267  enddo
268  else
269  xf(ipoint_out,:)=zero
270  endif
271 
272  ! pe and grid of next point
273  ipe_now=ipe_next
274  igrid_now=igrid_next
275  ipoint_in=ipoint_out
276  enddo
277 
278  if (tcondi/='TRAC') then
279  numrt=int(wl(1))
280  allocate(data_send(numrt,ndim+nwp),data_recv(numrt,ndim+nwp))
281  data_send(:,:)=zero
282  data_recv(:,:)=zero
283  data_send(1:numrt,1:ndim)=xf(1:numrt,1:ndim)
284  if (nwp>0) data_send(1:numrt,1+ndim:ndim+nwp)=wp(1:numrt,1:nwp)
285  call mpi_allreduce(data_send,data_recv,numrt*(ndim+nwp),mpi_double_precision,&
286  mpi_sum,icomm,ierrmpi)
287  xf(1:numrt,1:ndim)=data_recv(1:numrt,1:ndim)
288  if (nwp>0) wp(1:numrt,1:nwp)=data_recv(1:numrt,1+ndim:ndim+nwp)
289  deallocate(data_send,data_recv)
290  endif
291 
292  end subroutine trace_field_single
293 
294  subroutine find_points_in_pe(igrid,ipoint_in,xf,wP,wL,dL,numP,nwP,nwL,forward,ftype,tcondi,statusF)
295 
296  integer, intent(inout) :: igrid
297  integer, intent(in) :: ipoint_in,numP,nwP,nwL
298  double precision, intent(inout) :: xf(numP,ndim),wP(numP,nwP),wL(1+nwL)
299  double precision, intent(in) :: dL
300  logical, intent(in) :: forward
301  character(len=std_len), intent(in) :: ftype,tcondi
302  double precision, intent(inout) :: statusF(4+ndim)
303 
304  integer :: ipe_next,igrid_next,ip_in,ip_out,j,indomain
305  logical :: newpe,stopT
306  double precision :: xfout(ndim)
307 
308  ip_in=ipoint_in
309  newpe=.false.
310 
311  do while(newpe .eqv. .false.)
312  ! looking for points in given grid
313  call find_points_interp(igrid,ip_in,ip_out,xf,wp,wl,nump,nwp,nwl,dl,forward,ftype,tcondi)
314  ip_in=ip_out
315 
316  ! when next point is out of given grid, find next grid
317  indomain=0
318  {if (xf(ip_out,^db)>=xprobmin^db .and. xf(ip_out,^db)<=xprobmax^db) indomain=indomain+1\}
319  if (ip_out<nump .and. indomain==ndim) then
320  if (tcondi/='TRAC') then
321  stopt=.false.
322  xfout=xf(ip_out,:)
323  call find_next_grid(igrid,igrid_next,ipe_next,xfout,newpe,stopt)
324  else
325  if (xf(ip_out,ndim)>phys_trac_mask) then
326  newpe=.true.
327  stopt=.true.
328  else
329  stopt=.false.
330  xfout=xf(ip_out,:)
331  call find_next_grid(igrid,igrid_next,ipe_next,xfout,newpe,stopt)
332  endif
333  endif
334  else
335  newpe=.true.
336  stopt=.true.
337  endif
338 
339  if (newpe) then
340  statusf(1)=ip_out
341  statusf(2)=ipe_next
342  statusf(3)=igrid_next
343  statusf(4)=0
344  if (stopt) statusf(4)=1
345  do j=1,ndim
346  statusf(4+j)=xf(ip_out,j)
347  enddo
348  endif
349 
350  if (newpe .eqv. .false.) igrid=igrid_next
351  enddo
352 
353  end subroutine find_points_in_pe
354 
355  subroutine find_next_grid(igrid,igrid_next,ipe_next,xf1,newpe,stopT)
356  ! check the grid and pe of next point
357  use mod_usr_methods
359  use mod_forest
360 
361  integer, intent(inout) :: igrid,igrid_next,ipe_next
362  double precision, intent(in) :: xf1(ndim)
363  logical, intent(inout) :: newpe,stopT
364 
365  double precision :: dxb^D,xb^L,xbmid^D
366  integer :: idn^D,my_neighbor_type,inblock
367  integer :: ic^D,inc^D,ipe_neighbor,igrid_neighbor
368  double precision :: xbn^L
369 
370  ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
371  ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
372  inblock=0
373 
374  ! direction of next grid
375  idn^d=0\
376  {if (xf1(^d)<=xbmin^d) idn^d=-1\}
377  {if (xf1(^d)>=xbmax^d) idn^d=1\}
378  my_neighbor_type=neighbor_type(idn^d,igrid)
379  igrid_neighbor=neighbor(1,idn^d,igrid)
380  ipe_neighbor=neighbor(2,idn^d,igrid)
381 
382  ! ipe and igrid of next grid
383  select case(my_neighbor_type)
384  case (neighbor_boundary)
385  ! next point is not in simulation box
386  newpe=.true.
387  stopt=.true.
388 
389  case(neighbor_coarse)
390  ! neighbor grid has lower refinement level
391  igrid_next=igrid_neighbor
392  ipe_next=ipe_neighbor
393  if (mype==ipe_neighbor) then
394  newpe=.false.
395  else
396  newpe=.true.
397  endif
398 
399  case(neighbor_sibling)
400  ! neighbor grid has lower refinement level
401  igrid_next=igrid_neighbor
402  ipe_next=ipe_neighbor
403  if (mype==ipe_neighbor) then
404  newpe=.false.
405  else
406  newpe=.true.
407  endif
408 
409  case(neighbor_fine)
410  ! neighbor grid has higher refinement level
411  {xbmid^d=(xbmin^d+xbmax^d)/2.d0\}
412  ^d&inc^d=1\
413  {if (xf1(^d)<=xbmin^d) inc^d=0\}
414  {if (xf1(^d)>xbmin^d .and. xf1(^d)<=xbmid^d) inc^d=1\}
415  {if (xf1(^d)>xbmid^d .and. xf1(^d)<xbmax^d) inc^d=2\}
416  {if (xf1(^d)>=xbmax^d) inc^d=3\}
417  ipe_next=neighbor_child(2,inc^d,igrid)
418  igrid_next=neighbor_child(1,inc^d,igrid)
419  if (mype==ipe_next) then
420  newpe=.false.
421  else
422  newpe=.true.
423  endif
424  end select
425 
426  end subroutine find_next_grid
427 
428  subroutine find_points_interp(igrid,ip_in,ip_out,xf,wP,wL,numP,nwP,nwL,dL,forward,ftype,tcondi)
430  use mod_usr_methods
431 
432  integer, intent(in) :: igrid,ip_in,numP,nwP,nwL
433  integer, intent(inout) :: ip_out
434  double precision, intent(inout) :: xf(numP,ndim),wP(numP,nwP),wL(1+nwL)
435  double precision, intent(in) :: dL
436  logical, intent(in) :: forward
437  character(len=std_len), intent(in) :: ftype,tcondi
438 
439  double precision :: dxb^D,xb^L
440  integer :: ip,inblock,ixI^L,ixO^L,j
441  double precision :: field(ixg^T,ndir)
442  double precision :: xs1(ndim),xs2(ndim),K1(ndim),K2(ndim)
443  double precision :: xfpre(ndim),xfnow(ndim),xfnext(ndim)
444  double precision :: Tpre,Tnow,Tnext,dTds,Lt,Lr,ds,T_bott,trac_delta
445 
446  ixi^l=ixg^ll;
447  ixo^l=ixm^ll;
448  ^d&dxb^d=rnode(rpdx^d_,igrid);
449  ^d&xbmin^d=rnode(rpxmin^d_,igrid);
450  ^d&xbmax^d=rnode(rpxmax^d_,igrid);
451 
452  if (tcondi/='TRAC') then
453  ds=dl
454  else
455  ds=dxb^nd
456  t_bott=2.d4/unit_temperature
457  trac_delta=0.25d0
458  endif
459 
460  ! main loop
461  mainloop: do ip=ip_in,nump-1
462 
463  ! integrate magnetic field with Runge-Kutta method
464  xs1(:)=xf(ip,:)
465  call get_k(xs1,igrid,k1,ixi^l,dxb^d,ftype)
466  if (forward) then
467  xs2(:)=xf(ip,:)+ds*k1(:)
468  else
469  xs2(:)=xf(ip,:)-ds*k1(:)
470  endif
471  call get_k(xs2,igrid,k2,ixi^l,dxb^d,ftype)
472  if (forward) then
473  xf(ip+1,:)=xf(ip,:)+ds*(0.5*k1(:)+0.5*k2(:))
474  else
475  xf(ip+1,:)=xf(ip,:)-ds*(0.5*k1(:)+0.5*k2(:))
476  endif
477  ip_out=ip+1
478 
479  ! get local values for variable via interpolation
480  if (tcondi/='TRAC') then
481  if (associated(usr_set_field_w)) then
482  call usr_set_field_w(igrid,ip,xf,wp,wl,nump,nwp,nwl,dl,forward,ftype,tcondi)
483  endif
484  else ! get TRAC Tcoff
485  wp(ip,1)=mype
486  wp(ip,2)=igrid
487  if (ip==ip_in) then
488  if (forward) then
489  xfpre(:)=xf(ip,:)-ds*k1(:)
490  else
491  xfpre(:)=xf(ip,:)+ds*k1(:)
492  endif
493  xfnow(:)=xf(ip,:)
494  xfnext(:)=xf(ip+1,:)
495  call get_t_loc_trac(igrid,xfpre,tpre,ixi^l,dxb^d)
496  call get_t_loc_trac(igrid,xfnow,tnow,ixi^l,dxb^d)
497  call get_t_loc_trac(igrid,xfnext,tnext,ixi^l,dxb^d)
498  else
499  xfpre=xf(ip-1,:)
500  xfnow(:)=xf(ip,:)
501  xfnext(:)=xf(ip+1,:)
502  tpre=tnow
503  tnow=tnext
504  call get_t_loc_trac(igrid,xfnext,tnext,ixi^l,dxb^d)
505  endif
506  dtds=abs(tnext-tpre)/(2*ds)
507  if (ip==1) then
508  lt=0.d0
509  wl(2)=t_bott ! current Tcofl
510  wl(3)=tnow ! current Tlmax
511  else
512  lt=0.d0
513  if (dtds>0.d0) then
514  lt=tnow/dtds
515  lr=ds
516  ! renew cutoff temperature
517  if(lr>trac_delta*lt) then
518  if (tnow>wl(2)) wl(2)=tnow
519  endif
520  endif
521  if (tnow>wl(3)) wl(3)=tnow
522  endif
523  endif
524 
525  ! exit loop if next point is not in this block
526  inblock=0
527  {if (xf(ip+1,^db)>=xbmin^db .and. xf(ip+1,^db)<xbmax^db) inblock=inblock+1\}
528  if (tcondi=='TRAC' .and. xf(ip+1,ndim)>phys_trac_mask) inblock=0
529  if (inblock/=ndim) exit mainloop
530 
531  enddo mainloop
532 
533  end subroutine find_points_interp
534 
535  subroutine get_k(xfn,igrid,K,ixI^L,dxb^D,ftype)
536  use mod_usr_methods
537 
538  integer :: ixI^L,igrid
539  double precision :: dxb^D
540  double precision :: xfn(ndim),K(ndim)
541  character(len=std_len) :: ftype
542 
543  integer :: ixb^D,ix^D,ixbl^D,ixO^L,j
544  double precision :: xd^D
545  double precision :: field(0:1^D&,ndir),Fx(ndim),factor(0:1^D&)
546  double precision :: vector(ixI^S,1:ndir)
547  double precision :: Ftotal
548 
549  ^d&ixbl^d=floor((xfn(^d)-ps(igrid)%x(iximin^dd,^d))/dxb^d)+iximin^d;
550  ^d&xd^d=(xfn(^d)-ps(igrid)%x(ixbl^dd,^d))/dxb^d;
551  ^d&ixomin^d=ixbl^d;
552  ^d&ixomax^d=ixbl^d+1;
553  vector=zero
554 
555  field=zero
556  if(ftype=='Bfield') then
557  if(b0field) then
558  if(allocated(iw_mag)) then
559  vector(ixo^s,1:ndir)=ps(igrid)%w(ixo^s,iw_mag(1:ndir))+ps(igrid)%B0(ixo^s,1:ndir,0)
560  else
561  vector(ixo^s,1:ndir)=ps(igrid)%B0(ixo^s,1:ndir,0)
562  endif
563  else
564  vector(ixo^s,1:ndir)=ps(igrid)%w(ixo^s,iw_mag(1:ndir))
565  endif
566  else if (ftype=='Vfield') then
567  call phys_get_v(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,vector)
568  endif
569  {do ix^d=0,1\}
570  factor(ix^d)={abs(1-ix^d-xd^d)*}
571  field(ix^d,1:ndir)=vector(ixbl^d+ix^d,1:ndir)
572  {enddo\}
573 
574  fx=0.d0
575  {do ix^db=0,1\}
576  do j=1,ndim
577  fx(j)=fx(j)+field(ix^d,j)*factor(ix^d)
578  enddo
579  {enddo\}
580 
581  if (ftype=='Bfield' .or. ftype=='Vfield') then
582  fx=0.d0
583  {do ix^db=0,1\}
584  do j=1,ndim
585  fx(j)=fx(j)+field(ix^d,j)*factor(ix^d)
586  enddo
587  {enddo\}
588  else if (associated(usr_set_field)) then
589  call usr_set_field(xfn,igrid,fx,ftype)
590  else
591  call mpistop('Field tracing error: wrong field type!')
592  endif
593 
594  ftotal=zero
595  do j=1,ndim
596  ftotal=ftotal+(fx(j))**2
597  enddo
598  ftotal=dsqrt(ftotal)
599 
600  if (ftotal==0.d0) then
601  k=0.d0
602  k(1)=1.d0
603  else
604  k(1:ndim)=fx(1:ndim)/ftotal
605  endif
606 
607  end subroutine get_k
608 
609  subroutine get_t_loc_trac(igrid,xloc,Tloc,ixI^L,dxb^D)
610  ! grid T has been calculated and stored in wextra(ixI^S,Tcoff_)
611  ! see mod_trac
612  integer, intent(in) :: igrid,ixI^L
613  double precision, intent(inout) :: xloc(ndim)
614  double precision, intent(inout) :: Tloc
615  double precision, intent(in) :: dxb^D
616 
617  integer :: ixb^D,ix^D,ixbl^D,j,ixO^L
618  double precision :: xd^D
619  double precision :: factor(0:1^D&),Tnear(0:1^D&)
620 
621  ^d&ixbl^d=floor((xloc(^d)-ps(igrid)%x(iximin^dd,^d))/dxb^d)+iximin^d;
622  ^d&xd^d=(xloc(^d)-ps(igrid)%x(ixbl^dd,^d))/dxb^d;
623  ^d&ixomin^d=ixbl^d;
624  ^d&ixomax^d=ixomin^d+1;
625 
626  {do ix^d=0,1\}
627  factor(ix^d)={abs(1-ix^d-xd^d)*}
628  tnear(ix^d)=ps(igrid)%wextra(ixbl^d+ix^d,iw_tcoff)
629  {enddo\}
630 
631  tloc=0.d0
632  ! interpolation
633  {do ix^db=0,1\}
634  tloc=tloc+tnear(ix^d)*factor(ix^d)
635  {enddo\}
636 
637  end subroutine get_t_loc_trac
638 
639 end module mod_trace_field
Module with basic grid data structures.
Definition: mod_forest.t:2
This module contains definitions of global parameters and variables and some generic functions/subrou...
double precision phys_trac_mask
integer icomm
The MPI communicator.
integer mype
The rank of the current MPI task.
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
Module with shared functionality for all the particle movers.
subroutine find_particle_ipe(x, igrid_particle, ipe_particle)
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
procedure(sub_get_v), pointer phys_get_v
Definition: mod_physics.t:69
subroutine trace_field_single(xf, wP, wL, dL, numP, nwP, nwL, forward, ftype, tcondi)
subroutine get_k(xfn, igrid, K, ixIL, dxbD, ftype)
subroutine find_next_grid(igrid, igrid_next, ipe_next, xf1, newpe, stopT)
subroutine get_t_loc_trac(igrid, xloc, Tloc, ixIL, dxbD)
subroutine trace_field_multi(xfm, wPm, wLm, dL, numL, numP, nwP, nwL, forwardm, ftype, tcondi)
subroutine find_points_in_pe(igrid, ipoint_in, xf, wP, wL, dL, numP, nwP, nwL, forward, ftype, tcondi, statusF)
subroutine find_points_interp(igrid, ip_in, ip_out, xf, wP, wL, numP, nwP, nwL, dL, forward, ftype, tcondi)
Module with all the methods that users can customize in AMRVAC.
procedure(set_field_w), pointer usr_set_field_w
procedure(set_field), pointer usr_set_field