MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
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 
61 contains
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)
74 111 close(unitpar)
75  end do
76 
77  end subroutine mf_params_read
78 
79  !> Initialize the module
80  subroutine magnetofriction_init()
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 
113  subroutine magnetofriction
115  use mod_physics
117  use mod_input_output
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, &
191  icomm,ierrmpi)
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,&
306  icomm,ierrmpi)
307  call mpi_allreduce(f_i_ipe,f_i,1,mpi_double_precision,mpi_sum,&
308  icomm,ierrmpi)
309  call mpi_allreduce(sum_l_ipe,sum_l,1,mpi_double_precision,mpi_sum,&
310  icomm,ierrmpi)
311  call mpi_allreduce(volumepe,volume,1,mpi_double_precision,mpi_sum,&
312  icomm,ierrmpi)
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, &
491  icomm,ierrmpi)
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)
629  use mod_fix_conserve
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
687  use mod_fix_conserve
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
745  use mod_fix_conserve
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 
1441 end module mod_magnetofriction
subroutine metrics
subroutine mask_inner(ixIL, ixOL, w, x)
subroutine printlog_mf
double precision function integral_grid_mf(ixIL, ixOL, w, x, iw, patchwi)
subroutine, public resettree
reset AMR and (de)allocate boundary flux storage at level changes
Definition: mod_amr_grid.t:47
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
Module for flux conservation near refinement boundaries.
subroutine, public recvflux(idimLIM)
subroutine, public sendflux(idimLIM)
subroutine, public fix_conserve(psb, idimLIM, nw0, nwfluxin)
subroutine, public init_comm_fix_conserve(idimLIM, nwfluxin)
subroutine, public store_flux(igrid, fC, idimLIM, nwfluxin)
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
subroutine divvectors(qvec, ixIL, ixOL, divq)
Calculate divergence of a vector qvec within ixL using limited extrapolation to cell edges.
Definition: mod_geometry.t:609
integer coordinate
Definition: mod_geometry.t:7
integer, parameter spherical
Definition: mod_geometry.t:11
integer, parameter cylindrical
Definition: mod_geometry.t:10
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
Definition: mod_geometry.t:321
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...
Definition: mod_geometry.t:664
subroutine divvector(qvec, ixIL, ixOL, divq, fourthorder, sixthorder)
Calculate divergence of a vector qvec within ixL.
Definition: mod_geometry.t:516
update ghost cells of all blocks including physical boundaries
integer, dimension(-1:2^d &,-1:1^d &), target type_recv_srl_f
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_p_p1
integer, dimension(-1:2^d &,-1:1^d &), target type_recv_srl_p2
integer, dimension(-1:1^d &,-1:1^d &), target type_send_r_p2
integer, dimension(-1:1^d &, 0:3^d &), target type_send_p_p1
integer, dimension(-1:2^d &,-1:1^d &), target type_send_srl_p1
integer, dimension(-1:1^d &, 0:3^d &), target type_send_p_f
integer, dimension(:^d &,:^d &), pointer type_recv_srl
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_r_p2
subroutine create_bc_mpi_datatype(nwstart, nwbc)
integer, dimension(-1:2^d &,-1:1^d &), target type_recv_srl_p1
integer, dimension(-1:2^d &,-1:1^d &), target type_send_srl_p2
integer, dimension(:^d &,:^d &), pointer type_send_srl
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_p_p2
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_r_f
integer, dimension(:^d &,:^d &), pointer type_recv_r
integer, dimension(:^d &,:^d &), pointer type_send_r
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_r_p1
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_p_f
subroutine getbc(time, qdt, psb, nwstart, nwbc, req_diag)
do update ghost cells of all blocks including physical boundaries
integer, dimension(-1:1^d &, 0:3^d &), target type_send_p_p2
integer, dimension(-1:2^d &,-1:1^d &), target type_send_srl_f
integer, dimension(-1:1^d &,-1:1^d &), target type_send_r_p1
integer, dimension(:^d &,:^d &), pointer type_recv_p
integer, dimension(:^d &,:^d &), pointer type_send_p
integer, dimension(-1:1^d &,-1:1^d &), target type_send_r_f
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, parameter fs_tvdlf
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, parameter threestep
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
integer, parameter plevel_
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, parameter twostep
integer, parameter onestep
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
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...
Definition: mod_limiter.t:129
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
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 reconstructrmf(ixIL, iLL, idims, w, wRC)
subroutine advect1mf(method, dtin, dtfactor, idimLIM, qtC, psa, qt, psb)
subroutine getdtfff_courant(w, x, ixIL, ixOL, dtnew)
subroutine process1_gridmf(method, igrid, qdt, ixGL, idimLIM, qtC, wCT, qt, w)
subroutine divbclean(qdt, ixIL, ixOL, wCT, w, x)
Clean divergence of magnetic field by Janhunen's and Linde's source terms.
subroutine getfluxmf(w, x, ixIL, ixOL, idir, idim, f)
subroutine frictional_velocity(w, x, ixIL, ixOL, qvmax, qdt)
subroutine get_divb(w, ixIL, ixOL, divb)
Calculate div B within ixO.
double precision mf_tvdlfeps
TVDLF dissipation coefficient controls the dissipation term.
subroutine magnetofriction_init()
Initialize the module.
subroutine vhat(w, x, ixIL, ixOL, vhatmaxgrid)
subroutine advectmf(idimLIM, qt, qdt)
double precision, public mf_vmax
maximal limit of magnetofrictional velocity in cm s^-1 (Pomoell 2019 A&A)
double precision tmf
time in magnetofriction process
double precision mf_cy_max
subroutine tvdlfmf(qdt, ixIL, ixOL, idimLIM, qtC, wCT, qt, wnew, fC, dxD, x)
subroutine fdmf(qdt, ixIL, ixOL, idimLIM, qtC, wCT, qt, wnew, fC, dxD, x)
subroutine centdiff4mf(qdt, ixIL, ixOL, idimLIM, qtC, wCT, qt, w, fC, dxD, x)
double precision mf_cdivb
divb cleaning coefficient controls diffusion speed of divb
subroutine getcmaxfff(w, ixIL, ixOL, idims, cmax)
subroutine mf_velocity_update(dtfff)
double precision cmax_mype
maximal speed for fd scheme
double precision mf_cy
frictional velocity coefficient
double precision mf_tvdlfeps_min
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.
subroutine hancockmf(qdt, ixIL, ixOL, idimLIM, qtC, wCT, qt, wnew, dxD, x)
double precision cmax_global
maximal speed for fd scheme
subroutine upwindlrmf(ixIL, ixLL, ixRL, idim, w, wCT, wLC, wRC, x)
subroutine addgeometrymf(qdt, ixIL, ixOL, wCT, w, x)
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)...
subroutine reconstructlmf(ixIL, iLL, idims, w, wLC)
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:59
procedure(sub_convert), pointer phys_to_conserved
Definition: mod_physics.t:58
Module with all the methods that users can customize in AMRVAC.
procedure(p_no_args), pointer usr_before_main_loop