MPI-AMRVAC  2.0
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)
115 
116  time_in=mpi_wtime()
117  if(mype==0) write(*,*) 'Evolving to force-free field using magnetofricitonal method...'
118  if(prolongprimitive) call mpistop('use prolongprimitive=.false. in MF module')
119  mf_advance=.false.
120  dtfff=1.d-2
121  tmpt=global_time
122  tmpit=it
124  i=it
125  if(snapshotini==-1 .and. i==0) then
126  call saveamrfile(1)
127  call saveamrfile(2)
128  end if
129  mf_advance=.true.
130  ! point bc mpi datatype to partial type for magnetic field
137  ! create bc mpi datatype for ghostcells update
138  call create_bc_mpi_datatype(mag(1)-1,ndir)
139  ! point bc mpi datatype to partial type for velocity field
146  ! create bc mpi datatype for ghostcells update
147  call create_bc_mpi_datatype(mom(1)-1,ndir)
148  ! convert conservative variables to primitive ones which are used during MF
149  do iigrid=1,igridstail; igrid=igrids(iigrid);
150  call phys_to_primitive(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x)
151  end do
152  ! calculate magnetofrictional velocity
153  call mf_velocity_update(dtfff)
154  ! update velocity in ghost cells
155  call getbc(tmf,0.d0,ps,mom(1)-1,ndir)
156  ! calculate initial values of metrics
157  if(i==0) then
158  call metrics
159  call printlog_mf
160  end if
161  ! magnetofrictional loops
162  do
163  ! calculate time step based on Cmax= Alfven speed + abs(frictional speed)
164  dtfff_pe=bigdouble
165  cmax_mype=zero
166  do iigrid=1,igridstail; igrid=igrids(iigrid);
167  block=>ps(igrid)
168  pso(igrid)%w(ixg^t,mag(:))=ps(igrid)%w(ixg^t,mag(:))
169  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
170  call getdtfff_courant(ps(igrid)%w,ps(igrid)%x,ixg^ll,ixm^ll,dtnew)
171  dtfff_pe=min(dtfff_pe,dtnew)
172  end do
173  call mpi_allreduce(dtfff_pe,dtfff,1,mpi_double_precision,mpi_min, &
174  icomm,ierrmpi)
175  call mpi_allreduce(cmax_mype,cmax_global,1,mpi_double_precision,&
176  mpi_max,icomm,ierrmpi)
177 
178  ! =======
179  ! evolve
180  ! =======
181  call advectmf(1,ndim,tmf,dtfff)
182 
183  if(i>=10000) then
184  mf_tvdlfeps = 0.9998d0 * mf_tvdlfeps
186  end if
187  if(i<=60000) then
188  mf_cy=1.0001d0*mf_cy
189  mf_cy = min(mf_cy_max,mf_cy)
190  end if
191 
192  i=i+1
193  tmf=tmf+dtfff
194  if(mod(i,10)==0) then
195  ! calculate metrics
196  call metrics
197  call printlog_mf
198  end if
199  if(mod(i,mf_ditsave)==0) then
200  it=i
202  do iigrid=1,igridstail; igrid=igrids(iigrid);
203  call phys_to_conserved(ixg^ll,ixg^ll,ps(igrid)%w,ps(igrid)%x)
204  end do
205  mf_advance=.false.
206  call saveamrfile(1)
207  call saveamrfile(2)
208  do iigrid=1,igridstail; igrid=igrids(iigrid);
209  call phys_to_primitive(ixg^ll,ixg^ll,ps(igrid)%w,ps(igrid)%x)
210  end do
211  mf_advance=.true.
212  if(mype==0) then
213  write(*,*) "itmf=",i
214  write(*,*) '<CW sin theta>:',cwsin_theta_new
215  write(*,*) '<f_i>:',f_i
216  write(*,*) '----------------------------------------------------------'
217  end if
218  end if
219  ! reconstruct AMR grid every 10 step
220  if(mod(i,ditregrid)==0 .and. refine_max_level>1) call resettree
221  if (i>=mf_it_max) then
222  if(mod(i,10)/=0) then
223  ! calculate metrics
224  call metrics
225  call printlog_mf
226  end if
227  if(mype==0) then
228  write (*,*) 'Reach maximum iteration step!'
229  write (*,*) 'The total iteration step is:', i
230  end if
231  exit
232  end if
233  enddo
234  ! point bc mpi data type back to full type for MHD
241  bcphys=.true.
242  ! set velocity back to zero and convert primitive variables back to conservative ones
243  do iigrid=1,igridstail; igrid=igrids(iigrid);
244  ps(igrid)%w(ixg^t,mom(:))=zero
245  call phys_to_conserved(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x)
246  end do
247  global_time=tmpt
248  it=tmpit
249  if (mype==0) call mpi_file_close(fhmf,ierrmpi)
250  mf_advance=.false.
251  if(mype==0) write(*,*) 'Magnetofriction phase took : ',mpi_wtime()-time_in,' sec'
252  contains
253 
254  subroutine metrics
256  sum_jbb_ipe = 0.d0
257  sum_j_ipe = 0.d0
258  sum_l_ipe = 0.d0
259  f_i_ipe = 0.d0
260  volumepe=0.d0
261  do iigrid=1,igridstail; igrid=igrids(iigrid);
262  block=>ps(igrid)
263  if(slab) then
264  dvone={rnode(rpdx^d_,igrid)|*}
265  dvolume(ixm^t)=dvone
266  dsurface(ixm^t)=two*(^d&dvone/rnode(rpdx^d_,igrid)+)
267  else
268  dvolume(ixm^t)=block%dvolume(ixm^t)
269  dsurface(ixm^t)= sum(block%surfaceC(ixm^t,:),dim=ndim+1)
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  end if
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 
305  end subroutine metrics
306 
307  subroutine mask_inner(ixI^L,ixO^L,w,x)
309  integer, intent(in) :: ixI^L,ixO^L
310  double precision, intent(in):: w(ixi^s,nw),x(ixi^s,1:ndim)
311  double precision :: xO^L
312  integer :: ix^D
313 
314  {xomin^d = xprobmin^d + 0.05d0*(xprobmax^d-xprobmin^d)\}
315  {xomax^d = xprobmax^d - 0.05d0*(xprobmax^d-xprobmin^d)\}
316  if(slab .or. slab_stretched) then
317  xomin^nd = xprobmin^nd
318  else
319  xomin1 = xprobmin1
320  end if
321 
322  {do ix^db=ixomin^db,ixomax^db\}
323  if({ x(ix^dd,^d) > xomin^d .and. x(ix^dd,^d) < xomax^d | .and. }) then
324  patchwi(ix^d)=.true.
325  volumepe=volumepe+dvolume(ix^d)
326  else
327  patchwi(ix^d)=.false.
328  endif
329  {end do\}
330 
331  end subroutine mask_inner
332 
333  subroutine printlog_mf
334  integer :: amode, status(mpi_status_size)
335  character(len=800) :: filename,filehead
336  character(len=2048) :: line,datastr
337  logical, save :: logmfopened=.false.
338 
339  if(mype==0) then
340  if(.not.logmfopened) then
341  ! generate filename
342  write(filename,"(a,a)") trim(base_filename), "_mflog.csv"
343 
344  amode=ior(mpi_mode_create,mpi_mode_wronly)
345  amode=ior(amode,mpi_mode_append)
346  call mpi_file_open(mpi_comm_self,filename,amode,mpi_info_null,fhmf,ierrmpi)
347  logmfopened=.true.
348  filehead=" itmf, dt, <f_i>, <CW sin theta>, <Current>, <Lorenz force>"
349  call mpi_file_write(fhmf,filehead,len_trim(filehead), &
350  mpi_character,status,ierrmpi)
351  call mpi_file_write(fhmf,achar(10),1,mpi_character,status,ierrmpi)
352  end if
353  line=''
354  write(datastr,'(i6,a)') i,','
355  line=trim(line)//trim(datastr)
356  write(datastr,'(es13.6,a)') dtfff,','
357  line=trim(line)//trim(datastr)
358  write(datastr,'(es13.6,a)') f_i,','
359  line=trim(line)//trim(datastr)
360  write(datastr,'(es13.6,a)') cwsin_theta_new,','
361  line=trim(line)//trim(datastr)
362  write(datastr,'(es13.6,a)') sum_j,','
363  line=trim(line)//trim(datastr)
364  write(datastr,'(es13.6)') sum_l
365  line=trim(line)//trim(datastr)//new_line('A')
366  call mpi_file_write(fhmf,line,len_trim(line),mpi_character,status,ierrmpi)
367  end if
368 
369  end subroutine printlog_mf
370 
371  function integral_grid_mf(ixI^L,ixO^L,w,x,iw,patchwi)
373 
374  integer, intent(in) :: ixI^L,ixO^L,iw
375  double precision, intent(in) :: x(ixi^s,1:ndim)
376  double precision :: w(ixi^s,nw+nwauxio)
377  logical, intent(in) :: patchwi(ixi^s)
378 
379  double precision, dimension(ixI^S,1:ndir) :: bvec,qvec,current
380  double precision :: integral_grid_mf,tmp(ixi^s),b_mag(ixi^s)
381  integer :: ix^D,i,idirmin,idir,jdir,kdir
382 
383  integral_grid_mf=0.d0
384  select case(iw)
385  case(1)
386  ! Sum(dvolume*|JxB|/|B|)
387  if(b0field) then
388  bvec(ixi^s,:)=w(ixi^s,mag(:))+block%b0(ixi^s,mag(:),0)
389  else
390  bvec(ixi^s,:)=w(ixi^s,mag(:))
391  endif
392  call get_current(w,ixi^l,ixo^l,idirmin,current)
393  ! calculate Lorentz force
394  qvec(ixo^s,1:ndir)=zero
395  do idir=1,ndir; do jdir=1,ndir; do kdir=idirmin,3
396  if(lvc(idir,jdir,kdir)/=0)then
397  tmp(ixo^s)=current(ixo^s,jdir)*bvec(ixo^s,kdir)
398  if(lvc(idir,jdir,kdir)==1)then
399  qvec(ixo^s,idir)=qvec(ixo^s,idir)+tmp(ixo^s)
400  else
401  qvec(ixo^s,idir)=qvec(ixo^s,idir)-tmp(ixo^s)
402  endif
403  endif
404  enddo; enddo; enddo
405 
406  {do ix^db=ixomin^db,ixomax^db\}
407  if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+sqrt(sum(qvec(ix^d,:)**2)/&
408  sum(bvec(ix^d,:)**2))*dvolume(ix^d)
409  {end do\}
410  case(2)
411  ! Sum(dvolume*|J|)
412  call get_current(w,ixi^l,ixo^l,idirmin,current)
413  {do ix^db=ixomin^db,ixomax^db\}
414  if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+sqrt(sum(current(ix^d,:)**2))*&
415  dvolume(ix^d)
416  {end do\}
417  case(3)
418  ! f_i solenoidal property of B: (dvolume |div B|)/(dsurface |B|)
419  ! Sum(dvolume*f_i)
420  if(b0field) then
421  bvec(ixi^s,:)=w(ixi^s,mag(:))+block%b0(ixi^s,mag(:),0)
422  else
423  bvec(ixi^s,:)=w(ixi^s,mag(:))
424  endif
425  call divvector(bvec,ixi^l,ixo^l,tmp)
426  {do ix^db=ixomin^db,ixomax^db\}
427  if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+abs(tmp(ix^d))*&
428  dvolume(ix^d)**2/sqrt(sum(bvec(ix^d,:)**2))/dsurface(ix^d)
429  {end do\}
430  case(4)
431  ! Sum(|JxB|)
432  if(b0field) then
433  bvec(ixi^s,:)=w(ixi^s,mag(:))+block%b0(ixi^s,mag(:),0)
434  else
435  bvec(ixi^s,:)=w(ixi^s,mag(:))
436  endif
437  call curlvector(bvec,ixi^l,ixo^l,current,idirmin,1,ndir)
438  ! calculate Lorentz force
439  qvec(ixo^s,1:ndir)=zero
440  do idir=1,ndir; do jdir=1,ndir; do kdir=idirmin,3
441  if(lvc(idir,jdir,kdir)/=0)then
442  tmp(ixo^s)=current(ixo^s,jdir)*bvec(ixo^s,kdir)
443  if(lvc(idir,jdir,kdir)==1)then
444  qvec(ixo^s,idir)=qvec(ixo^s,idir)+tmp(ixo^s)
445  else
446  qvec(ixo^s,idir)=qvec(ixo^s,idir)-tmp(ixo^s)
447  endif
448  endif
449  enddo; enddo; enddo
450 
451  {do ix^db=ixomin^db,ixomax^db\}
452  if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+sqrt(sum(qvec(ix^d,:)**2))*dvolume(ix^d)
453  {end do\}
454  end select
455  return
456  end function integral_grid_mf
457 
458  end subroutine magnetofriction
459 
460  subroutine mf_velocity_update(dtfff)
462 
463  double precision, intent(in) :: dtfff
464  integer :: i,iigrid, igrid
465  double precision :: vhatmax,vhatmax_pe,vhatmaxgrid
466 
467  vhatmax_pe=smalldouble
468  do iigrid=1,igridstail; igrid=igrids(iigrid);
469  block=>ps(igrid)
470  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
471  call vhat(ps(igrid)%w,ps(igrid)%x,ixg^ll,ixm^ll,vhatmaxgrid)
472  vhatmax_pe=max(vhatmax_pe,vhatmaxgrid)
473  end do
474  call mpi_allreduce(vhatmax_pe,vhatmax,1,mpi_double_precision,mpi_max, &
475  icomm,ierrmpi)
476  do iigrid=1,igridstail; igrid=igrids(iigrid);
477  block=>ps(igrid)
478  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
479  ! calculate frictional velocity
480  call frictional_velocity(ps(igrid)%w,ps(igrid)%x,ixg^ll,ixm^ll,vhatmax,dtfff)
481  end do
482 
483  end subroutine mf_velocity_update
484 
485  subroutine vhat(w,x,ixI^L,ixO^L,vhatmaxgrid)
486  ! Calculate v_hat
488 
489  integer, intent(in) :: ixI^L, ixO^L
490  double precision, intent(inout) :: w(ixi^s,nw)
491  double precision, intent(in) :: x(ixi^s,1:ndim)
492  double precision, intent(out) :: vhatmaxgrid
493 
494  double precision :: current(ixi^s,7-2*ndir:3),tmp(ixi^s),dxhm
495  double precision :: dxhms(ixo^s)
496  integer :: idirmin,idir,jdir,kdir
497 
498  call get_current(w,ixi^l,ixo^l,idirmin,current)
499  w(ixi^s,mom(:))=0.d0
500  ! calculate Lorentz force
501  do idir=1,ndir; do jdir=1,ndir; do kdir=idirmin,3
502  if(lvc(idir,jdir,kdir)/=0)then
503  if(b0field) then
504  tmp(ixo^s)=current(ixo^s,jdir)*(w(ixo^s,mag(kdir))+block%b0(ixo^s,kdir,0))
505  else
506  tmp(ixo^s)=current(ixo^s,jdir)*w(ixo^s,mag(kdir))
507  endif
508  if(lvc(idir,jdir,kdir)==1)then
509  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))+tmp(ixo^s)
510  else
511  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-tmp(ixo^s)
512  endif
513  endif
514  enddo; enddo; enddo
515 
516  if(b0field) then
517  tmp(ixo^s)=sum((w(ixo^s,mag(:))+block%b0(ixo^s,:,0))**2,dim=ndim+1) ! |B|**2
518  else
519  tmp(ixo^s)=sum(w(ixo^s,mag(:))**2,dim=ndim+1) ! |B|**2
520  endif
521 
522  if(slab) then
523  dxhm=dble(ndim)/(^d&1.0d0/dxlevel(^d)+)
524  do idir=1,ndir
525  w(ixo^s,mom(idir))=dxhm*w(ixo^s,mom(idir))/tmp(ixo^s)
526  end do
527  else
528  dxhms(ixo^s)=dble(ndim)/sum(1.d0/block%dx(ixo^s,:),dim=ndim+1)
529  do idir=1,ndir
530  w(ixo^s,mom(idir))=dxhms(ixo^s)*w(ixo^s,mom(idir))/tmp(ixo^s)
531  end do
532  end if
533  vhatmaxgrid=maxval(sqrt(sum(w(ixo^s,mom(:))**2,dim=ndim+1)))
534 
535  end subroutine vhat
536 
537  subroutine frictional_velocity(w,x,ixI^L,ixO^L,qvmax,qdt)
539 
540  integer, intent(in) :: ixI^L, ixO^L
541  double precision, intent(in) :: x(ixi^s,1:ndim),qdt,qvmax
542  double precision, intent(inout) :: w(ixi^s,1:nw)
543 
544  double precision :: dxhm,disbd(6),bfzone^D
545  double precision :: dxhms(ixo^s)
546  integer :: ix^D, idir
547  logical :: buffer
548 
549  if(slab) then
550  dxhm=dble(ndim)/(^d&1.0d0/dxlevel(^d)+)
551  dxhm=mf_cc*mf_cy/qvmax*dxhm/qdt
552  ! dxhm=mf_cc*mf_cy/qvmax
553  w(ixo^s,mom(:))=w(ixo^s,mom(:))*dxhm
554  else
555  dxhms(ixo^s)=dble(ndim)/sum(1.d0/block%dx(ixo^s,:),dim=ndim+1)
556  dxhms(ixo^s)=mf_cc*mf_cy/qvmax*dxhms(ixo^s)/qdt
557  ! dxhm=mf_cc*mf_cy/qvmax
558  do idir=1,ndir
559  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))*dxhms(ixo^s)
560  end do
561  end if
562 
563 {^ifthreed
564  bfzone1=0.05d0*(xprobmax1-xprobmin1)
565  bfzone2=0.05d0*(xprobmax2-xprobmin2)
566  bfzone3=0.05d0*(xprobmax3-xprobmin3)
567  {do ix^db=ixomin^db,ixomax^db\}
568  disbd(1)=x(ix^d,1)-xprobmin1
569  disbd(2)=xprobmax1-x(ix^d,1)
570  disbd(3)=x(ix^d,2)-xprobmin2
571  disbd(4)=xprobmax2-x(ix^d,2)
572  disbd(5)=x(ix^d,3)-xprobmin1
573  disbd(6)=xprobmax3-x(ix^d,3)
574 
575  if(typeaxial=='slab'.or.typeaxial=='slabstretch') then
576  if(disbd(1)<bfzone1) then
577  w(ix^d,mom(:))=(1.d0-((bfzone1-disbd(1))/bfzone1)**2)*w(ix^d,mom(:))
578  endif
579  else
580  if(disbd(5)<bfzone3) then
581  w(ix^d,mom(:))=(1.d0-((bfzone3-disbd(5))/bfzone3)**2)*w(ix^d,mom(:))
582  endif
583  end if
584  if(disbd(2)<bfzone1) then
585  w(ix^d,mom(:))=(1.d0-((bfzone1-disbd(2))/bfzone1)**2)*w(ix^d,mom(:))
586  endif
587  if(disbd(3)<bfzone2) then
588  w(ix^d,mom(:))=(1.d0-((bfzone2-disbd(3))/bfzone2)**2)*w(ix^d,mom(:))
589  endif
590  if(disbd(4)<bfzone2) then
591  w(ix^d,mom(:))=(1.d0-((bfzone2-disbd(4))/bfzone2)**2)*w(ix^d,mom(:))
592  endif
593  if(disbd(6)<bfzone3) then
594  w(ix^d,mom(:))=(1.d0-((bfzone3-disbd(6))/bfzone3)**2)*w(ix^d,mom(:))
595  endif
596  {end do\}
597 }
598  end subroutine frictional_velocity
599 
600  subroutine advectmf(idim^LIM,qt,qdt)
601  ! integrate all grids by one step of its delta(global_time)
602  ! This subroutine is in VAC terminology equivalent to
603  ! `advect' (with the difference that it will `advect' all grids)
605  use mod_fix_conserve
606 
607  integer, intent(in) :: idim^LIM
608  double precision, intent(in) :: qt, qdt
609 
610  integer :: iigrid, igrid
611 
612  call init_comm_fix_conserve(idim^lim,ndir)
614 
615  ! copy w instead of wold because of potential use of dimsplit or sourcesplit
616  do iigrid=1,igridstail; igrid=igrids(iigrid);
617  ps1(igrid)%w=ps(igrid)%w
618  end do
619 
620  istep=0
621 
622  select case (time_integrator)
623  case ("onestep")
624  call advect1mf(flux_scheme,qdt,one,idim^lim,qt,ps1,qt,ps)
625  case ("twostep")
626  ! predictor step
627  fix_conserve_at_step = .false.
628  call advect1mf(typepred1,qdt,half,idim^lim,qt,ps,qt,ps1)
629  ! corrector step
631  call advect1mf(flux_scheme,qdt,one,idim^lim,qt+half*qdt,ps1,qt,ps)
632  case ("threestep")
633  ! three step Runge-Kutta in accordance with Gottlieb & Shu 1998
634  call advect1mf(flux_scheme,qdt,one,idim^lim,qt,ps,qt,ps1)
635 
636  do iigrid=1,igridstail; igrid=igrids(iigrid);
637  ps2(igrid)%w(ixg^t,1:nwflux)=0.75d0*ps(igrid)%w(ixg^t,1:nwflux)+0.25d0*&
638  ps1(igrid)%w(ixg^t,1:nwflux)
639  if (nw>nwflux) ps2(igrid)%w(ixg^t,nwflux+1:nw) = &
640  ps(igrid)%w(ixg^t,nwflux+1:nw)
641  end do
642 
643  call advect1mf(flux_scheme,qdt,0.25d0,idim^lim,qt+qdt,ps1,qt+dt*0.25d0,ps2)
644 
645  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
646  ps(igrid)%w(ixg^t,1:nwflux)=1.0d0/3.0d0*ps(igrid)%w(ixg^t,1:nwflux)+&
647  2.0d0/3.0d0*ps2(igrid)%w(ixg^t,1:nwflux)
648  end do
649 
650  call advect1mf(flux_scheme,qdt,2.0d0/3.0d0,idim^lim,qt+qdt/2.0d0,ps2,qt+qdt/3.0d0,ps)
651  case default
652  write(unitterm,*) "time_integrator=",time_integrator
653  write(unitterm,*) "Error in advectmf: Unknown time integration method"
654  call mpistop("Correct time_integrator")
655  end select
656 
657  end subroutine advectmf
658 
659  subroutine advect1mf(method,dtin,dtfactor,idim^LIM,qtC,psa,qt,psb)
660  ! Integrate all grids by one partial step
661  ! This subroutine is equivalent to VAC's `advect1', but does
662  ! the advection for all grids
665  use mod_fix_conserve
666 
667  integer, intent(in) :: idim^LIM
668  type(state) :: psa(max_blocks)! Compute fluxes based on this state
669  type(state) :: psb(max_blocks) ! update on this state
670  double precision, intent(in) :: dtin,dtfactor, qtC, qt
671  character(len=*), intent(in) :: method(nlevelshi)
672 
673  double precision :: qdt
674  integer :: iigrid, igrid, level, i^D
675  logical :: setigrid
676 
677  istep=istep+1
678 
679  ! loop over all grids to arrive at equivalent
680  qdt=dtfactor*dtin
681  do iigrid=1,igridstail; igrid=igrids(iigrid);
682  block=>ps(igrid)
683  level=node(plevel_,igrid)
684 
685  call process1_gridmf(method(level),igrid,qdt,ixg^ll,idim^lim,qtc,&
686  psa(igrid)%w,qt,psb(igrid)%w,pso(igrid)%w)
687  end do
688 
689  ! opedit: Send flux for all grids, expects sends for all
690  ! nsend_fc(^D), set in connectivity.t.
691 
692  if (fix_conserve_at_step) then
693  call recvflux(idim^lim)
694  call sendflux(idim^lim)
695  call fix_conserve(psb,idim^lim,mag(1),ndir)
696  end if
697  ! point bc mpi datatype to partial type for magnetic field
704  ! update B in ghost cells
705  call getbc(qt+qdt,qdt,psb,mag(1)-1,ndir)
706  ! calculate magnetofrictional velocity
707  call mf_velocity_update(qdt)
708  ! point bc mpi datatype to partial type for velocity field
715  ! update magnetofrictional velocity in ghost cells
716  call getbc(qt+qdt,qdt,psb,mom(1)-1,ndir)
717 
718  end subroutine advect1mf
719 
720  subroutine process1_gridmf(method,igrid,qdt,ixG^L,idim^LIM,qtC,wCT,qt,w,wold)
721  ! This subroutine is equivalent to VAC's `advect1' for one grid
723  use mod_fix_conserve
724 
725  character(len=*), intent(in) :: method
726  integer, intent(in) :: igrid, ixG^L, idim^LIM
727  double precision, intent(in) :: qdt, qtC, qt
728  double precision :: wCT(ixg^s,1:nw), w(ixg^s,1:nw), wold(ixg^s,1:nw)
729  double precision :: dx^D, fC(ixg^s,1:ndir,1:ndim)
730  integer :: ixO^L
731 
732  dx^d=rnode(rpdx^d_,igrid);
733  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
734  saveigrid=igrid
735  fc=0.d0
736 
739 
740  ixo^l=ixg^l^lsubnghostcells;
741  select case (method)
742  case ("cd4")
743  !================================
744  ! 4th order central difference
745  !================================
746  call centdiff4mf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,wold,fc,dx^d,ps(igrid)%x)
747  case ("tvdlf")
748  !================================
749  ! TVDLF
750  !================================
751  call tvdlfmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,wold,fc,dx^d,ps(igrid)%x)
752  case ('hancock')
753  ! hancock predict (first) step for twostep tvdlf and tvdmu scheme
754  call hancockmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,dx^d,ps(igrid)%x)
755  case ("fd")
756  !================================
757  ! finite difference
758  !================================
759  call fdmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,wold,fc,dx^d,ps(igrid)%x)
760  case default
761  if(mype==0) write(unitterm,*)'Error in advect1_gridmf:',method,' is not there!'
762  call mpistop("The scheme is not implemented.")
763  end select
764 
765  if (fix_conserve_at_step) then
766  call storeflux(igrid,fc,idim^lim,ndir)
767  end if
768 
769  end subroutine process1_gridmf
770 
771  subroutine upwindlrmf(ixI^L,ixL^L,ixR^L,idim,w,wCT,wLC,wRC,x)
772  ! Determine the upwinded wLC(ixL) and wRC(ixR) from w.
773  ! the wCT is only used when PPM is exploited.
775  use mod_limiter
776 
777  integer, intent(in) :: ixI^L, ixL^L, ixR^L, idim
778  double precision, dimension(ixI^S,1:nw) :: w, wCT
779  double precision, dimension(ixI^S,1:nw) :: wLC, wRC
780  double precision, dimension(ixI^S,1:ndim) :: x
781 
782  integer :: jxR^L, ixC^L, jxC^L, iw
783  double precision :: ldw(ixi^s), rdw(ixi^s), dwC(ixi^s)
784 
785  if (typelimiter == limiter_mp5) then
786  call mp5limiter(ixi^l,ixl^l,idim,w,wlc,wrc)
787  else if (typelimiter == limiter_ppm) then
788  call ppmlimiter(ixi^l,ixm^ll,idim,w,wct,wlc,wrc)
789  else
790  jxr^l=ixr^l+kr(idim,^d);
791  ixcmax^d=jxrmax^d; ixcmin^d=ixlmin^d-kr(idim,^d);
792  jxc^l=ixc^l+kr(idim,^d);
793 
794  do iw=1,nwflux
795  if (loglimit(iw)) then
796  w(ixcmin^d:jxcmax^d,iw)=dlog10(w(ixcmin^d:jxcmax^d,iw))
797  wlc(ixl^s,iw)=dlog10(wlc(ixl^s,iw))
798  wrc(ixr^s,iw)=dlog10(wrc(ixr^s,iw))
799  end if
800 
801  dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
802 
803  ! limit flux from left and/or right
804  call dwlimiter2(dwc,ixi^l,ixc^l,idim,typelimiter,ldw,rdw)
805  wlc(ixl^s,iw)=wlc(ixl^s,iw)+half*ldw(ixl^s)
806  wrc(ixr^s,iw)=wrc(ixr^s,iw)-half*rdw(jxr^s)
807 
808  if (loglimit(iw)) then
809  w(ixcmin^d:jxcmax^d,iw)=10.0d0**w(ixcmin^d:jxcmax^d,iw)
810  wlc(ixl^s,iw)=10.0d0**wlc(ixl^s,iw)
811  wrc(ixr^s,iw)=10.0d0**wrc(ixr^s,iw)
812  end if
813  end do
814 
815  endif
816 
817  end subroutine upwindlrmf
818 
819  subroutine getfluxmf(w,x,ixI^L,ixO^L,idir,idim,f)
820  ! Calculate lux f_idim[idir] within ixO^L.
822 
823  integer, intent(in) :: ixI^L, ixO^L, idir, idim
824  double precision, intent(in) :: w(ixi^s,nw)
825  double precision, intent(in) :: x(ixi^s,1:ndim)
826  double precision,intent(out) :: f(ixi^s)
827 
828  ! compute flux of magnetic field
829  ! f_i[b_k]=v_i*b_k-v_k*b_i
830  if (idim==idir) then
831  f(ixo^s)=zero
832  else
833  f(ixo^s)=w(ixo^s,mom(idim))*w(ixo^s,mag(idir))-w(ixo^s,mag(idim))*w(ixo^s,mom(idir))
834  if (b0field) then
835  f(ixo^s)=f(ixo^s)&
836  +w(ixo^s,mom(idim))*block%B0(ixo^s,idir,idim)&
837  -w(ixo^s,mom(idir))*block%B0(ixo^s,idim,idim)
838  end if
839  end if
840 
841  end subroutine getfluxmf
842 
843  subroutine tvdlfmf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,wnew,wold,fC,dx^D,x)
844  ! method=='tvdlf' --> 2nd order TVD-Lax-Friedrich scheme.
845  ! method=='tvdlf1' --> 1st order TVD-Lax-Friedrich scheme.
847 
848  double precision, intent(in) :: qdt, qtC, qt, dx^D
849  integer, intent(in) :: ixI^L, ixO^L, idim^LIM
850  double precision, dimension(ixI^S,1:ndim), intent(in) :: x
851  double precision, dimension(ixI^S,1:nw) :: wCT, wnew, wold
852  double precision, dimension(ixI^S,1:ndir,1:ndim) :: fC
853 
854  double precision, dimension(ixI^S,1:nw) :: wLC, wRC, wmean
855  double precision, dimension(ixI^S) :: fLC, fRC
856  double precision, dimension(ixI^S) :: cmaxC
857  double precision :: dxinv(1:ndim), inv_volume(ixo^s)
858  integer :: idims, idir, ix^L, hxO^L, ixC^L, ixCR^L, jxC^L, kxC^L, kxR^L
859 
860  ! The flux calculation contracts by one in the idim direction it is applied.
861  ! The limiter contracts the same directions by one more, so expand ixO by 2.
862  ix^l=ixo^l;
863  do idims= idim^lim
864  ix^l=ix^l^ladd2*kr(idims,^d);
865  end do
866  if (ixi^l^ltix^l|.or.|.or.) &
867  call mpistop("Error in tvdlfmf: Nonconforming input limits")
868 
869  ^d&dxinv(^d)=-qdt/dx^d;
870  fc=0.d0
871  do idims= idim^lim
872  block%iw0=idims
873 
874  hxo^l=ixo^l-kr(idims,^d);
875  ! ixC is centered index in the idim direction from ixOmin-1/2 to ixOmax+1/2
876  ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
877  ! Calculate wRC=uR_{j+1/2} and wLC=uL_j+1/2
878  jxc^l=ixc^l+kr(idims,^d);
879  kxcmin^d=iximin^d; kxcmax^d=iximax^d-kr(idims,^d);
880  kxr^l=kxc^l+kr(idims,^d);
881  ixcr^l=ixc^l;
882 
883  wrc(kxc^s,1:nwflux)=wct(kxr^s,1:nwflux)
884  wlc(kxc^s,1:nwflux)=wct(kxc^s,1:nwflux)
885 
886  call upwindlrmf(ixi^l,ixcr^l,ixcr^l,idims,wct,wct,wlc,wrc,x)
887 
888  ! For the high order Lax-Friedrich TVDLF scheme the limiter is based on
889  ! the maximum eigenvalue, it is calculated in advance.
890  ! determine mean state and store in wLC
891  wmean=0.5d0*(wlc+wrc)
892  call getcmaxfff(wmean,ixg^ll,ixc^l,idims,cmaxc)
893 
894  ! Calculate fLC=f(uL_j+1/2) and fRC=f(uR_j+1/2) for each idir
895  do idir=1,ndir
896  call getfluxmf(wlc,x,ixg^ll,ixc^l,idir,idims,flc)
897  call getfluxmf(wrc,x,ixg^ll,ixc^l,idir,idims,frc)
898  ! To save memory we use fLC to store (F_L+F_R)/2=half*(fLC+fRC)
899  flc(ixc^s)=half*(flc(ixc^s)+frc(ixc^s))
900 
901  ! Add TVDLF dissipation to the flux
902  if (idir==idims) then
903  flc(ixc^s)=flc(ixc^s)-mf_tvdlfeps*tvdlfeps*cmaxc(ixc^s)*half*(wrc(ixc^s,mag(idir))-wlc(ixc^s,mag(idir)))
904  end if
905  if (slab) then
906  fc(ixc^s,idir,idims)=flc(ixc^s)
907  else
908  fc(ixc^s,idir,idims)=block%surfaceC(ixc^s,idims)*flc(ixc^s)
909  end if
910 
911  end do ! Next idir
912  end do ! Next idims
913 
914  !Now update the state:
915  do idims= idim^lim
916  hxo^l=ixo^l-kr(idims,^d);
917  ! Multiply the fluxes by -dt/dx since Flux fixing expects this
918  if (slab) then
919  fc(ixi^s,:,idims)=dxinv(idims)*fc(ixi^s,:,idims)
920  wnew(ixo^s,mag(:))=wnew(ixo^s,mag(:)) &
921  + (fc(ixo^s,:,idims)-fc(hxo^s,:,idims))
922  else
923  inv_volume = 1.0d0/block%dvolume(ixo^s)
924  fc(ixi^s,:,idims)=-qdt*fc(ixi^s,:,idims)
925 
926  do idir = 1, ndir
927  wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir)) + (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims)) * &
928  inv_volume
929  end do
930  end if
931 
932  end do ! Next idims
933 
934  if (.not.slab) call addgeometrymf(qdt,ixi^l,ixo^l,wct,wnew,x)
935  call divbclean(qdt,ixi^l,ixo^l,wct,wnew,x)
936 
937  end subroutine tvdlfmf
938 
939  subroutine hancockmf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,wnew,dx^D,x)
940  ! The non-conservative Hancock predictor for TVDLFmf
941  ! on entry:
942  ! input available on ixI^L=ixG^L asks for output on ixO^L=ixG^L^LSUBnghostcells
943  ! one entry: (predictor): wCT -- w_n wnew -- w_n qdt=dt/2
944  ! on exit : (predictor): wCT -- w_n wnew -- w_n+1/2
946 
947  integer, intent(in) :: ixI^L, ixO^L, idim^LIM
948  double precision, intent(in) :: qdt, qtC, qt, dx^D, x(ixi^s,1:ndim)
949  double precision, intent(inout) :: wCT(ixi^s,1:nw), wnew(ixi^s,1:nw)
950 
951  double precision, dimension(ixI^S,1:nw) :: wLC, wRC
952  double precision, dimension(ixI^S) :: fLC, fRC
953  double precision :: dxinv(1:ndim)
954  integer :: idims, idir, ix^L, hxO^L, ixtest^L
955 
956  ! Expand limits in each idims direction in which fluxes are added
957  ix^l=ixo^l;
958  do idims= idim^lim
959  ix^l=ix^l^laddkr(idims,^d);
960  end do
961  if (ixi^l^ltix^l|.or.|.or.) &
962  call mpistop("Error in Hancockmf: Nonconforming input limits")
963 
964  ^d&dxinv(^d)=-qdt/dx^d;
965  do idims= idim^lim
966  block%iw0=idims
967  ! Calculate w_j+g_j/2 and w_j-g_j/2
968  ! First copy all variables, then upwind wLC and wRC.
969  ! wLC is to the left of ixO, wRC is to the right of wCT.
970  hxo^l=ixo^l-kr(idims,^d);
971 
972  wrc(hxo^s,1:nwflux)=wct(ixo^s,1:nwflux)
973  wlc(ixo^s,1:nwflux)=wct(ixo^s,1:nwflux)
974 
975  call upwindlrmf(ixi^l,ixo^l,hxo^l,idims,wct,wct,wlc,wrc,x)
976 
977  ! Advect mag(idir)
978  do idir=1,ndir
979  ! Calculate the fLC and fRC fluxes
980  call getfluxmf(wrc,x,ixi^l,hxo^l,idir,idims,frc)
981  call getfluxmf(wlc,x,ixi^l,ixo^l,idir,idims,flc)
982 
983  if (slab) then
984  wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))+dxinv(idims)* &
985  (flc(ixo^s)-frc(hxo^s))
986  else
987  wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))-qdt/block%dvolume(ixo^s) &
988  *(block%surfaceC(ixo^s,idims)*flc(ixo^s) &
989  -block%surfaceC(hxo^s,idims)*frc(hxo^s))
990  end if
991  end do
992  end do ! next idims
993 
994  if (.not.slab) call addgeometrymf(qdt,ixi^l,ixo^l,wct,wnew,x)
995 
996  end subroutine hancockmf
997 
998  subroutine fdmf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,wnew,wold,fC,dx^D,x)
1000  double precision, intent(in) :: qdt, qtC, qt, dx^D
1001  integer, intent(in) :: ixI^L, ixO^L, idim^LIM
1002  double precision, dimension(ixI^S,1:ndim), intent(in) :: x
1003 
1004  double precision, dimension(ixI^S,1:nw), intent(inout) :: wCT, wnew, wold
1005  double precision, dimension(ixI^S,1:ndir,1:ndim), intent(out) :: fC
1006 
1007  double precision, dimension(ixI^S) :: fCT
1008  double precision, dimension(ixI^S,1:nw) :: fm, fp, fmR, fpL
1009  double precision, dimension(ixI^S) :: v
1010  double precision :: dxinv(1:ndim)
1011  integer :: idims, idir, ixC^L, ix^L, hxO^L, ixCR^L
1012 
1013  ^d&dxinv(^d)=-qdt/dx^d;
1014  do idims= idim^lim
1015 
1016  ! Get fluxes for the whole grid (mesh+nghostcells)
1017  {^d& ixcmin^d = ixomin^d - nghostcells * kr(idims,^d)\}
1018  {^d& ixcmax^d = ixomax^d + nghostcells * kr(idims,^d)\}
1019 
1020  hxo^l=ixo^l-kr(idims,^d);
1021  ! ix is centered index in the idim direction from ixOmin-1/2 to ixOmax+1/2
1022  ixmax^d=ixomax^d; ixmin^d=hxomin^d;
1023  ixcr^l=ixc^l;
1024 
1025  do idir=1,ndir
1026  call getfluxmf(wct,x,ixg^ll,ixcr^l,idir,idims,fct)
1027  ! Lax-Friedrich splitting:
1028  fp(ixcr^s,mag(idir)) = half * (fct(ixcr^s) + mf_tvdlfeps * tvdlfeps * cmax_global * wct(ixcr^s,mag(idir)))
1029  fm(ixcr^s,mag(idir)) = half * (fct(ixcr^s) - mf_tvdlfeps * tvdlfeps * cmax_global * wct(ixcr^s,mag(idir)))
1030  end do ! iw loop
1031 
1032  ! now do the reconstruction of fp and fm:
1033  call reconstructlmf(ixi^l,ix^l,idims,fp,fpl)
1034  call reconstructrmf(ixi^l,ix^l,idims,fm,fmr)
1035 
1036  do idir=1,ndir
1037  if (slab) then
1038  fc(ix^s,idir,idims) = dxinv(idims) * (fpl(ix^s,mag(idir)) + fmr(ix^s,mag(idir)))
1039  wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))+ &
1040  (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))
1041  else
1042  fc(ix^s,idir,idims)=-qdt*block%surfaceC(ix^s,idims) * (fpl(ix^s,mag(idir)) + fmr(ix^s,mag(idir)))
1043  wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))+ &
1044  (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))/block%dvolume(ixo^s)
1045  end if
1046  end do ! iw loop
1047 
1048  end do !idims loop
1049 
1050  if (.not.slab) call addgeometrymf(qdt,ixi^l,ixo^l,wct,wnew,x)
1051  call divbclean(qdt,ixi^l,ixo^l,wct,wnew,x)
1052 
1053  end subroutine fdmf
1054 
1055  subroutine reconstructlmf(ixI^L,iL^L,idims,w,wLC)
1057  use mod_limiter
1058 
1059  integer, intent(in) :: ixI^L, iL^L, idims
1060  double precision, intent(in) :: w(ixi^s,1:nw)
1061 
1062  double precision, intent(out) :: wLC(ixi^s,1:nw)
1063 
1064  double precision :: ldw(ixi^s), dwC(ixi^s)
1065  integer :: jxR^L, ixC^L, jxC^L, kxC^L, iw
1066 
1067  select case (typelimiter)
1068  case (limiter_mp5)
1069  call mp5limiterl(ixi^l,il^l,idims,w,wlc)
1070  case default
1071 
1072  kxcmin^d=iximin^d; kxcmax^d=iximax^d-kr(idims,^d);
1073 
1074  wlc(kxc^s,1:nwflux) = w(kxc^s,1:nwflux)
1075 
1076  jxr^l=il^l+kr(idims,^d);
1077 
1078  ixcmax^d=jxrmax^d; ixcmin^d=ilmin^d-kr(idims,^d);
1079  jxc^l=ixc^l+kr(idims,^d);
1080 
1081  do iw=1,nwflux
1082  dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
1083 
1084  call dwlimiter2(dwc,ixi^l,ixc^l,idims,typelimiter,ldw=ldw)
1085 
1086  wlc(il^s,iw)=wlc(il^s,iw)+half*ldw(il^s)
1087  end do
1088  end select
1089 
1090  end subroutine reconstructlmf
1091 
1092  subroutine reconstructrmf(ixI^L,iL^L,idims,w,wRC)
1094  use mod_limiter
1095 
1096  integer, intent(in) :: ixI^L, iL^L, idims
1097  double precision, intent(in) :: w(ixi^s,1:nw)
1098 
1099  double precision, intent(out) :: wRC(ixi^s,1:nw)
1100 
1101  double precision :: rdw(ixi^s), dwC(ixi^s)
1102  integer :: jxR^L, ixC^L, jxC^L, kxC^L, kxR^L, iw
1103 
1104  select case (typelimiter)
1105  case (limiter_mp5)
1106  call mp5limiterr(ixi^l,il^l,idims,w,wrc)
1107  case default
1108 
1109  kxcmin^d=iximin^d; kxcmax^d=iximax^d-kr(idims,^d);
1110  kxr^l=kxc^l+kr(idims,^d);
1111 
1112  wrc(kxc^s,1:nwflux)=w(kxr^s,1:nwflux)
1113 
1114  jxr^l=il^l+kr(idims,^d);
1115  ixcmax^d=jxrmax^d; ixcmin^d=ilmin^d-kr(idims,^d);
1116  jxc^l=ixc^l+kr(idims,^d);
1117 
1118  do iw=1,nwflux
1119  dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
1120  call dwlimiter2(dwc,ixi^l,ixc^l,idims,typelimiter,rdw=rdw)
1121 
1122  wrc(il^s,iw)=wrc(il^s,iw)-half*rdw(jxr^s)
1123  end do
1124  end select
1125 
1126  end subroutine reconstructrmf
1127 
1128  subroutine centdiff4mf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,w,wold,fC,dx^D,x)
1129  ! Advance the flow variables from global_time to global_time+qdt within ixO^L by
1130  ! fourth order centered differencing in space
1131  ! for the dw/dt+dF_i(w)/dx_i=S type equation.
1132  ! wCT contains the time centered variables at time qtC for flux and source.
1133  ! w is the old value at qt on input and the new value at qt+qdt on output.
1135 
1136  integer, intent(in) :: ixI^L, ixO^L, idim^LIM
1137  double precision, intent(in) :: qdt, qtC, qt, dx^D
1138  double precision :: wCT(ixi^s,1:nw), w(ixi^s,1:nw), wold(ixi^s,1:nw)
1139  double precision, intent(in) :: x(ixi^s,1:ndim)
1140  double precision :: fC(ixi^s,1:ndir,1:ndim)
1141 
1142  double precision :: v(ixi^s,ndim), f(ixi^s)
1143  double precision, dimension(ixI^S,1:nw) :: wLC, wRC
1144  double precision, dimension(ixI^S) :: vLC, vRC,cmaxLC,cmaxRC
1145  double precision :: dxinv(1:ndim)
1146  integer :: idims, idir, idirmin,ix^D
1147  integer :: ix^L, hxO^L, ixC^L, jxC^L, hxC^L, kxC^L, kkxC^L, kkxR^L
1148 
1149  ! two extra layers are needed in each direction for which fluxes are added.
1150  ix^l=ixo^l;
1151  do idims= idim^lim
1152  ix^l=ix^l^ladd2*kr(idims,^d);
1153  end do
1154 
1155  if (ixi^l^ltix^l|.or.|.or.) then
1156  call mpistop("Error in evolve_CentDiff4: Non-conforming input limits")
1157  end if
1158  ^d&dxinv(^d)=-qdt/dx^d;
1159 
1160  ! Add fluxes to w
1161  do idims= idim^lim
1162  block%iw0=idims
1163  ix^l=ixo^l^ladd2*kr(idims,^d);
1164  hxo^l=ixo^l-kr(idims,^d);
1165 
1166  ixcmin^d=hxomin^d; ixcmax^d=ixomax^d;
1167  hxc^l=ixc^l-kr(idims,^d);
1168  jxc^l=ixc^l+kr(idims,^d);
1169  kxc^l=ixc^l+2*kr(idims,^d);
1170 
1171  kkxcmin^d=iximin^d; kkxcmax^d=iximax^d-kr(idims,^d);
1172  kkxr^l=kkxc^l+kr(idims,^d);
1173  wrc(kkxc^s,1:nwflux)=wct(kkxr^s,1:nwflux)
1174  wlc(kkxc^s,1:nwflux)=wct(kkxc^s,1:nwflux)
1175 
1176  call upwindlrmf(ixi^l,ixc^l,ixc^l,idims,wct,wct,wlc,wrc,x)
1177 
1178  ! Calculate velocities from upwinded values
1179  call getcmaxfff(wlc,ixg^ll,ixc^l,idims,cmaxlc)
1180  call getcmaxfff(wrc,ixg^ll,ixc^l,idims,cmaxrc)
1181  ! now take the maximum of left and right states
1182  vlc(ixc^s)=max(cmaxrc(ixc^s),cmaxlc(ixc^s))
1183 
1184  do idir=1,ndir
1185  ! Get non-transported flux
1186  call getfluxmf(wct,x,ixi^l,ix^l,idir,idims,f)
1187  ! Center flux to interface
1188  ! f_i+1/2= (-f_(i+2) +7 f_(i+1) + 7 f_i - f_(i-1))/12
1189  fc(ixc^s,idir,idims)=(-f(kxc^s)+7.0d0*(f(jxc^s)+f(ixc^s))-f(hxc^s))/12.0d0
1190  ! add rempel dissipative flux, only second order version for now
1191  ! one could gradually reduce the dissipative flux to improve solutions
1192  ! for computing steady states (Keppens et al. 2012, JCP)
1193  fc(ixc^s,idir,idims)=fc(ixc^s,idir,idims)-mf_tvdlfeps*tvdlfeps*half*vlc(ixc^s) &
1194  *(wrc(ixc^s,mag(idir))-wlc(ixc^s,mag(idir)))
1195 
1196  if (slab) then
1197  fc(ixc^s,idir,idims)=dxinv(idims)*fc(ixc^s,idir,idims)
1198  ! result: f_(i+1/2)-f_(i-1/2) = [-f_(i+2)+8(f_(i+1)+f_(i-1))-f_(i-2)]/12
1199  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+(fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))
1200  else
1201  fc(ixc^s,idir,idims)=-qdt*block%surfaceC(ixc^s,idims)*fc(ixc^s,idir,idims)
1202  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+ &
1203  (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))/block%dvolume(ixo^s)
1204  end if
1205  end do !next idir
1206  end do !next idims
1207 
1208  if (.not.slab) call addgeometrymf(qdt,ixi^l,ixo^l,wct,w,x)
1209  call divbclean(qdt,ixi^l,ixo^l,wct,w,x)
1210 
1211  end subroutine centdiff4mf
1212 
1213  subroutine getdtfff_courant(w,x,ixI^L,ixO^L,dtnew)
1214  ! compute CFL limited dt (for variable time stepping)
1216 
1217  integer, intent(in) :: ixI^L, ixO^L
1218  double precision, intent(in) :: x(ixi^s,1:ndim)
1219  double precision, intent(inout) :: w(ixi^s,1:nw), dtnew
1220 
1221  double precision :: courantmax, dxinv(1:ndim)
1222  double precision :: cmax(ixi^s),tmp(ixi^s),alfven(ixi^s)
1223  integer :: idims
1224 
1225  dtnew=bigdouble
1226  courantmax=0.d0
1227  ^d&dxinv(^d)=1.d0/dxlevel(^d);
1228 
1229  do idims=1,ndim
1230  call getcmaxfff(w,ixi^l,ixo^l,idims,cmax)
1231  cmax_mype = max(cmax_mype,maxval(cmax(ixo^s)))
1232  if (.not.slab) then
1233  tmp(ixo^s)=cmax(ixo^s)/block%dx(ixo^s,idims)
1234  courantmax=max(courantmax,maxval(tmp(ixo^s)))
1235  else
1236  tmp(ixo^s)=cmax(ixo^s)*dxinv(idims)
1237  courantmax=max(courantmax,maxval(tmp(ixo^s)))
1238  end if
1239  end do
1240  ! courantmax='max( c/dx)'
1241  if (courantmax>smalldouble) dtnew=min(dtnew,mf_cc/courantmax)
1242 
1243  end subroutine getdtfff_courant
1244 
1245  subroutine getcmaxfff(w,ixI^L,ixO^L,idims,cmax)
1247 
1248  logical :: new_cmax,needcmin
1249  integer, intent(in) :: ixI^L, ixO^L, idims
1250  double precision, intent(in) :: w(ixi^s,1:nw)
1251  double precision, intent(out) :: cmax(ixi^s)
1252 
1253  ! calculate alfven speed
1254  if(b0field) then
1255  cmax(ixo^s)=sqrt(sum((w(ixo^s,mag(:))+block%b0(ixo^s,:,0))**2,dim=ndim+1)/w(ixo^s,rho_))
1256  else
1257  cmax(ixo^s)=sqrt(sum(w(ixo^s,mag(:))**2,dim=ndim+1)/w(ixo^s,rho_))
1258  endif
1259  cmax(ixo^s)=cmax(ixo^s)+abs(w(ixo^s,mom(idims)))
1260 
1261  end subroutine getcmaxfff
1262 
1263  !> Clean divergence of magnetic field by Janhunen's and Linde's source terms
1264  subroutine divbclean(qdt,ixI^L,ixO^L,wCT,w,x)
1266  use mod_geometry
1267 
1268  integer, intent(in) :: ixI^L, ixO^L
1269  double precision, intent(in) :: x(ixi^s,1:ndim),wCT(ixi^s,1:nw),qdt
1270  double precision, intent(inout) :: w(ixi^s,1:nw)
1271  integer :: idims, ix^L, ixp^L, i^D, iside
1272  double precision :: divb(ixi^s),graddivb(ixi^s),bdivb(ixi^s,1:ndir)
1273 
1274  ! Calculate div B
1275  ix^l=ixo^l^ladd1;
1276  call get_divb(wct,ixi^l,ix^l,divb)
1277 
1278  ixp^l=ixo^l;
1279 
1280  ! Add Linde's diffusive terms
1281  do idims=1,ndim
1282  ! Calculate grad_idim(divb)
1283  call gradient(divb,ixi^l,ixp^l,idims,graddivb)
1284 
1285  ! Multiply by Linde's eta*dt = divbdiff*(c_max*dx)*dt = divbdiff*dx**2
1286  if (slab) then
1287  graddivb(ixp^s)=graddivb(ixp^s)*mf_cdivb/(^d&1.0d0/dxlevel(^d)**2+)
1288  else
1289  graddivb(ixp^s)=graddivb(ixp^s)*mf_cdivb &
1290  /(^d&1.0d0/block%dx(ixp^s,^d)**2+)
1291  end if
1292  ! B_idim += eta*grad_idim(divb)
1293  ! with Janhunen's term
1294  !w(ixp^S,mag(idims))=w(ixp^S,mag(idims))+&
1295  ! graddivb(ixp^S)-qdt*w(ixp^S,mom(idims))*divb(ixp^S)
1296  ! without Janjunen's term
1297  w(ixp^s,mag(idims))=w(ixp^s,mag(idims))+&
1298  graddivb(ixp^s)
1299  end do
1300 
1301  end subroutine divbclean
1302 
1303  subroutine addgeometrymf(qdt,ixI^L,ixO^L,wCT,w,x)
1304  ! Add geometrical source terms to w
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 (typeaxial)
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, 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.
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:255
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:430
integer, dimension(-1:2^d &,-1:1^d &), target type_recv_srl_p2
integer, parameter limiter_mp5
Definition: mod_limiter.t:23
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
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, 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)
character(len=std_len) typeaxial
subroutine resettree
reset AMR and (de)allocate boundary flux storage at level changes
Definition: amrgrid.t:44
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:22
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
procedure(sub_convert), pointer phys_to_conserved
Definition: mod_physics.t:36
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:485
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)
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:359
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, public storeflux(igrid, fC, idimLIM, nwfluxin)
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:160
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.
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:37
integer, dimension(:,:), allocatable node
integer, dimension(-1:1, 0:3) l
logical slab_stretched
Stretched Cartesian geometry or not.
integer, dimension(-1:1^d &,-1:1^d &), target type_send_r_f
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:80
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)