MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_magnetofriction.t
Go to the documentation of this file.
1!> module mod_magnetofriction.t
2!> Purpose: use magnetofrictional method to relax 3D magnetic field to
3!> force-free field
4!> 01.04.2016 developed by Chun Xia and Yang Guo
5!> 04.10.2017 modulized by Chun Xia
6!> Usage:
7!> in amrvac.par:
8!> &methodlist
9!> time_stepper='onestep' ! time marching scheme, or 'twostep','threestep'
10!> flux_method=13*'cd4' ! or 'tvdlf', 'fd'
11!> limiter= 13*'koren' ! or 'vanleer','cada3','mp5' so on
12!> /
13!> &meshlist
14!> ditregrid=20 ! set iteration interval for adjusting AMR
15!> /
16!> &mhd_list
17!> mhd_magnetofriction=.true.
18!> /
19!> &mf_list
20!> mf_it_max=60000 ! set the maximum iteration number
21!> mf_ditsave=20000 ! set iteration interval for data output
22!> mf_cc=0.3 ! stability coefficient controls numerical stability
23!> mf_cy=0.2 ! frictional velocity coefficient
24!> mf_cdivb=0.01 ! divb cleaning coefficient controls diffusion speed of divb
25!> /
27 implicit none
28
29 !> stability coefficient controls numerical stability
30 double precision :: mf_cc
31 !> frictional velocity coefficient
32 double precision :: mf_cy, mf_cy_max
33 !> divb cleaning coefficient controls diffusion speed of divb
34 double precision :: mf_cdivb, mf_cdivb_max
35 !> TVDLF dissipation coefficient controls the dissipation term
36 double precision :: mf_tvdlfeps, mf_tvdlfeps_min
37 !> time in magnetofriction process
38 double precision :: tmf
39 !> maximal speed for fd scheme
40 double precision :: cmax_mype
41 !> maximal speed for fd scheme
42 double precision :: cmax_global
43 !> maximal limit of magnetofrictional velocity in cm s^-1 (Pomoell 2019 A&A)
44 double precision, public :: mf_vmax = 3.d6
45
46 !> Index of the density (in the w array)
47 integer, private, protected :: rho_
48
49 !> Indices of the momentum density
50 integer, allocatable, private, protected :: mom(:)
51
52 !> Indices of the magnetic field
53 integer, allocatable, private, protected :: mag(:)
54
55 integer :: mf_ditsave
56 integer :: mf_it_max
57 integer :: mf_it
58 logical :: mf_advance
59 logical :: fix_conserve_at_step = .true.
60
61contains
62 !> Read this module"s parameters from a file
63 subroutine mf_params_read(files)
65 character(len=*), intent(in) :: files(:)
66 integer :: n
67
68 namelist /mf_list/ mf_ditsave, mf_it_max, mf_it, mf_cc, mf_cy, mf_cy_max, &
70
71 do n = 1, size(files)
72 open(unitpar, file=trim(files(n)), status="old")
73 read(unitpar, mf_list, end=111)
74111 close(unitpar)
75 end do
76
77 end subroutine mf_params_read
78
79 !> Initialize the module
83
84 rho_=iw_rho
85 allocate(mom(ndir))
86 mom=iw_mom
87 allocate(mag(ndir))
88 mag=iw_mag
89
90 mf_it=0 ! set the initial iteration number
91 mf_it_max=60000 ! set the maximum iteration number
92 mf_ditsave=20000 ! set iteration interval for data output
93 mf_cc=0.3d0 ! stability coefficient controls numerical stability
94 mf_cy=0.2d0 ! frictional velocity coefficient. The default value is mf_cy itself.
95 mf_cy_max=mf_cy ! maximum of the frictional velocity coefficient
96 mf_cdivb=0.01d0 ! divb cleaning coefficient controls diffusion speed of divb
97 mf_cdivb_max=mf_cdivb ! maximum of the divb cleaning coefficient
98 mf_tvdlfeps=1.d0 ! coefficient to control the TVDLF dissipation
99 mf_tvdlfeps_min = mf_tvdlfeps ! minimum of the TVDLF dissipation coefficient
100 ! get dimensionless maximal mf velocity limit
102
104
108
110
111 end subroutine magnetofriction_init
112
115 use mod_physics
118 use mod_amr_grid, only: resettree
119 use mod_comm_lib, only: mpistop
120
121 double precision :: dvolume(ixG^T),dsurface(ixG^T),dvone
122 double precision :: dtfff,dtfff_pe,dtnew,dx^D
123 double precision :: cwsin_theta_new,cwsin_theta_old
124 double precision :: sum_jbb,sum_jbb_ipe,sum_j,sum_j_ipe,sum_l_ipe,sum_l
125 double precision :: f_i_ipe,f_i,volumepe,volume,tmpt,time_in
126 double precision, external :: integral_grid
127 integer :: i,iigrid, igrid, idims,ix^D,hxM^LL,fhmf,tmpit,i^D
128 logical :: patchwi(ixG^T), stagger_flag
129
130 ! not do fix conserve and getbc for staggered values if stagger is used
131 stagger_flag=stagger_grid
132 stagger_grid=.false.
133
134 time_in=mpi_wtime()
135 if(mype==0) write(*,*) 'Evolving to force-free field using magnetofricitonal method...'
136 if(prolongprimitive) call mpistop('use prolongprimitive=.false. in MF module')
137 mf_advance=.false.
138 dtfff=1.d-2
139 tmpt=global_time
140 tmpit=it
142 i=mf_it
143 if(snapshotini==0 .and. i==0) then
144 call saveamrfile(1)
145 call saveamrfile(2)
146 end if
147 mf_advance=.true.
148 ! point bc mpi datatype to partial type for magnetic field
155 ! create bc mpi datatype for ghostcells update
156 call create_bc_mpi_datatype(mag(1),ndir)
157 ! point bc mpi datatype to partial type for velocity field
164 ! create bc mpi datatype for ghostcells update
165 call create_bc_mpi_datatype(mom(1),ndir)
166 ! convert conservative variables to primitive ones which are used during MF
167 do iigrid=1,igridstail; igrid=igrids(iigrid);
168 call phys_to_primitive(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x)
169 end do
170 ! calculate magnetofrictional velocity
171 call mf_velocity_update(dtfff)
172 ! update velocity in ghost cells
173 call getbc(tmf,0.d0,ps,mom(1),ndir)
174 ! calculate initial values of metrics
175 if(i==0) then
176 call metrics
177 call printlog_mf
178 end if
179 ! magnetofrictional loops
180 do
181 ! calculate time step based on Cmax= Alfven speed + abs(frictional speed)
182 dtfff_pe=bigdouble
183 cmax_mype=zero
184 do iigrid=1,igridstail; igrid=igrids(iigrid);
185 block=>ps(igrid)
186 ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
187 call getdtfff_courant(ps(igrid)%w,ps(igrid)%x,ixg^ll,ixm^ll,dtnew)
188 dtfff_pe=min(dtfff_pe,dtnew)
189 end do
190 call mpi_allreduce(dtfff_pe,dtfff,1,mpi_double_precision,mpi_min, &
192 call mpi_allreduce(cmax_mype,cmax_global,1,mpi_double_precision,&
193 mpi_max,icomm,ierrmpi)
194
195 ! =======
196 ! evolve
197 ! =======
198 call advectmf(1,ndim,tmf,dtfff)
199
200 if(i>=10000) then
201 mf_tvdlfeps = 0.9998d0 * mf_tvdlfeps
203 end if
204 if(i<=60000) then
205 mf_cy=1.0001d0*mf_cy
206 mf_cy = min(mf_cy_max,mf_cy)
207 end if
208 if(i>=100000) then
209 mf_cdivb=1.0000018d0*mf_cdivb
211 end if
212
213 i=i+1
214 tmf=tmf+dtfff
215 if(mod(i,10)==0) then
216 ! calculate metrics
217 call metrics
218 call printlog_mf
219 end if
220 if(mod(i,mf_ditsave)==0) then
221 it=i
223 do iigrid=1,igridstail; igrid=igrids(iigrid);
224 call phys_to_conserved(ixg^ll,ixg^ll,ps(igrid)%w,ps(igrid)%x)
225 end do
226 mf_advance=.false.
227 call saveamrfile(1)
228 call saveamrfile(2)
229 do iigrid=1,igridstail; igrid=igrids(iigrid);
230 call phys_to_primitive(ixg^ll,ixg^ll,ps(igrid)%w,ps(igrid)%x)
231 end do
232 mf_advance=.true.
233 if(mype==0) then
234 write(*,*) "itmf=",i
235 write(*,*) '<CW sin theta>:',cwsin_theta_new
236 write(*,*) '<f_i>:',f_i
237 write(*,*) '----------------------------------------------------------'
238 end if
239 end if
240 ! reconstruct AMR grid every 10 step
241 if(mod(i,ditregrid)==0 .and. refine_max_level>1) call resettree
242 if (i>=mf_it_max) then
243 if(mod(i,10)/=0) then
244 ! calculate metrics
245 call metrics
246 call printlog_mf
247 end if
248 if(mype==0) then
249 write (*,*) 'Reach maximum iteration step!'
250 write (*,*) 'The total iteration step is:', i
251 end if
252 exit
253 end if
254 enddo
255 ! point bc mpi data type back to full type for MHD
262 bcphys=.true.
263 ! set velocity back to zero and convert primitive variables back to conservative ones
264 do iigrid=1,igridstail; igrid=igrids(iigrid);
265 ps(igrid)%w(ixg^t,mom(:))=zero
266 call phys_to_conserved(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x)
267 end do
268 global_time=tmpt
269 it=tmpit
270 if (mype==0) call mpi_file_close(fhmf,ierrmpi)
271 mf_advance=.false.
272 ! restore stagger_grid value
273 stagger_grid=stagger_flag
274 if(mype==0) write(*,*) 'Magnetofriction phase took : ',mpi_wtime()-time_in,' sec'
275 contains
276
277 subroutine metrics
278
279 sum_jbb_ipe = 0.d0
280 sum_j_ipe = 0.d0
281 sum_l_ipe = 0.d0
282 f_i_ipe = 0.d0
283 volumepe=0.d0
284 dsurface=zero
285 do iigrid=1,igridstail; igrid=igrids(iigrid);
286 block=>ps(igrid)
287 dvolume(ixm^t)=block%dvolume(ixm^t)
288 do idims=1,ndim
289 hxm^ll=ixm^ll-kr(idims,^d);
290 dsurface(ixm^t)=dsurface(ixm^t)+block%surfaceC(hxm^t,idims)
291 end do
292 ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
293 call mask_inner(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x)
294 sum_jbb_ipe = sum_jbb_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
295 ps(igrid)%x,1,patchwi)
296 sum_j_ipe = sum_j_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
297 ps(igrid)%x,2,patchwi)
298 f_i_ipe=f_i_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
299 ps(igrid)%x,3,patchwi)
300 sum_l_ipe = sum_l_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
301 ps(igrid)%x,4,patchwi)
302 end do
303 call mpi_allreduce(sum_jbb_ipe,sum_jbb,1,mpi_double_precision,&
304 mpi_sum,icomm,ierrmpi)
305 call mpi_allreduce(sum_j_ipe,sum_j,1,mpi_double_precision,mpi_sum,&
307 call mpi_allreduce(f_i_ipe,f_i,1,mpi_double_precision,mpi_sum,&
309 call mpi_allreduce(sum_l_ipe,sum_l,1,mpi_double_precision,mpi_sum,&
311 call mpi_allreduce(volumepe,volume,1,mpi_double_precision,mpi_sum,&
313 ! current- and volume-weighted average of the sine of the angle
314 ! between the magnetic field B and the current density J
315 cwsin_theta_new = sum_jbb/sum_j
316 ! volume-weighted average of the absolute value of the fractional
317 ! magnetic flux change
318 f_i = f_i/volume
319 sum_j=sum_j/volume
320 sum_l=sum_l/volume
321 end subroutine metrics
322
323 subroutine mask_inner(ixI^L,ixO^L,w,x)
324
325 integer, intent(in) :: ixI^L,ixO^L
326 double precision, intent(in):: w(ixI^S,nw),x(ixI^S,1:ndim)
327 double precision :: xO^L
328 integer :: ix^D
329
330 {xomin^d = xprobmin^d + 0.05d0*(xprobmax^d-xprobmin^d)\}
331 {xomax^d = xprobmax^d - 0.05d0*(xprobmax^d-xprobmin^d)\}
332 if(slab) then
333 xomin^nd = xprobmin^nd
334 else
335 xomin1 = xprobmin1
336 end if
337
338 {do ix^db=ixomin^db,ixomax^db\}
339 if({ x(ix^dd,^d) > xomin^d .and. x(ix^dd,^d) < xomax^d | .and. }) then
340 patchwi(ix^d)=.true.
341 volumepe=volumepe+dvolume(ix^d)
342 else
343 patchwi(ix^d)=.false.
344 endif
345 {end do\}
346
347 end subroutine mask_inner
348
349 subroutine printlog_mf
350 integer :: amode, status(MPI_STATUS_SIZE)
351 character(len=800) :: filename,filehead
352 character(len=2048) :: line,datastr
353 logical, save :: logmfopened=.false.
354
355 if(mype==0) then
356 if(.not.logmfopened) then
357 ! generate filename
358 write(filename,"(a,a)") trim(base_filename), "_mflog.csv"
359
360 amode=ior(mpi_mode_create,mpi_mode_wronly)
361 amode=ior(amode,mpi_mode_append)
362 call mpi_file_open(mpi_comm_self,filename,amode,mpi_info_null,fhmf,ierrmpi)
363 logmfopened=.true.
364 filehead=" itmf, dt, <f_i>, <CW sin theta>, <Current>, <Lorenz force>"
365 call mpi_file_write(fhmf,filehead,len_trim(filehead), &
366 mpi_character,status,ierrmpi)
367 call mpi_file_write(fhmf,achar(10),1,mpi_character,status,ierrmpi)
368 end if
369 line=''
370 write(datastr,'(i6,a)') i,','
371 line=trim(line)//trim(datastr)
372 write(datastr,'(es13.6,a)') dtfff,','
373 line=trim(line)//trim(datastr)
374 write(datastr,'(es13.6,a)') f_i,','
375 line=trim(line)//trim(datastr)
376 write(datastr,'(es13.6,a)') cwsin_theta_new,','
377 line=trim(line)//trim(datastr)
378 write(datastr,'(es13.6,a)') sum_j,','
379 line=trim(line)//trim(datastr)
380 write(datastr,'(es13.6)') sum_l
381 line=trim(line)//trim(datastr)//new_line('A')
382 call mpi_file_write(fhmf,line,len_trim(line),mpi_character,status,ierrmpi)
383 end if
384
385 end subroutine printlog_mf
386
387 function integral_grid_mf(ixI^L,ixO^L,w,x,iw,patchwi)
388 use mod_geometry
389
390 integer, intent(in) :: ixI^L,ixO^L,iw
391 double precision, intent(in) :: x(ixI^S,1:ndim)
392 double precision, intent(in) :: w(ixI^S,nw+nwauxio)
393 logical, intent(in) :: patchwi(ixI^S)
394
395 double precision, dimension(ixI^S,1:ndir) :: bvec,qvec,current
396 double precision :: integral_grid_mf,tmp(ixI^S),b_mag(ixI^S)
397 integer :: ix^D,i,idirmin,idir,jdir,kdir
398
399 integral_grid_mf=0.d0
400 select case(iw)
401 case(1)
402 ! Sum(dvolume*|JxB|/|B|)
403 if(b0field) then
404 bvec(ixi^s,:)=w(ixi^s,mag(:))+block%b0(ixi^s,mag(:),0)
405 else
406 bvec(ixi^s,:)=w(ixi^s,mag(:))
407 endif
408 call get_current(w,ixi^l,ixo^l,idirmin,current)
409 ! calculate Lorentz force
410 qvec(ixo^s,1:ndir)=zero
411 do idir=1,ndir; do jdir=1,ndir; do kdir=idirmin,3
412 if(lvc(idir,jdir,kdir)/=0)then
413 tmp(ixo^s)=current(ixo^s,jdir)*bvec(ixo^s,kdir)
414 if(lvc(idir,jdir,kdir)==1)then
415 qvec(ixo^s,idir)=qvec(ixo^s,idir)+tmp(ixo^s)
416 else
417 qvec(ixo^s,idir)=qvec(ixo^s,idir)-tmp(ixo^s)
418 endif
419 endif
420 enddo; enddo; enddo
421
422 {do ix^db=ixomin^db,ixomax^db\}
423 if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+sqrt(sum(qvec(ix^d,:)**2)/&
424 sum(bvec(ix^d,:)**2))*dvolume(ix^d)
425 {end do\}
426 case(2)
427 ! Sum(dvolume*|J|)
428 call get_current(w,ixi^l,ixo^l,idirmin,current)
429 {do ix^db=ixomin^db,ixomax^db\}
430 if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+sqrt(sum(current(ix^d,:)**2))*&
431 dvolume(ix^d)
432 {end do\}
433 case(3)
434 ! f_i solenoidal property of B: (dvolume |div B|)/(dsurface |B|)
435 ! Sum(dvolume*f_i)
436 if(b0field) then
437 bvec(ixi^s,:)=w(ixi^s,mag(:))+block%b0(ixi^s,mag(:),0)
438 else
439 bvec(ixi^s,:)=w(ixi^s,mag(:))
440 endif
441 call divvector(bvec,ixi^l,ixo^l,tmp)
442 {do ix^db=ixomin^db,ixomax^db\}
443 if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+abs(tmp(ix^d))*&
444 dvolume(ix^d)**2/sqrt(sum(bvec(ix^d,:)**2))/dsurface(ix^d)
445 {end do\}
446 case(4)
447 ! Sum(|JxB|)
448 if(b0field) then
449 bvec(ixi^s,:)=w(ixi^s,mag(:))+block%b0(ixi^s,mag(:),0)
450 else
451 bvec(ixi^s,:)=w(ixi^s,mag(:))
452 endif
453 call curlvector(bvec,ixi^l,ixo^l,current,idirmin,1,ndir)
454 ! calculate Lorentz force
455 qvec(ixo^s,1:ndir)=zero
456 do idir=1,ndir; do jdir=1,ndir; do kdir=idirmin,3
457 if(lvc(idir,jdir,kdir)/=0)then
458 tmp(ixo^s)=current(ixo^s,jdir)*bvec(ixo^s,kdir)
459 if(lvc(idir,jdir,kdir)==1)then
460 qvec(ixo^s,idir)=qvec(ixo^s,idir)+tmp(ixo^s)
461 else
462 qvec(ixo^s,idir)=qvec(ixo^s,idir)-tmp(ixo^s)
463 endif
464 endif
465 enddo; enddo; enddo
466
467 {do ix^db=ixomin^db,ixomax^db\}
468 if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+sqrt(sum(qvec(ix^d,:)**2))*dvolume(ix^d)
469 {end do\}
470 end select
471 return
472 end function integral_grid_mf
473
474 end subroutine magnetofriction
475
476 subroutine mf_velocity_update(dtfff)
478
479 double precision, intent(in) :: dtfff
480 double precision :: vhatmax,vhatmax_pe,vhatmaxgrid
481 integer :: i,iigrid, igrid
482
483 vhatmax_pe=smalldouble
484 do iigrid=1,igridstail; igrid=igrids(iigrid);
485 block=>ps(igrid)
486 ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
487 call vhat(ps(igrid)%w,ps(igrid)%x,ixg^ll,ixm^ll,vhatmaxgrid)
488 vhatmax_pe=max(vhatmax_pe,vhatmaxgrid)
489 end do
490 call mpi_allreduce(vhatmax_pe,vhatmax,1,mpi_double_precision,mpi_max, &
492 do iigrid=1,igridstail; igrid=igrids(iigrid);
493 block=>ps(igrid)
494 ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
495 ! calculate frictional velocity
496 call frictional_velocity(ps(igrid)%w,ps(igrid)%x,ixg^ll,ixm^ll,vhatmax,dtfff)
497 end do
498
499 end subroutine mf_velocity_update
500
501 subroutine vhat(w,x,ixI^L,ixO^L,vhatmaxgrid)
502 ! Calculate v_hat
504
505 integer, intent(in) :: ixI^L, ixO^L
506 double precision, intent(inout) :: w(ixI^S,nw)
507 double precision, intent(in) :: x(ixI^S,1:ndim)
508 double precision, intent(out) :: vhatmaxgrid
509
510 double precision :: current(ixI^S,7-2*ndir:3),tmp(ixI^S),dxhm
511 double precision :: dxhms(ixO^S)
512 integer :: idirmin,idir,jdir,kdir
513
514 call get_current(w,ixi^l,ixo^l,idirmin,current)
515 w(ixi^s,mom(:))=0.d0
516 ! calculate Lorentz force
517 do idir=1,ndir; do jdir=1,ndir; do kdir=idirmin,3
518 if(lvc(idir,jdir,kdir)/=0)then
519 if(b0field) then
520 tmp(ixo^s)=current(ixo^s,jdir)*(w(ixo^s,mag(kdir))+block%b0(ixo^s,kdir,0))
521 else
522 tmp(ixo^s)=current(ixo^s,jdir)*w(ixo^s,mag(kdir))
523 endif
524 if(lvc(idir,jdir,kdir)==1)then
525 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))+tmp(ixo^s)
526 else
527 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-tmp(ixo^s)
528 endif
529 endif
530 enddo; enddo; enddo
531
532 ! 1/B**2
533 if(b0field) then
534 tmp(ixo^s)=1.d0/(sum((w(ixo^s,mag(:))+block%b0(ixo^s,:,0))**2,dim=ndim+1)+smalldouble)
535 else
536 tmp(ixo^s)=1.d0/(sum(w(ixo^s,mag(:))**2,dim=ndim+1)+smalldouble)
537 endif
538
539 if(slab_uniform) then
540 dxhm=dble(ndim)/(^d&1.0d0/dxlevel(^d)+)
541 do idir=1,ndir
542 w(ixo^s,mom(idir))=dxhm*w(ixo^s,mom(idir))*tmp(ixo^s)
543 end do
544 else
545 dxhms(ixo^s)=dble(ndim)/sum(1.d0/block%dx(ixo^s,:),dim=ndim+1)
546 do idir=1,ndir
547 w(ixo^s,mom(idir))=dxhms(ixo^s)*w(ixo^s,mom(idir))*tmp(ixo^s)
548 end do
549 end if
550 vhatmaxgrid=maxval(sqrt(sum(w(ixo^s,mom(:))**2,dim=ndim+1)))
551
552 end subroutine vhat
553
554 subroutine frictional_velocity(w,x,ixI^L,ixO^L,qvmax,qdt)
556
557 integer, intent(in) :: ixI^L, ixO^L
558 double precision, intent(in) :: x(ixI^S,1:ndim),qdt,qvmax
559 double precision, intent(inout) :: w(ixI^S,1:nw)
560
561 double precision :: dxhm,disbd(6),bfzone^D
562 double precision :: dxhms(ixO^S)
563 integer :: ix^D, idir
564 logical :: buffer
565
566 if(slab_uniform) then
567 dxhm=dble(ndim)/(^d&1.0d0/dxlevel(^d)+)
568 dxhm=mf_cc*mf_cy/qvmax*dxhm/qdt
569 ! dxhm=mf_cc*mf_cy/qvmax
570 w(ixo^s,mom(:))=w(ixo^s,mom(:))*dxhm
571 else
572 dxhms(ixo^s)=dble(ndim)/sum(1.d0/block%dx(ixo^s,:),dim=ndim+1)
573 dxhms(ixo^s)=mf_cc*mf_cy/qvmax*dxhms(ixo^s)/qdt
574 ! dxhm=mf_cc*mf_cy/qvmax
575 do idir=1,ndir
576 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))*dxhms(ixo^s)
577 end do
578 end if
579
580{^ifthreed
581 bfzone1=0.05d0*(xprobmax1-xprobmin1)
582 bfzone2=0.05d0*(xprobmax2-xprobmin2)
583 bfzone3=0.05d0*(xprobmax3-xprobmin3)
584 {do ix^db=ixomin^db,ixomax^db\}
585 disbd(1)=x(ix^d,1)-xprobmin1
586 disbd(2)=xprobmax1-x(ix^d,1)
587 disbd(3)=x(ix^d,2)-xprobmin2
588 disbd(4)=xprobmax2-x(ix^d,2)
589 disbd(5)=x(ix^d,3)-xprobmin1
590 disbd(6)=xprobmax3-x(ix^d,3)
591
592 if(slab) then
593 if(disbd(1)<bfzone1) then
594 w(ix^d,mom(:))=(1.d0-((bfzone1-disbd(1))/bfzone1)**2)*w(ix^d,mom(:))
595 endif
596 else
597 if(disbd(5)<bfzone3) then
598 w(ix^d,mom(:))=(1.d0-((bfzone3-disbd(5))/bfzone3)**2)*w(ix^d,mom(:))
599 endif
600 end if
601 if(disbd(2)<bfzone1) then
602 w(ix^d,mom(:))=(1.d0-((bfzone1-disbd(2))/bfzone1)**2)*w(ix^d,mom(:))
603 endif
604 if(disbd(3)<bfzone2) then
605 w(ix^d,mom(:))=(1.d0-((bfzone2-disbd(3))/bfzone2)**2)*w(ix^d,mom(:))
606 endif
607 if(disbd(4)<bfzone2) then
608 w(ix^d,mom(:))=(1.d0-((bfzone2-disbd(4))/bfzone2)**2)*w(ix^d,mom(:))
609 endif
610 if(disbd(6)<bfzone3) then
611 w(ix^d,mom(:))=(1.d0-((bfzone3-disbd(6))/bfzone3)**2)*w(ix^d,mom(:))
612 endif
613 {end do\}
614}
615
616 ! saturate mf velocity at mf_vmax
617 dxhms(ixo^s)=sqrt(sum(w(ixo^s,mom(:))**2,dim=ndim+1))/mf_vmax+1.d-12
618 dxhms(ixo^s)=dtanh(dxhms(ixo^s))/dxhms(ixo^s)
619 do idir=1,ndir
620 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))*dxhms(ixo^s)
621 end do
622 end subroutine frictional_velocity
623
624 subroutine advectmf(idim^LIM,qt,qdt)
625 ! integrate all grids by one step of its delta(global_time)
626 ! This subroutine is in VAC terminology equivalent to
627 ! `advect' (with the difference that it will `advect' all grids)
630 use mod_comm_lib, only: mpistop
631
632 integer, intent(in) :: idim^LIM
633 double precision, intent(in) :: qt, qdt
634
635 integer :: iigrid, igrid
636
637 call init_comm_fix_conserve(idim^lim,ndir)
639
640 do iigrid=1,igridstail; igrid=igrids(iigrid);
641 ps1(igrid)%w=ps(igrid)%w
642 end do
643
644 istep=0
645
646 select case (t_stepper)
647 case (onestep)
648 call advect1mf(flux_method,qdt,one,idim^lim,qt,ps1,qt,ps)
649 case (twostep)
650 ! predictor step
651 fix_conserve_at_step = .false.
652 call advect1mf(typepred1,qdt,half,idim^lim,qt,ps,qt,ps1)
653 ! corrector step
655 call advect1mf(flux_method,qdt,one,idim^lim,qt+half*qdt,ps1,qt,ps)
656 case (threestep)
657 ! three step Runge-Kutta in accordance with Gottlieb & Shu 1998
658 call advect1mf(flux_method,qdt,one,idim^lim,qt,ps,qt,ps1)
659
660 do iigrid=1,igridstail; igrid=igrids(iigrid);
661 ps2(igrid)%w(ixg^t,1:nwflux)=0.75d0*ps(igrid)%w(ixg^t,1:nwflux)+0.25d0*&
662 ps1(igrid)%w(ixg^t,1:nwflux)
663 if (nw>nwflux) ps2(igrid)%w(ixg^t,nwflux+1:nw) = &
664 ps(igrid)%w(ixg^t,nwflux+1:nw)
665 end do
666
667 call advect1mf(flux_method,qdt,0.25d0,idim^lim,qt+qdt,ps1,qt+dt*0.25d0,ps2)
668
669 do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
670 ps(igrid)%w(ixg^t,1:nwflux)=1.0d0/3.0d0*ps(igrid)%w(ixg^t,1:nwflux)+&
671 2.0d0/3.0d0*ps2(igrid)%w(ixg^t,1:nwflux)
672 end do
673
674 call advect1mf(flux_method,qdt,2.0d0/3.0d0,idim^lim,qt+qdt/2.0d0,ps2,qt+qdt/3.0d0,ps)
675 case default
676 call mpistop("unkown time_stepper in advectmf")
677 end select
678
679 end subroutine advectmf
680
681 subroutine advect1mf(method,dtin,dtfactor,idim^LIM,qtC,psa,qt,psb)
682 ! Integrate all grids by one partial step
683 ! This subroutine is equivalent to VAC's `advect1', but does
684 ! the advection for all grids
688
689 integer, intent(in) :: idim^LIM
690 type(state) :: psa(max_blocks)! Compute fluxes based on this state
691 type(state) :: psb(max_blocks) ! update on this state
692 double precision, intent(in) :: dtin,dtfactor, qtC, qt
693 integer, intent(in) :: method(nlevelshi)
694
695 double precision :: qdt
696 integer :: iigrid, igrid, level, i^D
697 logical :: setigrid
698
699 istep=istep+1
700
701 ! loop over all grids to arrive at equivalent
702 qdt=dtfactor*dtin
703 do iigrid=1,igridstail; igrid=igrids(iigrid);
704 block=>ps(igrid)
705 level=node(plevel_,igrid)
706
707 call process1_gridmf(method(level),igrid,qdt,ixg^ll,idim^lim,qtc,&
708 psa(igrid)%w,qt,psb(igrid)%w)
709 end do
710
711 ! opedit: Send flux for all grids, expects sends for all
712 ! nsend_fc(^D), set in connectivity.t.
713
714 if (fix_conserve_at_step) then
715 call recvflux(idim^lim)
716 call sendflux(idim^lim)
717 call fix_conserve(psb,idim^lim,mag(1),ndir)
718 end if
719 ! point bc mpi datatype to partial type for magnetic field
726 ! update B in ghost cells
727 call getbc(qt+qdt,qdt,psb,mag(1),ndir)
728 ! calculate magnetofrictional velocity
729 call mf_velocity_update(qdt)
730 ! point bc mpi datatype to partial type for velocity field
737 ! update magnetofrictional velocity in ghost cells
738 call getbc(qt+qdt,qdt,psb,mom(1),ndir)
739
740 end subroutine advect1mf
741
742 subroutine process1_gridmf(method,igrid,qdt,ixG^L,idim^LIM,qtC,wCT,qt,w)
743 ! This subroutine is equivalent to VAC's `advect1' for one grid
746 use mod_comm_lib, only: mpistop
747
748 integer, intent(in) :: method
749 integer, intent(in) :: igrid, ixG^L, idim^LIM
750 double precision, intent(in) :: qdt, qtC, qt
751 double precision :: wCT(ixG^S,1:nw), w(ixG^S,1:nw)
752 double precision :: dx^D, fC(ixG^S,1:ndir,1:ndim)
753 integer :: ixO^L
754
755 dx^d=rnode(rpdx^d_,igrid);
756 ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
757 fc=0.d0
758
759 ixo^l=ixg^l^lsubnghostcells;
760 select case (method)
761 case (fs_cd4)
762 !================================
763 ! 4th order central difference
764 !================================
765 call centdiff4mf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,fc,dx^d,ps(igrid)%x)
766 case (fs_tvdlf)
767 !================================
768 ! TVDLF
769 !================================
770 call tvdlfmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,fc,dx^d,ps(igrid)%x)
771 case (fs_hancock)
772 ! hancock predict (first) step for twostep tvdlf and tvdmu scheme
773 call hancockmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,dx^d,ps(igrid)%x)
774 case (fs_fd)
775 !================================
776 ! finite difference
777 !================================
778 call fdmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,fc,dx^d,ps(igrid)%x)
779 case default
780 call mpistop("unknown flux scheme in advect1_gridmf")
781 end select
782
783 if (fix_conserve_at_step) then
784 call store_flux(igrid,fc,idim^lim,ndir)
785 end if
786
787 end subroutine process1_gridmf
788
789 subroutine upwindlrmf(ixI^L,ixL^L,ixR^L,idim,w,wCT,wLC,wRC,x)
790 ! Determine the upwinded wLC(ixL) and wRC(ixR) from w.
791 ! the wCT is only used when PPM is exploited.
793 use mod_limiter
794
795 integer, intent(in) :: ixI^L, ixL^L, ixR^L, idim
796 double precision, dimension(ixI^S,1:nw) :: w, wCT
797 double precision, dimension(ixI^S,1:nw) :: wLC, wRC
798 double precision, dimension(ixI^S,1:ndim) :: x
799
800 double precision :: ldw(ixI^S), rdw(ixI^S), dwC(ixI^S)
801 integer :: jxR^L, ixC^L, jxC^L, iw
802
803 if (type_limiter(block%level) == limiter_mp5) then
804 call mp5limiter(ixi^l,ixl^l,idim,w,wlc,wrc)
805 else if (type_limiter(block%level) == limiter_ppm) then
806 call ppmlimiter(ixi^l,ixm^ll,idim,w,wct,wlc,wrc)
807 else
808 jxr^l=ixr^l+kr(idim,^d);
809 ixcmax^d=jxrmax^d; ixcmin^d=ixlmin^d-kr(idim,^d);
810 jxc^l=ixc^l+kr(idim,^d);
811
812 do iw=1,nwflux
813 if (loglimit(iw)) then
814 w(ixcmin^d:jxcmax^d,iw)=dlog10(w(ixcmin^d:jxcmax^d,iw))
815 wlc(ixl^s,iw)=dlog10(wlc(ixl^s,iw))
816 wrc(ixr^s,iw)=dlog10(wrc(ixr^s,iw))
817 end if
818
819 dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
820
821 ! limit flux from left and/or right
822 call dwlimiter2(dwc,ixi^l,ixc^l,idim,type_limiter(block%level),ldw,rdw)
823 wlc(ixl^s,iw)=wlc(ixl^s,iw)+half*ldw(ixl^s)
824 wrc(ixr^s,iw)=wrc(ixr^s,iw)-half*rdw(jxr^s)
825
826 if (loglimit(iw)) then
827 w(ixcmin^d:jxcmax^d,iw)=10.0d0**w(ixcmin^d:jxcmax^d,iw)
828 wlc(ixl^s,iw)=10.0d0**wlc(ixl^s,iw)
829 wrc(ixr^s,iw)=10.0d0**wrc(ixr^s,iw)
830 end if
831 end do
832
833 endif
834
835 end subroutine upwindlrmf
836
837 subroutine getfluxmf(w,x,ixI^L,ixO^L,idir,idim,f)
838 ! Calculate lux f_idim[idir] within ixO^L.
840
841 integer, intent(in) :: ixI^L, ixO^L, idir, idim
842 double precision, intent(in) :: w(ixI^S,nw)
843 double precision, intent(in) :: x(ixI^S,1:ndim)
844 double precision,intent(out) :: f(ixI^S)
845
846 ! compute flux of magnetic field
847 ! f_i[b_k]=v_i*b_k-v_k*b_i
848 if (idim==idir) then
849 f(ixo^s)=zero
850 else
851 f(ixo^s)=w(ixo^s,mom(idim))*w(ixo^s,mag(idir))-w(ixo^s,mag(idim))*w(ixo^s,mom(idir))
852 if (b0field) then
853 f(ixo^s)=f(ixo^s)&
854 +w(ixo^s,mom(idim))*block%B0(ixo^s,idir,idim)&
855 -w(ixo^s,mom(idir))*block%B0(ixo^s,idim,idim)
856 end if
857 end if
858
859 end subroutine getfluxmf
860
861 subroutine tvdlfmf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,wnew,fC,dx^D,x)
862 ! method=='tvdlf' --> 2nd order TVD-Lax-Friedrich scheme.
863 ! method=='tvdlf1' --> 1st order TVD-Lax-Friedrich scheme.
865 use mod_comm_lib, only: mpistop
866
867 double precision, intent(in) :: qdt, qtC, qt, dx^D
868 integer, intent(in) :: ixI^L, ixO^L, idim^LIM
869 double precision, dimension(ixI^S,1:ndim), intent(in) :: x
870 double precision, dimension(ixI^S,1:nw) :: wCT, wnew
871 double precision, dimension(ixI^S,1:ndir,1:ndim) :: fC
872
873 double precision, dimension(ixI^S,1:nw) :: wLC, wRC, wmean
874 double precision, dimension(ixI^S) :: fLC, fRC
875 double precision, dimension(ixI^S) :: cmaxC
876 double precision :: dxinv(1:ndim), inv_volume(ixO^S)
877 integer :: idims, idir, ix^L, hxO^L, ixC^L, ixCR^L, jxC^L, kxC^L, kxR^L
878
879 ! The flux calculation contracts by one in the idim direction it is applied.
880 ! The limiter contracts the same directions by one more, so expand ixO by 2.
881 ix^l=ixo^l;
882 do idims= idim^lim
883 ix^l=ix^l^ladd2*kr(idims,^d);
884 end do
885 if (ixi^l^ltix^l|.or.|.or.) &
886 call mpistop("Error in tvdlfmf: Nonconforming input limits")
887
888 ^d&dxinv(^d)=-qdt/dx^d;
889 fc=0.d0
890 do idims= idim^lim
891 b0i=idims
892
893 hxo^l=ixo^l-kr(idims,^d);
894 ! ixC is centered index in the idim direction from ixOmin-1/2 to ixOmax+1/2
895 ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
896 ! Calculate wRC=uR_{j+1/2} and wLC=uL_j+1/2
897 jxc^l=ixc^l+kr(idims,^d);
898 kxcmin^d=iximin^d; kxcmax^d=iximax^d-kr(idims,^d);
899 kxr^l=kxc^l+kr(idims,^d);
900 ixcr^l=ixc^l;
901
902 wrc(kxc^s,1:nwflux)=wct(kxr^s,1:nwflux)
903 wlc(kxc^s,1:nwflux)=wct(kxc^s,1:nwflux)
904
905 call upwindlrmf(ixi^l,ixcr^l,ixcr^l,idims,wct,wct,wlc,wrc,x)
906
907 ! For the high order Lax-Friedrich TVDLF scheme the limiter is based on
908 ! the maximum eigenvalue, it is calculated in advance.
909 ! determine mean state and store in wLC
910 wmean=0.5d0*(wlc+wrc)
911 call getcmaxfff(wmean,ixg^ll,ixc^l,idims,cmaxc)
912
913 ! Calculate fLC=f(uL_j+1/2) and fRC=f(uR_j+1/2) for each idir
914 do idir=1,ndir
915 call getfluxmf(wlc,x,ixg^ll,ixc^l,idir,idims,flc)
916 call getfluxmf(wrc,x,ixg^ll,ixc^l,idir,idims,frc)
917 ! To save memory we use fLC to store (F_L+F_R)/2=half*(fLC+fRC)
918 flc(ixc^s)=half*(flc(ixc^s)+frc(ixc^s))
919
920 ! Add TVDLF dissipation to the flux
921 if (idir==idims) then
922 flc(ixc^s)=flc(ixc^s)-mf_tvdlfeps*tvdlfeps*cmaxc(ixc^s)*half*(wrc(ixc^s,mag(idir))-wlc(ixc^s,mag(idir)))
923 end if
924 if (slab_uniform) then
925 fc(ixc^s,idir,idims)=flc(ixc^s)
926 else
927 fc(ixc^s,idir,idims)=block%surfaceC(ixc^s,idims)*flc(ixc^s)
928 end if
929
930 end do ! Next idir
931 end do ! Next idims
932 b0i=0
933
934 !Now update the state:
935 do idims= idim^lim
936 hxo^l=ixo^l-kr(idims,^d);
937 ! Multiply the fluxes by -dt/dx since Flux fixing expects this
938 if (slab_uniform) then
939 fc(ixi^s,:,idims)=dxinv(idims)*fc(ixi^s,:,idims)
940 wnew(ixo^s,mag(:))=wnew(ixo^s,mag(:)) &
941 + (fc(ixo^s,:,idims)-fc(hxo^s,:,idims))
942 else
943 inv_volume = 1.0d0/block%dvolume(ixo^s)
944 fc(ixi^s,:,idims)=-qdt*fc(ixi^s,:,idims)
945
946 do idir = 1, ndir
947 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir)) + (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims)) * &
948 inv_volume
949 end do
950 end if
951
952 end do ! Next idims
953
954 if (.not.slab) call addgeometrymf(qdt,ixi^l,ixo^l,wct,wnew,x)
955 call divbclean(qdt,ixi^l,ixo^l,wct,wnew,x)
956
957 end subroutine tvdlfmf
958
959 subroutine hancockmf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,wnew,dx^D,x)
960 ! The non-conservative Hancock predictor for TVDLFmf
961 ! on entry:
962 ! input available on ixI^L=ixG^L asks for output on ixO^L=ixG^L^LSUBnghostcells
963 ! one entry: (predictor): wCT -- w_n wnew -- w_n qdt=dt/2
964 ! on exit : (predictor): wCT -- w_n wnew -- w_n+1/2
966 use mod_comm_lib, only: mpistop
967
968 integer, intent(in) :: ixI^L, ixO^L, idim^LIM
969 double precision, intent(in) :: qdt, qtC, qt, dx^D, x(ixI^S,1:ndim)
970 double precision, intent(inout) :: wCT(ixI^S,1:nw), wnew(ixI^S,1:nw)
971
972 double precision, dimension(ixI^S,1:nw) :: wLC, wRC
973 double precision, dimension(ixI^S) :: fLC, fRC
974 double precision :: dxinv(1:ndim)
975 integer :: idims, idir, ix^L, hxO^L, ixtest^L
976
977 ! Expand limits in each idims direction in which fluxes are added
978 ix^l=ixo^l;
979 do idims= idim^lim
980 ix^l=ix^l^laddkr(idims,^d);
981 end do
982 if (ixi^l^ltix^l|.or.|.or.) &
983 call mpistop("Error in Hancockmf: Nonconforming input limits")
984
985 ^d&dxinv(^d)=-qdt/dx^d;
986 do idims= idim^lim
987 b0i=idims
988 ! Calculate w_j+g_j/2 and w_j-g_j/2
989 ! First copy all variables, then upwind wLC and wRC.
990 ! wLC is to the left of ixO, wRC is to the right of wCT.
991 hxo^l=ixo^l-kr(idims,^d);
992
993 wrc(hxo^s,1:nwflux)=wct(ixo^s,1:nwflux)
994 wlc(ixo^s,1:nwflux)=wct(ixo^s,1:nwflux)
995
996 call upwindlrmf(ixi^l,ixo^l,hxo^l,idims,wct,wct,wlc,wrc,x)
997
998 ! Advect mag(idir)
999 do idir=1,ndir
1000 ! Calculate the fLC and fRC fluxes
1001 call getfluxmf(wrc,x,ixi^l,hxo^l,idir,idims,frc)
1002 call getfluxmf(wlc,x,ixi^l,ixo^l,idir,idims,flc)
1003
1004 if (slab_uniform) then
1005 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))+dxinv(idims)* &
1006 (flc(ixo^s)-frc(hxo^s))
1007 else
1008 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))-qdt/block%dvolume(ixo^s) &
1009 *(block%surfaceC(ixo^s,idims)*flc(ixo^s) &
1010 -block%surfaceC(hxo^s,idims)*frc(hxo^s))
1011 end if
1012 end do
1013 end do ! next idims
1014 b0i=0
1015
1016 if (.not.slab) call addgeometrymf(qdt,ixi^l,ixo^l,wct,wnew,x)
1017
1018 end subroutine hancockmf
1019
1020 subroutine fdmf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,wnew,fC,dx^D,x)
1022 double precision, intent(in) :: qdt, qtC, qt, dx^D
1023 integer, intent(in) :: ixI^L, ixO^L, idim^LIM
1024 double precision, dimension(ixI^S,1:ndim), intent(in) :: x
1025
1026 double precision, dimension(ixI^S,1:nw), intent(inout) :: wCT, wnew
1027 double precision, dimension(ixI^S,1:ndir,1:ndim), intent(out) :: fC
1028
1029 double precision, dimension(ixI^S) :: fCT
1030 double precision, dimension(ixI^S,1:nw) :: fm, fp, fmR, fpL
1031 double precision, dimension(ixI^S) :: v
1032 double precision :: dxinv(1:ndim)
1033 integer :: idims, idir, ixC^L, ix^L, hxO^L, ixCR^L
1034
1035 ^d&dxinv(^d)=-qdt/dx^d;
1036 do idims= idim^lim
1037
1038 ! Get fluxes for the whole grid (mesh+nghostcells)
1039 {^d& ixcmin^d = ixomin^d - nghostcells * kr(idims,^d)\}
1040 {^d& ixcmax^d = ixomax^d + nghostcells * kr(idims,^d)\}
1041
1042 hxo^l=ixo^l-kr(idims,^d);
1043 ! ix is centered index in the idim direction from ixOmin-1/2 to ixOmax+1/2
1044 ixmax^d=ixomax^d; ixmin^d=hxomin^d;
1045 ixcr^l=ixc^l;
1046
1047 do idir=1,ndir
1048 call getfluxmf(wct,x,ixg^ll,ixcr^l,idir,idims,fct)
1049 ! Lax-Friedrich splitting:
1050 fp(ixcr^s,mag(idir)) = half * (fct(ixcr^s) + mf_tvdlfeps * tvdlfeps * cmax_global * wct(ixcr^s,mag(idir)))
1051 fm(ixcr^s,mag(idir)) = half * (fct(ixcr^s) - mf_tvdlfeps * tvdlfeps * cmax_global * wct(ixcr^s,mag(idir)))
1052 end do ! iw loop
1053
1054 ! now do the reconstruction of fp and fm:
1055 call reconstructlmf(ixi^l,ix^l,idims,fp,fpl)
1056 call reconstructrmf(ixi^l,ix^l,idims,fm,fmr)
1057
1058 do idir=1,ndir
1059 if (slab_uniform) then
1060 fc(ix^s,idir,idims) = dxinv(idims) * (fpl(ix^s,mag(idir)) + fmr(ix^s,mag(idir)))
1061 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))+ &
1062 (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))
1063 else
1064 fc(ix^s,idir,idims)=-qdt*block%surfaceC(ix^s,idims) * (fpl(ix^s,mag(idir)) + fmr(ix^s,mag(idir)))
1065 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))+ &
1066 (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))/block%dvolume(ixo^s)
1067 end if
1068 end do ! iw loop
1069
1070 end do !idims loop
1071
1072 if (.not.slab) call addgeometrymf(qdt,ixi^l,ixo^l,wct,wnew,x)
1073 call divbclean(qdt,ixi^l,ixo^l,wct,wnew,x)
1074
1075 end subroutine fdmf
1076
1077 subroutine reconstructlmf(ixI^L,iL^L,idims,w,wLC)
1079 use mod_limiter
1080
1081 integer, intent(in) :: ixI^L, iL^L, idims
1082 double precision, intent(in) :: w(ixI^S,1:nw)
1083
1084 double precision, intent(out) :: wLC(ixI^S,1:nw)
1085
1086 double precision :: ldw(ixI^S), dwC(ixI^S)
1087 integer :: jxR^L, ixC^L, jxC^L, kxC^L, iw
1088
1089 select case (type_limiter(block%level))
1090 case (limiter_mp5)
1091 call mp5limiterl(ixi^l,il^l,idims,w,wlc)
1092 case (limiter_weno5)
1093 call weno5limiterl(ixi^l,il^l,idims,w,wlc,1)
1094 case (limiter_wenoz5)
1095 call weno5limiterl(ixi^l,il^l,idims,w,wlc,2)
1096 case default
1097
1098 kxcmin^d=iximin^d; kxcmax^d=iximax^d-kr(idims,^d);
1099
1100 wlc(kxc^s,1:nwflux) = w(kxc^s,1:nwflux)
1101
1102 jxr^l=il^l+kr(idims,^d);
1103
1104 ixcmax^d=jxrmax^d; ixcmin^d=ilmin^d-kr(idims,^d);
1105 jxc^l=ixc^l+kr(idims,^d);
1106
1107 do iw=1,nwflux
1108 dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
1109
1110 call dwlimiter2(dwc,ixi^l,ixc^l,idims,type_limiter(block%level),ldw=ldw)
1111
1112 wlc(il^s,iw)=wlc(il^s,iw)+half*ldw(il^s)
1113 end do
1114 end select
1115
1116 end subroutine reconstructlmf
1117
1118 subroutine reconstructrmf(ixI^L,iL^L,idims,w,wRC)
1120 use mod_limiter
1121
1122 integer, intent(in) :: ixI^L, iL^L, idims
1123 double precision, intent(in) :: w(ixI^S,1:nw)
1124
1125 double precision, intent(out) :: wRC(ixI^S,1:nw)
1126
1127 double precision :: rdw(ixI^S), dwC(ixI^S)
1128 integer :: jxR^L, ixC^L, jxC^L, kxC^L, kxR^L, iw
1129
1130 select case (type_limiter(block%level))
1131 case (limiter_mp5)
1132 call mp5limiterr(ixi^l,il^l,idims,w,wrc)
1133 case (limiter_weno5)
1134 call weno5limiterr(ixi^l,il^l,idims,w,wrc,1)
1135 case (limiter_wenoz5)
1136 call weno5limiterr(ixi^l,il^l,idims,w,wrc,2)
1137 case default
1138
1139 kxcmin^d=iximin^d; kxcmax^d=iximax^d-kr(idims,^d);
1140 kxr^l=kxc^l+kr(idims,^d);
1141
1142 wrc(kxc^s,1:nwflux)=w(kxr^s,1:nwflux)
1143
1144 jxr^l=il^l+kr(idims,^d);
1145 ixcmax^d=jxrmax^d; ixcmin^d=ilmin^d-kr(idims,^d);
1146 jxc^l=ixc^l+kr(idims,^d);
1147
1148 do iw=1,nwflux
1149 dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
1150 call dwlimiter2(dwc,ixi^l,ixc^l,idims,type_limiter(block%level),rdw=rdw)
1151
1152 wrc(il^s,iw)=wrc(il^s,iw)-half*rdw(jxr^s)
1153 end do
1154 end select
1155
1156 end subroutine reconstructrmf
1157
1158 subroutine centdiff4mf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,w,fC,dx^D,x)
1159 ! Advance the flow variables from global_time to global_time+qdt within ixO^L by
1160 ! fourth order centered differencing in space
1161 ! for the dw/dt+dF_i(w)/dx_i=S type equation.
1162 ! wCT contains the time centered variables at time qtC for flux and source.
1163 ! w is the old value at qt on input and the new value at qt+qdt on output.
1165 use mod_comm_lib, only: mpistop
1166
1167 integer, intent(in) :: ixI^L, ixO^L, idim^LIM
1168 double precision, intent(in) :: qdt, qtC, qt, dx^D
1169 double precision :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
1170 double precision, intent(in) :: x(ixI^S,1:ndim)
1171 double precision :: fC(ixI^S,1:ndir,1:ndim)
1172
1173 double precision :: v(ixI^S,ndim), f(ixI^S)
1174 double precision, dimension(ixI^S,1:nw) :: wLC, wRC
1175 double precision, dimension(ixI^S) :: vLC, vRC,cmaxLC,cmaxRC
1176 double precision :: dxinv(1:ndim)
1177 integer :: idims, idir, idirmin,ix^D
1178 integer :: ix^L, hxO^L, ixC^L, jxC^L, hxC^L, kxC^L, kkxC^L, kkxR^L
1179
1180 ! two extra layers are needed in each direction for which fluxes are added.
1181 ix^l=ixo^l;
1182 do idims= idim^lim
1183 ix^l=ix^l^ladd2*kr(idims,^d);
1184 end do
1185
1186 if (ixi^l^ltix^l|.or.|.or.) then
1187 call mpistop("Error in evolve_CentDiff4: Non-conforming input limits")
1188 end if
1189 ^d&dxinv(^d)=-qdt/dx^d;
1190
1191 ! Add fluxes to w
1192 do idims= idim^lim
1193 b0i=idims
1194 ix^l=ixo^l^ladd2*kr(idims,^d);
1195 hxo^l=ixo^l-kr(idims,^d);
1196
1197 ixcmin^d=hxomin^d; ixcmax^d=ixomax^d;
1198 hxc^l=ixc^l-kr(idims,^d);
1199 jxc^l=ixc^l+kr(idims,^d);
1200 kxc^l=ixc^l+2*kr(idims,^d);
1201
1202 kkxcmin^d=iximin^d; kkxcmax^d=iximax^d-kr(idims,^d);
1203 kkxr^l=kkxc^l+kr(idims,^d);
1204 wrc(kkxc^s,1:nwflux)=wct(kkxr^s,1:nwflux)
1205 wlc(kkxc^s,1:nwflux)=wct(kkxc^s,1:nwflux)
1206
1207 call upwindlrmf(ixi^l,ixc^l,ixc^l,idims,wct,wct,wlc,wrc,x)
1208
1209 ! Calculate velocities from upwinded values
1210 call getcmaxfff(wlc,ixg^ll,ixc^l,idims,cmaxlc)
1211 call getcmaxfff(wrc,ixg^ll,ixc^l,idims,cmaxrc)
1212 ! now take the maximum of left and right states
1213 vlc(ixc^s)=max(cmaxrc(ixc^s),cmaxlc(ixc^s))
1214
1215 do idir=1,ndir
1216 ! Get non-transported flux
1217 call getfluxmf(wct,x,ixi^l,ix^l,idir,idims,f)
1218 ! Center flux to interface
1219 ! f_i+1/2= (-f_(i+2) +7 f_(i+1) + 7 f_i - f_(i-1))/12
1220 fc(ixc^s,idir,idims)=(-f(kxc^s)+7.0d0*(f(jxc^s)+f(ixc^s))-f(hxc^s))/12.0d0
1221 ! add rempel dissipative flux, only second order version for now
1222 ! one could gradually reduce the dissipative flux to improve solutions
1223 ! for computing steady states (Keppens et al. 2012, JCP)
1224 fc(ixc^s,idir,idims)=fc(ixc^s,idir,idims)-mf_tvdlfeps*tvdlfeps*half*vlc(ixc^s) &
1225 *(wrc(ixc^s,mag(idir))-wlc(ixc^s,mag(idir)))
1226
1227 if (slab_uniform) then
1228 fc(ixc^s,idir,idims)=dxinv(idims)*fc(ixc^s,idir,idims)
1229 ! result: f_(i+1/2)-f_(i-1/2) = [-f_(i+2)+8(f_(i+1)+f_(i-1))-f_(i-2)]/12
1230 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+(fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))
1231 else
1232 fc(ixc^s,idir,idims)=-qdt*block%surfaceC(ixc^s,idims)*fc(ixc^s,idir,idims)
1233 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+ &
1234 (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))/block%dvolume(ixo^s)
1235 end if
1236 end do !next idir
1237 end do !next idims
1238 b0i=0
1239
1240 if (.not.slab) call addgeometrymf(qdt,ixi^l,ixo^l,wct,w,x)
1241 call divbclean(qdt,ixi^l,ixo^l,wct,w,x)
1242
1243 end subroutine centdiff4mf
1244
1245 subroutine getdtfff_courant(w,x,ixI^L,ixO^L,dtnew)
1246 ! compute CFL limited dt (for variable time stepping)
1248
1249 integer, intent(in) :: ixI^L, ixO^L
1250 double precision, intent(in) :: x(ixI^S,1:ndim)
1251 double precision, intent(inout) :: w(ixI^S,1:nw), dtnew
1252
1253 double precision :: courantmax, dxinv(1:ndim)
1254 double precision :: cmax(ixI^S),tmp(ixI^S),alfven(ixI^S)
1255 integer :: idims
1256
1257 dtnew=bigdouble
1258 courantmax=0.d0
1259 ^d&dxinv(^d)=1.d0/dxlevel(^d);
1260
1261 do idims=1,ndim
1262 call getcmaxfff(w,ixi^l,ixo^l,idims,cmax)
1263 cmax_mype = max(cmax_mype,maxval(cmax(ixo^s)))
1264 if (.not.slab_uniform) then
1265 tmp(ixo^s)=cmax(ixo^s)/block%dx(ixo^s,idims)
1266 courantmax=max(courantmax,maxval(tmp(ixo^s)))
1267 else
1268 tmp(ixo^s)=cmax(ixo^s)*dxinv(idims)
1269 courantmax=max(courantmax,maxval(tmp(ixo^s)))
1270 end if
1271 end do
1272 ! courantmax='max( c/dx)'
1273 if (courantmax>smalldouble) dtnew=min(dtnew,mf_cc/courantmax)
1274
1275 end subroutine getdtfff_courant
1276
1277 subroutine getcmaxfff(w,ixI^L,ixO^L,idims,cmax)
1279
1280 logical :: new_cmax,needcmin
1281 integer, intent(in) :: ixI^L, ixO^L, idims
1282 double precision, intent(in) :: w(ixI^S,1:nw)
1283 double precision, intent(out) :: cmax(ixI^S)
1284
1285 ! calculate alfven speed
1286 if(b0field) then
1287 cmax(ixo^s)=sqrt(sum((w(ixo^s,mag(:))+block%b0(ixo^s,:,0))**2,dim=ndim+1)/w(ixo^s,rho_))
1288 else
1289 cmax(ixo^s)=sqrt(sum(w(ixo^s,mag(:))**2,dim=ndim+1)/w(ixo^s,rho_))
1290 endif
1291 cmax(ixo^s)=cmax(ixo^s)+abs(w(ixo^s,mom(idims)))
1292
1293 end subroutine getcmaxfff
1294
1295 !> Clean divergence of magnetic field by Janhunen's and Linde's source terms
1296 subroutine divbclean(qdt,ixI^L,ixO^L,wCT,w,x)
1298 use mod_geometry
1299
1300 integer, intent(in) :: ixI^L, ixO^L
1301 double precision, intent(in) :: x(ixI^S,1:ndim),wCT(ixI^S,1:nw),qdt
1302 double precision, intent(inout) :: w(ixI^S,1:nw)
1303 double precision :: divb(ixI^S),graddivb(ixI^S),bdivb(ixI^S,1:ndir)
1304 integer :: idims, ix^L, ixp^L, i^D, iside
1305
1306 ! Calculate div B
1307 ix^l=ixo^l^ladd1;
1308 call get_divb(wct,ixi^l,ix^l,divb)
1309
1310 ixp^l=ixo^l;
1311
1312 ! Add Linde's diffusive terms
1313 do idims=1,ndim
1314 ! Calculate grad_idim(divb)
1315 call gradient(divb,ixi^l,ixp^l,idims,graddivb)
1316
1317 ! Multiply by Linde's eta*dt = divbdiff*(c_max*dx)*dt = divbdiff*dx**2
1318 if (slab_uniform) then
1319 graddivb(ixp^s)=graddivb(ixp^s)*mf_cdivb/(^d&1.0d0/dxlevel(^d)**2+)
1320 else
1321 graddivb(ixp^s)=graddivb(ixp^s)*mf_cdivb &
1322 /(^d&1.0d0/block%dx(ixp^s,^d)**2+)
1323 end if
1324 ! B_idim += eta*grad_idim(divb)
1325 ! with Janhunen's term
1326 !w(ixp^S,mag(idims))=w(ixp^S,mag(idims))+&
1327 ! graddivb(ixp^S)-qdt*w(ixp^S,mom(idims))*divb(ixp^S)
1328 ! without Janjunen's term
1329 w(ixp^s,mag(idims))=w(ixp^s,mag(idims))+&
1330 graddivb(ixp^s)
1331 end do
1332
1333 end subroutine divbclean
1334
1335 subroutine addgeometrymf(qdt,ixI^L,ixO^L,wCT,w,x)
1336 ! Add geometrical source terms to w
1338 use mod_geometry
1339
1340 integer, intent(in) :: ixI^L, ixO^L
1341 double precision, intent(in) :: qdt, x(ixI^S,1:ndim)
1342 double precision, intent(inout) :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
1343 !.. local ..
1344 double precision :: tmp(ixI^S)
1345 integer :: iw
1346 integer :: mr_,mphi_ ! Polar var. names
1347 integer :: br_,bphi_
1348
1349 mr_=mom(1); mphi_=mom(1)-1+phi_ ! Polar var. names
1350 br_=mag(1); bphi_=mag(1)-1+phi_
1351
1352 select case (coordinate)
1353 case (cylindrical)
1354 if(phi_>0) then
1355 ! s[Bphi]=(Bphi*vr-Br*vphi)/radius
1356 tmp(ixo^s)=(wct(ixo^s,bphi_)*wct(ixo^s,mom(1)) &
1357 -wct(ixo^s,br_)*wct(ixo^s,mom(3)))
1358 w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt*tmp(ixo^s)/x(ixo^s,1)
1359 end if
1360 case (spherical)
1361 {^nooned
1362 ! s[b2]=(vr*Btheta-vtheta*Br)/r
1363 ! + cot(theta)*psi/r
1364 tmp(ixo^s)= wct(ixo^s,mom(1))*wct(ixo^s,mag(2)) &
1365 -wct(ixo^s,mom(2))*wct(ixo^s,mag(1))
1366 if (b0field) then
1367 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,mom(1))*block%b0(ixo^s,2,0) &
1368 -wct(ixo^s,mom(2))*block%b0(ixo^s,1,0)
1369 end if
1370 ! Divide by radius and add to w
1371 w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
1372 }
1373 if(ndir==3) then
1374 ! s[b3]=(vr*Bphi-vphi*Br)/r
1375 ! -cot(theta)*(vphi*Btheta-vtheta*Bphi)/r
1376 tmp(ixo^s)=wct(ixo^s,mom(1))*wct(ixo^s,mag(3)) &
1377 -wct(ixo^s,mom(3))*wct(ixo^s,mag(1)){^nooned &
1378 -(wct(ixo^s,mom(3))*wct(ixo^s,mag(2)) &
1379 -wct(ixo^s,mom(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
1380 /dsin(x(ixo^s,2)) }
1381 if (b0field) then
1382 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,mom(1))*block%b0(ixo^s,3,0) &
1383 -wct(ixo^s,mom(3))*block%b0(ixo^s,1,0){^nooned &
1384 -(wct(ixo^s,mom(3))*block%b0(ixo^s,2,0) &
1385 -wct(ixo^s,mom(2))*block%b0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
1386 /dsin(x(ixo^s,2)) }
1387 end if
1388 ! Divide by radius and add to w
1389 w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
1390 end if
1391 end select
1392
1393 end subroutine addgeometrymf
1394
1395 !> Calculate idirmin and the idirmin:3 components of the common current array
1396 !> make sure that dxlevel(^D) is set correctly.
1397 subroutine get_current(w,ixI^L,ixO^L,idirmin,current)
1399 use mod_geometry
1400
1401 integer :: idirmin0
1402 integer :: ixO^L, idirmin, ixI^L
1403 double precision :: w(ixI^S,1:nw)
1404
1405 ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
1406 double precision :: current(ixI^S,7-2*ndir:3),bvec(ixI^S,1:ndir)
1407 integer :: idir
1408
1409 idirmin0 = 7-2*ndir
1410
1411 bvec(ixi^s,1:ndir)=w(ixi^s,mag(1:ndir))
1412
1413 call curlvector(bvec,ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
1414
1415 if(b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
1416 block%J0(ixo^s,idirmin0:3)
1417
1418 end subroutine get_current
1419
1420 !> Calculate div B within ixO
1421 subroutine get_divb(w,ixI^L,ixO^L,divb)
1422 use mod_geometry
1424
1425 integer, intent(in) :: ixI^L, ixO^L
1426 double precision, intent(in) :: w(ixI^S,1:nw)
1427 double precision :: divb(ixI^S)
1428
1429 double precision :: bvec(ixI^S,1:ndir)
1430
1431 bvec(ixi^s,:)=w(ixi^s,mag(:))
1432
1433 select case(typediv)
1434 case("central")
1435 call divvector(bvec,ixi^l,ixo^l,divb)
1436 case("limited")
1437 call divvectors(bvec,ixi^l,ixo^l,divb)
1438 end select
1439 end subroutine get_divb
1440
1441end module mod_magnetofriction
subroutine metrics
subroutine, public resettree
reset AMR and (de)allocate boundary flux storage at level changes
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for flux conservation near refinement boundaries.
subroutine, public init_comm_fix_conserve(idimlim, nwfluxin)
subroutine, public recvflux(idimlim)
subroutine, public sendflux(idimlim)
subroutine, public store_flux(igrid, fc, idimlim, nwfluxin)
subroutine, public fix_conserve(psb, idimlim, nw0, nwfluxin)
Module with geometry-related routines (e.g., divergence, curl)
Definition mod_geometry.t:2
subroutine divvector(qvec, ixil, ixol, divq, fourthorder, sixthorder)
Calculate divergence of a vector qvec within ixL.
integer coordinate
Definition mod_geometry.t:7
integer, parameter spherical
subroutine gradient(q, ixil, ixol, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
integer, parameter cylindrical
subroutine curlvector(qvec, ixil, ixol, curlvec, idirmin, idirmin0, ndir0, fourthorder)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
subroutine divvectors(qvec, ixil, ixol, divq)
Calculate divergence of a vector qvec within ixL using limited extrapolation to cell edges.
update ghost cells of all blocks including physical boundaries
integer, dimension(0:3^d &), target type_recv_p_p1
integer, dimension( :^d &), pointer type_recv_r
integer, dimension(-1:1^d &), target type_send_srl_p1
integer, dimension(-1:1^d &), target type_send_r_p2
subroutine getbc(time, qdt, psb, nwstart, nwbc)
do update ghost cells of all blocks including physical boundaries
integer, dimension(0:3^d &), target type_send_p_f
integer, dimension( :^d &), pointer type_send_p
integer, dimension(0:3^d &), target type_recv_p_p2
integer, dimension( :^d &), pointer type_send_srl
integer, dimension(-1:1^d &), target type_send_r_p1
integer, dimension(-1:1^d &), target type_send_srl_p2
subroutine create_bc_mpi_datatype(nwstart, nwbc)
integer, dimension(0:3^d &), target type_send_p_p2
integer, dimension(0:3^d &), target type_recv_r_f
integer, dimension(-1:1^d &), target type_recv_srl_f
integer, dimension(-1:1^d &), target type_send_r_f
integer, dimension(-1:1^d &), target type_send_srl_f
integer, dimension(0:3^d &), target type_send_p_p1
integer, dimension( :^d &), pointer type_send_r
integer, dimension(0:3^d &), target type_recv_r_p1
integer, dimension(0:3^d &), target type_recv_r_p2
integer, dimension(-1:1^d &), target type_recv_srl_p2
integer, dimension( :^d &), pointer type_recv_p
integer, dimension(-1:1^d &), target type_recv_srl_p1
integer, dimension(0:3^d &), target type_recv_p_f
integer, dimension( :^d &), pointer type_recv_srl
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
integer, dimension(:), allocatable typepred1
The spatial discretization for the predictor step when using a two step PC method.
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
integer, parameter unitpar
file handle for IO
double precision global_time
The global simulation time.
integer istep
Index of the sub-step in a multi-step time integrator.
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer snapshotini
Resume from the snapshot with this index.
integer it
Number of time steps taken.
logical, dimension(:), allocatable loglimit
integer ditregrid
Reconstruct the AMR grid once every ditregrid iteration(s)
integer, parameter ndim
Number of spatial dimensions for grid variables.
logical stagger_grid
True for using stagger grid.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer icomm
The MPI communicator.
integer b0i
background magnetic field location indicator
integer mype
The rank of the current MPI task.
character(len=std_len) typediv
double precision, dimension(:), allocatable, parameter d
double precision dt
global time step
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
integer, dimension(:), allocatable flux_method
Which flux scheme of spatial discretization to use (per grid level)
logical slab
Cartesian geometry or not.
double precision unit_velocity
Physical scaling factor for velocity.
logical prolongprimitive
prolongate primitive variables in level-jump ghost cells
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
integer, dimension(:), allocatable type_limiter
Type of slope limiter used for reconstructing variables on cell edges.
integer nghostcells
Number of ghost cells surrounding a grid.
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
integer t_stepper
time stepper type
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
integer refine_max_level
Maximal number of AMR levels.
integer, parameter fs_hancock
integer, dimension(:,:), allocatable node
Module for reading input and writing output.
subroutine saveamrfile(ifile)
Module with slope/flux limiters.
Definition mod_limiter.t:2
integer, parameter limiter_ppm
Definition mod_limiter.t:27
double precision d
Definition mod_limiter.t:14
integer, parameter limiter_weno5
Definition mod_limiter.t:31
integer, parameter limiter_wenoz5
Definition mod_limiter.t:33
subroutine dwlimiter2(dwc, ixil, ixcl, idims, typelim, ldw, rdw, a2max)
Limit the centered dwC differences within ixC for iw in direction idim. The limiter is chosen accordi...
integer, parameter limiter_mp5
Definition mod_limiter.t:28
module mod_magnetofriction.t Purpose: use magnetofrictional method to relax 3D magnetic field to forc...
subroutine fdmf(qdt, ixil, ixol, idimlim, qtc, wct, qt, wnew, fc, dxd, x)
subroutine getdtfff_courant(w, x, ixil, ixol, dtnew)
subroutine hancockmf(qdt, ixil, ixol, idimlim, qtc, wct, qt, wnew, dxd, x)
subroutine getfluxmf(w, x, ixil, ixol, idir, idim, f)
subroutine reconstructrmf(ixil, ill, idims, w, wrc)
subroutine getcmaxfff(w, ixil, ixol, idims, cmax)
subroutine divbclean(qdt, ixil, ixol, wct, w, x)
Clean divergence of magnetic field by Janhunen's and Linde's source terms.
subroutine centdiff4mf(qdt, ixil, ixol, idimlim, qtc, wct, qt, w, fc, dxd, x)
double precision mf_tvdlfeps
TVDLF dissipation coefficient controls the dissipation term.
subroutine magnetofriction_init()
Initialize the module.
double precision, public mf_vmax
maximal limit of magnetofrictional velocity in cm s^-1 (Pomoell 2019 A&A)
subroutine vhat(w, x, ixil, ixol, vhatmaxgrid)
double precision tmf
time in magnetofriction process
subroutine get_current(w, ixil, ixol, idirmin, current)
Calculate idirmin and the idirmin:3 components of the common current array make sure that dxlevel(^D)...
double precision mf_cy_max
subroutine addgeometrymf(qdt, ixil, ixol, wct, w, x)
subroutine frictional_velocity(w, x, ixil, ixol, qvmax, qdt)
double precision mf_cdivb
divb cleaning coefficient controls diffusion speed of divb
subroutine upwindlrmf(ixil, ixll, ixrl, idim, w, wct, wlc, wrc, x)
subroutine mf_velocity_update(dtfff)
double precision cmax_mype
maximal speed for fd scheme
subroutine process1_gridmf(method, igrid, qdt, ixgl, idimlim, qtc, wct, qt, w)
subroutine advect1mf(method, dtin, dtfactor, idimlim, qtc, psa, qt, psb)
double precision mf_cy
frictional velocity coefficient
subroutine advectmf(idimlim, qt, qdt)
double precision mf_tvdlfeps_min
subroutine tvdlfmf(qdt, ixil, ixol, idimlim, qtc, wct, qt, wnew, fc, dxd, x)
double precision mf_cdivb_max
double precision mf_cc
stability coefficient controls numerical stability
subroutine mf_params_read(files)
Read this module"s parameters from a file.
double precision cmax_global
maximal speed for fd scheme
subroutine reconstructlmf(ixil, ill, idims, w, wlc)
subroutine get_divb(w, ixil, ixol, divb)
Calculate div B within ixO.
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition mod_physics.t:4
procedure(sub_convert), pointer phys_to_primitive
Definition mod_physics.t:55
procedure(sub_convert), pointer phys_to_conserved
Definition mod_physics.t:54
Module with all the methods that users can customize in AMRVAC.
procedure(p_no_args), pointer usr_before_main_loop