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