MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_mf_phys.t
Go to the documentation of this file.
1 !> Magnetofriction module
2 module mod_mf_phys
3  use mod_global_parameters, only: std_len
5  use mod_comm_lib, only: mpistop
6 
7 
8  implicit none
9  private
10 
11  !> viscosity coefficient in s cm^-2 for solar corona (Cheung 2012 ApJ)
12  double precision, public :: mf_nu = 1.d-15
13 
14  !> maximal limit of magnetofrictional velocity in cm s^-1 (Pomoell 2019 A&A)
15  double precision, public :: mf_vmax = 3.d6
16 
17  !> decay scale of frictional velocity near boundaries
18  double precision, public :: mf_decay_scale(2*^nd)=0.d0
19 
20  !> GLM-MHD parameter: ratio of the diffusive and advective time scales for div b
21  !> taking values within [0, 1]
22  double precision, public :: mf_glm_alpha = 0.5d0
23 
24  !> The resistivity
25  double precision, public :: mf_eta = 0.0d0
26 
27  !> The hyper-resistivity
28  double precision, public :: mf_eta_hyper = 0.0d0
29 
30  !> Coefficient of diffusive divB cleaning
31  double precision :: divbdiff = 0.8d0
32 
33  !> Helium abundance over Hydrogen
34  double precision, public, protected :: he_abundance=0.1d0
35 
36  !> Method type in a integer for good performance
37  integer :: type_divb
38 
39  !> Indices of the momentum density
40  integer, allocatable, public, protected :: mom(:)
41 
42  !> Indices of the GLM psi
43  integer, public, protected :: psi_
44 
45  !> To skip * layer of ghost cells during divB=0 fix for boundary
46  integer, public, protected :: boundary_divbfix_skip(2*^nd)=0
47 
48  ! DivB cleaning methods
49  integer, parameter :: divb_none = 0
50  integer, parameter :: divb_multigrid = -1
51  integer, parameter :: divb_glm = 1
52  integer, parameter :: divb_powel = 2
53  integer, parameter :: divb_janhunen = 3
54  integer, parameter :: divb_linde = 4
55  integer, parameter :: divb_lindejanhunen = 5
56  integer, parameter :: divb_lindepowel = 6
57  integer, parameter :: divb_lindeglm = 7
58  integer, parameter :: divb_ct = 8
59 
60  !> Whether particles module is added
61  logical, public, protected :: mf_particles = .false.
62 
63  !> Whether GLM-MHD is used
64  logical, public, protected :: mf_glm = .false.
65 
66  !> Whether divB cleaning sources are added splitting from fluid solver
67  logical, public, protected :: source_split_divb = .false.
68 
69  !> MHD fourth order
70  logical, public, protected :: mf_4th_order = .false.
71 
72  !> set to true if need to record electric field on cell edges
73  logical, public, protected :: mf_record_electric_field = .false.
74 
75  !> Whether divB is computed with a fourth order approximation
76  logical, public, protected :: mf_divb_4thorder = .false.
77 
78  !> To control divB=0 fix for boundary
79  logical, public, protected :: boundary_divbfix(2*^nd)=.true.
80 
81  !> Use a compact way to add resistivity
82  logical :: compactres = .false.
83 
84  !> Add divB wave in Roe solver
85  logical, public :: divbwave = .true.
86 
87  !> clean divb in the initial condition
88  logical, public, protected :: clean_initial_divb=.false.
89 
90  !> Method type to clean divergence of B
91  character(len=std_len), public, protected :: typedivbfix = 'ct'
92 
93  !> Method type of constrained transport
94  character(len=std_len), public, protected :: type_ct = 'average'
95 
96  !> Update all equations due to divB cleaning
97  character(len=std_len) :: typedivbdiff = 'all'
98 
99 
100  ! Public methods
101  public :: mf_phys_init
102  public :: mf_to_conserved
103  public :: mf_to_primitive
104  public :: mf_face_to_center
105  public :: get_divb
106  public :: get_current
107  public :: get_normalized_divb
108  public :: b_from_vector_potential
109  public :: record_force_free_metrics
110  {^nooned
111  public :: mf_clean_divb_multigrid
112  }
113 
114 contains
115 
116  !> Read this module's parameters from a file
117  subroutine mf_read_params(files)
119  use mod_particles, only: particles_eta
120  character(len=*), intent(in) :: files(:)
121  integer :: n
122 
123  namelist /mf_list/ mf_nu, mf_vmax, mf_decay_scale, &
125  particles_eta, mf_record_electric_field,&
127  typedivbdiff, type_ct, compactres, divbwave, he_abundance, si_unit, &
130 
131  do n = 1, size(files)
132  open(unitpar, file=trim(files(n)), status="old")
133  read(unitpar, mf_list, end=111)
134 111 close(unitpar)
135  end do
136 
137  end subroutine mf_read_params
138 
139  !> Write this module's parameters to a snapsoht
140  subroutine mf_write_info(fh)
142  integer, intent(in) :: fh
143  integer, parameter :: n_par = 1
144  double precision :: values(n_par)
145  character(len=name_len) :: names(n_par)
146  integer, dimension(MPI_STATUS_SIZE) :: st
147  integer :: er
148 
149  call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
150 
151  names(1) = "nu"
152  values(1) = mf_nu
153  call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
154  call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
155  end subroutine mf_write_info
156 
157  subroutine mf_phys_init()
159  use mod_physics
160  use mod_particles, only: particles_init, particles_eta, particles_etah
161  {^nooned
163  }
164 
165  integer :: itr, idir
166 
167  call mf_read_params(par_files)
168 
169  physics_type = "mf"
170  phys_energy = .false.
172 
173  if(ndim==1) typedivbfix='none'
174  select case (typedivbfix)
175  case ('none')
176  type_divb = divb_none
177  {^nooned
178  case ('multigrid')
179  type_divb = divb_multigrid
180  use_multigrid = .true.
181  mg%operator_type = mg_laplacian
183  }
184  case ('glm')
185  mf_glm = .true.
186  need_global_cmax = .true.
187  type_divb = divb_glm
188  case ('powel', 'powell')
189  type_divb = divb_powel
190  case ('janhunen')
191  type_divb = divb_janhunen
192  case ('linde')
193  type_divb = divb_linde
194  case ('lindejanhunen')
195  type_divb = divb_lindejanhunen
196  case ('lindepowel')
197  type_divb = divb_lindepowel
198  case ('lindeglm')
199  mf_glm = .true.
200  need_global_cmax = .true.
201  type_divb = divb_lindeglm
202  case ('ct')
203  type_divb = divb_ct
204  stagger_grid = .true.
205  case default
206  call mpistop('Unknown divB fix')
207  end select
208 
209  ! Whether diagonal ghost cells are required for the physics
210  if(type_divb < divb_linde) phys_req_diagonal = .false.
211 
212  ! set magnetic field as flux variables
213  allocate(mag(ndir))
214  mag(:) = var_set_bfield(ndir)
215 
216  ! start with magnetic field and skip velocity when update ghostcells
217  iwstart=mag(1)
218  allocate(start_indices(number_species),stop_indices(number_species))
219  ! set the index of the first flux variable for species 1
220  start_indices(1)=mag(1)
221 
222  if (mf_glm) then
223  psi_ = var_set_fluxvar('psi', 'psi', need_bc=.false.)
224  else
225  psi_ = -1
226  end if
227 
228  ! set velocity field as flux variables
229  allocate(mom(ndir))
230  mom(:) = var_set_momentum(ndir)
231 
232  ! Number of variables need reconstruction in w
233  nw_recon=nwflux
234 
235  ! set number of variables which need update ghostcells
236  nwgc=nwflux-ndir
237 
238  ! remove velocity field from flux variables
239  nwflux=nwflux-ndir
240 
241  ! set the index of the last flux variable for species 1
242  stop_indices(1)=nwflux
243 
244  ! determine number of stagger variables
245  nws=ndim
246 
247  nvector = 2 ! No. vector vars
248  allocate(iw_vector(nvector))
249  iw_vector(1) = mom(1) - 1 ! TODO: why like this?
250  iw_vector(2) = mag(1) - 1 ! TODO: why like this?
251 
252  ! Check whether custom flux types have been defined
253  if (.not. allocated(flux_type)) then
254  allocate(flux_type(ndir, nw))
256  else if (any(shape(flux_type) /= [ndir, nw])) then
257  call mpistop("phys_check error: flux_type has wrong shape")
258  end if
259  do idir=1,ndir
260  if(ndim>1) flux_type(idir,mag(idir))=flux_tvdlf
261  end do
262  if(mf_glm .and. ndim>1) flux_type(:,psi_)=flux_tvdlf
263 
275 
276  if(type_divb==divb_glm) then
278  end if
279 
280  ! pass to global variable to record electric field
282 
283  ! Initialize particles module
284  if(mf_particles) then
285  call particles_init()
286  phys_req_diagonal = .true.
287  ! never allow Hall effects in particles when doing magnetofrictional
288  particles_etah=0.0d0
289  if(mype==0) then
290  write(*,*) '*****Using particles: with mf_eta, mf_eta_hyper :', mf_eta, mf_eta_hyper
291  write(*,*) '*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
292  end if
293  end if
294 
295  ! if using ct stagger grid, boundary divb=0 is not done here
296  if(stagger_grid) then
300  else if(ndim>1) then
302  end if
303 
304  {^nooned
305  ! clean initial divb
307  }
308 
309  ! derive units from basic units
310  call mf_physical_units()
311 
312  end subroutine mf_phys_init
313 
314  subroutine mf_check_params
316 
317  end subroutine mf_check_params
318 
319  subroutine mf_physical_units()
321  double precision :: mp,kB,miu0,c_lightspeed
322  ! Derive scaling units
323  if(si_unit) then
324  mp=mp_si
325  kb=kb_si
326  miu0=miu0_si
327  c_lightspeed=c_si
328  else
329  mp=mp_cgs
330  kb=kb_cgs
331  miu0=4.d0*dpi
332  c_lightspeed=const_c
333  end if
334  if(unit_density/=1.d0) then
336  else
337  ! unit of numberdensity is independent by default
339  end if
340  if(unit_velocity/=1.d0) then
344  else if(unit_pressure/=1.d0) then
348  else if(unit_magneticfield/=1.d0) then
352  else
353  ! unit of temperature is independent by default
357  end if
358  if(unit_time/=1.d0) then
360  else
361  ! unit of length is independent by default
363  end if
364  ! Additional units needed for the particles
365  c_norm=c_lightspeed/unit_velocity
367  if (.not. si_unit) unit_charge = unit_charge*const_c
369 
370  ! get dimensionless mf nu
372 
373  ! get dimensionless maximal mf velocity limit
375 
376  end subroutine mf_physical_units
377 
378  !> Transform primitive variables into conservative ones
379  subroutine mf_to_conserved(ixI^L,ixO^L,w,x)
381  integer, intent(in) :: ixi^l, ixo^l
382  double precision, intent(inout) :: w(ixi^s, nw)
383  double precision, intent(in) :: x(ixi^s, 1:ndim)
384 
385  ! nothing to do for mf
386  end subroutine mf_to_conserved
387 
388  !> Transform conservative variables into primitive ones
389  subroutine mf_to_primitive(ixI^L,ixO^L,w,x)
391  integer, intent(in) :: ixi^l, ixo^l
392  double precision, intent(inout) :: w(ixi^s, nw)
393  double precision, intent(in) :: x(ixi^s, 1:ndim)
394 
395  ! nothing to do for mf
396  end subroutine mf_to_primitive
397 
398  !> Calculate cmax_idim=csound+abs(v_idim) within ixO^L
399  subroutine mf_get_cmax(w,x,ixI^L,ixO^L,idim,cmax)
401 
402  integer, intent(in) :: ixI^L, ixO^L, idim
403  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
404  double precision, intent(inout) :: cmax(ixI^S)
405 
406  cmax(ixo^s)=abs(w(ixo^s,mom(idim)))+one
407 
408  end subroutine mf_get_cmax
409 
410  !> Estimating bounds for the minimum and maximum signal velocities
411  subroutine mf_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
413  use mod_variables
414 
415  integer, intent(in) :: ixI^L, ixO^L, idim
416  double precision, intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
417  double precision, intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
418  double precision, intent(in) :: x(ixI^S,1:ndim)
419  double precision, intent(inout) :: cmax(ixI^S,1:number_species)
420  double precision, intent(inout), optional :: cmin(ixI^S,1:number_species)
421  double precision, intent(in) :: Hspeed(ixI^S,1:number_species)
422 
423  double precision, dimension(ixI^S) :: tmp1
424 
425  tmp1(ixo^s)=0.5d0*(wlc(ixo^s,mom(idim))+wrc(ixo^s,mom(idim)))
426  if(present(cmin)) then
427  cmax(ixo^s,1)=max(tmp1(ixo^s)+one,zero)
428  cmin(ixo^s,1)=min(tmp1(ixo^s)-one,zero)
429  else
430  cmax(ixo^s,1)=abs(tmp1(ixo^s))+one
431  end if
432 
433  end subroutine mf_get_cbounds
434 
435  !> prepare velocities for ct methods
436  subroutine mf_get_ct_velocity(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
438 
439  integer, intent(in) :: ixI^L, ixO^L, idim
440  double precision, intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
441  double precision, intent(in) :: cmax(ixI^S)
442  double precision, intent(in), optional :: cmin(ixI^S)
443  type(ct_velocity), intent(inout):: vcts
444 
445  integer :: idimE,idimN
446 
447  ! calculate velocities related to different UCT schemes
448  select case(type_ct)
449  case('average')
450  case('uct_contact')
451  if(.not.allocated(vcts%vnorm)) allocate(vcts%vnorm(ixi^s,1:ndim))
452  ! get average normal velocity at cell faces
453  vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,mom(idim))+wrp(ixo^s,mom(idim)))
454  case('uct_hll')
455  if(.not.allocated(vcts%vbarC)) then
456  allocate(vcts%vbarC(ixi^s,1:ndir,2),vcts%vbarLC(ixi^s,1:ndir,2),vcts%vbarRC(ixi^s,1:ndir,2))
457  allocate(vcts%cbarmin(ixi^s,1:ndim),vcts%cbarmax(ixi^s,1:ndim))
458  end if
459  ! Store magnitude of characteristics
460  if(present(cmin)) then
461  vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
462  vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
463  else
464  vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
465  vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
466  end if
467 
468  idimn=mod(idim,ndir)+1 ! 'Next' direction
469  idime=mod(idim+1,ndir)+1 ! Electric field direction
470  ! Store velocities
471  vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,mom(idimn))
472  vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,mom(idimn))
473  vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
474  +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
475  /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
476 
477  vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,mom(idime))
478  vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,mom(idime))
479  vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
480  +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
481  /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
482  case default
483  call mpistop('choose average, uct_contact,or uct_hll for type_ct!')
484  end select
485 
486  end subroutine mf_get_ct_velocity
487 
488  !> Calculate fluxes within ixO^L.
489  subroutine mf_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
491  use mod_usr_methods
492  use mod_geometry
493 
494  integer, intent(in) :: ixI^L, ixO^L, idim
495  ! conservative w
496  double precision, intent(in) :: wC(ixI^S,nw)
497  ! primitive w
498  double precision, intent(in) :: w(ixI^S,nw)
499  double precision, intent(in) :: x(ixI^S,1:ndim)
500  double precision,intent(out) :: f(ixI^S,nwflux)
501 
502  integer :: idir, ix^D
503 
504  do idir=1,ndir
505  {do ix^db=ixomin^db,ixomax^db\}
506  ! compute flux of magnetic field
507  ! f_i[b_k]=v_i*b_k-v_k*b_i
508  f(ix^d,mag(idir))=w(ix^d,mom(idim))*w(ix^d,mag(idir))-w(ix^d,mag(idim))*w(ix^d,mom(idir))
509  {end do\}
510  end do
511  if(mf_glm) then
512  {do ix^db=ixomin^db,ixomax^db\}
513  f(ix^d,mag(idim))=w(ix^d,psi_)
514  !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645
515  f(ix^d,psi_)=cmax_global**2*w(ix^d,mag(idim))
516  {end do\}
517  end if
518 
519  end subroutine mf_get_flux
520 
521  !> Add global source terms to update frictional velocity on complete domain
522  subroutine mf_velocity_update(qt,psa)
525  double precision, intent(in) :: qt !< Current time
526  type(state), target :: psa(max_blocks) !< Compute based on this state
527 
528  integer :: iigrid,igrid
529  logical :: stagger_flag
530  logical :: firstmf=.true.
531 
532  if(firstmf) then
533  ! point bc mpi datatype to partial type for velocity field
541  firstmf=.false.
542  end if
543 
544  !$OMP PARALLEL DO PRIVATE(igrid)
545  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
546  block=>psa(igrid)
547  call frictional_velocity(psa(igrid)%w,psa(igrid)%x,ixg^ll,ixm^ll)
548  end do
549  !$OMP END PARALLEL DO
550 
551  ! only update mf velocity in ghost cells
552  stagger_flag=stagger_grid
559  stagger_grid=.false.
560  bcphys=.false.
561  call getbc(qt,0.d0,psa,mom(1),ndir,.true.)
562  bcphys=.true.
569  stagger_grid=stagger_flag
570 
571  end subroutine mf_velocity_update
572 
573  !> w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
574  subroutine mf_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
576 
577  integer, intent(in) :: ixI^L, ixO^L
578  double precision, intent(in) :: qdt, dtfactor
579  double precision, intent(in) :: wCT(ixI^S,1:nw),wCTprim(ixI^S,1:nw), x(ixI^S,1:ndim)
580  double precision, intent(inout) :: w(ixI^S,1:nw)
581  logical, intent(in) :: qsourcesplit
582  logical, intent(inout) :: active
583 
584  if (.not. qsourcesplit) then
585 
586  ! Sources for resistivity in eqs. for e, B1, B2 and B3
587  if (abs(mf_eta)>smalldouble)then
588  active = .true.
589  call add_source_res2(qdt,ixi^l,ixo^l,wct,w,x)
590  end if
591 
592  if (mf_eta_hyper>0.d0)then
593  active = .true.
594  call add_source_hyperres(qdt,ixi^l,ixo^l,wct,w,x)
595  end if
596 
597  end if
598 
599  {^nooned
600  if(source_split_divb .eqv. qsourcesplit) then
601  ! Sources related to div B
602  select case (type_divb)
603  case (divb_none)
604  ! Do nothing
605  case (divb_glm)
606  active = .true.
607  call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
608  case (divb_powel)
609  active = .true.
610  call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
611  case (divb_janhunen)
612  active = .true.
613  call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
614  case (divb_linde)
615  active = .true.
616  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
617  case (divb_lindejanhunen)
618  active = .true.
619  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
620  call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
621  case (divb_lindepowel)
622  active = .true.
623  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
624  call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
625  case (divb_lindeglm)
626  active = .true.
627  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
628  call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
629  case (divb_ct)
630  continue ! Do nothing
631  case (divb_multigrid)
632  continue ! Do nothing
633  case default
634  call mpistop('Unknown divB fix')
635  end select
636  end if
637  }
638 
639  end subroutine mf_add_source
640 
641  subroutine frictional_velocity(w,x,ixI^L,ixO^L)
643 
644  integer, intent(in) :: ixI^L, ixO^L
645  double precision, intent(in) :: x(ixI^S,1:ndim)
646  double precision, intent(inout) :: w(ixI^S,1:nw)
647 
648  double precision :: xmin(ndim),xmax(ndim)
649  double precision :: decay(ixO^S)
650  double precision :: current(ixI^S,7-2*ndir:3),tmp(ixO^S)
651  integer :: ix^D, idirmin,idir,jdir,kdir
652 
653  call get_current(w,ixi^l,ixo^l,idirmin,current)
654  ! extrapolate current for the outmost layer,
655  ! because extrapolation of B at boundaries introduces artificial current
656 {
657  if(block%is_physical_boundary(2*^d-1)) then
658  if(slab) then
659  ! exclude boundary 5 in 3D Cartesian
660  if(2*^d-1==5.and.ndim==3) then
661  else
662  current(ixomin^d^d%ixO^s,:)=current(ixomin^d+1^d%ixO^s,:)
663  end if
664  else
665  ! exclude boundary 1 spherical/cylindrical
666  if(2*^d-1>1) then
667  current(ixomin^d^d%ixO^s,:)=current(ixomin^d+1^d%ixO^s,:)
668  end if
669  end if
670  end if
671  if(block%is_physical_boundary(2*^d)) then
672  current(ixomax^d^d%ixO^s,:)=current(ixomax^d-1^d%ixO^s,:)
673  end if
674 \}
675  w(ixo^s,mom(:))=0.d0
676  ! calculate Lorentz force
677  do idir=1,ndir; do jdir=idirmin,3; do kdir=1,ndir
678  if(lvc(idir,jdir,kdir)/=0) then
679  if(lvc(idir,jdir,kdir)==1) then
680  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))+current(ixo^s,jdir)*w(ixo^s,mag(kdir))
681  else
682  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-current(ixo^s,jdir)*w(ixo^s,mag(kdir))
683  end if
684  end if
685  end do; end do; end do
686 
687  tmp(ixo^s)=sum(w(ixo^s,mag(:))**2,dim=ndim+1)
688  ! frictional coefficient
689  where(tmp(ixo^s)/=0.d0)
690  tmp(ixo^s)=1.d0/(tmp(ixo^s)*mf_nu)
691  endwhere
692 
693  do idir=1,ndir
694  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))*tmp(ixo^s)
695  end do
696 
697  ! decay frictional velocity near selected boundaries
698  ^d&xmin(^d)=xprobmin^d\
699  ^d&xmax(^d)=xprobmax^d\
700  decay(ixo^s)=1.d0
701  do idir=1,ndim
702  if(mf_decay_scale(2*idir-1)>0.d0) decay(ixo^s)=decay(ixo^s)*&
703  (1.d0-exp(-(x(ixo^s,idir)-xmin(idir))/mf_decay_scale(2*idir-1)))
704  if(mf_decay_scale(2*idir)>0.d0) decay(ixo^s)=decay(ixo^s)*&
705  (1.d0-exp((x(ixo^s,idir)-xmax(idir))/mf_decay_scale(2*idir)))
706  end do
707 
708  ! saturate mf velocity at mf_vmax
709  tmp(ixo^s)=sqrt(sum(w(ixo^s,mom(:))**2,dim=ndim+1))/mf_vmax+1.d-12
710  tmp(ixo^s)=dtanh(tmp(ixo^s))/tmp(ixo^s)
711  do idir=1,ndir
712  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))*tmp(ixo^s)*decay(ixo^s)
713  end do
714 
715  end subroutine frictional_velocity
716 
717  !> Add resistive source to w within ixO Uses 3 point stencil (1 neighbour) in
718  !> each direction, non-conservative. If the fourthorder precompiler flag is
719  !> set, uses fourth order central difference for the laplacian. Then the
720  !> stencil is 5 (2 neighbours).
721  subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
723  use mod_usr_methods
724  use mod_geometry
725 
726  integer, intent(in) :: ixI^L, ixO^L
727  double precision, intent(in) :: qdt
728  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
729  double precision, intent(inout) :: w(ixI^S,1:nw)
730 
731  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
732  double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
733  double precision :: gradeta(ixI^S,1:ndim), Bf(ixI^S,1:ndir)
734  double precision :: tmp(ixI^S),tmp2(ixI^S)
735  integer :: ixA^L,idir,jdir,kdir,idirmin,idim,jxO^L,hxO^L,ix
736  integer :: lxO^L, kxO^L
737 
738  ! Calculating resistive sources involve one extra layer
739  if (mf_4th_order) then
740  ixa^l=ixo^l^ladd2;
741  else
742  ixa^l=ixo^l^ladd1;
743  end if
744 
745  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
746  call mpistop("Error in add_source_res1: Non-conforming input limits")
747 
748  ! Calculate current density and idirmin
749  call get_current(wct,ixi^l,ixo^l,idirmin,current)
750 
751  if (mf_eta>zero)then
752  eta(ixa^s)=mf_eta
753  gradeta(ixo^s,1:ndim)=zero
754  else
755  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
756  ! assumes that eta is not function of current?
757  do idim=1,ndim
758  call gradient(eta,ixi^l,ixo^l,idim,tmp)
759  gradeta(ixo^s,idim)=tmp(ixo^s)
760  end do
761  end if
762 
763  bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
764 
765  do idir=1,ndir
766  ! Put B_idir into tmp2 and eta*Laplace B_idir into tmp
767  if (mf_4th_order) then
768  tmp(ixo^s)=zero
769  tmp2(ixi^s)=bf(ixi^s,idir)
770  do idim=1,ndim
771  lxo^l=ixo^l+2*kr(idim,^d);
772  jxo^l=ixo^l+kr(idim,^d);
773  hxo^l=ixo^l-kr(idim,^d);
774  kxo^l=ixo^l-2*kr(idim,^d);
775  tmp(ixo^s)=tmp(ixo^s)+&
776  (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
777  /(12.0d0 * dxlevel(idim)**2)
778  end do
779  else
780  tmp(ixo^s)=zero
781  tmp2(ixi^s)=bf(ixi^s,idir)
782  do idim=1,ndim
783  jxo^l=ixo^l+kr(idim,^d);
784  hxo^l=ixo^l-kr(idim,^d);
785  tmp(ixo^s)=tmp(ixo^s)+&
786  (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/dxlevel(idim)**2
787  end do
788  end if
789 
790  ! Multiply by eta
791  tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
792 
793  ! Subtract grad(eta) x J = eps_ijk d_j eta J_k if eta is non-constant
794  if (mf_eta<zero)then
795  do jdir=1,ndim; do kdir=idirmin,3
796  if (lvc(idir,jdir,kdir)/=0)then
797  if (lvc(idir,jdir,kdir)==1)then
798  tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
799  else
800  tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
801  end if
802  end if
803  end do; end do
804  end if
805 
806  ! Add sources related to eta*laplB-grad(eta) x J to B and e
807  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
808 
809  end do ! idir
810 
811  end subroutine add_source_res1
812 
813  !> Add resistive source to w within ixO
814  !> Uses 5 point stencil (2 neighbours) in each direction, conservative
815  subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
817  use mod_usr_methods
818  use mod_geometry
819 
820  integer, intent(in) :: ixI^L, ixO^L
821  double precision, intent(in) :: qdt
822  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
823  double precision, intent(inout) :: w(ixI^S,1:nw)
824 
825  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
826  double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S),curlj(ixI^S,1:3)
827  double precision :: tmpvec(ixI^S,1:3)
828  integer :: ixA^L,idir,idirmin,idirmin1
829 
830  ! resistive term already added as an electric field component
831  if(stagger_grid .and. ndim==ndir) return
832 
833  ixa^l=ixo^l^ladd2;
834 
835  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
836  call mpistop("Error in add_source_res2: Non-conforming input limits")
837 
838  ixa^l=ixo^l^ladd1;
839  ! Calculate current density within ixL: J=curl B, thus J_i=eps_ijk*d_j B_k
840  ! Determine exact value of idirmin while doing the loop.
841  call get_current(wct,ixi^l,ixa^l,idirmin,current)
842 
843  if (mf_eta>zero)then
844  eta(ixa^s)=mf_eta
845  else
846  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
847  end if
848 
849  ! dB/dt= -curl(J*eta), thus B_i=B_i-eps_ijk d_j Jeta_k
850  tmpvec(ixa^s,1:ndir)=zero
851  do idir=idirmin,3
852  tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
853  end do
854  curlj=0.d0
855  call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
856  if(stagger_grid) then
857  if(ndim==2.and.ndir==3) then
858  ! if 2.5D
859  w(ixo^s,mag(ndir)) = w(ixo^s,mag(ndir))-qdt*curlj(ixo^s,ndir)
860  end if
861  else
862  w(ixo^s,mag(1:ndir)) = w(ixo^s,mag(1:ndir))-qdt*curlj(ixo^s,1:ndir)
863  end if
864 
865  end subroutine add_source_res2
866 
867  !> Add Hyper-resistive source to w within ixO
868  !> Uses 9 point stencil (4 neighbours) in each direction.
869  subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
871  use mod_geometry
872 
873  integer, intent(in) :: ixI^L, ixO^L
874  double precision, intent(in) :: qdt
875  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
876  double precision, intent(inout) :: w(ixI^S,1:nw)
877  !.. local ..
878  double precision :: current(ixI^S,7-2*ndir:3)
879  double precision :: tmpvec(ixI^S,1:3),tmpvec2(ixI^S,1:3),tmp(ixI^S),ehyper(ixI^S,1:3)
880  integer :: ixA^L,idir,jdir,kdir,idirmin,idirmin1
881 
882  ixa^l=ixo^l^ladd3;
883  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
884  call mpistop("Error in add_source_hyperres: Non-conforming input limits")
885 
886  call get_current(wct,ixi^l,ixa^l,idirmin,current)
887  tmpvec(ixa^s,1:ndir)=zero
888  do jdir=idirmin,3
889  tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
890  end do
891 
892  ixa^l=ixo^l^ladd2;
893  call curlvector(tmpvec,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
894 
895  ixa^l=ixo^l^ladd1;
896  tmpvec(ixa^s,1:ndir)=zero
897  call curlvector(tmpvec2,ixi^l,ixa^l,tmpvec,idirmin1,1,3)
898  ehyper(ixa^s,1:ndir) = - tmpvec(ixa^s,1:ndir)*mf_eta_hyper
899 
900  ixa^l=ixo^l;
901  tmpvec2(ixa^s,1:ndir)=zero
902  call curlvector(ehyper,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
903 
904  do idir=1,ndir
905  w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
906  end do
907 
908  end subroutine add_source_hyperres
909 
910  subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
911  ! Add divB related sources to w within ixO
912  ! corresponding to Dedner JCP 2002, 175, 645 _equation 24_
913  ! giving the EGLM-MHD scheme
915  use mod_geometry
916 
917  integer, intent(in) :: ixI^L, ixO^L
918  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
919  double precision, intent(inout) :: w(ixI^S,1:nw)
920  double precision:: divb(ixI^S)
921  double precision :: gradPsi(ixI^S)
922  integer :: idim,idir
923 
924  ! We calculate now div B
925  call get_divb(wct,ixi^l,ixo^l,divb, mf_divb_4thorder)
926 
927  ! dPsi/dt = - Ch^2/Cp^2 Psi
928  if (mf_glm_alpha < zero) then
929  w(ixo^s,psi_) = abs(mf_glm_alpha)*wct(ixo^s,psi_)
930  else
931  ! implicit update of Psi variable
932  ! equation (27) in Mignone 2010 J. Com. Phys. 229, 2117
933  if(slab_uniform) then
934  w(ixo^s,psi_) = dexp(-qdt*cmax_global*mf_glm_alpha/minval(dxlevel(:)))*w(ixo^s,psi_)
935  else
936  w(ixo^s,psi_) = dexp(-qdt*cmax_global*mf_glm_alpha/minval(block%ds(ixo^s,:),dim=ndim+1))*w(ixo^s,psi_)
937  end if
938  end if
939 
940  ! gradient of Psi
941  do idim=1,ndim
942  select case(typegrad)
943  case("central")
944  call gradient(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
945  case("limited")
946  call gradients(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
947  end select
948  end do
949 
950  end subroutine add_source_glm
951 
952  !> Add divB related sources to w within ixO corresponding to Powel
953  subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
955 
956  integer, intent(in) :: ixI^L, ixO^L
957  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
958  double precision, intent(inout) :: w(ixI^S,1:nw)
959 
960  double precision :: divb(ixI^S)
961  integer :: idir
962 
963  ! We calculate now div B
964  call get_divb(wct,ixi^l,ixo^l,divb, mf_divb_4thorder)
965 
966  ! b = b - qdt v * div b
967  do idir=1,ndir
968  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*wct(ixo^s,mom(idir))*divb(ixo^s)
969  end do
970 
971  end subroutine add_source_powel
972 
973  subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
974  ! Add divB related sources to w within ixO
975  ! corresponding to Janhunen, just the term in the induction equation.
977 
978  integer, intent(in) :: ixI^L, ixO^L
979  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
980  double precision, intent(inout) :: w(ixI^S,1:nw)
981  double precision :: divb(ixI^S)
982  integer :: idir
983 
984  ! We calculate now div B
985  call get_divb(wct,ixi^l,ixo^l,divb, mf_divb_4thorder)
986 
987  ! b = b - qdt v * div b
988  do idir=1,ndir
989  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*wct(ixo^s,mom(idir))*divb(ixo^s)
990  end do
991 
992  end subroutine add_source_janhunen
993 
994  subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
995  ! Add Linde's divB related sources to wnew within ixO
997  use mod_geometry
998 
999  integer, intent(in) :: ixI^L, ixO^L
1000  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
1001  double precision, intent(inout) :: w(ixI^S,1:nw)
1002  double precision :: divb(ixI^S),graddivb(ixI^S)
1003  integer :: idim, idir, ixp^L, i^D, iside
1004  logical, dimension(-1:1^D&) :: leveljump
1005 
1006  ! Calculate div B
1007  ixp^l=ixo^l^ladd1;
1008  call get_divb(wct,ixi^l,ixp^l,divb, mf_divb_4thorder)
1009 
1010  ! for AMR stability, retreat one cell layer from the boarders of level jump
1011  {do i^db=-1,1\}
1012  if(i^d==0|.and.) cycle
1013  if(neighbor_type(i^d,block%igrid)==2 .or. neighbor_type(i^d,block%igrid)==4) then
1014  leveljump(i^d)=.true.
1015  else
1016  leveljump(i^d)=.false.
1017  end if
1018  {end do\}
1019 
1020  ixp^l=ixo^l;
1021  do idim=1,ndim
1022  select case(idim)
1023  {case(^d)
1024  do iside=1,2
1025  i^dd=kr(^dd,^d)*(2*iside-3);
1026  if (leveljump(i^dd)) then
1027  if (iside==1) then
1028  ixpmin^d=ixomin^d-i^d
1029  else
1030  ixpmax^d=ixomax^d-i^d
1031  end if
1032  end if
1033  end do
1034  \}
1035  end select
1036  end do
1037 
1038  ! Add Linde's diffusive terms
1039  do idim=1,ndim
1040  ! Calculate grad_idim(divb)
1041  select case(typegrad)
1042  case("central")
1043  call gradient(divb,ixi^l,ixp^l,idim,graddivb)
1044  case("limited")
1045  call gradients(divb,ixi^l,ixp^l,idim,graddivb)
1046  end select
1047 
1048  ! Multiply by Linde's eta*dt = divbdiff*(c_max*dx)*dt = divbdiff*dx**2
1049  if (slab_uniform) then
1050  graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
1051  else
1052  graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
1053  /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
1054  end if
1055 
1056  w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
1057  end do
1058 
1059  end subroutine add_source_linde
1060 
1061 
1062  !> get dimensionless div B = |divB| * volume / area / |B|
1063  subroutine get_normalized_divb(w,ixI^L,ixO^L,divb)
1064 
1066 
1067  integer, intent(in) :: ixi^l, ixo^l
1068  double precision, intent(in) :: w(ixi^s,1:nw)
1069  double precision :: divb(ixi^s), dsurface(ixi^s)
1070 
1071  double precision :: invb(ixo^s)
1072  integer :: ixa^l,idims
1073 
1074  call get_divb(w,ixi^l,ixo^l,divb)
1075  invb(ixo^s)=sqrt(sum(w(ixo^s, mag(:))**2, dim=ndim+1))
1076  where(invb(ixo^s)/=0.d0)
1077  invb(ixo^s)=1.d0/invb(ixo^s)
1078  end where
1079  if(slab_uniform) then
1080  divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/dxlevel(:))
1081  else
1082  ixamin^d=ixomin^d-1;
1083  ixamax^d=ixomax^d-1;
1084  dsurface(ixo^s)= sum(block%surfaceC(ixo^s,:),dim=ndim+1)
1085  do idims=1,ndim
1086  ixa^l=ixo^l-kr(idims,^d);
1087  dsurface(ixo^s)=dsurface(ixo^s)+block%surfaceC(ixa^s,idims)
1088  end do
1089  divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
1090  block%dvolume(ixo^s)/dsurface(ixo^s)
1091  end if
1092 
1093  end subroutine get_normalized_divb
1094 
1095  !> Calculate idirmin and the idirmin:3 components of the common current array
1096  !> make sure that dxlevel(^D) is set correctly.
1097  subroutine get_current(w,ixI^L,ixO^L,idirmin,current,fourthorder)
1099  use mod_geometry
1100 
1101  integer, intent(in) :: ixo^l, ixi^l
1102  double precision, intent(in) :: w(ixi^s,1:nw)
1103  integer, intent(out) :: idirmin
1104  logical, intent(in), optional :: fourthorder
1105 
1106  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
1107  double precision :: current(ixi^s,7-2*ndir:3)
1108  integer :: idir, idirmin0
1109 
1110  idirmin0 = 7-2*ndir
1111 
1112  if(present(fourthorder)) then
1113  call curlvector(w(ixi^s,mag(1:ndir)),ixi^l,ixo^l,current,idirmin,idirmin0,ndir,fourthorder)
1114  else
1115  call curlvector(w(ixi^s,mag(1:ndir)),ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
1116  end if
1117 
1118  end subroutine get_current
1119 
1120  !> If resistivity is not zero, check diffusion time limit for dt
1121  subroutine mf_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
1123  use mod_usr_methods
1124 
1125  integer, intent(in) :: ixI^L, ixO^L
1126  double precision, intent(inout) :: dtnew
1127  double precision, intent(in) :: dx^D
1128  double precision, intent(in) :: w(ixI^S,1:nw)
1129  double precision, intent(in) :: x(ixI^S,1:ndim)
1130 
1131  double precision :: dxarr(ndim)
1132  double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
1133  integer :: idirmin,idim
1134 
1135  dtnew = bigdouble
1136 
1137  ^d&dxarr(^d)=dx^d;
1138  if (mf_eta>zero)then
1139  dtnew=dtdiffpar*minval(dxarr(1:ndim))**2/mf_eta
1140  else if (mf_eta<zero)then
1141  call get_current(w,ixi^l,ixo^l,idirmin,current)
1142  call usr_special_resistivity(w,ixi^l,ixo^l,idirmin,x,current,eta)
1143  dtnew=bigdouble
1144  do idim=1,ndim
1145  if(slab_uniform) then
1146  dtnew=min(dtnew,&
1147  dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
1148  else
1149  dtnew=min(dtnew,&
1150  dtdiffpar/(smalldouble+maxval(eta(ixo^s)/block%ds(ixo^s,idim)**2)))
1151  end if
1152  end do
1153  end if
1154 
1155  if(mf_eta_hyper>zero) then
1156  if(slab_uniform) then
1157  dtnew=min(dtdiffpar*minval(dxarr(1:ndim))**4/mf_eta_hyper,dtnew)
1158  else
1159  dtnew=min(dtdiffpar*minval(block%ds(ixo^s,1:ndim))**4/mf_eta_hyper,dtnew)
1160  end if
1161  end if
1162  end subroutine mf_get_dt
1163 
1164  ! Add geometrical source terms to w
1165  subroutine mf_add_source_geom(qdt,dtfactor, ixI^L,ixO^L,wCT,wprim,w,x)
1167  use mod_geometry
1168 
1169  integer, intent(in) :: ixI^L, ixO^L
1170  double precision, intent(in) :: qdt, dtfactor, x(ixI^S,1:ndim)
1171  double precision, intent(inout) :: wCT(ixI^S,1:nw), wprim(ixI^S,1:nw), w(ixI^S,1:nw)
1172 
1173  double precision :: tmp,invr,cs2
1174  integer :: iw,idir,ix^D
1175 
1176  integer :: mr_,mphi_ ! Polar var. names
1177  integer :: br_,bphi_
1178 
1179  mr_=mom(1); mphi_=mom(1)-1+phi_ ! Polar var. names
1180  br_=mag(1); bphi_=mag(1)-1+phi_
1181 
1182  select case (coordinate)
1183  case (cylindrical)
1184  if(phi_>0) then
1185  if(.not.stagger_grid) then
1186  {do ix^db=ixomin^db,ixomax^db\}
1187  w(ix^d,bphi_)=w(ix^d,bphi_)+qdt/x(ix^d,1)*&
1188  (wct(ix^d,bphi_)*wct(ix^d,mr_) &
1189  -wct(ix^d,br_)*wct(ix^d,mphi_))
1190  {end do\}
1191  end if
1192  end if
1193  if(mf_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,psi_)/x(ixo^s,1)
1194  case (spherical)
1195  {do ix^db=ixomin^db,ixomax^db\}
1196  invr=qdt/x(ix^d,1)
1197  ! b1
1198  if(mf_glm) then
1199  w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wct(ix^d,psi_)
1200  end if
1201 
1202  {^nooned
1203  ! b2
1204  if(.not.stagger_grid) then
1205  tmp=wct(ix^d,mom(1))*wct(ix^d,mag(2))-wct(ix^d,mom(2))*wct(ix^d,mag(1))
1206  cs2=1.d0/tan(x(ix^d,2))
1207  if(mf_glm) then
1208  tmp=tmp+cs2*wct(ix^d,psi_)
1209  end if
1210  w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
1211  end if
1212  }
1213 
1214  if(ndir==3) then
1215  {^ifoned
1216  ! b3
1217  if(.not.stagger_grid) then
1218  w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
1219  (wct(ix^d,mom(1))*wct(ix^d,mag(3)) &
1220  -wct(ix^d,mom(3))*wct(ix^d,mag(1)))
1221  end if
1222  }
1223  {^nooned
1224  ! b3
1225  if(.not.stagger_grid) then
1226  w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
1227  (wct(ix^d,mom(1))*wct(ix^d,mag(3)) &
1228  -wct(ix^d,mom(3))*wct(ix^d,mag(1)) &
1229  -(wct(ix^d,mom(3))*wct(ix^d,mag(2)) &
1230  -wct(ix^d,mom(2))*wct(ix^d,mag(3)))*cs2)
1231  end if
1232  }
1233  end if
1234  {end do\}
1235  end select
1236 
1237  end subroutine mf_add_source_geom
1238 
1239  subroutine mf_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
1241  use mod_usr_methods
1242  integer, intent(in) :: ixI^L, ixO^L, idir
1243  double precision, intent(in) :: qt
1244  double precision, intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
1245  double precision, intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
1246  type(state) :: s
1247  double precision :: dB(ixI^S), dPsi(ixI^S)
1248 
1249  ! Solve the Riemann problem for the linear 2x2 system for normal
1250  ! B-field and GLM_Psi according to Dedner 2002:
1251  ! This implements eq. (42) in Dedner et al. 2002 JcP 175
1252  ! Gives the Riemann solution on the interface
1253  ! for the normal B component and Psi in the GLM-MHD system.
1254  ! 23/04/2013 Oliver Porth
1255  db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
1256  dpsi(ixo^s) = wrp(ixo^s,psi_) - wlp(ixo^s,psi_)
1257 
1258  wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
1259  - 0.5d0/cmax_global * dpsi(ixo^s)
1260  wlp(ixo^s,psi_) = 0.5d0 * (wrp(ixo^s,psi_) + wlp(ixo^s,psi_)) &
1261  - 0.5d0*cmax_global * db(ixo^s)
1262 
1263  wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
1264  wrp(ixo^s,psi_) = wlp(ixo^s,psi_)
1265 
1266  end subroutine mf_modify_wlr
1267 
1268  subroutine mf_boundary_adjust(igrid,psb)
1270  integer, intent(in) :: igrid
1271  type(state), target :: psb(max_blocks)
1272 
1273  integer :: iB, idims, iside, ixO^L, i^D
1274 
1275  block=>ps(igrid)
1276  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
1277  do idims=1,ndim
1278  ! to avoid using as yet unknown corner info in more than 1D, we
1279  ! fill only interior mesh ranges of the ghost cell ranges at first,
1280  ! and progressively enlarge the ranges to include corners later
1281  do iside=1,2
1282  i^d=kr(^d,idims)*(2*iside-3);
1283  if (neighbor_type(i^d,igrid)/=1) cycle
1284  ib=(idims-1)*2+iside
1285  if(.not.boundary_divbfix(ib)) cycle
1286  if(any(typeboundary(:,ib)==bc_special)) then
1287  ! MF nonlinear force-free B field extrapolation and data driven
1288  ! require normal B of the first ghost cell layer to be untouched by
1289  ! fixdivB=0 process, set boundary_divbfix_skip(iB)=1 in par file
1290  select case (idims)
1291  {case (^d)
1292  if (iside==2) then
1293  ! maximal boundary
1294  ixomin^dd=ixghi^d+1-nghostcells+boundary_divbfix_skip(2*^d)^d%ixOmin^dd=ixglo^dd;
1295  ixomax^dd=ixghi^dd;
1296  else
1297  ! minimal boundary
1298  ixomin^dd=ixglo^dd;
1299  ixomax^dd=ixglo^d-1+nghostcells-boundary_divbfix_skip(2*^d-1)^d%ixOmax^dd=ixghi^dd;
1300  end if \}
1301  end select
1302  call fixdivb_boundary(ixg^ll,ixo^l,psb(igrid)%w,psb(igrid)%x,ib)
1303  end if
1304  end do
1305  end do
1306 
1307  end subroutine mf_boundary_adjust
1308 
1309  subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
1311 
1312  integer, intent(in) :: ixG^L,ixO^L,iB
1313  double precision, intent(inout) :: w(ixG^S,1:nw)
1314  double precision, intent(in) :: x(ixG^S,1:ndim)
1315 
1316  double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
1317  integer :: ix^D,ixF^L
1318 
1319  select case(ib)
1320  case(1)
1321  ! 2nd order CD for divB=0 to set normal B component better
1322  {^iftwod
1323  ixfmin1=ixomin1+1
1324  ixfmax1=ixomax1+1
1325  ixfmin2=ixomin2+1
1326  ixfmax2=ixomax2-1
1327  if(slab_uniform) then
1328  dx1x2=dxlevel(1)/dxlevel(2)
1329  do ix1=ixfmax1,ixfmin1,-1
1330  w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
1331  +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
1332  w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
1333  enddo
1334  else
1335  do ix1=ixfmax1,ixfmin1,-1
1336  w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
1337  w(ix1,ixfmin2:ixfmax2,mag(1)))*block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
1338  +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
1339  block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
1340  -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
1341  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
1342  /block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
1343  end do
1344  end if
1345  }
1346  {^ifthreed
1347  ixfmin1=ixomin1+1
1348  ixfmax1=ixomax1+1
1349  ixfmin2=ixomin2+1
1350  ixfmax2=ixomax2-1
1351  ixfmin3=ixomin3+1
1352  ixfmax3=ixomax3-1
1353  if(slab_uniform) then
1354  dx1x2=dxlevel(1)/dxlevel(2)
1355  dx1x3=dxlevel(1)/dxlevel(3)
1356  do ix1=ixfmax1,ixfmin1,-1
1357  w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1358  w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
1359  +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
1360  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
1361  +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
1362  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
1363  end do
1364  else
1365  do ix1=ixfmax1,ixfmin1,-1
1366  w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1367  ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
1368  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
1369  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
1370  +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
1371  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
1372  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
1373  -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
1374  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
1375  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
1376  +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
1377  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
1378  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
1379  -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
1380  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1381  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
1382  /block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
1383  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
1384  end do
1385  end if
1386  }
1387  case(2)
1388  {^iftwod
1389  ixfmin1=ixomin1-1
1390  ixfmax1=ixomax1-1
1391  ixfmin2=ixomin2+1
1392  ixfmax2=ixomax2-1
1393  if(slab_uniform) then
1394  dx1x2=dxlevel(1)/dxlevel(2)
1395  do ix1=ixfmin1,ixfmax1
1396  w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
1397  -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
1398  w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
1399  enddo
1400  else
1401  do ix1=ixfmin1,ixfmax1
1402  w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
1403  w(ix1,ixfmin2:ixfmax2,mag(1)))*block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
1404  -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
1405  block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
1406  +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
1407  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
1408  /block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
1409  end do
1410  end if
1411  }
1412  {^ifthreed
1413  ixfmin1=ixomin1-1
1414  ixfmax1=ixomax1-1
1415  ixfmin2=ixomin2+1
1416  ixfmax2=ixomax2-1
1417  ixfmin3=ixomin3+1
1418  ixfmax3=ixomax3-1
1419  if(slab_uniform) then
1420  dx1x2=dxlevel(1)/dxlevel(2)
1421  dx1x3=dxlevel(1)/dxlevel(3)
1422  do ix1=ixfmin1,ixfmax1
1423  w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1424  w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
1425  -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
1426  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
1427  -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
1428  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
1429  end do
1430  else
1431  do ix1=ixfmin1,ixfmax1
1432  w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1433  ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
1434  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
1435  block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
1436  -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
1437  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
1438  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
1439  +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
1440  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
1441  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
1442  -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
1443  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
1444  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
1445  +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
1446  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1447  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
1448  /block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
1449  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
1450  end do
1451  end if
1452  }
1453  case(3)
1454  {^iftwod
1455  ixfmin1=ixomin1+1
1456  ixfmax1=ixomax1-1
1457  ixfmin2=ixomin2+1
1458  ixfmax2=ixomax2+1
1459  if(slab_uniform) then
1460  dx2x1=dxlevel(2)/dxlevel(1)
1461  do ix2=ixfmax2,ixfmin2,-1
1462  w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
1463  +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
1464  w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
1465  enddo
1466  else
1467  do ix2=ixfmax2,ixfmin2,-1
1468  w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
1469  w(ixfmin1:ixfmax1,ix2,mag(2)))*block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
1470  +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
1471  block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
1472  -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
1473  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
1474  /block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
1475  end do
1476  end if
1477  }
1478  {^ifthreed
1479  ixfmin1=ixomin1+1
1480  ixfmax1=ixomax1-1
1481  ixfmin3=ixomin3+1
1482  ixfmax3=ixomax3-1
1483  ixfmin2=ixomin2+1
1484  ixfmax2=ixomax2+1
1485  if(slab_uniform) then
1486  dx2x1=dxlevel(2)/dxlevel(1)
1487  dx2x3=dxlevel(2)/dxlevel(3)
1488  do ix2=ixfmax2,ixfmin2,-1
1489  w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
1490  ix2+1,ixfmin3:ixfmax3,mag(2)) &
1491  +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
1492  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
1493  +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
1494  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
1495  end do
1496  else
1497  do ix2=ixfmax2,ixfmin2,-1
1498  w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
1499  ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
1500  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
1501  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
1502  +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
1503  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1504  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
1505  -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
1506  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1507  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
1508  +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
1509  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
1510  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
1511  -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
1512  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1513  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
1514  /block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
1515  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
1516  end do
1517  end if
1518  }
1519  case(4)
1520  {^iftwod
1521  ixfmin1=ixomin1+1
1522  ixfmax1=ixomax1-1
1523  ixfmin2=ixomin2-1
1524  ixfmax2=ixomax2-1
1525  if(slab_uniform) then
1526  dx2x1=dxlevel(2)/dxlevel(1)
1527  do ix2=ixfmin2,ixfmax2
1528  w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
1529  -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
1530  w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
1531  end do
1532  else
1533  do ix2=ixfmin2,ixfmax2
1534  w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
1535  w(ixfmin1:ixfmax1,ix2,mag(2)))*block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
1536  -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
1537  block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
1538  +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
1539  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
1540  /block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
1541  end do
1542  end if
1543  }
1544  {^ifthreed
1545  ixfmin1=ixomin1+1
1546  ixfmax1=ixomax1-1
1547  ixfmin3=ixomin3+1
1548  ixfmax3=ixomax3-1
1549  ixfmin2=ixomin2-1
1550  ixfmax2=ixomax2-1
1551  if(slab_uniform) then
1552  dx2x1=dxlevel(2)/dxlevel(1)
1553  dx2x3=dxlevel(2)/dxlevel(3)
1554  do ix2=ixfmin2,ixfmax2
1555  w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
1556  ix2-1,ixfmin3:ixfmax3,mag(2)) &
1557  -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
1558  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
1559  -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
1560  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
1561  end do
1562  else
1563  do ix2=ixfmin2,ixfmax2
1564  w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
1565  ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
1566  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
1567  block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
1568  -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
1569  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1570  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
1571  +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
1572  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1573  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
1574  -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
1575  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
1576  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
1577  +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
1578  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1579  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
1580  /block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
1581  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
1582  end do
1583  end if
1584  }
1585  {^ifthreed
1586  case(5)
1587  ixfmin1=ixomin1+1
1588  ixfmax1=ixomax1-1
1589  ixfmin2=ixomin2+1
1590  ixfmax2=ixomax2-1
1591  ixfmin3=ixomin3+1
1592  ixfmax3=ixomax3+1
1593  if(slab_uniform) then
1594  dx3x1=dxlevel(3)/dxlevel(1)
1595  dx3x2=dxlevel(3)/dxlevel(2)
1596  do ix3=ixfmax3,ixfmin3,-1
1597  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
1598  ixfmin2:ixfmax2,ix3+1,mag(3)) &
1599  +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
1600  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
1601  +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
1602  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
1603  end do
1604  else
1605  do ix3=ixfmax3,ixfmin3,-1
1606  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
1607  ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
1608  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
1609  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
1610  +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
1611  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1612  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
1613  -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
1614  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1615  block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
1616  +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
1617  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
1618  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
1619  -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
1620  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
1621  block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
1622  /block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
1623  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
1624  end do
1625  end if
1626  case(6)
1627  ixfmin1=ixomin1+1
1628  ixfmax1=ixomax1-1
1629  ixfmin2=ixomin2+1
1630  ixfmax2=ixomax2-1
1631  ixfmin3=ixomin3-1
1632  ixfmax3=ixomax3-1
1633  if(slab_uniform) then
1634  dx3x1=dxlevel(3)/dxlevel(1)
1635  dx3x2=dxlevel(3)/dxlevel(2)
1636  do ix3=ixfmin3,ixfmax3
1637  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
1638  ixfmin2:ixfmax2,ix3-1,mag(3)) &
1639  -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
1640  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
1641  -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
1642  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
1643  end do
1644  else
1645  do ix3=ixfmin3,ixfmax3
1646  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
1647  ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
1648  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
1649  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
1650  -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
1651  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1652  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
1653  +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
1654  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1655  block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
1656  -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
1657  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
1658  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
1659  +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
1660  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
1661  block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
1662  /block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
1663  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
1664  end do
1665  end if
1666  }
1667  case default
1668  call mpistop("Special boundary is not defined for this region")
1669  end select
1670 
1671  end subroutine fixdivb_boundary
1672 
1673  {^nooned
1674  subroutine mf_clean_divb_multigrid(qdt, qt, active)
1675  use mod_forest
1678  use mod_geometry
1679 
1680  double precision, intent(in) :: qdt !< Current time step
1681  double precision, intent(in) :: qt !< Current time
1682  logical, intent(inout) :: active !< Output if the source is active
1683 
1684  double precision :: tmp(ixg^t), grad(ixg^t, ndim)
1685  double precision :: res
1686  double precision, parameter :: max_residual = 1d-3
1687  double precision, parameter :: residual_reduction = 1d-10
1688  integer :: id
1689  integer, parameter :: max_its = 50
1690  double precision :: residual_it(max_its), max_divb
1691  integer :: iigrid, igrid
1692  integer :: n, nc, lvl, ix^l, ixc^l, idim
1693  type(tree_node), pointer :: pnode
1694 
1695  mg%operator_type = mg_laplacian
1696 
1697  ! Set boundary conditions
1698  do n = 1, 2*ndim
1699  idim = (n+1)/2
1700  select case (typeboundary(mag(idim), n))
1701  case (bc_symm)
1702  ! d/dx B = 0, take phi = 0
1703  mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
1704  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1705  case (bc_asymm)
1706  ! B = 0, so grad(phi) = 0
1707  mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
1708  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1709  case (bc_cont)
1710  mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
1711  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1712  case (bc_special)
1713  ! Assume Dirichlet boundary conditions, derivative zero
1714  mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
1715  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1716  case (bc_periodic)
1717  ! Nothing to do here
1718  case default
1719  write(*,*) "mf_clean_divb_multigrid warning: unknown boundary type"
1720  mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
1721  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1722  end select
1723  end do
1724 
1725  ix^l=ixm^ll^ladd1;
1726  max_divb = 0.0d0
1727 
1728  ! Store divergence of B as right-hand side
1729  do iigrid = 1, igridstail
1730  igrid = igrids(iigrid);
1731  pnode => igrid_to_node(igrid, mype)%node
1732  id = pnode%id
1733  lvl = mg%boxes(id)%lvl
1734  nc = mg%box_size_lvl(lvl)
1735 
1736  ! Geometry subroutines expect this to be set
1737  block => ps(igrid)
1738  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
1739 
1740  call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^ll, ixm^ll, tmp, &
1742  mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(ixm^t)
1743  max_divb = max(max_divb, maxval(abs(tmp(ixm^t))))
1744  end do
1745 
1746  ! Solve laplacian(phi) = divB
1747  if(stagger_grid) then
1748  call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
1749  mpi_max, icomm, ierrmpi)
1750 
1751  if (mype == 0) print *, "Performing multigrid divB cleaning"
1752  if (mype == 0) print *, "iteration vs residual"
1753  ! Solve laplacian(phi) = divB
1754  do n = 1, max_its
1755  call mg_fas_fmg(mg, n>1, max_res=residual_it(n))
1756  if (mype == 0) write(*, "(I4,E11.3)") n, residual_it(n)
1757  if (residual_it(n) < residual_reduction * max_divb) exit
1758  end do
1759  if (mype == 0 .and. n > max_its) then
1760  print *, "divb_multigrid warning: not fully converged"
1761  print *, "current amplitude of divb: ", residual_it(max_its)
1762  print *, "multigrid smallest grid: ", &
1763  mg%domain_size_lvl(:, mg%lowest_lvl)
1764  print *, "note: smallest grid ideally has <= 8 cells"
1765  print *, "multigrid dx/dy/dz ratio: ", mg%dr(:, 1)/mg%dr(1, 1)
1766  print *, "note: dx/dy/dz should be similar"
1767  end if
1768  else
1769  do n = 1, max_its
1770  call mg_fas_vcycle(mg, max_res=res)
1771  if (res < max_residual) exit
1772  end do
1773  if (res > max_residual) call mpistop("divb_multigrid: no convergence")
1774  end if
1775 
1776 
1777  ! Correct the magnetic field
1778  do iigrid = 1, igridstail
1779  igrid = igrids(iigrid);
1780  pnode => igrid_to_node(igrid, mype)%node
1781  id = pnode%id
1782 
1783  ! Geometry subroutines expect this to be set
1784  block => ps(igrid)
1785  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
1786 
1787  ! Compute the gradient of phi
1788  tmp(ix^s) = mg%boxes(id)%cc({:,}, mg_iphi)
1789 
1790  if(stagger_grid) then
1791  do idim =1, ndim
1792  ixcmin^d=ixmlo^d-kr(idim,^d);
1793  ixcmax^d=ixmhi^d;
1794  call gradientx(tmp,ps(igrid)%x,ixg^ll,ixc^l,idim,grad(ixg^t,idim),.false.)
1795  ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
1796  end do
1797  call mf_face_to_center(ixm^ll,ps(igrid))
1798  else
1799  do idim = 1, ndim
1800  call gradient(tmp,ixg^ll,ixm^ll,idim,grad(ixg^t, idim))
1801  end do
1802  ! Apply the correction B* = B - gradient(phi)
1803  ps(igrid)%w(ixm^t, mag(1:ndim)) = &
1804  ps(igrid)%w(ixm^t, mag(1:ndim)) - grad(ixm^t, :)
1805  end if
1806 
1807  end do
1808 
1809  active = .true.
1810 
1811  end subroutine mf_clean_divb_multigrid
1812  }
1813 
1814  subroutine mf_update_faces(ixI^L,ixO^L,qt,qdt,wprim,fC,fE,sCT,s,vcts)
1816 
1817  integer, intent(in) :: ixI^L, ixO^L
1818  double precision, intent(in) :: qt, qdt
1819  ! cell-center primitive variables
1820  double precision, intent(in) :: wprim(ixI^S,1:nw)
1821  type(state) :: sCT, s
1822  type(ct_velocity) :: vcts
1823  double precision, intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
1824  double precision, intent(inout) :: fE(ixI^S,sdim:3)
1825 
1826  select case(type_ct)
1827  case('average')
1828  call update_faces_average(ixi^l,ixo^l,qt,qdt,fc,fe,sct,s)
1829  case('uct_contact')
1830  call update_faces_contact(ixi^l,ixo^l,qt,qdt,wprim,fc,fe,sct,s,vcts)
1831  case('uct_hll')
1832  call update_faces_hll(ixi^l,ixo^l,qt,qdt,fe,sct,s,vcts)
1833  case default
1834  call mpistop('choose average, uct_contact,or uct_hll for type_ct!')
1835  end select
1836 
1837  end subroutine mf_update_faces
1838 
1839  !> get electric field though averaging neighors to update faces in CT
1840  subroutine update_faces_average(ixI^L,ixO^L,qt,qdt,fC,fE,sCT,s)
1842  use mod_usr_methods
1843 
1844  integer, intent(in) :: ixI^L, ixO^L
1845  double precision, intent(in) :: qt, qdt
1846  type(state) :: sCT, s
1847  double precision, intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
1848  double precision, intent(inout) :: fE(ixI^S,sdim:3)
1849 
1850  double precision :: circ(ixI^S,1:ndim)
1851  double precision :: E_resi(ixI^S,sdim:3)
1852  integer :: ix^D,ixC^L,ixA^L,i1kr^D,i2kr^D
1853  integer :: idim1,idim2,idir,iwdim1,iwdim2
1854 
1855  associate(bfaces=>s%ws,x=>s%x)
1856 
1857  ! Calculate contribution to FEM of each edge,
1858  ! that is, estimate value of line integral of
1859  ! electric field in the positive idir direction.
1860 
1861  ! if there is resistivity, get eta J
1862  if(mf_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,e_resi)
1863 
1864  do idim1=1,ndim
1865  iwdim1 = mag(idim1)
1866  i1kr^d=kr(idim1,^d);
1867  do idim2=1,ndim
1868  iwdim2 = mag(idim2)
1869  i2kr^d=kr(idim2,^d);
1870  do idir=sdim,3! Direction of line integral
1871  ! Allow only even permutations
1872  if (lvc(idim1,idim2,idir)==1) then
1873  ixcmax^d=ixomax^d;
1874  ixcmin^d=ixomin^d+kr(idir,^d)-1;
1875  ! average cell-face electric field to cell edges
1876  {do ix^db=ixcmin^db,ixcmax^db\}
1877  fe(ix^d,idir)=quarter*&
1878  (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
1879  -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
1880  ! add resistive electric field at cell edges E=-vxB+eta J
1881  if(mf_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
1882  if(record_electric_field) s%we(ix^d,idir)=fe(ix^d,idir)
1883  ! times time step and edge length
1884  fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
1885  {end do\}
1886  end if
1887  end do
1888  end do
1889  end do
1890 
1891  ! allow user to change inductive electric field, especially for boundary driven applications
1892  if(associated(usr_set_electric_field)) &
1893  call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
1894 
1895  circ(ixi^s,1:ndim)=zero
1896 
1897  ! Calculate circulation on each face
1898 
1899  do idim1=1,ndim ! Coordinate perpendicular to face
1900  ixcmax^d=ixomax^d;
1901  ixcmin^d=ixomin^d-kr(idim1,^d);
1902  do idim2=1,ndim
1903  ixa^l=ixc^l-kr(idim2,^d);
1904  do idir=sdim,3 ! Direction of line integral
1905  ! Assemble indices
1906  if(lvc(idim1,idim2,idir)==1) then
1907  ! Add line integrals in direction idir
1908  circ(ixc^s,idim1)=circ(ixc^s,idim1)&
1909  +(fe(ixc^s,idir)&
1910  -fe(ixa^s,idir))
1911  else if(lvc(idim1,idim2,idir)==-1) then
1912  ! Add line integrals in direction idir
1913  circ(ixc^s,idim1)=circ(ixc^s,idim1)&
1914  -(fe(ixc^s,idir)&
1915  -fe(ixa^s,idir))
1916  end if
1917  end do
1918  end do
1919  ! Divide by the area of the face to get dB/dt
1920  circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
1921  ! Time update
1922  bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
1923  end do
1924 
1925  end associate
1926 
1927  end subroutine update_faces_average
1928 
1929  !> update faces using UCT contact mode by Gardiner and Stone 2005 JCP 205, 509
1930  subroutine update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
1932  use mod_usr_methods
1933 
1934  integer, intent(in) :: ixI^L, ixO^L
1935  double precision, intent(in) :: qt, qdt
1936  ! cell-center primitive variables
1937  double precision, intent(in) :: wp(ixI^S,1:nw)
1938  type(state) :: sCT, s
1939  type(ct_velocity) :: vcts
1940  double precision, intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
1941  double precision, intent(inout) :: fE(ixI^S,sdim:3)
1942 
1943  double precision :: circ(ixI^S,1:ndim)
1944  ! electric field at cell centers
1945  double precision :: ECC(ixI^S,sdim:3)
1946  ! gradient of E at left and right side of a cell face
1947  double precision :: EL(ixI^S),ER(ixI^S)
1948  ! gradient of E at left and right side of a cell corner
1949  double precision :: ELC(ixI^S),ERC(ixI^S)
1950  ! current on cell edges
1951  double precision :: jce(ixI^S,sdim:3)
1952  integer :: hxC^L,ixC^L,jxC^L,ixA^L,ixB^L
1953  integer :: idim1,idim2,idir,iwdim1,iwdim2
1954 
1955  associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm)
1956 
1957  ecc=0.d0
1958  ! Calculate electric field at cell centers
1959  do idim1=1,ndim; do idim2=1,ndim; do idir=sdim,3
1960  if(lvc(idim1,idim2,idir)==1)then
1961  ecc(ixi^s,idir)=ecc(ixi^s,idir)+wp(ixi^s,mag(idim1))*wp(ixi^s,mom(idim2))
1962  else if(lvc(idim1,idim2,idir)==-1) then
1963  ecc(ixi^s,idir)=ecc(ixi^s,idir)-wp(ixi^s,mag(idim1))*wp(ixi^s,mom(idim2))
1964  endif
1965  enddo; enddo; enddo
1966 
1967  ! if there is resistivity, get eta J
1968  if(mf_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,jce)
1969 
1970  ! Calculate contribution to FEM of each edge,
1971  ! that is, estimate value of line integral of
1972  ! electric field in the positive idir direction.
1973  ! evaluate electric field along cell edges according to equation (41)
1974  do idim1=1,ndim
1975  iwdim1 = mag(idim1)
1976  do idim2=1,ndim
1977  iwdim2 = mag(idim2)
1978  do idir=sdim,3 ! Direction of line integral
1979  ! Allow only even permutations
1980  if (lvc(idim1,idim2,idir)==1) then
1981  ixcmax^d=ixomax^d;
1982  ixcmin^d=ixomin^d+kr(idir,^d)-1;
1983  ! Assemble indices
1984  jxc^l=ixc^l+kr(idim1,^d);
1985  hxc^l=ixc^l+kr(idim2,^d);
1986  ! average cell-face electric field to cell edges
1987  fe(ixc^s,idir)=quarter*&
1988  (fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
1989  -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
1990 
1991  ! add slope in idim2 direction from equation (50)
1992  ixamin^d=ixcmin^d;
1993  ixamax^d=ixcmax^d+kr(idim1,^d);
1994  el(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(ixa^s,idir)
1995  hxc^l=ixa^l+kr(idim2,^d);
1996  er(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(hxc^s,idir)
1997  where(vnorm(ixc^s,idim1)>0.d0)
1998  elc(ixc^s)=el(ixc^s)
1999  else where(vnorm(ixc^s,idim1)<0.d0)
2000  elc(ixc^s)=el(jxc^s)
2001  else where
2002  elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
2003  end where
2004  hxc^l=ixc^l+kr(idim2,^d);
2005  where(vnorm(hxc^s,idim1)>0.d0)
2006  erc(ixc^s)=er(ixc^s)
2007  else where(vnorm(hxc^s,idim1)<0.d0)
2008  erc(ixc^s)=er(jxc^s)
2009  else where
2010  erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
2011  end where
2012  fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
2013 
2014  ! add slope in idim1 direction from equation (50)
2015  jxc^l=ixc^l+kr(idim2,^d);
2016  ixamin^d=ixcmin^d;
2017  ixamax^d=ixcmax^d+kr(idim2,^d);
2018  el(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(ixa^s,idir)
2019  hxc^l=ixa^l+kr(idim1,^d);
2020  er(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(hxc^s,idir)
2021  where(vnorm(ixc^s,idim2)>0.d0)
2022  elc(ixc^s)=el(ixc^s)
2023  else where(vnorm(ixc^s,idim2)<0.d0)
2024  elc(ixc^s)=el(jxc^s)
2025  else where
2026  elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
2027  end where
2028  hxc^l=ixc^l+kr(idim1,^d);
2029  where(vnorm(hxc^s,idim2)>0.d0)
2030  erc(ixc^s)=er(ixc^s)
2031  else where(vnorm(hxc^s,idim2)<0.d0)
2032  erc(ixc^s)=er(jxc^s)
2033  else where
2034  erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
2035  end where
2036  fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
2037 
2038  ! add current component of electric field at cell edges E=-vxB+eta J
2039  if(mf_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+jce(ixc^s,idir)
2040 
2041  if(record_electric_field) s%we(ixc^s,idir)=fe(ixc^s,idir)
2042  ! times time step and edge length
2043  fe(ixc^s,idir)=fe(ixc^s,idir)*qdt*s%dsC(ixc^s,idir)
2044  if (.not.slab) then
2045  where(abs(x(ixc^s,r_)+half*dxlevel(r_))<1.0d-9)
2046  fe(ixc^s,idir)=zero
2047  end where
2048  end if
2049  end if
2050  end do
2051  end do
2052  end do
2053 
2054  ! allow user to change inductive electric field, especially for boundary driven applications
2055  if(associated(usr_set_electric_field)) &
2056  call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
2057 
2058  circ(ixi^s,1:ndim)=zero
2059 
2060  ! Calculate circulation on each face
2061  do idim1=1,ndim ! Coordinate perpendicular to face
2062  ixcmax^d=ixomax^d;
2063  ixcmin^d=ixomin^d-kr(idim1,^d);
2064  do idim2=1,ndim
2065  do idir=sdim,3 ! Direction of line integral
2066  ! Assemble indices
2067  if(lvc(idim1,idim2,idir)/=0) then
2068  hxc^l=ixc^l-kr(idim2,^d);
2069  ! Add line integrals in direction idir
2070  circ(ixc^s,idim1)=circ(ixc^s,idim1)&
2071  +lvc(idim1,idim2,idir)&
2072  *(fe(ixc^s,idir)&
2073  -fe(hxc^s,idir))
2074  end if
2075  end do
2076  end do
2077  ! Divide by the area of the face to get dB/dt
2078  where(s%surfaceC(ixc^s,idim1) > 1.0d-9*s%dvolume(ixc^s))
2079  circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
2080  elsewhere
2081  circ(ixc^s,idim1)=zero
2082  end where
2083  ! Time update cell-face magnetic field component
2084  bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
2085  end do
2086 
2087  end associate
2088 
2089  end subroutine update_faces_contact
2090 
2091  !> update faces
2092  subroutine update_faces_hll(ixI^L,ixO^L,qt,qdt,fE,sCT,s,vcts)
2095  use mod_usr_methods
2096 
2097  integer, intent(in) :: ixI^L, ixO^L
2098  double precision, intent(in) :: qt, qdt
2099  double precision, intent(inout) :: fE(ixI^S,sdim:3)
2100  type(state) :: sCT, s
2101  type(ct_velocity) :: vcts
2102 
2103  double precision :: vtilL(ixI^S,2)
2104  double precision :: vtilR(ixI^S,2)
2105  double precision :: btilL(ixI^S,ndim)
2106  double precision :: btilR(ixI^S,ndim)
2107  double precision :: cp(ixI^S,2)
2108  double precision :: cm(ixI^S,2)
2109  double precision :: circ(ixI^S,1:ndim)
2110  ! current on cell edges
2111  double precision :: jce(ixI^S,sdim:3)
2112  integer :: hxC^L,ixC^L,ixCp^L,jxC^L,ixCm^L
2113  integer :: idim1,idim2,idir
2114 
2115  associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
2116  cbarmax=>vcts%cbarmax)
2117 
2118  ! Calculate contribution to FEM of each edge,
2119  ! that is, estimate value of line integral of
2120  ! electric field in the positive idir direction.
2121 
2122  ! Loop over components of electric field
2123 
2124  ! idir: electric field component we need to calculate
2125  ! idim1: directions in which we already performed the reconstruction
2126  ! idim2: directions in which we perform the reconstruction
2127 
2128  ! if there is resistivity, get eta J
2129  if(mf_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,jce)
2130 
2131  do idir=sdim,3
2132  ! Indices
2133  ! idir: electric field component
2134  ! idim1: one surface
2135  ! idim2: the other surface
2136  ! cyclic permutation: idim1,idim2,idir=1,2,3
2137  ! Velocity components on the surface
2138  ! follow cyclic premutations:
2139  ! Sx(1),Sx(2)=y,z ; Sy(1),Sy(2)=z,x ; Sz(1),Sz(2)=x,y
2140 
2141  ixcmax^d=ixomax^d;
2142  ixcmin^d=ixomin^d-1+kr(idir,^d);
2143 
2144  ! Set indices and directions
2145  idim1=mod(idir,3)+1
2146  idim2=mod(idir+1,3)+1
2147 
2148  jxc^l=ixc^l+kr(idim1,^d);
2149  ixcp^l=ixc^l+kr(idim2,^d);
2150 
2151  ! Reconstruct transverse transport velocities
2152  call reconstruct(ixi^l,ixc^l,idim2,vbarc(ixi^s,idim1,1),&
2153  vtill(ixi^s,2),vtilr(ixi^s,2))
2154 
2155  call reconstruct(ixi^l,ixc^l,idim1,vbarc(ixi^s,idim2,2),&
2156  vtill(ixi^s,1),vtilr(ixi^s,1))
2157 
2158  ! Reconstruct magnetic fields
2159  ! Eventhough the arrays are larger, reconstruct works with
2160  ! the limits ixG.
2161  call reconstruct(ixi^l,ixc^l,idim2,bfacesct(ixi^s,idim1),&
2162  btill(ixi^s,idim1),btilr(ixi^s,idim1))
2163 
2164  call reconstruct(ixi^l,ixc^l,idim1,bfacesct(ixi^s,idim2),&
2165  btill(ixi^s,idim2),btilr(ixi^s,idim2))
2166 
2167  ! Take the maximum characteristic
2168 
2169  cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
2170  cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
2171 
2172  cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
2173  cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
2174 
2175 
2176  ! Calculate eletric field
2177  fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
2178  + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
2179  - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
2180  /(cp(ixc^s,1)+cm(ixc^s,1)) &
2181  +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
2182  + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
2183  - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
2184  /(cp(ixc^s,2)+cm(ixc^s,2))
2185 
2186  ! add current component of electric field at cell edges E=-vxB+eta J
2187  if(mf_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+jce(ixc^s,idir)
2188 
2189  if(record_electric_field) s%we(ixc^s,idir)=fe(ixc^s,idir)
2190  ! times time step and edge length
2191  fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
2192 
2193  if (.not.slab) then
2194  where(abs(x(ixc^s,r_)+half*dxlevel(r_)).lt.1.0d-9)
2195  fe(ixc^s,idir)=zero
2196  end where
2197  end if
2198 
2199  end do
2200 
2201  ! allow user to change inductive electric field, especially for boundary driven applications
2202  if(associated(usr_set_electric_field)) &
2203  call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
2204 
2205  circ(ixi^s,1:ndim)=zero
2206 
2207  ! Calculate circulation on each face: interal(fE dot dl)
2208 
2209  do idim1=1,ndim ! Coordinate perpendicular to face
2210  ixcmax^d=ixomax^d;
2211  ixcmin^d=ixomin^d-kr(idim1,^d);
2212  do idim2=1,ndim
2213  do idir=sdim,3 ! Direction of line integral
2214  ! Assemble indices
2215  if(lvc(idim1,idim2,idir)/=0) then
2216  hxc^l=ixc^l-kr(idim2,^d);
2217  ! Add line integrals in direction idir
2218  circ(ixc^s,idim1)=circ(ixc^s,idim1)&
2219  +lvc(idim1,idim2,idir)&
2220  *(fe(ixc^s,idir)&
2221  -fe(hxc^s,idir))
2222  end if
2223  end do
2224  end do
2225  ! Divide by the area of the face to get dB/dt
2226  where(s%surfaceC(ixc^s,idim1) > 1.0d-9*s%dvolume(ixc^s))
2227  circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
2228  elsewhere
2229  circ(ixc^s,idim1)=zero
2230  end where
2231  ! Time update
2232  bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
2233  end do
2234 
2235  end associate
2236  end subroutine update_faces_hll
2237 
2238  !> calculate eta J at cell edges
2239  subroutine get_resistive_electric_field(ixI^L,ixO^L,sCT,s,jce)
2241  use mod_usr_methods
2242  use mod_geometry
2243 
2244  integer, intent(in) :: ixI^L, ixO^L
2245  type(state), intent(in) :: sCT, s
2246  ! current on cell edges
2247  double precision :: jce(ixI^S,sdim:3)
2248 
2249  ! current on cell centers
2250  double precision :: jcc(ixI^S,7-2*ndir:3)
2251  ! location at cell faces
2252  double precision :: xs(ixGs^T,1:ndim)
2253  ! resistivity
2254  double precision :: eta(ixI^S)
2255  double precision :: gradi(ixGs^T)
2256  integer :: ix^D,ixC^L,ixA^L,ixB^L,idir,idirmin,idim1,idim2
2257 
2258  associate(x=>s%x,dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
2259  ! calculate current density at cell edges
2260  jce=0.d0
2261  do idim1=1,ndim
2262  do idim2=1,ndim
2263  do idir=sdim,3
2264  if (lvc(idim1,idim2,idir)==0) cycle
2265  ixcmax^d=ixomax^d+1;
2266  ixcmin^d=ixomin^d+kr(idir,^d)-2;
2267  ixbmax^d=ixcmax^d-kr(idir,^d)+1;
2268  ixbmin^d=ixcmin^d;
2269  ! current at transverse faces
2270  xs(ixb^s,:)=x(ixb^s,:)
2271  xs(ixb^s,idim2)=x(ixb^s,idim2)+half*dx(ixb^s,idim2)
2272  call gradientx(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi,.false.)
2273  if (lvc(idim1,idim2,idir)==1) then
2274  jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
2275  else
2276  jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
2277  end if
2278  end do
2279  end do
2280  end do
2281  ! get resistivity
2282  if(mf_eta>zero)then
2283  jce(ixi^s,:)=jce(ixi^s,:)*mf_eta
2284  else
2285  ixa^l=ixo^l^ladd1;
2286  call get_current(wct,ixi^l,ixa^l,idirmin,jcc)
2287  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,jcc,eta)
2288  ! calcuate eta on cell edges
2289  do idir=sdim,3
2290  ixcmax^d=ixomax^d;
2291  ixcmin^d=ixomin^d+kr(idir,^d)-1;
2292  jcc(ixc^s,idir)=0.d0
2293  {do ix^db=0,1\}
2294  if({ ix^d==1 .and. ^d==idir | .or.}) cycle
2295  ixamin^d=ixcmin^d+ix^d;
2296  ixamax^d=ixcmax^d+ix^d;
2297  jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
2298  {end do\}
2299  jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
2300  jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
2301  enddo
2302  end if
2303 
2304  end associate
2305  end subroutine get_resistive_electric_field
2306 
2307  !> calculate cell-center values from face-center values
2308  subroutine mf_face_to_center(ixO^L,s)
2310  ! Non-staggered interpolation range
2311  integer, intent(in) :: ixo^l
2312  type(state) :: s
2313 
2314  integer :: ix^d
2315 
2316  ! calculate cell-center values from face-center values in 2nd order
2317  {!dir$ ivdep
2318  do ix^db=ixomin^db,ixomax^db\}
2319  {^ifthreed
2320  s%w(ix^d,mag(1))=half/s%surface(ix^d,1)*(s%ws(ix^d,1)*s%surfaceC(ix^d,1)&
2321  +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
2322  s%w(ix^d,mag(2))=half/s%surface(ix^d,2)*(s%ws(ix^d,2)*s%surfaceC(ix^d,2)&
2323  +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
2324  s%w(ix^d,mag(3))=half/s%surface(ix^d,3)*(s%ws(ix^d,3)*s%surfaceC(ix^d,3)&
2325  +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
2326  }
2327  {^iftwod
2328  s%w(ix^d,mag(1))=half/s%surface(ix^d,1)*(s%ws(ix^d,1)*s%surfaceC(ix^d,1)&
2329  +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
2330  s%w(ix^d,mag(2))=half/s%surface(ix^d,2)*(s%ws(ix^d,2)*s%surfaceC(ix^d,2)&
2331  +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
2332  }
2333  {end do\}
2334 
2335  ! calculate cell-center values from face-center values in 4th order
2336  !do idim=1,ndim
2337  ! gxO^L=ixO^L-2*kr(idim,^D);
2338  ! hxO^L=ixO^L-kr(idim,^D);
2339  ! jxO^L=ixO^L+kr(idim,^D);
2340 
2341  ! ! Interpolate to cell barycentre using fourth order central formula
2342  ! w(ixO^S,mag(idim))=(0.0625d0/s%surface(ixO^S,idim))*&
2343  ! ( -ws(gxO^S,idim)*s%surfaceC(gxO^S,idim) &
2344  ! +9.0d0*ws(hxO^S,idim)*s%surfaceC(hxO^S,idim) &
2345  ! +9.0d0*ws(ixO^S,idim)*s%surfaceC(ixO^S,idim) &
2346  ! -ws(jxO^S,idim)*s%surfaceC(jxO^S,idim) )
2347  !end do
2348 
2349  ! calculate cell-center values from face-center values in 6th order
2350  !do idim=1,ndim
2351  ! fxO^L=ixO^L-3*kr(idim,^D);
2352  ! gxO^L=ixO^L-2*kr(idim,^D);
2353  ! hxO^L=ixO^L-kr(idim,^D);
2354  ! jxO^L=ixO^L+kr(idim,^D);
2355  ! kxO^L=ixO^L+2*kr(idim,^D);
2356 
2357  ! ! Interpolate to cell barycentre using sixth order central formula
2358  ! w(ixO^S,mag(idim))=(0.00390625d0/s%surface(ixO^S,idim))* &
2359  ! ( +3.0d0*ws(fxO^S,idim)*s%surfaceC(fxO^S,idim) &
2360  ! -25.0d0*ws(gxO^S,idim)*s%surfaceC(gxO^S,idim) &
2361  ! +150.0d0*ws(hxO^S,idim)*s%surfaceC(hxO^S,idim) &
2362  ! +150.0d0*ws(ixO^S,idim)*s%surfaceC(ixO^S,idim) &
2363  ! -25.0d0*ws(jxO^S,idim)*s%surfaceC(jxO^S,idim) &
2364  ! +3.0d0*ws(kxO^S,idim)*s%surfaceC(kxO^S,idim) )
2365  !end do
2366 
2367  end subroutine mf_face_to_center
2368 
2369  !> calculate magnetic field from vector potential
2370  subroutine b_from_vector_potential(ixIs^L, ixI^L, ixO^L, ws, x)
2373 
2374  integer, intent(in) :: ixis^l, ixi^l, ixo^l
2375  double precision, intent(inout) :: ws(ixis^s,1:nws)
2376  double precision, intent(in) :: x(ixi^s,1:ndim)
2377 
2378  double precision :: adummy(ixis^s,1:3)
2379 
2380  call b_from_vector_potentiala(ixis^l, ixi^l, ixo^l, ws, x, adummy)
2381 
2382  end subroutine b_from_vector_potential
2383 
2386 
2387  double precision :: sum_jbb_ipe, sum_j_ipe, sum_l_ipe, f_i_ipe, volume_pe
2388  double precision :: sum_jbb, sum_j, sum_l, f_i, volume, cw_sin_theta
2389  integer :: iigrid, igrid, ix^d
2390  integer :: amode, istatus(mpi_status_size)
2391  integer, save :: fhmf
2392  logical :: patchwi(ixg^t)
2393  logical, save :: logmfopened=.false.
2394  character(len=800) :: filename,filehead
2395  character(len=800) :: line,datastr
2396 
2397  sum_jbb_ipe = 0.d0
2398  sum_j_ipe = 0.d0
2399  sum_l_ipe = 0.d0
2400  f_i_ipe = 0.d0
2401  volume_pe=0.d0
2402  do iigrid=1,igridstail; igrid=igrids(iigrid);
2403  block=>ps(igrid)
2404  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
2405  call mask_inner(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x,patchwi,volume_pe)
2406  sum_jbb_ipe = sum_jbb_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
2407  ps(igrid)%x,1,patchwi)
2408  sum_j_ipe = sum_j_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
2409  ps(igrid)%x,2,patchwi)
2410  f_i_ipe=f_i_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
2411  ps(igrid)%x,3,patchwi)
2412  sum_l_ipe = sum_l_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
2413  ps(igrid)%x,4,patchwi)
2414  end do
2415  call mpi_reduce(sum_jbb_ipe,sum_jbb,1,mpi_double_precision,&
2416  mpi_sum,0,icomm,ierrmpi)
2417  call mpi_reduce(sum_j_ipe,sum_j,1,mpi_double_precision,mpi_sum,0,&
2418  icomm,ierrmpi)
2419  call mpi_reduce(f_i_ipe,f_i,1,mpi_double_precision,mpi_sum,0,&
2420  icomm,ierrmpi)
2421  call mpi_reduce(sum_l_ipe,sum_l,1,mpi_double_precision,mpi_sum,0,&
2422  icomm,ierrmpi)
2423  call mpi_reduce(volume_pe,volume,1,mpi_double_precision,mpi_sum,0,&
2424  icomm,ierrmpi)
2425 
2426  if(mype==0) then
2427  ! current- and volume-weighted average of the sine of the angle
2428  ! between the magnetic field B and the current density J
2429  cw_sin_theta = sum_jbb/sum_j
2430  ! volume-weighted average of the absolute value of the fractional
2431  ! magnetic flux change
2432  f_i = f_i/volume
2433  sum_j=sum_j/volume
2434  sum_l=sum_l/volume
2435  if(.not.logmfopened) then
2436  ! generate filename
2437  write(filename,"(a,a)") trim(base_filename), "_mf_metrics.csv"
2438 
2439  ! Delete the log when not doing a restart run
2440  if (restart_from_file == undefined) then
2441  open(unit=20,file=trim(filename),status='replace')
2442  close(20, status='delete')
2443  end if
2444 
2445  amode=ior(mpi_mode_create,mpi_mode_wronly)
2446  amode=ior(amode,mpi_mode_append)
2447  call mpi_file_open(mpi_comm_self,filename,amode,mpi_info_null,fhmf,ierrmpi)
2448  logmfopened=.true.
2449  filehead=" it,time,CW_sin_theta,current,Lorenz_force,f_i"
2450  if (it == 0) then
2451  call mpi_file_write(fhmf,filehead,len_trim(filehead), &
2452  mpi_character,istatus,ierrmpi)
2453  call mpi_file_write(fhmf,achar(10),1,mpi_character,istatus,ierrmpi)
2454  end if
2455  end if
2456  line=''
2457  write(datastr,'(i6,a)') it,','
2458  line=trim(line)//trim(datastr)
2459  write(datastr,'(es13.6,a)') global_time,','
2460  line=trim(line)//trim(datastr)
2461  write(datastr,'(es13.6,a)') cw_sin_theta,','
2462  line=trim(line)//trim(datastr)
2463  write(datastr,'(es13.6,a)') sum_j,','
2464  line=trim(line)//trim(datastr)
2465  write(datastr,'(es13.6,a)') sum_l,','
2466  line=trim(line)//trim(datastr)
2467  write(datastr,'(es13.6)') f_i
2468  line=trim(line)//trim(datastr)//new_line('A')
2469  call mpi_file_write(fhmf,line,len_trim(line),mpi_character,istatus,ierrmpi)
2470  if(.not.time_advance) then
2471  call mpi_file_close(fhmf,ierrmpi)
2472  logmfopened=.false.
2473  end if
2474  end if
2475 
2476  end subroutine record_force_free_metrics
2477 
2478  subroutine mask_inner(ixI^L,ixO^L,w,x,patchwi,volume_pe)
2480 
2481  integer, intent(in) :: ixI^L,ixO^L
2482  double precision, intent(in):: w(ixI^S,nw),x(ixI^S,1:ndim)
2483  double precision, intent(inout) :: volume_pe
2484  logical, intent(out) :: patchwi(ixI^S)
2485 
2486  double precision :: xO^L
2487  integer :: ix^D
2488 
2489  {xomin^d = xprobmin^d + 0.05d0*(xprobmax^d-xprobmin^d)\}
2490  {xomax^d = xprobmax^d - 0.05d0*(xprobmax^d-xprobmin^d)\}
2491  if(slab) then
2492  xomin^nd = xprobmin^nd
2493  else
2494  xomin1 = xprobmin1
2495  end if
2496 
2497  {do ix^db=ixomin^db,ixomax^db\}
2498  if({ x(ix^dd,^d) > xomin^d .and. x(ix^dd,^d) < xomax^d | .and. }) then
2499  patchwi(ix^d)=.true.
2500  volume_pe=volume_pe+block%dvolume(ix^d)
2501  else
2502  patchwi(ix^d)=.false.
2503  endif
2504  {end do\}
2505 
2506  end subroutine mask_inner
2507 
2508  function integral_grid_mf(ixI^L,ixO^L,w,x,iw,patchwi)
2510 
2511  integer, intent(in) :: ixi^l,ixo^l,iw
2512  double precision, intent(in) :: x(ixi^s,1:ndim)
2513  double precision :: w(ixi^s,nw+nwauxio)
2514  logical, intent(in) :: patchwi(ixi^s)
2515 
2516  double precision, dimension(ixI^S,7-2*ndir:3) :: current
2517  double precision, dimension(ixI^S,1:ndir) :: bvec,qvec
2518  double precision :: integral_grid_mf,tmp(ixi^s),bm2
2519  integer :: ix^d,idirmin0,idirmin,idir,jdir,kdir
2520 
2521  integral_grid_mf=0.d0
2522  select case(iw)
2523  case(1)
2524  ! Sum(|JxB|/|B|*dvolume)
2525  bvec(ixo^s,:)=w(ixo^s,mag(:))
2526  call get_current(w,ixi^l,ixo^l,idirmin,current)
2527  ! calculate Lorentz force
2528  qvec(ixo^s,1:ndir)=zero
2529  do idir=1,ndir; do jdir=idirmin,3; do kdir=1,ndir
2530  if(lvc(idir,jdir,kdir)/=0)then
2531  if(lvc(idir,jdir,kdir)==1)then
2532  qvec(ixo^s,idir)=qvec(ixo^s,idir)+current(ixo^s,jdir)*bvec(ixo^s,kdir)
2533  else
2534  qvec(ixo^s,idir)=qvec(ixo^s,idir)-current(ixo^s,jdir)*bvec(ixo^s,kdir)
2535  end if
2536  end if
2537  end do; end do; end do
2538 
2539  {do ix^db=ixomin^db,ixomax^db\}
2540  if(patchwi(ix^d)) then
2541  bm2=sum(bvec(ix^d,:)**2)
2542  if(bm2/=0.d0) bm2=1.d0/bm2
2543  integral_grid_mf=integral_grid_mf+sqrt(sum(qvec(ix^d,:)**2)*&
2544  bm2)*block%dvolume(ix^d)
2545  end if
2546  {end do\}
2547  case(2)
2548  ! Sum(|J|*dvolume)
2549  call get_current(w,ixi^l,ixo^l,idirmin,current)
2550  {do ix^db=ixomin^db,ixomax^db\}
2551  if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+sqrt(sum(current(ix^d,:)**2))*&
2552  block%dvolume(ix^d)
2553  {end do\}
2554  case(3)
2555  ! f_i solenoidal property of B: (dvolume |div B|)/(dsurface |B|)
2556  ! Sum(f_i*dvolume)
2557  call get_normalized_divb(w,ixi^l,ixo^l,tmp)
2558  {do ix^db=ixomin^db,ixomax^db\}
2559  if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+tmp(ix^d)*block%dvolume(ix^d)
2560  {end do\}
2561  case(4)
2562  ! Sum(|JxB|*dvolume)
2563  bvec(ixo^s,:)=w(ixo^s,mag(:))
2564  call get_current(w,ixi^l,ixo^l,idirmin,current)
2565  ! calculate Lorentz force
2566  qvec(ixo^s,1:ndir)=zero
2567  do idir=1,ndir; do jdir=idirmin,3; do kdir=1,ndir
2568  if(lvc(idir,jdir,kdir)/=0)then
2569  if(lvc(idir,jdir,kdir)==1)then
2570  qvec(ixo^s,idir)=qvec(ixo^s,idir)+current(ixo^s,jdir)*bvec(ixo^s,kdir)
2571  else
2572  qvec(ixo^s,idir)=qvec(ixo^s,idir)-current(ixo^s,jdir)*bvec(ixo^s,kdir)
2573  end if
2574  end if
2575  end do; end do; end do
2576 
2577  {do ix^db=ixomin^db,ixomax^db\}
2578  if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+&
2579  sqrt(sum(qvec(ix^d,:)**2))*block%dvolume(ix^d)
2580  {end do\}
2581  end select
2582  return
2583  end function integral_grid_mf
2584 
2585 end module mod_mf_phys
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
subroutine b_from_vector_potentiala(ixIsL, ixIL, ixOL, ws, x, A)
calculate magnetic field from vector potential A at cell edges
subroutine reconstruct(ixIL, ixCL, idir, q, qL, qR)
Reconstruct scalar q within ixO^L to 1/2 dx in direction idir Return both left and right reconstructe...
Module with basic grid data structures.
Definition: mod_forest.t:2
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
Definition: mod_forest.t:32
integer, dimension(:), allocatable, public mag
Indices of the magnetic field.
subroutine, public get_divb(w, ixIL, ixOL, divb, fourthorder)
Calculate div B within ixO.
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
integer coordinate
Definition: mod_geometry.t:7
integer, parameter cylindrical
Definition: mod_geometry.t:10
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
Definition: mod_geometry.t:321
subroutine curlvector(qvec, ixIL, ixOL, curlvec, idirmin, idirmin0, ndir0, fourthorder)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
Definition: mod_geometry.t:664
subroutine gradients(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir first use limiter to go from cell cente...
Definition: mod_geometry.t:458
subroutine gradientx(q, x, ixIL, ixOL, idir, gradq, fourth_order)
Calculate gradient of a scalar q in direction idir at cell interfaces.
Definition: mod_geometry.t:363
update ghost cells of all blocks including physical boundaries
integer, dimension(0:3^d &), target type_recv_p_p1
integer, dimension( :^d &), pointer type_recv_r
integer, dimension(-1:1^d &), target type_send_srl_p1
integer, dimension(0:3^d &), target type_send_p_f
integer, dimension( :^d &), pointer type_send_p
integer, dimension( :^d &), pointer type_send_srl
integer, dimension(-1:1^d &), target type_send_r_p1
subroutine create_bc_mpi_datatype(nwstart, nwbc)
integer, dimension(0:3^d &), target type_recv_r_f
integer, dimension(-1:1^d &), target type_recv_srl_f
integer, dimension(-1:1^d &), target type_send_r_f
integer, dimension(-1:1^d &), target type_send_srl_f
integer, dimension(0:3^d &), target type_send_p_p1
integer, dimension( :^d &), pointer type_send_r
integer, dimension(0:3^d &), target type_recv_r_p1
subroutine getbc(time, qdt, psb, nwstart, nwbc, req_diag)
do update ghost cells of all blocks including physical boundaries
integer, dimension( :^d &), pointer type_recv_p
integer, dimension(-1:1^d &), target type_recv_srl_p1
integer, dimension(0:3^d &), target type_recv_p_f
integer, dimension( :^d &), pointer type_recv_srl
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.
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
character(len=std_len) typegrad
double precision unit_charge
Physical scaling factor for charge.
integer ixghi
Upper index of grid block arrays.
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
double precision unit_time
Physical scaling factor for time.
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
integer, parameter bc_asymm
double precision global_time
The global simulation time.
double precision unit_mass
Physical scaling factor for mass.
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer it
Number of time steps taken.
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
double precision unit_numberdensity
Physical scaling factor for number density.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision unit_length
Physical scaling factor for length.
logical stagger_grid
True for using stagger grid.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
logical use_particles
Use particles module or not.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer icomm
The MPI communicator.
double precision bdip
amplitude of background dipolar, quadrupolar, octupolar, user's field
integer mype
The rank of the current MPI task.
double precision, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
logical slab
Cartesian geometry or not.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision unit_magneticfield
Physical scaling factor for magnetic field.
integer nwauxio
Number of auxiliary variables that are only included in output.
double precision unit_velocity
Physical scaling factor for velocity.
logical time_advance
do time evolving
character(len=std_len) restart_from_file
If not 'unavailable', resume from snapshot with this base file name.
double precision c_norm
Normalised speed of light.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
integer, parameter bc_cont
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter bc_symm
character(len= *), parameter undefined
logical need_global_cmax
need global maximal wave speed
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
logical use_multigrid
Use multigrid (only available in 2D and 3D)
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
character(len=std_len) base_filename
Base file name for simulation output, which will be followed by a number.
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
logical record_electric_field
True for record electric field.
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Magnetofriction module.
Definition: mod_mf_phys.t:2
double precision, public mf_nu
viscosity coefficient in s cm^-2 for solar corona (Cheung 2012 ApJ)
Definition: mod_mf_phys.t:12
logical, public, protected mf_divb_4thorder
Whether divB is computed with a fourth order approximation.
Definition: mod_mf_phys.t:76
subroutine frictional_velocity(w, x, ixIL, ixOL)
Definition: mod_mf_phys.t:642
subroutine mf_get_ct_velocity(vcts, wLp, wRp, ixIL, ixOL, idim, cmax, cmin)
prepare velocities for ct methods
Definition: mod_mf_phys.t:437
subroutine add_source_hyperres(qdt, ixIL, ixOL, wCT, w, x)
Add Hyper-resistive source to w within ixO Uses 9 point stencil (4 neighbours) in each direction.
Definition: mod_mf_phys.t:870
subroutine, public mf_to_primitive(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
Definition: mod_mf_phys.t:390
subroutine, public mf_phys_init()
Definition: mod_mf_phys.t:158
subroutine mf_get_flux(wC, w, x, ixIL, ixOL, idim, f)
Calculate fluxes within ixO^L.
Definition: mod_mf_phys.t:490
subroutine add_source_glm(qdt, ixIL, ixOL, wCT, w, x)
Definition: mod_mf_phys.t:911
subroutine mf_check_params
Definition: mod_mf_phys.t:315
subroutine mf_add_source_geom(qdt, dtfactor, ixIL, ixOL, wCT, wprim, w, x)
Definition: mod_mf_phys.t:1166
subroutine add_source_janhunen(qdt, ixIL, ixOL, wCT, w, x)
Definition: mod_mf_phys.t:974
subroutine, public mf_face_to_center(ixOL, s)
calculate cell-center values from face-center values
Definition: mod_mf_phys.t:2309
subroutine mf_modify_wlr(ixIL, ixOL, qt, wLC, wRC, wLp, wRp, s, idir)
Definition: mod_mf_phys.t:1240
logical, public, protected mf_record_electric_field
set to true if need to record electric field on cell edges
Definition: mod_mf_phys.t:73
subroutine add_source_linde(qdt, ixIL, ixOL, wCT, w, x)
Definition: mod_mf_phys.t:995
subroutine update_faces_hll(ixIL, ixOL, qt, qdt, fE, sCT, s, vcts)
update faces
Definition: mod_mf_phys.t:2093
subroutine, public b_from_vector_potential(ixIsL, ixIL, ixOL, ws, x)
calculate magnetic field from vector potential
Definition: mod_mf_phys.t:2371
logical, public, protected clean_initial_divb
clean divb in the initial condition
Definition: mod_mf_phys.t:88
subroutine fixdivb_boundary(ixGL, ixOL, w, x, iB)
Definition: mod_mf_phys.t:1310
logical, public divbwave
Add divB wave in Roe solver.
Definition: mod_mf_phys.t:85
subroutine mf_get_cmax(w, x, ixIL, ixOL, idim, cmax)
Calculate cmax_idim=csound+abs(v_idim) within ixO^L.
Definition: mod_mf_phys.t:400
double precision, public mf_vmax
maximal limit of magnetofrictional velocity in cm s^-1 (Pomoell 2019 A&A)
Definition: mod_mf_phys.t:15
logical, public, protected mf_glm
Whether GLM-MHD is used.
Definition: mod_mf_phys.t:64
subroutine get_resistive_electric_field(ixIL, ixOL, sCT, s, jce)
calculate eta J at cell edges
Definition: mod_mf_phys.t:2240
double precision, public mf_eta_hyper
The hyper-resistivity.
Definition: mod_mf_phys.t:28
subroutine add_source_res2(qdt, ixIL, ixOL, wCT, w, x)
Add resistive source to w within ixO Uses 5 point stencil (2 neighbours) in each direction,...
Definition: mod_mf_phys.t:816
subroutine, public get_normalized_divb(w, ixIL, ixOL, divb)
get dimensionless div B = |divB| * volume / area / |B|
Definition: mod_mf_phys.t:1064
subroutine add_source_res1(qdt, ixIL, ixOL, wCT, w, x)
Add resistive source to w within ixO Uses 3 point stencil (1 neighbour) in each direction,...
Definition: mod_mf_phys.t:722
subroutine, public get_current(w, ixIL, ixOL, idirmin, current, fourthorder)
Calculate idirmin and the idirmin:3 components of the common current array make sure that dxlevel(^D)...
Definition: mod_mf_phys.t:1098
subroutine mf_get_cbounds(wLC, wRC, wLp, wRp, x, ixIL, ixOL, idim, Hspeed, cmax, cmin)
Estimating bounds for the minimum and maximum signal velocities.
Definition: mod_mf_phys.t:412
double precision function integral_grid_mf(ixIL, ixOL, w, x, iw, patchwi)
Definition: mod_mf_phys.t:2509
subroutine mf_velocity_update(qt, psa)
Add global source terms to update frictional velocity on complete domain.
Definition: mod_mf_phys.t:523
double precision, public mf_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
Definition: mod_mf_phys.t:22
subroutine mf_add_source(qdt, dtfactor, ixIL, ixOL, wCT, wCTprim, w, x, qsourcesplit, active)
w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
Definition: mod_mf_phys.t:575
subroutine mf_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
If resistivity is not zero, check diffusion time limit for dt.
Definition: mod_mf_phys.t:1122
subroutine mask_inner(ixIL, ixOL, w, x, patchwi, volume_pe)
Definition: mod_mf_phys.t:2479
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
Definition: mod_mf_phys.t:79
subroutine update_faces_average(ixIL, ixOL, qt, qdt, fC, fE, sCT, s)
get electric field though averaging neighors to update faces in CT
Definition: mod_mf_phys.t:1841
character(len=std_len), public, protected type_ct
Method type of constrained transport.
Definition: mod_mf_phys.t:94
character(len=std_len), public, protected typedivbfix
Method type to clean divergence of B.
Definition: mod_mf_phys.t:91
subroutine, public mf_to_conserved(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
Definition: mod_mf_phys.t:380
subroutine mf_boundary_adjust(igrid, psb)
Definition: mod_mf_phys.t:1269
integer, dimension(2 *^nd), public, protected boundary_divbfix_skip
To skip * layer of ghost cells during divB=0 fix for boundary.
Definition: mod_mf_phys.t:46
subroutine add_source_powel(qdt, ixIL, ixOL, wCT, w, x)
Add divB related sources to w within ixO corresponding to Powel.
Definition: mod_mf_phys.t:954
double precision, dimension(2 *^nd), public mf_decay_scale
decay scale of frictional velocity near boundaries
Definition: mod_mf_phys.t:18
subroutine, public record_force_free_metrics()
Definition: mod_mf_phys.t:2385
subroutine mf_write_info(fh)
Write this module's parameters to a snapsoht.
Definition: mod_mf_phys.t:141
logical, public, protected mf_4th_order
MHD fourth order.
Definition: mod_mf_phys.t:70
integer, public, protected psi_
Indices of the GLM psi.
Definition: mod_mf_phys.t:43
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
Definition: mod_mf_phys.t:40
subroutine, public mf_clean_divb_multigrid(qdt, qt, active)
Definition: mod_mf_phys.t:1675
logical, public, protected mf_particles
Whether particles module is added.
Definition: mod_mf_phys.t:61
double precision, public mf_eta
The resistivity.
Definition: mod_mf_phys.t:25
subroutine update_faces_contact(ixIL, ixOL, qt, qdt, wp, fC, fE, sCT, s, vcts)
update faces using UCT contact mode by Gardiner and Stone 2005 JCP 205, 509
Definition: mod_mf_phys.t:1931
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
Definition: mod_mf_phys.t:67
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
Definition: mod_mf_phys.t:34
subroutine mf_update_faces(ixIL, ixOL, qt, qdt, wprim, fC, fE, sCT, s, vcts)
Definition: mod_mf_phys.t:1815
subroutine mf_physical_units()
Definition: mod_mf_phys.t:320
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
Module containing all the particle routines.
Definition: mod_particles.t:2
subroutine particles_init()
Initialize particle data and parameters.
Definition: mod_particles.t:15
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
procedure(sub_get_ct_velocity), pointer phys_get_ct_velocity
Definition: mod_physics.t:83
procedure(sub_convert), pointer phys_to_primitive
Definition: mod_physics.t:59
integer, parameter flux_tvdlf
Indicates the flux should be treated with tvdlf.
Definition: mod_physics.t:26
procedure(sub_write_info), pointer phys_write_info
Definition: mod_physics.t:81
procedure(sub_get_flux), pointer phys_get_flux
Definition: mod_physics.t:68
logical phys_req_diagonal
Whether the physics routines require diagonal ghost cells, for example for computing a curl.
Definition: mod_physics.t:36
procedure(sub_get_cbounds), pointer phys_get_cbounds
Definition: mod_physics.t:67
procedure(sub_get_dt), pointer phys_get_dt
Definition: mod_physics.t:71
procedure(sub_add_source_geom), pointer phys_add_source_geom
Definition: mod_physics.t:72
procedure(sub_check_params), pointer phys_check_params
Definition: mod_physics.t:56
integer, parameter flux_default
Indicates a normal flux.
Definition: mod_physics.t:24
integer, dimension(:, :), allocatable flux_type
Array per direction per variable, which can be used to specify that certain fluxes have to be treated...
Definition: mod_physics.t:21
procedure(sub_update_faces), pointer phys_update_faces
Definition: mod_physics.t:84
procedure(sub_convert), pointer phys_to_conserved
Definition: mod_physics.t:58
character(len=name_len) physics_type
String describing the physics type of the simulation.
Definition: mod_physics.t:54
procedure(sub_clean_divb), pointer phys_clean_divb
Definition: mod_physics.t:88
procedure(sub_boundary_adjust), pointer phys_boundary_adjust
Definition: mod_physics.t:80
procedure(sub_global_source), pointer phys_global_source_after
Definition: mod_physics.t:74
procedure(sub_add_source), pointer phys_add_source
Definition: mod_physics.t:73
procedure(sub_face_to_center), pointer phys_face_to_center
Definition: mod_physics.t:85
procedure(sub_modify_wlr), pointer phys_modify_wlr
Definition: mod_physics.t:60
procedure(sub_get_cmax), pointer phys_get_cmax
Definition: mod_physics.t:61
logical phys_energy
Solve energy equation or not.
Definition: mod_physics.t:39
procedure(sub_special_advance), pointer phys_special_advance
Definition: mod_physics.t:75
Module with all the methods that users can customize in AMRVAC.
procedure(special_resistivity), pointer usr_special_resistivity
procedure(set_electric_field), pointer usr_set_electric_field
The data structure that contains information about a tree node/grid block.
Definition: mod_forest.t:11