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==0 .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),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),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  dsurface=zero
268  do iigrid=1,igridstail; igrid=igrids(iigrid);
269  block=>ps(igrid)
270  dvolume(ixm^t)=block%dvolume(ixm^t)
271  do idims=1,ndim
272  hxm^ll=ixm^ll-kr(idims,^d);
273  dsurface(ixm^t)=dsurface(ixm^t)+block%surfaceC(hxm^t,idims)
274  end do
275  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
276  call mask_inner(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x)
277  sum_jbb_ipe = sum_jbb_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
278  ps(igrid)%x,1,patchwi)
279  sum_j_ipe = sum_j_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
280  ps(igrid)%x,2,patchwi)
281  f_i_ipe=f_i_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
282  ps(igrid)%x,3,patchwi)
283  sum_l_ipe = sum_l_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
284  ps(igrid)%x,4,patchwi)
285  end do
286  call mpi_allreduce(sum_jbb_ipe,sum_jbb,1,mpi_double_precision,&
287  mpi_sum,icomm,ierrmpi)
288  call mpi_allreduce(sum_j_ipe,sum_j,1,mpi_double_precision,mpi_sum,&
289  icomm,ierrmpi)
290  call mpi_allreduce(f_i_ipe,f_i,1,mpi_double_precision,mpi_sum,&
291  icomm,ierrmpi)
292  call mpi_allreduce(sum_l_ipe,sum_l,1,mpi_double_precision,mpi_sum,&
293  icomm,ierrmpi)
294  call mpi_allreduce(volumepe,volume,1,mpi_double_precision,mpi_sum,&
295  icomm,ierrmpi)
296  ! current- and volume-weighted average of the sine of the angle
297  ! between the magnetic field B and the current density J
298  cwsin_theta_new = sum_jbb/sum_j
299  ! volume-weighted average of the absolute value of the fractional
300  ! magnetic flux change
301  f_i = f_i/volume
302  sum_j=sum_j/volume
303  sum_l=sum_l/volume
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, intent(in) :: 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 (limiter_weno5)
1070  call weno5limiterl(ixi^l,il^l,idims,w,wlc,1)
1071  case (limiter_wenoz5)
1072  call weno5limiterl(ixi^l,il^l,idims,w,wlc,2)
1073  case default
1074 
1075  kxcmin^d=iximin^d; kxcmax^d=iximax^d-kr(idims,^d);
1076 
1077  wlc(kxc^s,1:nwflux) = w(kxc^s,1:nwflux)
1078 
1079  jxr^l=il^l+kr(idims,^d);
1080 
1081  ixcmax^d=jxrmax^d; ixcmin^d=ilmin^d-kr(idims,^d);
1082  jxc^l=ixc^l+kr(idims,^d);
1083 
1084  do iw=1,nwflux
1085  dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
1086 
1087  call dwlimiter2(dwc,ixi^l,ixc^l,idims,typelimiter,ldw=ldw)
1088 
1089  wlc(il^s,iw)=wlc(il^s,iw)+half*ldw(il^s)
1090  end do
1091  end select
1092 
1093  end subroutine reconstructlmf
1094 
1095  subroutine reconstructrmf(ixI^L,iL^L,idims,w,wRC)
1097  use mod_limiter
1098 
1099  integer, intent(in) :: ixI^L, iL^L, idims
1100  double precision, intent(in) :: w(ixi^s,1:nw)
1101 
1102  double precision, intent(out) :: wRC(ixi^s,1:nw)
1103 
1104  double precision :: rdw(ixi^s), dwC(ixi^s)
1105  integer :: jxR^L, ixC^L, jxC^L, kxC^L, kxR^L, iw
1106 
1107  select case (typelimiter)
1108  case (limiter_mp5)
1109  call mp5limiterr(ixi^l,il^l,idims,w,wrc)
1110  case (limiter_weno5)
1111  call weno5limiterr(ixi^l,il^l,idims,w,wrc,1)
1112  case (limiter_wenoz5)
1113  call weno5limiterr(ixi^l,il^l,idims,w,wrc,2)
1114  case default
1115 
1116  kxcmin^d=iximin^d; kxcmax^d=iximax^d-kr(idims,^d);
1117  kxr^l=kxc^l+kr(idims,^d);
1118 
1119  wrc(kxc^s,1:nwflux)=w(kxr^s,1:nwflux)
1120 
1121  jxr^l=il^l+kr(idims,^d);
1122  ixcmax^d=jxrmax^d; ixcmin^d=ilmin^d-kr(idims,^d);
1123  jxc^l=ixc^l+kr(idims,^d);
1124 
1125  do iw=1,nwflux
1126  dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
1127  call dwlimiter2(dwc,ixi^l,ixc^l,idims,typelimiter,rdw=rdw)
1128 
1129  wrc(il^s,iw)=wrc(il^s,iw)-half*rdw(jxr^s)
1130  end do
1131  end select
1132 
1133  end subroutine reconstructrmf
1134 
1135  subroutine centdiff4mf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,w,wold,fC,dx^D,x)
1136  ! Advance the flow variables from global_time to global_time+qdt within ixO^L by
1137  ! fourth order centered differencing in space
1138  ! for the dw/dt+dF_i(w)/dx_i=S type equation.
1139  ! wCT contains the time centered variables at time qtC for flux and source.
1140  ! w is the old value at qt on input and the new value at qt+qdt on output.
1142 
1143  integer, intent(in) :: ixI^L, ixO^L, idim^LIM
1144  double precision, intent(in) :: qdt, qtC, qt, dx^D
1145  double precision :: wCT(ixi^s,1:nw), w(ixi^s,1:nw), wold(ixi^s,1:nw)
1146  double precision, intent(in) :: x(ixi^s,1:ndim)
1147  double precision :: fC(ixi^s,1:ndir,1:ndim)
1148 
1149  double precision :: v(ixi^s,ndim), f(ixi^s)
1150  double precision, dimension(ixI^S,1:nw) :: wLC, wRC
1151  double precision, dimension(ixI^S) :: vLC, vRC,cmaxLC,cmaxRC
1152  double precision :: dxinv(1:ndim)
1153  integer :: idims, idir, idirmin,ix^D
1154  integer :: ix^L, hxO^L, ixC^L, jxC^L, hxC^L, kxC^L, kkxC^L, kkxR^L
1155 
1156  ! two extra layers are needed in each direction for which fluxes are added.
1157  ix^l=ixo^l;
1158  do idims= idim^lim
1159  ix^l=ix^l^ladd2*kr(idims,^d);
1160  end do
1161 
1162  if (ixi^l^ltix^l|.or.|.or.) then
1163  call mpistop("Error in evolve_CentDiff4: Non-conforming input limits")
1164  end if
1165  ^d&dxinv(^d)=-qdt/dx^d;
1166 
1167  ! Add fluxes to w
1168  do idims= idim^lim
1169  block%iw0=idims
1170  ix^l=ixo^l^ladd2*kr(idims,^d);
1171  hxo^l=ixo^l-kr(idims,^d);
1172 
1173  ixcmin^d=hxomin^d; ixcmax^d=ixomax^d;
1174  hxc^l=ixc^l-kr(idims,^d);
1175  jxc^l=ixc^l+kr(idims,^d);
1176  kxc^l=ixc^l+2*kr(idims,^d);
1177 
1178  kkxcmin^d=iximin^d; kkxcmax^d=iximax^d-kr(idims,^d);
1179  kkxr^l=kkxc^l+kr(idims,^d);
1180  wrc(kkxc^s,1:nwflux)=wct(kkxr^s,1:nwflux)
1181  wlc(kkxc^s,1:nwflux)=wct(kkxc^s,1:nwflux)
1182 
1183  call upwindlrmf(ixi^l,ixc^l,ixc^l,idims,wct,wct,wlc,wrc,x)
1184 
1185  ! Calculate velocities from upwinded values
1186  call getcmaxfff(wlc,ixg^ll,ixc^l,idims,cmaxlc)
1187  call getcmaxfff(wrc,ixg^ll,ixc^l,idims,cmaxrc)
1188  ! now take the maximum of left and right states
1189  vlc(ixc^s)=max(cmaxrc(ixc^s),cmaxlc(ixc^s))
1190 
1191  do idir=1,ndir
1192  ! Get non-transported flux
1193  call getfluxmf(wct,x,ixi^l,ix^l,idir,idims,f)
1194  ! Center flux to interface
1195  ! f_i+1/2= (-f_(i+2) +7 f_(i+1) + 7 f_i - f_(i-1))/12
1196  fc(ixc^s,idir,idims)=(-f(kxc^s)+7.0d0*(f(jxc^s)+f(ixc^s))-f(hxc^s))/12.0d0
1197  ! add rempel dissipative flux, only second order version for now
1198  ! one could gradually reduce the dissipative flux to improve solutions
1199  ! for computing steady states (Keppens et al. 2012, JCP)
1200  fc(ixc^s,idir,idims)=fc(ixc^s,idir,idims)-mf_tvdlfeps*tvdlfeps*half*vlc(ixc^s) &
1201  *(wrc(ixc^s,mag(idir))-wlc(ixc^s,mag(idir)))
1202 
1203  if (slab_uniform) then
1204  fc(ixc^s,idir,idims)=dxinv(idims)*fc(ixc^s,idir,idims)
1205  ! result: f_(i+1/2)-f_(i-1/2) = [-f_(i+2)+8(f_(i+1)+f_(i-1))-f_(i-2)]/12
1206  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+(fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))
1207  else
1208  fc(ixc^s,idir,idims)=-qdt*block%surfaceC(ixc^s,idims)*fc(ixc^s,idir,idims)
1209  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+ &
1210  (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))/block%dvolume(ixo^s)
1211  end if
1212  end do !next idir
1213  end do !next idims
1214 
1215  if (.not.slab) call addgeometrymf(qdt,ixi^l,ixo^l,wct,w,x)
1216  call divbclean(qdt,ixi^l,ixo^l,wct,w,x)
1217 
1218  end subroutine centdiff4mf
1219 
1220  subroutine getdtfff_courant(w,x,ixI^L,ixO^L,dtnew)
1221  ! compute CFL limited dt (for variable time stepping)
1223 
1224  integer, intent(in) :: ixI^L, ixO^L
1225  double precision, intent(in) :: x(ixi^s,1:ndim)
1226  double precision, intent(inout) :: w(ixi^s,1:nw), dtnew
1227 
1228  double precision :: courantmax, dxinv(1:ndim)
1229  double precision :: cmax(ixi^s),tmp(ixi^s),alfven(ixi^s)
1230  integer :: idims
1231 
1232  dtnew=bigdouble
1233  courantmax=0.d0
1234  ^d&dxinv(^d)=1.d0/dxlevel(^d);
1235 
1236  do idims=1,ndim
1237  call getcmaxfff(w,ixi^l,ixo^l,idims,cmax)
1238  cmax_mype = max(cmax_mype,maxval(cmax(ixo^s)))
1239  if (.not.slab_uniform) then
1240  tmp(ixo^s)=cmax(ixo^s)/block%dx(ixo^s,idims)
1241  courantmax=max(courantmax,maxval(tmp(ixo^s)))
1242  else
1243  tmp(ixo^s)=cmax(ixo^s)*dxinv(idims)
1244  courantmax=max(courantmax,maxval(tmp(ixo^s)))
1245  end if
1246  end do
1247  ! courantmax='max( c/dx)'
1248  if (courantmax>smalldouble) dtnew=min(dtnew,mf_cc/courantmax)
1249 
1250  end subroutine getdtfff_courant
1251 
1252  subroutine getcmaxfff(w,ixI^L,ixO^L,idims,cmax)
1254 
1255  logical :: new_cmax,needcmin
1256  integer, intent(in) :: ixI^L, ixO^L, idims
1257  double precision, intent(in) :: w(ixi^s,1:nw)
1258  double precision, intent(out) :: cmax(ixi^s)
1259 
1260  ! calculate alfven speed
1261  if(b0field) then
1262  cmax(ixo^s)=sqrt(sum((w(ixo^s,mag(:))+block%b0(ixo^s,:,0))**2,dim=ndim+1)/w(ixo^s,rho_))
1263  else
1264  cmax(ixo^s)=sqrt(sum(w(ixo^s,mag(:))**2,dim=ndim+1)/w(ixo^s,rho_))
1265  endif
1266  cmax(ixo^s)=cmax(ixo^s)+abs(w(ixo^s,mom(idims)))
1267 
1268  end subroutine getcmaxfff
1269 
1270  !> Clean divergence of magnetic field by Janhunen's and Linde's source terms
1271  subroutine divbclean(qdt,ixI^L,ixO^L,wCT,w,x)
1273  use mod_geometry
1274 
1275  integer, intent(in) :: ixI^L, ixO^L
1276  double precision, intent(in) :: x(ixi^s,1:ndim),wCT(ixi^s,1:nw),qdt
1277  double precision, intent(inout) :: w(ixi^s,1:nw)
1278  integer :: idims, ix^L, ixp^L, i^D, iside
1279  double precision :: divb(ixi^s),graddivb(ixi^s),bdivb(ixi^s,1:ndir)
1280 
1281  ! Calculate div B
1282  ix^l=ixo^l^ladd1;
1283  call get_divb(wct,ixi^l,ix^l,divb)
1284 
1285  ixp^l=ixo^l;
1286 
1287  ! Add Linde's diffusive terms
1288  do idims=1,ndim
1289  ! Calculate grad_idim(divb)
1290  call gradient(divb,ixi^l,ixp^l,idims,graddivb)
1291 
1292  ! Multiply by Linde's eta*dt = divbdiff*(c_max*dx)*dt = divbdiff*dx**2
1293  if (slab_uniform) then
1294  graddivb(ixp^s)=graddivb(ixp^s)*mf_cdivb/(^d&1.0d0/dxlevel(^d)**2+)
1295  else
1296  graddivb(ixp^s)=graddivb(ixp^s)*mf_cdivb &
1297  /(^d&1.0d0/block%dx(ixp^s,^d)**2+)
1298  end if
1299  ! B_idim += eta*grad_idim(divb)
1300  ! with Janhunen's term
1301  !w(ixp^S,mag(idims))=w(ixp^S,mag(idims))+&
1302  ! graddivb(ixp^S)-qdt*w(ixp^S,mom(idims))*divb(ixp^S)
1303  ! without Janjunen's term
1304  w(ixp^s,mag(idims))=w(ixp^s,mag(idims))+&
1305  graddivb(ixp^s)
1306  end do
1307 
1308  end subroutine divbclean
1309 
1310  subroutine addgeometrymf(qdt,ixI^L,ixO^L,wCT,w,x)
1311  ! Add geometrical source terms to w
1313  use mod_geometry
1314 
1315  integer, intent(in) :: ixI^L, ixO^L
1316  double precision, intent(in) :: qdt, x(ixi^s,1:ndim)
1317  double precision, intent(inout) :: wCT(ixi^s,1:nw), w(ixi^s,1:nw)
1318  !.. local ..
1319  double precision :: tmp(ixi^s)
1320  integer :: iw
1321  integer :: mr_,mphi_ ! Polar var. names
1322  integer :: br_,bphi_
1323 
1324  mr_=mom(1); mphi_=mom(1)-1+phi_ ! Polar var. names
1325  br_=mag(1); bphi_=mag(1)-1+phi_
1326 
1327  select case (coordinate)
1328  case (cylindrical)
1329  if(phi_>0) then
1330  ! s[Bphi]=(Bphi*vr-Br*vphi)/radius
1331  tmp(ixo^s)=(wct(ixo^s,bphi_)*wct(ixo^s,mom(1)) &
1332  -wct(ixo^s,br_)*wct(ixo^s,mom(3)))
1333  w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt*tmp(ixo^s)/x(ixo^s,1)
1334  end if
1335  case (spherical)
1336  {^nooned
1337  ! s[b2]=(vr*Btheta-vtheta*Br)/r
1338  ! + cot(theta)*psi/r
1339  tmp(ixo^s)= wct(ixo^s,mom(1))*wct(ixo^s,mag(2)) &
1340  -wct(ixo^s,mom(2))*wct(ixo^s,mag(1))
1341  if (b0field) then
1342  tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,mom(1))*block%b0(ixo^s,2,0) &
1343  -wct(ixo^s,mom(2))*block%b0(ixo^s,1,0)
1344  end if
1345  ! Divide by radius and add to w
1346  w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
1347  }
1348  if(ndir==3) then
1349  ! s[b3]=(vr*Bphi-vphi*Br)/r
1350  ! -cot(theta)*(vphi*Btheta-vtheta*Bphi)/r
1351  tmp(ixo^s)=wct(ixo^s,mom(1))*wct(ixo^s,mag(3)) &
1352  -wct(ixo^s,mom(3))*wct(ixo^s,mag(1)){^nooned &
1353  -(wct(ixo^s,mom(3))*wct(ixo^s,mag(2)) &
1354  -wct(ixo^s,mom(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
1355  /dsin(x(ixo^s,2)) }
1356  if (b0field) then
1357  tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,mom(1))*block%b0(ixo^s,3,0) &
1358  -wct(ixo^s,mom(3))*block%b0(ixo^s,1,0){^nooned &
1359  -(wct(ixo^s,mom(3))*block%b0(ixo^s,2,0) &
1360  -wct(ixo^s,mom(2))*block%b0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
1361  /dsin(x(ixo^s,2)) }
1362  end if
1363  ! Divide by radius and add to w
1364  w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
1365  end if
1366  end select
1367 
1368  end subroutine addgeometrymf
1369 
1370  !> Calculate idirmin and the idirmin:3 components of the common current array
1371  !> make sure that dxlevel(^D) is set correctly.
1372  subroutine get_current(w,ixI^L,ixO^L,idirmin,current)
1374  use mod_geometry
1375 
1376  integer :: idirmin0
1377  integer :: ixO^L, idirmin, ixI^L
1378  double precision :: w(ixi^s,1:nw)
1379  integer :: idir
1380 
1381  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
1382  double precision :: current(ixi^s,7-2*ndir:3),bvec(ixi^s,1:ndir)
1383 
1384  idirmin0 = 7-2*ndir
1385 
1386  bvec(ixi^s,1:ndir)=w(ixi^s,mag(1:ndir))
1387 
1388  call curlvector(bvec,ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
1389 
1390  if(b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
1391  block%J0(ixo^s,idirmin0:3)
1392 
1393  end subroutine get_current
1394 
1395  !> Calculate div B within ixO
1396  subroutine get_divb(w,ixI^L,ixO^L,divb)
1399 
1400  integer, intent(in) :: ixI^L, ixO^L
1401  double precision, intent(in) :: w(ixi^s,1:nw)
1402  double precision :: divb(ixi^s)
1403 
1404  double precision :: bvec(ixi^s,1:ndir)
1405 
1406  bvec(ixi^s,:)=w(ixi^s,mag(:))
1407 
1408  select case(typediv)
1409  case("central")
1410  call divvector(bvec,ixi^l,ixo^l,divb)
1411  case("limited")
1412  call divvectors(bvec,ixi^l,ixo^l,divb)
1413  end select
1414  end subroutine get_divb
1415 
1416 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
integer, parameter limiter_wenoz5
Definition: mod_limiter.t:31
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
integer, parameter limiter_weno5
Definition: mod_limiter.t:29
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:119
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)