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