MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_mhd_phys.t
Go to the documentation of this file.
1 !> Magneto-hydrodynamics module
3  use mod_global_parameters, only: std_len
4  implicit none
5  private
6 
7  !> Whether an energy equation is used
8  logical, public, protected :: mhd_energy = .true.
9 
10  !> Whether thermal conduction is used
11  logical, public, protected :: mhd_thermal_conduction = .false.
12 
13  !> Whether radiative cooling is added
14  logical, public, protected :: mhd_radiative_cooling = .false.
15 
16  !> Whether viscosity is added
17  logical, public, protected :: mhd_viscosity = .false.
18 
19  !> Whether gravity is added
20  logical, public, protected :: mhd_gravity = .false.
21 
22  !> Whether Hall-MHD is used
23  logical, public, protected :: mhd_hall = .false.
24 
25  !> Whether particles module is added
26  logical, public, protected :: mhd_particles = .false.
27 
28  !> Whether magnetofriction is added
29  logical, public, protected :: mhd_magnetofriction = .false.
30 
31  !> Whether GLM-MHD is used
32  logical, public, protected :: mhd_glm = .false.
33 
34  !> Whether auxiliary internal energy is solved
35  logical, public, protected :: mhd_solve_eaux = .false.
36 
37  !> Whether internal energy is solved instead of total energy
38  logical, public, protected :: mhd_internal_e = .false.
39 
40  !> Whether divB cleaning sources are added splitting from fluid solver
41  logical, public, protected :: source_split_divb = .false.
42 
43  !> GLM-MHD parameter: ratio of the diffusive and advective time scales for div b
44  !> taking values within [0, 1]
45  double precision, public :: mhd_glm_alpha = 0.5d0
46 
47  !> Use Boris approximation
48  character(len=20) :: mhd_boris_method = "none"
49 
50  integer, parameter :: boris_none = 0
51  integer, parameter :: boris_reduced_force = 1
52  integer, parameter :: boris_simplification = 2
53  integer :: mhd_boris_type = boris_none
54 
55  !> Speed of light for Boris' approximation. If negative, test changes to the
56  !> momentum equation with gamma_A = 1
57  double precision :: mhd_boris_c = 0.0d0
58 
59  !> MHD fourth order
60  logical, public, protected :: mhd_4th_order = .false.
61 
62  !> Number of tracer species
63  integer, public, protected :: mhd_n_tracer = 0
64 
65  !> Index of the density (in the w array)
66  integer, public, protected :: rho_
67 
68  !> Indices of the momentum density
69  integer, allocatable, public, protected :: mom(:)
70 
71  !> Index of the energy density (-1 if not present)
72  integer, public, protected :: e_
73 
74  !> Index of the gas pressure (-1 if not present) should equal e_
75  integer, public, protected :: p_
76 
77  !> Indices of the magnetic field
78  integer, allocatable, public, protected :: mag(:)
79 
80  !> Indices of the GLM psi
81  integer, public, protected :: psi_
82 
83  !> Indices of auxiliary internal energy
84  integer, public, protected :: eaux_
85  integer, public, protected :: paux_
86 
87  !> Indices of the tracers
88  integer, allocatable, public, protected :: tracer(:)
89 
90  !> The adiabatic index
91  double precision, public :: mhd_gamma = 5.d0/3.0d0
92 
93  !> The adiabatic constant
94  double precision, public :: mhd_adiab = 1.0d0
95 
96  !> The MHD resistivity
97  double precision, public :: mhd_eta = 0.0d0
98 
99  !> The MHD hyper-resistivity
100  double precision, public :: mhd_eta_hyper = 0.0d0
101 
102  !> TODO: what is this?
103  double precision, public :: mhd_etah = 0.0d0
104 
105  !> The small_est allowed energy
106  double precision, protected :: small_e
107 
108  !> The number of waves
109  integer :: nwwave=8
110 
111  !> Method type to clean divergence of B
112  character(len=std_len), public, protected :: typedivbfix = 'linde'
113 
114  !> Method type of constrained transport
115  character(len=std_len), public, protected :: type_ct = 'uct_contact'
116 
117  !> Whether divB is computed with a fourth order approximation
118  logical, public, protected :: mhd_divb_4thorder = .false.
119 
120  !> Method type in a integer for good performance
121  integer :: type_divb
122 
123  !> Coefficient of diffusive divB cleaning
124  double precision :: divbdiff = 0.8d0
125 
126  !> Update all equations due to divB cleaning
127  character(len=std_len) :: typedivbdiff = 'all'
128 
129  !> Use a compact way to add resistivity
130  logical :: compactres = .false.
131 
132  !> Add divB wave in Roe solver
133  logical, public :: divbwave = .true.
134 
135  !> clean initial divB
136  logical, public :: clean_initial_divb = .false.
137 
138  !> Helium abundance over Hydrogen
139  double precision, public, protected :: he_abundance=0.1d0
140 
141  !> To control divB=0 fix for boundary
142  logical, public, protected :: boundary_divbfix(2*^nd)=.true.
143 
144  !> To skip * layer of ghost cells during divB=0 fix for boundary
145  integer, public, protected :: boundary_divbfix_skip(2*^nd)=0
146 
147  !> B0 field is force-free
148  logical, public, protected :: b0field_forcefree=.true.
149 
150  !> Whether an total energy equation is used
151  logical :: total_energy = .true.
152 
153  !> gamma minus one and its inverse
154  double precision :: gamma_1, inv_gamma_1
155 
156  ! DivB cleaning methods
157  integer, parameter :: divb_none = 0
158  integer, parameter :: divb_multigrid = -1
159  integer, parameter :: divb_glm = 1
160  integer, parameter :: divb_powel = 2
161  integer, parameter :: divb_janhunen = 3
162  integer, parameter :: divb_linde = 4
163  integer, parameter :: divb_lindejanhunen = 5
164  integer, parameter :: divb_lindepowel = 6
165  integer, parameter :: divb_lindeglm = 7
166  integer, parameter :: divb_ct = 8
167 
168 
169  ! Public methods
170  public :: mhd_phys_init
171  public :: mhd_kin_en
172  public :: mhd_get_pthermal
173  public :: mhd_get_v
174  public :: mhd_get_v_idim
175  public :: mhd_to_conserved
176  public :: mhd_to_primitive
177  public :: mhd_get_csound2
178  public :: mhd_face_to_center
179  public :: get_divb
180  public :: get_current
181  public :: get_normalized_divb
182  public :: b_from_vector_potential
183  public :: mhd_mag_en_all
184  {^nooned
185  public :: mhd_clean_divb_multigrid
186  }
187 
188 contains
189 
190  !> Read this module"s parameters from a file
191  subroutine mhd_read_params(files)
193  character(len=*), intent(in) :: files(:)
194  integer :: n
195 
196  namelist /mhd_list/ mhd_energy, mhd_n_tracer, mhd_gamma, mhd_adiab,&
200  typedivbdiff, type_ct, compactres, divbwave, he_abundance, si_unit, b0field,&
203  mhd_boris_method, mhd_boris_c, clean_initial_divb, mhd_solve_eaux, mhd_internal_e
204 
205  do n = 1, size(files)
206  open(unitpar, file=trim(files(n)), status="old")
207  read(unitpar, mhd_list, end=111)
208 111 close(unitpar)
209  end do
210 
211  end subroutine mhd_read_params
212 
213  !> Write this module's parameters to a snapsoht
214  subroutine mhd_write_info(fh)
216  integer, intent(in) :: fh
217  integer, parameter :: n_par = 1
218  double precision :: values(n_par)
219  character(len=name_len) :: names(n_par)
220  integer, dimension(MPI_STATUS_SIZE) :: st
221  integer :: er
222 
223  call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
224 
225  names(1) = "gamma"
226  values(1) = mhd_gamma
227  call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
228  call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
229  end subroutine mhd_write_info
230 
231  subroutine mhd_angmomfix(fC,x,wnew,ixI^L,ixO^L,idim)
233  double precision, intent(in) :: x(ixi^s,1:ndim)
234  double precision, intent(inout) :: fC(ixi^s,1:nwflux,1:ndim), wnew(ixi^s,1:nw)
235  integer, intent(in) :: ixI^L, ixO^L
236  integer, intent(in) :: idim
237  integer :: hxO^L, kxC^L, iw
238  double precision :: inv_volume(ixi^s)
239 
240  call mpistop("to do")
241 
242  end subroutine mhd_angmomfix
243 
244  subroutine mhd_phys_init()
248  use mod_viscosity, only: viscosity_init
249  use mod_gravity, only: gravity_init
250  use mod_particles, only: particles_init
252  use mod_physics
253  {^nooned
255  }
256 
257  integer :: itr, idir
258 
259  call mhd_read_params(par_files)
260 
261  physics_type = "mhd"
265 
266  if(mhd_energy.and..not.mhd_internal_e) then
267  total_energy=.true.
268  else
269  total_energy=.false.
270  end if
271  phys_total_energy=total_energy
272 
273  if(.not. mhd_energy) then
274  if(mhd_internal_e) then
275  mhd_internal_e=.false.
276  if(mype==0) write(*,*) 'WARNING: set mhd_internal_e=F when mhd_energy=F'
277  end if
278  if(mhd_solve_eaux) then
279  mhd_solve_eaux=.false.
280  if(mype==0) write(*,*) 'WARNING: set mhd_solve_eaux=F when mhd_energy=F'
281  end if
282  if(mhd_thermal_conduction) then
283  mhd_thermal_conduction=.false.
284  if(mype==0) write(*,*) 'WARNING: set mhd_thermal_conduction=F when mhd_energy=F'
285  end if
286  if(mhd_radiative_cooling) then
287  mhd_radiative_cooling=.false.
288  if(mype==0) write(*,*) 'WARNING: set mhd_radiative_cooling=F when mhd_energy=F'
289  end if
290  end if
291 
293 
294  ! set default gamma for polytropic/isothermal process
295  if(.not.mhd_energy) mhd_gamma=1.d0
297  if(ndim==1) typedivbfix='none'
298  select case (typedivbfix)
299  case ('none')
300  type_divb = divb_none
301  {^nooned
302  case ('multigrid')
303  type_divb = divb_multigrid
304  use_multigrid = .true.
305  mg%operator_type = mg_laplacian
307  }
308  case ('glm')
309  mhd_glm = .true.
310  need_global_cmax = .true.
311  type_divb = divb_glm
312  case ('powel', 'powell')
313  type_divb = divb_powel
314  case ('janhunen')
315  type_divb = divb_janhunen
316  case ('linde')
317  type_divb = divb_linde
318  case ('lindejanhunen')
319  type_divb = divb_lindejanhunen
320  case ('lindepowel')
321  type_divb = divb_lindepowel
322  case ('lindeglm')
323  mhd_glm = .true.
324  need_global_cmax = .true.
325  type_divb = divb_lindeglm
326  case ('ct')
327  type_divb = divb_ct
328  stagger_grid = .true.
329  case default
330  call mpistop('Unknown divB fix')
331  end select
332 
333  select case (mhd_boris_method)
334  case ("none")
335  mhd_boris_type = boris_none
336  case ("reduced_force")
337  mhd_boris_type = boris_reduced_force
338  case ("simplification")
339  mhd_boris_type = boris_simplification
340  case default
341  call mpistop("Unknown mhd_boris_method (none, reduced_force, simplification)")
342  end select
343 
344  ! Determine flux variables
345  rho_ = var_set_rho()
346 
347  allocate(mom(ndir))
348  mom(:) = var_set_momentum(ndir)
349 
350  ! Set index of energy variable
351  if (mhd_energy) then
352  nwwave = 8
353  e_ = var_set_energy() ! energy density
354  p_ = e_ ! gas pressure
355  else
356  nwwave = 7
357  e_ = -1
358  p_ = -1
359  end if
360 
361  allocate(mag(ndir))
362  mag(:) = var_set_bfield(ndir)
363 
364  if (mhd_glm) then
365  psi_ = var_set_fluxvar('psi', 'psi', need_bc=.false.)
366  else
367  psi_ = -1
368  end if
369 
370  ! set auxiliary internal energy variable
371  if(mhd_energy .and. mhd_solve_eaux) then
372  eaux_ = var_set_internal_energy()
373  paux_ = eaux_
374  else
375  eaux_ = -1
376  paux_ = -1
377  end if
378 
379  allocate(tracer(mhd_n_tracer))
380 
381  ! Set starting index of tracers
382  do itr = 1, mhd_n_tracer
383  tracer(itr) = var_set_fluxvar("trc", "trp", itr, need_bc=.false.)
384  end do
385 
386  ! determine number of stagger variables
387  if(stagger_grid) nws=ndim
388 
389  nvector = 2 ! No. vector vars
390  allocate(iw_vector(nvector))
391  iw_vector(1) = mom(1) - 1 ! TODO: why like this?
392  iw_vector(2) = mag(1) - 1 ! TODO: why like this?
393 
394  ! Check whether custom flux types have been defined
395  if (.not. allocated(flux_type)) then
396  allocate(flux_type(ndir, nw))
398  else if (any(shape(flux_type) /= [ndir, nw])) then
399  call mpistop("phys_check error: flux_type has wrong shape")
400  end if
401 
402  if(ndim>1) then
403  if(mhd_glm) then
405  do idir=1,ndir
406  flux_type(idir,mag(idir))=flux_special
407  end do
408  else
409  do idir=1,ndir
410  flux_type(idir,mag(idir))=flux_tvdlf
411  end do
412  end if
413  end if
414 
415  select case (mhd_boris_method)
416  case ("none")
417  mhd_boris_type = boris_none
418  case ("reduced_force")
419  mhd_boris_type = boris_reduced_force
420  case ("simplification")
421  mhd_boris_type = boris_simplification
422  do idir = 1, ndir
423  phys_iw_methods(mom(idir))%inv_capacity => mhd_gamma2_alfven
424  end do
425  case default
426  call mpistop("Unknown mhd_boris_method (none, reduced_force, simplification)")
427  end select
428 
449 
450  if(type_divb==divb_glm) then
452  end if
453 
454  ! if using ct stagger grid, boundary divb=0 is not done here
455  if(stagger_grid) then
459  else if(ndim>1) then
461  end if
462 
463  ! Whether diagonal ghost cells are required for the physics
464  if(type_divb < divb_linde) phys_req_diagonal = .false.
465 
466  ! derive units from basic units
467  call mhd_physical_units()
468 
469  if(.not. mhd_energy .and. mhd_thermal_conduction) then
470  call mpistop("thermal conduction needs mhd_energy=T")
471  end if
472  if(.not. mhd_energy .and. mhd_radiative_cooling) then
473  call mpistop("radiative cooling needs mhd_energy=T")
474  end if
475 
476  ! initialize thermal conduction module
477  if (mhd_thermal_conduction) then
478  phys_req_diagonal = .true.
480  end if
481 
482  ! Initialize radiative cooling module
483  if (mhd_radiative_cooling) then
485  end if
486 
487  ! Initialize viscosity module
489 
490  ! Initialize gravity module
491  if(mhd_gravity) then
492  call gravity_init()
493  end if
494 
495  ! Initialize particles module
496  if(mhd_particles) then
497  call particles_init()
498  phys_req_diagonal = .true.
499  end if
500 
501  ! initialize magnetofriction module
502  if(mhd_magnetofriction) then
503  phys_req_diagonal = .true.
504  call magnetofriction_init()
505  end if
506 
507  ! For Hall, we need one more reconstructed layer since currents are computed
508  ! in getflux: assuming one additional ghost layer (two for FOURTHORDER) was
509  ! added in nghostcells.
510  if (mhd_hall) then
511  phys_req_diagonal = .true.
512  if (mhd_4th_order) then
514  else
516  end if
517  end if
518 
519  end subroutine mhd_phys_init
520 
521  subroutine mhd_check_params
523 
524  ! after user parameter setting
525  gamma_1=mhd_gamma-1.d0
526  if (.not. mhd_energy) then
527  if (mhd_gamma <= 0.0d0) call mpistop ("Error: mhd_gamma <= 0")
528  if (mhd_adiab < 0.0d0) call mpistop ("Error: mhd_adiab < 0")
530  else
531  if (mhd_gamma <= 0.0d0 .or. mhd_gamma == 1.0d0) &
532  call mpistop ("Error: mhd_gamma <= 0 or mhd_gamma == 1")
533  inv_gamma_1=1.d0/gamma_1
534  small_e = small_pressure * inv_gamma_1
535  end if
536 
537  if (mhd_boris_type > 0 .and. abs(mhd_boris_c) <= 0.0d0) then
538  call mpistop("You have not specified mhd_boris_c")
539  end if
540 
541  end subroutine mhd_check_params
542 
543  subroutine mhd_physical_units()
545  double precision :: mp,kB,miu0,c_lightspeed
546  ! Derive scaling units
547  if(si_unit) then
548  mp=mp_si
549  kb=kb_si
550  miu0=miu0_si
551  c_lightspeed=c_si
552  else
553  mp=mp_cgs
554  kb=kb_cgs
555  miu0=4.d0*dpi
556  c_lightspeed=const_c
557  end if
558  if(unit_velocity==0) then
564  else
570  end if
571  ! Additional units needed for the particles
572  c_norm=c_lightspeed/unit_velocity
574  if (.not. si_unit) unit_charge = unit_charge*const_c
576 
577  end subroutine mhd_physical_units
578 
579  subroutine mhd_check_w(primitive,ixI^L,ixO^L,w,flag,smallw)
581 
582  logical, intent(in) :: primitive
583  integer, intent(in) :: ixI^L, ixO^L
584  double precision, intent(in) :: w(ixi^s,nw)
585  integer, intent(inout) :: flag(ixi^s)
586  double precision, intent(out) :: smallw(1:nw)
587  double precision :: tmp(ixi^s)
588 
589  smallw=1.d0
590  flag(ixo^s)=0
591  where(w(ixo^s, rho_) < small_density) flag(ixo^s) = rho_
592  if(any(flag(ixo^s)==rho_)) smallw(rho_)=minval(w(ixo^s,rho_))
593 
594  if(mhd_energy) then
595  if(primitive) then
596  where(w(ixo^s,e_) < small_pressure) flag(ixo^s) = e_
597  if(any(flag(ixo^s)==e_)) smallw(e_)=minval(w(ixo^s,e_))
598  else
599  if(mhd_internal_e) then
600  where(w(ixo^s, e_) < small_pressure*inv_gamma_1) flag(ixo^s) = e_
601  if(any(flag(ixo^s)==e_)) smallw(e_)=minval(w(ixo^s,e_))
602  else
603  ! Calculate pressure=(gamma-1)*(e-0.5*(2ek+2eb))
604  tmp(ixo^s)=w(ixo^s,e_)-&
605  mhd_kin_en(w,ixi^l,ixo^l)-mhd_mag_en(w,ixi^l,ixo^l)
606  where(tmp(ixo^s) < small_pressure*inv_gamma_1) flag(ixo^s) = e_
607  if(any(flag(ixo^s)==e_)) smallw(e_)=minval(tmp(ixo^s))
608  end if
609  end if
610  end if
611 
612  end subroutine mhd_check_w
613 
614  !> Transform primitive variables into conservative ones
615  subroutine mhd_to_conserved(ixI^L,ixO^L,w,x)
617  integer, intent(in) :: ixI^L, ixO^L
618  double precision, intent(inout) :: w(ixi^s, nw)
619  double precision, intent(in) :: x(ixi^s, 1:ndim)
620  integer :: idir, itr
621 
623  call mhd_handle_small_values(.true., w, x, ixi^l, ixo^l, 'mhd_to_conserved')
624  end if
625 
626  ! Calculate total energy from pressure, kinetic and magnetic energy
627  if(mhd_energy) then
628  if(mhd_internal_e) then
629  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1
630  else
631  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1&
632  +half*sum(w(ixo^s,mom(:))**2,dim=ndim+1)*w(ixo^s,rho_)&
633  +mhd_mag_en(w, ixi^l, ixo^l)
634  if(mhd_solve_eaux) w(ixo^s,eaux_)=w(ixo^s,paux_)*inv_gamma_1
635  end if
636  end if
637 
638  ! Convert velocity to momentum
639  do idir = 1, ndir
640  w(ixo^s, mom(idir)) = w(ixo^s, rho_) * w(ixo^s, mom(idir))
641  end do
642 
643  if (check_small_values .and. .not. small_values_use_primitive) then
644  call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_conserved')
645  end if
646  end subroutine mhd_to_conserved
647 
648  !> Transform conservative variables into primitive ones
649  subroutine mhd_to_primitive(ixI^L,ixO^L,w,x)
651  integer, intent(in) :: ixI^L, ixO^L
652  double precision, intent(inout) :: w(ixi^s, nw)
653  double precision, intent(in) :: x(ixi^s, 1:ndim)
654  double precision :: inv_rho(ixo^s)
655  integer :: itr, idir
656 
657  if (check_small_values .and. .not. small_values_use_primitive) then
658  call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive')
659  end if
660 
661  inv_rho=1.0d0/w(ixo^s,rho_)
662 
663  ! Calculate pressure = (gamma-1) * (e-ek-eb)
664  if(mhd_energy) then
665  if(mhd_internal_e) then
666  w(ixo^s,p_)=w(ixo^s,e_)*gamma_1
667  else
668  w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)&
669  -mhd_kin_en(w,ixi^l,ixo^l,inv_rho)&
670  -mhd_mag_en(w,ixi^l,ixo^l))
671  if(mhd_solve_eaux) w(ixo^s,paux_)=w(ixo^s,eaux_)*gamma_1
672  end if
673  end if
674 
675  ! Convert momentum to velocity
676  do idir = 1, ndir
677  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))*inv_rho
678  end do
679 
681  call mhd_handle_small_values(.true., w, x, ixi^l, ixo^l, 'mhd_to_primitive')
682  end if
683  end subroutine mhd_to_primitive
684 
685  !> Transform internal energy to total energy
686  subroutine mhd_ei_to_e(ixI^L,ixO^L,w,x)
688  integer, intent(in) :: ixI^L, ixO^L
689  double precision, intent(inout) :: w(ixi^s, nw)
690  double precision, intent(in) :: x(ixi^s, 1:ndim)
691 
692  ! Calculate total energy from internal, kinetic and magnetic energy
693  if(total_energy) then
694  w(ixo^s,e_)=w(ixo^s,e_)&
695  +mhd_kin_en(w,ixi^l,ixo^l)&
696  +mhd_mag_en(w,ixi^l,ixo^l)
697  end if
698 
699  end subroutine mhd_ei_to_e
700 
701  !> Transform total energy to internal energy
702  subroutine mhd_e_to_ei(ixI^L,ixO^L,w,x)
704  integer, intent(in) :: ixI^L, ixO^L
705  double precision, intent(inout) :: w(ixi^s, nw)
706  double precision, intent(in) :: x(ixi^s, 1:ndim)
707 
708  ! Calculate ei = e - ek - eb
709  if(total_energy) then
710  w(ixo^s,e_)=w(ixo^s,e_)&
711  -mhd_kin_en(w,ixi^l,ixo^l)&
712  -mhd_mag_en(w,ixi^l,ixo^l)
713  end if
714 
715  end subroutine mhd_e_to_ei
716 
717  subroutine mhd_energy_synchro(ixI^L,ixO^L,w,x)
719  integer, intent(in) :: ixI^L,ixO^L
720  double precision, intent(in) :: x(ixi^s,1:ndim)
721  double precision, intent(inout) :: w(ixi^s,1:nw)
722 
723  double precision :: pth1(ixi^s),pth2(ixi^s),alfa(ixi^s),beta(ixi^s)
724  double precision, parameter :: beta_low=0.005d0,beta_high=0.05d0
725 ! double precision :: vtot(ixI^S),cs2(ixI^S),mach(ixI^S)
726 ! double precision, parameter :: mach_low=20.d0,mach_high=200.d0
727 
728  ! get magnetic energy
729  alfa(ixo^s)=mhd_mag_en(w,ixi^l,ixo^l)
730  pth1(ixo^s)=gamma_1*(w(ixo^s,e_)-mhd_kin_en(w,ixi^l,ixo^l)-alfa(ixo^s))
731  pth2(ixo^s)=w(ixo^s,eaux_)*gamma_1
732  ! get plasma beta
733  beta(ixo^s)=min(pth1(ixo^s),pth2(ixo^s))/alfa(ixo^s)
734 
735  ! whether Mach number should be another criterion ?
736 ! vtot(ixO^S)=sum(w(ixO^S,mom(:))**2,dim=ndim+1)
737 ! call mhd_get_csound2(w,x,ixI^L,ixO^L,cs2)
738 ! mach(ixO^S)=sqrt(vtot(ixO^S)/cs2(ixO^S))/w(ixO^S,rho_)
739  where(beta(ixo^s) .ge. beta_high)
740 ! where(beta(ixO^S) .ge. beta_high .and. mach(ixO^S) .le. mach_low)
741  w(ixo^s,eaux_)=pth1(ixo^s)*inv_gamma_1
742  else where(beta(ixo^s) .le. beta_low)
743 ! else where(beta(ixO^S) .le. beta_low .or. mach(ixO^S) .ge. mach_high)
744  w(ixo^s,e_)=w(ixo^s,e_)-pth1(ixo^s)*inv_gamma_1+w(ixo^s,eaux_)
745  else where
746  alfa(ixo^s)=dlog(beta(ixo^s)/beta_low)/dlog(beta_high/beta_low)
747 ! alfa(ixO^S)=min(dlog(beta(ixO^S)/beta_low)/dlog(beta_high/beta_low),
748 ! dlog(mach_high(ixO^S)/mach(ixO^S))/dlog(mach_high/mach_low))
749  w(ixo^s,eaux_)=(pth2(ixo^s)*(one-alfa(ixo^s))&
750  +pth1(ixo^s)*alfa(ixo^s))*inv_gamma_1
751  w(ixo^s,e_)=w(ixo^s,e_)-pth1(ixo^s)*inv_gamma_1+w(ixo^s,eaux_)
752  end where
753  end subroutine mhd_energy_synchro
754 
755  subroutine mhd_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
757  use mod_small_values
758  logical, intent(in) :: primitive
759  integer, intent(in) :: ixI^L,ixO^L
760  double precision, intent(inout) :: w(ixi^s,1:nw)
761  double precision, intent(in) :: x(ixi^s,1:ndim)
762  character(len=*), intent(in) :: subname
763 
764  double precision :: smallw(1:nw)
765  integer :: idir, flag(ixi^s)
766 
767  if (small_values_method == "ignore") return
768 
769  call mhd_check_w(primitive, ixi^l, ixo^l, w, flag, smallw)
770 
771  if (any(flag(ixo^s) /= 0)) then
772  select case (small_values_method)
773  case ("replace")
774  if (small_values_fix_iw(rho_)) then
775  where(flag(ixo^s) /= 0) w(ixo^s,rho_) = small_density
776  end if
777 
778  do idir = 1, ndir
779  if (small_values_fix_iw(mom(idir))) then
780  where(flag(ixo^s) /= 0) w(ixo^s, mom(idir)) = 0.0d0
781  end if
782  end do
783 
784  if (mhd_energy) then
785  if (small_values_fix_iw(e_)) then
786  if(mhd_solve_eaux) then
787  if(primitive) then
788  where(flag(ixo^s) /= 0) w(ixo^s,p_)=w(ixo^s,paux_)
789  else
790  where(flag(ixo^s) /= 0)
791  w(ixo^s,e_) = w(ixo^s,eaux_) + 0.5d0 * &
792  sum(w(ixo^s, mom(:))**2, dim=ndim+1) / w(ixo^s, rho_) + &
793  mhd_mag_en(w, ixi^l, ixo^l)
794  end where
795  end if
796  else
797  if(primitive) then
798  where(flag(ixo^s) /= 0) w(ixo^s,p_) = small_pressure
799  else
800  where(flag(ixo^s) /= 0)
801  w(ixo^s,e_) = small_e + 0.5d0 * &
802  sum(w(ixo^s, mom(:))**2, dim=ndim+1) / w(ixo^s, rho_) + &
803  mhd_mag_en(w, ixi^l, ixo^l)
804  end where
805  end if
806  end if
807  end if
808  end if
809  case ("average")
810  call small_values_average(ixi^l, ixo^l, w, x, flag)
811  case default
812  call small_values_error(w, x, ixi^l, ixo^l, flag, subname, smallw)
813  end select
814  end if
815  end subroutine mhd_handle_small_values
816 
817  !> Convert energy to entropy
818  subroutine e_to_rhos(ixI^L,ixO^L,w,x)
820  integer, intent(in) :: ixI^L, ixO^L
821  double precision,intent(inout) :: w(ixi^s,nw)
822  double precision, intent(in) :: x(ixi^s,1:ndim)
823 
824  if(mhd_energy) then
825  if(.not.mhd_internal_e) &
826  w(ixo^s, e_)=w(ixo^s, e_)-mhd_kin_en(w, ixi^l, ixo^l) &
827  -mhd_mag_en(w, ixi^l, ixo^l)
828  w(ixo^s, e_)=gamma_1*w(ixo^s, rho_)**(1.0d0 - mhd_gamma)&
829  *w(ixo^s, e_)
830  else
831  call mpistop("e_to_rhos can not be used without energy equation!")
832  end if
833  end subroutine e_to_rhos
834 
835  !> Convert entropy to energy
836  subroutine rhos_to_e(ixI^L,ixO^L,w,x)
838  integer, intent(in) :: ixI^L, ixO^L
839  double precision :: w(ixi^s,nw)
840  double precision, intent(in) :: x(ixi^s,1:ndim)
841 
842  if(mhd_energy) then
843  w(ixo^s, e_)=w(ixo^s, rho_)**gamma_1 * w(ixo^s, e_) &
844  * inv_gamma_1
845  if(.not.mhd_internal_e) &
846  w(ixo^s, e_)=w(ixo^s, e_)+mhd_kin_en(w, ixi^l, ixo^l) + &
847  mhd_mag_en(w, ixi^l, ixo^l)
848  else
849  call mpistop("rhos_to_e can not be used without energy equation!")
850  end if
851  end subroutine rhos_to_e
852 
853  !> Calculate v vector
854  subroutine mhd_get_v(w,x,ixI^L,ixO^L,v)
856 
857  integer, intent(in) :: ixI^L, ixO^L
858  double precision, intent(in) :: w(ixi^s,nw), x(ixi^s,1:ndim)
859  double precision, intent(out) :: v(ixi^s,ndir)
860 
861  integer :: idir
862 
863  do idir=1,ndir
864  v(ixo^s,idir) = w(ixo^s, mom(idir)) / w(ixo^s,rho_)
865  end do
866 
867  end subroutine mhd_get_v
868 
869  !> Calculate v component
870  subroutine mhd_get_v_idim(w,x,ixI^L,ixO^L,idim,v)
872 
873  integer, intent(in) :: ixI^L, ixO^L, idim
874  double precision, intent(in) :: w(ixi^s,nw), x(ixi^s,1:ndim)
875  double precision, intent(out) :: v(ixi^s)
876 
877  v(ixo^s) = w(ixo^s, mom(idim)) / w(ixo^s,rho_)
878 
879  end subroutine mhd_get_v_idim
880 
881  !> Calculate cmax_idim=csound+abs(v_idim) within ixO^L
882  subroutine mhd_get_cmax(w,x,ixI^L,ixO^L,idim,cmax)
884 
885  integer, intent(in) :: ixI^L, ixO^L, idim
886  double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
887  double precision, intent(inout) :: cmax(ixi^s)
888 
889  call mhd_get_csound(w,x,ixi^l,ixo^l,idim,cmax)
890 
891  cmax(ixo^s)=abs(w(ixo^s,mom(idim))/w(ixo^s,rho_))+cmax(ixo^s)
892 
893  end subroutine mhd_get_cmax
894 
895  subroutine mhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
897 
898  integer, intent(in) :: ixI^L, ixO^L
899  double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
900  double precision, intent(inout) :: a2max(ndim)
901  double precision :: a2(ixi^s,ndim,nw)
902  integer :: gxO^L,hxO^L,jxO^L,kxO^L,i,j
903 
904  a2=zero
905  do i = 1,ndim
906  !> 4th order
907  hxo^l=ixo^l-kr(i,^d);
908  gxo^l=hxo^l-kr(i,^d);
909  jxo^l=ixo^l+kr(i,^d);
910  kxo^l=jxo^l+kr(i,^d);
911  a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
912  -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
913  a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/dxlevel(i)**2
914  end do
915  end subroutine mhd_get_a2max
916 
917  !> get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
918  subroutine mhd_get_tcutoff(ixI^L,ixO^L,w,x,tco_local,Tmax_local)
920  use mod_geometry
921  integer, intent(in) :: ixI^L,ixO^L
922  double precision, intent(in) :: x(ixi^s,1:ndim),w(ixi^s,1:nw)
923  double precision, intent(out) :: tco_local, Tmax_local
924 
925  double precision, parameter :: delta=0.5d0
926  double precision :: tmp1(ixi^s),Te(ixi^s),lts(ixi^s)
927  double precision, dimension(ixI^S,1:ndir) :: bunitvec
928  double precision, dimension(ixI^S,1:ndim) :: gradT
929  integer :: idims
930  logical :: lrlt(ixi^s)
931 
932  if(mhd_internal_e) then
933  tmp1(ixi^s)=w(ixi^s,e_)
934  else
935  tmp1(ixi^s)=w(ixi^s,e_)-0.5d0*(sum(w(ixi^s,iw_mom(:))**2,dim=ndim+1)/&
936  w(ixi^s,rho_)+sum(w(ixi^s,iw_mag(:))**2,dim=ndim+1))
937  end if
938  te(ixi^s)=tmp1(ixi^s)/w(ixi^s,rho_)*(mhd_gamma-1.d0)
939 
940  tmax_local=maxval(te(ixo^s))
941 
942  ! temperature gradient at cell centers
943  do idims=1,ndim
944  call gradient(te,ixi^l,ixo^l,idims,tmp1)
945  gradt(ixo^s,idims)=tmp1(ixo^s)
946  end do
947  ! B vector
948  if(b0field) then
949  bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+block%B0(ixo^s,:,0)
950  else
951  bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
952  end if
953  ! |B|
954  tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
955  where(tmp1(ixo^s)/=0.d0)
956  tmp1(ixo^s)=1.d0/tmp1(ixo^s)
957  elsewhere
958  tmp1(ixo^s)=bigdouble
959  end where
960  ! b unit vector: magnetic field direction vector
961  do idims=1,ndim
962  bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
963  end do
964  ! temperature length scale inversed
965  lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
966  ! fraction of cells size to temperature length scale
967  lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
968  lrlt=.false.
969  where(lts(ixo^s) > delta)
970  lrlt(ixo^s)=.true.
971  end where
972  block%special_values(1)=0.d0
973  if(any(lrlt(ixo^s))) then
974  block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
975  end if
976 
977  end subroutine mhd_get_tcutoff
978 
979  !> Estimating bounds for the minimum and maximum signal velocities
980  subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,cmax,cmin)
983 
984  integer, intent(in) :: ixI^L, ixO^L, idim
985  double precision, intent(in) :: wLC(ixi^s, nw), wRC(ixi^s, nw)
986  double precision, intent(in) :: wLp(ixi^s, nw), wRp(ixi^s, nw)
987  double precision, intent(in) :: x(ixi^s,1:ndim)
988  double precision, intent(inout) :: cmax(ixi^s)
989  double precision, intent(inout), optional :: cmin(ixi^s)
990 
991  double precision :: wmean(ixi^s,nw)
992  double precision, dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
993  integer :: idimE,idimN
994 
995  if (typeboundspeed=='cmaxmean') then
996  wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
997  tmp1(ixo^s)=wmean(ixo^s,mom(idim))/wmean(ixo^s,rho_)
998  call mhd_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
999  if(present(cmin)) then
1000  cmax(ixo^s)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
1001  cmin(ixo^s)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
1002  else
1003  cmax(ixo^s)=abs(tmp1(ixo^s))+csoundr(ixo^s)
1004  end if
1005  else
1006  ! This implements formula (10.52) from "Riemann Solvers and Numerical
1007  ! Methods for Fluid Dynamics" by Toro.
1008  tmp1(ixo^s)=sqrt(abs(wlp(ixo^s,rho_)))
1009  tmp2(ixo^s)=sqrt(abs(wrp(ixo^s,rho_)))
1010  tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
1011  umean(ixo^s)=(wlp(ixo^s,mom(idim))*tmp1(ixo^s)+wrp(ixo^s,mom(idim))*tmp2(ixo^s))*tmp3(ixo^s)
1012  call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
1013  call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
1014  dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
1015  0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
1016  (wrp(ixo^s,mom(idim))-wlp(ixo^s,mom(idim)))**2
1017  dmean(ixo^s)=sqrt(dmean(ixo^s))
1018  if(present(cmin)) then
1019  cmin(ixo^s)=umean(ixo^s)-dmean(ixo^s)
1020  cmax(ixo^s)=umean(ixo^s)+dmean(ixo^s)
1021  else
1022  cmax(ixo^s)=abs(umean(ixo^s))+dmean(ixo^s)
1023  end if
1024  end if
1025 
1026  if(stagger_grid) then
1027  ! calculate velocities related to different UCT schemes
1028  select case(type_ct)
1029  case('average')
1030  case('uct_contact')
1031  if(.not.allocated(vcts%vnorm)) allocate(vcts%vnorm(ixi^s,1:ndim))
1032  ! get average normal velocity at cell faces
1033  vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,mom(idim))+wrp(ixo^s,mom(idim)))
1034  case('uct_hll')
1035  if(.not.allocated(vcts%vbarC)) then
1036  allocate(vcts%vbarC(ixi^s,1:ndir,2),vcts%vbarLC(ixi^s,1:ndir,2),vcts%vbarRC(ixi^s,1:ndir,2))
1037  allocate(vcts%cbarmin(ixi^s,1:ndim),vcts%cbarmax(ixi^s,1:ndim))
1038  end if
1039  ! Store magnitude of characteristics
1040  if(present(cmin)) then
1041  vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
1042  vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
1043  else
1044  vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
1045  vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
1046  end if
1047 
1048  idimn=mod(idim,ndir)+1 ! 'Next' direction
1049  idime=mod(idim+1,ndir)+1 ! Electric field direction
1050  ! Store velocities
1051  vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,mom(idimn))
1052  vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,mom(idimn))
1053  vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
1054  +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
1055  /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
1056 
1057  vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,mom(idime))
1058  vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,mom(idime))
1059  vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
1060  +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
1061  /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
1062  case default
1063  call mpistop('choose average, uct_contact,or uct_hll for type_ct!')
1064  end select
1065  end if
1066 
1067  end subroutine mhd_get_cbounds
1068 
1069  !> Calculate fast magnetosonic wave speed
1070  subroutine mhd_get_csound(w,x,ixI^L,ixO^L,idim,csound)
1072 
1073  integer, intent(in) :: ixI^L, ixO^L, idim
1074  double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
1075  double precision, intent(out):: csound(ixi^s)
1076  double precision :: cfast2(ixi^s), AvMinCs2(ixi^s), b2(ixi^s), kmax
1077  double precision :: inv_rho(ixo^s), gamma2(ixo^s)
1078 
1079  inv_rho=1.d0/w(ixo^s,rho_)
1080 
1081  if (mhd_boris_type == boris_reduced_force) then
1082  call mhd_gamma2_alfven(ixi^l, ixo^l, w, gamma2)
1083  else
1084  gamma2 = 1.0d0
1085  end if
1086 
1087  call mhd_get_csound2(w,x,ixi^l,ixo^l,csound)
1088 
1089  ! store |B|^2 in v
1090  b2(ixo^s) = mhd_mag_en_all(w,ixi^l,ixo^l) * gamma2
1091 
1092  cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
1093  avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
1094  * mhd_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
1095  * inv_rho * gamma2
1096 
1097  where(avmincs2(ixo^s)<zero)
1098  avmincs2(ixo^s)=zero
1099  end where
1100 
1101  avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
1102 
1103  if (.not. mhd_hall) then
1104  csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
1105  if (mhd_boris_type == boris_simplification) then
1106  csound(ixo^s) = mhd_gamma_alfven(w, ixi^l,ixo^l) * csound(ixo^s)
1107  end if
1108  else
1109  ! take the Hall velocity into account:
1110  ! most simple estimate, high k limit:
1111  ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
1112  kmax = dpi/min({dxlevel(^d)},bigdouble)*half
1113  csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
1114  mhd_etah * sqrt(b2(ixo^s))*inv_rho*kmax)
1115  end if
1116 
1117  end subroutine mhd_get_csound
1118 
1119  !> Calculate fast magnetosonic wave speed
1120  subroutine mhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
1122 
1123  integer, intent(in) :: ixI^L, ixO^L, idim
1124  double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
1125  double precision, intent(out):: csound(ixi^s)
1126  double precision :: cfast2(ixi^s), AvMinCs2(ixi^s), b2(ixi^s), kmax
1127  double precision :: inv_rho(ixo^s), gamma_A2(ixo^s)
1128 
1129  inv_rho=1.d0/w(ixo^s,rho_)
1130 
1131  if (mhd_boris_type == boris_reduced_force) then
1132  call mhd_gamma2_alfven(ixi^l, ixo^l, w, gamma_a2)
1133  else
1134  gamma_a2 = 1.0d0
1135  end if
1136 
1137  if(mhd_energy) then
1138  csound(ixo^s)=mhd_gamma*w(ixo^s,p_)*inv_rho
1139  else
1140  csound(ixo^s)=mhd_gamma*mhd_adiab*w(ixo^s,rho_)**gamma_1
1141  end if
1142  ! store |B|^2 in v
1143  b2(ixo^s) = mhd_mag_en_all(w,ixi^l,ixo^l) * gamma_a2
1144  cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
1145  avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
1146  * mhd_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
1147  * inv_rho * gamma_a2
1148 
1149  where(avmincs2(ixo^s)<zero)
1150  avmincs2(ixo^s)=zero
1151  end where
1152 
1153  avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
1154 
1155  if (.not. mhd_hall) then
1156  csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
1157  if (mhd_boris_type == boris_simplification) then
1158  csound(ixo^s) = mhd_gamma_alfven(w, ixi^l,ixo^l) * csound(ixo^s)
1159  end if
1160  else
1161  ! take the Hall velocity into account:
1162  ! most simple estimate, high k limit:
1163  ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
1164  kmax = dpi/min({dxlevel(^d)},bigdouble)*half
1165  csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
1166  mhd_etah * sqrt(b2(ixo^s))*inv_rho*kmax)
1167  end if
1168 
1169  end subroutine mhd_get_csound_prim
1170 
1171  !> Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho-b**2/2) within ixO^L
1172  subroutine mhd_get_pthermal(w,x,ixI^L,ixO^L,pth)
1174 
1175  integer, intent(in) :: ixI^L, ixO^L
1176  double precision, intent(in) :: w(ixi^s,nw)
1177  double precision, intent(in) :: x(ixi^s,1:ndim)
1178  double precision, intent(out):: pth(ixi^s)
1179 
1180  if (check_small_values) then
1181  call mhd_find_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_get_pthermal')
1182  end if
1183 
1184  if(mhd_energy) then
1185  if(mhd_internal_e) then
1186  pth(ixo^s)=gamma_1*w(ixo^s,e_)
1187  else
1188  pth(ixo^s)=gamma_1*(w(ixo^s,e_)&
1189  - mhd_kin_en(w,ixi^l,ixo^l)&
1190  - mhd_mag_en(w,ixi^l,ixo^l))
1191  end if
1192  else
1193  pth(ixo^s)=mhd_adiab*w(ixo^s,rho_)**mhd_gamma
1194  end if
1195  end subroutine mhd_get_pthermal
1196 
1197  subroutine mhd_find_small_values(primitive, w, x, ixI^L, ixO^L, subname)
1199  use mod_small_values
1200  logical, intent(in) :: primitive
1201  integer, intent(in) :: ixI^L,ixO^L
1202  double precision, intent(in) :: w(ixi^s,1:nw)
1203  double precision, intent(in) :: x(ixi^s,1:ndim)
1204  character(len=*), intent(in) :: subname
1205 
1206  double precision :: smallw(1:nw)
1207  integer :: idir, flag(ixi^s)
1208 
1209  call mhd_check_w(primitive, ixi^l, ixo^l, w, flag, smallw)
1210 
1211  if (any(flag(ixo^s) /= 0)) call small_values_error(w, x, ixi^l, ixo^l, flag, subname, smallw)
1212 
1213  end subroutine mhd_find_small_values
1214 
1215  !> Calculate the square of the thermal sound speed csound2 within ixO^L.
1216  !> csound2=gamma*p/rho
1217  subroutine mhd_get_csound2(w,x,ixI^L,ixO^L,csound2)
1219  integer, intent(in) :: ixI^L, ixO^L
1220  double precision, intent(in) :: w(ixi^s,nw)
1221  double precision, intent(in) :: x(ixi^s,1:ndim)
1222  double precision, intent(out) :: csound2(ixi^s)
1223 
1224  if(mhd_energy) then
1225  call mhd_get_pthermal(w,x,ixi^l,ixo^l,csound2)
1226  csound2(ixo^s)=mhd_gamma*csound2(ixo^s)/w(ixo^s,rho_)
1227  else
1228  csound2(ixo^s)=mhd_gamma*mhd_adiab*w(ixo^s,rho_)**gamma_1
1229  end if
1230  end subroutine mhd_get_csound2
1231 
1232  !> Calculate total pressure within ixO^L including magnetic pressure
1233  subroutine mhd_get_p_total(w,x,ixI^L,ixO^L,p)
1235 
1236  integer, intent(in) :: ixI^L, ixO^L
1237  double precision, intent(in) :: w(ixi^s,nw)
1238  double precision, intent(in) :: x(ixi^s,1:ndim)
1239  double precision, intent(out) :: p(ixi^s)
1240 
1241  call mhd_get_pthermal(w,x,ixi^l,ixo^l,p)
1242 
1243  p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
1244 
1245  end subroutine mhd_get_p_total
1246 
1247  !> Calculate fluxes within ixO^L.
1248  subroutine mhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
1250  use mod_geometry
1251 
1252  integer, intent(in) :: ixI^L, ixO^L, idim
1253  ! conservative w
1254  double precision, intent(in) :: wC(ixi^s,nw)
1255  ! primitive w
1256  double precision, intent(in) :: w(ixi^s,nw)
1257  double precision, intent(in) :: x(ixi^s,1:ndim)
1258  double precision,intent(out) :: f(ixi^s,nwflux)
1259 
1260  double precision :: pgas(ixo^s), ptotal(ixo^s),tmp(ixi^s)
1261  double precision, allocatable:: vHall(:^d&,:)
1262  integer :: idirmin, iw, idir, jdir, kdir
1263 
1264  if (mhd_hall) then
1265  allocate(vhall(ixi^s,1:ndir))
1266  call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
1267  end if
1268 
1269  if(b0field) tmp(ixo^s)=sum(block%B0(ixo^s,:,idim)*w(ixo^s,mag(:)),dim=ndim+1)
1270 
1271  if(mhd_energy) then
1272  pgas=w(ixo^s,p_)
1273  else
1274  pgas=mhd_adiab*w(ixo^s,rho_)**mhd_gamma
1275  end if
1276 
1277  ptotal = pgas + 0.5d0*sum(w(ixo^s, mag(:))**2, dim=ndim+1)
1278 
1279  ! Get flux of density
1280  f(ixo^s,rho_)=w(ixo^s,mom(idim))*w(ixo^s,rho_)
1281 
1282  ! Get flux of tracer
1283  do iw=1,mhd_n_tracer
1284  f(ixo^s,tracer(iw))=w(ixo^s,mom(idim))*w(ixo^s,tracer(iw))
1285  end do
1286 
1287  ! Get flux of momentum
1288  ! f_i[m_k]=v_i*m_k-b_k*b_i [+ptotal if i==k]
1289  if (mhd_boris_type == boris_reduced_force) then
1290  do idir=1,ndir
1291  if(idim==idir) then
1292  f(ixo^s,mom(idir)) = pgas(ixo^s)
1293  else
1294  f(ixo^s,mom(idir)) = 0.0d0
1295  end if
1296  f(ixo^s,mom(idir))=f(ixo^s,mom(idir))+w(ixo^s,mom(idim))*wc(ixo^s,mom(idir))
1297  end do
1298  else
1299  ! Normal case (no Boris approximation)
1300  do idir=1,ndir
1301  if(idim==idir) then
1302  f(ixo^s,mom(idir))=ptotal(ixo^s)-w(ixo^s,mag(idim))*w(ixo^s,mag(idir))
1303  if(b0field) f(ixo^s,mom(idir))=f(ixo^s,mom(idir))+tmp(ixo^s)
1304  else
1305  f(ixo^s,mom(idir))= -w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
1306  end if
1307  if (b0field) then
1308  f(ixo^s,mom(idir))=f(ixo^s,mom(idir))&
1309  -w(ixo^s,mag(idir))*block%B0(ixo^s,idim,idim)&
1310  -w(ixo^s,mag(idim))*block%B0(ixo^s,idir,idim)
1311  end if
1312  f(ixo^s,mom(idir))=f(ixo^s,mom(idir))+w(ixo^s,mom(idim))*wc(ixo^s,mom(idir))
1313  end do
1314  end if
1315 
1316  ! Get flux of energy
1317  ! f_i[e]=v_i*e+v_i*ptotal-b_i*(b_k*v_k)
1318  if(mhd_energy) then
1319  if (mhd_internal_e) then
1320  f(ixo^s,e_)=w(ixo^s,mom(idim))*w(ixo^s,p_)
1321  if (mhd_hall) then
1322  call mpistop("solve internal energy not implemented for Hall MHD")
1323  endif
1324  else
1325  f(ixo^s,e_)=w(ixo^s,mom(idim))*(wc(ixo^s,e_)+ptotal(ixo^s))&
1326  -w(ixo^s,mag(idim))*sum(w(ixo^s,mag(:))*w(ixo^s,mom(:)),dim=ndim+1)
1327  if(mhd_solve_eaux) f(ixo^s,eaux_)=w(ixo^s,mom(idim))*wc(ixo^s,eaux_)
1328 
1329  if (b0field) then
1330  f(ixo^s,e_) = f(ixo^s,e_) &
1331  + w(ixo^s,mom(idim)) * tmp(ixo^s) &
1332  - sum(w(ixo^s,mom(:))*w(ixo^s,mag(:)),dim=ndim+1) * block%B0(ixo^s,idim,idim)
1333  end if
1334 
1335  if (mhd_hall) then
1336  ! f_i[e]= f_i[e] + vHall_i*(b_k*b_k) - b_i*(vHall_k*b_k)
1337  if (mhd_etah>zero) then
1338  f(ixo^s,e_) = f(ixo^s,e_) + vhall(ixo^s,idim) * &
1339  sum(w(ixo^s, mag(:))**2,dim=ndim+1) &
1340  - w(ixo^s,mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1)
1341  if (b0field) then
1342  f(ixo^s,e_) = f(ixo^s,e_) &
1343  + vhall(ixo^s,idim) * tmp(ixo^s) &
1344  - sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1) * block%B0(ixo^s,idim,idim)
1345  end if
1346  end if
1347  end if
1348  end if
1349  end if
1350 
1351  ! compute flux of magnetic field
1352  ! f_i[b_k]=v_i*b_k-v_k*b_i
1353  do idir=1,ndir
1354  if (idim==idir) then
1355  ! f_i[b_i] should be exactly 0, so we do not use the transport flux
1356  if (mhd_glm) then
1357  f(ixo^s,mag(idir))=w(ixo^s,psi_)
1358  else
1359  f(ixo^s,mag(idir))=zero
1360  end if
1361  else
1362  f(ixo^s,mag(idir))=w(ixo^s,mom(idim))*w(ixo^s,mag(idir))-w(ixo^s,mag(idim))*w(ixo^s,mom(idir))
1363 
1364  if (b0field) then
1365  f(ixo^s,mag(idir))=f(ixo^s,mag(idir))&
1366  +w(ixo^s,mom(idim))*block%B0(ixo^s,idir,idim)&
1367  -w(ixo^s,mom(idir))*block%B0(ixo^s,idim,idim)
1368  end if
1369 
1370  if (mhd_hall) then
1371  ! f_i[b_k] = f_i[b_k] + vHall_i*b_k - vHall_k*b_i
1372  if (mhd_etah>zero) then
1373  if (b0field) then
1374  f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
1375  - vhall(ixo^s,idir)*(w(ixo^s,mag(idim))+block%B0(ixo^s,idim,idim)) &
1376  + vhall(ixo^s,idim)*(w(ixo^s,mag(idir))+block%B0(ixo^s,idir,idim))
1377  else
1378  f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
1379  - vhall(ixo^s,idir)*w(ixo^s,mag(idim)) &
1380  + vhall(ixo^s,idim)*w(ixo^s,mag(idir))
1381  end if
1382  end if
1383  end if
1384 
1385  end if
1386  end do
1387 
1388  if (mhd_glm) then
1389  !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645
1390  f(ixo^s,psi_) = cmax_global**2*w(ixo^s,mag(idim))
1391  end if
1392 
1393  end subroutine mhd_get_flux
1394 
1395  !> w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
1396  subroutine mhd_add_source(qdt,ixI^L,ixO^L,wCT,w,x,qsourcesplit,active)
1400  use mod_gravity, only: gravity_add_source
1401 
1402  integer, intent(in) :: ixI^L, ixO^L
1403  double precision, intent(in) :: qdt
1404  double precision, intent(in) :: wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1405  double precision, intent(inout) :: w(ixi^s,1:nw)
1406  logical, intent(in) :: qsourcesplit
1407  logical, intent(inout) :: active
1408 
1409  if (.not. qsourcesplit) then
1410  ! Source for solving internal energy
1411  if(mhd_energy) then
1412  if(mhd_internal_e) then
1413  active = .true.
1414  call internal_energy_add_source(qdt,ixi^l,ixo^l,wct,w,x,e_)
1415  else if(mhd_solve_eaux) then
1416  active = .true.
1417  call internal_energy_add_source(qdt,ixi^l,ixo^l,wct,w,x,eaux_)
1418  end if
1419  endif
1420 
1421  ! Source for B0 splitting
1422  if (b0field) then
1423  active = .true.
1424  call add_source_b0split(qdt,ixi^l,ixo^l,wct,w,x)
1425  end if
1426 
1427  ! Sources for resistivity in eqs. for e, B1, B2 and B3
1428  if (abs(mhd_eta)>smalldouble)then
1429  active = .true.
1430  call add_source_res2(qdt,ixi^l,ixo^l,wct,w,x)
1431  end if
1432 
1433  if (mhd_eta_hyper>0.d0)then
1434  active = .true.
1435  call add_source_hyperres(qdt,ixi^l,ixo^l,wct,w,x)
1436  end if
1437  end if
1438 
1439  {^nooned
1440  if(.not.source_split_divb .and. .not.qsourcesplit .and. istep==nstep) then
1441  ! Sources related to div B
1442  select case (type_divb)
1443  case (divb_none)
1444  ! Do nothing
1445  case (divb_glm)
1446  active = .true.
1447  call add_source_glm(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1448  case (divb_powel)
1449  active = .true.
1450  call add_source_powel(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1451  case (divb_janhunen)
1452  active = .true.
1453  call add_source_janhunen(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1454  case (divb_linde)
1455  active = .true.
1456  call add_source_linde(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1457  case (divb_lindejanhunen)
1458  active = .true.
1459  call add_source_linde(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1460  call add_source_janhunen(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1461  case (divb_lindepowel)
1462  active = .true.
1463  call add_source_linde(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1464  call add_source_powel(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1465  case (divb_lindeglm)
1466  active = .true.
1467  call add_source_linde(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1468  call add_source_glm(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1469  case (divb_ct)
1470  continue ! Do nothing
1471  case (divb_multigrid)
1472  continue ! Do nothing
1473  case default
1474  call mpistop('Unknown divB fix')
1475  end select
1476  else if(source_split_divb .and. qsourcesplit) then
1477  ! Sources related to div B
1478  select case (type_divb)
1479  case (divb_none)
1480  ! Do nothing
1481  case (divb_glm)
1482  active = .true.
1483  call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
1484  case (divb_powel)
1485  active = .true.
1486  call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
1487  case (divb_janhunen)
1488  active = .true.
1489  call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
1490  case (divb_linde)
1491  active = .true.
1492  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
1493  case (divb_lindejanhunen)
1494  active = .true.
1495  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
1496  call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
1497  case (divb_lindepowel)
1498  active = .true.
1499  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
1500  call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
1501  case (divb_lindeglm)
1502  active = .true.
1503  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
1504  call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
1505  case (divb_ct)
1506  continue ! Do nothing
1507  case (divb_multigrid)
1508  continue ! Do nothing
1509  case default
1510  call mpistop('Unknown divB fix')
1511  end select
1512  end if
1513  }
1514 
1515  if(mhd_radiative_cooling) then
1516  call radiative_cooling_add_source(qdt,ixi^l,ixo^l,wct,&
1517  w,x,qsourcesplit,active)
1518  end if
1519 
1520  if(mhd_viscosity) then
1521  call viscosity_add_source(qdt,ixi^l,ixo^l,wct,&
1522  w,x,mhd_energy,qsourcesplit,active)
1523  end if
1524 
1525  if(mhd_gravity) then
1526  call gravity_add_source(qdt,ixi^l,ixo^l,wct,&
1527  w,x,total_energy,qsourcesplit,active)
1528  end if
1529 
1530  if (mhd_boris_type == boris_reduced_force) then
1531  call boris_add_source(qdt,ixi^l,ixo^l,wct,&
1532  w,x,qsourcesplit,active)
1533  end if
1534 
1535  end subroutine mhd_add_source
1536 
1537  subroutine boris_add_source(qdt,ixI^L,ixO^L,wCT,w,x,&
1538  qsourcesplit,active)
1540 
1541  integer, intent(in) :: ixI^L, ixO^L
1542  double precision, intent(in) :: qdt, x(ixi^s,1:ndim)
1543  double precision, intent(in) :: wCT(ixi^s,1:nw)
1544  double precision, intent(inout) :: w(ixi^s,1:nw)
1545  logical, intent(in) :: qsourcesplit
1546  logical, intent(inout) :: active
1547 
1548  double precision :: JxB(ixi^s,3)
1549  double precision :: gamma_A2(ixo^s)
1550  integer :: idir
1551 
1552  ! Boris source term is always unsplit
1553  if (qsourcesplit) return
1554 
1555  call get_lorentz(ixi^l,ixo^l,w,jxb)
1556  call mhd_gamma2_alfven(ixi^l, ixo^l, wct, gamma_a2)
1557 
1558  do idir = 1, ndir
1559  w(ixo^s,mom(idir)) = w(ixo^s,mom(idir)) + &
1560  qdt * gamma_a2 * jxb(ixo^s, idir)
1561  end do
1562 
1563  end subroutine boris_add_source
1564 
1565  !> Compute the Lorentz force (JxB)
1566  subroutine get_lorentz(ixI^L,ixO^L,w,JxB)
1568  integer, intent(in) :: ixI^L, ixO^L
1569  double precision, intent(in) :: w(ixi^s,1:nw)
1570  double precision, intent(inout) :: JxB(ixi^s,3)
1571  double precision :: a(ixi^s,3), b(ixi^s,3)
1572  integer :: idir, idirmin
1573  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
1574  double precision :: current(ixi^s,7-2*ndir:3)
1575 
1576  b=0.0d0
1577  do idir = 1, ndir
1578  b(ixo^s, idir) = mhd_mag_i_all(w, ixi^l, ixo^l,idir)
1579  end do
1580 
1581  ! store J current in a
1582  call get_current(w,ixi^l,ixo^l,idirmin,current)
1583 
1584  a=0.0d0
1585  do idir=7-2*ndir,3
1586  a(ixo^s,idir)=current(ixo^s,idir)
1587  end do
1588 
1589  call cross_product(ixi^l,ixo^l,a,b,jxb)
1590  end subroutine get_lorentz
1591 
1592  !> Compute 1/(1+v_A^2/c^2) for Boris' approximation, where v_A is the Alfven
1593  !> velocity
1594  subroutine mhd_gamma2_alfven(ixI^L, ixO^L, w, gamma_A2)
1596  integer, intent(in) :: ixI^L, ixO^L
1597  double precision, intent(in) :: w(ixi^s, nw)
1598  double precision, intent(out) :: gamma_A2(ixo^s)
1599 
1600  if (mhd_boris_c < 0.0d0) then
1601  ! Good for testing the non-conservative momentum treatment
1602  gamma_a2 = 1.0d0
1603  else
1604  ! Compute the inverse of 1 + B^2/(rho * c^2)
1605  gamma_a2 = 1.0d0 / (1.0d0 + mhd_mag_en_all(w, ixi^l, ixo^l) / &
1606  (w(ixo^s, rho_) * mhd_boris_c**2))
1607  end if
1608  end subroutine mhd_gamma2_alfven
1609 
1610  !> Compute 1/sqrt(1+v_A^2/c^2) for Boris simplification, where v_A is the
1611  !> Alfven velocity
1612  function mhd_gamma_alfven(w, ixI^L, ixO^L) result(gamma_A)
1614  integer, intent(in) :: ixI^L, ixO^L
1615  double precision, intent(in) :: w(ixi^s, nw)
1616  double precision :: gamma_A(ixo^s)
1617 
1618  call mhd_gamma2_alfven(ixi^l, ixo^l, w, gamma_a)
1619  gamma_a = sqrt(gamma_a)
1620  end function mhd_gamma_alfven
1621 
1622  subroutine internal_energy_add_source(qdt,ixI^L,ixO^L,wCT,w,x,ie)
1624  use mod_geometry
1625 
1626  integer, intent(in) :: ixI^L, ixO^L,ie
1627  double precision, intent(in) :: qdt
1628  double precision, intent(in) :: wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1629  double precision, intent(inout) :: w(ixi^s,1:nw)
1630  double precision :: pth(ixi^s),v(ixi^s,1:ndir),divv(ixi^s)
1631 
1632  call mhd_get_v(wct,x,ixi^l,ixi^l,v)
1633  if(slab_uniform) then
1634  if(nghostcells .gt. 2) then
1635  call divvector(v,ixi^l,ixo^l,divv,sixthorder=.true.)
1636  else
1637  call divvector(v,ixi^l,ixo^l,divv,fourthorder=.true.)
1638  end if
1639  else
1640  call divvector(v,ixi^l,ixo^l,divv)
1641  end if
1642  call mhd_get_pthermal(wct,x,ixi^l,ixo^l,pth)
1643  w(ixo^s,ie)=w(ixo^s,ie)-qdt*pth(ixo^s)*divv(ixo^s)
1644  end subroutine internal_energy_add_source
1645 
1646  !> Source terms after split off time-independent magnetic field
1647  subroutine add_source_b0split(qdt,ixI^L,ixO^L,wCT,w,x)
1649 
1650  integer, intent(in) :: ixI^L, ixO^L
1651  double precision, intent(in) :: qdt, wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1652  double precision, intent(inout) :: w(ixi^s,1:nw)
1653 
1654  double precision :: a(ixi^s,3), b(ixi^s,3), axb(ixi^s,3)
1655  integer :: idir
1656 
1657  a=0.d0
1658  b=0.d0
1659  ! for force-free field J0xB0 =0
1660  if(.not.b0field_forcefree) then
1661  ! store B0 magnetic field in b
1662  b(ixo^s,1:ndir)=block%B0(ixo^s,1:ndir,0)
1663 
1664  ! store J0 current in a
1665  do idir=7-2*ndir,3
1666  a(ixo^s,idir)=block%J0(ixo^s,idir)
1667  end do
1668  call cross_product(ixi^l,ixo^l,a,b,axb)
1669  axb(ixo^s,:)=axb(ixo^s,:)*qdt
1670  ! add J0xB0 source term in momentum equations
1671  w(ixo^s,mom(1:ndir))=w(ixo^s,mom(1:ndir))+axb(ixo^s,1:ndir)
1672  end if
1673 
1674  if(total_energy) then
1675  a=0.d0
1676  ! for free-free field -(vxB0) dot J0 =0
1677  b(ixo^s,:)=wct(ixo^s,mag(:))
1678  ! store full magnetic field B0+B1 in b
1679  if(.not.b0field_forcefree) b(ixo^s,:)=b(ixo^s,:)+block%B0(ixo^s,:,0)
1680  ! store velocity in a
1681  do idir=1,ndir
1682  a(ixo^s,idir)=wct(ixo^s,mom(idir))/wct(ixo^s,rho_)
1683  end do
1684  call cross_product(ixi^l,ixo^l,a,b,axb)
1685  axb(ixo^s,:)=axb(ixo^s,:)*qdt
1686  ! add -(vxB) dot J0 source term in energy equation
1687  do idir=7-2*ndir,3
1688  w(ixo^s,e_)=w(ixo^s,e_)-axb(ixo^s,idir)*block%J0(ixo^s,idir)
1689  end do
1690  end if
1691 
1692  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_B0')
1693 
1694  end subroutine add_source_b0split
1695 
1696  !> Add resistive source to w within ixO Uses 3 point stencil (1 neighbour) in
1697  !> each direction, non-conservative. If the fourthorder precompiler flag is
1698  !> set, uses fourth order central difference for the laplacian. Then the
1699  !> stencil is 5 (2 neighbours).
1700  subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
1702  use mod_usr_methods
1703  use mod_geometry
1704 
1705  integer, intent(in) :: ixI^L, ixO^L
1706  double precision, intent(in) :: qdt
1707  double precision, intent(in) :: wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1708  double precision, intent(inout) :: w(ixi^s,1:nw)
1709  integer :: ixA^L,idir,jdir,kdir,idirmin,idim,jxO^L,hxO^L,ix
1710  integer :: lxO^L, kxO^L
1711 
1712  double precision :: tmp(ixi^s),tmp2(ixi^s)
1713 
1714  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
1715  double precision :: current(ixi^s,7-2*ndir:3),eta(ixi^s)
1716  double precision :: gradeta(ixi^s,1:ndim), Bf(ixi^s,1:ndir)
1717 
1718  ! Calculating resistive sources involve one extra layer
1719  if (mhd_4th_order) then
1720  ixa^l=ixo^l^ladd2;
1721  else
1722  ixa^l=ixo^l^ladd1;
1723  end if
1724 
1725  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
1726  call mpistop("Error in add_source_res1: Non-conforming input limits")
1727 
1728  ! Calculate current density and idirmin
1729  call get_current(wct,ixi^l,ixo^l,idirmin,current)
1730 
1731  if (mhd_eta>zero)then
1732  eta(ixa^s)=mhd_eta
1733  gradeta(ixo^s,1:ndim)=zero
1734  else
1735  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
1736  ! assumes that eta is not function of current?
1737  do idim=1,ndim
1738  call gradient(eta,ixi^l,ixo^l,idim,tmp)
1739  gradeta(ixo^s,idim)=tmp(ixo^s)
1740  end do
1741  end if
1742 
1743  if(b0field) then
1744  bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))+block%B0(ixi^s,1:ndir,0)
1745  else
1746  bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
1747  end if
1748 
1749  do idir=1,ndir
1750  ! Put B_idir into tmp2 and eta*Laplace B_idir into tmp
1751  if (mhd_4th_order) then
1752  tmp(ixo^s)=zero
1753  tmp2(ixi^s)=bf(ixi^s,idir)
1754  do idim=1,ndim
1755  lxo^l=ixo^l+2*kr(idim,^d);
1756  jxo^l=ixo^l+kr(idim,^d);
1757  hxo^l=ixo^l-kr(idim,^d);
1758  kxo^l=ixo^l-2*kr(idim,^d);
1759  tmp(ixo^s)=tmp(ixo^s)+&
1760  (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
1761  /(12.0d0 * dxlevel(idim)**2)
1762  end do
1763  else
1764  tmp(ixo^s)=zero
1765  tmp2(ixi^s)=bf(ixi^s,idir)
1766  do idim=1,ndim
1767  jxo^l=ixo^l+kr(idim,^d);
1768  hxo^l=ixo^l-kr(idim,^d);
1769  tmp(ixo^s)=tmp(ixo^s)+&
1770  (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/dxlevel(idim)**2
1771  end do
1772  end if
1773 
1774  ! Multiply by eta
1775  tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
1776 
1777  ! Subtract grad(eta) x J = eps_ijk d_j eta J_k if eta is non-constant
1778  if (mhd_eta<zero)then
1779  do jdir=1,ndim; do kdir=idirmin,3
1780  if (lvc(idir,jdir,kdir)/=0)then
1781  if (lvc(idir,jdir,kdir)==1)then
1782  tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
1783  else
1784  tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
1785  end if
1786  end if
1787  end do; end do
1788  end if
1789 
1790  ! Add sources related to eta*laplB-grad(eta) x J to B and e
1791  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
1792  if (mhd_energy) then
1793  w(ixo^s,e_)=w(ixo^s,e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
1794  if(mhd_solve_eaux) then
1795  w(ixo^s,eaux_)=w(ixo^s,eaux_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
1796  end if
1797  end if
1798  end do ! idir
1799 
1800  if (mhd_energy) then
1801  ! de/dt+=eta*J**2
1802  tmp(ixo^s)=zero
1803  do idir=idirmin,3
1804  tmp(ixo^s)=tmp(ixo^s)+current(ixo^s,idir)**2
1805  end do
1806  w(ixo^s,e_)=w(ixo^s,e_)+qdt*eta(ixo^s)*tmp(ixo^s)
1807  if(mhd_solve_eaux) then
1808  w(ixo^s,eaux_)=w(ixo^s,eaux_)+qdt*eta(ixo^s)*tmp(ixo^s)
1809  end if
1810  end if
1811 
1812  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_res1')
1813 
1814  end subroutine add_source_res1
1815 
1816  !> Add resistive source to w within ixO
1817  !> Uses 5 point stencil (2 neighbours) in each direction, conservative
1818  subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
1820  use mod_usr_methods
1821  use mod_geometry
1822 
1823  integer, intent(in) :: ixI^L, ixO^L
1824  double precision, intent(in) :: qdt
1825  double precision, intent(in) :: wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1826  double precision, intent(inout) :: w(ixi^s,1:nw)
1827 
1828  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
1829  double precision :: current(ixi^s,7-2*ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
1830  double precision :: tmpvec(ixi^s,1:3)
1831  integer :: ixA^L,idir,idirmin,idirmin1
1832 
1833  ixa^l=ixo^l^ladd2;
1834 
1835  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
1836  call mpistop("Error in add_source_res2: Non-conforming input limits")
1837 
1838  ixa^l=ixo^l^ladd1;
1839  ! Calculate current density within ixL: J=curl B, thus J_i=eps_ijk*d_j B_k
1840  ! Determine exact value of idirmin while doing the loop.
1841  call get_current(wct,ixi^l,ixa^l,idirmin,current)
1842 
1843  if (mhd_eta>zero)then
1844  eta(ixa^s)=mhd_eta
1845  else
1846  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
1847  end if
1848 
1849  ! dB/dt= -curl(J*eta), thus B_i=B_i-eps_ijk d_j Jeta_k
1850  tmpvec(ixa^s,1:ndir)=zero
1851  do idir=idirmin,3
1852  tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
1853  end do
1854  curlj=0.d0
1855  call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
1856  if(stagger_grid.and.ndim==2.and.ndir==3) then
1857  ! if 2.5D
1858  w(ixo^s,mag(ndir)) = w(ixo^s,mag(ndir))-qdt*curlj(ixo^s,ndir)
1859  else
1860  w(ixo^s,mag(1:ndir)) = w(ixo^s,mag(1:ndir))-qdt*curlj(ixo^s,1:ndir)
1861  end if
1862 
1863  if(mhd_energy) then
1864  ! de/dt= +div(B x Jeta) = eta J^2 - B dot curl(eta J)
1865  ! de1/dt= eta J^2 - B1 dot curl(eta J)
1866  w(ixo^s,e_)=w(ixo^s,e_)+qdt*(sum(current(ixo^s,:)**2,dim=ndim+1)*eta(ixo^s)-&
1867  sum(wct(ixo^s,mag(1:ndir))*curlj(ixo^s,1:ndir),dim=ndim+1))
1868  if(mhd_solve_eaux) then
1869  w(ixo^s,eaux_)=w(ixo^s,eaux_)+qdt*(sum(current(ixo^s,:)**2,dim=ndim+1)*eta(ixo^s)-&
1870  sum(wct(ixo^s,mag(1:ndir))*curlj(ixo^s,1:ndir),dim=ndim+1))
1871  end if
1872  end if
1873 
1874  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_res2')
1875  end subroutine add_source_res2
1876 
1877  !> Add Hyper-resistive source to w within ixO
1878  !> Uses 9 point stencil (4 neighbours) in each direction.
1879  subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
1881  use mod_geometry
1882 
1883  integer, intent(in) :: ixI^L, ixO^L
1884  double precision, intent(in) :: qdt
1885  double precision, intent(in) :: wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1886  double precision, intent(inout) :: w(ixi^s,1:nw)
1887  !.. local ..
1888  double precision :: current(ixi^s,7-2*ndir:3)
1889  double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
1890  integer :: ixA^L,idir,jdir,kdir,idirmin,idirmin1
1891 
1892  ixa^l=ixo^l^ladd3;
1893  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
1894  call mpistop("Error in add_source_hyperres: Non-conforming input limits")
1895 
1896  call get_current(wct,ixi^l,ixa^l,idirmin,current)
1897  tmpvec(ixa^s,1:ndir)=zero
1898  do jdir=idirmin,3
1899  tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
1900  end do
1901 
1902  ixa^l=ixo^l^ladd2;
1903  call curlvector(tmpvec,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
1904 
1905  ixa^l=ixo^l^ladd1;
1906  tmpvec(ixa^s,1:ndir)=zero
1907  call curlvector(tmpvec2,ixi^l,ixa^l,tmpvec,idirmin1,1,3)
1908  ehyper(ixa^s,1:ndir) = - tmpvec(ixa^s,1:ndir)*mhd_eta_hyper
1909 
1910  ixa^l=ixo^l;
1911  tmpvec2(ixa^s,1:ndir)=zero
1912  call curlvector(ehyper,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
1913 
1914  do idir=1,ndir
1915  w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
1916  end do
1917 
1918  if (mhd_energy) then
1919  ! de/dt= +div(B x Ehyper)
1920  ixa^l=ixo^l^ladd1;
1921  tmpvec2(ixa^s,1:ndir)=zero
1922  do idir=1,ndir; do jdir=1,ndir; do kdir=idirmin,3
1923  tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
1924  + lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
1925  end do; end do; end do
1926  tmp(ixo^s)=zero
1927  call divvector(tmpvec2,ixi^l,ixo^l,tmp)
1928  w(ixo^s,e_)=w(ixo^s,e_)+tmp(ixo^s)*qdt
1929  if(mhd_solve_eaux) then
1930  w(ixo^s,eaux_)=w(ixo^s,eaux_)+tmp(ixo^s)*qdt
1931  end if
1932  end if
1933 
1934  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_hyperres')
1935 
1936  end subroutine add_source_hyperres
1937 
1938  subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
1939  ! Add divB related sources to w within ixO
1940  ! corresponding to Dedner JCP 2002, 175, 645 _equation 24_
1941  ! giving the EGLM-MHD scheme
1943  use mod_geometry
1944 
1945  integer, intent(in) :: ixI^L, ixO^L
1946  double precision, intent(in) :: qdt, wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1947  double precision, intent(inout) :: w(ixi^s,1:nw)
1948  double precision:: divb(ixi^s)
1949  integer :: idim,idir
1950  double precision :: gradPsi(ixi^s)
1951 
1952  ! We calculate now div B
1953  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
1954 
1955  ! dPsi/dt = - Ch^2/Cp^2 Psi
1956  if (mhd_glm_alpha < zero) then
1957  w(ixo^s,psi_) = abs(mhd_glm_alpha)*wct(ixo^s,psi_)
1958  else
1959  ! implicit update of Psi variable
1960  ! equation (27) in Mignone 2010 J. Com. Phys. 229, 2117
1961  if(slab_uniform) then
1962  w(ixo^s,psi_) = dexp(-qdt*cmax_global*mhd_glm_alpha/minval(dxlevel(:)))*w(ixo^s,psi_)
1963  else
1964  w(ixo^s,psi_) = dexp(-qdt*cmax_global*mhd_glm_alpha/minval(block%ds(ixo^s,:),dim=ndim+1))*w(ixo^s,psi_)
1965  end if
1966  end if
1967 
1968  ! gradient of Psi
1969  do idim=1,ndim
1970  select case(typegrad)
1971  case("central")
1972  call gradient(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
1973  case("limited")
1974  call gradients(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
1975  end select
1976  if (total_energy) then
1977  ! e = e -qdt (b . grad(Psi))
1978  w(ixo^s,e_) = w(ixo^s,e_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
1979  end if
1980  end do
1981 
1982  ! m = m - qdt b div b
1983  do idir=1,ndir
1984  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-qdt*mhd_mag_i_all(w,ixi^l,ixo^l,idir)*divb(ixo^s)
1985  end do
1986 
1987  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_glm')
1988 
1989  end subroutine add_source_glm
1990 
1991  !> Add divB related sources to w within ixO corresponding to Powel
1992  subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
1994 
1995  integer, intent(in) :: ixI^L, ixO^L
1996  double precision, intent(in) :: qdt, wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1997  double precision, intent(inout) :: w(ixi^s,1:nw)
1998  double precision :: divb(ixi^s),v(ixi^s,1:ndir)
1999  integer :: idir
2000 
2001  ! We calculate now div B
2002  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
2003 
2004  ! calculate velocity
2005  call mhd_get_v(wct,x,ixi^l,ixo^l,v)
2006 
2007  if (total_energy) then
2008  ! e = e - qdt (v . b) * div b
2009  w(ixo^s,e_)=w(ixo^s,e_)-&
2010  qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=ndim+1)*divb(ixo^s)
2011  end if
2012 
2013  ! b = b - qdt v * div b
2014  do idir=1,ndir
2015  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
2016  end do
2017 
2018  ! m = m - qdt b div b
2019  do idir=1,ndir
2020  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-qdt*mhd_mag_i_all(w,ixi^l,ixo^l,idir)*divb(ixo^s)
2021  end do
2022 
2023  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_powel')
2024 
2025  end subroutine add_source_powel
2026 
2027  subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
2028  ! Add divB related sources to w within ixO
2029  ! corresponding to Janhunen, just the term in the induction equation.
2031 
2032  integer, intent(in) :: ixI^L, ixO^L
2033  double precision, intent(in) :: qdt, wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
2034  double precision, intent(inout) :: w(ixi^s,1:nw)
2035  double precision :: divb(ixi^s)
2036  integer :: idir
2037 
2038  ! We calculate now div B
2039  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
2040 
2041  ! b = b - qdt v * div b
2042  do idir=1,ndir
2043  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*wct(ixo^s,mom(idir))/wct(ixo^s,rho_)*divb(ixo^s)
2044  end do
2045 
2046  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_janhunen')
2047 
2048  end subroutine add_source_janhunen
2049 
2050  subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
2051  ! Add Linde's divB related sources to wnew within ixO
2053  use mod_geometry
2054 
2055  integer, intent(in) :: ixI^L, ixO^L
2056  double precision, intent(in) :: qdt, wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
2057  double precision, intent(inout) :: w(ixi^s,1:nw)
2058  integer :: idim, idir, ixp^L, i^D, iside
2059  double precision :: divb(ixi^s),graddivb(ixi^s)
2060  logical, dimension(-1:1^D&) :: leveljump
2061 
2062  ! Calculate div B
2063  ixp^l=ixo^l^ladd1;
2064  call get_divb(wct,ixi^l,ixp^l,divb, mhd_divb_4thorder)
2065 
2066  ! for AMR stability, retreat one cell layer from the boarders of level jump
2067  {do i^db=-1,1\}
2068  if(i^d==0|.and.) cycle
2069  if(neighbor_type(i^d,saveigrid)==2 .or. neighbor_type(i^d,saveigrid)==4) then
2070  leveljump(i^d)=.true.
2071  else
2072  leveljump(i^d)=.false.
2073  end if
2074  {end do\}
2075 
2076  ixp^l=ixo^l;
2077  do idim=1,ndim
2078  select case(idim)
2079  {case(^d)
2080  do iside=1,2
2081  i^dd=kr(^dd,^d)*(2*iside-3);
2082  if (leveljump(i^dd)) then
2083  if (iside==1) then
2084  ixpmin^d=ixomin^d-i^d
2085  else
2086  ixpmax^d=ixomax^d-i^d
2087  end if
2088  end if
2089  end do
2090  \}
2091  end select
2092  end do
2093 
2094  ! Add Linde's diffusive terms
2095  do idim=1,ndim
2096  ! Calculate grad_idim(divb)
2097  select case(typegrad)
2098  case("central")
2099  call gradient(divb,ixi^l,ixp^l,idim,graddivb)
2100  case("limited")
2101  call gradients(divb,ixi^l,ixp^l,idim,graddivb)
2102  end select
2103 
2104  ! Multiply by Linde's eta*dt = divbdiff*(c_max*dx)*dt = divbdiff*dx**2
2105  if (slab_uniform) then
2106  graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
2107  else
2108  graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
2109  /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
2110  end if
2111 
2112  w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
2113 
2114  if (typedivbdiff=='all' .and. total_energy) then
2115  ! e += B_idim*eta*grad_idim(divb)
2116  w(ixp^s,e_)=w(ixp^s,e_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
2117  end if
2118  end do
2119 
2120  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_linde')
2121 
2122  end subroutine add_source_linde
2123 
2124  !> Calculate div B within ixO
2125  subroutine get_divb(w,ixI^L,ixO^L,divb, fourthorder)
2128  use mod_geometry
2129 
2130  integer, intent(in) :: ixI^L, ixO^L
2131  double precision, intent(in) :: w(ixi^s,1:nw)
2132  double precision, intent(inout) :: divb(ixi^s)
2133  logical, intent(in), optional :: fourthorder
2134 
2135  double precision :: bvec(ixi^s,1:ndir)
2136  double precision :: divb_corner(ixi^s), sign
2137  double precision :: aux_vol(ixi^s)
2138  integer :: ixC^L, idir, ic^D, ix^L
2139 
2140  if(stagger_grid) then
2141  divb=0.d0
2142  do idir=1,ndim
2143  ixc^l=ixo^l-kr(idir,^d);
2144  divb(ixo^s)=divb(ixo^s)+block%ws(ixo^s,idir)*block%surfaceC(ixo^s,idir)-&
2145  block%ws(ixc^s,idir)*block%surfaceC(ixc^s,idir)
2146  end do
2147  divb(ixo^s)=divb(ixo^s)/block%dvolume(ixo^s)
2148  else
2149  bvec(ixi^s,:)=w(ixi^s,mag(:))
2150  select case(typediv)
2151  case("central")
2152  call divvector(bvec,ixi^l,ixo^l,divb,fourthorder)
2153  case("limited")
2154  call divvectors(bvec,ixi^l,ixo^l,divb)
2155  end select
2156  end if
2157 
2158  end subroutine get_divb
2159 
2160  !> get dimensionless div B = |divB| * volume / area / |B|
2161  subroutine get_normalized_divb(w,ixI^L,ixO^L,divb)
2164 
2165  integer, intent(in) :: ixI^L, ixO^L
2166  double precision, intent(in) :: w(ixi^s,1:nw)
2167  double precision :: divb(ixi^s), dsurface(ixi^s)
2168 
2169  integer :: ixA^L,idims
2170 
2171  call get_divb(w,ixi^l,ixo^l,divb)
2172  if(slab_uniform) then
2173  divb(ixo^s)=0.5d0*abs(divb(ixo^s))/sqrt(mhd_mag_en_all(w,ixi^l,ixo^l))/sum(1.d0/dxlevel(:))
2174  else
2175  ixamin^d=ixomin^d-1;
2176  ixamax^d=ixomax^d-1;
2177  dsurface(ixo^s)= sum(block%surfaceC(ixo^s,:),dim=ndim+1)
2178  do idims=1,ndim
2179  ixa^l=ixo^l-kr(idims,^d);
2180  dsurface(ixo^s)=dsurface(ixo^s)+block%surfaceC(ixa^s,idims)
2181  end do
2182  divb(ixo^s)=abs(divb(ixo^s))/sqrt(mhd_mag_en_all(w,ixi^l,ixo^l))*&
2183  block%dvolume(ixo^s)/dsurface(ixo^s)
2184  end if
2185 
2186  end subroutine get_normalized_divb
2187 
2188  !> Calculate idirmin and the idirmin:3 components of the common current array
2189  !> make sure that dxlevel(^D) is set correctly.
2190  subroutine get_current(w,ixI^L,ixO^L,idirmin,current)
2192  use mod_geometry
2193 
2194  integer, intent(in) :: ixO^L, ixI^L
2195  double precision, intent(in) :: w(ixi^s,1:nw)
2196  integer, intent(out) :: idirmin
2197  integer :: idir, idirmin0
2198 
2199  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
2200  double precision :: current(ixi^s,7-2*ndir:3),bvec(ixi^s,1:ndir)
2201 
2202  idirmin0 = 7-2*ndir
2203 
2204  bvec(ixi^s,1:ndir)=w(ixi^s,mag(1:ndir))
2205 
2206  call curlvector(bvec,ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
2207 
2208  if(b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
2209  block%J0(ixo^s,idirmin0:3)
2210 
2211  end subroutine get_current
2212 
2213  !> If resistivity is not zero, check diffusion time limit for dt
2214  subroutine mhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
2216  use mod_usr_methods
2218  use mod_viscosity, only: viscosity_get_dt
2219  use mod_gravity, only: gravity_get_dt
2220 
2221  integer, intent(in) :: ixI^L, ixO^L
2222  double precision, intent(inout) :: dtnew
2223  double precision, intent(in) :: dx^D
2224  double precision, intent(in) :: w(ixi^s,1:nw)
2225  double precision, intent(in) :: x(ixi^s,1:ndim)
2226 
2227  integer :: idirmin,idim
2228  double precision :: dxarr(ndim)
2229  double precision :: current(ixi^s,7-2*ndir:3),eta(ixi^s)
2230 
2231  dtnew = bigdouble
2232 
2233  ^d&dxarr(^d)=dx^d;
2234  if (mhd_eta>zero)then
2235  dtnew=dtdiffpar*minval(dxarr(1:ndim))**2/mhd_eta
2236  else if (mhd_eta<zero)then
2237  call get_current(w,ixi^l,ixo^l,idirmin,current)
2238  call usr_special_resistivity(w,ixi^l,ixo^l,idirmin,x,current,eta)
2239  dtnew=bigdouble
2240  do idim=1,ndim
2241  if(slab_uniform) then
2242  dtnew=min(dtnew,&
2243  dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
2244  else
2245  dtnew=min(dtnew,&
2246  dtdiffpar/(smalldouble+maxval(eta(ixo^s)/block%ds(ixo^s,idim)**2)))
2247  end if
2248  end do
2249  end if
2250 
2251  if(mhd_eta_hyper>zero) then
2252  if(slab_uniform) then
2253  dtnew=min(dtdiffpar*minval(dxarr(1:ndim))**4/mhd_eta_hyper,dtnew)
2254  else
2255  dtnew=min(dtdiffpar*minval(block%ds(ixo^s,1:ndim))**4/mhd_eta_hyper,dtnew)
2256  end if
2257  end if
2258 
2259  if(mhd_radiative_cooling) then
2260  call cooling_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
2261  end if
2262 
2263  if(mhd_viscosity) then
2264  call viscosity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
2265  end if
2266 
2267  if(mhd_gravity) then
2268  call gravity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
2269  end if
2270 
2271  end subroutine mhd_get_dt
2272 
2273  ! Add geometrical source terms to w
2274  subroutine mhd_add_source_geom(qdt,ixI^L,ixO^L,wCT,w,x)
2276  use mod_geometry
2277 
2278  integer, intent(in) :: ixI^L, ixO^L
2279  double precision, intent(in) :: qdt, x(ixi^s,1:ndim)
2280  double precision, intent(inout) :: wCT(ixi^s,1:nw), w(ixi^s,1:nw)
2281 
2282  integer :: iw,idir, h1x^L{^NOONED, h2x^L}
2283  double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s)
2284 
2285  integer :: mr_,mphi_ ! Polar var. names
2286  integer :: br_,bphi_
2287 
2288  mr_=mom(1); mphi_=mom(1)-1+phi_ ! Polar var. names
2289  br_=mag(1); bphi_=mag(1)-1+phi_
2290 
2291  select case (coordinate)
2292  case (cylindrical)
2293  if (angmomfix) then
2294  call mpistop("angmomfix not implemented yet in MHD")
2295  endif
2296  call mhd_get_p_total(wct,x,ixi^l,ixo^l,tmp)
2297  if(phi_>0) then
2298  w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*(tmp(ixo^s)-&
2299  wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2/wct(ixo^s,rho_))
2300  w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt/x(ixo^s,1)*(&
2301  -wct(ixo^s,mphi_)*wct(ixo^s,mr_)/wct(ixo^s,rho_) &
2302  +wct(ixo^s,bphi_)*wct(ixo^s,br_))
2303  if(.not.stagger_grid) then
2304  w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt/x(ixo^s,1)*&
2305  (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
2306  -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
2307  /wct(ixo^s,rho_)
2308  end if
2309  else
2310  w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*tmp(ixo^s)
2311  end if
2312  if(mhd_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,psi_)/x(ixo^s,1)
2313  case (spherical)
2314  h1x^l=ixo^l-kr(1,^d); {^nooned h2x^l=ixo^l-kr(2,^d);}
2315  call mhd_get_p_total(wct,x,ixi^l,ixo^l,tmp1)
2316  tmp(ixo^s)=tmp1(ixo^s)
2317  if(b0field) then
2318  tmp2(ixo^s)=sum(block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=ndim+1)
2319  tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
2320  end if
2321  ! m1
2322  tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
2323  *(block%surfaceC(ixo^s,1)-block%surfaceC(h1x^s,1))/block%dvolume(ixo^s)
2324  if(ndir>1) then
2325  do idir=2,ndir
2326  tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,mom(idir))**2/wct(ixo^s,rho_)-wct(ixo^s,mag(idir))**2
2327  if(b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
2328  end do
2329  end if
2330  w(ixo^s,mom(1))=w(ixo^s,mom(1))+qdt*tmp(ixo^s)/x(ixo^s,1)
2331  ! b1
2332  if(mhd_glm) then
2333  w(ixo^s,mag(1))=w(ixo^s,mag(1))+qdt/x(ixo^s,1)*2.0d0*wct(ixo^s,psi_)
2334  end if
2335 
2336  {^nooned
2337  ! m2
2338  tmp(ixo^s)=tmp1(ixo^s)
2339  if(b0field) then
2340  tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
2341  end if
2342  ! This will make hydrostatic p=const an exact solution
2343  w(ixo^s,mom(2))=w(ixo^s,mom(2))+qdt*tmp(ixo^s) &
2344  *(block%surfaceC(ixo^s,2)-block%surfaceC(h2x^s,2)) &
2345  /block%dvolume(ixo^s)
2346  tmp(ixo^s)=-(wct(ixo^s,mom(1))*wct(ixo^s,mom(2))/wct(ixo^s,rho_) &
2347  -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
2348  if (b0field) then
2349  tmp(ixo^s)=tmp(ixo^s)+block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
2350  +wct(ixo^s,mag(1))*block%B0(ixo^s,2,0)
2351  end if
2352  if(ndir==3) then
2353  tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom(3))**2/wct(ixo^s,rho_) &
2354  -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
2355  if (b0field) then
2356  tmp(ixo^s)=tmp(ixo^s)-2.0d0*block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
2357  *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
2358  end if
2359  end if
2360  w(ixo^s,mom(2))=w(ixo^s,mom(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
2361  ! b2
2362  if(.not.stagger_grid) then
2363  tmp(ixo^s)=(wct(ixo^s,mom(1))*wct(ixo^s,mag(2)) &
2364  -wct(ixo^s,mom(2))*wct(ixo^s,mag(1)))/wct(ixo^s,rho_)
2365  if(b0field) then
2366  tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom(1))*block%B0(ixo^s,2,0) &
2367  -wct(ixo^s,mom(2))*block%B0(ixo^s,1,0))/wct(ixo^s,rho_)
2368  end if
2369  if(mhd_glm) then
2370  tmp(ixo^s)=tmp(ixo^s) &
2371  + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,psi_)
2372  end if
2373  w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
2374  end if
2375  }
2376 
2377  if(ndir==3) then
2378  ! m3
2379  if(.not.angmomfix) then
2380  tmp(ixo^s)=-(wct(ixo^s,mom(3))*wct(ixo^s,mom(1))/wct(ixo^s,rho_) &
2381  -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
2382  -(wct(ixo^s,mom(2))*wct(ixo^s,mom(3))/wct(ixo^s,rho_) &
2383  -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
2384  *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
2385  if (b0field) then
2386  tmp(ixo^s)=tmp(ixo^s)+block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
2387  +wct(ixo^s,mag(1))*block%B0(ixo^s,3,0) {^nooned &
2388  +(block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
2389  +wct(ixo^s,mag(2))*block%B0(ixo^s,3,0)) &
2390  *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
2391  end if
2392  w(ixo^s,mom(3))=w(ixo^s,mom(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
2393  else
2394  call mpistop("angmomfix not implemented yet in MHD")
2395  end if
2396  ! b3
2397  if(.not.stagger_grid) then
2398  tmp(ixo^s)=(wct(ixo^s,mom(1))*wct(ixo^s,mag(3)) &
2399  -wct(ixo^s,mom(3))*wct(ixo^s,mag(1)))/wct(ixo^s,rho_) {^nooned &
2400  -(wct(ixo^s,mom(3))*wct(ixo^s,mag(2)) &
2401  -wct(ixo^s,mom(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
2402  /(wct(ixo^s,rho_)*dsin(x(ixo^s,2))) }
2403  if (b0field) then
2404  tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom(1))*block%B0(ixo^s,3,0) &
2405  -wct(ixo^s,mom(3))*block%B0(ixo^s,1,0))/wct(ixo^s,rho_){^nooned &
2406  -(wct(ixo^s,mom(3))*block%B0(ixo^s,2,0) &
2407  -wct(ixo^s,mom(2))*block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
2408  /(wct(ixo^s,rho_)*dsin(x(ixo^s,2))) }
2409  end if
2410  w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
2411  end if
2412  end if
2413  end select
2414  end subroutine mhd_add_source_geom
2415 
2416  !> Compute 2 times total magnetic energy
2417  function mhd_mag_en_all(w, ixI^L, ixO^L) result(mge)
2419  integer, intent(in) :: ixI^L, ixO^L
2420  double precision, intent(in) :: w(ixi^s, nw)
2421  double precision :: mge(ixo^s)
2422 
2423  if (b0field) then
2424  mge = sum((w(ixo^s, mag(:))+block%B0(ixo^s,:,block%iw0))**2, dim=ndim+1)
2425  else
2426  mge = sum(w(ixo^s, mag(:))**2, dim=ndim+1)
2427  end if
2428  end function mhd_mag_en_all
2429 
2430  !> Compute full magnetic field by direction
2431  function mhd_mag_i_all(w, ixI^L, ixO^L,idir) result(mgf)
2433  integer, intent(in) :: ixI^L, ixO^L, idir
2434  double precision, intent(in) :: w(ixi^s, nw)
2435  double precision :: mgf(ixo^s)
2436 
2437  if (b0field) then
2438  mgf = w(ixo^s, mag(idir))+block%B0(ixo^s,idir,block%iw0)
2439  else
2440  mgf = w(ixo^s, mag(idir))
2441  end if
2442  end function mhd_mag_i_all
2443 
2444  !> Compute evolving magnetic energy
2445  function mhd_mag_en(w, ixI^L, ixO^L) result(mge)
2447  integer, intent(in) :: ixI^L, ixO^L
2448  double precision, intent(in) :: w(ixi^s, nw)
2449  double precision :: mge(ixo^s)
2450 
2451  mge = 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
2452  end function mhd_mag_en
2453 
2454  !> compute kinetic energy
2455  function mhd_kin_en(w, ixI^L, ixO^L, inv_rho) result(ke)
2457  integer, intent(in) :: ixI^L, ixO^L
2458  double precision, intent(in) :: w(ixi^s, nw)
2459  double precision :: ke(ixo^s)
2460  double precision, intent(in), optional :: inv_rho(ixo^s)
2461 
2462  if (present(inv_rho)) then
2463  ke = 0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) * inv_rho
2464  else
2465  ke = 0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) / w(ixo^s, rho_)
2466  end if
2467  end function mhd_kin_en
2468 
2469  subroutine mhd_getv_hall(w,x,ixI^L,ixO^L,vHall)
2471 
2472  integer, intent(in) :: ixI^L, ixO^L
2473  double precision, intent(in) :: w(ixi^s,nw)
2474  double precision, intent(in) :: x(ixi^s,1:ndim)
2475  double precision, intent(inout) :: vHall(ixi^s,1:3)
2476 
2477  integer :: idir, idirmin
2478  double precision :: current(ixi^s,7-2*ndir:3)
2479 
2480  ! Calculate current density and idirmin
2481  call get_current(w,ixi^l,ixo^l,idirmin,current)
2482  vhall(ixo^s,1:3) = zero
2483  vhall(ixo^s,idirmin:3) = - mhd_etah*current(ixo^s,idirmin:3)
2484  do idir = idirmin, 3
2485  vhall(ixo^s,idir) = vhall(ixo^s,idir)/w(ixo^s,rho_)
2486  end do
2487 
2488  end subroutine mhd_getv_hall
2489 
2490  subroutine mhd_getdt_hall(w,x,ixI^L,ixO^L,dx^D,dthall)
2492 
2493  integer, intent(in) :: ixI^L, ixO^L
2494  double precision, intent(in) :: dx^D
2495  double precision, intent(in) :: w(ixi^s,1:nw)
2496  double precision, intent(in) :: x(ixi^s,1:ndim)
2497  double precision, intent(out) :: dthall
2498  !.. local ..
2499  double precision :: dxarr(ndim)
2500  double precision :: bmag(ixi^s)
2501 
2502  dthall=bigdouble
2503 
2504  ! because we have that in cmax now:
2505  return
2506 
2507  ^d&dxarr(^d)=dx^d;
2508 
2509  if (.not. b0field) then
2510  bmag(ixo^s)=sqrt(sum(w(ixo^s,mag(:))**2, dim=ndim+1))
2511  bmag(ixo^s)=sqrt(sum((w(ixo^s,mag(:)) + block%B0(ixo^s,1:ndir,block%iw0))**2))
2512  end if
2513 
2514  if(slab_uniform) then
2515  dthall=dtdiffpar*minval(dxarr(1:ndim))**2.0d0/(mhd_etah*maxval(bmag(ixo^s)/w(ixo^s,rho_)))
2516  else
2517  dthall=dtdiffpar*minval(block%ds(ixo^s,1:ndim))**2.0d0/(mhd_etah*maxval(bmag(ixo^s)/w(ixo^s,rho_)))
2518  end if
2519 
2520  end subroutine mhd_getdt_hall
2521 
2522  subroutine mhd_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
2524  use mod_usr_methods
2525  integer, intent(in) :: ixI^L, ixO^L, idir
2526  double precision, intent(in) :: qt
2527  double precision, intent(inout) :: wLC(ixi^s,1:nw), wRC(ixi^s,1:nw)
2528  double precision, intent(inout) :: wLp(ixi^s,1:nw), wRp(ixi^s,1:nw)
2529  type(state) :: s
2530  double precision :: dB(ixi^s), dPsi(ixi^s)
2531 
2532  if(stagger_grid) then
2533  wlc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
2534  wrc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
2535  wlp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
2536  wrp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
2537  else
2538  ! Solve the Riemann problem for the linear 2x2 system for normal
2539  ! B-field and GLM_Psi according to Dedner 2002:
2540  ! This implements eq. (42) in Dedner et al. 2002 JcP 175
2541  ! Gives the Riemann solution on the interface
2542  ! for the normal B component and Psi in the GLM-MHD system.
2543  ! 23/04/2013 Oliver Porth
2544  db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
2545  dpsi(ixo^s) = wrp(ixo^s,psi_) - wlp(ixo^s,psi_)
2546 
2547  wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
2548  - 0.5d0/cmax_global * dpsi(ixo^s)
2549  wlp(ixo^s,psi_) = 0.5d0 * (wrp(ixo^s,psi_) + wlp(ixo^s,psi_)) &
2550  - 0.5d0*cmax_global * db(ixo^s)
2551 
2552  wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
2553  wrp(ixo^s,psi_) = wlp(ixo^s,psi_)
2554  wrc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
2555  wrc(ixo^s,psi_) = wlp(ixo^s,psi_)
2556  wlc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
2557  wlc(ixo^s,psi_) = wlp(ixo^s,psi_)
2558  end if
2559 
2560  if(associated(usr_set_wlr)) call usr_set_wlr(ixi^l,ixo^l,qt,wlc,wrc,wlp,wrp,s,idir)
2561 
2562  end subroutine mhd_modify_wlr
2563 
2564  subroutine mhd_boundary_adjust
2566  integer :: iB, idim, iside, iigrid, igrid
2567  integer :: ixG^L, ixO^L, i^D
2568 
2569  ixg^l=ixg^ll;
2570  do iigrid=1,igridstail; igrid=igrids(iigrid);
2571  if(.not.phyboundblock(igrid)) cycle
2572  saveigrid=igrid
2573  block=>ps(igrid)
2574  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
2575  do idim=1,ndim
2576  ! to avoid using as yet unknown corner info in more than 1D, we
2577  ! fill only interior mesh ranges of the ghost cell ranges at first,
2578  ! and progressively enlarge the ranges to include corners later
2579  do iside=1,2
2580  i^d=kr(^d,idim)*(2*iside-3);
2581  if (neighbor_type(i^d,igrid)/=1) cycle
2582  ib=(idim-1)*2+iside
2583  if(.not.boundary_divbfix(ib)) cycle
2584  if(any(typeboundary(:,ib)=="special")) then
2585  ! MF nonlinear force-free B field extrapolation and data driven
2586  ! require normal B of the first ghost cell layer to be untouched by
2587  ! fixdivB=0 process, set boundary_divbfix_skip(iB)=1 in par file
2588  select case (idim)
2589  {case (^d)
2590  if (iside==2) then
2591  ! maximal boundary
2592  ixomin^dd=ixgmax^d+1-nghostcells+boundary_divbfix_skip(2*^d)^d%ixOmin^dd=ixgmin^dd;
2593  ixomax^dd=ixgmax^dd;
2594  else
2595  ! minimal boundary
2596  ixomin^dd=ixgmin^dd;
2597  ixomax^dd=ixgmin^d-1+nghostcells-boundary_divbfix_skip(2*^d-1)^d%ixOmax^dd=ixgmax^dd;
2598  end if \}
2599  end select
2600  call fixdivb_boundary(ixg^l,ixo^l,ps(igrid)%w,ps(igrid)%x,ib)
2601  end if
2602  end do
2603  end do
2604  end do
2605 
2606  end subroutine mhd_boundary_adjust
2607 
2608  subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
2610 
2611  integer, intent(in) :: ixG^L,ixO^L,iB
2612  double precision, intent(inout) :: w(ixg^s,1:nw)
2613  double precision, intent(in) :: x(ixg^s,1:ndim)
2614 
2615  double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
2616  integer :: ix^D,ixF^L
2617 
2618  select case(ib)
2619  case(1)
2620  ! 2nd order CD for divB=0 to set normal B component better
2621  if(total_energy) call mhd_to_primitive(ixg^l,ixo^l,w,x)
2622  {^iftwod
2623  ixfmin1=ixomin1+1
2624  ixfmax1=ixomax1+1
2625  ixfmin2=ixomin2+1
2626  ixfmax2=ixomax2-1
2627  if(slab_uniform) then
2628  dx1x2=dxlevel(1)/dxlevel(2)
2629  do ix1=ixfmax1,ixfmin1,-1
2630  w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
2631  +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
2632  w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
2633  enddo
2634  else
2635  do ix1=ixfmax1,ixfmin1,-1
2636  w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
2637  w(ix1,ixfmin2:ixfmax2,mag(1)))*block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
2638  +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
2639  block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
2640  -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
2641  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
2642  /block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
2643  end do
2644  end if
2645  }
2646  {^ifthreed
2647  ixfmin1=ixomin1+1
2648  ixfmax1=ixomax1+1
2649  ixfmin2=ixomin2+1
2650  ixfmax2=ixomax2-1
2651  ixfmin3=ixomin3+1
2652  ixfmax3=ixomax3-1
2653  if(slab_uniform) then
2654  dx1x2=dxlevel(1)/dxlevel(2)
2655  dx1x3=dxlevel(1)/dxlevel(3)
2656  do ix1=ixfmax1,ixfmin1,-1
2657  w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
2658  w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
2659  +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
2660  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
2661  +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
2662  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
2663  end do
2664  else
2665  do ix1=ixfmax1,ixfmin1,-1
2666  w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
2667  ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
2668  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
2669  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
2670  +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
2671  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
2672  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
2673  -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
2674  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
2675  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
2676  +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
2677  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
2678  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
2679  -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
2680  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
2681  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
2682  /block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
2683  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
2684  end do
2685  end if
2686  }
2687  if(total_energy) call mhd_to_conserved(ixg^l,ixo^l,w,x)
2688  case(2)
2689  if(total_energy) call mhd_to_primitive(ixg^l,ixo^l,w,x)
2690  {^iftwod
2691  ixfmin1=ixomin1-1
2692  ixfmax1=ixomax1-1
2693  ixfmin2=ixomin2+1
2694  ixfmax2=ixomax2-1
2695  if(slab_uniform) then
2696  dx1x2=dxlevel(1)/dxlevel(2)
2697  do ix1=ixfmin1,ixfmax1
2698  w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
2699  -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
2700  w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
2701  enddo
2702  else
2703  do ix1=ixfmin1,ixfmax1
2704  w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
2705  w(ix1,ixfmin2:ixfmax2,mag(1)))*block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
2706  -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
2707  block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
2708  +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
2709  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
2710  /block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
2711  end do
2712  end if
2713  }
2714  {^ifthreed
2715  ixfmin1=ixomin1-1
2716  ixfmax1=ixomax1-1
2717  ixfmin2=ixomin2+1
2718  ixfmax2=ixomax2-1
2719  ixfmin3=ixomin3+1
2720  ixfmax3=ixomax3-1
2721  if(slab_uniform) then
2722  dx1x2=dxlevel(1)/dxlevel(2)
2723  dx1x3=dxlevel(1)/dxlevel(3)
2724  do ix1=ixfmin1,ixfmax1
2725  w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
2726  w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
2727  -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
2728  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
2729  -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
2730  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
2731  end do
2732  else
2733  do ix1=ixfmin1,ixfmax1
2734  w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
2735  ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
2736  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
2737  block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
2738  -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
2739  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
2740  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
2741  +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
2742  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
2743  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
2744  -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
2745  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
2746  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
2747  +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
2748  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
2749  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
2750  /block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
2751  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
2752  end do
2753  end if
2754  }
2755  if(total_energy) call mhd_to_conserved(ixg^l,ixo^l,w,x)
2756  case(3)
2757  if(total_energy) call mhd_to_primitive(ixg^l,ixo^l,w,x)
2758  {^iftwod
2759  ixfmin1=ixomin1+1
2760  ixfmax1=ixomax1-1
2761  ixfmin2=ixomin2+1
2762  ixfmax2=ixomax2+1
2763  if(slab_uniform) then
2764  dx2x1=dxlevel(2)/dxlevel(1)
2765  do ix2=ixfmax2,ixfmin2,-1
2766  w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
2767  +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
2768  w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
2769  enddo
2770  else
2771  do ix2=ixfmax2,ixfmin2,-1
2772  w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
2773  w(ixfmin1:ixfmax1,ix2,mag(2)))*block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
2774  +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
2775  block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
2776  -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
2777  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
2778  /block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
2779  end do
2780  end if
2781  }
2782  {^ifthreed
2783  ixfmin1=ixomin1+1
2784  ixfmax1=ixomax1-1
2785  ixfmin3=ixomin3+1
2786  ixfmax3=ixomax3-1
2787  ixfmin2=ixomin2+1
2788  ixfmax2=ixomax2+1
2789  if(slab_uniform) then
2790  dx2x1=dxlevel(2)/dxlevel(1)
2791  dx2x3=dxlevel(2)/dxlevel(3)
2792  do ix2=ixfmax2,ixfmin2,-1
2793  w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
2794  ix2+1,ixfmin3:ixfmax3,mag(2)) &
2795  +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
2796  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
2797  +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
2798  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
2799  end do
2800  else
2801  do ix2=ixfmax2,ixfmin2,-1
2802  w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
2803  ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
2804  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
2805  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
2806  +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
2807  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
2808  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
2809  -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
2810  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
2811  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
2812  +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
2813  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
2814  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
2815  -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
2816  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
2817  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
2818  /block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
2819  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
2820  end do
2821  end if
2822  }
2823  if(total_energy) call mhd_to_conserved(ixg^l,ixo^l,w,x)
2824  case(4)
2825  if(total_energy) call mhd_to_primitive(ixg^l,ixo^l,w,x)
2826  {^iftwod
2827  ixfmin1=ixomin1+1
2828  ixfmax1=ixomax1-1
2829  ixfmin2=ixomin2-1
2830  ixfmax2=ixomax2-1
2831  if(slab_uniform) then
2832  dx2x1=dxlevel(2)/dxlevel(1)
2833  do ix2=ixfmin2,ixfmax2
2834  w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
2835  -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
2836  w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
2837  end do
2838  else
2839  do ix2=ixfmin2,ixfmax2
2840  w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
2841  w(ixfmin1:ixfmax1,ix2,mag(2)))*block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
2842  -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
2843  block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
2844  +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
2845  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
2846  /block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
2847  end do
2848  end if
2849  }
2850  {^ifthreed
2851  ixfmin1=ixomin1+1
2852  ixfmax1=ixomax1-1
2853  ixfmin3=ixomin3+1
2854  ixfmax3=ixomax3-1
2855  ixfmin2=ixomin2-1
2856  ixfmax2=ixomax2-1
2857  if(slab_uniform) then
2858  dx2x1=dxlevel(2)/dxlevel(1)
2859  dx2x3=dxlevel(2)/dxlevel(3)
2860  do ix2=ixfmin2,ixfmax2
2861  w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
2862  ix2-1,ixfmin3:ixfmax3,mag(2)) &
2863  -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
2864  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
2865  -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
2866  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
2867  end do
2868  else
2869  do ix2=ixfmin2,ixfmax2
2870  w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
2871  ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
2872  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
2873  block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
2874  -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
2875  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
2876  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
2877  +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
2878  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
2879  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
2880  -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
2881  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
2882  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
2883  +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
2884  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
2885  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
2886  /block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
2887  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
2888  end do
2889  end if
2890  }
2891  if(total_energy) call mhd_to_conserved(ixg^l,ixo^l,w,x)
2892  {^ifthreed
2893  case(5)
2894  if(total_energy) call mhd_to_primitive(ixg^l,ixo^l,w,x)
2895  ixfmin1=ixomin1+1
2896  ixfmax1=ixomax1-1
2897  ixfmin2=ixomin2+1
2898  ixfmax2=ixomax2-1
2899  ixfmin3=ixomin3+1
2900  ixfmax3=ixomax3+1
2901  if(slab_uniform) then
2902  dx3x1=dxlevel(3)/dxlevel(1)
2903  dx3x2=dxlevel(3)/dxlevel(2)
2904  do ix3=ixfmax3,ixfmin3,-1
2905  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
2906  ixfmin2:ixfmax2,ix3+1,mag(3)) &
2907  +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
2908  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
2909  +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
2910  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
2911  end do
2912  else
2913  do ix3=ixfmax3,ixfmin3,-1
2914  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
2915  ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
2916  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
2917  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
2918  +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
2919  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
2920  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
2921  -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
2922  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
2923  block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
2924  +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
2925  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
2926  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
2927  -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
2928  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
2929  block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
2930  /block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
2931  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
2932  end do
2933  end if
2934  if(total_energy) call mhd_to_conserved(ixg^l,ixo^l,w,x)
2935  case(6)
2936  if(total_energy) call mhd_to_primitive(ixg^l,ixo^l,w,x)
2937  ixfmin1=ixomin1+1
2938  ixfmax1=ixomax1-1
2939  ixfmin2=ixomin2+1
2940  ixfmax2=ixomax2-1
2941  ixfmin3=ixomin3-1
2942  ixfmax3=ixomax3-1
2943  if(slab_uniform) then
2944  dx3x1=dxlevel(3)/dxlevel(1)
2945  dx3x2=dxlevel(3)/dxlevel(2)
2946  do ix3=ixfmin3,ixfmax3
2947  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
2948  ixfmin2:ixfmax2,ix3-1,mag(3)) &
2949  -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
2950  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
2951  -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
2952  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
2953  end do
2954  else
2955  do ix3=ixfmin3,ixfmax3
2956  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
2957  ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
2958  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
2959  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
2960  -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
2961  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
2962  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
2963  +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
2964  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
2965  block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
2966  -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
2967  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
2968  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
2969  +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
2970  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
2971  block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
2972  /block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
2973  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
2974  end do
2975  end if
2976  if(total_energy) call mhd_to_conserved(ixg^l,ixo^l,w,x)
2977  }
2978  case default
2979  call mpistop("Special boundary is not defined for this region")
2980  end select
2981 
2982  end subroutine fixdivb_boundary
2983 
2984  {^nooned
2985  subroutine mhd_clean_divb_multigrid(qdt, qt, active)
2989  use mod_geometry
2990 
2991  double precision, intent(in) :: qdt !< Current time step
2992  double precision, intent(in) :: qt !< Current time
2993  logical, intent(inout) :: active !< Output if the source is active
2994  integer :: iigrid, igrid, id
2995  integer :: n, nc, lvl, ix^L, ixC^L, idim
2996  type(tree_node), pointer :: pnode
2997  double precision :: tmp(ixg^t), grad(ixg^t, ndim)
2998  double precision :: res
2999  double precision, parameter :: max_residual = 1d-3
3000  double precision, parameter :: residual_reduction = 1d-10
3001  integer, parameter :: max_its = 50
3002  double precision :: residual_it(max_its), max_divb
3003 
3004  mg%operator_type = mg_laplacian
3005 
3006  ! Set boundary conditions
3007  do n = 1, 2*ndim
3008  idim = (n+1)/2
3009  select case (typeboundary(mag(idim), n))
3010  case ('symm')
3011  ! d/dx B = 0, take phi = 0
3012  mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
3013  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
3014  case ('asymm')
3015  ! B = 0, so grad(phi) = 0
3016  mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
3017  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
3018  case ('cont')
3019  mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
3020  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
3021  case ('special')
3022  ! Assume Dirichlet boundary conditions, derivative zero
3023  mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
3024  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
3025  case ('periodic')
3026  ! Nothing to do here
3027  case default
3028  print *, "divb_multigrid warning: unknown b.c.: ", &
3029  trim(typeboundary(mag(idim), n))
3030  mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
3031  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
3032  end select
3033  end do
3034 
3035  ix^l=ixm^ll^ladd1;
3036  max_divb = 0.0d0
3037 
3038  ! Store divergence of B as right-hand side
3039  do iigrid = 1, igridstail
3040  igrid = igrids(iigrid);
3041  pnode => igrid_to_node(igrid, mype)%node
3042  id = pnode%id
3043  lvl = mg%boxes(id)%lvl
3044  nc = mg%box_size_lvl(lvl)
3045 
3046  ! Geometry subroutines expect this to be set
3047  block => ps(igrid)
3048  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
3049 
3050  call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^ll, ixm^ll, tmp, &
3052  mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(ixm^t)
3053  max_divb = max(max_divb, maxval(abs(tmp(ixm^t))))
3054  end do
3055 
3056  ! Solve laplacian(phi) = divB
3057  if(stagger_grid) then
3058  call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
3059  mpi_max, icomm, ierrmpi)
3060 
3061  if (mype == 0) print *, "Performing multigrid divB cleaning"
3062  if (mype == 0) print *, "iteration vs residual"
3063  ! Solve laplacian(phi) = divB
3064  do n = 1, max_its
3065  call mg_fas_fmg(mg, n>1, max_res=residual_it(n))
3066  if (mype == 0) write(*, "(I4,E11.3)") n, residual_it(n)
3067  if (residual_it(n) < residual_reduction * max_divb) exit
3068  end do
3069  if (mype == 0 .and. n > max_its) then
3070  print *, "divb_multigrid warning: not fully converged"
3071  print *, "current amplitude of divb: ", residual_it(max_its)
3072  print *, "multigrid smallest grid: ", &
3073  mg%domain_size_lvl(:, mg%lowest_lvl)
3074  print *, "note: smallest grid ideally has <= 8 cells"
3075  print *, "multigrid dx/dy/dz ratio: ", mg%dr(:, 1)/mg%dr(1, 1)
3076  print *, "note: dx/dy/dz should be similar"
3077  end if
3078  else
3079  do n = 1, max_its
3080  call mg_fas_vcycle(mg, max_res=res)
3081  if (res < max_residual) exit
3082  end do
3083  if (res > max_residual) call mpistop("divb_multigrid: no convergence")
3084  end if
3085 
3086 
3087  ! Correct the magnetic field
3088  do iigrid = 1, igridstail
3089  igrid = igrids(iigrid);
3090  pnode => igrid_to_node(igrid, mype)%node
3091  id = pnode%id
3092 
3093  ! Geometry subroutines expect this to be set
3094  block => ps(igrid)
3095  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
3096 
3097  ! Compute the gradient of phi
3098  tmp(ix^s) = mg%boxes(id)%cc({:,}, mg_iphi)
3099 
3100  if(stagger_grid) then
3101  do idim =1, ndim
3102  ixcmin^d=ixmlo^d-kr(idim,^d);
3103  ixcmax^d=ixmhi^d;
3104  call gradientx(tmp,ps(igrid)%x,ixg^ll,ixc^l,idim,grad(ixg^t,idim),.false.)
3105  ! Apply the correction B* = B - gradient(phi)
3106  ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
3107  end do
3108  ! store cell-center magnetic energy
3109  tmp(ixm^t) = sum(ps(igrid)%w(ixm^t, mag(1:ndim))**2, dim=ndim+1)
3110  ! change cell-center magnetic field
3111  call mhd_face_to_center(ixm^ll,ps(igrid))
3112  else
3113  do idim = 1, ndim
3114  call gradient(tmp,ixg^ll,ixm^ll,idim,grad(ixg^t, idim))
3115  end do
3116  ! store cell-center magnetic energy
3117  tmp(ixm^t) = sum(ps(igrid)%w(ixm^t, mag(1:ndim))**2, dim=ndim+1)
3118  ! Apply the correction B* = B - gradient(phi)
3119  ps(igrid)%w(ixm^t, mag(1:ndim)) = &
3120  ps(igrid)%w(ixm^t, mag(1:ndim)) - grad(ixm^t, :)
3121  end if
3122 
3123  if(total_energy) then
3124  ! Determine magnetic energy difference
3125  tmp(ixm^t) = 0.5_dp * (sum(ps(igrid)%w(ixm^t, &
3126  mag(1:ndim))**2, dim=ndim+1) - tmp(ixm^t))
3127  ! Keep thermal pressure the same
3128  ps(igrid)%w(ixm^t, e_) = ps(igrid)%w(ixm^t, e_) + tmp(ixm^t)
3129  end if
3130  end do
3131 
3132  active = .true.
3133 
3134  end subroutine mhd_clean_divb_multigrid
3135  }
3136 
3137  subroutine mhd_update_faces(ixI^L,ixO^L,qt,qdt,wprim,fC,fE,sCT,s)
3139 
3140  integer, intent(in) :: ixI^L, ixO^L
3141  double precision, intent(in) :: qt,qdt
3142  ! cell-center primitive variables
3143  double precision, intent(in) :: wprim(ixi^s,1:nw)
3144  type(state) :: sCT, s
3145  double precision, intent(in) :: fC(ixi^s,1:nwflux,1:ndim)
3146  double precision, intent(inout) :: fE(ixi^s,7-2*ndim:3)
3147 
3148  select case(type_ct)
3149  case('average')
3150  call update_faces_average(ixi^l,ixo^l,qt,qdt,fc,fe,sct,s)
3151  case('uct_contact')
3152  call update_faces_contact(ixi^l,ixo^l,qt,qdt,wprim,fc,fe,sct,s)
3153  case('uct_hll')
3154  call update_faces_hll(ixi^l,ixo^l,qt,qdt,fe,sct,s)
3155  case default
3156  call mpistop('choose average, uct_contact,or uct_hll for type_ct!')
3157  end select
3158 
3159  end subroutine mhd_update_faces
3160 
3161  !> get electric field though averaging neighors to update faces in CT
3162  subroutine update_faces_average(ixI^L,ixO^L,qt,qdt,fC,fE,sCT,s)
3165  use mod_usr_methods
3166 
3167  integer, intent(in) :: ixI^L, ixO^L
3168  double precision, intent(in) :: qt, qdt
3169  type(state) :: sCT, s
3170  double precision, intent(in) :: fC(ixi^s,1:nwflux,1:ndim)
3171  double precision, intent(inout) :: fE(ixi^s,7-2*ndim:3)
3172 
3173  integer :: hxC^L,ixC^L,jxC^L,ixCm^L
3174  integer :: idim1,idim2,idir,iwdim1,iwdim2
3175  double precision :: circ(ixi^s,1:ndim)
3176  ! current on cell edges
3177  double precision :: jce(ixi^s,7-2*ndim:3)
3178 
3179  associate(bfaces=>s%ws,x=>s%x)
3180 
3181  ! Calculate contribution to FEM of each edge,
3182  ! that is, estimate value of line integral of
3183  ! electric field in the positive idir direction.
3184  ixcmax^d=ixomax^d;
3185  ixcmin^d=ixomin^d-1;
3186 
3187  ! if there is resistivity, get eta J
3188  if(mhd_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,jce)
3189 
3190  fe=zero
3191 
3192  do idim1=1,ndim
3193  iwdim1 = mag(idim1)
3194  do idim2=1,ndim
3195  iwdim2 = mag(idim2)
3196  do idir=7-2*ndim,3! Direction of line integral
3197  ! Allow only even permutations
3198  if (lvc(idim1,idim2,idir)==1) then
3199  ! Assemble indices
3200  jxc^l=ixc^l+kr(idim1,^d);
3201  hxc^l=ixc^l+kr(idim2,^d);
3202  ! Interpolate to edges
3203  fe(ixc^s,idir)=quarter*(fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
3204  -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
3205 
3206  ! add current component of electric field at cell edges E=-vxB+eta J
3207  if(mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+jce(ixc^s,idir)
3208 
3209  fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
3210 
3211  if (.not.slab) then
3212  where(abs(x(ixc^s,r_)+half*dxlevel(r_))<1.0d-9)
3213  fe(ixc^s,idir)=zero
3214  end where
3215  end if
3216  end if
3217  end do
3218  end do
3219  end do
3220 
3221  ! allow user to change inductive electric field, especially for boundary driven applications
3222  if(associated(usr_set_electric_field)) &
3223  call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
3224 
3225  circ(ixi^s,1:ndim)=zero
3226 
3227  ! Calculate circulation on each face
3228 
3229  do idim1=1,ndim ! Coordinate perpendicular to face
3230  do idim2=1,ndim
3231  do idir=7-2*ndim,3 ! Direction of line integral
3232  ! Assemble indices
3233  hxc^l=ixc^l-kr(idim2,^d);
3234  ! Add line integrals in direction idir
3235  circ(ixc^s,idim1)=circ(ixc^s,idim1)&
3236  +lvc(idim1,idim2,idir)&
3237  *(fe(ixc^s,idir)&
3238  -fe(hxc^s,idir))
3239  end do
3240  end do
3241  end do
3242 
3243  ! Divide by the area of the face to get dB/dt
3244  do idim1=1,ndim
3245  ixcmax^d=ixomax^d;
3246  ixcmin^d=ixomin^d-kr(idim1,^d);
3247  where(s%surfaceC(ixc^s,idim1) > 1.0d-9*s%dvolume(ixc^s))
3248  circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
3249  elsewhere
3250  circ(ixc^s,idim1)=zero
3251  end where
3252  ! Time update
3253  bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
3254  end do
3255 
3256  end associate
3257 
3258  end subroutine update_faces_average
3259 
3260  !> update faces using UCT contact mode by Gardiner and Stone 2005 JCP 205, 509
3261  subroutine update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s)
3264  use mod_usr_methods
3265 
3266  integer, intent(in) :: ixI^L, ixO^L
3267  double precision, intent(in) :: qt, qdt
3268  ! cell-center primitive variables
3269  double precision, intent(in) :: wp(ixi^s,1:nw)
3270  type(state) :: sCT, s
3271  double precision, intent(in) :: fC(ixi^s,1:nwflux,1:ndim)
3272  double precision, intent(inout) :: fE(ixi^s,7-2*ndim:3)
3273 
3274  double precision :: circ(ixi^s,1:ndim)
3275  ! electric field at cell centers
3276  double precision :: ECC(ixi^s,7-2*ndim:3)
3277  ! gradient of E at left and right side of a cell face
3278  double precision :: EL(ixi^s),ER(ixi^s)
3279  ! gradient of E at left and right side of a cell corner
3280  double precision :: ELC(ixi^s),ERC(ixi^s)
3281  ! current on cell edges
3282  double precision :: jce(ixi^s,7-2*ndim:3)
3283  ! total magnetic field at cell centers
3284  double precision :: Btot(ixi^s,1:ndim)
3285  integer :: hxC^L,ixC^L,jxC^L,ixA^L,ixB^L
3286  integer :: idim1,idim2,idir,iwdim1,iwdim2
3287 
3288  associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm)
3289 
3290  if(b0field) then
3291  btot(ixi^s,1:ndim)=wp(ixi^s,mag(1:ndim))+block%B0(ixi^s,1:ndim,0)
3292  else
3293  btot(ixi^s,1:ndim)=wp(ixi^s,mag(1:ndim))
3294  end if
3295  ecc=0.d0
3296  ! Calculate electric field at cell centers
3297  do idim1=1,ndim; do idim2=1,ndim; do idir=7-2*ndim,3
3298  if(lvc(idim1,idim2,idir)==1)then
3299  ecc(ixi^s,idir)=ecc(ixi^s,idir)+btot(ixi^s,idim1)*wp(ixi^s,mom(idim2))
3300  else if(lvc(idim1,idim2,idir)==-1) then
3301  ecc(ixi^s,idir)=ecc(ixi^s,idir)-btot(ixi^s,idim1)*wp(ixi^s,mom(idim2))
3302  endif
3303  enddo; enddo; enddo
3304 
3305  ! if there is resistivity, get eta J
3306  if(mhd_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,jce)
3307 
3308  ! Calculate contribution to FEM of each edge,
3309  ! that is, estimate value of line integral of
3310  ! electric field in the positive idir direction.
3311  fe=zero
3312  ! evaluate electric field along cell edges according to equation (41)
3313  do idim1=1,ndim
3314  iwdim1 = mag(idim1)
3315  do idim2=1,ndim
3316  iwdim2 = mag(idim2)
3317  do idir=7-2*ndim,3 ! Direction of line integral
3318  ! Allow only even permutations
3319  if (lvc(idim1,idim2,idir)==1) then
3320  ixcmax^d=ixomax^d;
3321  ixcmin^d=ixomin^d+kr(idir,^d)-1;
3322  ! Assemble indices
3323  jxc^l=ixc^l+kr(idim1,^d);
3324  hxc^l=ixc^l+kr(idim2,^d);
3325  ! average cell-face electric field to cell edges
3326  fe(ixc^s,idir)=quarter*&
3327  (fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
3328  -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
3329 
3330  ! add slope in idim2 direction from equation (50)
3331  ixamin^d=ixcmin^d;
3332  ixamax^d=ixcmax^d+kr(idim1,^d);
3333  el(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(ixa^s,idir)
3334  hxc^l=ixa^l+kr(idim2,^d);
3335  er(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(hxc^s,idir)
3336  where(vnorm(ixc^s,idim1)>0.d0)
3337  elc(ixc^s)=el(ixc^s)
3338  else where(vnorm(ixc^s,idim1)<0.d0)
3339  elc(ixc^s)=el(jxc^s)
3340  else where
3341  elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
3342  end where
3343  hxc^l=ixc^l+kr(idim2,^d);
3344  where(vnorm(hxc^s,idim1)>0.d0)
3345  erc(ixc^s)=er(ixc^s)
3346  else where(vnorm(hxc^s,idim1)<0.d0)
3347  erc(ixc^s)=er(jxc^s)
3348  else where
3349  erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
3350  end where
3351  fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
3352 
3353  ! add slope in idim1 direction from equation (50)
3354  jxc^l=ixc^l+kr(idim2,^d);
3355  ixamin^d=ixcmin^d;
3356  ixamax^d=ixcmax^d+kr(idim2,^d);
3357  el(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(ixa^s,idir)
3358  hxc^l=ixa^l+kr(idim1,^d);
3359  er(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(hxc^s,idir)
3360  where(vnorm(ixc^s,idim2)>0.d0)
3361  elc(ixc^s)=el(ixc^s)
3362  else where(vnorm(ixc^s,idim2)<0.d0)
3363  elc(ixc^s)=el(jxc^s)
3364  else where
3365  elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
3366  end where
3367  hxc^l=ixc^l+kr(idim1,^d);
3368  where(vnorm(hxc^s,idim2)>0.d0)
3369  erc(ixc^s)=er(ixc^s)
3370  else where(vnorm(hxc^s,idim2)<0.d0)
3371  erc(ixc^s)=er(jxc^s)
3372  else where
3373  erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
3374  end where
3375  fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
3376 
3377  ! add current component of electric field at cell edges E=-vxB+eta J
3378  if(mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+jce(ixc^s,idir)
3379 
3380  ! times time step and edge length
3381  fe(ixc^s,idir)=fe(ixc^s,idir)*qdt*s%dsC(ixc^s,idir)
3382  if (.not.slab) then
3383  where(abs(x(ixc^s,r_)+half*dxlevel(r_))<1.0d-9)
3384  fe(ixc^s,idir)=zero
3385  end where
3386  end if
3387  end if
3388  end do
3389  end do
3390  end do
3391 
3392  ! allow user to change inductive electric field, especially for boundary driven applications
3393  if(associated(usr_set_electric_field)) &
3394  call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
3395 
3396  circ(ixi^s,1:ndim)=zero
3397 
3398  ! Calculate circulation on each face
3399  do idim1=1,ndim ! Coordinate perpendicular to face
3400  ixcmax^d=ixomax^d;
3401  ixcmin^d=ixomin^d-kr(idim1,^d);
3402  do idim2=1,ndim
3403  do idir=7-2*ndim,3 ! Direction of line integral
3404  ! Assemble indices
3405  hxc^l=ixc^l-kr(idim2,^d);
3406  ! Add line integrals in direction idir
3407  circ(ixc^s,idim1)=circ(ixc^s,idim1)&
3408  +lvc(idim1,idim2,idir)&
3409  *(fe(ixc^s,idir)&
3410  -fe(hxc^s,idir))
3411  end do
3412  end do
3413  ! Divide by the area of the face to get dB/dt
3414  ixcmax^d=ixomax^d;
3415  ixcmin^d=ixomin^d-kr(idim1,^d);
3416  where(s%surfaceC(ixc^s,idim1) > 1.0d-9*s%dvolume(ixc^s))
3417  circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
3418  elsewhere
3419  circ(ixc^s,idim1)=zero
3420  end where
3421  ! Time update cell-face magnetic field component
3422  bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
3423  end do
3424 
3425  end associate
3426 
3427  end subroutine update_faces_contact
3428 
3429  !> update faces
3430  subroutine update_faces_hll(ixI^L,ixO^L,qt,qdt,fE,sCT,s)
3433  use mod_usr_methods
3434 
3435  integer, intent(in) :: ixI^L, ixO^L
3436  double precision, intent(in) :: qt, qdt
3437  double precision, intent(inout) :: fE(ixi^s,7-2*ndim:3)
3438  type(state) :: sCT, s
3439 
3440  double precision :: vtilL(ixi^s,2)
3441  double precision :: vtilR(ixi^s,2)
3442  double precision :: bfacetot(ixi^s,ndim)
3443  double precision :: btilL(s%ixgs^s,ndim)
3444  double precision :: btilR(s%ixgs^s,ndim)
3445  double precision :: cp(ixi^s,2)
3446  double precision :: cm(ixi^s,2)
3447  double precision :: circ(ixi^s,1:ndim)
3448  ! current on cell edges
3449  double precision :: jce(ixi^s,7-2*ndim:3)
3450  integer :: hxC^L,ixC^L,ixCp^L,jxC^L,ixCm^L
3451  integer :: idim1,idim2,idir
3452 
3453  associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
3454  cbarmax=>vcts%cbarmax)
3455 
3456  ! Calculate contribution to FEM of each edge,
3457  ! that is, estimate value of line integral of
3458  ! electric field in the positive idir direction.
3459 
3460  ! Loop over components of electric field
3461 
3462  ! idir: electric field component we need to calculate
3463  ! idim1: directions in which we already performed the reconstruction
3464  ! idim2: directions in which we perform the reconstruction
3465 
3466  ! if there is resistivity, get eta J
3467  if(mhd_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,jce)
3468 
3469  fe=zero
3470 
3471  do idir=7-2*ndim,3
3472  ! Indices
3473  ! idir: electric field component
3474  ! idim1: one surface
3475  ! idim2: the other surface
3476  ! cyclic permutation: idim1,idim2,idir=1,2,3
3477  ! Velocity components on the surface
3478  ! follow cyclic premutations:
3479  ! Sx(1),Sx(2)=y,z ; Sy(1),Sy(2)=z,x ; Sz(1),Sz(2)=x,y
3480 
3481  ixcmax^d=ixomax^d;
3482  ixcmin^d=ixomin^d-1+kr(idir,^d);
3483 
3484  ! Set indices and directions
3485  idim1=mod(idir,3)+1
3486  idim2=mod(idir+1,3)+1
3487 
3488  jxc^l=ixc^l+kr(idim1,^d);
3489  ixcp^l=ixc^l+kr(idim2,^d);
3490 
3491  ! Reconstruct transverse transport velocities
3492  call reconstruct(ixi^l,ixc^l,idim2,vbarc(ixi^s,idim1,1),&
3493  vtill(ixi^s,2),vtilr(ixi^s,2))
3494 
3495  call reconstruct(ixi^l,ixc^l,idim1,vbarc(ixi^s,idim2,2),&
3496  vtill(ixi^s,1),vtilr(ixi^s,1))
3497 
3498  ! Reconstruct magnetic fields
3499  ! Eventhough the arrays are larger, reconstruct works with
3500  ! the limits ixG.
3501  if(b0field) then
3502  bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)+block%B0(ixi^s,idim1,idim1)
3503  bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)+block%B0(ixi^s,idim2,idim2)
3504  else
3505  bfacetot(ixi^s,idim1)=bfacesct(ixi^s,idim1)
3506  bfacetot(ixi^s,idim2)=bfacesct(ixi^s,idim2)
3507  end if
3508  call reconstruct(ixi^l,ixc^l,idim2,bfacetot(ixi^s,idim1),&
3509  btill(ixi^s,idim1),btilr(ixi^s,idim1))
3510 
3511  call reconstruct(ixi^l,ixc^l,idim1,bfacetot(ixi^s,idim2),&
3512  btill(ixi^s,idim2),btilr(ixi^s,idim2))
3513 
3514  ! Take the maximum characteristic
3515 
3516  cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
3517  cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
3518 
3519  cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
3520  cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
3521 
3522 
3523  ! Calculate eletric field
3524  fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
3525  + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
3526  - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
3527  /(cp(ixc^s,1)+cm(ixc^s,1)) &
3528  +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
3529  + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
3530  - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
3531  /(cp(ixc^s,2)+cm(ixc^s,2))
3532 
3533  ! add current component of electric field at cell edges E=-vxB+eta J
3534  if(mhd_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+jce(ixc^s,idir)
3535 
3536  fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
3537 
3538  if (.not.slab) then
3539  where(abs(x(ixc^s,r_)+half*dxlevel(r_)).lt.1.0d-9)
3540  fe(ixc^s,idir)=zero
3541  end where
3542  end if
3543 
3544  end do
3545 
3546  ! allow user to change inductive electric field, especially for boundary driven applications
3547  if(associated(usr_set_electric_field)) &
3548  call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
3549 
3550  circ(ixi^s,1:ndim)=zero
3551 
3552  ! Calculate circulation on each face: interal(fE dot dl)
3553 
3554  do idim1=1,ndim ! Coordinate perpendicular to face
3555  ixcmax^d=ixomax^d;
3556  ixcmin^d=ixomin^d-kr(idim1,^d);
3557  do idim2=1,ndim
3558  do idir=7-2*ndim,3 ! Direction of line integral
3559  ! Assemble indices
3560  hxc^l=ixc^l-kr(idim2,^d);
3561  ! Add line integrals in direction idir
3562  circ(ixc^s,idim1)=circ(ixc^s,idim1)&
3563  +lvc(idim1,idim2,idir)&
3564  *(fe(ixc^s,idir)&
3565  -fe(hxc^s,idir))
3566  end do
3567  end do
3568  end do
3569 
3570  ! Divide by the area of the face to get dB/dt
3571  do idim1=1,ndim
3572  ixcmax^d=ixomax^d;
3573  ixcmin^d=ixomin^d-kr(idim1,^d);
3574  where(s%surfaceC(ixc^s,idim1) > 1.0d-9*s%dvolume(ixc^s))
3575  circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
3576  elsewhere
3577  circ(ixc^s,idim1)=zero
3578  end where
3579  ! Time update
3580  bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
3581  end do
3582 
3583  end associate
3584  end subroutine update_faces_hll
3585 
3586  !> calculate eta J at cell edges
3587  subroutine get_resistive_electric_field(ixI^L,ixO^L,sCT,s,jce)
3589  use mod_usr_methods
3590  use mod_geometry
3591 
3592  integer, intent(in) :: ixI^L, ixO^L
3593  type(state), intent(in) :: sCT, s
3594  ! current on cell edges
3595  double precision :: jce(ixi^s,7-2*ndim:3)
3596 
3597  ! current on cell centers
3598  double precision :: jcc(ixi^s,7-2*ndim:3)
3599  ! location at cell faces
3600  double precision :: xs(ixgs^t,1:ndim)
3601  ! resistivity
3602  double precision :: eta(ixi^s)
3603  double precision :: gradi(ixgs^t)
3604  integer :: ix^D,ixC^L,ixA^L,ixB^L,idir,idirmin,idim1,idim2
3605 
3606  associate(x=>s%x,dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
3607  ! calculate current density at cell edges
3608  jce=0.d0
3609  do idim1=1,ndim
3610  do idim2=1,ndim
3611  do idir=7-2*ndim,3
3612  if (lvc(idim1,idim2,idir)==0) cycle
3613  ixcmax^d=ixomax^d;
3614  ixcmin^d=ixomin^d+kr(idir,^d)-1;
3615  ixbmax^d=ixcmax^d-kr(idir,^d)+1;
3616  ixbmin^d=ixcmin^d;
3617  ! current at transverse faces
3618  xs(ixb^s,:)=x(ixb^s,:)
3619  xs(ixb^s,idim2)=x(ixb^s,idim2)+half*dx(ixb^s,idim2)
3620  call gradientx(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi,.true.)
3621  if (lvc(idim1,idim2,idir)==1) then
3622  jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
3623  else
3624  jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
3625  end if
3626  end do
3627  end do
3628  end do
3629  ! get resistivity
3630  if(mhd_eta>zero)then
3631  jce(ixi^s,:)=jce(ixi^s,:)*mhd_eta
3632  else
3633  ixa^l=ixo^l^ladd1;
3634  call get_current(wct,ixi^l,ixo^l,idirmin,jcc)
3635  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,jcc,eta)
3636  ! calcuate eta on cell edges
3637  do idir=7-2*ndim,3
3638  ixcmax^d=ixomax^d;
3639  ixcmin^d=ixomin^d+kr(idir,^d)-1;
3640  jcc(ixc^s,idir)=0.d0
3641  {do ix^db=0,1\}
3642  if({ ix^d==1 .and. ^d==idir | .or.}) cycle
3643  ixamin^d=ixcmin^d+ix^d;
3644  ixamax^d=ixcmax^d+ix^d;
3645  jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
3646  {end do\}
3647  jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
3648  jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
3649  enddo
3650  end if
3651 
3652  end associate
3653  end subroutine get_resistive_electric_field
3654 
3655  !> calculate cell-center values from face-center values
3656  subroutine mhd_face_to_center(ixO^L,s)
3658  ! Non-staggered interpolation range
3659  integer, intent(in) :: ixO^L
3660  type(state) :: s
3661 
3662  integer :: fxO^L, gxO^L, hxO^L, jxO^L, kxO^L, idim
3663 
3664  associate(w=>s%w, ws=>s%ws)
3665 
3666  ! calculate cell-center values from face-center values in 2nd order
3667  do idim=1,ndim
3668  ! Displace index to the left
3669  ! Even if ixI^L is the full size of the w arrays, this is ok
3670  ! because the staggered arrays have an additional place to the left.
3671  hxo^l=ixo^l-kr(idim,^d);
3672  ! Interpolate to cell barycentre using arithmetic average
3673  ! This might be done better later, to make the method less diffusive.
3674  w(ixo^s,mag(idim))=half/s%surface(ixo^s,idim)*&
3675  (ws(ixo^s,idim)*s%surfaceC(ixo^s,idim)&
3676  +ws(hxo^s,idim)*s%surfaceC(hxo^s,idim))
3677  end do
3678 
3679  ! calculate cell-center values from face-center values in 4th order
3680  !do idim=1,ndim
3681  ! gxO^L=ixO^L-2*kr(idim,^D);
3682  ! hxO^L=ixO^L-kr(idim,^D);
3683  ! jxO^L=ixO^L+kr(idim,^D);
3684 
3685  ! ! Interpolate to cell barycentre using fourth order central formula
3686  ! w(ixO^S,mag(idim))=(0.0625d0/s%surface(ixO^S,idim))*&
3687  ! ( -ws(gxO^S,idim)*s%surfaceC(gxO^S,idim) &
3688  ! +9.0d0*ws(hxO^S,idim)*s%surfaceC(hxO^S,idim) &
3689  ! +9.0d0*ws(ixO^S,idim)*s%surfaceC(ixO^S,idim) &
3690  ! -ws(jxO^S,idim)*s%surfaceC(jxO^S,idim) )
3691  !end do
3692 
3693  ! calculate cell-center values from face-center values in 6th order
3694  !do idim=1,ndim
3695  ! fxO^L=ixO^L-3*kr(idim,^D);
3696  ! gxO^L=ixO^L-2*kr(idim,^D);
3697  ! hxO^L=ixO^L-kr(idim,^D);
3698  ! jxO^L=ixO^L+kr(idim,^D);
3699  ! kxO^L=ixO^L+2*kr(idim,^D);
3700 
3701  ! ! Interpolate to cell barycentre using sixth order central formula
3702  ! w(ixO^S,mag(idim))=(0.00390625d0/s%surface(ixO^S,idim))* &
3703  ! ( +3.0d0*ws(fxO^S,idim)*s%surfaceC(fxO^S,idim) &
3704  ! -25.0d0*ws(gxO^S,idim)*s%surfaceC(gxO^S,idim) &
3705  ! +150.0d0*ws(hxO^S,idim)*s%surfaceC(hxO^S,idim) &
3706  ! +150.0d0*ws(ixO^S,idim)*s%surfaceC(ixO^S,idim) &
3707  ! -25.0d0*ws(jxO^S,idim)*s%surfaceC(jxO^S,idim) &
3708  ! +3.0d0*ws(kxO^S,idim)*s%surfaceC(kxO^S,idim) )
3709  !end do
3710 
3711  end associate
3712 
3713  end subroutine mhd_face_to_center
3714 
3715  !> calculate magnetic field from vector potential
3716  subroutine b_from_vector_potential(ixIs^L, ixI^L, ixO^L, ws, x)
3719 
3720  integer, intent(in) :: ixIs^L, ixI^L, ixO^L
3721  double precision, intent(inout) :: ws(ixis^s,1:nws)
3722  double precision, intent(in) :: x(ixi^s,1:ndim)
3723 
3724  double precision :: Adummy(ixis^s,1:3)
3725 
3726  call b_from_vector_potentiala(ixis^l, ixi^l, ixo^l, ws, x, adummy)
3727 
3728  end subroutine b_from_vector_potential
3729 
3730 end module mod_mhd_phys
integer, parameter cylindrical
Definition: mod_geometry.t:9
type(iw_methods), dimension(max_nw) phys_iw_methods
Special methods defined per variable.
Definition: mod_physics.t:57
procedure(special_resistivity), pointer usr_special_resistivity
double precision function, dimension(ixo^s) mhd_mag_en(w, ixIL, ixOL)
Compute evolving magnetic energy.
This module contains definitions of global parameters and variables and some generic functions/subrou...
subroutine viscosity_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
subroutine b_from_vector_potentiala(ixIsL, ixIL, ixOL, ws, x, A)
calculate magnetic field from vector potential A at cell edges
procedure(sub_global_source), pointer phys_global_source_after
Definition: mod_physics.t:76
subroutine, public mhd_clean_divb_multigrid(qdt, qt, active)
subroutine gravity_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
Definition: mod_gravity.t:74
logical, public, protected mhd_magnetofriction
Whether magnetofriction is added.
Definition: mod_mhd_phys.t:29
subroutine, public small_values_error(w, x, ixIL, ixOL, w_flag, subname, smallw)
integer, public, protected e_
Index of the energy density (-1 if not present)
Definition: mod_mhd_phys.t:72
integer, public, protected p_
Index of the gas pressure (-1 if not present) should equal e_.
Definition: mod_mhd_phys.t:75
procedure(sub_energy_synchro), pointer phys_energy_synchro
Definition: mod_physics.t:70
subroutine, public get_normalized_divb(w, ixIL, ixOL, divb)
get dimensionless div B = |divB| * volume / area / |B|
logical small_values_use_primitive
Use primitive variables when correcting small values.
subroutine viscosity_add_source(qdt, ixIL, ixOL, wCT, w, x, energy, qsourcesplit, active)
Definition: mod_viscosity.t:86
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, conservative.
subroutine fixdivb_boundary(ixGL, ixOL, w, x, iB)
subroutine mhd_angmomfix(fC, x, wnew, ixIL, ixOL, idim)
Definition: mod_mhd_phys.t:232
double precision, public mhd_eta
The MHD resistivity.
Definition: mod_mhd_phys.t:97
double precision unit_length
Physical scaling factor for length.
The data structure that contains information about a tree node/grid block.
Definition: mod_forest.t:11
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
double precision unit_density
Physical scaling factor for density.
character(len=name_len) physics_type
String describing the physics type of the simulation.
Definition: mod_physics.t:14
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
Definition: mod_geometry.t:310
subroutine mhd_modify_wlr(ixIL, ixOL, qt, wLC, wRC, wLp, wRp, s, idir)
subroutine divvectors(qvec, ixIL, ixOL, divq)
Calculate divergence of a vector qvec within ixL using limited extrapolation to cell edges...
Definition: mod_geometry.t:561
logical phys_total_energy
Solve total energy equation or not.
Definition: mod_physics.t:28
subroutine radiative_cooling_add_source(qdt, ixIL, ixOL, wCT, w, x, qsourcesplit, active)
subroutine mhd_add_source_geom(qdt, ixIL, ixOL, wCT, w, x)
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:354
logical need_global_cmax
need global maximal wave speed
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
integer ndir
Number of spatial dimensions (components) for vector variables.
logical, public, protected mhd_viscosity
Whether viscosity is added.
Definition: mod_mhd_phys.t:17
integer nstep
How many sub-steps the time integrator takes.
subroutine e_to_rhos(ixIL, ixOL, w, x)
Convert energy to entropy.
Definition: mod_mhd_phys.t:819
integer coordinate
Definition: mod_geometry.t:6
procedure(sub_check_w), pointer phys_check_w
Definition: mod_physics.t:78
procedure(sub_face_to_center), pointer phys_face_to_center
Definition: mod_physics.t:85
logical, public, protected mhd_gravity
Whether gravity is added.
Definition: mod_mhd_phys.t:20
procedure(sub_small_values), pointer phys_handle_small_values
Definition: mod_physics.t:83
double precision unit_mass
Physical scaling factor for mass.
Module with basic grid data structures.
Definition: mod_forest.t:2
procedure(sub_get_a2max), pointer phys_get_a2max
Definition: mod_physics.t:66
subroutine cooling_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
integer istep
Index of the sub-step in a multi-step time integrator.
character(len=std_len) typeboundspeed
Which type of TVDLF method to use.
integer, dimension(2 *^nd), public, protected boundary_divbfix_skip
To skip * layer of ghost cells during divB=0 fix for boundary.
Definition: mod_mhd_phys.t:145
subroutine mhd_getv_hall(w, x, ixIL, ixOL, vHall)
procedure(sub_boundary_adjust), pointer phys_boundary_adjust
Definition: mod_physics.t:80
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...
logical phys_req_diagonal
Whether the physics routines require diagonal ghost cells, for example for computing a curl...
Definition: mod_physics.t:22
type(mg_t) mg
Data structure containing the multigrid tree.
subroutine update_faces_average(ixIL, ixOL, qt, qdt, fC, fE, sCT, s)
get electric field though averaging neighors to update faces in CT
logical b0field
split magnetic field as background B0 field
subroutine mhd_get_csound(w, x, ixIL, ixOL, idim, csound)
Calculate fast magnetosonic wave speed.
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor, but its use is kept to a minimum.
subroutine mhd_get_a2max(w, x, ixIL, ixOL, a2max)
Definition: mod_mhd_phys.t:896
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
Definition: mod_mhd_phys.t:139
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
double precision small_density
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, non-conservative. If the fourthorder precompiler flag is set, uses fourth order central difference for the laplacian. Then the stencil is 5 (2 neighbours).
subroutine particles_init()
Initialize particle data and parameters.
Definition: mod_particles.t:15
logical, public, protected mhd_energy
Whether an energy equation is used.
Definition: mod_mhd_phys.t:8
logical, public, protected mhd_divb_4thorder
Whether divB is computed with a fourth order approximation.
Definition: mod_mhd_phys.t:118
subroutine, public small_values_average(ixIL, ixOL, w, x, w_flag)
double precision, public mhd_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
Definition: mod_mhd_phys.t:45
procedure(sub_update_faces), pointer phys_update_faces
Definition: mod_physics.t:84
integer, public, protected psi_
Indices of the GLM psi.
Definition: mod_mhd_phys.t:81
subroutine mhd_boundary_adjust
integer nghostcells
Number of ghost cells surrounding a grid.
subroutine radiative_cooling_init(phys_gamma, He_abund)
Radiative cooling initialization.
double precision, public mhd_gamma
The adiabatic index.
Definition: mod_mhd_phys.t:91
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
Module with all the methods that users can customize in AMRVAC.
subroutine add_source_powel(qdt, ixIL, ixOL, wCT, w, x)
Add divB related sources to w within ixO corresponding to Powel.
procedure(sub_convert), pointer phys_e_to_ei
Definition: mod_physics.t:63
procedure(sub_get_tcutoff), pointer phys_get_tcutoff
Definition: mod_physics.t:67
logical, public, protected mhd_radiative_cooling
Whether radiative cooling is added.
Definition: mod_mhd_phys.t:14
subroutine add_source_b0split(qdt, ixIL, ixOL, wCT, w, x)
Source terms after split off time-independent magnetic field.
subroutine mhd_check_params
Definition: mod_mhd_phys.t:522
procedure(set_wlr), pointer usr_set_wlr
logical angmomfix
Enable to strictly conserve the angular momentum (works both in cylindrical and spherical coordinates...
integer, dimension(:), allocatable, public, protected mag
Indices of the magnetic field.
Definition: mod_mhd_phys.t:78
logical, dimension(:), allocatable small_values_fix_iw
Whether to apply small value fixes to certain variables.
double precision small_pressure
subroutine mhd_get_p_total(w, x, ixIL, ixOL, p)
Calculate total pressure within ixO^L including magnetic pressure.
subroutine mhd_add_source(qdt, ixIL, ixOL, wCT, w, x, qsourcesplit, active)
w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
logical stagger_grid
True for using stagger grid.
procedure(sub_angmomfix), pointer phys_angmomfix
Definition: mod_physics.t:82
Thermal conduction for HD and MHD.
procedure(sub_convert), pointer phys_to_conserved
Definition: mod_physics.t:60
integer ierrmpi
A global MPI error return code.
integer ixm
the mesh range (within a block with ghost cells)
subroutine, public mhd_get_v(w, x, ixIL, ixOL, v)
Calculate v vector.
Definition: mod_mhd_phys.t:855
character(len=std_len), dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
double precision unit_velocity
Physical scaling factor for velocity.
subroutine mhd_write_info(fh)
Write this module&#39;s parameters to a snapsoht.
Definition: mod_mhd_phys.t:215
subroutine mhd_get_cbounds(wLC, wRC, wLp, wRp, x, ixIL, ixOL, idim, cmax, cmin)
Estimating bounds for the minimum and maximum signal velocities.
Definition: mod_mhd_phys.t:981
procedure(set_electric_field), pointer usr_set_electric_field
subroutine mhd_getdt_hall(w, x, ixIL, ixOL, dxD, dthall)
subroutine, public get_current(w, ixIL, ixOL, idirmin, current)
Calculate idirmin and the idirmin:3 components of the common current array make sure that dxlevel(^D)...
logical, public, protected mhd_internal_e
Whether internal energy is solved instead of total energy.
Definition: mod_mhd_phys.t:38
double precision unit_charge
Physical scaling factor for charge.
subroutine, public mhd_to_conserved(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
Definition: mod_mhd_phys.t:616
subroutine mhd_get_csound_prim(w, x, ixIL, ixOL, idim, csound)
Calculate fast magnetosonic wave speed.
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
Definition: mod_mhd_phys.t:41
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
Definition: mod_mhd_phys.t:69
double precision, public mhd_etah
TODO: what is this?
Definition: mod_mhd_phys.t:103
subroutine divvector(qvec, ixIL, ixOL, divq, fourthorder, sixthorder)
Calculate divergence of a vector qvec within ixL.
Definition: mod_geometry.t:469
integer, parameter unitpar
integer, dimension(3, 3) kr
Kronecker delta tensor.
double precision, dimension(:,:), allocatable dx
integer, parameter flux_default
Indicates a normal flux.
Definition: mod_physics.t:41
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
Definition: mod_mhd_phys.t:142
character(len=std_len) typegrad
character(len=std_len), public, protected typedivbfix
Method type to clean divergence of B.
Definition: mod_mhd_phys.t:112
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
integer, dimension(:), allocatable, parameter d
procedure(sub_get_flux), pointer phys_get_flux
Definition: mod_physics.t:69
integer mype
The rank of the current MPI task.
subroutine, public get_divb(w, ixIL, ixOL, divb, fourthorder)
Calculate div B within ixO.
subroutine mhd_get_cmax(w, x, ixIL, ixOL, idim, cmax)
Calculate cmax_idim=csound+abs(v_idim) within ixO^L.
Definition: mod_mhd_phys.t:883
subroutine gravity_add_source(qdt, ixIL, ixOL, wCT, w, x, energy, qsourcesplit, active)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
Definition: mod_gravity.t:37
logical, public, protected mhd_solve_eaux
Whether auxiliary internal energy is solved.
Definition: mod_mhd_phys.t:35
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 for including gravity in (magneto)hydrodynamics simulations.
Definition: mod_gravity.t:2
subroutine mhd_physical_units()
Definition: mod_mhd_phys.t:544
subroutine rhos_to_e(ixIL, ixOL, w, x)
Convert entropy to energy.
Definition: mod_mhd_phys.t:837
subroutine add_source_janhunen(qdt, ixIL, ixOL, wCT, w, x)
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
procedure(sub_get_dt), pointer phys_get_dt
Definition: mod_physics.t:72
subroutine get_lorentz(ixIL, ixOL, w, JxB)
Compute the Lorentz force (JxB)
subroutine mhd_find_small_values(primitive, w, x, ixIL, ixOL, subname)
subroutine update_faces_hll(ixIL, ixOL, qt, qdt, fE, sCT, s)
update faces
logical phys_internal_e
Solve internal enery instead of total energy.
Definition: mod_physics.t:31
character(len=std_len), public, protected type_ct
Method type of constrained transport.
Definition: mod_mhd_phys.t:115
Module for handling problematic values in simulations, such as negative pressures.
double precision unit_pressure
Physical scaling factor for pressure.
double precision unit_magneticfield
Physical scaling factor for magnetic field.
logical use_particles
Use particles module or not.
double precision unit_temperature
Physical scaling factor for temperature.
integer, public, protected rho_
Index of the density (in the w array)
Definition: mod_mhd_phys.t:66
procedure(sub_add_source_geom), pointer phys_add_source_geom
Definition: mod_physics.t:73
subroutine get_resistive_electric_field(ixIL, ixOL, sCT, s, jce)
calculate eta J at cell edges
subroutine, public mhd_get_csound2(w, x, ixIL, ixOL, csound2)
Calculate the square of the thermal sound speed csound2 within ixO^L. csound2=gamma*p/rho.
Module containing all the particle routines.
Definition: mod_particles.t:2
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:198
subroutine, public mhd_face_to_center(ixOL, s)
calculate cell-center values from face-center values
double precision function, dimension(ixo^s), public mhd_kin_en(w, ixIL, ixOL, inv_rho)
compute kinetic energy
procedure(sub_check_params), pointer phys_check_params
Definition: mod_physics.t:59
integer, parameter flux_special
Indicates the flux should be specially treated.
Definition: mod_physics.t:47
subroutine add_source_glm(qdt, ixIL, ixOL, wCT, w, x)
double precision, public mhd_eta_hyper
The MHD hyper-resistivity.
Definition: mod_mhd_phys.t:100
double precision, dimension(ndim) dxlevel
double precision bdip
amplitude of background dipolar, quadrupolar, octupolar, user&#39;s field
subroutine mhd_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
If resistivity is not zero, check diffusion time limit for dt.
integer, parameter ndim
Number of spatial dimensions for grid variables.
integer, dimension(:), allocatable, public, protected tracer
Indices of the tracers.
Definition: mod_mhd_phys.t:88
subroutine gravity_init()
Initialize the module.
Definition: mod_gravity.t:27
logical, public, protected mhd_particles
Whether particles module is added.
Definition: mod_mhd_phys.t:26
integer phys_wider_stencil
To use wider stencils in flux calculations. A value of 1 will extend it by one cell in both direction...
Definition: mod_physics.t:18
subroutine, public mhd_get_pthermal(w, x, ixIL, ixOL, pth)
Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho-b**2/2) within ixO^L.
logical phys_solve_eaux
Solve internal energy and total energy equations.
Definition: mod_physics.t:34
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
subroutine add_source_linde(qdt, ixIL, ixOL, wCT, w, x)
integer, public, protected paux_
Definition: mod_mhd_phys.t:85
subroutine mhd_e_to_ei(ixIL, ixOL, w, x)
Transform total energy to internal energy.
Definition: mod_mhd_phys.t:703
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:411
logical use_multigrid
Use multigrid (only available in 2D and 3D)
subroutine mhd_handle_small_values(primitive, w, x, ixIL, ixOL, subname)
Definition: mod_mhd_phys.t:756
subroutine mhd_check_w(primitive, ixIL, ixOL, w, flag, smallw)
Definition: mod_mhd_phys.t:580
logical slab
Cartesian geometry or not.
subroutine, public mhd_phys_init()
Definition: mod_mhd_phys.t:245
integer icomm
The MPI communicator.
double precision c_norm
Normalised speed of light.
double precision function, dimension(ixo^s), public mhd_mag_en_all(w, ixIL, ixOL)
Compute 2 times total magnetic energy.
subroutine mhd_energy_synchro(ixIL, ixOL, w, x)
Definition: mod_mhd_phys.t:718
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
logical, public, protected mhd_4th_order
MHD fourth order.
Definition: mod_mhd_phys.t:60
procedure(sub_get_v_idim), pointer phys_get_v_idim
Definition: mod_physics.t:71
subroutine mhd_gamma2_alfven(ixIL, ixOL, w, gamma_A2)
Compute 1/(1+v_A^2/c^2) for Boris&#39; approximation, where v_A is the Alfven velocity.
character(len=std_len) typediv
logical, public divbwave
Add divB wave in Roe solver.
Definition: mod_mhd_phys.t:133
integer, public, protected mhd_n_tracer
Number of tracer species.
Definition: mod_mhd_phys.t:63
procedure(sub_get_cmax), pointer phys_get_cmax
Definition: mod_physics.t:65
logical phys_energy
Solve energy equation or not.
Definition: mod_physics.t:25
subroutine boris_add_source(qdt, ixIL, ixOL, wCT, w, x, qsourcesplit, active)
The module add viscous source terms and check time step.
Definition: mod_viscosity.t:10
subroutine, public b_from_vector_potential(ixIsL, ixIL, ixOL, ws, x)
calculate magnetic field from vector potential
procedure(sub_convert), pointer phys_ei_to_e
Definition: mod_physics.t:62
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
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:616
procedure(sub_add_source), pointer phys_add_source
Definition: mod_physics.t:74
subroutine update_faces_contact(ixIL, ixOL, qt, qdt, wp, fC, fE, sCT, s)
update faces using UCT contact mode by Gardiner and Stone 2005 JCP 205, 509
procedure(sub_convert), pointer phys_to_primitive
Definition: mod_physics.t:61
subroutine viscosity_init(phys_wider_stencil, phys_req_diagonal)
Initialize the module.
Definition: mod_viscosity.t:57
subroutine thermal_conduction_init(phys_gamma)
Initialize the module.
integer, parameter spherical
Definition: mod_geometry.t:10
module radiative cooling – add optically thin radiative cooling for HD and MHD
double precision unit_time
Physical scaling factor for time.
logical, public, protected mhd_thermal_conduction
Whether thermal conduction is used.
Definition: mod_mhd_phys.t:11
logical, public, protected mhd_glm
Whether GLM-MHD is used.
Definition: mod_mhd_phys.t:32
subroutine mhd_get_tcutoff(ixIL, ixOL, w, x, tco_local, Tmax_local)
get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
Definition: mod_mhd_phys.t:919
subroutine mhd_ei_to_e(ixIL, ixOL, w, x)
Transform internal energy to total energy.
Definition: mod_mhd_phys.t:687
logical, public, protected b0field_forcefree
B0 field is force-free.
Definition: mod_mhd_phys.t:148
subroutine internal_energy_add_source(qdt, ixIL, ixOL, wCT, w, x, ie)
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:38
logical, dimension(:), allocatable phyboundblock
True if a block has any physical boundary.
subroutine magnetofriction_init()
Initialize the module.
Magneto-hydrodynamics module.
Definition: mod_mhd_phys.t:2
integer, public, protected eaux_
Indices of auxiliary internal energy.
Definition: mod_mhd_phys.t:84
character(len=20), public small_values_method
How to handle small values.
module mod_magnetofriction.t Purpose: use magnetofrictional method to relax 3D magnetic field to forc...
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision function, dimension(ixo^s) mhd_gamma_alfven(w, ixIL, ixOL)
Compute 1/sqrt(1+v_A^2/c^2) for Boris simplification, where v_A is the Alfven velocity.
double precision, public mhd_adiab
The adiabatic constant.
Definition: mod_mhd_phys.t:94
double precision unit_numberdensity
Physical scaling factor for number density.
subroutine mhd_update_faces(ixIL, ixOL, qt, qdt, wprim, fC, fE, sCT, s)
logical, public, protected mhd_hall
Whether Hall-MHD is used.
Definition: mod_mhd_phys.t:23
procedure(sub_get_pthermal), pointer phys_get_pthermal
Definition: mod_physics.t:79
logical, public clean_initial_divb
clean initial divB
Definition: mod_mhd_phys.t:136
subroutine, public mhd_to_primitive(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
Definition: mod_mhd_phys.t:650
procedure(sub_get_cbounds), pointer phys_get_cbounds
Definition: mod_physics.t:68
procedure(sub_modify_wlr), pointer phys_modify_wlr
Definition: mod_physics.t:64
procedure(sub_write_info), pointer phys_write_info
Definition: mod_physics.t:81
integer, parameter flux_tvdlf
Indicates the flux should be treated with tvdlf.
Definition: mod_physics.t:43
subroutine mhd_get_flux(wC, w, x, ixIL, ixOL, idim, f)
Calculate fluxes within ixO^L.
pure subroutine cross_product(ixIL, ixOL, a, b, axb)
Cross product of two vectors.
double precision function, dimension(ixo^s) mhd_mag_i_all(w, ixIL, ixOL, idir)
Compute full magnetic field by direction.
subroutine, public mhd_get_v_idim(w, x, ixIL, ixOL, idim, v)
Calculate v component.
Definition: mod_mhd_phys.t:871