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