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