MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_trace_field.t
Go to the documentation of this file.
4 implicit none
5
6contains
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 double precision :: x3d(3),statusF(4+ndim),statusL(numL,4+ndim+nwL),statusS(numL,4+ndim+nwL)
43 double precision :: xf(numP,ndim),wP(numP,nwP),wL(1+nwL)
44 double precision, allocatable :: data_send(:,:,:),data_recv(:,:,:)
45 integer :: indomain,ipe_now,igrid_now,igrid,j,iL
46 integer :: ipoint_in,ipoint_out,numSend,nRT,nRTmax
47 integer :: ipointm(numL),igridm(numL)
48 logical :: continueL(numL),myL(numL)
49 logical :: stopT,forward
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
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 double precision :: x3d(3),statusF(4+ndim),status_bcast(4+ndim+nwL)
208 double precision, allocatable :: data_send(:,:),data_recv(:,:)
209 integer :: indomain,ipoint_in,ipe_now,igrid_now,igrid,j
210 integer :: ipoint_out,ipe_next,igrid_next,numRT
211 logical :: stopT
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 double precision :: xfout(ndim)
305 integer :: ipe_next,igrid_next,ip_in,ip_out,j,indomain
306 logical :: newpe,stopT
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
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 double precision :: xbn^L
367 integer :: idn^D,my_neighbor_type,inblock
368 integer :: ic^D,inc^D,ipe_neighbor,igrid_neighbor
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)
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 double precision :: field(ixg^T,ndir)
441 double precision :: xs1(ndim),xs2(ndim),K1(ndim),K2(ndim)
442 double precision :: xfpre(ndim),xfnow(ndim),xfnext(ndim)
443 double precision :: Tpre,Tnow,Tnext,dTds,Lt,Lr,ds,T_bott,trac_delta
444 integer :: ip,inblock,ixI^L,ixO^L,j
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)
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 double precision :: xd^D
544 double precision :: field(0:1^D&,ndir),Fx(ndim),factor(0:1^D&)
545 double precision :: vector(ixI^S,1:ndir)
546 double precision :: Ftotal
547 integer :: ixb^D,ix^D,ixbl^D,ixO^L,j
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 double precision :: xd^D
618 double precision :: factor(0:1^D&),Tnear(0:1^D&)
619 integer :: ixb^D,ix^D,ixbl^D,j,ixO^L
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
639end 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...
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.
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:65
subroutine trace_field_single(xf, wp, wl, dl, nump, nwp, nwl, forward, ftype, tcondi)
subroutine trace_field_multi(xfm, wpm, wlm, dl, numl, nump, nwp, nwl, forwardm, ftype, tcondi)
subroutine get_k(xfn, igrid, k, ixil, dxbd, ftype)
subroutine find_points_in_pe(igrid, ipoint_in, xf, wp, wl, dl, nump, nwp, nwl, forward, ftype, tcondi, statusf)
subroutine find_next_grid(igrid, igrid_next, ipe_next, xf1, newpe, stopt)
subroutine get_t_loc_trac(igrid, xloc, tloc, ixil, dxbd)
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