MPI-AMRVAC  2.0
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  !> MHD fourth order
42  logical, public, protected :: mhd_4th_order = .false.
43 
44  !> Number of tracer species
45  integer, public, protected :: mhd_n_tracer = 0
46 
47  !> Index of the density (in the w array)
48  integer, public, protected :: rho_
49 
50  !> Indices of the momentum density
51  integer, allocatable, public, protected :: mom(:)
52 
53  !> Index of the energy density (-1 if not present)
54  integer, public, protected :: e_
55 
56  !> Index of the gas pressure (-1 if not present) should equal e_
57  integer, public, protected :: p_
58 
59  !> Indices of the magnetic field
60  integer, allocatable, public, protected :: mag(:)
61 
62  !> Indices of the GLM psi
63  integer, public, protected :: psi_
64 
65  !> Indices of the tracers
66  integer, allocatable, public, protected :: tracer(:)
67 
68  !> The adiabatic index
69  double precision, public :: mhd_gamma = 5.d0/3.0d0
70 
71  !> The adiabatic constant
72  double precision, public :: mhd_adiab = 1.0d0
73 
74  !> The MHD resistivity
75  double precision, public :: mhd_eta = 0.0d0
76 
77  !> The MHD hyper-resistivity
78  double precision, public :: mhd_eta_hyper = 0.0d0
79 
80  !> TODO: what is this?
81  double precision, public :: mhd_etah = 0.0d0
82 
83  !> The small_est allowed energy
84  double precision, protected :: small_e
85 
86  !> The number of waves
87  integer :: nwwave=8
88 
89  !> Method type to clean divergence of B
90  character(len=std_len), public, protected :: typedivbfix = 'linde'
91 
92  !> Whether divB is computed with a fourth order approximation
93  logical, public, protected :: mhd_divb_4thorder = .false.
94 
95  !> Method type in a integer for good performance
96  integer :: type_divb
97 
98  !> Coefficient of diffusive divB cleaning
99  double precision :: divbdiff = 0.8d0
100 
101  !> Update all equations due to divB cleaning
102  character(len=std_len) :: typedivbdiff = 'all'
103 
104  !> Use a compact way to add resistivity
105  logical :: compactres = .false.
106 
107  !> Add divB wave in Roe solver
108  logical, public :: divbwave = .true.
109 
110  !> Helium abundance over Hydrogen
111  double precision, public, protected :: he_abundance=0.1d0
112 
113  !> To control divB=0 fix for boundary
114  logical, public, protected :: boundary_divbfix(2*^nd)=.true.
115 
116  !> To skip * layer of ghost cells during divB=0 fix for boundary
117  integer, public, protected :: boundary_divbfix_skip(2*^nd)=0
118 
119  !> B0 field is force-free
120  logical, public, protected :: b0field_forcefree=.true.
121 
122  !> gamma minus one and its inverse
123  double precision :: gamma_1, inv_gamma_1
124 
125  ! DivB cleaning methods
126  integer, parameter :: divb_none = 0
127  integer, parameter :: divb_multigrid = -1
128  integer, parameter :: divb_glm1 = 1
129  integer, parameter :: divb_glm2 = 2
130  integer, parameter :: divb_powel = 3
131  integer, parameter :: divb_janhunen = 4
132  integer, parameter :: divb_linde = 5
133  integer, parameter :: divb_lindejanhunen = 6
134  integer, parameter :: divb_lindepowel = 7
135  integer, parameter :: divb_lindeglm = 8
136 
137 
138  ! Public methods
139  public :: mhd_phys_init
140  public :: mhd_kin_en
141  public :: mhd_get_pthermal
142  public :: mhd_get_v
143  public :: mhd_to_conserved
144  public :: mhd_to_primitive
145  public :: mhd_get_csound2
146  public :: get_divb
147  public :: get_current
148  public :: get_normalized_divb
149 
150 contains
151 
152  !> Read this module"s parameters from a file
153  subroutine mhd_read_params(files)
155  character(len=*), intent(in) :: files(:)
156  integer :: n
157 
158  namelist /mhd_list/ mhd_energy, mhd_n_tracer, mhd_gamma, mhd_adiab,&
162  typedivbdiff, compactres, divbwave, he_abundance, si_unit, b0field,&
165 
166  do n = 1, size(files)
167  open(unitpar, file=trim(files(n)), status="old")
168  read(unitpar, mhd_list, end=111)
169 111 close(unitpar)
170  end do
171 
172  end subroutine mhd_read_params
173 
174  !> Write this module's parameters to a snapsoht
175  subroutine mhd_write_info(fh)
177  integer, intent(in) :: fh
178  integer, parameter :: n_par = 1
179  double precision :: values(n_par)
180  character(len=name_len) :: names(n_par)
181  integer, dimension(MPI_STATUS_SIZE) :: st
182  integer :: er
183 
184  call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
185 
186  names(1) = "gamma"
187  values(1) = mhd_gamma
188  call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
189  call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
190  end subroutine mhd_write_info
191 
192  subroutine mhd_angmomfix(fC,x,wnew,ixI^L,ixO^L,idim)
194  double precision, intent(in) :: x(ixi^s,1:ndim)
195  double precision, intent(inout) :: fC(ixi^s,1:nwflux,1:ndim), wnew(ixi^s,1:nw)
196  integer, intent(in) :: ixI^L, ixO^L
197  integer, intent(in) :: idim
198  integer :: hxO^L, kxC^L, iw
199  double precision :: inv_volume(ixi^s)
200 
201  call mpistop("to do")
202  ! ! shifted indexes
203  ! hxO^L=ixO^L-kr(idim,^D);
204  ! ! all the indexes
205  ! kxCmin^D=hxOmin^D;
206  ! kxCmax^D=ixOmax^D;
207  !
208  ! inv_volume = 1.0d0/block%dvolume(ixO^S)
209  !
210  ! select case(typeaxial)
211  ! case ("cylindrical")
212  ! do iw=1,nwflux
213  ! if (idim==r_ .and. iw==iw_mom(phi_)) then
214  ! fC(kxC^S,iw,idim)= fC(kxC^S,iw,idim)*(x(kxC^S,r_)+half*block%dx(kxC^S,r_))
215  ! wnew(ixO^S,iw)=wnew(ixO^S,iw) + (fC(ixO^S,iw,idim)-fC(hxO^S,iw,idim)) * &
216  ! (inv_volume/x(ixO^S,r_))
217  ! else
218  ! wnew(ixO^S,iw)=wnew(ixO^S,iw) + (fC(ixO^S,iw,idim)-fC(hxO^S,iw,idim)) * &
219  ! inv_volume
220  ! endif
221  ! enddo
222  ! case ("spherical")
223  ! do iw=1,nwflux
224  ! if (idim==r_ .and. (iw==iw_mom(2) .or. iw==iw_mom(phi_))) then
225  ! fC(kxC^S,iw,idim)= fC(kxC^S,iw,idim)*(x(kxC^S,r_)+half*block%dx(kxC^S,r_))
226  ! wnew(ixO^S,iw)=wnew(ixO^S,iw) + (fC(ixO^S,iw,idim)-fC(hxO^S,iw,idim)) * &
227  ! (inv_volume/x(ixO^S,r_))
228  ! elseif (idim==2 .and. iw==iw_mom(phi_)) then
229  ! fC(kxC^S,iw,idim)=fC(kxC^S,iw,idim)*dsin(x(kxC^S,2)+half*block%dx(kxC^S,2)) ! (x(4,3,1)-x(3,3,1)))
230  ! wnew(ixO^S,iw)=wnew(ixO^S,iw) + (fC(ixO^S,iw,idim)-fC(hxO^S,iw,idim)) * &
231  ! (inv_volume/dsin(x(ixO^S,2)))
232  ! else
233  ! wnew(ixO^S,iw)=wnew(ixO^S,iw) + (fC(ixO^S,iw,idim)-fC(hxO^S,iw,idim)) * &
234  ! inv_volume
235  ! endif
236  ! enddo
237  !
238  ! ! if (idim==r_) then
239  ! ! fC(kxC^S,iw_mom(phi_),idim)= fC(kxC^S,iw_mom(phi_),idim)*(x(kxC^S,r_)+half*block%dx(kxC^S,r_))
240  ! ! fC(kxC^S,iw_mom(phi_),idim)= fC(kxC^S,iw_mom(phi_),idim)*(x(kxC^S,r_)+half*block%dx(kxC^S,r_))
241  ! ! wnew(ixO^S,iw_mom(phi_))=wnew(ixO^S,iw_mom(phi_)) + (fC(ixO^S,iw_mom(phi_),idim)-fC(hxO^S,iw_mom(phi_),idim)) * &
242  ! ! (inv_volume/x(ixO^S,r_))
243  ! !
244  ! ! elseif (idim==2) then
245  ! ! fC(hxOmin1:hxOmax1,hxOmin2,hxOmin3:hxOmax3,iw,idim)=fC(hxOmin1:hxOmax1,hxOmin2,hxOmin3:hxOmax3,iw,idim)*dsin(x(hxOmin1:hxOmax1,hxOmin2,hxOmin3:hxOmax3,2)+half*block%dx(hxOmin1:hxOmax1,hxOmin2,hxOmin3:hxOmax3,2)) ! (x(4,3,1)-x(3,3,1)))
246  ! ! fC(ixO^S,iw,idim)=fC(ixO^S,iw,idim)*dsin(x(ixO^S,2)+half*block%dx(ixO^S,2)) ! (x(4,3,1)-x(3,3,1)))
247  ! ! wnew(ixO^S,iw)=wnew(ixO^S,iw) + (fC(ixO^S,iw,idim)-fC(hxO^S,iw,idim)) * &
248  ! ! (inv_volume/dsin(x(ixO^S,2)))
249  ! ! endif
250  !
251  ! end select
252 
253  end subroutine mhd_angmomfix
254 
255  subroutine mhd_phys_init()
259  use mod_viscosity, only: viscosity_init
260  use mod_gravity, only: gravity_init
261  use mod_particles, only: particles_init
263  use mod_physics
264  {^nooned
266  }
267 
268  integer :: itr, idir
269 
270  call mhd_read_params(par_files)
271 
272  physics_type = "mhd"
274  ! set default gamma for polytropic/isothermal process
275  if(.not.mhd_energy) mhd_gamma=1.d0
277  if(ndim==1) typedivbfix='none'
278  select case (typedivbfix)
279  case ('none')
280  type_divb = divb_none
281  {^nooned
282  case ('multigrid')
283  type_divb = divb_multigrid
284  use_multigrid = .true.
285  mg%operator_type = mg_laplacian
287  }
288  case ('glm1')
289  mhd_glm = .true.
290  need_global_cmax = .true.
291  type_divb = divb_glm1
292  case ('glm2')
293  mhd_glm = .true.
294  need_global_cmax = .true.
295  need_global_vmax = .true.
296  type_divb = divb_glm2
297  case ('powel', 'powell')
298  type_divb = divb_powel
299  case ('janhunen')
300  type_divb = divb_janhunen
301  case ('linde')
302  type_divb = divb_linde
303  case ('lindejanhunen')
304  type_divb = divb_lindejanhunen
305  case ('lindepowel')
306  type_divb = divb_lindepowel
307  case ('lindeglm')
308  mhd_glm = .true.
309  need_global_cmax = .true.
310  need_global_vmax = .true.
311  type_divb = divb_lindeglm
312  case default
313  call mpistop('Unknown divB fix')
314  end select
315 
316  ! Determine flux variables
317  rho_ = var_set_rho()
318 
319  allocate(mom(ndir))
320  mom(:) = var_set_momentum(ndir)
321 
322  ! Set index of energy variable
323  if (mhd_energy) then
324  nwwave = 8
325  e_ = var_set_energy() ! energy density
326  p_ = e_ ! gas pressure
327  else
328  nwwave = 7
329  e_ = -1
330  p_ = -1
331  end if
332 
333  allocate(mag(ndir))
334  mag(:) = var_set_bfield(ndir)
335 
336  if (mhd_glm) then
337  psi_ = var_set_fluxvar('psi', 'psi', need_bc=.false.)
338  else
339  psi_ = -1
340  end if
341 
342  allocate(tracer(mhd_n_tracer))
343 
344  ! Set starting index of tracers
345  do itr = 1, mhd_n_tracer
346  tracer(itr) = var_set_fluxvar("trc", "trp", itr, need_bc=.false.)
347  end do
348 
349  nvector = 2 ! No. vector vars
350  allocate(iw_vector(nvector))
351  iw_vector(1) = mom(1) - 1 ! TODO: why like this?
352  iw_vector(2) = mag(1) - 1 ! TODO: why like this?
353 
354  ! Check whether custom flux types have been defined
355  if (.not. allocated(flux_type)) then
356  allocate(flux_type(ndir, nw))
358  else if (any(shape(flux_type) /= [ndir, nw])) then
359  call mpistop("phys_check error: flux_type has wrong shape")
360  end if
361  do idir=1,ndir
362  if(ndim>1) flux_type(idir,mag(idir))=flux_tvdlf
363  end do
364  if(mhd_glm .and. ndim>1) flux_type(:,psi_)=flux_tvdlf
365 
382 
383  ! Whether diagonal ghost cells are required for the physics
384  if(type_divb < divb_linde) phys_req_diagonal = .false.
385 
386  ! derive units from basic units
387  call mhd_physical_units()
388 
389  if(.not. mhd_energy .and. mhd_thermal_conduction) then
390  call mpistop("thermal conduction needs mhd_energy=T")
391  end if
392  if(.not. mhd_energy .and. mhd_radiative_cooling) then
393  call mpistop("radiative cooling needs mhd_energy=T")
394  end if
395 
396  ! initialize thermal conduction module
397  if (mhd_thermal_conduction) then
398  phys_req_diagonal = .true.
400  end if
401 
402  ! Initialize radiative cooling module
403  if (mhd_radiative_cooling) then
405  end if
406 
407  ! Initialize viscosity module
409 
410  ! Initialize gravity module
411  if(mhd_gravity) then
412  call gravity_init()
413  end if
414 
415  ! Initialize particles module
416  if(mhd_particles) then
417  call particles_init()
418  phys_req_diagonal = .true.
419  end if
420 
421  ! initialize magnetofriction module
422  if(mhd_magnetofriction) then
423  phys_req_diagonal = .true.
424  call magnetofriction_init()
425  end if
426 
427  if(type_divb==divb_glm1) then
428  ! Solve the Riemann problem for the linear 2x2 system for normal
429  ! B-field and GLM_Psi according to Dedner 2002:
431  end if
432 
433  ! For Hall, we need one more reconstructed layer since currents are computed
434  ! in getflux: assuming one additional ghost layer (two for FOURTHORDER) was
435  ! added in nghostcells.
436  if (mhd_hall) then
437  phys_req_diagonal = .true.
438  if (mhd_4th_order) then
440  else
442  end if
443  end if
444 
445  end subroutine mhd_phys_init
446 
447  subroutine mhd_check_params
449 
450  ! after user parameter setting
451  gamma_1=mhd_gamma-1.d0
452 
453  if (.not. mhd_energy) then
454  if (mhd_gamma <= 0.0d0) call mpistop ("Error: mhd_gamma <= 0")
455  if (mhd_adiab < 0.0d0) call mpistop ("Error: mhd_adiab < 0")
457  else
458  if (mhd_gamma <= 0.0d0 .or. mhd_gamma == 1.0d0) &
459  call mpistop ("Error: mhd_gamma <= 0 or mhd_gamma == 1")
460  inv_gamma_1=1.d0/gamma_1
461  small_e = small_pressure * inv_gamma_1
462  end if
463 
464  end subroutine mhd_check_params
465 
466  subroutine mhd_physical_units()
468  double precision :: mp,kB,miu0
469  ! Derive scaling units
470  if(si_unit) then
471  mp=mp_si
472  kb=kb_si
473  miu0=miu0_si
474  else
475  mp=mp_cgs
476  kb=kb_cgs
477  miu0=4.d0*dpi
478  end if
479  if(unit_velocity==0) then
485  else
491  end if
492 
493  end subroutine mhd_physical_units
494 
495  subroutine mhd_check_w(primitive,ixI^L,ixO^L,w,flag)
497 
498  logical, intent(in) :: primitive
499  integer, intent(in) :: ixI^L, ixO^L
500  double precision, intent(in) :: w(ixi^s,nw)
501  integer, intent(inout) :: flag(ixi^s)
502  double precision :: tmp(ixi^s)
503 
504  flag(ixo^s)=0
505  where(w(ixo^s, rho_) < small_density) flag(ixo^s) = rho_
506 
507  if (mhd_energy) then
508  if (block%e_is_internal) then
509  where(w(ixo^s, e_) < small_pressure*inv_gamma_1) flag(ixo^s) = e_
510  else
511  if (primitive)then
512  where(w(ixo^s, e_) < small_pressure) flag(ixo^s) = e_
513  else
514  ! Calculate pressure=(gamma-1)*(e-0.5*(2ek+2eb))
515  tmp(ixo^s)=w(ixo^s,e_) - &
516  mhd_kin_en(w,ixi^l,ixo^l)-mhd_mag_en(w,ixi^l,ixo^l)
517  if(type_divb==divb_glm2) tmp(ixo^s)=tmp(ixo^s)-0.5d0*w(ixo^s,psi_)**2
518  tmp(ixo^s)=gamma_1*tmp(ixo^s)
519  where(tmp(ixo^s) < small_pressure) flag(ixo^s) = e_
520  end if
521  end if
522  end if
523  end subroutine mhd_check_w
524 
525  !> Transform primitive variables into conservative ones
526  subroutine mhd_to_conserved(ixI^L,ixO^L,w,x)
528  integer, intent(in) :: ixI^L, ixO^L
529  double precision, intent(inout) :: w(ixi^s, nw)
530  double precision, intent(in) :: x(ixi^s, 1:ndim)
531  integer :: idir, itr
532 
534  call mhd_handle_small_values(.true., w, x, ixi^l, ixo^l, 'mhd_to_conserved')
535  end if
536 
537  if (mhd_energy) then
538  ! Calculate total energy from pressure, kinetic and magnetic energy
539  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1
540  if(.not.block%e_is_internal) w(ixo^s,e_)=w(ixo^s,e_) + &
541  0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) * w(ixo^s, rho_) + &
542  mhd_mag_en(w, ixi^l, ixo^l)
543  if(type_divb==divb_glm2) w(ixo^s,e_)=w(ixo^s,e_) + 0.5d0*w(ixo^s,psi_)**2
544  end if
545 
546  ! Convert velocity to momentum
547  do idir = 1, ndir
548  w(ixo^s, mom(idir)) = w(ixo^s, rho_) * w(ixo^s, mom(idir))
549  end do
550 
551  if (check_small_values .and. .not. small_values_use_primitive) then
552  call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_conserved')
553  end if
554  end subroutine mhd_to_conserved
555 
556  !> Transform conservative variables into primitive ones
557  subroutine mhd_to_primitive(ixI^L,ixO^L,w,x)
559  integer, intent(in) :: ixI^L, ixO^L
560  double precision, intent(inout) :: w(ixi^s, nw)
561  double precision, intent(in) :: x(ixi^s, 1:ndim)
562  double precision :: inv_rho(ixo^s)
563  integer :: itr, idir
564 
565  if (check_small_values .and. .not. small_values_use_primitive) then
566  call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive')
567  end if
568 
569  inv_rho = 1.0d0 / w(ixo^s, rho_)
570 
571  if (mhd_energy) then
572  ! Calculate pressure = (gamma-1) * (e-ek-eb)
573  if(.not.block%e_is_internal) then
574  w(ixo^s, p_) = w(ixo^s, e_) &
575  - mhd_kin_en(w, ixi^l, ixo^l, inv_rho) &
576  - mhd_mag_en(w, ixi^l, ixo^l)
577  ! Calculate pressure = (gamma-1) * (e-ek-eb-epsi)
578  if(type_divb==divb_glm2) w(ixo^s, p_)=w(ixo^s, p_)-0.5d0*w(ixo^s,psi_)**2
579  end if
580  w(ixo^s, p_) = gamma_1*w(ixo^s, p_)
581  end if
582 
583  ! Convert momentum to velocity
584  do idir = 1, ndir
585  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))*inv_rho
586  end do
587 
589  call mhd_handle_small_values(.true., w, x, ixi^l, ixo^l, 'mhd_to_primitive')
590  end if
591  end subroutine mhd_to_primitive
592 
593  subroutine mhd_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
595  use mod_small_values
596  logical, intent(in) :: primitive
597  integer, intent(in) :: ixI^L,ixO^L
598  double precision, intent(inout) :: w(ixi^s,1:nw)
599  double precision, intent(in) :: x(ixi^s,1:ndim)
600  character(len=*), intent(in) :: subname
601 
602  double precision :: smallone
603  integer :: idir, flag(ixi^s)
604 
605  if (small_values_method == "ignore") return
606 
607  call mhd_check_w(primitive, ixi^l, ixo^l, w, flag)
608 
609  if (any(flag(ixo^s) /= 0)) then
610  select case (small_values_method)
611  case ("replace")
612  if (small_values_fix_iw(rho_)) then
613  where(flag(ixo^s) /= 0) w(ixo^s,rho_) = small_density
614  end if
615 
616  do idir = 1, ndir
617  if (small_values_fix_iw(mom(idir))) then
618  where(flag(ixo^s) /= 0) w(ixo^s, mom(idir)) = 0.0d0
619  end if
620  end do
621 
622  if (mhd_energy) then
623  if (small_values_fix_iw(e_)) then
624  if(primitive) then
625  where(flag(ixo^s) /= 0) w(ixo^s,e_) = small_pressure
626  else
627  where(flag(ixo^s) /= 0)
628  w(ixo^s,e_) = small_e + 0.5d0 * &
629  sum(w(ixo^s, mom(:))**2, dim=ndim+1) / w(ixo^s, rho_) + &
630  mhd_mag_en(w, ixi^l, ixo^l)
631  end where
632  end if
633  end if
634  end if
635  case ("average")
636  call small_values_average(ixi^l, ixo^l, w, x, flag)
637  case default
638  call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
639  end select
640  end if
641  end subroutine mhd_handle_small_values
642 
643  !> Convert energy to entropy
644  subroutine e_to_rhos(ixI^L,ixO^L,w,x)
646  integer, intent(in) :: ixI^L, ixO^L
647  double precision,intent(inout) :: w(ixi^s,nw)
648  double precision, intent(in) :: x(ixi^s,1:ndim)
649 
650  if (mhd_energy) then
651  if(.not.block%e_is_internal) &
652  w(ixo^s, e_) = w(ixo^s, e_) - mhd_kin_en(w, ixi^l, ixo^l) &
653  - mhd_mag_en(w, ixi^l, ixo^l)
654  w(ixo^s, e_) = gamma_1* w(ixo^s, rho_)**(1.0d0 - mhd_gamma) * &
655  w(ixo^s, e_)
656  else
657  call mpistop("e_to_rhos can not be used without energy equation!")
658  end if
659  end subroutine e_to_rhos
660 
661  !> Convert entropy to energy
662  subroutine rhos_to_e(ixI^L,ixO^L,w,x)
664  integer, intent(in) :: ixI^L, ixO^L
665  double precision :: w(ixi^s,nw)
666  double precision, intent(in) :: x(ixi^s,1:ndim)
667 
668  if (mhd_energy) then
669  w(ixo^s, e_) = w(ixo^s, rho_)**gamma_1 * w(ixo^s, e_) &
670  * inv_gamma_1
671  if(.not.block%e_is_internal) &
672  w(ixo^s, e_) =w(ixo^s, e_) + mhd_kin_en(w, ixi^l, ixo^l) + &
673  mhd_mag_en(w, ixi^l, ixo^l)
674  else
675  call mpistop("rhos_to_e can not be used without energy equation!")
676  end if
677  end subroutine rhos_to_e
678 
679  !> Calculate v vector
680  subroutine mhd_get_v(w,x,ixI^L,ixO^L,v)
682 
683  integer, intent(in) :: ixI^L, ixO^L
684  double precision, intent(in) :: w(ixi^s,nw), x(ixi^s,1:ndim)
685  double precision, intent(out) :: v(ixi^s,ndir)
686 
687  integer :: idir
688 
689  do idir=1,ndir
690  v(ixo^s,idir) = w(ixo^s, mom(idir)) / w(ixo^s,rho_)
691  end do
692 
693  end subroutine mhd_get_v
694 
695  !> Calculate v component
696  subroutine mhd_get_v_idim(w,x,ixI^L,ixO^L,idim,v)
698 
699  integer, intent(in) :: ixI^L, ixO^L, idim
700  double precision, intent(in) :: w(ixi^s,nw), x(ixi^s,1:ndim)
701  double precision, intent(out) :: v(ixi^s)
702 
703  v(ixo^s) = w(ixo^s, mom(idim)) / w(ixo^s,rho_)
704 
705  end subroutine mhd_get_v_idim
706 
707  !> Calculate cmax_idim=csound+abs(v_idim) within ixO^L
708  subroutine mhd_get_cmax(w,x,ixI^L,ixO^L,idim,cmax)
710 
711  integer, intent(in) :: ixI^L, ixO^L, idim
712  double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
713  double precision, intent(inout) :: cmax(ixi^s)
714 
715  call mhd_get_csound(w,x,ixi^l,ixo^l,idim,cmax)
716 
717  cmax(ixo^s)=abs(w(ixo^s,mom(idim))/w(ixo^s,rho_))+cmax(ixo^s)
718 
719  end subroutine mhd_get_cmax
720 
721  !> Estimating bounds for the minimum and maximum signal velocities
722  subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,cmax,cmin)
724 
725  integer, intent(in) :: ixI^L, ixO^L, idim
726  double precision, intent(in) :: wLC(ixi^s, nw), wRC(ixi^s, nw)
727  double precision, intent(in) :: wLp(ixi^s, nw), wRp(ixi^s, nw)
728  double precision, intent(in) :: x(ixi^s,1:ndim)
729  double precision, intent(inout) :: cmax(ixi^s)
730  double precision, intent(inout), optional :: cmin(ixi^s)
731 
732  double precision :: wmean(ixi^s,nw)
733  double precision, dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
734 
735  if (typeboundspeed/='cmaxmean') then
736  ! This implements formula (10.52) from "Riemann Solvers and Numerical
737  ! Methods for Fluid Dynamics" by Toro.
738 
739  tmp1(ixo^s)=sqrt(wlp(ixo^s,rho_))
740  tmp2(ixo^s)=sqrt(wrp(ixo^s,rho_))
741  tmp3(ixo^s)=1.d0/(sqrt(wlp(ixo^s,rho_))+sqrt(wrp(ixo^s,rho_)))
742  umean(ixo^s)=(wlp(ixo^s,mom(idim))*tmp1(ixo^s)+wrp(ixo^s,mom(idim))*tmp2(ixo^s))*tmp3(ixo^s)
743  call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
744  call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
745  dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
746  0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
747  (wrp(ixo^s,mom(idim))-wlp(ixo^s,mom(idim)))**2
748  dmean(ixo^s)=sqrt(dmean(ixo^s))
749  if(present(cmin)) then
750  cmin(ixo^s)=umean(ixo^s)-dmean(ixo^s)
751  cmax(ixo^s)=umean(ixo^s)+dmean(ixo^s)
752  else
753  cmax(ixo^s)=abs(umean(ixo^s))+dmean(ixo^s)
754  end if
755  else
756  wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
757  tmp1(ixo^s)=wmean(ixo^s,mom(idim))/wmean(ixo^s,rho_)
758  call mhd_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
759  if(present(cmin)) then
760  cmax(ixo^s)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
761  cmin(ixo^s)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
762  else
763  cmax(ixo^s)=abs(tmp1(ixo^s))+csoundr(ixo^s)
764  end if
765  end if
766 
767  end subroutine mhd_get_cbounds
768 
769  !> Calculate fast magnetosonic wave speed
770  subroutine mhd_get_csound(w,x,ixI^L,ixO^L,idim,csound)
772 
773  integer, intent(in) :: ixI^L, ixO^L, idim
774  double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
775  double precision, intent(out):: csound(ixi^s)
776  double precision :: cfast2(ixi^s), AvMinCs2(ixi^s), b2(ixi^s), kmax
777  double precision :: inv_rho(ixo^s)
778 
779  inv_rho=1.d0/w(ixo^s,rho_)
780 
781  call mhd_get_csound2(w,x,ixi^l,ixo^l,csound)
782  ! store |B|^2 in v
783  b2(ixo^s) = mhd_mag_en_all(w,ixi^l,ixo^l)
784  cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
785  avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
786  * mhd_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
787  * inv_rho
788 
789  where(avmincs2(ixo^s)<zero)
790  avmincs2(ixo^s)=zero
791  end where
792 
793  avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
794 
795  if (.not. mhd_hall) then
796  csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
797  else
798  ! take the Hall velocity into account:
799  ! most simple estimate, high k limit:
800  ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
801  kmax = dpi/min({dxlevel(^d)},bigdouble)*half
802  csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
803  mhd_etah * sqrt(b2(ixo^s))*inv_rho*kmax)
804  end if
805 
806  end subroutine mhd_get_csound
807 
808  !> Calculate fast magnetosonic wave speed
809  subroutine mhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
811 
812  integer, intent(in) :: ixI^L, ixO^L, idim
813  double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
814  double precision, intent(out):: csound(ixi^s)
815  double precision :: cfast2(ixi^s), AvMinCs2(ixi^s), b2(ixi^s), kmax
816  double precision :: inv_rho(ixo^s)
817 
818  inv_rho=1.d0/w(ixo^s,rho_)
819 
820  if(mhd_energy) then
821  csound(ixo^s)=mhd_gamma*w(ixo^s,p_)*inv_rho
822  else
823  csound(ixo^s)=mhd_gamma*mhd_adiab*w(ixo^s,rho_)**gamma_1
824  end if
825  ! store |B|^2 in v
826  b2(ixo^s) = mhd_mag_en_all(w,ixi^l,ixo^l)
827  cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
828  avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
829  * mhd_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
830  * inv_rho
831 
832  where(avmincs2(ixo^s)<zero)
833  avmincs2(ixo^s)=zero
834  end where
835 
836  avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
837 
838  if (.not. mhd_hall) then
839  csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
840  else
841  ! take the Hall velocity into account:
842  ! most simple estimate, high k limit:
843  ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
844  kmax = dpi/min({dxlevel(^d)},bigdouble)*half
845  csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
846  mhd_etah * sqrt(b2(ixo^s))*inv_rho*kmax)
847  end if
848 
849  end subroutine mhd_get_csound_prim
850 
851  !> Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho-b**2/2) within ixO^L
852  subroutine mhd_get_pthermal(w,x,ixI^L,ixO^L,pth)
854 
855  integer, intent(in) :: ixI^L, ixO^L
856  double precision, intent(in) :: w(ixi^s,nw)
857  double precision, intent(in) :: x(ixi^s,1:ndim)
858  double precision, intent(out):: pth(ixi^s)
859 
860  if(mhd_energy) then
861  if(block%e_is_internal) then
862  pth(ixo^s)=gamma_1*w(ixo^s,e_)
863  else
864  pth(ixo^s)=gamma_1*(w(ixo^s,e_)&
865  - mhd_kin_en(w,ixi^l,ixo^l)&
866  - mhd_mag_en(w,ixi^l,ixo^l))
867  end if
868  else
869  pth(ixo^s)=mhd_adiab*w(ixo^s,rho_)**mhd_gamma
870  end if
871  end subroutine mhd_get_pthermal
872 
873  !> Calculate the square of the thermal sound speed csound2 within ixO^L.
874  !> csound2=gamma*p/rho
875  subroutine mhd_get_csound2(w,x,ixI^L,ixO^L,csound2)
877  integer, intent(in) :: ixI^L, ixO^L
878  double precision, intent(in) :: w(ixi^s,nw)
879  double precision, intent(in) :: x(ixi^s,1:ndim)
880  double precision, intent(out) :: csound2(ixi^s)
881 
882  if(mhd_energy) then
883  call mhd_get_pthermal(w,x,ixi^l,ixo^l,csound2)
884  csound2(ixo^s)=mhd_gamma*csound2(ixo^s)/w(ixo^s,rho_)
885  else
886  csound2(ixo^s)=mhd_gamma*mhd_adiab*w(ixo^s,rho_)**gamma_1
887  end if
888  end subroutine mhd_get_csound2
889 
890  !> Calculate total pressure within ixO^L including magnetic pressure
891  subroutine mhd_get_p_total(w,x,ixI^L,ixO^L,p)
893 
894  integer, intent(in) :: ixI^L, ixO^L
895  double precision, intent(in) :: w(ixi^s,nw)
896  double precision, intent(in) :: x(ixi^s,1:ndim)
897  double precision, intent(out) :: p(ixi^s)
898 
899  call mhd_get_pthermal(w,x,ixi^l,ixo^l,p)
900 
901  p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
902 
903  end subroutine mhd_get_p_total
904 
905  !> Calculate fluxes within ixO^L.
906  subroutine mhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
908  use mod_usr_methods
909 
910  integer, intent(in) :: ixI^L, ixO^L, idim
911  ! conservative w
912  double precision, intent(in) :: wC(ixi^s,nw)
913  ! primitive w
914  double precision, intent(in) :: w(ixi^s,nw)
915  double precision, intent(in) :: x(ixi^s,1:ndim)
916  double precision,intent(out) :: f(ixi^s,nwflux)
917 
918  double precision :: ptotal(ixo^s),tmp(ixi^s)
919  double precision, allocatable:: vHall(:^d&,:)
920  integer :: idirmin, iw, idir
921 
922  if (mhd_hall) then
923  allocate(vhall(ixi^s,1:ndir))
924  call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
925  end if
926 
927  if(b0field) tmp(ixo^s)=sum(block%B0(ixo^s,:,idim)*w(ixo^s,mag(:)),dim=ndim+1)
928 
929  if(mhd_energy) then
930  ptotal=w(ixo^s,p_)+0.5d0*sum(w(ixo^s, mag(:))**2, dim=ndim+1)
931  else
932  ptotal(ixo^s)=mhd_adiab*w(ixo^s,rho_)**mhd_gamma+0.5d0*sum(w(ixo^s, mag(:))**2, dim=ndim+1)
933  end if
934 
935  ! Get flux of density
936  f(ixo^s,rho_)=w(ixo^s,mom(idim))*w(ixo^s,rho_)
937 
938  ! Get flux of tracer
939  do iw=1,mhd_n_tracer
940  f(ixo^s,tracer(iw))=w(ixo^s,mom(idim))*w(ixo^s,tracer(iw))
941  end do
942 
943  ! Get flux of momentum
944  ! f_i[m_k]=v_i*m_k-b_k*b_i [+ptotal if i==k]
945  do idir=1,ndir
946  if(idim==idir) then
947  f(ixo^s,mom(idir))=ptotal(ixo^s)-w(ixo^s,mag(idim))*w(ixo^s,mag(idir))
948  if(b0field) f(ixo^s,mom(idir))=f(ixo^s,mom(idir))+tmp(ixo^s)
949  else
950  f(ixo^s,mom(idir))= -w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
951  end if
952  if (b0field) then
953  f(ixo^s,mom(idir))=f(ixo^s,mom(idir))&
954  -w(ixo^s,mag(idir))*block%B0(ixo^s,idim,idim)&
955  -w(ixo^s,mag(idim))*block%B0(ixo^s,idir,idim)
956  end if
957  f(ixo^s,mom(idir))=f(ixo^s,mom(idir))+w(ixo^s,mom(idim))*wc(ixo^s,mom(idir))
958  end do
959 
960  ! Get flux of energy
961  ! f_i[e]=v_i*e+v_i*ptotal-b_i*(b_k*v_k)
962  if (mhd_energy) then
963  if (block%e_is_internal) then
964  f(ixo^s,e_)=w(ixo^s,mom(idim))*w(ixo^s,p_)
965  if (mhd_hall) then
966  call mpistop("solve pthermal not designed for Hall MHD")
967  endif
968  else
969  f(ixo^s,e_)=w(ixo^s,mom(idim))*(wc(ixo^s,e_) + ptotal(ixo^s))- &
970  w(ixo^s,mag(idim))*sum(w(ixo^s,mag(:))*w(ixo^s,mom(:)),dim=ndim+1)
971 
972  if(type_divb==divb_glm2) then
973  f(ixo^s,e_) = f(ixo^s,e_) + vmax_global*w(ixo^s,psi_)*w(ixo^s,mag(idim))
974  end if
975 
976  if (b0field) then
977  f(ixo^s,e_) = f(ixo^s,e_) &
978  + w(ixo^s,mom(idim)) * tmp(ixo^s) &
979  - sum(w(ixo^s,mom(:))*w(ixo^s,mag(:)),dim=ndim+1) * block%B0(ixo^s,idim,idim)
980  end if
981 
982  if (mhd_hall) then
983  ! f_i[e]= f_i[e] + vHall_i*(b_k*b_k) - b_i*(vHall_k*b_k)
984  if (mhd_etah>zero) then
985  f(ixo^s,e_) = f(ixo^s,e_) + vhall(ixo^s,idim) * &
986  sum(w(ixo^s, mag(:))**2,dim=ndim+1) &
987  - w(ixo^s,mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1)
988  if (b0field) then
989  f(ixo^s,e_) = f(ixo^s,e_) &
990  + vhall(ixo^s,idim) * tmp(ixo^s) &
991  - sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1) * block%B0(ixo^s,idim,idim)
992  end if
993  end if
994  end if
995  end if
996  end if
997 
998  ! compute flux of magnetic field
999  ! f_i[b_k]=v_i*b_k-v_k*b_i
1000  do idir=1,ndir
1001  if (idim==idir) then
1002  ! f_i[b_i] should be exactly 0, so we do not use the transport flux
1003  if (mhd_glm) then
1004  if(type_divb==divb_glm1) then
1005  f(ixo^s,mag(idir))=w(ixo^s,psi_)
1006  else
1007  f(ixo^s,mag(idir))=vmax_global*w(ixo^s,psi_)
1008  end if
1009  else
1010  f(ixo^s,mag(idir))=zero
1011  end if
1012  else
1013  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))
1014 
1015  if (b0field) then
1016  f(ixo^s,mag(idir))=f(ixo^s,mag(idir))&
1017  +w(ixo^s,mom(idim))*block%B0(ixo^s,idir,idim)&
1018  -w(ixo^s,mom(idir))*block%B0(ixo^s,idim,idim)
1019  end if
1020 
1021  if (mhd_hall) then
1022  ! f_i[b_k] = f_i[b_k] + vHall_i*b_k - vHall_k*b_i
1023  if (mhd_etah>zero) then
1024  if (b0field) then
1025  f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
1026  - vhall(ixo^s,idir)*(w(ixo^s,mag(idim))+block%B0(ixo^s,idim,idim)) &
1027  + vhall(ixo^s,idim)*(w(ixo^s,mag(idir))+block%B0(ixo^s,idir,idim))
1028  else
1029  f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
1030  - vhall(ixo^s,idir)*w(ixo^s,mag(idim)) &
1031  + vhall(ixo^s,idim)*w(ixo^s,mag(idir))
1032  end if
1033  end if
1034  end if
1035 
1036  end if
1037  end do
1038 
1039  if (mhd_glm) then
1040  if(type_divb==divb_glm1) then
1041  !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645
1042  f(ixo^s,psi_) = cmax_global**2*w(ixo^s,mag(idim))
1043  else
1044  !f_i[psi]=Ch*b_{i} Eq. 3.16e Derigs et al 2018 JCP, 364, 420
1045  f(ixo^s,psi_) = vmax_global*w(ixo^s,mag(idim))
1046  end if
1047  end if
1048 
1049  end subroutine mhd_get_flux
1050 
1051  !> w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
1052  subroutine mhd_add_source(qdt,ixI^L,ixO^L,wCT,w,x,qsourcesplit,active)
1056  use mod_gravity, only: gravity_add_source
1057 
1058  integer, intent(in) :: ixI^L, ixO^L
1059  double precision, intent(in) :: qdt
1060  double precision, intent(in) :: wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1061  double precision, intent(inout) :: w(ixi^s,1:nw)
1062  logical, intent(in) :: qsourcesplit
1063  logical, intent(inout) :: active
1064 
1065  if (.not. qsourcesplit) then
1066  ! Source for solving internal energy
1067  if (mhd_energy .and. block%e_is_internal) then
1068  active = .true.
1069  call internal_energy_add_source(qdt,ixi^l,ixo^l,wct,w,x)
1070  endif
1071 
1072  ! Source for B0 splitting
1073  if (b0field) then
1074  active = .true.
1075  call add_source_b0split(qdt,ixi^l,ixo^l,wct,w,x)
1076  end if
1077 
1078  ! Sources for resistivity in eqs. for e, B1, B2 and B3
1079  if (abs(mhd_eta)>smalldouble)then
1080  active = .true.
1081  call add_source_res2(qdt,ixi^l,ixo^l,wct,w,x)
1082  end if
1083 
1084  if (mhd_eta_hyper>0.d0)then
1085  active = .true.
1086  call add_source_hyperres(qdt,ixi^l,ixo^l,wct,w,x)
1087  end if
1088  end if
1089 
1090  {^nooned
1091  if(.not.source_split_divb .and. .not.qsourcesplit .and. istep==nstep) then
1092  ! Sources related to div B
1093  select case (type_divb)
1094  case (divb_none)
1095  ! Do nothing
1096  case (divb_glm1)
1097  active = .true.
1098  call add_source_glm1(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1099  case (divb_glm2)
1100  active = .true.
1101  call add_source_glm2(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1102  case (divb_powel)
1103  active = .true.
1104  call add_source_powel(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1105  case (divb_janhunen)
1106  active = .true.
1107  call add_source_janhunen(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1108  case (divb_linde)
1109  active = .true.
1110  call add_source_linde(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1111  case (divb_lindejanhunen)
1112  active = .true.
1113  call add_source_linde(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1114  call add_source_janhunen(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1115  case (divb_lindepowel)
1116  active = .true.
1117  call add_source_linde(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1118  call add_source_powel(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1119  case (divb_lindeglm)
1120  active = .true.
1121  call add_source_linde(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1122  call add_source_glm2(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1123  case (divb_multigrid)
1124  continue ! Do nothing
1125  case default
1126  call mpistop('Unknown divB fix')
1127  end select
1128  else if(source_split_divb .and. qsourcesplit) then
1129  ! Sources related to div B
1130  select case (type_divb)
1131  case (divb_none)
1132  ! Do nothing
1133  case (divb_glm1)
1134  active = .true.
1135  call add_source_glm1(qdt,ixi^l,ixo^l,wct,w,x)
1136  case (divb_glm2)
1137  active = .true.
1138  call add_source_glm2(dt,ixi^l,ixo^l,pso(saveigrid)%w,w,x)
1139  case (divb_powel)
1140  active = .true.
1141  call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
1142  case (divb_janhunen)
1143  active = .true.
1144  call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
1145  case (divb_linde)
1146  active = .true.
1147  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
1148  case (divb_lindejanhunen)
1149  active = .true.
1150  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
1151  call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
1152  case (divb_lindepowel)
1153  active = .true.
1154  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
1155  call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
1156  case (divb_lindeglm)
1157  active = .true.
1158  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
1159  call add_source_glm2(qdt,ixi^l,ixo^l,wct,w,x)
1160  case (divb_multigrid)
1161  continue ! Do nothing
1162  case default
1163  call mpistop('Unknown divB fix')
1164  end select
1165  end if
1166  }
1167 
1168  if(mhd_radiative_cooling) then
1169  call radiative_cooling_add_source(qdt,ixi^l,ixo^l,wct,&
1170  w,x,qsourcesplit,active)
1171  end if
1172 
1173  if(mhd_viscosity) then
1174  call viscosity_add_source(qdt,ixi^l,ixo^l,wct,&
1175  w,x,mhd_energy,qsourcesplit,active)
1176  end if
1177 
1178  if(mhd_gravity) then
1179  call gravity_add_source(qdt,ixi^l,ixo^l,wct,&
1180  w,x,mhd_energy,qsourcesplit,active)
1181  end if
1182 
1183  end subroutine mhd_add_source
1184 
1185  subroutine internal_energy_add_source(qdt,ixI^L,ixO^L,wCT,w,x)
1187  use mod_geometry
1188 
1189  integer, intent(in) :: ixI^L, ixO^L
1190  double precision, intent(in) :: qdt
1191  double precision, intent(in) :: wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1192  double precision, intent(inout) :: w(ixi^s,1:nw)
1193  double precision :: pth(ixi^s),v(ixi^s,1:ndir),divv(ixi^s)
1194 
1195  call mhd_get_v(wct,x,ixi^l,ixi^l,v)
1196  call divvector(v,ixi^l,ixo^l,divv)
1197  call mhd_get_pthermal(wct,x,ixi^l,ixo^l,pth)
1198  w(ixo^s,e_)=w(ixo^s,e_)-qdt*pth(ixo^s)*divv(ixo^s)
1199 
1200  end subroutine internal_energy_add_source
1201 
1202  !> Source terms after split off time-independent magnetic field
1203  subroutine add_source_b0split(qdt,ixI^L,ixO^L,wCT,w,x)
1205 
1206  integer, intent(in) :: ixI^L, ixO^L
1207  double precision, intent(in) :: qdt, wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1208  double precision, intent(inout) :: w(ixi^s,1:nw)
1209 
1210  double precision :: a(ixi^s,3), b(ixi^s,3), axb(ixi^s,3)
1211  integer :: idir
1212 
1213  a=0.d0
1214  b=0.d0
1215  ! for force-free field J0xB0 =0
1216  if(.not.b0field_forcefree) then
1217  ! store B0 magnetic field in b
1218  b(ixo^s,1:ndir)=block%B0(ixo^s,1:ndir,0)
1219 
1220  ! store J0 current in a
1221  do idir=7-2*ndir,3
1222  a(ixo^s,idir)=block%J0(ixo^s,idir)
1223  end do
1224  call cross_product(ixi^l,ixo^l,a,b,axb)
1225  axb(ixo^s,:)=axb(ixo^s,:)*qdt
1226  ! add J0xB0 source term in momentum equations
1227  w(ixo^s,mom(1:ndir))=w(ixo^s,mom(1:ndir))+axb(ixo^s,1:ndir)
1228  end if
1229 
1230  if(mhd_energy) then
1231  if(.not.block%e_is_internal) then
1232  a=0.d0
1233  ! for free-free field -(vxB0) dot J0 =0
1234  b(ixo^s,:)=wct(ixo^s,mag(:))
1235  ! store full magnetic field B0+B1 in b
1236  if(.not.b0field_forcefree) b(ixo^s,:)=b(ixo^s,:)+block%B0(ixo^s,:,0)
1237  ! store velocity in a
1238  do idir=1,ndir
1239  a(ixo^s,idir)=wct(ixo^s,mom(idir))/wct(ixo^s,rho_)
1240  end do
1241  call cross_product(ixi^l,ixo^l,a,b,axb)
1242  axb(ixo^s,:)=axb(ixo^s,:)*qdt
1243  ! add -(vxB) dot J0 source term in energy equation
1244  do idir=7-2*ndir,3
1245  w(ixo^s,e_)=w(ixo^s,e_)-axb(ixo^s,idir)*block%J0(ixo^s,idir)
1246  end do
1247  end if
1248  end if
1249 
1250  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_B0')
1251 
1252  end subroutine add_source_b0split
1253 
1254  !> Add resistive source to w within ixO Uses 3 point stencil (1 neighbour) in
1255  !> each direction, non-conservative. If the fourthorder precompiler flag is
1256  !> set, uses fourth order central difference for the laplacian. Then the
1257  !> stencil is 5 (2 neighbours).
1258  subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
1260  use mod_usr_methods
1261  use mod_geometry
1262 
1263  integer, intent(in) :: ixI^L, ixO^L
1264  double precision, intent(in) :: qdt
1265  double precision, intent(in) :: wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1266  double precision, intent(inout) :: w(ixi^s,1:nw)
1267  integer :: ixA^L,idir,jdir,kdir,idirmin,idim,jxO^L,hxO^L,ix
1268  integer :: lxO^L, kxO^L
1269 
1270  double precision :: tmp(ixi^s),tmp2(ixi^s)
1271 
1272  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
1273  double precision :: current(ixi^s,7-2*ndir:3),eta(ixi^s)
1274  double precision :: gradeta(ixi^s,1:ndim), Bf(ixi^s,1:ndir)
1275 
1276  ! Calculating resistive sources involve one extra layer
1277  if (mhd_4th_order) then
1278  ixa^l=ixo^l^ladd2;
1279  else
1280  ixa^l=ixo^l^ladd1;
1281  end if
1282 
1283  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
1284  call mpistop("Error in add_source_res1: Non-conforming input limits")
1285 
1286  ! Calculate current density and idirmin
1287  call get_current(wct,ixi^l,ixo^l,idirmin,current)
1288 
1289  if (mhd_eta>zero)then
1290  eta(ixa^s)=mhd_eta
1291  gradeta(ixo^s,1:ndim)=zero
1292  else
1293  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
1294  ! assumes that eta is not function of current?
1295  do idim=1,ndim
1296  call gradient(eta,ixi^l,ixo^l,idim,tmp)
1297  gradeta(ixo^s,idim)=tmp(ixo^s)
1298  end do
1299  end if
1300 
1301  if(b0field) then
1302  bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))+block%B0(ixi^s,1:ndir,0)
1303  else
1304  bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
1305  end if
1306 
1307  do idir=1,ndir
1308  ! Put B_idir into tmp2 and eta*Laplace B_idir into tmp
1309  if (mhd_4th_order) then
1310  tmp(ixo^s)=zero
1311  tmp2(ixi^s)=bf(ixi^s,idir)
1312  do idim=1,ndim
1313  lxo^l=ixo^l+2*kr(idim,^d);
1314  jxo^l=ixo^l+kr(idim,^d);
1315  hxo^l=ixo^l-kr(idim,^d);
1316  kxo^l=ixo^l-2*kr(idim,^d);
1317  tmp(ixo^s)=tmp(ixo^s)+&
1318  (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
1319  /(12.0d0 * dxlevel(idim)**2)
1320  end do
1321  else
1322  tmp(ixo^s)=zero
1323  tmp2(ixi^s)=bf(ixi^s,idir)
1324  do idim=1,ndim
1325  jxo^l=ixo^l+kr(idim,^d);
1326  hxo^l=ixo^l-kr(idim,^d);
1327  tmp(ixo^s)=tmp(ixo^s)+&
1328  (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/dxlevel(idim)**2
1329  end do
1330  end if
1331 
1332  ! Multiply by eta
1333  tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
1334 
1335  ! Subtract grad(eta) x J = eps_ijk d_j eta J_k if eta is non-constant
1336  if (mhd_eta<zero)then
1337  do jdir=1,ndim; do kdir=idirmin,3
1338  if (lvc(idir,jdir,kdir)/=0)then
1339  if (lvc(idir,jdir,kdir)==1)then
1340  tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
1341  else
1342  tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
1343  end if
1344  end if
1345  end do; end do
1346  end if
1347 
1348  ! Add sources related to eta*laplB-grad(eta) x J to B and e
1349  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
1350  if (mhd_energy) then
1351  w(ixo^s,e_)=w(ixo^s,e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
1352  end if
1353 
1354  end do ! idir
1355 
1356  if (mhd_energy) then
1357  ! de/dt+=eta*J**2
1358  tmp(ixo^s)=zero
1359  do idir=idirmin,3
1360  tmp(ixo^s)=tmp(ixo^s)+current(ixo^s,idir)**2
1361  end do
1362  w(ixo^s,e_)=w(ixo^s,e_)+qdt*eta(ixo^s)*tmp(ixo^s)
1363  end if
1364 
1365  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_res1')
1366 
1367  end subroutine add_source_res1
1368 
1369  !> Add resistive source to w within ixO
1370  !> Uses 5 point stencil (2 neighbours) in each direction, conservative
1371  subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
1373  use mod_usr_methods
1374  use mod_geometry
1375 
1376  integer, intent(in) :: ixI^L, ixO^L
1377  double precision, intent(in) :: qdt
1378  double precision, intent(in) :: wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1379  double precision, intent(inout) :: w(ixi^s,1:nw)
1380  integer :: ixA^L,idir,jdir,kdir,idirmin,iw,idim,idirmin1
1381 
1382  double precision :: tmp(ixi^s),tmp2(ixi^s)
1383 
1384  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
1385  double precision :: current(ixi^s,7-2*ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
1386  double precision :: tmpvec(ixi^s,1:3)
1387 
1388  ixa^l=ixo^l^ladd2;
1389 
1390  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
1391  call mpistop("Error in add_source_res2: Non-conforming input limits")
1392 
1393  ixa^l=ixo^l^ladd1;
1394  ! Calculate current density within ixL: J=curl B, thus J_i=eps_ijk*d_j B_k
1395  ! Determine exact value of idirmin while doing the loop.
1396  call get_current(wct,ixi^l,ixa^l,idirmin,current)
1397 
1398  if (mhd_eta>zero)then
1399  eta(ixa^s)=mhd_eta
1400  else
1401  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
1402  end if
1403 
1404  ! dB/dt= -curl(J*eta), thus B_i=B_i-eps_ijk d_j Jeta_k
1405  tmpvec(ixa^s,1:ndir)=zero
1406  do jdir=idirmin,3
1407  tmpvec(ixa^s,jdir)=current(ixa^s,jdir)*eta(ixa^s)
1408  end do
1409  call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
1410  do idir=1,ndir
1411  w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-qdt*curlj(ixo^s,idir)
1412  end do
1413 
1414  if(mhd_energy) then
1415  ! de/dt= +div(B x Jeta) = eta J^2 - B dot curl(eta J)
1416  ! de1/dt= eta J^2 - B1 dot curl(eta J)
1417  w(ixo^s,e_)=w(ixo^s,e_)+qdt*(sum(current(ixo^s,:)**2,dim=ndim+1)*eta(ixo^s)-&
1418  sum(wct(ixo^s,mag(1:ndir))*curlj(ixo^s,1:ndir),dim=ndim+1))
1419  end if
1420 
1421  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_res2')
1422 
1423  end subroutine add_source_res2
1424 
1425  !> Add Hyper-resistive source to w within ixO
1426  !> Uses 9 point stencil (4 neighbours) in each direction.
1427  subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
1429  use mod_geometry
1430  use mod_geometry
1431 
1432  integer, intent(in) :: ixI^L, ixO^L
1433  double precision, intent(in) :: qdt
1434  double precision, intent(in) :: wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1435  double precision, intent(inout) :: w(ixi^s,1:nw)
1436  !.. local ..
1437  double precision :: current(ixi^s,7-2*ndir:3)
1438  double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
1439  integer :: ixA^L,idir,jdir,kdir,idirmin,idirmin1
1440 
1441  ixa^l=ixo^l^ladd3;
1442  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
1443  call mpistop("Error in add_source_hyperres: Non-conforming input limits")
1444 
1445  call get_current(wct,ixi^l,ixa^l,idirmin,current)
1446  tmpvec(ixa^s,1:ndir)=zero
1447  do jdir=idirmin,3
1448  tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
1449  end do
1450 
1451  ixa^l=ixo^l^ladd2;
1452  call curlvector(tmpvec,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
1453 
1454  ixa^l=ixo^l^ladd1;
1455  tmpvec(ixa^s,1:ndir)=zero
1456  call curlvector(tmpvec2,ixi^l,ixa^l,tmpvec,idirmin1,1,3)
1457  ehyper(ixa^s,1:ndir) = - tmpvec(ixa^s,1:ndir)*mhd_eta_hyper
1458 
1459  ixa^l=ixo^l;
1460  tmpvec2(ixa^s,1:ndir)=zero
1461  call curlvector(ehyper,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
1462 
1463  do idir=1,ndir
1464  w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
1465  end do
1466 
1467  if (mhd_energy) then
1468  ! de/dt= +div(B x Ehyper)
1469  ixa^l=ixo^l^ladd1;
1470  tmpvec2(ixa^s,1:ndir)=zero
1471  do idir=1,ndir; do jdir=1,ndir; do kdir=idirmin,3
1472  tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
1473  + lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
1474  end do; end do; end do
1475  tmp(ixo^s)=zero
1476  call divvector(tmpvec2,ixi^l,ixo^l,tmp)
1477  w(ixo^s,e_)=w(ixo^s,e_)+tmp(ixo^s)*qdt
1478  end if
1479 
1480  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_hyperres')
1481 
1482  end subroutine add_source_hyperres
1483 
1484  subroutine add_source_glm1(qdt,ixI^L,ixO^L,wCT,w,x)
1485  ! Add divB related sources to w within ixO
1486  ! corresponding to Dedner JCP 2002, 175, 645 _equation 24_
1487  ! giving the EGLM-MHD scheme
1489  use mod_geometry
1490 
1491  integer, intent(in) :: ixI^L, ixO^L
1492  double precision, intent(in) :: qdt, wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1493  double precision, intent(inout) :: w(ixi^s,1:nw)
1494  double precision:: divb(ixi^s)
1495  integer :: idim,idir
1496  double precision :: gradPsi(ixi^s)
1497 
1498  ! We calculate now div B
1499  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
1500 
1501  ! dPsi/dt = - Ch^2/Cp^2 Psi
1502  if (mhd_glm_alpha < zero) then
1503  w(ixo^s,psi_) = abs(mhd_glm_alpha)*wct(ixo^s,psi_)
1504  else
1505  ! implicit update of Psi variable
1506  ! equation (27) in Mignone 2010 J. Com. Phys. 229, 2117
1507  if(slab) then
1508  w(ixo^s,psi_) = dexp(-qdt*cmax_global*mhd_glm_alpha/minval(dxlevel(:)))*w(ixo^s,psi_)
1509  else
1510  w(ixo^s,psi_) = dexp(-qdt*cmax_global*mhd_glm_alpha/minval(block%ds(ixo^s,:),dim=ndim+1))*w(ixo^s,psi_)
1511  end if
1512  end if
1513 
1514  ! gradient of Psi
1515  do idim=1,ndim
1516  select case(typegrad)
1517  case("central")
1518  call gradient(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
1519  case("limited")
1520  call gradients(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
1521  end select
1522  if (mhd_energy .and. .not.block%e_is_internal) then
1523  ! e = e -qdt (b . grad(Psi))
1524  w(ixo^s,e_) = w(ixo^s,e_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
1525  end if
1526  end do
1527 
1528  ! m = m - qdt b div b
1529  do idir=1,ndir
1530  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-qdt*mhd_mag_i_all(w,ixi^l,ixo^l,idir)*divb(ixo^s)
1531  end do
1532 
1533  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_glm1')
1534 
1535  end subroutine add_source_glm1
1536 
1537  subroutine add_source_glm2(qdt,ixI^L,ixO^L,wCT,w,x)
1538  ! Add divB related sources to w within ixO
1539  ! corresponding to Eq. 3.17 Derigs et al 2018 JCP, 364, 420
1541  use mod_geometry
1542 
1543  integer, intent(in) :: ixI^L, ixO^L
1544  double precision, intent(in) :: qdt, wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1545  double precision, intent(inout) :: w(ixi^s,1:nw)
1546  double precision:: divb(ixi^s),v(ixi^s,1:ndir)
1547  integer :: idim,idir
1548  double precision :: gradPsi(ixi^s,ndim), Bf(ixi^s,1:ndir)
1549 
1550  ! calculate now div B
1551  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
1552 
1553  ! calculate velocity
1554  call mhd_get_v(wct,x,ixi^l,ixo^l,v)
1555 
1556  ! gradient of Psi
1557  do idim=1,ndim
1558  select case(typegrad)
1559  case("central")
1560  call gradient(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi(ixi^s,idim))
1561  case("limited")
1562  call gradients(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi(ixi^s,idim))
1563  end select
1564  end do
1565 
1566  if(b0field) then
1567  bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))+block%B0(ixi^s,1:ndir,0)
1568  else
1569  bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
1570  end if
1571 
1572  if (mhd_energy .and. .not.block%e_is_internal) then
1573  ! e = e - qdt ( (v . b) * div b + (grad psi . v) * psi)
1574  w(ixo^s,e_)=w(ixo^s,e_) - qdt * (divb(ixo^s) * &
1575  sum(v(ixo^s,:)*bf(ixo^s,:),dim=ndim+1) + wct(ixo^s,psi_)*&
1576  sum(v(ixo^s,1:ndim)*gradpsi(ixo^s,1:ndim),dim=ndim+1))
1577  end if
1578 
1579  ! b_i = b_i - qdt * v_i * div b
1580  do idir=1,ndir
1581  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
1582  end do
1583 
1584  ! m_i = m_i - qdt * b_i * div b
1585  do idir=1,ndir
1586  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-qdt*bf(ixo^s,idir)*divb(ixo^s)
1587  end do
1588 
1589  ! psi = psi - qdt * (v . grad(psi) + alpha * psi)
1590  w(ixo^s,psi_) = w(ixo^s,psi_)-qdt*(sum(v(ixo^s,1:ndim)*gradpsi(ixo^s,1:ndim),dim=ndim+1)+&
1591  vmax_global/0.18d0*wct(ixo^s,psi_))
1592 
1593  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_glm')
1594 
1595  end subroutine add_source_glm2
1596 
1597  !> Add divB related sources to w within ixO corresponding to Powel
1598  subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
1600 
1601  integer, intent(in) :: ixI^L, ixO^L
1602  double precision, intent(in) :: qdt, wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1603  double precision, intent(inout) :: w(ixi^s,1:nw)
1604  double precision :: divb(ixi^s),v(ixi^s,1:ndir)
1605  integer :: idir
1606 
1607  ! We calculate now div B
1608  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
1609 
1610  ! calculate velocity
1611  call mhd_get_v(wct,x,ixi^l,ixo^l,v)
1612 
1613  if (mhd_energy .and. .not.block%e_is_internal) then
1614  ! e = e - qdt (v . b) * div b
1615  w(ixo^s,e_)=w(ixo^s,e_)-&
1616  qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=ndim+1)*divb(ixo^s)
1617  end if
1618 
1619  ! b = b - qdt v * div b
1620  do idir=1,ndir
1621  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
1622  end do
1623 
1624  ! m = m - qdt b div b
1625  do idir=1,ndir
1626  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-qdt*mhd_mag_i_all(w,ixi^l,ixo^l,idir)*divb(ixo^s)
1627  end do
1628 
1629  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_powel')
1630 
1631  end subroutine add_source_powel
1632 
1633  subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
1634  ! Add divB related sources to w within ixO
1635  ! corresponding to Janhunen, just the term in the induction equation.
1637 
1638  integer, intent(in) :: ixI^L, ixO^L
1639  double precision, intent(in) :: qdt, wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1640  double precision, intent(inout) :: w(ixi^s,1:nw)
1641  double precision :: divb(ixi^s)
1642  integer :: idir
1643 
1644  ! We calculate now div B
1645  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
1646 
1647  ! b = b - qdt v * div b
1648  do idir=1,ndir
1649  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*wct(ixo^s,mom(idir))/wct(ixo^s,rho_)*divb(ixo^s)
1650  end do
1651 
1652  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_janhunen')
1653 
1654  end subroutine add_source_janhunen
1655 
1656  subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
1657  ! Add Linde's divB related sources to wnew within ixO
1659  use mod_geometry
1660 
1661  integer, intent(in) :: ixI^L, ixO^L
1662  double precision, intent(in) :: qdt, wCT(ixi^s,1:nw), x(ixi^s,1:ndim)
1663  double precision, intent(inout) :: w(ixi^s,1:nw)
1664  integer :: idim, idir, ixp^L, i^D, iside
1665  double precision :: divb(ixi^s),graddivb(ixi^s)
1666  logical, dimension(-1:1^D&) :: leveljump
1667 
1668  ! Calculate div B
1669  ixp^l=ixo^l^ladd1;
1670  call get_divb(wct,ixi^l,ixp^l,divb, mhd_divb_4thorder)
1671 
1672  ! for AMR stability, retreat one cell layer from the boarders of level jump
1673  {do i^db=-1,1\}
1674  if(i^d==0|.and.) cycle
1675  if(neighbor_type(i^d,saveigrid)==2 .or. neighbor_type(i^d,saveigrid)==4) then
1676  leveljump(i^d)=.true.
1677  else
1678  leveljump(i^d)=.false.
1679  end if
1680  {end do\}
1681 
1682  ixp^l=ixo^l;
1683  do idim=1,ndim
1684  select case(idim)
1685  {case(^d)
1686  do iside=1,2
1687  i^dd=kr(^dd,^d)*(2*iside-3);
1688  if (leveljump(i^dd)) then
1689  if (iside==1) then
1690  ixpmin^d=ixomin^d-i^d
1691  else
1692  ixpmax^d=ixomax^d-i^d
1693  end if
1694  end if
1695  end do
1696  \}
1697  end select
1698  end do
1699 
1700  ! Add Linde's diffusive terms
1701  do idim=1,ndim
1702  ! Calculate grad_idim(divb)
1703  select case(typegrad)
1704  case("central")
1705  call gradient(divb,ixi^l,ixp^l,idim,graddivb)
1706  case("limited")
1707  call gradients(divb,ixi^l,ixp^l,idim,graddivb)
1708  end select
1709 
1710  ! Multiply by Linde's eta*dt = divbdiff*(c_max*dx)*dt = divbdiff*dx**2
1711  if (slab) then
1712  graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
1713  else
1714  graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
1715  /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
1716  end if
1717 
1718  w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
1719 
1720  if (mhd_energy .and. typedivbdiff=='all' .and. .not.block%e_is_internal) then
1721  ! e += B_idim*eta*grad_idim(divb)
1722  w(ixp^s,e_)=w(ixp^s,e_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
1723  end if
1724  end do
1725 
1726  if (check_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_linde')
1727 
1728  end subroutine add_source_linde
1729 
1730  !> Calculate div B within ixO
1731  subroutine get_divb(w,ixI^L,ixO^L,divb, fourthorder)
1734  use mod_geometry
1735 
1736  integer, intent(in) :: ixI^L, ixO^L
1737  double precision, intent(in) :: w(ixi^s,1:nw)
1738  double precision, intent(inout) :: divb(ixi^s)
1739  logical, intent(in), optional :: fourthorder
1740 
1741  double precision :: bvec(ixi^s,1:ndir)
1742 
1743  bvec(ixi^s,:)=w(ixi^s,mag(:))
1744 
1745  select case(typediv)
1746  case("central")
1747  call divvector(bvec,ixi^l,ixo^l,divb,fourthorder)
1748  case("limited")
1749  call divvectors(bvec,ixi^l,ixo^l,divb)
1750  end select
1751  end subroutine get_divb
1752 
1753  !> get dimensionless div B = |divB| * volume / area / |B|
1754  subroutine get_normalized_divb(w,ixI^L,ixO^L,divb)
1757 
1758  integer, intent(in) :: ixI^L, ixO^L
1759  double precision, intent(in) :: w(ixi^s,1:nw)
1760  double precision :: divb(ixi^s), dsurface(ixi^s)
1761 
1762  integer :: ixA^L,idims
1763 
1764  call get_divb(w,ixi^l,ixo^l,divb)
1765  if(slab) then
1766  divb(ixo^s)=0.5d0*abs(divb(ixo^s))/sqrt(mhd_mag_en_all(w,ixi^l,ixo^l))/sum(1.d0/dxlevel(:))
1767  else
1768  ixamin^d=ixomin^d-1;
1769  ixamax^d=ixomax^d-1;
1770  dsurface(ixo^s)= sum(block%surfaceC(ixo^s,:),dim=ndim+1)
1771  do idims=1,ndim
1772  ixa^l=ixo^l-kr(idims,^d);
1773  dsurface(ixo^s)=dsurface(ixo^s)+block%surfaceC(ixa^s,idims)
1774  end do
1775  divb(ixo^s)=abs(divb(ixo^s))/sqrt(mhd_mag_en_all(w,ixi^l,ixo^l))*&
1776  block%dvolume(ixo^s)/dsurface(ixo^s)
1777  end if
1778 
1779  end subroutine get_normalized_divb
1780 
1781  !> Calculate idirmin and the idirmin:3 components of the common current array
1782  !> make sure that dxlevel(^D) is set correctly.
1783  subroutine get_current(w,ixI^L,ixO^L,idirmin,current)
1785  use mod_geometry
1786 
1787  integer :: idirmin0
1788  integer :: ixO^L, idirmin, ixI^L
1789  double precision :: w(ixi^s,1:nw)
1790  integer :: idir
1791 
1792  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
1793  double precision :: current(ixi^s,7-2*ndir:3),bvec(ixi^s,1:ndir)
1794 
1795  idirmin0 = 7-2*ndir
1796 
1797  bvec(ixi^s,1:ndir)=w(ixi^s,mag(1:ndir))
1798 
1799  call curlvector(bvec,ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
1800 
1801  if(b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
1802  block%J0(ixo^s,idirmin0:3)
1803 
1804  end subroutine get_current
1805 
1806  !> If resistivity is not zero, check diffusion time limit for dt
1807  subroutine mhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
1809  use mod_usr_methods
1811  use mod_viscosity, only: viscosity_get_dt
1812  use mod_gravity, only: gravity_get_dt
1813 
1814  integer, intent(in) :: ixI^L, ixO^L
1815  double precision, intent(inout) :: dtnew
1816  double precision, intent(in) :: dx^D
1817  double precision, intent(in) :: w(ixi^s,1:nw)
1818  double precision, intent(in) :: x(ixi^s,1:ndim)
1819 
1820  integer :: idirmin,idim
1821  double precision :: dxarr(ndim)
1822  double precision :: current(ixi^s,7-2*ndir:3),eta(ixi^s)
1823 
1824  dtnew = bigdouble
1825 
1826  ^d&dxarr(^d)=dx^d;
1827  if (mhd_eta>zero)then
1828  dtnew=dtdiffpar*minval(dxarr(1:ndim))**2/mhd_eta
1829  else if (mhd_eta<zero)then
1830  call get_current(w,ixi^l,ixo^l,idirmin,current)
1831  call usr_special_resistivity(w,ixi^l,ixo^l,idirmin,x,current,eta)
1832  dtnew=bigdouble
1833  do idim=1,ndim
1834  if(slab) then
1835  dtnew=min(dtnew,&
1836  dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
1837  else
1838  dtnew=min(dtnew,&
1839  dtdiffpar/(smalldouble+maxval(eta(ixo^s)/block%ds(ixo^s,idim)**2)))
1840  end if
1841  end do
1842  end if
1843 
1844  if(mhd_eta_hyper>zero) then
1845  if(slab) then
1846  dtnew=min(dtdiffpar*minval(dxarr(1:ndim))**4/mhd_eta_hyper,dtnew)
1847  else
1848  dtnew=min(dtdiffpar*minval(block%ds(ixo^s,1:ndim))**4/mhd_eta_hyper,dtnew)
1849  end if
1850  end if
1851 
1852  if(mhd_radiative_cooling) then
1853  call cooling_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
1854  end if
1855 
1856  if(mhd_viscosity) then
1857  call viscosity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
1858  end if
1859 
1860  if(mhd_gravity) then
1861  call gravity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
1862  end if
1863 
1864  end subroutine mhd_get_dt
1865 
1866  ! Add geometrical source terms to w
1867  subroutine mhd_add_source_geom(qdt,ixI^L,ixO^L,wCT,w,x)
1869 
1870  integer, intent(in) :: ixI^L, ixO^L
1871  double precision, intent(in) :: qdt, x(ixi^s,1:ndim)
1872  double precision, intent(inout) :: wCT(ixi^s,1:nw), w(ixi^s,1:nw)
1873 
1874  integer :: iw,idir, h1x^L{^NOONED, h2x^L}
1875  double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s)
1876 
1877  integer :: mr_,mphi_ ! Polar var. names
1878  integer :: br_,bphi_
1879 
1880  mr_=mom(1); mphi_=mom(1)-1+phi_ ! Polar var. names
1881  br_=mag(1); bphi_=mag(1)-1+phi_
1882 
1883  select case (typeaxial)
1884  case ('cylindrical')
1885  if (angmomfix) then
1886  call mpistop("angmomfix not implemented yet in MHD")
1887  endif
1888  call mhd_get_p_total(wct,x,ixi^l,ixo^l,tmp)
1889  if(phi_>0) then
1890  w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*(tmp(ixo^s)-&
1891  wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2/wct(ixo^s,rho_))
1892  w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt/x(ixo^s,1)*(&
1893  -wct(ixo^s,mphi_)*wct(ixo^s,mr_)/wct(ixo^s,rho_) &
1894  +wct(ixo^s,bphi_)*wct(ixo^s,br_))
1895  w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt/x(ixo^s,1)*&
1896  (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
1897  -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
1898  /wct(ixo^s,rho_)
1899  else
1900  w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*tmp(ixo^s)
1901  end if
1902  if(mhd_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,psi_)/x(ixo^s,1)
1903  case ('spherical')
1904  h1x^l=ixo^l-kr(1,^d); {^nooned h2x^l=ixo^l-kr(2,^d);}
1905  call mhd_get_p_total(wct,x,ixi^l,ixo^l,tmp1)
1906  tmp(ixo^s)=tmp1(ixo^s)
1907  if(b0field) then
1908  tmp2(ixo^s)=sum(block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=ndim+1)
1909  tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
1910  end if
1911  ! m1
1912  tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
1913  *(block%surfaceC(ixo^s,1)-block%surfaceC(h1x^s,1))/block%dvolume(ixo^s)
1914  if(ndir>1) then
1915  do idir=2,ndir
1916  tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,mom(idir))**2/wct(ixo^s,rho_)-wct(ixo^s,mag(idir))**2
1917  if(b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
1918  end do
1919  end if
1920  w(ixo^s,mom(1))=w(ixo^s,mom(1))+qdt*tmp(ixo^s)/x(ixo^s,1)
1921  ! b1
1922  if(mhd_glm) then
1923  w(ixo^s,mag(1))=w(ixo^s,mag(1))+qdt/x(ixo^s,1)*2.0d0*wct(ixo^s,psi_)
1924  end if
1925 
1926  {^nooned
1927  ! m2
1928  tmp(ixo^s)=tmp1(ixo^s)
1929  if(b0field) then
1930  tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
1931  end if
1932  ! This will make hydrostatic p=const an exact solution
1933  w(ixo^s,mom(2))=w(ixo^s,mom(2))+qdt*tmp(ixo^s) &
1934  *(block%surfaceC(ixo^s,2)-block%surfaceC(h2x^s,2)) &
1935  /block%dvolume(ixo^s)
1936  tmp(ixo^s)=-(wct(ixo^s,mom(1))*wct(ixo^s,mom(2))/wct(ixo^s,rho_) &
1937  -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
1938  if (b0field) then
1939  tmp(ixo^s)=tmp(ixo^s)+block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
1940  +wct(ixo^s,mag(1))*block%B0(ixo^s,2,0)
1941  end if
1942  if(ndir==3) then
1943  tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom(3))**2/wct(ixo^s,rho_) &
1944  -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
1945  if (b0field) then
1946  tmp(ixo^s)=tmp(ixo^s)-2.0d0*block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
1947  *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
1948  end if
1949  end if
1950  w(ixo^s,mom(2))=w(ixo^s,mom(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
1951  ! b2
1952  tmp(ixo^s)=(wct(ixo^s,mom(1))*wct(ixo^s,mag(2)) &
1953  -wct(ixo^s,mom(2))*wct(ixo^s,mag(1)))/wct(ixo^s,rho_)
1954  if(b0field) then
1955  tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom(1))*block%B0(ixo^s,2,0) &
1956  -wct(ixo^s,mom(2))*block%B0(ixo^s,1,0))/wct(ixo^s,rho_)
1957  end if
1958  if(mhd_glm) then
1959  tmp(ixo^s)=tmp(ixo^s) &
1960  + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,psi_)
1961  end if
1962  w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
1963  }
1964 
1965  if(ndir==3) then
1966  ! m3
1967  if(.not.angmomfix) then
1968  tmp(ixo^s)=-(wct(ixo^s,mom(3))*wct(ixo^s,mom(1))/wct(ixo^s,rho_) &
1969  -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
1970  -(wct(ixo^s,mom(2))*wct(ixo^s,mom(3))/wct(ixo^s,rho_) &
1971  -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
1972  *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
1973  if (b0field) then
1974  tmp(ixo^s)=tmp(ixo^s)+block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
1975  +wct(ixo^s,mag(1))*block%B0(ixo^s,3,0) {^nooned &
1976  +(block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
1977  +wct(ixo^s,mag(2))*block%B0(ixo^s,3,0)) &
1978  *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
1979  end if
1980  w(ixo^s,mom(3))=w(ixo^s,mom(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
1981  else
1982  call mpistop("angmomfix not implemented yet in MHD")
1983  end if
1984  ! b3
1985  tmp(ixo^s)=(wct(ixo^s,mom(1))*wct(ixo^s,mag(3)) &
1986  -wct(ixo^s,mom(3))*wct(ixo^s,mag(1)))/wct(ixo^s,rho_) {^nooned &
1987  -(wct(ixo^s,mom(3))*wct(ixo^s,mag(2)) &
1988  -wct(ixo^s,mom(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
1989  /(wct(ixo^s,rho_)*dsin(x(ixo^s,2))) }
1990  if (b0field) then
1991  tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom(1))*block%B0(ixo^s,3,0) &
1992  -wct(ixo^s,mom(3))*block%B0(ixo^s,1,0))/wct(ixo^s,rho_){^nooned &
1993  -(wct(ixo^s,mom(3))*block%B0(ixo^s,2,0) &
1994  -wct(ixo^s,mom(2))*block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
1995  /(wct(ixo^s,rho_)*dsin(x(ixo^s,2))) }
1996  end if
1997  w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
1998  end if
1999  end select
2000  end subroutine mhd_add_source_geom
2001 
2002  !> Compute 2 times total magnetic energy
2003  function mhd_mag_en_all(w, ixI^L, ixO^L) result(mge)
2005  integer, intent(in) :: ixI^L, ixO^L
2006  double precision, intent(in) :: w(ixi^s, nw)
2007  double precision :: mge(ixo^s)
2008 
2009  if (b0field) then
2010  mge = sum((w(ixo^s, mag(:))+block%B0(ixo^s,:,block%iw0))**2, dim=ndim+1)
2011  else
2012  mge = sum(w(ixo^s, mag(:))**2, dim=ndim+1)
2013  end if
2014  end function mhd_mag_en_all
2015 
2016  !> Compute full magnetic field by direction
2017  function mhd_mag_i_all(w, ixI^L, ixO^L,idir) result(mgf)
2019  integer, intent(in) :: ixI^L, ixO^L, idir
2020  double precision, intent(in) :: w(ixi^s, nw)
2021  double precision :: mgf(ixo^s)
2022 
2023  if (b0field) then
2024  mgf = w(ixo^s, mag(idir))+block%B0(ixo^s,idir,block%iw0)
2025  else
2026  mgf = w(ixo^s, mag(idir))
2027  end if
2028  end function mhd_mag_i_all
2029 
2030  !> Compute evolving magnetic energy
2031  function mhd_mag_en(w, ixI^L, ixO^L) result(mge)
2033  integer, intent(in) :: ixI^L, ixO^L
2034  double precision, intent(in) :: w(ixi^s, nw)
2035  double precision :: mge(ixo^s)
2036 
2037  mge = 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
2038  end function mhd_mag_en
2039 
2040  !> compute kinetic energy
2041  function mhd_kin_en(w, ixI^L, ixO^L, inv_rho) result(ke)
2043  integer, intent(in) :: ixI^L, ixO^L
2044  double precision, intent(in) :: w(ixi^s, nw)
2045  double precision :: ke(ixo^s)
2046  double precision, intent(in), optional :: inv_rho(ixo^s)
2047 
2048  if (present(inv_rho)) then
2049  ke = 0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) * inv_rho
2050  else
2051  ke = 0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) / w(ixo^s, rho_)
2052  end if
2053  end function mhd_kin_en
2054 
2055  subroutine mhd_getv_hall(w,x,ixI^L,ixO^L,vHall)
2057 
2058  integer, intent(in) :: ixI^L, ixO^L
2059  double precision, intent(in) :: w(ixi^s,nw)
2060  double precision, intent(in) :: x(ixi^s,1:ndim)
2061  double precision, intent(inout) :: vHall(ixi^s,1:3)
2062 
2063  integer :: idir, idirmin
2064  double precision :: current(ixi^s,7-2*ndir:3)
2065 
2066  ! Calculate current density and idirmin
2067  call get_current(w,ixi^l,ixo^l,idirmin,current)
2068  vhall(ixo^s,1:3) = zero
2069  vhall(ixo^s,idirmin:3) = - mhd_etah*current(ixo^s,idirmin:3)
2070  do idir = idirmin, 3
2071  vhall(ixo^s,idir) = vhall(ixo^s,idir)/w(ixo^s,rho_)
2072  end do
2073 
2074  end subroutine mhd_getv_hall
2075 
2076  subroutine mhd_getdt_hall(w,x,ixI^L,ixO^L,dx^D,dthall)
2078 
2079  integer, intent(in) :: ixI^L, ixO^L
2080  double precision, intent(in) :: dx^D
2081  double precision, intent(in) :: w(ixi^s,1:nw)
2082  double precision, intent(in) :: x(ixi^s,1:ndim)
2083  double precision, intent(out) :: dthall
2084  !.. local ..
2085  double precision :: dxarr(ndim)
2086  double precision :: bmag(ixi^s)
2087 
2088  dthall=bigdouble
2089 
2090  ! because we have that in cmax now:
2091  return
2092 
2093  ^d&dxarr(^d)=dx^d;
2094 
2095  if (.not. b0field) then
2096  bmag(ixo^s)=sqrt(sum(w(ixo^s,mag(:))**2, dim=ndim+1))
2097  bmag(ixo^s)=sqrt(sum((w(ixo^s,mag(:)) + block%B0(ixo^s,1:ndir,block%iw0))**2))
2098  end if
2099 
2100  if(slab) then
2101  dthall=dtdiffpar*minval(dxarr(1:ndim))**2.0d0/(mhd_etah*maxval(bmag(ixo^s)/w(ixo^s,rho_)))
2102  else
2103  dthall=dtdiffpar*minval(block%ds(ixo^s,1:ndim))**2.0d0/(mhd_etah*maxval(bmag(ixo^s)/w(ixo^s,rho_)))
2104  end if
2105 
2106  end subroutine mhd_getdt_hall
2107 
2108  !> This implements eq. (42) in Dedner et al. 2002 JcP 175
2109  !> Gives the Riemann solution on the interface
2110  !> for the normal B component and Psi in the GLM-MHD system.
2111  !> 23/04/2013 Oliver Porth
2112  subroutine glmsolve(wLC,wRC,ixI^L,ixO^L,idir)
2114  double precision, intent(inout) :: wLC(ixi^s,1:nw), wRC(ixi^s,1:nw)
2115  integer, intent(in) :: ixI^L, ixO^L, idir
2116  double precision :: dB(ixi^s), dPsi(ixi^s)
2117 
2118  db(ixo^s) = wrc(ixo^s,mag(idir)) - wlc(ixo^s,mag(idir))
2119  dpsi(ixo^s) = wrc(ixo^s,psi_) - wlc(ixo^s,psi_)
2120 
2121  wlc(ixo^s,mag(idir)) = 0.5d0 * (wrc(ixo^s,mag(idir)) + wlc(ixo^s,mag(idir))) &
2122  - 0.5d0/cmax_global * dpsi(ixo^s)
2123  wlc(ixo^s,psi_) = 0.5d0 * (wrc(ixo^s,psi_) + wlc(ixo^s,psi_)) &
2124  - 0.5d0*cmax_global * db(ixo^s)
2125 
2126  wrc(ixo^s,mag(idir)) = wlc(ixo^s,mag(idir))
2127  wrc(ixo^s,psi_) = wlc(ixo^s,psi_)
2128 
2129  end subroutine glmsolve
2130 
2131  subroutine mhd_boundary_adjust
2133  integer :: iB, idim, iside, iigrid, igrid
2134  integer :: ixG^L, ixO^L, i^D
2135 
2136  ixg^l=ixg^ll;
2137  do iigrid=1,igridstail; igrid=igrids(iigrid);
2138  if(.not.phyboundblock(igrid)) cycle
2139  saveigrid=igrid
2140  block=>ps(igrid)
2141  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
2142  do idim=1,ndim
2143  ! to avoid using as yet unknown corner info in more than 1D, we
2144  ! fill only interior mesh ranges of the ghost cell ranges at first,
2145  ! and progressively enlarge the ranges to include corners later
2146  do iside=1,2
2147  i^d=kr(^d,idim)*(2*iside-3);
2148  if (neighbor_type(i^d,igrid)/=1) cycle
2149  ib=(idim-1)*2+iside
2150  if(.not.boundary_divbfix(ib)) cycle
2151  if(any(typeboundary(:,ib)=="special")) then
2152  ! MF nonlinear force-free B field extrapolation and data driven
2153  ! require normal B of the first ghost cell layer to be untouched by
2154  ! fixdivB=0 process, set boundary_divbfix_skip(iB)=1 in par file
2155  select case (idim)
2156  {case (^d)
2157  if (iside==2) then
2158  ! maximal boundary
2159  ixomin^dd=ixgmax^d+1-nghostcells+boundary_divbfix_skip(2*^d)^d%ixOmin^dd=ixgmin^dd;
2160  ixomax^dd=ixgmax^dd;
2161  else
2162  ! minimal boundary
2163  ixomin^dd=ixgmin^dd;
2164  ixomax^dd=ixgmin^d-1+nghostcells-boundary_divbfix_skip(2*^d-1)^d%ixOmax^dd=ixgmax^dd;
2165  end if \}
2166  end select
2167  call fixdivb_boundary(ixg^l,ixo^l,ps(igrid)%w,ps(igrid)%x,ib)
2168  end if
2169  end do
2170  end do
2171  end do
2172 
2173  end subroutine mhd_boundary_adjust
2174 
2175  subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
2177 
2178  integer, intent(in) :: ixG^L,ixO^L,iB
2179  double precision, intent(inout) :: w(ixg^s,1:nw)
2180  double precision, intent(in) :: x(ixg^s,1:ndim)
2181 
2182  double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
2183  integer :: ix^D,ixF^L
2184 
2185  select case(ib)
2186  case(1)
2187  ! 2nd order CD for divB=0 to set normal B component better
2188  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_primitive(ixg^l,ixo^l,w,x)
2189  {^iftwod
2190  ixfmin1=ixomin1+1
2191  ixfmax1=ixomax1+1
2192  ixfmin2=ixomin2+1
2193  ixfmax2=ixomax2-1
2194  if(slab) then
2195  dx1x2=dxlevel(1)/dxlevel(2)
2196  do ix1=ixfmax1,ixfmin1,-1
2197  w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
2198  +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
2199  w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
2200  enddo
2201  else
2202  do ix1=ixfmax1,ixfmin1,-1
2203  w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
2204  w(ix1,ixfmin2:ixfmax2,mag(1)))*block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
2205  +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
2206  block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
2207  -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
2208  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
2209  /block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
2210  end do
2211  end if
2212  }
2213  {^ifthreed
2214  ixfmin1=ixomin1+1
2215  ixfmax1=ixomax1+1
2216  ixfmin2=ixomin2+1
2217  ixfmax2=ixomax2-1
2218  ixfmin3=ixomin3+1
2219  ixfmax3=ixomax3-1
2220  if(slab) then
2221  dx1x2=dxlevel(1)/dxlevel(2)
2222  dx1x3=dxlevel(1)/dxlevel(3)
2223  do ix1=ixfmax1,ixfmin1,-1
2224  w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
2225  w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
2226  +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
2227  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
2228  +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
2229  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
2230  end do
2231  else
2232  do ix1=ixfmax1,ixfmin1,-1
2233  w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
2234  ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
2235  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
2236  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
2237  +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
2238  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
2239  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
2240  -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
2241  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
2242  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
2243  +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
2244  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
2245  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
2246  -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
2247  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
2248  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
2249  /block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
2250  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
2251  end do
2252  end if
2253  }
2254  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_conserved(ixg^l,ixo^l,w,x)
2255  case(2)
2256  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_primitive(ixg^l,ixo^l,w,x)
2257  {^iftwod
2258  ixfmin1=ixomin1-1
2259  ixfmax1=ixomax1-1
2260  ixfmin2=ixomin2+1
2261  ixfmax2=ixomax2-1
2262  if(slab) then
2263  dx1x2=dxlevel(1)/dxlevel(2)
2264  do ix1=ixfmin1,ixfmax1
2265  w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
2266  -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
2267  w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
2268  enddo
2269  else
2270  do ix1=ixfmin1,ixfmax1
2271  w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
2272  w(ix1,ixfmin2:ixfmax2,mag(1)))*block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
2273  -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
2274  block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
2275  +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
2276  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
2277  /block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
2278  end do
2279  end if
2280  }
2281  {^ifthreed
2282  ixfmin1=ixomin1-1
2283  ixfmax1=ixomax1-1
2284  ixfmin2=ixomin2+1
2285  ixfmax2=ixomax2-1
2286  ixfmin3=ixomin3+1
2287  ixfmax3=ixomax3-1
2288  if(slab) then
2289  dx1x2=dxlevel(1)/dxlevel(2)
2290  dx1x3=dxlevel(1)/dxlevel(3)
2291  do ix1=ixfmin1,ixfmax1
2292  w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
2293  w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
2294  -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
2295  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
2296  -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
2297  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
2298  end do
2299  else
2300  do ix1=ixfmin1,ixfmax1
2301  w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
2302  ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
2303  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
2304  block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
2305  -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
2306  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
2307  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
2308  +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
2309  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
2310  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
2311  -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
2312  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
2313  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
2314  +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
2315  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
2316  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
2317  /block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
2318  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
2319  end do
2320  end if
2321  }
2322  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_conserved(ixg^l,ixo^l,w,x)
2323  case(3)
2324  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_primitive(ixg^l,ixo^l,w,x)
2325  {^iftwod
2326  ixfmin1=ixomin1+1
2327  ixfmax1=ixomax1-1
2328  ixfmin2=ixomin2+1
2329  ixfmax2=ixomax2+1
2330  if(slab) then
2331  dx2x1=dxlevel(2)/dxlevel(1)
2332  do ix2=ixfmax2,ixfmin2,-1
2333  w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
2334  +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
2335  w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
2336  enddo
2337  else
2338  do ix2=ixfmax2,ixfmin2,-1
2339  w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
2340  w(ixfmin1:ixfmax1,ix2,mag(2)))*block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
2341  +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
2342  block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
2343  -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
2344  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
2345  /block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
2346  end do
2347  end if
2348  }
2349  {^ifthreed
2350  ixfmin1=ixomin1+1
2351  ixfmax1=ixomax1-1
2352  ixfmin3=ixomin3+1
2353  ixfmax3=ixomax3-1
2354  ixfmin2=ixomin2+1
2355  ixfmax2=ixomax2+1
2356  if(slab) then
2357  dx2x1=dxlevel(2)/dxlevel(1)
2358  dx2x3=dxlevel(2)/dxlevel(3)
2359  do ix2=ixfmax2,ixfmin2,-1
2360  w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
2361  ix2+1,ixfmin3:ixfmax3,mag(2)) &
2362  +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
2363  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
2364  +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
2365  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
2366  end do
2367  else
2368  do ix2=ixfmax2,ixfmin2,-1
2369  w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
2370  ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
2371  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
2372  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
2373  +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
2374  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
2375  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
2376  -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
2377  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
2378  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
2379  +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
2380  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
2381  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
2382  -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
2383  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
2384  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
2385  /block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
2386  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
2387  end do
2388  end if
2389  }
2390  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_conserved(ixg^l,ixo^l,w,x)
2391  case(4)
2392  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_primitive(ixg^l,ixo^l,w,x)
2393  {^iftwod
2394  ixfmin1=ixomin1+1
2395  ixfmax1=ixomax1-1
2396  ixfmin2=ixomin2-1
2397  ixfmax2=ixomax2-1
2398  if(slab) then
2399  dx2x1=dxlevel(2)/dxlevel(1)
2400  do ix2=ixfmin2,ixfmax2
2401  w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
2402  -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
2403  w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
2404  end do
2405  else
2406  do ix2=ixfmin2,ixfmax2
2407  w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
2408  w(ixfmin1:ixfmax1,ix2,mag(2)))*block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
2409  -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
2410  block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
2411  +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
2412  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
2413  /block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
2414  end do
2415  end if
2416  }
2417  {^ifthreed
2418  ixfmin1=ixomin1+1
2419  ixfmax1=ixomax1-1
2420  ixfmin3=ixomin3+1
2421  ixfmax3=ixomax3-1
2422  ixfmin2=ixomin2-1
2423  ixfmax2=ixomax2-1
2424  if(slab) then
2425  dx2x1=dxlevel(2)/dxlevel(1)
2426  dx2x3=dxlevel(2)/dxlevel(3)
2427  do ix2=ixfmin2,ixfmax2
2428  w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
2429  ix2-1,ixfmin3:ixfmax3,mag(2)) &
2430  -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
2431  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
2432  -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
2433  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
2434  end do
2435  else
2436  do ix2=ixfmin2,ixfmax2
2437  w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
2438  ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
2439  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
2440  block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
2441  -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
2442  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
2443  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
2444  +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
2445  w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
2446  block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
2447  -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
2448  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
2449  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
2450  +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
2451  w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
2452  block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
2453  /block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
2454  w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
2455  end do
2456  end if
2457  }
2458  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_conserved(ixg^l,ixo^l,w,x)
2459  {^ifthreed
2460  case(5)
2461  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_primitive(ixg^l,ixo^l,w,x)
2462  ixfmin1=ixomin1+1
2463  ixfmax1=ixomax1-1
2464  ixfmin2=ixomin2+1
2465  ixfmax2=ixomax2-1
2466  ixfmin3=ixomin3+1
2467  ixfmax3=ixomax3+1
2468  if(slab) then
2469  dx3x1=dxlevel(3)/dxlevel(1)
2470  dx3x2=dxlevel(3)/dxlevel(2)
2471  do ix3=ixfmax3,ixfmin3,-1
2472  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
2473  ixfmin2:ixfmax2,ix3+1,mag(3)) &
2474  +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
2475  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
2476  +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
2477  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
2478  end do
2479  else
2480  do ix3=ixfmax3,ixfmin3,-1
2481  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
2482  ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
2483  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
2484  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
2485  +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
2486  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
2487  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
2488  -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
2489  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
2490  block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
2491  +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
2492  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
2493  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
2494  -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
2495  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
2496  block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
2497  /block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
2498  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
2499  end do
2500  end if
2501  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_conserved(ixg^l,ixo^l,w,x)
2502  case(6)
2503  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_primitive(ixg^l,ixo^l,w,x)
2504  ixfmin1=ixomin1+1
2505  ixfmax1=ixomax1-1
2506  ixfmin2=ixomin2+1
2507  ixfmax2=ixomax2-1
2508  ixfmin3=ixomin3-1
2509  ixfmax3=ixomax3-1
2510  if(slab) then
2511  dx3x1=dxlevel(3)/dxlevel(1)
2512  dx3x2=dxlevel(3)/dxlevel(2)
2513  do ix3=ixfmin3,ixfmax3
2514  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
2515  ixfmin2:ixfmax2,ix3-1,mag(3)) &
2516  -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
2517  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
2518  -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
2519  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
2520  end do
2521  else
2522  do ix3=ixfmin3,ixfmax3
2523  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
2524  ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
2525  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
2526  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
2527  -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
2528  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
2529  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
2530  +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
2531  w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
2532  block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
2533  -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
2534  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
2535  block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
2536  +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
2537  w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
2538  block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
2539  /block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
2540  w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
2541  end do
2542  end if
2543  if(mhd_energy.and..not.block%e_is_internal) call mhd_to_conserved(ixg^l,ixo^l,w,x)
2544  }
2545  case default
2546  call mpistop("Special boundary is not defined for this region")
2547  end select
2548 
2549  end subroutine fixdivb_boundary
2550 
2551  {^nooned
2552  subroutine mhd_clean_divb_multigrid(qdt, qt, active)
2556  use mod_geometry
2557  double precision, intent(in) :: qdt !< Current time step
2558  double precision, intent(in) :: qt !< Current time
2559  logical, intent(inout) :: active !< Output if the source is active
2560  integer :: iigrid, igrid, id
2561  integer :: n, nc, lvl, ix^L, idim
2562  type(tree_node), pointer :: pnode
2563  double precision :: tmp(ixg^t), grad(ixg^t, ndim)
2564  double precision :: res
2565  double precision, parameter :: max_residual = 1d-3
2566  integer, parameter :: max_its = 10
2567 
2568  mg%operator_type = mg_laplacian
2569 
2570  ! Set boundary conditions
2571  do n = 1, 2*ndim
2572  idim = (n+1)/2
2573  select case (typeboundary(mag(idim), n))
2574  case ('symm')
2575  ! d/dx B = 0, take phi = 0
2576  mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
2577  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
2578  case ('asymm')
2579  ! B = 0, so grad(phi) = 0
2580  mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
2581  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
2582  case ('cont')
2583  mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
2584  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
2585  case ('special')
2586  ! Assume Dirichlet boundary conditions, derivative zero
2587  mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
2588  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
2589  case ('periodic')
2590  ! Nothing to do here
2591  case default
2592  print *, "divb_multigrid warning: unknown b.c.: ", &
2593  trim(typeboundary(mag(idim), n))
2594  mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
2595  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
2596  end select
2597  end do
2598 
2599  ix^l=ixm^ll^ladd1;
2600 
2601  ! Store divergence of B as right-hand side
2602  do iigrid = 1, igridstail
2603  igrid = igrids(iigrid);
2604  pnode => igrid_to_node(igrid, mype)%node
2605  id = pnode%id
2606  lvl = mg%boxes(id)%lvl
2607  nc = mg%box_size_lvl(lvl)
2608 
2609  ! Geometry subroutines expect this to be set
2610  block => ps(igrid)
2611  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
2612 
2613  call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^ll, ixm^ll, tmp, &
2615  mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(ixm^t)
2616  end do
2617 
2618  ! Solve laplacian(phi) = divB
2619  do n = 1, max_its
2620  call mg_fas_vcycle(mg, max_res=res)
2621  if (res < max_residual) exit
2622  end do
2623 
2624  if (res > max_residual) call mpistop("divb_multigrid: no convergence")
2625 
2626  ! Correct the magnetic field
2627  do iigrid = 1, igridstail
2628  igrid = igrids(iigrid);
2629  pnode => igrid_to_node(igrid, mype)%node
2630  id = pnode%id
2631 
2632  ! Geometry subroutines expect this to be set
2633  block => ps(igrid)
2634  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
2635 
2636  ! Compute the gradient of phi
2637  tmp(ix^s) = mg%boxes(id)%cc({:,}, mg_iphi)
2638  do idim = 1, ndim
2639  call gradient(tmp,ixg^ll,ixm^ll,idim,grad(ixg^t, idim))
2640  end do
2641 
2642  ! Apply the correction B* = B - gradient(phi)
2643  tmp(ixm^t) = sum(ps(igrid)%w(ixm^t, mag(1:ndim))**2, dim=ndim+1)
2644  ps(igrid)%w(ixm^t, mag(1:ndim)) = &
2645  ps(igrid)%w(ixm^t, mag(1:ndim)) - grad(ixm^t, :)
2646 
2647  ! Determine magnetic energy difference
2648  tmp(ixm^t) = 0.5_dp * (sum(ps(igrid)%w(ixm^t, &
2649  mag(1:ndim))**2, dim=ndim+1) - tmp(ixm^t))
2650 
2651  if (mhd_energy) then
2652  ! Keep thermal pressure the same
2653  ps(igrid)%w(ixm^t, e_) = ps(igrid)%w(ixm^t, e_) + tmp(ixm^t)
2654 
2655  ! Possible alternative: keep total pressure the same
2656  ! ps(igrid)%w(ixM^T, e_) = ps(igrid)%w(ixM^T, e_) + &
2657  ! (mhd_gamma-2) * inv_gamma_1 * tmp(ixM^T)
2658  end if
2659  end do
2660 
2661  active = .true.
2662 
2663  end subroutine mhd_clean_divb_multigrid
2664  }
2665 
2666 end module mod_mhd_phys
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 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:46
subroutine mhd_clean_divb_multigrid(qdt, qt, active)
double precision vmax_global
global fastest flow speed needed in glm method
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:54
integer, public, protected p_
Index of the gas pressure (-1 if not present) should equal e_.
Definition: mod_mhd_phys.t:57
subroutine add_source_glm2(qdt, ixIL, ixOL, wCT, w, x)
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:193
double precision, public mhd_eta
The MHD resistivity.
Definition: mod_mhd_phys.t:75
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
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
Definition: mod_geometry.t:255
subroutine divvectors(qvec, ixIL, ixOL, divq)
Calculate divergence of a vector qvec within ixL using limited extrapolation to cell edges...
Definition: mod_geometry.t:430
subroutine glmsolve(wLC, wRC, ixIL, ixOL, idir)
This implements eq. (42) in Dedner et al. 2002 JcP 175 Gives the Riemann solution on the interface fo...
subroutine radiative_cooling_add_source(qdt, ixIL, ixOL, wCT, w, x, qsourcesplit, active)
subroutine mhd_add_source_geom(qdt, ixIL, ixOL, wCT, w, x)
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:645
procedure(sub_check_w), pointer phys_check_w
Definition: mod_physics.t:48
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:53
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:117
subroutine mhd_getv_hall(w, x, ixIL, ixOL, vHall)
character(len=std_len) typeaxial
procedure(sub_boundary_adjust), pointer phys_boundary_adjust
Definition: mod_physics.t:50
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:771
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor, but its use is kept to a minimum.
subroutine add_source_glm1(qdt, ixIL, ixOL, wCT, w, x)
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
Definition: mod_mhd_phys.t:111
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:496
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:93
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
integer, public, protected psi_
Indices of the GLM psi.
Definition: mod_mhd_phys.t:63
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:69
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:448
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:60
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:892
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
procedure(sub_angmomfix), pointer phys_angmomfix
Definition: mod_physics.t:52
Thermal conduction for HD and MHD.
procedure(sub_convert), pointer phys_to_conserved
Definition: mod_physics.t:36
integer ixm
the mesh range (within a block with ghost cells)
subroutine curlvector(qvec, ixIL, ixOL, curlvec, idirmin, idirmin0, ndir0)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
Definition: mod_geometry.t:485
subroutine, public mhd_get_v(w, x, ixIL, ixOL, v)
Calculate v vector.
Definition: mod_mhd_phys.t:681
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:176
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:723
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:527
subroutine mhd_get_csound_prim(w, x, ixIL, ixOL, idim, csound)
Calculate fast magnetosonic wave speed.
Definition: mod_mhd_phys.t:810
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:51
double precision, public mhd_etah
TODO: what is this?
Definition: mod_mhd_phys.t:81
integer, parameter unitpar
integer, dimension(3, 3) kr
Kronecker delta tensor.
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:114
character(len=std_len) typegrad
character(len=std_len), public, protected typedivbfix
Method type to clean divergence of B.
Definition: mod_mhd_phys.t:90
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
integer, dimension(:), allocatable, parameter d
subroutine divvector(qvec, ixIL, ixOL, divq, fourthorder)
Calculate divergence of a vector qvec within ixL.
Definition: mod_geometry.t:359
procedure(sub_get_flux), pointer phys_get_flux
Definition: mod_physics.t:41
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:709
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
Module for including gravity in (magneto)hydrodynamics simulations.
Definition: mod_gravity.t:2
subroutine mhd_physical_units()
Definition: mod_mhd_phys.t:467
subroutine rhos_to_e(ixIL, ixOL, w, x)
Convert entropy to energy.
Definition: mod_mhd_phys.t:663
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:42
subroutine internal_energy_add_source(qdt, ixIL, ixOL, wCT, w, x)
subroutine mhd_get_v_idim(w, x, ixIL, ixOL, idim, v)
Calculate v component.
Definition: mod_mhd_phys.t:697
Module for handling problematic values in simulations, such as negative pressures.
double precision unit_pressure
Physical scaling factor for pressure.
double precision unit_magneticfield
Physical scaling factor for magnetic field.
logical use_particles
Use particles module or not.
double precision unit_temperature
Physical scaling factor for temperature.
integer, public, protected rho_
Index of the density (in the w array)
Definition: mod_mhd_phys.t:48
procedure(sub_add_source_geom), pointer phys_add_source_geom
Definition: mod_physics.t:44
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:876
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:160
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:35
double precision, public mhd_eta_hyper
The MHD hyper-resistivity.
Definition: mod_mhd_phys.t:78
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:66
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:853
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:302
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:594
logical slab
Cartesian geometry or not.
subroutine, public mhd_phys_init()
Definition: mod_mhd_phys.t:256
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:42
procedure(sub_get_v_idim), pointer phys_get_v_idim
Definition: mod_physics.t:43
character(len=std_len) typediv
logical, public divbwave
Add divB wave in Roe solver.
Definition: mod_mhd_phys.t:108
integer, public, protected mhd_n_tracer
Number of tracer species.
Definition: mod_mhd_phys.t:45
procedure(sub_get_cmax), pointer phys_get_cmax
Definition: mod_physics.t:39
The module add viscous source terms and check time step.
Definition: mod_viscosity.t:10
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:45
procedure(sub_convert), pointer phys_to_primitive
Definition: mod_physics.t:37
subroutine viscosity_init(phys_wider_stencil, phys_req_diagonal)
Initialize the module.
Definition: mod_viscosity.t:57
double precision function, dimension(ixo^s) mhd_mag_en_all(w, ixIL, ixOL)
Compute 2 times total magnetic energy.
subroutine thermal_conduction_init(phys_gamma)
Initialize the module.
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:120
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
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, public mhd_adiab
The adiabatic constant.
Definition: mod_mhd_phys.t:72
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:49
subroutine, public mhd_to_primitive(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
Definition: mod_mhd_phys.t:558
procedure(sub_get_cbounds), pointer phys_get_cbounds
Definition: mod_physics.t:40
procedure(sub_modify_wlr), pointer phys_modify_wlr
Definition: mod_physics.t:38
procedure(sub_write_info), pointer phys_write_info
Definition: mod_physics.t:51
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:907
pure subroutine cross_product(ixIL, ixOL, a, b, axb)
Cross product of two vectors.
logical need_global_vmax
need global maximal flow speed
double precision function, dimension(ixo^s) mhd_mag_i_all(w, ixIL, ixOL, idir)
Compute full magnetic field by direction.