MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_trac.t
Go to the documentation of this file.
1 module mod_trac
3  use mod_mhd
4  use mod_functions_bfield, only: mag
5 
6  implicit none
7  private
8 
9  ! common
10  integer :: numFL,numLP
11  double precision :: dL,Tmax,trac_delta,T_bott
12  double precision, allocatable :: xFi(:,:)
13  ! TRACB
14  integer :: numxT^D
15  double precision :: dxT^D
16  double precision :: xTmin(ndim),xTmax(ndim)
17  double precision, allocatable :: xT(:^D&,:)
18  ! TRACL mask
19  integer, allocatable :: trac_grid(:),ground_grid(:)
20  integer :: ngrid_trac,ngrid_ground
21  logical, allocatable :: trac_pe(:)
22 
24 
25 contains
26 
27  subroutine init_trac_line(mask)
28  logical, intent(in) :: mask
29  integer :: refine_factor,ix^D,ix(ndim),j,iFL,numL(ndim),finegrid
30  double precision :: lengthFL
31  double precision :: xprobmin(ndim),xprobmax(ndim),domain_nx(ndim)
32  integer :: memxFi
33 
34  trac_delta=0.25d0
35  !----------------- init magnetic field -------------------!
36  refine_factor=2**(refine_max_level-1)
37  ^d&xprobmin(^d)=xprobmin^d\
38  ^d&xprobmax(^d)=xprobmax^d\
39  ^d&domain_nx(^d)=domain_nx^d\
40  dl=(xprobmax(ndim)-xprobmin(ndim))/(domain_nx(ndim)*refine_factor)
41  ! max length of a field line
42  if (.not. mask) phys_trac_mask=xprobmax^nd
43  {^iftwod
44  lengthfl=(phys_trac_mask-xprobmin(2))*3.d0
45  }
46  {^ifthreed
47  lengthfl=(phys_trac_mask-xprobmin(3))*3.d0
48  }
49 
50  numlp=floor(lengthfl/dl)
51  numl=1
52  numfl=1
53  do j=1,ndim-1
54  ! number of field lines, every 4 finest grid, every direction
55  finegrid=mhd_trac_finegrid
56  numl(j)=floor((xprobmax(j)-xprobmin(j))/dl/finegrid)
57  numfl=numfl*numl(j)
58  end do
59  allocate(xfi(numfl,ndim))
60  xfi(:,ndim)=xprobmin(ndim)+dl/50.d0
61  {do ix^db=1,numl(^db)\}
62  ^d&ix(^d)=ix^d\
63  ifl=0
64  do j=ndim-1,1,-1
65  ifl=ifl+(ix(j)-(ndim-1-j))*(numl(j))**(ndim-1-j)
66  end do
67  xfi(ifl,1:ndim-1)=xprobmin(1:ndim-1)+finegrid*dl*ix(1:ndim-1)-finegrid*dl/2.d0
68  {end do\}
69  {^iftwod
70  if(mype .eq. 0) write(*,*) 'NOTE: 2D TRAC method take the y-dir == grav-dir'
71  }
72  {^ifthreed
73  if(mype .eq. 0) write(*,*) 'NOTE: 3D TRAC method take the z-dir == grav-dir'
74  }
75 
76  memxfi=floor(8*numfl*numlp*ndim/1e6)
77  if (mype==0) write(*,*) 'Memory requirement for each processor in TRAC:'
78  if (mype==0) write(*,*) memxfi,' MB'
79 
80  allocate(trac_pe(0:npe-1),trac_grid(max_blocks),ground_grid(max_blocks))
81 
82  end subroutine init_trac_line
83 
84  subroutine init_trac_block(mask)
85  logical, intent(in) :: mask
86  integer :: refine_factor,finegrid,iFL,j
87  integer :: ix^D,ixT^L
88  integer :: numL(ndim),ix(ndim)
89  double precision :: lengthFL
90  double precision :: ration,a0
91  double precision :: xprobmin(ndim),xprobmax(ndim),dxT(ndim)
92 
93  refine_factor=2**(refine_max_level-1)
94  ^d&xprobmin(^d)=xprobmin^d\
95  ^d&xprobmax(^d)=xprobmax^d\
96  ^d&dxt^d=(xprobmax^d-xprobmin^d)/(domain_nx^d*refine_factor/block_nx^d)\
97  ^d&dxt(ndim)=dxt^d\
98  finegrid=mhd_trac_finegrid
99  {^ifoned
100  dl=dxt1/finegrid
101  }
102  {^nooned
103  dl=min(dxt^d)/finegrid
104  }
105  ! table for interpolation
106  ^d&xtmin(^d)=xprobmin^d\
107  ^d&xtmax(^d)=xprobmax^d\
108  if(mask) xtmax(ndim)=phys_trac_mask
109  ! max length of a field line
110  if(mask) then
111  lengthfl=maxval(xtmax-xprobmin)*3.d0
112  else
113  lengthfl=maxval(xprobmax-xprobmin)*3.d0
114  end if
115  numlp=floor(lengthfl/dl)
116  ^d&numxt^d=ceiling((xtmax(^d)-xtmin(^d)-smalldouble)/dxt^d)\
117  allocate(xt(numxt^d,ndim))
118  ixtmin^d=1\
119  ixtmax^d=numxt^d\
120  {do j=1,numxt^d
121  xt(j^d%ixT^s,^d)=(j-0.5d0)*dxt^d+xtmin(^d)
122  end do\}
123  if(mask) xtmax(ndim)=maxval(xt(:^d&,ndim))+half*dxt(ndim)
124  numl=1
125  numfl=1
126  do j=1,ndim-1
127  ! number of field lines, every 4 finest grid, every direction
128  numl(j)=floor((xprobmax(j)-xprobmin(j))/dl)
129  numfl=numfl*numl(j)
130  end do
131  allocate(xfi(numfl,ndim))
132  xfi(:,ndim)=xprobmin(ndim)+dl/50.d0
133  {do ix^db=1,numl(^db)\}
134  ^d&ix(^d)=ix^d\
135  ifl=0
136  do j=ndim-1,1,-1
137  ifl=ifl+(ix(j)-(ndim-1-j))*(numl(j))**(ndim-1-j)
138  end do
139  xfi(ifl,1:ndim-1)=xprobmin(1:ndim-1)+dl*ix(1:ndim-1)-dl/2.d0
140  {end do\}
141  {^iftwod
142  if(mype .eq. 0) write(*,*) 'NOTE: 2D TRAC method take the y-dir == grav-dir'
143  }
144  {^ifthreed
145  if(mype .eq. 0) write(*,*) 'NOTE: 3D TRAC method take the z-dir == grav-dir'
146  }
147  end subroutine init_trac_block
148 
149  subroutine trac_simple(tco_global,trac_alfa,T_peak)
150  double precision, intent(in) :: tco_global, trac_alfa,T_peak
151  integer :: iigrid, igrid
152 
153  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
154  {^ifoned
155  ps(igrid)%special_values(1)=tco_global
156  }
157  if(ps(igrid)%special_values(1)<trac_alfa*ps(igrid)%special_values(2)) then
158  ps(igrid)%special_values(1)=trac_alfa*ps(igrid)%special_values(2)
159  end if
160  if(ps(igrid)%special_values(1) .lt. t_bott) then
161  ps(igrid)%special_values(1)=t_bott
162  else if(ps(igrid)%special_values(1) .gt. 0.2d0*t_peak) then
163  ps(igrid)%special_values(1)=0.2d0*t_peak
164  end if
165  ps(igrid)%wextra(ixg^t,iw_tcoff)=ps(igrid)%special_values(1)
166  !> special values(2) to save old tcutoff
167  ps(igrid)%special_values(2)=ps(igrid)%special_values(1)
168  end do
169  end subroutine trac_simple
170 
171  subroutine ltrac(T_peak)
172  double precision, intent(in) :: T_peak
173  integer :: iigrid, igrid
174  integer :: ixO^L,trac_tcoff
175 
176  ixo^l=ixm^ll^ladd1;
177  trac_tcoff=iw_tcoff
178  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
179  where(ps(igrid)%wextra(ixo^s,trac_tcoff) .lt. t_bott)
180  ps(igrid)%wextra(ixo^s,trac_tcoff)=t_bott
181  else where(ps(igrid)%wextra(ixo^s,trac_tcoff) .gt. 0.2d0*t_peak)
182  ps(igrid)%wextra(ixo^s,trac_tcoff)=0.2d0*t_peak
183  end where
184  end do
185  end subroutine ltrac
186 
187  subroutine tracl(mask,T_peak)
188  logical, intent(in) :: mask
189  double precision, intent(in) :: T_peak
190  integer :: ix^D,j
191  integer :: iigrid, igrid
192  double precision :: xF(numFL,numLP,ndim)
193  integer :: numR(numFL),ix^L
194  double precision :: Tlcoff(numFL)
195  integer :: ipel(numFL,numLP),igridl(numFL,numLP)
196  logical :: forwardl(numFL)
197 
198  tmax=t_peak
199  xf=zero
200 
201  do j=1,ndim
202  xf(:,1,j)=xfi(:,j)
203  end do
204 
205  call mpi_barrier(icomm,ierrmpi)
206  ! record pe and grid for mask region
207  call update_pegrid()
208  ! get temperature for the calculation of Te gradient,
209  ! which will be stored in wextra(ixI^S,Tcoff_) temporarily
210  call get_te_grid()
211  ! find out direction for tracing B field
212  call get_btracing_dir(ipel,igridl,forwardl)
213  ! trace Bfield and get Tcoff for each field line
214  call get_tcoff_line(xf,numr,tlcoff,ipel,igridl,forwardl,mask)
215  ! init vairable Tcoff_ and Tweight_
216  call init_trac_tcoff()
217  ! get cell Tcoff via interpolation
218  call interp_tcoff(xf,ipel,igridl,numr,tlcoff)
219  ! convert primitive data back to conserved data
220  call mpi_barrier(icomm,ierrmpi)
221 
222 
223  end subroutine tracl
224 
225  subroutine tracb(mask,T_peak)
226  logical, intent(in) :: mask
227  double precision, intent(in) :: T_peak
228  integer :: peArr(numxT^D),gdArr(numxT^D),numR(numFL)
229  double precision :: Tcoff(numxT^D),Tcmax(numxT^D),Bdir(numxT^D,ndim)
230  double precision :: xF(numFL,numLP,ndim),Tcoff_line(numFL)
231  integer :: xpe(numFL,numLP,2**ndim)
232  integer :: xgd(numFL,numLP,2**ndim)
233 
234  tmax=t_peak
235  tcoff=zero
236  tcmax=zero
237  bdir=zero
238  pearr=-1
239  gdarr=-1
240  call block_estable(mask,tcoff,tcmax,bdir,pearr,gdarr)
241  xf=zero
242  numr=0
243  tcoff_line=zero
244  call block_trace_mfl(mask,tcoff,tcoff_line,tcmax,bdir,pearr,gdarr,xf,numr,xpe,xgd)
245  call block_interp_grid(mask,xf,numr,xpe,xgd,tcoff_line)
246  end subroutine tracb
247 
248  subroutine block_estable(mask,Tcoff,Tcmax,Bdir,peArr,gdArr)
249  logical :: mask
250  double precision :: Tcoff(numxT^D),Tcoff_recv(numxT^D)
251  double precision :: Tcmax(numxT^D),Tcmax_recv(numxT^D)
252  double precision :: Bdir(numxT^D,ndim),Bdir_recv(numxT^D,ndim)
253  integer :: peArr(numxT^D),peArr_recv(numxT^D)
254  integer :: gdArr(numxT^D),gdArr_recv(numxT^D)
255  integer :: xc^L,xd^L,ix^D
256  integer :: iigrid,igrid,numxT,intab
257  double precision :: xb^L
258 
259  tcoff_recv=zero
260  tcmax_recv=zero
261  bdir_recv=zero
262  pearr_recv=-1
263  gdarr_recv=-1
264  !> combine table from different processors
265  xcmin^d=nghostcells+1\
266  xcmax^d=block_nx^d+nghostcells\
267  do iigrid=1,igridstail; igrid=igrids(iigrid);
268  ps(igrid)%wextra(:^d&,tweight_)=zero
269  ps(igrid)%wextra(:^d&,tcoff_)=zero
270  ^d&xbmin^d=rnode(rpxmin^d_,igrid)-xtmin(^d)\
271  ^d&xbmax^d=rnode(rpxmax^d_,igrid)-xtmin(^d)\
272  xdmin^d=nint(xbmin^d/dxt^d)+1\
273  xdmax^d=ceiling((xbmax^d-smalldouble)/dxt^d)\
274  {do ix^d=xdmin^d,xdmax^d \}
275  intab=0
276  {if (ix^d .le. numxt^d) intab=intab+1 \}
277  if(intab .eq. ndim) then
278  !> in principle, no overlap will happen here
279  tcoff(ix^d)=max(tcoff(ix^d),ps(igrid)%special_values(1))
280  tcmax(ix^d)=ps(igrid)%special_values(2)
281  !> abs(Bdir) <= 1, so that Bdir+2 should always be positive
282  bdir(ix^d,1:ndim)=ps(igrid)%special_values(3:3+ndim-1)+2.d0
283  pearr(ix^d)=mype
284  gdarr(ix^d)=igrid
285  end if
286  {end do\}
287  end do
288  call mpi_barrier(icomm,ierrmpi)
289  numxt={numxt^d*}
290  call mpi_allreduce(pearr,pearr_recv,numxt,mpi_integer,&
291  mpi_max,icomm,ierrmpi)
292  call mpi_allreduce(gdarr,gdarr_recv,numxt,mpi_integer,&
293  mpi_max,icomm,ierrmpi)
294  call mpi_allreduce(tcoff,tcoff_recv,numxt,mpi_double_precision,&
295  mpi_max,icomm,ierrmpi)
296  call mpi_allreduce(bdir,bdir_recv,numxt*ndim,mpi_double_precision,&
297  mpi_max,icomm,ierrmpi)
298  if(.not. mask) then
299  call mpi_allreduce(tcmax,tcmax_recv,numxt,mpi_double_precision,&
300  mpi_max,icomm,ierrmpi)
301  end if
302  pearr=pearr_recv
303  gdarr=gdarr_recv
304  tcoff=tcoff_recv
305  bdir=bdir_recv-2.d0
306  if(.not. mask) tcmax=tcmax_recv
307  end subroutine block_estable
308 
309  subroutine block_trace_mfl(mask,Tcoff,Tcoff_line,Tcmax,Bdir,peArr,gdArr,xF,numR,xpe,xgd)
310  integer :: i,j,k,k^D,ix_next^D
311  logical :: mask,flag,first
312  double precision :: Tcoff(numxT^D),Tcoff_line(numFL)
313  double precision :: Tcmax(numxT^D),Tcmax_line(numFL)
314  double precision :: xF(numFL,numLP,ndim)
315  integer :: ix_mod(ndim,2),numR(numFL)
316  double precision :: alfa_mod(ndim,2)
317  double precision :: nowpoint(ndim),nowgridc(ndim)
318  double precision :: Bdir(numxT^D,ndim)
319  double precision :: init_dir,now_dir1(ndim),now_dir2(ndim)
320  integer :: peArr(numxT^D),xpe(numFL,numLP,2**ndim)
321  integer :: gdArr(numxT^D),xgd(numFL,numLP,2**ndim)
322 
323  do i=1,numfl
324  flag=.true.
325  ^d&k^d=ceiling((xfi(i,^d)-xtmin(^d)-smalldouble)/dxt^d)\
326  tcoff_line(i)=tcoff(k^d)
327  if(.not. mask) tcmax_line(i)=tcmax(k^d)
328  ix_next^d=k^d\
329  j=1
330  xf(i,j,:)=xfi(i,:)
331  do while(flag)
332  nowpoint(:)=xf(i,j,:)
333  nowgridc(:)=xt(ix_next^d,:)
334  first=.true.
335  if(j .eq. 1) then
336  call rk_bdir(nowgridc,nowpoint,ix_next^d,now_dir1,bdir,&
337  ix_mod,first,init_dir)
338  else
339  call rk_bdir(nowgridc,nowpoint,ix_next^d,now_dir1,bdir,&
340  ix_mod,first)
341  end if
342  {^iftwod
343  xgd(i,j,1)=gdarr(ix_mod(1,1),ix_mod(2,1))
344  xgd(i,j,2)=gdarr(ix_mod(1,2),ix_mod(2,1))
345  xgd(i,j,3)=gdarr(ix_mod(1,1),ix_mod(2,2))
346  xgd(i,j,4)=gdarr(ix_mod(1,2),ix_mod(2,2))
347  xpe(i,j,1)=pearr(ix_mod(1,1),ix_mod(2,1))
348  xpe(i,j,2)=pearr(ix_mod(1,2),ix_mod(2,1))
349  xpe(i,j,3)=pearr(ix_mod(1,1),ix_mod(2,2))
350  xpe(i,j,4)=pearr(ix_mod(1,2),ix_mod(2,2))
351  }
352  {^ifthreed
353  xgd(i,j,1)=gdarr(ix_mod(1,1),ix_mod(2,1),ix_mod(3,1))
354  xgd(i,j,2)=gdarr(ix_mod(1,2),ix_mod(2,1),ix_mod(3,1))
355  xgd(i,j,3)=gdarr(ix_mod(1,1),ix_mod(2,2),ix_mod(3,1))
356  xgd(i,j,4)=gdarr(ix_mod(1,2),ix_mod(2,2),ix_mod(3,1))
357  xgd(i,j,5)=gdarr(ix_mod(1,1),ix_mod(2,1),ix_mod(3,2))
358  xgd(i,j,6)=gdarr(ix_mod(1,2),ix_mod(2,1),ix_mod(3,2))
359  xgd(i,j,7)=gdarr(ix_mod(1,1),ix_mod(2,2),ix_mod(3,2))
360  xgd(i,j,8)=gdarr(ix_mod(1,2),ix_mod(2,2),ix_mod(3,2))
361  xpe(i,j,1)=pearr(ix_mod(1,1),ix_mod(2,1),ix_mod(3,1))
362  xpe(i,j,2)=pearr(ix_mod(1,2),ix_mod(2,1),ix_mod(3,1))
363  xpe(i,j,3)=pearr(ix_mod(1,1),ix_mod(2,2),ix_mod(3,1))
364  xpe(i,j,4)=pearr(ix_mod(1,2),ix_mod(2,2),ix_mod(3,1))
365  xpe(i,j,5)=pearr(ix_mod(1,1),ix_mod(2,1),ix_mod(3,2))
366  xpe(i,j,6)=pearr(ix_mod(1,2),ix_mod(2,1),ix_mod(3,2))
367  xpe(i,j,7)=pearr(ix_mod(1,1),ix_mod(2,2),ix_mod(3,2))
368  xpe(i,j,8)=pearr(ix_mod(1,2),ix_mod(2,2),ix_mod(3,2))
369  }
370  nowpoint(:)=nowpoint(:)+init_dir*now_dir1*dl
371  {if(nowpoint(^d) .gt. xtmax(^d) .or. nowpoint(^d) .lt. xtmin(^d)) then
372  flag=.false.
373  end if\}
374  if(mask .and. nowpoint(ndim) .gt. phys_trac_mask) then
375  flag=.false.
376  end if
377  if(flag) then
378  first=.false.
379  ^d&ix_next^d=ceiling((nowpoint(^d)-xtmin(^d)-smalldouble)/dxt^d)\
380  nowgridc(:)=xt(ix_next^d,:)
381  call rk_bdir(nowgridc,nowpoint,ix_next^d,now_dir2,bdir,&
382  ix_mod,first)
383  xf(i,j+1,:)=xf(i,j,:)+init_dir*dl*half*(now_dir1+now_dir2)
384  {if(xf(i,j+1,^d) .gt. xtmax(^d) .or. xf(i,j+1,^d) .lt. xtmin(^d)) then
385  flag=.false.
386  end if\}
387  if(mask .and. xf(i,j+1,ndim) .gt. phys_trac_mask) then
388  flag=.false.
389  end if
390  if(flag) then
391  ^d&ix_next^d=ceiling((xf(i,j+1,^d)-xtmin(^d)-smalldouble)/dxt^d)\
392  j=j+1
393  tcoff_line(i)=max(tcoff_line(i),tcoff(ix_next^d))
394  if(.not.mask) tcmax_line(i)=max(tcmax_line(i),tcmax(ix_next^d))
395  end if
396  end if
397  end do
398  numr(i)=j
399  if(mask) then
400  if(tcoff_line(i) .gt. tmax*0.2d0) then
401  tcoff_line(i)=tmax*0.2d0
402  end if
403  else
404  if(tcoff_line(i) .gt. tcmax_line(i)*0.2d0) then
405  tcoff_line(i)=tcmax_line(i)*0.2d0
406  end if
407  end if
408  end do
409  end subroutine block_trace_mfl
410 
411  subroutine rk_bdir(nowgridc,nowpoint,ix_next^D,now_dir,Bdir,ix_mod,first,init_dir)
412  double precision :: nowpoint(ndim),nowgridc(ndim)
413  integer :: ix_mod(ndim,2)
414  double precision :: alfa_mod(ndim,2)
415  integer :: ix_next^D,k^D
416  double precision :: now_dir(ndim)
417  double precision :: Bdir(numxT^D,ndim)
418  logical :: first
419  double precision, optional :: init_dir
420 
421  {if(nowpoint(^d) .gt. xtmin(^d)+half*dxt^d .and. nowpoint(^d) .lt. xtmax(^d)-half*dxt^d) then
422  if(nowpoint(^d) .le. nowgridc(^d)) then
423  ix_mod(^d,1)=ix_next^d-1
424  ix_mod(^d,2)=ix_next^d
425  alfa_mod(^d,1)=abs(nowgridc(^d)-nowpoint(^d))/dxt^d
426  alfa_mod(^d,2)=one-alfa_mod(^d,1)
427  else
428  ix_mod(^d,1)=ix_next^d
429  ix_mod(^d,2)=ix_next^d+1
430  alfa_mod(^d,2)=abs(nowgridc(^d)-nowpoint(^d))/dxt^d
431  alfa_mod(^d,1)=one-alfa_mod(^d,2)
432  end if
433  else
434  ix_mod(^d,:)=ix_next^d
435  alfa_mod(^d,:)=half
436  end if\}
437  now_dir=zero
438  {^iftwod
439  do k1=1,2
440  do k2=1,2
441  now_dir=now_dir + bdir(ix_mod(1,k1),ix_mod(2,k2),:)*alfa_mod(1,k1)*alfa_mod(2,k2)
442  end do
443  end do
444  }
445  {^ifthreed
446  do k1=1,2
447  do k2=1,2
448  do k3=1,2
449  now_dir=now_dir + bdir(ix_mod(1,k1),ix_mod(2,k2),ix_mod(3,k3),:)&
450  *alfa_mod(1,k1)*alfa_mod(2,k2)*alfa_mod(3,k3)
451  end do
452  end do
453  end do
454  }
455  if(present(init_dir)) then
456  init_dir=sign(one,now_dir(ndim))
457  end if
458  end subroutine rk_bdir
459 
460  subroutine block_interp_grid(mask,xF,numR,xpe,xgd,Tcoff_line)
461  logical :: mask
462  double precision :: xF(numFL,numLP,ndim)
463  integer :: numR(numFL)
464  integer :: xpe(numFL,numLP,2**ndim)
465  integer :: xgd(numFL,numLP,2**ndim)
466  double precision :: Tcoff_line(numFL)
467  double precision :: weightIndex,weight,ds
468  integer :: i,j,k,igrid,iigrid,ixO^L,ixc^L,ixc^D
469  double precision :: dxMax^D,dxb^D
470 
471  !> interpolate lines into grid cells
472  weightindex=2.d0
473  dxmax^d=2.d0*dl;
474  ixo^l=ixm^ll;
475  do i=1,numfl
476  do j=1,numr(i)
477  do k=1,2**ndim
478  if(mype .eq. xpe(i,j,k)) then
479  igrid=xgd(i,j,k)
480  if(igrid .le. igrids(igridstail)) then
481  ^d&dxb^d=rnode(rpdx^d_,igrid)\
482  ^d&ixcmin^d=floor((xf(i,j,^d)-dxmax^d-ps(igrid)%x(ixomin^dd,^d))/dxb^d)+ixomin^d\
483  ^d&ixcmax^d=floor((xf(i,j,^d)+dxmax^d-ps(igrid)%x(ixomin^dd,^d))/dxb^d)+ixomin^d\
484  {if (ixcmin^d<ixomin^d) ixcmin^d=ixomin^d\}
485  {if (ixcmax^d>ixomax^d) ixcmax^d=ixomax^d\}
486  {do ixc^d=ixcmin^d,ixcmax^d\}
487  ds=0.d0
488  {ds=ds+(xf(i,j,^d)-ps(igrid)%x(ixc^dd,^d))**2\}
489  ds=sqrt(ds)
490  if(ds .le. 0.099d0*dl) then
491  weight=(1/(0.099d0*dl))**weightindex
492  else
493  weight=(1/ds)**weightindex
494  endif
495  ps(igrid)%wextra(ixc^d,tweight_)=ps(igrid)%wextra(ixc^d,tweight_)+weight
496  ps(igrid)%wextra(ixc^d,tcoff_)=ps(igrid)%wextra(ixc^d,tcoff_)+weight*tcoff_line(i)
497  {enddo\}
498  else
499  call mpistop("we need to check here 366Line in mod_trac.t")
500  end if
501  end if
502  end do
503  end do
504  end do
505  ! finish interpolation
506  do iigrid=1,igridstail; igrid=igrids(iigrid);
507  where (ps(igrid)%wextra(ixo^s,tweight_)>0.d0)
508  ps(igrid)%wextra(ixo^s,tcoff_)=ps(igrid)%wextra(ixo^s,tcoff_)/ps(igrid)%wextra(ixo^s,tweight_)
509  elsewhere
510  ps(igrid)%wextra(ixo^s,tcoff_)=0.2d0*tmax
511  endwhere
512  enddo
513  end subroutine block_interp_grid
514 
515  ! TODO remove, not used?
516 ! subroutine get_Tmax_grid(x,w,ixI^L,ixO^L,Tmax_grid)
517 ! integer, intent(in) :: ixI^L,ixO^L
518 ! double precision, intent(in) :: x(ixI^S,1:ndim)
519 ! double precision, intent(out) :: w(ixI^S,1:nw)
520 ! double precision :: Tmax_grid
521 ! double precision :: pth(ixI^S),Te(ixI^S)
522 !
523 ! call mhd_get_pthermal(w,x,ixI^L,ixO^L,pth)
524 ! Te(ixO^S)=pth(ixO^S)/w(ixO^S,rho_)
525 ! Tmax_grid=maxval(Te(ixO^S))
526 ! end subroutine get_Tmax_grid
527 
528  subroutine init_trac_tcoff()
529  integer :: ixI^L,ixO^L,igrid,iigrid
530 
531  ixi^l=ixg^ll;
532  ixo^l=ixm^ll;
533 
534  do iigrid=1,igridstail; igrid=igrids(iigrid);
535  ps(igrid)%wextra(ixi^s,tcoff_)=0.d0
536  ps(igrid)%wextra(ixi^s,tweight_)=0.d0
537  end do
538 
539  end subroutine init_trac_tcoff
540 
541  subroutine update_pegrid()
542 
543  ngrid_trac=0
544  ngrid_ground=0
545  trac_pe(:)=.false.
546  trac_grid(:)=0
547  ground_grid(:)=0
548 
549  ! traverse gridtable to information of grid inside mask region
550  call traverse_gridtable()
551 
552  end subroutine update_pegrid
553 
554  subroutine traverse_gridtable()
556 
557  double precision :: dxb^D,xb^L
558  integer :: iigrid,igrid,j
559  logical, allocatable :: trac_pe_recv(:)
560  double precision :: hcmax_bt
561 
562  allocate(trac_pe_recv(npe))
563 
564  hcmax_bt=xprobmin^nd+(xprobmax^nd-xprobmin^nd)/(domain_nx^nd*2**(refine_max_level-1))
565 
566  do iigrid=1,igridstail; igrid=igrids(iigrid);
567  xbmin^nd=rnode(rpxmin^nd_,igrid)
568  if (xbmin^nd<phys_trac_mask) then
569  trac_pe(mype)=.true.
570  ngrid_trac=ngrid_trac+1
571  trac_grid(ngrid_trac)=igrid
572  if (xbmin^nd<hcmax_bt) then
573  ngrid_ground=ngrid_ground+1
574  ground_grid(ngrid_ground)=igrid
575  endif
576  endif
577  enddo
578 
579  call mpi_allreduce(trac_pe,trac_pe_recv,npe,mpi_logical,mpi_lor,icomm,ierrmpi)
580  trac_pe=trac_pe_recv
581 
582  deallocate(trac_pe_recv)
583 
584  end subroutine traverse_gridtable
585 
586  subroutine get_te_grid()
587  integer :: ixI^L,ixO^L,igrid,iigrid,j
588 
589  ixi^l=ixg^ll;
590  ixo^l=ixm^ll;
591 
592  do j=1,ngrid_trac
593  igrid=trac_grid(j)
594 
595  call mhd_get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixi^l,ps(igrid)%wextra(ixi^s,tcoff_))
596  !TODO move outside loop for optimziation?
597  if(has_equi_rho0) then
598  ps(igrid)%wextra(ixi^s,tcoff_)=ps(igrid)%wextra(ixi^s,tcoff_)/&
599  (ps(igrid)%w(ixi^s,rho_) + ps(igrid)%equi_vars(ixi^s,equi_rho0_,0))
600  else
601  ps(igrid)%wextra(ixi^s,tcoff_)=ps(igrid)%wextra(ixi^s,tcoff_)/ps(igrid)%w(ixi^s,rho_)
602  endif
603  enddo
604 
605  end subroutine get_te_grid
606 
607  subroutine get_btracing_dir(ipel,igridl,forwardl)
608  integer :: ipel(numFL,numLP),igridl(numFL,numLP)
609  logical :: forwardl(numFL)
610 
611  integer :: igrid,ixO^L,iFL,j,ix^D,idir,ixb^D,ixbb^D
612  double precision :: xb^L,dxb^D,xd^D,factor,Bh
613  integer :: numL(ndim),ixmin(ndim),ixmax(ndim),ix(ndim)
614  logical :: forwardRC(numFL)
615 
616  ipel=-1
617  igridl=0
618  forwardl=.true.
619 
620  ^d&ixomin^d=ixmlo^d;
621  ^d&ixomax^d=ixmhi^d;
622 
623  do j=1,ngrid_ground
624  igrid=ground_grid(j)
625  ^d&dxb^d=rnode(rpdx^d_,igrid);
626  ^d&xbmin^d=rnode(rpxmin^d_,igrid);
627  ^d&xbmax^d=rnode(rpxmax^d_,igrid);
628  ^d&ixmin(^d)=floor((xbmin^d-xprobmin^d)/(mhd_trac_finegrid*dl))+1;
629  ^d&ixmax(^d)=floor((xbmax^d-xprobmin^d)/(mhd_trac_finegrid*dl));
630  ^d&numl(^d)=floor((xprobmax^d-xprobmin^d)/(mhd_trac_finegrid*dl));
631  ixmin(^nd)=1
632  ixmax(^nd)=1
633  numl(^nd)=1
634 
635  {do ix^db=ixmin(^db),ixmax(^db)\}
636  ^d&ix(^d)=ix^d;
637  ifl=0
638  do idir=ndim-1,1,-1
639  ifl=ifl+(ix(idir)-(ndim-1-idir))*(numl(idir))**(ndim-1-idir)
640  enddo
641  ipel(ifl,1)=mype
642  igridl(ifl,1)=igrid
643 
644  ^d&ixb^d=floor((xfi(ifl,^d)-ps(igrid)%x(ixomin^dd,^d))/dxb^d)+ixomin^d;
645  ^d&xd^d=(xfi(ifl,^d)-ps(igrid)%x(ixb^dd,^d))/dxb^d;
646  bh=0.d0
647  {do ixbb^d=0,1\}
648  factor={abs(1-ix^d-xd^d)*}
649  if (b0field) then
650  bh=bh+factor*(ps(igrid)%w(ixb^d+ixbb^d,mag(^nd))+ps(igrid)%B0(ixb^d+ixbb^d,^nd,0))
651  else
652  bh=bh+factor*ps(igrid)%w(ixb^d+ixbb^d,mag(^nd))
653  endif
654  {enddo\}
655  if (bh>0) then
656  forwardl(ifl)=.true.
657  else
658  forwardl(ifl)=.false.
659  endif
660  {enddo\}
661  enddo
662 
663  call mpi_allreduce(forwardl,forwardrc,numfl,mpi_logical,&
664  mpi_land,icomm,ierrmpi)
665  forwardl=forwardrc
666 
667  end subroutine get_btracing_dir
668 
669  subroutine get_tcoff_line(xFL,numR,TcoffFL,ipeFL,igridFL,forwardFL,mask)
670  use mod_trace_field
671 
672  double precision :: xFL(numFL,numLP,ndim)
673  integer :: numR(numFL)
674  double precision :: TcoffFL(numFL),TmaxFL(numFL)
675  integer :: ipeFL(numFL,numLP),igridFL(numFL,numLP)
676  logical :: forwardFL(numFL)
677  logical, intent(in) :: mask
678 
679  integer :: nwP,nwL,iFL,iLP
680  double precision :: wPm(numFL,numLP,2),wLm(numFL,1+2)
681  character(len=std_len) :: ftype,tcondi
682 
683  nwp=2
684  nwl=2
685  ftype='Bfield'
686  tcondi='TRAC'
687  call trace_field_multi(xfl,wpm,wlm,dl,numfl,numlp,nwp,nwl,forwardfl,ftype,tcondi)
688  do ifl=1,numfl
689  numr(ifl)=int(wlm(ifl,1))
690  tcofffl(ifl)=wlm(ifl,2)
691  tmaxfl(ifl)=wlm(ifl,3)
692  if(mask) then
693  if(tcofffl(ifl)>0.2d0*tmax) tcofffl(ifl)=0.2d0*tmax
694  else
695  tmaxfl(ifl)=wlm(ifl,3)
696  if(tcofffl(ifl)>0.2d0*tmaxfl(ifl)) tcofffl(ifl)=0.2d0*tmaxfl(ifl)
697  end if
698 
699  if(tcofffl(ifl)<t_bott) tcofffl(ifl)=t_bott
700  enddo
701 
702  do ifl=1,numfl
703  if (numr(ifl)>0) then
704  do ilp=1,numr(ifl)
705  ipefl(ifl,ilp)=int(wpm(ifl,ilp,1))
706  igridfl(ifl,ilp)=int(wpm(ifl,ilp,2))
707  enddo
708  endif
709  enddo
710 
711 
712  end subroutine get_tcoff_line
713 
714  subroutine interp_tcoff(xF,ipel,igridl,numR,Tlcoff)
715  double precision :: xF(numFL,numLP,ndim)
716  integer :: numR(numFL),ipel(numFL,numLP),igridl(numFL,numLP)
717  double precision :: Tlcoff(numFL)
718 
719  integer :: iFL,iLP,ixO^L,ixI^L,ixc^L,ixb^L,ixc^D
720  integer :: igrid,j,ipmin,ipmax,igrid_nb
721  double precision :: dxb^D,dxMax^D,xb^L,Tcnow
722  double precision :: xFnow(ndim)
723  integer :: weightIndex,idn^D,ixmax^ND
724  double precision :: ds,weight
725 
726  weightindex=2
727 
728  ^d&iximin^d=ixglo^d;
729  ^d&iximax^d=ixghi^d;
730  ^d&ixomin^d=ixmlo^d;
731  ^d&ixomax^d=ixmhi^d;
732 
733  do ifl=1,numfl
734  ilp=1
735  tcnow=tlcoff(ifl)
736  do while (ilp<=numr(ifl))
737  ! find points in the same grid
738  ipmin=ilp
739  do while (ipel(ifl,ipmin)/=mype .and. ipmin<=numr(ifl))
740  ipmin=ipmin+1
741  enddo
742  igrid=igridl(ifl,ipmin)
743  ipmax=ipmin
744  do while (ipel(ifl,ipmax)==mype .and. igridl(ifl,ipmax+1)==igrid .and. ipmax<numr(ifl))
745  ipmax=ipmax+1
746  enddo
747 
748  ! parameters for the grid
749  ^d&dxb^d=rnode(rpdx^d_,igrid);
750  ^d&dxmax^d=4*dxb^d;
751 
752  do ilp=ipmin,ipmax
753  xfnow(:)=xf(ifl,ilp,:)
754  ^d&ixbmin^d=floor((xfnow(^d)-dxmax^d-ps(igrid)%x(ixomin^dd,^d))/dxb^d)+ixomin^d;
755  ^d&ixbmax^d=floor((xfnow(^d)+dxmax^d-ps(igrid)%x(ixomin^dd,^d))/dxb^d)+ixomin^d;
756 
757  ! do interpolation inside the grid
758  {ixcmin^d=max(ixbmin^d,ixomin^d)\}
759  {ixcmax^d=min(ixbmax^d,ixomax^d)\}
760  xbmin^nd=rnode(rpxmin^nd_,igrid)
761  xbmax^nd=rnode(rpxmax^nd_,igrid)
762  ixmax^nd=floor((phys_trac_mask-xbmin^nd)/dxb^nd)+ixomin^nd
763  if (xbmax^nd>phys_trac_mask) ixcmax^nd=min(ixmax^nd,ixcmax^nd)
764  {do ixc^d=ixcmin^d,ixcmax^d\}
765  ds=0.d0
766  {ds=ds+(xfnow(^d)-ps(igrid)%x(ixc^dd,^d))**2\}
767  ds=sqrt(ds)
768  if(ds<1.0d-2*dxb1) then
769  weight=(1/(1.0d-2*dxb1))**weightindex
770  else
771  weight=(1/ds)**weightindex
772  endif
773  ps(igrid)%wextra(ixc^d,tweight_)=ps(igrid)%wextra(ixc^d,tweight_)+weight
774  ps(igrid)%wextra(ixc^d,tcoff_)=ps(igrid)%wextra(ixc^d,tcoff_)+weight*tcnow
775  {enddo\}
776 
777  ! do interpolation in neighbor grids that have the same level and pe
778  {
779  if (ixbmin^d<ixomin^d) then
780  idn^dd=0;
781  idn^d=-1
782  if (neighbor(2,idn^dd,igrid)==mype .and. neighbor_type(idn^dd,igrid)==neighbor_sibling) then
783  igrid_nb=neighbor(1,idn^dd,igrid)
784  ixcmin^dd=max(ixbmin^dd,ixomin^dd);
785  ixcmax^dd=min(ixbmax^dd,ixomax^dd);
786  ixcmin^d=ixomax^d+(ixbmin^d-ixomin^d)
787  ixcmax^d=ixomax^d
788  xbmin^nd=rnode(rpxmin^nd_,igrid_nb)
789  xbmax^nd=rnode(rpxmax^nd_,igrid_nb)
790  ixmax^nd=floor((phys_trac_mask-xbmin^nd)/dxb^nd)+ixomin^nd
791  if (xbmax^nd>phys_trac_mask) ixcmax^nd=min(ixmax^nd,ixcmax^nd)
792 
793  {do ixc^dd=ixcmin^dd,ixcmax^dd;}
794  ds=0.d0
795  {ds=ds+(xfnow(^dd)-ps(igrid_nb)%x({ixc^dd},^dd))**2;}
796  ds=sqrt(ds)
797  if(ds<1.0d-2*dxb1) then
798  weight=(1/(1.0d-2*dxb1))**weightindex
799  else
800  weight=(1/ds)**weightindex
801  endif
802  ps(igrid_nb)%wextra(ixc^dd,tweight_)=ps(igrid_nb)%wextra(ixc^dd,tweight_)+weight
803  ps(igrid_nb)%wextra(ixc^dd,tcoff_)=ps(igrid_nb)%wextra(ixc^dd,tcoff_)+weight*tcnow
804  {enddo;}
805  endif
806  endif
807 
808  if (ixbmax^d>ixomin^d) then
809  idn^dd=0;
810  idn^d=1
811  if (neighbor(2,idn^dd,igrid)==mype .and. neighbor_type(idn^dd,igrid)==neighbor_sibling) then
812  igrid_nb=neighbor(1,idn^dd,igrid)
813  xbmin^nd=rnode(rpxmin^nd_,igrid_nb)
814  if (xbmin^nd<phys_trac_mask) then
815  ixcmin^dd=max(ixbmin^dd,ixomin^dd);
816  ixcmax^dd=min(ixbmax^dd,ixomax^dd);
817  ixcmin^d=ixomin^d
818  ixcmax^d=ixomin^d+(ixbmax^d-ixomax^d)
819  xbmax^nd=rnode(rpxmax^nd_,igrid_nb)
820  ixmax^nd=floor((phys_trac_mask-xbmin^nd)/dxb^nd)+ixomin^nd
821  if (xbmax^nd>phys_trac_mask) ixcmax^nd=min(ixmax^nd,ixcmax^nd)
822 
823  {do ixc^dd=ixcmin^dd,ixcmax^dd;}
824  ds=0.d0
825  {ds=ds+(xfnow(^dd)-ps(igrid_nb)%x({ixc^dd},^dd))**2;}
826  ds=sqrt(ds)
827  if(ds<1.0d-2*dxb1) then
828  weight=(1/(1.0d-2*dxb1))**weightindex
829  else
830  weight=(1/ds)**weightindex
831  endif
832  ps(igrid_nb)%wextra(ixc^dd,tweight_)=ps(igrid_nb)%wextra(ixc^dd,tweight_)+weight
833  ps(igrid_nb)%wextra(ixc^dd,tcoff_)=ps(igrid_nb)%wextra(ixc^dd,tcoff_)+weight*tcnow
834  {enddo;}
835  endif
836  endif
837  endif
838  \}
839  enddo
840 
841  ilp=ipmax+1
842  enddo
843  enddo
844 
845 
846  do j=1,ngrid_trac
847  igrid=trac_grid(j)
848  where(ps(igrid)%wextra(ixo^s,tweight_)>0.d0)
849  ps(igrid)%wextra(ixo^s,tcoff_)=ps(igrid)%wextra(ixo^s,tcoff_)/ps(igrid)%wextra(ixo^s,tweight_)
850  elsewhere
851  ps(igrid)%wextra(ixo^s,tcoff_)=t_bott
852  endwhere
853  enddo
854 
855  end subroutine interp_tcoff
856 
857 
858  subroutine trac_after_setdt(tco,trac_alfa,T_peak, T_bott_in)
859  double precision, intent(in) :: trac_alfa,tco,T_peak, T_bott_in
860  ! default lower limit of cutoff temperature
861  t_bott = t_bott_in
862  select case(phys_trac_type)
863  case(0)
864  !> Test case, fixed cutoff temperature
865  !> do nothing here
866  case(1)
867  !> 1D TRAC method
868  call trac_simple(tco,trac_alfa,t_peak)
869  case(2)
870  !> LTRAC method, by iijima et al. 2021
871  !> do nothing here
872  call ltrac(t_peak)
873  case(3)
874  !> 2D or 3D TRACL(ine) method
875  call tracl(.false.,t_peak)
876  case(4)
877  !> 2D or 3D TRACB(lock) method
878  call tracb(.false.,t_peak)
879  case(5)
880  !> 2D or 3D TRACL(ine) method with mask
881  call tracl(.true.,t_peak)
882  case(6)
883  !> 2D or 3D TRACB(lock) method with mask
884  call tracb(.true.,t_peak)
885  case default
886  call mpistop("undefined TRAC method type")
887  end select
888 
889  end subroutine trac_after_setdt
890 
894 
895  if(phys_trac) then
896  phys_trac_after_setdt => trac_after_setdt
897  if(phys_trac_type .eq. 3) then
898  if(mype .eq. 0) write(*,*) 'Using TRACL(ine) global method'
899  if(mype .eq. 0) write(*,*) 'By default, magnetic field lines are traced every 4 grid cells'
900  call init_trac_line(.false.)
901  end if
902  if(phys_trac_type .eq. 4) then
903  if(mype .eq. 0) write(*,*) 'Using TRACB(lock) global method'
904  if(mype .eq. 0) write(*,*) 'Currently, only valid in Cartesian uniform settings'
905  if(mype .eq. 0) write(*,*) 'By default, magnetic field lines are traced every 4 grid cells'
906  call init_trac_block(.false.)
907  end if
908  if(phys_trac_type .eq. 5) then
909  if(mype .eq. 0) write(*,*) 'Using TRACL(ine) method with a mask'
910  if(mype .eq. 0) write(*,*) 'By default, magnetic field lines are traced every 4 grid cells'
911  call init_trac_line(.true.)
912  end if
913  if(phys_trac_type .eq. 6) then
914  if(mype .eq. 0) write(*,*) 'Using TRACB(lock) method with a mask'
915  if(mype .eq. 0) write(*,*) 'Currently, only valid in Cartesian uniform settings'
916  if(mype .eq. 0) write(*,*) 'By default, magnetic field lines are traced every 4 grid cells'
917  call init_trac_block(.true.)
918  end if
919  end if
920 
921  end subroutine initialize_trac_after_settree
922 
923 
924 
925 end module mod_trac
integer, dimension(:), allocatable, public mag
Indices of the magnetic field.
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer domain_nx
number of cells for each dimension in level-one mesh
double precision phys_trac_mask
integer icomm
The MPI communicator.
integer mype
The rank of the current MPI task.
integer ierrmpi
A global MPI error return code.
integer npe
The number of MPI tasks.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
logical phys_trac
Use TRAC (Johnston 2019 ApJL, 873, L22) for MHD or 1D HD.
integer refine_max_level
Maximal number of AMR levels.
integer max_blocks
The maximum number of grid blocks in a processor.
Definition: mod_mhd.t:1
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
procedure(sub_trac_after_setdt), pointer phys_trac_after_setdt
Definition: mod_physics.t:64
subroutine get_te_grid()
Definition: mod_trac.t:587
subroutine interp_tcoff(xF, ipel, igridl, numR, Tlcoff)
Definition: mod_trac.t:715
subroutine, public initialize_trac_after_settree
Definition: mod_trac.t:892
subroutine trace_field_multi(xfm, wPm, wLm, dL, numL, numP, nwP, nwL, forwardm, ftype, tcondi)