MPI-AMRVAC  3.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 
4 #include "amrvac.h"
5 
6  use mod_global_parameters, only: std_len, const_c
10  use mod_physics
11 
12  implicit none
13  private
14 
15  !> Whether an energy equation is used
16  logical, public, protected :: mhd_energy = .true.
17 
18  !> Whether thermal conduction is used
19  logical, public, protected :: mhd_thermal_conduction = .false.
20  !> type of fluid for thermal conduction
21  type(tc_fluid), public, allocatable :: tc_fl
22  type(te_fluid), public, allocatable :: te_fl_mhd
23 
24  !> Whether radiative cooling is added
25  logical, public, protected :: mhd_radiative_cooling = .false.
26  !> type of fluid for radiative cooling
27  type(rc_fluid), public, allocatable :: rc_fl
28 
29  !> Whether viscosity is added
30  logical, public, protected :: mhd_viscosity = .false.
31 
32  !> Whether gravity is added
33  logical, public, protected :: mhd_gravity = .false.
34 
35  !> Whether Hall-MHD is used
36  logical, public, protected :: mhd_hall = .false.
37 
38  !> Whether Ambipolar term is used
39  logical, public, protected :: mhd_ambipolar = .false.
40 
41  !> Whether Ambipolar term is implemented using supertimestepping
42  logical, public, protected :: mhd_ambipolar_sts = .false.
43 
44  !> Whether Ambipolar term is implemented explicitly
45  logical, public, protected :: mhd_ambipolar_exp = .false.
46 
47  !> Whether particles module is added
48  logical, public, protected :: mhd_particles = .false.
49 
50  !> Whether magnetofriction is added
51  logical, public, protected :: mhd_magnetofriction = .false.
52 
53  !> Whether GLM-MHD is used to control div B
54  logical, public, protected :: mhd_glm = .false.
55 
56  !> Whether extended GLM-MHD is used with additional sources
57  logical, public, protected :: mhd_glm_extended = .true.
58 
59  !> Whether TRAC method is used
60  logical, public, protected :: mhd_trac = .false.
61 
62  !> Which TRAC method is used
63  integer, public, protected :: mhd_trac_type=1
64 
65  !> Height of the mask used in the TRAC method
66  double precision, public, protected :: mhd_trac_mask = 0.d0
67 
68  !> Distance between two adjacent traced magnetic field lines (in finest cell size)
69  integer, public, protected :: mhd_trac_finegrid=4
70 
71  !> Whether auxiliary internal energy is solved
72  logical, public, protected :: mhd_solve_eaux = .false.
73 
74  !> Whether internal energy is solved instead of total energy
75  logical, public, protected :: mhd_internal_e = .false.
76 
77  !TODO this does not work with the splitting: check mhd_check_w_hde and mhd_handle_small_values_hde
78  !> Whether hydrodynamic energy is solved instead of total energy
79  logical, public, protected :: mhd_hydrodynamic_e = .false.
80 
81  !> Whether divB cleaning sources are added splitting from fluid solver
82  logical, public, protected :: source_split_divb = .false.
83 
84  !> GLM-MHD parameter: ratio of the diffusive and advective time scales for div b
85  !> taking values within [0, 1]
86  double precision, public :: mhd_glm_alpha = 0.5d0
87 
88  !TODO this does not work with the splitting: check mhd_check_w_semirelati and mhd_handle_small_values_semirelati
89  !> Whether semirelativistic MHD equations (Gombosi 2002 JCP) are solved
90  logical, public, protected :: mhd_semirelativistic = .false.
91 
92  !> Whether boris simplified semirelativistic MHD equations (Gombosi 2002 JCP) are solved
93  logical, public, protected :: mhd_boris_simplification = .false.
94 
95  !> Reduced speed of light for semirelativistic MHD
96  double precision, public, protected :: mhd_reduced_c = const_c
97 
98  !> Whether CAK radiation line force is activated
99  logical, public, protected :: mhd_cak_force = .false.
100 
101  !> MHD fourth order
102  logical, public, protected :: mhd_4th_order = .false.
103 
104  !> whether split off equilibrium density
105  logical, public :: has_equi_rho0 = .false.
106  !> whether split off equilibrium thermal pressure
107  logical, public :: has_equi_pe0 = .false.
108  logical, public :: mhd_equi_thermal = .false.
109 
110  !> equi vars indices in the state%equi_vars array
111  integer, public :: equi_rho0_ = -1
112  integer, public :: equi_pe0_ = -1
113 
114  !> whether dump full variables (when splitting is used) in a separate dat file
115  logical, public, protected :: mhd_dump_full_vars = .false.
116 
117  !> Number of tracer species
118  integer, public, protected :: mhd_n_tracer = 0
119 
120  !> Index of the density (in the w array)
121  integer, public, protected :: rho_
122 
123  !> Indices of the momentum density
124  integer, allocatable, public, protected :: mom(:)
125 
126  !> Index of the energy density (-1 if not present)
127  integer, public, protected :: e_
128 
129  !> Index of the gas pressure (-1 if not present) should equal e_
130  integer, public, protected :: p_
131 
132  !> Indices of the magnetic field
133  integer, allocatable, public, protected :: mag(:)
134 
135  !> Indices of the GLM psi
136  integer, public, protected :: psi_
137 
138  !> Indices of auxiliary internal energy
139  integer, public, protected :: eaux_
140  integer, public, protected :: paux_
141 
142  !> Index of the cutoff temperature for the TRAC method
143  integer, public, protected :: tcoff_
144  integer, public, protected :: tweight_
145 
146  !> Indices of the tracers
147  integer, allocatable, public, protected :: tracer(:)
148 
149  !> The adiabatic index
150  double precision, public :: mhd_gamma = 5.d0/3.0d0
151 
152  !> The adiabatic constant
153  double precision, public :: mhd_adiab = 1.0d0
154 
155  !> The MHD resistivity
156  double precision, public :: mhd_eta = 0.0d0
157 
158  !> The MHD hyper-resistivity
159  double precision, public :: mhd_eta_hyper = 0.0d0
160 
161  !> TODO: what is this?
162  double precision, public :: mhd_etah = 0.0d0
163 
164  !> The MHD ambipolar coefficient
165  double precision, public :: mhd_eta_ambi = 0.0d0
166 
167  !> The small_est allowed energy
168  double precision, protected :: small_e
169 
170  !> The number of waves
171  integer :: nwwave=8
172 
173  !> Method type to clean divergence of B
174  character(len=std_len), public, protected :: typedivbfix = 'linde'
175 
176  !> Method type of constrained transport
177  character(len=std_len), public, protected :: type_ct = 'uct_contact'
178 
179  !> Whether divB is computed with a fourth order approximation
180  logical, public, protected :: mhd_divb_4thorder = .false.
181 
182  !> Method type in a integer for good performance
183  integer :: type_divb
184 
185  !> Coefficient of diffusive divB cleaning
186  double precision :: divbdiff = 0.8d0
187 
188  !> Update all equations due to divB cleaning
189  character(len=std_len) :: typedivbdiff = 'all'
190 
191  !> Use a compact way to add resistivity
192  logical :: compactres = .false.
193 
194  !> Add divB wave in Roe solver
195  logical, public :: divbwave = .true.
196 
197  !> clean initial divB
198  logical, public :: clean_initial_divb = .false.
199 
200  !> Helium abundance over Hydrogen
201  double precision, public, protected :: he_abundance=0.1d0
202  !> Ionization fraction of H
203  !> H_ion_fr = H+/(H+ + H)
204  double precision, public, protected :: h_ion_fr=1d0
205  !> Ionization fraction of He
206  !> He_ion_fr = (He2+ + He+)/(He2+ + He+ + He)
207  double precision, public, protected :: he_ion_fr=1d0
208  !> Ratio of number He2+ / number He+ + He2+
209  !> He_ion_fr2 = He2+/(He2+ + He+)
210  double precision, public, protected :: he_ion_fr2=1d0
211  ! used for eq of state when it is not defined by units,
212  ! the units do not contain terms related to ionization fraction
213  ! and it is p = RR * rho * T
214  double precision, public, protected :: rr=1d0
215  ! remove the below flag and assume default value = .false.
216  ! when eq state properly implemented everywhere
217  ! and not anymore through units
218  logical, public, protected :: eq_state_units = .true.
219 
220  !> To control divB=0 fix for boundary
221  logical, public, protected :: boundary_divbfix(2*^nd)=.true.
222 
223  !> To skip * layer of ghost cells during divB=0 fix for boundary
224  integer, public, protected :: boundary_divbfix_skip(2*^nd)=0
225 
226  !> B0 field is force-free
227  logical, public, protected :: b0field_forcefree=.true.
228 
229  !> Whether an total energy equation is used
230  logical :: total_energy = .true.
231 
232  !> Whether unsplit semirelativistic MHD is solved
233  logical :: unsplit_semirelativistic=.false.
234 
235  !> Whether gravity work is included in energy equation
236  logical :: gravity_energy
237 
238  !> gamma minus one and its inverse
239  double precision :: gamma_1, inv_gamma_1
240 
241  !> inverse of squared speed of light c0 and reduced speed of light c
242  double precision :: inv_squared_c0, inv_squared_c
243 
244  ! DivB cleaning methods
245  integer, parameter :: divb_none = 0
246  integer, parameter :: divb_multigrid = -1
247  integer, parameter :: divb_glm = 1
248  integer, parameter :: divb_powel = 2
249  integer, parameter :: divb_janhunen = 3
250  integer, parameter :: divb_linde = 4
251  integer, parameter :: divb_lindejanhunen = 5
252  integer, parameter :: divb_lindepowel = 6
253  integer, parameter :: divb_lindeglm = 7
254  integer, parameter :: divb_ct = 8
255 
256  !define the subroutine interface for the ambipolar mask
257  abstract interface
258 
259  subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
261  integer, intent(in) :: ixi^l, ixo^l
262  double precision, intent(in) :: x(ixi^s,1:ndim)
263  double precision, intent(in) :: w(ixi^s,1:nw)
264  double precision, intent(inout) :: res(ixi^s)
265  end subroutine mask_subroutine
266 
267  function fun_kin_en(w, ixI^L, ixO^L, inv_rho) result(ke)
268  use mod_global_parameters, only: nw, ndim,block
269  integer, intent(in) :: ixI^L, ixO^L
270  double precision, intent(in) :: w(ixI^S, nw)
271  double precision :: ke(ixO^S)
272  double precision, intent(in), optional :: inv_rho(ixO^S)
273  end function fun_kin_en
274 
275  end interface
276 
277  procedure(mask_subroutine), pointer :: usr_mask_ambipolar => null()
278  procedure(sub_get_pthermal), pointer :: usr_rfactor => null()
279  procedure(sub_convert), pointer :: mhd_to_primitive => null()
280  procedure(sub_convert), pointer :: mhd_to_conserved => null()
281  procedure(sub_small_values), pointer :: mhd_handle_small_values => null()
282  procedure(sub_get_pthermal), pointer :: mhd_get_pthermal => null()
283  procedure(sub_get_v), pointer :: mhd_get_v => null()
284  procedure(fun_kin_en), pointer :: mhd_kin_en => null()
285  ! Public methods
286  public :: usr_mask_ambipolar
287  public :: usr_rfactor
288  public :: mhd_phys_init
289  public :: mhd_kin_en
290  public :: mhd_get_pthermal
291  public :: mhd_get_v
292  public :: mhd_get_rho
293  public :: mhd_get_v_idim
294  public :: mhd_to_conserved
295  public :: mhd_to_primitive
296  public :: mhd_get_csound2
297  public :: mhd_e_to_ei
298  public :: mhd_ei_to_e
299  public :: mhd_face_to_center
300  public :: get_divb
301  public :: get_current
302  !> needed public if we want to use the ambipolar coefficient in the user file
303  public :: multiplyambicoef
304  public :: get_normalized_divb
305  public :: b_from_vector_potential
306  public :: mhd_mag_en_all
307  {^nooned
308  public :: mhd_clean_divb_multigrid
309  }
310 
311 contains
312 
313  !> Read this module"s parameters from a file
314  subroutine mhd_read_params(files)
316  use mod_particles, only: particles_eta, particles_etah
317  character(len=*), intent(in) :: files(:)
318  integer :: n
319 
320  namelist /mhd_list/ mhd_energy, mhd_n_tracer, mhd_gamma, mhd_adiab,&
324  typedivbdiff, type_ct, compactres, divbwave, he_abundance, &
327  particles_eta, particles_etah,has_equi_rho0, has_equi_pe0,mhd_equi_thermal,&
331 
332  do n = 1, size(files)
333  open(unitpar, file=trim(files(n)), status="old")
334  read(unitpar, mhd_list, end=111)
335 111 close(unitpar)
336  end do
337 
338  end subroutine mhd_read_params
339 
340  !> Write this module's parameters to a snapsoht
341  subroutine mhd_write_info(fh)
343  integer, intent(in) :: fh
344  integer, parameter :: n_par = 1
345  double precision :: values(n_par)
346  character(len=name_len) :: names(n_par)
347  integer, dimension(MPI_STATUS_SIZE) :: st
348  integer :: er
349 
350  call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
351 
352  names(1) = "gamma"
353  values(1) = mhd_gamma
354  call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
355  call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
356  end subroutine mhd_write_info
357 
358  subroutine mhd_angmomfix(fC,x,wnew,ixI^L,ixO^L,idim)
360  double precision, intent(in) :: x(ixI^S,1:ndim)
361  double precision, intent(inout) :: fC(ixI^S,1:nwflux,1:ndim), wnew(ixI^S,1:nw)
362  integer, intent(in) :: ixI^L, ixO^L
363  integer, intent(in) :: idim
364  integer :: hxO^L, kxC^L, iw
365  double precision :: inv_volume(ixI^S)
366 
367  call mpistop("to do")
368 
369  end subroutine mhd_angmomfix
370 
371  subroutine mhd_phys_init()
375  use mod_viscosity, only: viscosity_init
376  use mod_gravity, only: gravity_init
377  use mod_particles, only: particles_init, particles_eta, particles_etah
381  use mod_cak_force, only: cak_init
382  {^nooned
384  }
385 
386  integer :: itr, idir
387 
388  call mhd_read_params(par_files)
389 
390  if(mhd_internal_e) then
391  if(mhd_solve_eaux) then
392  mhd_solve_eaux=.false.
393  if(mype==0) write(*,*) 'WARNING: set mhd_solve_eaux=F when mhd_internal_e=T'
394  end if
395  if(mhd_hydrodynamic_e) then
396  mhd_hydrodynamic_e=.false.
397  if(mype==0) write(*,*) 'WARNING: set mhd_hydrodynamic_e=F when mhd_internal_e=T'
398  end if
399  end if
400 
401  if(mhd_semirelativistic) then
402  unsplit_semirelativistic=.true.
403  if(mhd_internal_e) then
404  mhd_internal_e=.false.
405  if(mype==0) write(*,*) 'WARNING: set mhd_internal_e=F when mhd_semirelativistic=T'
406  end if
407  if(mhd_hydrodynamic_e) then
408  ! use split semirelativistic MHD
409  unsplit_semirelativistic=.false.
410  end if
411  if(mhd_boris_simplification) then
413  if(mype==0) write(*,*) 'WARNING: set mhd_boris_simplification=F when mhd_semirelativistic=T'
414  end if
415  end if
416 
417  if(.not. mhd_energy) then
418  if(mhd_internal_e) then
419  mhd_internal_e=.false.
420  if(mype==0) write(*,*) 'WARNING: set mhd_internal_e=F when mhd_energy=F'
421  end if
422  if(mhd_solve_eaux) then
423  mhd_solve_eaux=.false.
424  if(mype==0) write(*,*) 'WARNING: set mhd_solve_eaux=F when mhd_energy=F'
425  end if
426  if(mhd_hydrodynamic_e) then
427  mhd_hydrodynamic_e=.false.
428  if(mype==0) write(*,*) 'WARNING: set mhd_hydrodynamic_e=F when mhd_energy=F'
429  end if
430  if(mhd_thermal_conduction) then
431  mhd_thermal_conduction=.false.
432  if(mype==0) write(*,*) 'WARNING: set mhd_thermal_conduction=F when mhd_energy=F'
433  end if
434  if(mhd_radiative_cooling) then
435  mhd_radiative_cooling=.false.
436  if(mype==0) write(*,*) 'WARNING: set mhd_radiative_cooling=F when mhd_energy=F'
437  end if
438  if(mhd_trac) then
439  mhd_trac=.false.
440  if(mype==0) write(*,*) 'WARNING: set mhd_trac=F when mhd_energy=F'
441  end if
442  end if
443 
444  physics_type = "mhd"
445  phys_energy=mhd_energy
446  phys_internal_e=mhd_internal_e
447  phys_solve_eaux=mhd_solve_eaux
450 
451  phys_gamma = mhd_gamma
452 
453  if(mhd_energy.and..not.mhd_internal_e.and..not.mhd_hydrodynamic_e) then
454  total_energy=.true.
455  else
456  total_energy=.false.
457  end if
458  phys_total_energy=total_energy
459  if(mhd_energy) then
460  if(mhd_internal_e) then
461  gravity_energy=.false.
462  else
463  gravity_energy=.true.
464  end if
465  else
466  gravity_energy=.false.
467  end if
468 
469  {^ifoned
470  if(mhd_trac .and. mhd_trac_type .gt. 2) then
471  mhd_trac_type=1
472  if(mype==0) write(*,*) 'WARNING: reset mhd_trac_type=1 for 1D simulation'
473  end if
474  }
475  if(mhd_trac .and. mhd_trac_type .le. 4) then
476  mhd_trac_mask=bigdouble
477  if(mype==0) write(*,*) 'WARNING: set mhd_trac_mask==bigdouble for global TRAC method'
478  end if
480 
482 
483  ! set default gamma for polytropic/isothermal process
485  if(ndim==1) typedivbfix='none'
486  select case (typedivbfix)
487  case ('none')
488  type_divb = divb_none
489  {^nooned
490  case ('multigrid')
491  type_divb = divb_multigrid
492  use_multigrid = .true.
493  mg%operator_type = mg_laplacian
494  phys_global_source_after => mhd_clean_divb_multigrid
495  }
496  case ('glm')
497  mhd_glm = .true.
498  need_global_cmax = .true.
499  type_divb = divb_glm
500  case ('powel', 'powell')
501  type_divb = divb_powel
502  case ('janhunen')
503  type_divb = divb_janhunen
504  case ('linde')
505  type_divb = divb_linde
506  case ('lindejanhunen')
507  type_divb = divb_lindejanhunen
508  case ('lindepowel')
509  type_divb = divb_lindepowel
510  case ('lindeglm')
511  mhd_glm = .true.
512  need_global_cmax = .true.
513  type_divb = divb_lindeglm
514  case ('ct')
515  type_divb = divb_ct
516  stagger_grid = .true.
517  case default
518  call mpistop('Unknown divB fix')
519  end select
520 
521  allocate(start_indices(number_species),stop_indices(number_species))
522  ! set the index of the first flux variable for species 1
523  start_indices(1)=1
524  ! Determine flux variables
525  rho_ = var_set_rho()
526 
527  allocate(mom(ndir))
528  mom(:) = var_set_momentum(ndir)
529 
530  ! Set index of energy variable
531  if (mhd_energy) then
532  nwwave = 8
533  e_ = var_set_energy() ! energy density
534  p_ = e_ ! gas pressure
535  else
536  nwwave = 7
537  e_ = -1
538  p_ = -1
539  end if
540 
541  allocate(mag(ndir))
542  mag(:) = var_set_bfield(ndir)
543 
544  ! set auxiliary internal energy variable
545  if(mhd_energy .and. mhd_solve_eaux) then
546  eaux_ = var_set_internal_energy()
547  paux_ = eaux_
548  else
549  eaux_ = -1
550  paux_ = -1
551  end if
552 
553  if (mhd_glm) then
554  psi_ = var_set_fluxvar('psi', 'psi', need_bc=.false.)
555  else
556  psi_ = -1
557  end if
558 
559  allocate(tracer(mhd_n_tracer))
560 
561  ! Set starting index of tracers
562  do itr = 1, mhd_n_tracer
563  tracer(itr) = var_set_fluxvar("trc", "trp", itr, need_bc=.false.)
564  end do
565 
566  ! set number of variables which need update ghostcells
567  nwgc=nwflux
568 
569  ! set the index of the last flux variable for species 1
570  stop_indices(1)=nwflux
571 
572  ! set cutoff temperature when using the TRAC method, as well as an auxiliary weight
573  tweight_ = -1
574  if(mhd_trac) then
575  tcoff_ = var_set_wextra()
576  iw_tcoff=tcoff_
577  if(mhd_trac_type .ge. 3) then
578  tweight_ = var_set_wextra()
579  endif
580  else
581  tcoff_ = -1
582  end if
583 
584  ! set indices of equi vars and update number_equi_vars
585  number_equi_vars = 0
586  if(has_equi_rho0) then
589  iw_equi_rho = equi_rho0_
590  endif
591  if(has_equi_pe0) then
594  iw_equi_p = equi_pe0_
595  endif
596  ! determine number of stagger variables
597  if(stagger_grid) nws=ndim
598 
599  nvector = 2 ! No. vector vars
600  allocate(iw_vector(nvector))
601  iw_vector(1) = mom(1) - 1 ! TODO: why like this?
602  iw_vector(2) = mag(1) - 1 ! TODO: why like this?
603 
604  ! Check whether custom flux types have been defined
605  if (.not. allocated(flux_type)) then
606  allocate(flux_type(ndir, nwflux))
607  flux_type = flux_default
608  else if (any(shape(flux_type) /= [ndir, nwflux])) then
609  call mpistop("phys_check error: flux_type has wrong shape")
610  end if
611 
612  if(nwflux>mag(ndir)) then
613  ! for flux of eaux and tracers, using hll flux
614  flux_type(:,mag(ndir)+1:nwflux)=flux_hll
615  end if
616 
617  if(ndim>1) then
618  if(mhd_glm) then
619  flux_type(:,psi_)=flux_special
620  do idir=1,ndir
621  flux_type(idir,mag(idir))=flux_special
622  end do
623  else
624  do idir=1,ndir
625  flux_type(idir,mag(idir))=flux_tvdlf
626  end do
627  end if
628  end if
629 
630  phys_get_dt => mhd_get_dt
631  if(mhd_semirelativistic) then
632  phys_get_cmax => mhd_get_cmax_semirelati
633  else
634  phys_get_cmax => mhd_get_cmax_origin
635  end if
636  phys_get_a2max => mhd_get_a2max
637  phys_get_tcutoff => mhd_get_tcutoff
638  phys_get_h_speed => mhd_get_h_speed
639  if(has_equi_rho0) then
640  phys_get_cbounds => mhd_get_cbounds_split_rho
641  else if(mhd_semirelativistic) then
642  phys_get_cbounds => mhd_get_cbounds_semirelati
643  else
644  phys_get_cbounds => mhd_get_cbounds
645  end if
646  if(has_equi_rho0) then
647  phys_to_primitive => mhd_to_primitive_split_rho
649  phys_to_conserved => mhd_to_conserved_split_rho
651  else if(mhd_internal_e) then
652  phys_to_primitive => mhd_to_primitive_inte
654  phys_to_conserved => mhd_to_conserved_inte
656  else if(unsplit_semirelativistic) then
657  phys_to_primitive => mhd_to_primitive_semirelati
659  phys_to_conserved => mhd_to_conserved_semirelati
661  else if(mhd_hydrodynamic_e) then
662  phys_to_primitive => mhd_to_primitive_hde
664  phys_to_conserved => mhd_to_conserved_hde
666  else
667  phys_to_primitive => mhd_to_primitive_origin
669  phys_to_conserved => mhd_to_conserved_origin
671  end if
672  if(unsplit_semirelativistic) then
673  phys_get_flux => mhd_get_flux_semirelati
674  else
675  if(b0field.or.has_equi_rho0.or.has_equi_pe0) then
676  phys_get_flux => mhd_get_flux_split
677  else if(mhd_hydrodynamic_e) then
678  phys_get_flux => mhd_get_flux_hde
679  else
680  phys_get_flux => mhd_get_flux
681  end if
682  end if
683  if(mhd_boris_simplification) then
684  phys_get_v => mhd_get_v_boris
687  else
688  phys_get_v => mhd_get_v_origin
691  end if
692  if(b0field.or.has_equi_rho0) then
693  phys_add_source_geom => mhd_add_source_geom_split
694  else
695  phys_add_source_geom => mhd_add_source_geom
696  end if
697  phys_add_source => mhd_add_source
698  phys_check_params => mhd_check_params
699  phys_write_info => mhd_write_info
700  phys_angmomfix => mhd_angmomfix
701  if(unsplit_semirelativistic) then
702  phys_handle_small_values => mhd_handle_small_values_semirelati
703  mhd_handle_small_values => mhd_handle_small_values_semirelati
704  phys_check_w => mhd_check_w_semirelati
705  phys_get_pthermal => mhd_get_pthermal_semirelati
707  else if(mhd_hydrodynamic_e) then
708  phys_handle_small_values => mhd_handle_small_values_hde
709  mhd_handle_small_values => mhd_handle_small_values_hde
710  phys_check_w => mhd_check_w_hde
711  phys_get_pthermal => mhd_get_pthermal_hde
713  else
714  phys_handle_small_values => mhd_handle_small_values_origin
715  mhd_handle_small_values => mhd_handle_small_values_origin
716  phys_check_w => mhd_check_w_origin
717  phys_get_pthermal => mhd_get_pthermal_origin
719  end if
720  phys_energy_synchro => mhd_energy_synchro
721  if(number_equi_vars>0) then
722  phys_set_equi_vars => set_equi_vars_grid
723  endif
724 
725  if(type_divb==divb_glm) then
726  phys_modify_wlr => mhd_modify_wlr
727  end if
728 
729  ! if using ct stagger grid, boundary divb=0 is not done here
730  if(stagger_grid) then
731  phys_get_ct_velocity => mhd_get_ct_velocity
732  phys_update_faces => mhd_update_faces
733  phys_face_to_center => mhd_face_to_center
734  phys_modify_wlr => mhd_modify_wlr
735  else if(ndim>1) then
736  phys_boundary_adjust => mhd_boundary_adjust
737  end if
738 
739  {^nooned
740  ! clean initial divb
741  if(clean_initial_divb) phys_clean_divb => mhd_clean_divb_multigrid
742  }
743 
744  ! Whether diagonal ghost cells are required for the physics
745  if(type_divb < divb_linde) phys_req_diagonal = .false.
746 
747  ! derive units from basic units
748  call mhd_physical_units()
749 
750  if(.not. mhd_energy .and. mhd_thermal_conduction) then
751  call mpistop("thermal conduction needs mhd_energy=T")
752  end if
753  if(.not. mhd_energy .and. mhd_radiative_cooling) then
754  call mpistop("radiative cooling needs mhd_energy=T")
755  end if
756 
757  ! resistive MHD needs diagonal ghost cells
758  if(mhd_eta/=0.d0) phys_req_diagonal = .true.
759 
760  ! initialize thermal conduction module
761  if (mhd_thermal_conduction) then
762  phys_req_diagonal = .true.
763 
764  call sts_init()
766 
767  allocate(tc_fl)
770  if(phys_internal_e) then
771  if(has_equi_pe0 .and. has_equi_rho0) then
772  tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint_with_equi
773  else
774  tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint
775  end if
776  else if(mhd_hydrodynamic_e) then
777  tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_hde
778  else
779  if(has_equi_pe0 .and. has_equi_rho0) then
780  tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_etot_with_equi
781  else
782  tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_etot
783  end if
784  end if
785  if(mhd_solve_eaux) then
787  else if(unsplit_semirelativistic) then
789  else if(mhd_hydrodynamic_e) then
791  else if(.not. mhd_internal_e) then
793  end if
794  if(has_equi_pe0 .and. has_equi_rho0) then
795  tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint_with_equi
796  if(mhd_equi_thermal) then
797  tc_fl%has_equi = .true.
798  tc_fl%get_temperature_equi => mhd_get_temperature_equi
799  tc_fl%get_rho_equi => mhd_get_rho_equi
800  else
801  tc_fl%has_equi = .false.
802  endif
803  else
804  tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint
805  endif
807  tc_fl%get_rho => mhd_get_rho
808  tc_fl%e_ = e_
809  tc_fl%Tcoff_ = tcoff_
810  end if
811 
812  ! Initialize radiative cooling module
813  if (mhd_radiative_cooling) then
815  allocate(rc_fl)
817  rc_fl%get_rho => mhd_get_rho
818  rc_fl%get_pthermal => mhd_get_pthermal
819  rc_fl%e_ = e_
820  rc_fl%eaux_ = eaux_
821  rc_fl%Tcoff_ = tcoff_
822  if(associated(usr_rfactor)) then
823  rc_fl%get_var_Rfactor => usr_rfactor
824  endif
825  if(has_equi_pe0 .and. has_equi_rho0 .and. mhd_equi_thermal) then
826  rc_fl%has_equi = .true.
827  rc_fl%get_rho_equi => mhd_get_rho_equi
828  rc_fl%get_pthermal_equi => mhd_get_pe_equi
829  else
830  rc_fl%has_equi = .false.
831  end if
832  end if
833  allocate(te_fl_mhd)
834  te_fl_mhd%get_rho=> mhd_get_rho
835  te_fl_mhd%get_pthermal=> mhd_get_pthermal
836  te_fl_mhd%Rfactor = rr
837  if(associated(usr_rfactor)) then
838  te_fl_mhd%get_var_Rfactor => usr_rfactor
839  endif
840 {^ifthreed
841  phys_te_images => mhd_te_images
842 }
843  ! Initialize viscosity module
844  if (mhd_viscosity) call viscosity_init(phys_wider_stencil,phys_req_diagonal)
845 
846  ! Initialize gravity module
847  if(mhd_gravity) then
848  call gravity_init()
849  end if
850 
851  ! Initialize particles module
852  if(mhd_particles) then
853  call particles_init()
854  if (particles_eta < zero) particles_eta = mhd_eta
855  if (particles_etah < zero) particles_eta = mhd_etah
856  phys_req_diagonal = .true.
857  if(mype==0) then
858  write(*,*) '*****Using particles: with mhd_eta, mhd_etah :', mhd_eta, mhd_etah
859  write(*,*) '*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
860  end if
861  end if
862 
863  ! initialize magnetofriction module
864  if(mhd_magnetofriction) then
865  phys_req_diagonal = .true.
866  call magnetofriction_init()
867  end if
868 
869  ! For Hall, we need one more reconstructed layer since currents are computed
870  ! in mhd_get_flux: assuming one additional ghost layer (two for FOURTHORDER) was
871  ! added in nghostcells.
872  if(mhd_hall) then
873  phys_req_diagonal = .true.
874  if(mhd_4th_order) then
875  phys_wider_stencil = 2
876  else
877  phys_wider_stencil = 1
878  end if
879  end if
880 
881  if(mhd_ambipolar) then
882  phys_req_diagonal = .true.
883  if(mhd_ambipolar_sts) then
884  call sts_init()
885  if(mhd_internal_e) then
887  ndir,mag(1),ndir,.true.)
888  else
890  mag(ndir)-mom(ndir),mag(1),ndir,.true.)
891  end if
892  else
893  mhd_ambipolar_exp=.true.
894  ! For flux ambipolar term, we need one more reconstructed layer since currents are computed
895  ! in mhd_get_flux: assuming one additional ghost layer (two for FOURTHORDER) was
896  ! added in nghostcells.
897  if(mhd_4th_order) then
898  phys_wider_stencil = 2
899  else
900  phys_wider_stencil = 1
901  end if
902  end if
903  end if
904 
905  ! Initialize CAK radiation force module
906  if (mhd_cak_force) call cak_init(mhd_gamma)
907 
908  end subroutine mhd_phys_init
909 
910 {^ifthreed
911  subroutine mhd_te_images
914 
915  select case(convert_type)
916  case('EIvtiCCmpi','EIvtuCCmpi')
918  case('ESvtiCCmpi','ESvtuCCmpi')
920  case('SIvtiCCmpi','SIvtuCCmpi')
922  case default
923  call mpistop("Error in synthesize emission: Unknown convert_type")
924  end select
925  end subroutine mhd_te_images
926 }
927 
928 !!start th cond
929  ! wrappers for STS functions in thermal_conductivity module
930  ! which take as argument the tc_fluid (defined in the physics module)
931  subroutine mhd_sts_set_source_tc_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
933  use mod_fix_conserve
935  integer, intent(in) :: ixI^L, ixO^L, igrid, nflux
936  double precision, intent(in) :: x(ixI^S,1:ndim)
937  double precision, intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
938  double precision, intent(in) :: my_dt
939  logical, intent(in) :: fix_conserve_at_step
940  call sts_set_source_tc_mhd(ixi^l,ixo^l,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux,tc_fl)
941  end subroutine mhd_sts_set_source_tc_mhd
942 
943  function mhd_get_tc_dt_mhd(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
944  !Check diffusion time limit dt < dx_i**2/((gamma-1)*tc_k_para_i/rho)
945  !where tc_k_para_i=tc_k_para*B_i**2/B**2
946  !and T=p/rho
949 
950  integer, intent(in) :: ixi^l, ixo^l
951  double precision, intent(in) :: dx^d, x(ixi^s,1:ndim)
952  double precision, intent(in) :: w(ixi^s,1:nw)
953  double precision :: dtnew
954 
955  dtnew=get_tc_dt_mhd(w,ixi^l,ixo^l,dx^d,x,tc_fl)
956  end function mhd_get_tc_dt_mhd
957 
958  subroutine mhd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
960 
961  integer, intent(in) :: ixI^L,ixO^L
962  double precision, intent(inout) :: w(ixI^S,1:nw)
963  double precision, intent(in) :: x(ixI^S,1:ndim)
964  integer, intent(in) :: step
965  character(len=140) :: error_msg
966 
967  write(error_msg,"(a,i3)") "Thermal conduction step ", step
968  call mhd_handle_small_ei(w,x,ixi^l,ixo^l,e_,error_msg)
969  end subroutine mhd_tc_handle_small_e
970 
971  ! fill in tc_fluid fields from namelist
972  subroutine tc_params_read_mhd(fl)
974  type(tc_fluid), intent(inout) :: fl
975 
976  integer :: n
977 
978  ! list parameters
979  logical :: tc_perpendicular=.true.
980  logical :: tc_saturate=.false.
981  double precision :: tc_k_para=0d0
982  double precision :: tc_k_perp=0d0
983  character(len=std_len) :: tc_slope_limiter="MC"
984 
985  namelist /tc_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
986 
987  do n = 1, size(par_files)
988  open(unitpar, file=trim(par_files(n)), status="old")
989  read(unitpar, tc_list, end=111)
990 111 close(unitpar)
991  end do
992 
993  fl%tc_perpendicular = tc_perpendicular
994  fl%tc_saturate = tc_saturate
995  fl%tc_k_para = tc_k_para
996  fl%tc_k_perp = tc_k_perp
997  fl%tc_slope_limiter = tc_slope_limiter
998  end subroutine tc_params_read_mhd
999 !!end th cond
1000 
1001 !!rad cool
1002  subroutine rc_params_read(fl)
1004  use mod_constants, only: bigdouble
1005  type(rc_fluid), intent(inout) :: fl
1006  integer :: n
1007  ! list parameters
1008  integer :: ncool = 4000
1009  double precision :: cfrac=0.1d0
1010 
1011  !> Name of cooling curve
1012  character(len=std_len) :: coolcurve='JCcorona'
1013 
1014  !> Name of cooling method
1015  character(len=std_len) :: coolmethod='exact'
1016 
1017  !> Fixed temperature not lower than tlow
1018  logical :: Tfix=.false.
1019 
1020  !> Lower limit of temperature
1021  double precision :: tlow=bigdouble
1022 
1023  !> Add cooling source in a split way (.true.) or un-split way (.false.)
1024  logical :: rc_split=.false.
1025 
1026 
1027  namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
1028 
1029  do n = 1, size(par_files)
1030  open(unitpar, file=trim(par_files(n)), status="old")
1031  read(unitpar, rc_list, end=111)
1032 111 close(unitpar)
1033  end do
1034 
1035  fl%ncool=ncool
1036  fl%coolcurve=coolcurve
1037  fl%coolmethod=coolmethod
1038  fl%tlow=tlow
1039  fl%Tfix=tfix
1040  fl%rc_split=rc_split
1041  fl%cfrac=cfrac
1042 
1043  end subroutine rc_params_read
1044 !! end rad cool
1045 
1046  !> sets the equilibrium variables
1047  subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1049  use mod_usr_methods
1050  integer, intent(in) :: igrid, ixI^L, ixO^L
1051  double precision, intent(in) :: x(ixI^S,1:ndim)
1052 
1053  double precision :: delx(ixI^S,1:ndim)
1054  double precision :: xC(ixI^S,1:ndim),xshift^D
1055  integer :: idims, ixC^L, hxO^L, ix, idims2
1056 
1057  if(slab_uniform)then
1058  ^d&delx(ixi^s,^d)=rnode(rpdx^d_,igrid)\
1059  else
1060  ! for all non-cartesian and stretched cartesian coordinates
1061  delx(ixi^s,1:ndim)=ps(igrid)%dx(ixi^s,1:ndim)
1062  endif
1063 
1064  do idims=1,ndim
1065  hxo^l=ixo^l-kr(idims,^d);
1066  if(stagger_grid) then
1067  ! ct needs all transverse cells
1068  ixcmax^d=ixomax^d+nghostcells-nghostcells*kr(idims,^d); ixcmin^d=hxomin^d-nghostcells+nghostcells*kr(idims,^d);
1069  else
1070  ! ixC is centered index in the idims direction from ixOmin-1/2 to ixOmax+1/2
1071  ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
1072  end if
1073  ! always xshift=0 or 1/2
1074  xshift^d=half*(one-kr(^d,idims));
1075  do idims2=1,ndim
1076  select case(idims2)
1077  {case(^d)
1078  do ix = ixc^lim^d
1079  ! xshift=half: this is the cell center coordinate
1080  ! xshift=0: this is the cell edge i+1/2 coordinate
1081  xc(ix^d%ixC^s,^d)=x(ix^d%ixC^s,^d)+(half-xshift^d)*delx(ix^d%ixC^s,^d)
1082  end do\}
1083  end select
1084  end do
1085  call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1086  end do
1087 
1088  end subroutine set_equi_vars_grid_faces
1089 
1090  !> sets the equilibrium variables
1091  subroutine set_equi_vars_grid(igrid)
1093  use mod_usr_methods
1094 
1095  integer, intent(in) :: igrid
1096 
1097  !values at the center
1098  call usr_set_equi_vars(ixg^ll,ixg^ll,ps(igrid)%x,ps(igrid)%equi_vars(ixg^t,1:number_equi_vars,0))
1099 
1100  !values at the interfaces
1101  call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^ll,ixm^ll)
1102 
1103  end subroutine set_equi_vars_grid
1104 
1105  ! w, wnew conserved, add splitted variables back to wnew
1106  function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc) result(wnew)
1108  integer, intent(in) :: ixi^l,ixo^l, nwc
1109  double precision, intent(in) :: w(ixi^s, 1:nw)
1110  double precision, intent(in) :: x(ixi^s,1:ndim)
1111  double precision :: wnew(ixo^s, 1:nwc)
1112  double precision :: rho(ixi^s)
1113 
1114  call mhd_get_rho(w,x,ixi^l,ixo^l,rho(ixi^s))
1115  wnew(ixo^s,rho_) = rho(ixo^s)
1116  wnew(ixo^s,mom(:)) = w(ixo^s,mom(:))
1117 
1118  if (b0field) then
1119  ! add background magnetic field B0 to B
1120  wnew(ixo^s,mag(:))=w(ixo^s,mag(:))+block%B0(ixo^s,:,0)
1121  else
1122  wnew(ixo^s,mag(:))=w(ixo^s,mag(:))
1123  end if
1124 
1125  if(mhd_energy) then
1126  wnew(ixo^s,e_) = w(ixo^s,e_)
1127  if(has_equi_pe0) then
1128  wnew(ixo^s,e_) = wnew(ixo^s,e_) + block%equi_vars(ixo^s,equi_pe0_,0)* inv_gamma_1
1129  end if
1130  if(b0field .and. .not. mhd_internal_e) then
1131  wnew(ixo^s,e_)=wnew(ixo^s,e_)+0.5d0*sum(block%B0(ixo^s,:,0)**2,dim=ndim+1) &
1132  + sum(w(ixo^s,mag(:))*block%B0(ixo^s,:,0),dim=ndim+1)
1133  end if
1134  end if
1135 
1136  end function convert_vars_splitting
1137 
1138  subroutine mhd_check_params
1140  use mod_usr_methods
1141  use mod_convert, only: add_convert_method
1142 
1143  ! after user parameter setting
1144  gamma_1=mhd_gamma-1.d0
1145  if (.not. mhd_energy) then
1146  if (mhd_gamma <= 0.0d0) call mpistop ("Error: mhd_gamma <= 0")
1147  if (mhd_adiab < 0.0d0) call mpistop ("Error: mhd_adiab < 0")
1149  else
1150  if (mhd_gamma <= 0.0d0 .or. mhd_gamma == 1.0d0) &
1151  call mpistop ("Error: mhd_gamma <= 0 or mhd_gamma == 1")
1152  inv_gamma_1=1.d0/gamma_1
1153  small_e = small_pressure * inv_gamma_1
1154  end if
1155 
1156  if (number_equi_vars > 0 .and. .not. associated(usr_set_equi_vars)) then
1157  call mpistop("usr_set_equi_vars has to be implemented in the user file")
1158  endif
1159  if(convert .or. autoconvert) then
1160  if(convert_type .eq. 'dat_generic_mpi') then
1161  if(mhd_dump_full_vars) then
1162  if(mype .eq. 0) print*, " add conversion method: split -> full "
1163  call add_convert_method(convert_vars_splitting, nw, cons_wnames, "new")
1164  endif
1165  endif
1166  endif
1167  end subroutine mhd_check_params
1168 
1169  subroutine mhd_physical_units()
1171  double precision :: mp,kB,miu0,c_lightspeed
1172  double precision :: a,b
1173  ! Derive scaling units
1174  if(si_unit) then
1175  mp=mp_si
1176  kb=kb_si
1177  miu0=miu0_si
1178  c_lightspeed=c_si
1179  else
1180  mp=mp_cgs
1181  kb=kb_cgs
1182  miu0=4.d0*dpi ! G^2 cm^2 dyne^-1
1183  c_lightspeed=const_c
1184  end if
1185  if(eq_state_units) then
1186  a = 1d0 + 4d0 * he_abundance
1187  b = 1d0 + h_ion_fr + he_abundance*(he_ion_fr*(he_ion_fr2 + 1d0)+1d0)
1188  rr = 1d0
1189  else
1190  a = 1d0
1191  b = 1d0
1192  rr = (1d0 + h_ion_fr + he_abundance*(he_ion_fr*(he_ion_fr2 + 1d0)+1d0))/(1d0 + 4d0 * he_abundance)
1193  endif
1194  if(unit_density/=1.d0) then
1196  else
1197  ! unit of numberdensity is independent by default
1199  end if
1200  if(unit_velocity/=1.d0) then
1204  else if(unit_pressure/=1.d0) then
1208  else if(unit_magneticfield/=1.d0) then
1212  else
1213  ! unit of temperature is independent by default
1217  end if
1218  if(unit_time/=1.d0) then
1220  else
1221  ! unit of length is independent by default
1223  end if
1224  ! Additional units needed for the particles
1225  c_norm=c_lightspeed/unit_velocity
1227  if (.not. si_unit) unit_charge = unit_charge*const_c
1229 
1231  if(mhd_reduced_c<1.d0) then
1232  ! dimensionless speed
1233  inv_squared_c0=1.d0
1234  inv_squared_c=1.d0/mhd_reduced_c**2
1235  else
1236  inv_squared_c0=(unit_velocity/c_lightspeed)**2
1237  inv_squared_c=(unit_velocity/mhd_reduced_c)**2
1238  end if
1239  end if
1240 
1241  end subroutine mhd_physical_units
1242 
1243  subroutine mhd_check_w_semirelati(primitive,ixI^L,ixO^L,w,flag)
1245 
1246  logical, intent(in) :: primitive
1247  integer, intent(in) :: ixI^L, ixO^L
1248  double precision, intent(in) :: w(ixI^S,nw)
1249  double precision :: tmp(ixO^S),b2(ixO^S),b(ixO^S,1:ndir),Ba(ixO^S,1:ndir)
1250  double precision :: v(ixO^S,1:ndir),gamma2(ixO^S),inv_rho(ixO^S)
1251  logical, intent(inout) :: flag(ixI^S,1:nw)
1252 
1253  integer :: idir, jdir, kdir
1254 
1255  flag=.false.
1256  where(w(ixo^s,rho_) < small_density) flag(ixo^s,rho_) = .true.
1257 
1258  if(mhd_energy) then
1259  if(primitive) then
1260  where(w(ixo^s,p_) < small_pressure) flag(ixo^s,e_) = .true.
1261  else
1262  if(b0field) then
1263  ba(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+block%B0(ixo^s,1:ndir,b0i)
1264  else
1265  ba(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
1266  end if
1267  inv_rho(ixo^s) = 1d0/w(ixo^s,rho_)
1268  b2(ixo^s)=sum(ba(ixo^s,:)**2,dim=ndim+1)
1269  tmp(ixo^s)=sqrt(b2(ixo^s))
1270  where(tmp(ixo^s)>smalldouble)
1271  tmp(ixo^s)=1.d0/tmp(ixo^s)
1272  else where
1273  tmp(ixo^s)=0.d0
1274  end where
1275  do idir=1,ndir
1276  b(ixo^s,idir)=ba(ixo^s,idir)*tmp(ixo^s)
1277  end do
1278  tmp(ixo^s)=sum(b(ixo^s,:)*w(ixo^s,mom(:)),dim=ndim+1)
1279  ! Va^2/c^2
1280  b2(ixo^s)=b2(ixo^s)*inv_rho(ixo^s)*inv_squared_c
1281  ! equation (15)
1282  gamma2(ixo^s)=1.d0/(1.d0+b2(ixo^s))
1283  ! Convert momentum to velocity
1284  do idir = 1, ndir
1285  v(ixo^s,idir) = gamma2*(w(ixo^s, mom(idir))+b2*b(ixo^s,idir)*tmp)*inv_rho
1286  end do
1287  ! E=Bxv
1288  b=0.d0
1289  do idir=1,ndir; do jdir=1,ndir; do kdir=1,ndir
1290  if(lvc(idir,jdir,kdir)==1)then
1291  b(ixo^s,idir)=b(ixo^s,idir)+ba(ixo^s,jdir)*v(ixo^s,kdir)
1292  else if(lvc(idir,jdir,kdir)==-1)then
1293  b(ixo^s,idir)=b(ixo^s,idir)-ba(ixo^s,jdir)*v(ixo^s,kdir)
1294  end if
1295  end do; end do; end do
1296  ! Calculate internal e = e-eK-eB-eE
1297  tmp(ixo^s)=w(ixo^s,e_)&
1298  -half*(sum(v(ixo^s,:)**2,dim=ndim+1)*w(ixo^s,rho_)&
1299  +sum(w(ixo^s,mag(:))**2,dim=ndim+1)&
1300  +sum(b(ixo^s,:)**2,dim=ndim+1)*inv_squared_c)
1301  where(tmp(ixo^s) < small_e) flag(ixo^s,e_) = .true.
1302  end if
1303  end if
1304 
1305  end subroutine mhd_check_w_semirelati
1306 
1307  subroutine mhd_check_w_origin(primitive,ixI^L,ixO^L,w,flag)
1309 
1310  logical, intent(in) :: primitive
1311  integer, intent(in) :: ixI^L, ixO^L
1312  double precision, intent(in) :: w(ixI^S,nw)
1313  double precision :: tmp(ixI^S)
1314  logical, intent(inout) :: flag(ixI^S,1:nw)
1315 
1316  flag=.false.
1317  if(has_equi_rho0) then
1318  tmp(ixo^s) = w(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,0)
1319  else
1320  tmp(ixo^s) = w(ixo^s,rho_)
1321  endif
1322  where(tmp(ixo^s) < small_density) flag(ixo^s,rho_) = .true.
1323 
1324  if(mhd_energy) then
1325  if(primitive) then
1326  tmp(ixo^s) = w(ixo^s,e_)
1327  if(has_equi_pe0) then
1328  tmp(ixo^s) = tmp(ixo^s)+block%equi_vars(ixo^s,equi_pe0_,0)
1329  endif
1330  where(tmp(ixo^s) < small_pressure) flag(ixo^s,e_) = .true.
1331  else
1332  if(mhd_internal_e) then
1333  tmp(ixo^s)=w(ixo^s,e_)
1334  if(has_equi_pe0) then
1335  tmp(ixo^s) = tmp(ixo^s)+block%equi_vars(ixo^s,equi_pe0_,0)*inv_gamma_1
1336  endif
1337  where(tmp(ixo^s) < small_e) flag(ixo^s,e_) = .true.
1338  else
1339  tmp(ixo^s)=w(ixo^s,e_)-&
1340  mhd_kin_en(w,ixi^l,ixo^l)-mhd_mag_en(w,ixi^l,ixo^l)
1341  if(has_equi_pe0) then
1342  tmp(ixo^s) = tmp(ixo^s)+block%equi_vars(ixo^s,equi_pe0_,0)*inv_gamma_1
1343  endif
1344  where(tmp(ixo^s) < small_e) flag(ixo^s,e_) = .true.
1345  end if
1346  end if
1347  end if
1348 
1349  end subroutine mhd_check_w_origin
1350 
1351  subroutine mhd_check_w_hde(primitive,ixI^L,ixO^L,w,flag)
1353 
1354  logical, intent(in) :: primitive
1355  integer, intent(in) :: ixI^L, ixO^L
1356  double precision, intent(in) :: w(ixI^S,nw)
1357  double precision :: tmp(ixI^S)
1358  logical, intent(inout) :: flag(ixI^S,1:nw)
1359 
1360  flag=.false.
1361  where(w(ixo^s,rho_) < small_density) flag(ixo^s,rho_) = .true.
1362 
1363  if(mhd_energy) then
1364  if(primitive) then
1365  where(w(ixo^s,e_) < small_pressure) flag(ixo^s,e_) = .true.
1366  else
1367  tmp(ixo^s)=w(ixo^s,e_)-mhd_kin_en(w,ixi^l,ixo^l)
1368  where(tmp(ixo^s) < small_e) flag(ixo^s,e_) = .true.
1369  end if
1370  end if
1371 
1372  end subroutine mhd_check_w_hde
1373 
1374  !> Transform primitive variables into conservative ones
1375  subroutine mhd_to_conserved_origin(ixI^L,ixO^L,w,x)
1377  integer, intent(in) :: ixI^L, ixO^L
1378  double precision, intent(inout) :: w(ixI^S, nw)
1379  double precision, intent(in) :: x(ixI^S, 1:ndim)
1380 
1381  double precision :: inv_gamma2(ixO^S)
1382  integer :: idir
1383 
1384  !if (fix_small_values) then
1385  ! call mhd_handle_small_values(.true., w, x, ixI^L, ixO^L, 'mhd_to_conserved')
1386  !end if
1387 
1388  ! Calculate total energy from pressure, kinetic and magnetic energy
1389  if(mhd_energy) then
1390  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1&
1391  +half*sum(w(ixo^s,mom(:))**2,dim=ndim+1)*w(ixo^s,rho_)&
1392  +mhd_mag_en(w, ixi^l, ixo^l)
1393  if(mhd_solve_eaux) w(ixo^s,eaux_)=w(ixo^s,paux_)*inv_gamma_1
1394  end if
1395 
1396  if(mhd_boris_simplification) then
1397  inv_gamma2=1.d0+sum(w(ixo^s,mag(:))**2,dim=ndim+1)/w(ixo^s,rho_)*inv_squared_c
1398  ! Convert velocity to momentum
1399  do idir = 1, ndir
1400  w(ixo^s, mom(idir)) = inv_gamma2*w(ixo^s,rho_)*w(ixo^s, mom(idir))
1401  end do
1402  else
1403  ! Convert velocity to momentum
1404  do idir = 1, ndir
1405  w(ixo^s, mom(idir)) = w(ixo^s,rho_)*w(ixo^s, mom(idir))
1406  end do
1407  end if
1408  end subroutine mhd_to_conserved_origin
1409 
1410  !> Transform primitive variables into conservative ones
1411  subroutine mhd_to_conserved_hde(ixI^L,ixO^L,w,x)
1413  integer, intent(in) :: ixI^L, ixO^L
1414  double precision, intent(inout) :: w(ixI^S, nw)
1415  double precision, intent(in) :: x(ixI^S, 1:ndim)
1416 
1417  integer :: idir
1418 
1419  !if (fix_small_values) then
1420  ! call mhd_handle_small_values(.true., w, x, ixI^L, ixO^L, 'mhd_to_conserved')
1421  !end if
1422 
1423  ! Calculate total energy from pressure, kinetic and magnetic energy
1424  if(mhd_energy) then
1425  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1&
1426  +half*sum(w(ixo^s,mom(:))**2,dim=ndim+1)*w(ixo^s,rho_)
1427  end if
1428 
1429  ! Convert velocity to momentum
1430  do idir = 1, ndir
1431  w(ixo^s, mom(idir)) = w(ixo^s,rho_)*w(ixo^s, mom(idir))
1432  end do
1433  end subroutine mhd_to_conserved_hde
1434 
1435  !> Transform primitive variables into conservative ones
1436  subroutine mhd_to_conserved_inte(ixI^L,ixO^L,w,x)
1438  integer, intent(in) :: ixI^L, ixO^L
1439  double precision, intent(inout) :: w(ixI^S, nw)
1440  double precision, intent(in) :: x(ixI^S, 1:ndim)
1441 
1442  double precision :: inv_gamma2(ixO^S)
1443  integer :: idir
1444 
1445  !if (fix_small_values) then
1446  ! call mhd_handle_small_values(.true., w, x, ixI^L, ixO^L, 'mhd_to_conserved')
1447  !end if
1448 
1449  ! Calculate total energy from pressure, kinetic and magnetic energy
1450  if(mhd_energy) then
1451  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1
1452  end if
1453 
1454  if(mhd_boris_simplification) then
1455  inv_gamma2=1.d0+sum(w(ixo^s,mag(:))**2,dim=ndim+1)/w(ixo^s,rho_)*inv_squared_c
1456  ! Convert velocity to momentum
1457  do idir = 1, ndir
1458  w(ixo^s, mom(idir)) = inv_gamma2*w(ixo^s,rho_)*w(ixo^s, mom(idir))
1459  end do
1460  else
1461  ! Convert velocity to momentum
1462  do idir = 1, ndir
1463  w(ixo^s, mom(idir)) = w(ixo^s,rho_)*w(ixo^s, mom(idir))
1464  end do
1465  end if
1466  end subroutine mhd_to_conserved_inte
1467 
1468  !> Transform primitive variables into conservative ones
1469  subroutine mhd_to_conserved_split_rho(ixI^L,ixO^L,w,x)
1471  integer, intent(in) :: ixI^L, ixO^L
1472  double precision, intent(inout) :: w(ixI^S, nw)
1473  double precision, intent(in) :: x(ixI^S, 1:ndim)
1474  integer :: idir
1475  double precision :: rho(ixI^S)
1476 
1477  !if (fix_small_values) then
1478  ! call mhd_handle_small_values(.true., w, x, ixI^L, ixO^L, 'mhd_to_conserved')
1479  !end if
1480 
1481  rho(ixo^s) = w(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,b0i)
1482  ! Calculate total energy from pressure, kinetic and magnetic energy
1483  if(mhd_energy) then
1484  if(mhd_internal_e) then
1485  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1
1486  else
1487  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1&
1488  +half*sum(w(ixo^s,mom(:))**2,dim=ndim+1)*rho(ixo^s)&
1489  +mhd_mag_en(w, ixi^l, ixo^l)
1490  if(mhd_solve_eaux) w(ixo^s,eaux_)=w(ixo^s,paux_)*inv_gamma_1
1491  end if
1492  end if
1493 
1494  ! Convert velocity to momentum
1495  do idir = 1, ndir
1496  w(ixo^s, mom(idir)) = rho(ixo^s) * w(ixo^s, mom(idir))
1497  end do
1498  end subroutine mhd_to_conserved_split_rho
1499 
1500  !> Transform primitive variables into conservative ones
1501  subroutine mhd_to_conserved_semirelati(ixI^L,ixO^L,w,x)
1503  integer, intent(in) :: ixI^L, ixO^L
1504  double precision, intent(inout) :: w(ixI^S, nw)
1505  double precision, intent(in) :: x(ixI^S, 1:ndim)
1506 
1507  double precision :: E(ixO^S,1:ndir), B(ixO^S,1:ndir), S(ixO^S,1:ndir),tmp(ixO^S), b2(ixO^S)
1508  integer :: idir, jdir, kdir
1509 
1510  !if (fix_small_values) then
1511  ! call mhd_handle_small_values(.true., w, x, ixI^L, ixO^L, 'mhd_to_conserved')
1512  !end if
1513 
1514  if(b0field) then
1515  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+block%B0(ixo^s,1:ndir,b0i)
1516  else
1517  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
1518  end if
1519  ! Calculate total energy from pressure, kinetic and magnetic energy
1520  if(mhd_energy) then
1521  e=0.d0
1522  do idir=1,ndir; do jdir=1,ndir; do kdir=1,ndir
1523  if(lvc(idir,jdir,kdir)==1)then
1524  e(ixo^s,idir)=e(ixo^s,idir)+b(ixo^s,jdir)*w(ixo^s,mom(kdir))
1525  else if(lvc(idir,jdir,kdir)==-1)then
1526  e(ixo^s,idir)=e(ixo^s,idir)-b(ixo^s,jdir)*w(ixo^s,mom(kdir))
1527  end if
1528  end do; end do; end do
1529  ! equation (9)
1530  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1&
1531  +half*(sum(w(ixo^s,mom(:))**2,dim=ndim+1)*w(ixo^s,rho_)&
1532  +sum(w(ixo^s,mag(:))**2,dim=ndim+1)&
1533  +sum(e(ixo^s,:)**2,dim=ndim+1)*inv_squared_c)
1534  end if
1535 
1536  ! Convert velocity to momentum
1537  do idir = 1, ndir
1538  w(ixo^s, mom(idir)) = w(ixo^s,rho_) * w(ixo^s, mom(idir))
1539  end do
1540  !b2(ixO^S)=sum(w(ixO^S,mag(:))**2,dim=ndim+1)
1541  !tmp(ixO^S)=sqrt(b2(ixO^S))
1542  !where(tmp(ixO^S)>smalldouble)
1543  ! tmp(ixO^S)=1.d0/tmp(ixO^S)
1544  !else where
1545  ! tmp(ixO^S)=0.d0
1546  !end where
1547  !do idir=1,ndir
1548  ! b(ixO^S,idir)=w(ixO^S,mag(idir))*tmp(ixO^S)
1549  !end do
1550  !tmp(ixO^S)=sum(b(ixO^S,:)*w(ixO^S,mom(:)),dim=ndim+1)
1551  !do idir = 1, ndir
1552  ! w(ixO^S, mom(idir)) = w(ixO^S, mom(idir))+b2(ixO^S)/w(ixO^S,rho_)*inv_squared_c*(w(ixO^S, mom(idir))-b(ixO^S,idir)*tmp(ixO^S))
1553  !end do
1554  ! equation (5) Poynting vector
1555  s=0.d0
1556  do idir=1,ndir; do jdir=1,ndir; do kdir=1,ndir
1557  if(lvc(idir,jdir,kdir)==1)then
1558  s(ixo^s,idir)=s(ixo^s,idir)+e(ixo^s,jdir)*b(ixo^s,kdir)
1559  else if(lvc(idir,jdir,kdir)==-1)then
1560  s(ixo^s,idir)=s(ixo^s,idir)-e(ixo^s,jdir)*b(ixo^s,kdir)
1561  end if
1562  end do; end do; end do
1563  ! equation (9)
1564  do idir = 1, ndir
1565  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))+s(ixo^s,idir)*inv_squared_c
1566  end do
1567  end subroutine mhd_to_conserved_semirelati
1568 
1569  !> Transform conservative variables into primitive ones
1570  subroutine mhd_to_primitive_origin(ixI^L,ixO^L,w,x)
1572  integer, intent(in) :: ixI^L, ixO^L
1573  double precision, intent(inout) :: w(ixI^S, nw)
1574  double precision, intent(in) :: x(ixI^S, 1:ndim)
1575 
1576  double precision :: inv_rho(ixO^S), gamma2(ixO^S)
1577  integer :: idir
1578 
1579  if (fix_small_values) then
1580  ! fix small values preventing NaN numbers in the following converting
1581  call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive_origin')
1582  end if
1583 
1584  inv_rho(ixo^s) = 1d0/w(ixo^s,rho_)
1585 
1586  ! Calculate pressure = (gamma-1) * (e-ek-eb)
1587  if(mhd_energy) then
1588  w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)&
1589  -mhd_kin_en(w,ixi^l,ixo^l,inv_rho)&
1590  -mhd_mag_en(w,ixi^l,ixo^l))
1591  if(mhd_solve_eaux) w(ixo^s,paux_)=w(ixo^s,eaux_)*gamma_1
1592  end if
1593 
1594  if(mhd_boris_simplification) then
1595  gamma2=1.d0/(1.d0+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*inv_rho*inv_squared_c)
1596  ! Convert momentum to velocity
1597  do idir = 1, ndir
1598  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))*inv_rho*gamma2
1599  end do
1600  else
1601  ! Convert momentum to velocity
1602  do idir = 1, ndir
1603  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))*inv_rho
1604  end do
1605  end if
1606 
1607  end subroutine mhd_to_primitive_origin
1608 
1609  !> Transform conservative variables into primitive ones
1610  subroutine mhd_to_primitive_hde(ixI^L,ixO^L,w,x)
1612  integer, intent(in) :: ixI^L, ixO^L
1613  double precision, intent(inout) :: w(ixI^S, nw)
1614  double precision, intent(in) :: x(ixI^S, 1:ndim)
1615 
1616  double precision :: inv_rho(ixO^S)
1617  integer :: idir
1618 
1619  if (fix_small_values) then
1620  ! fix small values preventing NaN numbers in the following converting
1621  call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive_hde')
1622  end if
1623 
1624  inv_rho(ixo^s) = 1d0/w(ixo^s,rho_)
1625 
1626  ! Calculate pressure = (gamma-1) * (e-ek)
1627  if(mhd_energy) then
1628  w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)-mhd_kin_en(w,ixi^l,ixo^l,inv_rho))
1629  end if
1630 
1631  ! Convert momentum to velocity
1632  do idir = 1, ndir
1633  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))*inv_rho
1634  end do
1635 
1636  end subroutine mhd_to_primitive_hde
1637 
1638  !> Transform conservative variables into primitive ones
1639  subroutine mhd_to_primitive_inte(ixI^L,ixO^L,w,x)
1641  integer, intent(in) :: ixI^L, ixO^L
1642  double precision, intent(inout) :: w(ixI^S, nw)
1643  double precision, intent(in) :: x(ixI^S, 1:ndim)
1644 
1645  double precision :: inv_rho(ixO^S), gamma2(ixO^S)
1646  integer :: idir
1647 
1648  if (fix_small_values) then
1649  ! fix small values preventing NaN numbers in the following converting
1650  call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive_inte')
1651  end if
1652 
1653  inv_rho(ixo^s) = 1d0/w(ixo^s,rho_)
1654 
1655  ! Calculate pressure = (gamma-1) * e_internal
1656  if(mhd_energy) then
1657  w(ixo^s,p_)=w(ixo^s,e_)*gamma_1
1658  end if
1659 
1660  if(mhd_boris_simplification) then
1661  gamma2=1.d0/(1.d0+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*inv_rho*inv_squared_c)
1662  ! Convert momentum to velocity
1663  do idir = 1, ndir
1664  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))*inv_rho*gamma2
1665  end do
1666  else
1667  ! Convert momentum to velocity
1668  do idir = 1, ndir
1669  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))*inv_rho
1670  end do
1671  end if
1672 
1673  end subroutine mhd_to_primitive_inte
1674 
1675  !> Transform conservative variables into primitive ones
1676  subroutine mhd_to_primitive_split_rho(ixI^L,ixO^L,w,x)
1678  integer, intent(in) :: ixI^L, ixO^L
1679  double precision, intent(inout) :: w(ixI^S, nw)
1680  double precision, intent(in) :: x(ixI^S, 1:ndim)
1681  double precision :: inv_rho(ixO^S)
1682  integer :: idir
1683 
1684  if (fix_small_values) then
1685  ! fix small values preventing NaN numbers in the following converting
1686  call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive_split_rho')
1687  end if
1688 
1689  inv_rho(ixo^s) = 1d0/(w(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,b0i))
1690 
1691  ! Calculate pressure = (gamma-1) * (e-ek-eb)
1692  if(mhd_energy) then
1693  if(mhd_internal_e) then
1694  w(ixo^s,p_)=w(ixo^s,e_)*gamma_1
1695  else
1696  w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)&
1697  -mhd_kin_en(w,ixi^l,ixo^l,inv_rho)&
1698  -mhd_mag_en(w,ixi^l,ixo^l))
1699  if(mhd_solve_eaux) w(ixo^s,paux_)=w(ixo^s,eaux_)*gamma_1
1700  end if
1701  end if
1702 
1703  ! Convert momentum to velocity
1704  do idir = 1, ndir
1705  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))*inv_rho
1706  end do
1707 
1708  end subroutine mhd_to_primitive_split_rho
1709 
1710  !> Transform conservative variables into primitive ones
1711  subroutine mhd_to_primitive_semirelati(ixI^L,ixO^L,w,x)
1713  integer, intent(in) :: ixI^L, ixO^L
1714  double precision, intent(inout) :: w(ixI^S, nw)
1715  double precision, intent(in) :: x(ixI^S, 1:ndim)
1716 
1717  double precision :: inv_rho(ixO^S)
1718  double precision :: b(ixO^S,1:ndir), Ba(ixO^S,1:ndir),tmp(ixO^S), b2(ixO^S), gamma2(ixO^S)
1719  integer :: idir, jdir, kdir
1720 
1721  if (fix_small_values) then
1722  ! fix small values preventing NaN numbers in the following converting
1723  call mhd_handle_small_values_semirelati(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive_semirelati')
1724  end if
1725 
1726  if(b0field) then
1727  ba(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+block%B0(ixo^s,1:ndir,b0i)
1728  else
1729  ba(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
1730  end if
1731 
1732  inv_rho(ixo^s) = 1d0/w(ixo^s,rho_)
1733 
1734  b2(ixo^s)=sum(ba(ixo^s,:)**2,dim=ndim+1)
1735  tmp(ixo^s)=sqrt(b2(ixo^s))
1736  where(tmp(ixo^s)>smalldouble)
1737  tmp(ixo^s)=1.d0/tmp(ixo^s)
1738  else where
1739  tmp(ixo^s)=0.d0
1740  end where
1741  do idir=1,ndir
1742  b(ixo^s,idir)=ba(ixo^s,idir)*tmp(ixo^s)
1743  end do
1744  tmp(ixo^s)=sum(b(ixo^s,:)*w(ixo^s,mom(:)),dim=ndim+1)
1745 
1746  ! Va^2/c^2
1747  b2(ixo^s)=b2(ixo^s)*inv_rho(ixo^s)*inv_squared_c
1748  ! equation (15)
1749  gamma2(ixo^s)=1.d0/(1.d0+b2(ixo^s))
1750  ! Convert momentum to velocity
1751  do idir = 1, ndir
1752  w(ixo^s, mom(idir)) = gamma2*(w(ixo^s, mom(idir))+b2*b(ixo^s,idir)*tmp)*inv_rho
1753  end do
1754 
1755  if(mhd_energy) then
1756  ! E=Bxv
1757  b=0.d0
1758  do idir=1,ndir; do jdir=1,ndir; do kdir=1,ndir
1759  if(lvc(idir,jdir,kdir)==1)then
1760  b(ixo^s,idir)=b(ixo^s,idir)+ba(ixo^s,jdir)*w(ixo^s,mom(kdir))
1761  else if(lvc(idir,jdir,kdir)==-1)then
1762  b(ixo^s,idir)=b(ixo^s,idir)-ba(ixo^s,jdir)*w(ixo^s,mom(kdir))
1763  end if
1764  end do; end do; end do
1765  ! Calculate pressure = (gamma-1) * (e-eK-eB-eE)
1766  w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)&
1767  -half*(sum(w(ixo^s,mom(:))**2,dim=ndim+1)*w(ixo^s,rho_)&
1768  +sum(w(ixo^s,mag(:))**2,dim=ndim+1)&
1769  +sum(b(ixo^s,:)**2,dim=ndim+1)*inv_squared_c))
1770  end if
1771 
1772  end subroutine mhd_to_primitive_semirelati
1773 
1774  !> Transform internal energy to total energy
1775  subroutine mhd_ei_to_e(ixI^L,ixO^L,w,x)
1777  integer, intent(in) :: ixi^l, ixo^l
1778  double precision, intent(inout) :: w(ixi^s, nw)
1779  double precision, intent(in) :: x(ixi^s, 1:ndim)
1780 
1781  ! Calculate total energy from internal, kinetic and magnetic energy
1782  w(ixi^s,e_)=w(ixi^s,e_)&
1783  +mhd_kin_en(w,ixi^l,ixi^l)&
1784  +mhd_mag_en(w,ixi^l,ixi^l)
1785 
1786  end subroutine mhd_ei_to_e
1787 
1788  !> Transform internal energy to hydrodynamic energy
1789  subroutine mhd_ei_to_e_hde(ixI^L,ixO^L,w,x)
1791  integer, intent(in) :: ixI^L, ixO^L
1792  double precision, intent(inout) :: w(ixI^S, nw)
1793  double precision, intent(in) :: x(ixI^S, 1:ndim)
1794 
1795  ! Calculate hydrodynamic energy from internal and kinetic
1796  w(ixi^s,e_)=w(ixi^s,e_)+mhd_kin_en(w,ixi^l,ixi^l)
1797 
1798  end subroutine mhd_ei_to_e_hde
1799 
1800  !> Transform internal energy to total energy and velocity to momentum
1801  subroutine mhd_ei_to_e_semirelati(ixI^L,ixO^L,w,x)
1803  integer, intent(in) :: ixI^L, ixO^L
1804  double precision, intent(inout) :: w(ixI^S, nw)
1805  double precision, intent(in) :: x(ixI^S, 1:ndim)
1806 
1807  w(ixi^s,p_)=w(ixi^s,e_)*gamma_1
1808  call mhd_to_conserved_semirelati(ixi^l,ixi^l,w,x)
1809 
1810  end subroutine mhd_ei_to_e_semirelati
1811 
1812  !> Transform total energy to internal energy
1813  subroutine mhd_e_to_ei(ixI^L,ixO^L,w,x)
1815  integer, intent(in) :: ixi^l, ixo^l
1816  double precision, intent(inout) :: w(ixi^s, nw)
1817  double precision, intent(in) :: x(ixi^s, 1:ndim)
1818 
1819  ! Calculate ei = e - ek - eb
1820  w(ixi^s,e_)=w(ixi^s,e_)&
1821  -mhd_kin_en(w,ixi^l,ixi^l)&
1822  -mhd_mag_en(w,ixi^l,ixi^l)
1823 
1824  if(fix_small_values) then
1825  call mhd_handle_small_ei(w,x,ixi^l,ixi^l,e_,'mhd_e_to_ei')
1826  end if
1827 
1828  end subroutine mhd_e_to_ei
1829 
1830  !> Transform hydrodynamic energy to internal energy
1831  subroutine mhd_e_to_ei_hde(ixI^L,ixO^L,w,x)
1833  integer, intent(in) :: ixI^L, ixO^L
1834  double precision, intent(inout) :: w(ixI^S, nw)
1835  double precision, intent(in) :: x(ixI^S, 1:ndim)
1836 
1837  ! Calculate ei = e - ek
1838  w(ixi^s,e_)=w(ixi^s,e_)-mhd_kin_en(w,ixi^l,ixi^l)
1839 
1840  if(fix_small_values) then
1841  call mhd_handle_small_ei(w,x,ixi^l,ixi^l,e_,'mhd_e_to_ei_hde')
1842  end if
1843 
1844  end subroutine mhd_e_to_ei_hde
1845 
1846  !> Transform total energy to internal energy and momentum to velocity
1847  subroutine mhd_e_to_ei_semirelati(ixI^L,ixO^L,w,x)
1849  integer, intent(in) :: ixI^L, ixO^L
1850  double precision, intent(inout) :: w(ixI^S, nw)
1851  double precision, intent(in) :: x(ixI^S, 1:ndim)
1852 
1853  call mhd_to_primitive_semirelati(ixi^l,ixi^l,w,x)
1854  w(ixi^s,e_)=w(ixi^s,p_)*inv_gamma_1
1855 
1856  end subroutine mhd_e_to_ei_semirelati
1857 
1858  !> Update eaux and transform internal energy to total energy
1859  subroutine mhd_ei_to_e_aux(ixI^L,ixO^L,w,x)
1861  integer, intent(in) :: ixI^L, ixO^L
1862  double precision, intent(inout) :: w(ixI^S, nw)
1863  double precision, intent(in) :: x(ixI^S, 1:ndim)
1864 
1865  w(ixi^s,eaux_)=w(ixi^s,e_)
1866  ! Calculate total energy from internal, kinetic and magnetic energy
1867  w(ixi^s,e_)=w(ixi^s,e_)&
1868  +mhd_kin_en(w,ixi^l,ixi^l)&
1869  +mhd_mag_en(w,ixi^l,ixi^l)
1870 
1871  end subroutine mhd_ei_to_e_aux
1872 
1873  !> Transform total energy to internal energy via eaux as internal energy
1874  subroutine mhd_e_to_ei_aux(ixI^L,ixO^L,w,x)
1876  integer, intent(in) :: ixI^L, ixO^L
1877  double precision, intent(inout) :: w(ixI^S, nw)
1878  double precision, intent(in) :: x(ixI^S, 1:ndim)
1879 
1880  w(ixi^s,e_)=w(ixi^s,eaux_)
1881 
1882  end subroutine mhd_e_to_ei_aux
1883 
1884  subroutine mhd_energy_synchro(ixI^L,ixO^L,w,x)
1886  integer, intent(in) :: ixI^L,ixO^L
1887  double precision, intent(in) :: x(ixI^S,1:ndim)
1888  double precision, intent(inout) :: w(ixI^S,1:nw)
1889 
1890  double precision :: pth1(ixI^S),pth2(ixI^S),alfa(ixI^S),beta(ixI^S)
1891  double precision, parameter :: beta_low=0.005d0,beta_high=0.05d0
1892 
1893 ! double precision :: vtot(ixI^S),cs2(ixI^S),mach(ixI^S)
1894 ! double precision, parameter :: mach_low=20.d0,mach_high=200.d0
1895 
1896  ! get magnetic energy
1897  alfa(ixo^s)=mhd_mag_en(w,ixi^l,ixo^l)
1898  pth1(ixo^s)=gamma_1*(w(ixo^s,e_)-mhd_kin_en(w,ixi^l,ixo^l)-alfa(ixo^s))
1899  pth2(ixo^s)=w(ixo^s,eaux_)*gamma_1
1900  ! get plasma beta
1901  beta(ixo^s)=min(pth1(ixo^s),pth2(ixo^s))/alfa(ixo^s)
1902 
1903  ! whether Mach number should be another criterion ?
1904 ! vtot(ixO^S)=sum(w(ixO^S,mom(:))**2,dim=ndim+1)
1905 ! call mhd_get_csound2(w,x,ixI^L,ixO^L,cs2)
1906 ! mach(ixO^S)=sqrt(vtot(ixO^S)/cs2(ixO^S))/w(ixO^S,rho_)
1907  where(beta(ixo^s) .ge. beta_high)
1908 ! where(beta(ixO^S) .ge. beta_high .and. mach(ixO^S) .le. mach_low)
1909  w(ixo^s,eaux_)=pth1(ixo^s)*inv_gamma_1
1910  else where(beta(ixo^s) .le. beta_low)
1911 ! else where(beta(ixO^S) .le. beta_low .or. mach(ixO^S) .ge. mach_high)
1912  w(ixo^s,e_)=w(ixo^s,e_)-pth1(ixo^s)*inv_gamma_1+w(ixo^s,eaux_)
1913  else where
1914  alfa(ixo^s)=dlog(beta(ixo^s)/beta_low)/dlog(beta_high/beta_low)
1915 ! alfa(ixO^S)=min(dlog(beta(ixO^S)/beta_low)/dlog(beta_high/beta_low),
1916 ! dlog(mach_high(ixO^S)/mach(ixO^S))/dlog(mach_high/mach_low))
1917  w(ixo^s,eaux_)=(pth2(ixo^s)*(one-alfa(ixo^s))&
1918  +pth1(ixo^s)*alfa(ixo^s))*inv_gamma_1
1919  w(ixo^s,e_)=w(ixo^s,e_)-pth1(ixo^s)*inv_gamma_1+w(ixo^s,eaux_)
1920  end where
1921  end subroutine mhd_energy_synchro
1922 
1923  subroutine mhd_handle_small_values_semirelati(primitive, w, x, ixI^L, ixO^L, subname)
1925  use mod_small_values
1926  logical, intent(in) :: primitive
1927  integer, intent(in) :: ixI^L,ixO^L
1928  double precision, intent(inout) :: w(ixI^S,1:nw)
1929  double precision, intent(in) :: x(ixI^S,1:ndim)
1930  character(len=*), intent(in) :: subname
1931 
1932  double precision :: pressure(ixI^S), inv_rho(ixI^S), b2(ixI^S), tmp(ixI^S), gamma2(ixI^S)
1933  double precision :: b(ixI^S,1:ndir), v(ixI^S,1:ndir), Ba(ixI^S,1:ndir)
1934  integer :: idir, jdir, kdir, ix^D
1935  logical :: flag(ixI^S,1:nw)
1936 
1937  if(small_values_method == "ignore") return
1938 
1939  flag=.false.
1940  where(w(ixo^s,rho_) < small_density) flag(ixo^s,rho_) = .true.
1941 
1942  if(mhd_energy) then
1943  if(primitive) then
1944  where(w(ixo^s,p_) < small_pressure) flag(ixo^s,e_) = .true.
1945  else
1946  if(b0field) then
1947  ba(ixi^s,1:ndir)=w(ixi^s,mag(1:ndir))+block%B0(ixi^s,1:ndir,b0i)
1948  else
1949  ba(ixi^s,1:ndir)=w(ixi^s,mag(1:ndir))
1950  end if
1951  inv_rho(ixi^s) = 1d0/w(ixi^s,rho_)
1952  b2(ixi^s)=sum(ba(ixi^s,:)**2,dim=ndim+1)
1953  tmp(ixi^s)=sqrt(b2(ixi^s))
1954  where(tmp(ixi^s)>smalldouble)
1955  tmp(ixi^s)=1.d0/tmp(ixi^s)
1956  else where
1957  tmp(ixi^s)=0.d0
1958  end where
1959  do idir=1,ndir
1960  b(ixi^s,idir)=ba(ixi^s,idir)*tmp(ixi^s)
1961  end do
1962  tmp(ixi^s)=sum(b(ixi^s,:)*w(ixi^s,mom(:)),dim=ndim+1)
1963  ! Va^2/c^2
1964  b2(ixi^s)=b2(ixi^s)*inv_rho(ixi^s)*inv_squared_c
1965  ! equation (15)
1966  gamma2(ixi^s)=1.d0/(1.d0+b2(ixi^s))
1967  ! Convert momentum to velocity
1968  do idir = 1, ndir
1969  v(ixi^s,idir) = gamma2*(w(ixi^s, mom(idir))+b2*b(ixi^s,idir)*tmp(ixi^s))*inv_rho(ixi^s)
1970  end do
1971  ! E=Bxv
1972  b=0.d0
1973  do idir=1,ndir; do jdir=1,ndir; do kdir=1,ndir
1974  if(lvc(idir,jdir,kdir)==1)then
1975  b(ixi^s,idir)=b(ixi^s,idir)+ba(ixi^s,jdir)*v(ixi^s,kdir)
1976  else if(lvc(idir,jdir,kdir)==-1)then
1977  b(ixi^s,idir)=b(ixi^s,idir)-ba(ixi^s,jdir)*v(ixi^s,kdir)
1978  end if
1979  end do; end do; end do
1980  ! Calculate pressure p = (gamma-1) e-eK-eB-eE
1981  pressure(ixi^s)=gamma_1*(w(ixi^s,e_)&
1982  -half*(sum(v(ixi^s,:)**2,dim=ndim+1)*w(ixi^s,rho_)&
1983  +sum(w(ixi^s,mag(:))**2,dim=ndim+1)&
1984  +sum(b(ixi^s,:)**2,dim=ndim+1)*inv_squared_c))
1985  where(pressure(ixo^s) < small_pressure) flag(ixo^s,p_) = .true.
1986  end if
1987  end if
1988 
1989  if(any(flag)) then
1990  select case (small_values_method)
1991  case ("replace")
1992  where(flag(ixo^s,rho_)) w(ixo^s,rho_) = small_density
1993 
1994  if(mhd_energy) then
1995  if(primitive) then
1996  where(flag(ixo^s,e_)) w(ixo^s,p_) = small_pressure
1997  else
1998  {do ix^db=ixomin^db,ixomax^db\}
1999  if(flag(ix^d,e_)) then
2000  w(ix^d,e_)=small_pressure*inv_gamma_1+half*(sum(v(ix^d,:)**2)*w(ix^d,rho_)&
2001  +sum(w(ix^d,mag(:))**2)+sum(b(ix^d,:)**2)*inv_squared_c)
2002  end if
2003  {end do\}
2004  end if
2005  end if
2006  case ("average")
2007  ! do averaging of density
2008  call small_values_average(ixi^l, ixo^l, w, x, flag, rho_)
2009  if(mhd_energy) then
2010  if(primitive) then
2011  call small_values_average(ixi^l, ixo^l, w, x, flag, p_)
2012  else
2013  w(ixi^s,e_)=pressure(ixi^s)
2014  call small_values_average(ixi^l, ixo^l, w, x, flag, p_)
2015  w(ixi^s,e_)=w(ixi^s,p_)*inv_gamma_1&
2016  +half*(sum(v(ixi^s,:)**2,dim=ndim+1)*w(ixi^s,rho_)&
2017  +sum(w(ixi^s,mag(:))**2,dim=ndim+1)&
2018  +sum(b(ixi^s,:)**2,dim=ndim+1)*inv_squared_c)
2019  end if
2020  end if
2021  case default
2022  call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2023  end select
2024  end if
2025  end subroutine mhd_handle_small_values_semirelati
2026 
2027  subroutine mhd_handle_small_values_origin(primitive, w, x, ixI^L, ixO^L, subname)
2029  use mod_small_values
2030  logical, intent(in) :: primitive
2031  integer, intent(in) :: ixI^L,ixO^L
2032  double precision, intent(inout) :: w(ixI^S,1:nw)
2033  double precision, intent(in) :: x(ixI^S,1:ndim)
2034  character(len=*), intent(in) :: subname
2035 
2036  integer :: idir
2037  logical :: flag(ixI^S,1:nw)
2038  double precision :: tmp2(ixI^S)
2039 
2040  if(small_values_method == "ignore") return
2041 
2042  call phys_check_w(primitive, ixi^l, ixo^l, w, flag)
2043 
2044  if(any(flag)) then
2045  select case (small_values_method)
2046  case ("replace")
2047  if(has_equi_rho0) then
2048  where(flag(ixo^s,rho_)) w(ixo^s,rho_) = &
2049  small_density-block%equi_vars(ixo^s,equi_rho0_,0)
2050  else
2051  where(flag(ixo^s,rho_)) w(ixo^s,rho_) = small_density
2052  endif
2053  do idir = 1, ndir
2054  if(small_values_fix_iw(mom(idir))) then
2055  where(flag(ixo^s,rho_)) w(ixo^s, mom(idir)) = 0.0d0
2056  end if
2057  end do
2058 
2059  if(mhd_energy) then
2060  if(primitive) then
2061  if(has_equi_pe0) then
2062  tmp2(ixo^s) = small_pressure - &
2063  block%equi_vars(ixo^s,equi_pe0_,0)
2064  else
2065  tmp2(ixo^s) = small_pressure
2066  endif
2067  where(flag(ixo^s,e_)) w(ixo^s,p_) = tmp2(ixo^s)
2068  else
2069  ! conserved
2070  if(has_equi_pe0) then
2071  tmp2(ixo^s) = small_e - &
2072  block%equi_vars(ixo^s,equi_pe0_,0)*inv_gamma_1
2073  else
2074  tmp2(ixo^s) = small_e
2075  endif
2076  if(mhd_internal_e) then
2077  where(flag(ixo^s,e_))
2078  w(ixo^s,e_)=tmp2(ixo^s)
2079  end where
2080  else
2081  where(flag(ixo^s,e_))
2082  w(ixo^s,e_) = tmp2(ixo^s)+&
2083  mhd_kin_en(w,ixi^l,ixo^l)+&
2084  mhd_mag_en(w,ixi^l,ixo^l)
2085  end where
2086  if(mhd_solve_eaux) then
2087  where(flag(ixo^s,e_))
2088  w(ixo^s,eaux_)=tmp2(ixo^s)
2089  end where
2090  end if
2091  end if
2092  end if
2093  end if
2094  case ("average")
2095  if(primitive)then
2096  call small_values_average(ixi^l, ixo^l, w, x, flag, rho_)
2097  if(mhd_energy) then
2098  call small_values_average(ixi^l, ixo^l, w, x, flag, p_)
2099  end if
2100  else
2101  ! do averaging of density
2102  call small_values_average(ixi^l, ixo^l, w, x, flag, rho_)
2103  if(mhd_energy) then
2104  ! do averaging of internal energy
2105  if(.not.mhd_internal_e) then
2106  w(ixi^s,e_)=w(ixi^s,e_)&
2107  -mhd_kin_en(w,ixi^l,ixi^l)&
2108  -mhd_mag_en(w,ixi^l,ixi^l)
2109  end if
2110  call small_values_average(ixi^l, ixo^l, w, x, flag, e_)
2111  ! convert back
2112  if(.not.mhd_internal_e) then
2113  w(ixi^s,e_)=w(ixi^s,e_)&
2114  +mhd_kin_en(w,ixi^l,ixi^l)&
2115  +mhd_mag_en(w,ixi^l,ixi^l)
2116  end if
2117  ! eaux
2118  if(mhd_solve_eaux) then
2119  call small_values_average(ixi^l, ixo^l, w, x, flag, paux_)
2120  end if
2121  end if
2122  endif
2123  case default
2124  if(.not.primitive) then
2125  !convert w to primitive
2126  ! Calculate pressure = (gamma-1) * (e-ek-eb)
2127  if(mhd_energy) then
2128  if(mhd_internal_e) then
2129  w(ixo^s,p_)=w(ixo^s,e_)*gamma_1
2130  else
2131  w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)&
2132  -mhd_kin_en(w,ixi^l,ixo^l)&
2133  -mhd_mag_en(w,ixi^l,ixo^l))
2134  if(mhd_solve_eaux) w(ixo^s,paux_)=w(ixo^s,eaux_)*gamma_1
2135  end if
2136  end if
2137  ! Convert momentum to velocity
2138  if(has_equi_rho0) then
2139  tmp2(ixo^s) = w(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,0)
2140  else
2141  tmp2(ixo^s) = w(ixo^s,rho_)
2142  endif
2143  do idir = 1, ndir
2144  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))/tmp2(ixo^s)
2145  end do
2146  end if
2147  call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2148  end select
2149  end if
2150  end subroutine mhd_handle_small_values_origin
2151 
2152  subroutine mhd_handle_small_values_hde(primitive, w, x, ixI^L, ixO^L, subname)
2154  use mod_small_values
2155  logical, intent(in) :: primitive
2156  integer, intent(in) :: ixI^L,ixO^L
2157  double precision, intent(inout) :: w(ixI^S,1:nw)
2158  double precision, intent(in) :: x(ixI^S,1:ndim)
2159  character(len=*), intent(in) :: subname
2160 
2161  integer :: idir
2162  logical :: flag(ixI^S,1:nw)
2163  double precision :: tmp2(ixI^S)
2164 
2165  if(small_values_method == "ignore") return
2166 
2167  call phys_check_w(primitive, ixi^l, ixo^l, w, flag)
2168 
2169  if(any(flag)) then
2170  select case (small_values_method)
2171  case ("replace")
2172  where(flag(ixo^s,rho_)) w(ixo^s,rho_) = small_density
2173  do idir = 1, ndir
2174  if(small_values_fix_iw(mom(idir))) then
2175  where(flag(ixo^s,rho_)) w(ixo^s, mom(idir)) = 0.0d0
2176  end if
2177  end do
2178 
2179  if(mhd_energy) then
2180  where(flag(ixo^s,e_))
2181  w(ixo^s,e_) = small_e+mhd_kin_en(w,ixi^l,ixo^l)
2182  end where
2183  end if
2184  case ("average")
2185  ! do averaging of density
2186  call small_values_average(ixi^l, ixo^l, w, x, flag, rho_)
2187  if(mhd_energy) then
2188  if(primitive) then
2189  call small_values_average(ixi^l, ixo^l, w, x, flag, p_)
2190  else
2191  ! do averaging of internal energy
2192  w(ixi^s,e_)=w(ixi^s,e_)-mhd_kin_en(w,ixi^l,ixi^l)
2193  call small_values_average(ixi^l, ixo^l, w, x, flag, e_)
2194  ! convert back
2195  w(ixi^s,e_)=w(ixi^s,e_)+mhd_kin_en(w,ixi^l,ixi^l)
2196  end if
2197  end if
2198  case default
2199  if(.not.primitive) then
2200  !convert w to primitive
2201  ! Calculate pressure = (gamma-1) * (e-ek)
2202  if(mhd_energy) then
2203  w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)-mhd_kin_en(w,ixi^l,ixo^l))
2204  end if
2205  ! Convert momentum to velocity
2206  do idir = 1, ndir
2207  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))/w(ixo^s,rho_)
2208  end do
2209  end if
2210  call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2211  end select
2212  end if
2213 
2214  end subroutine mhd_handle_small_values_hde
2215 
2216  !> Calculate v vector
2217  subroutine mhd_get_v_origin(w,x,ixI^L,ixO^L,v)
2219 
2220  integer, intent(in) :: ixI^L, ixO^L
2221  double precision, intent(in) :: w(ixI^S,nw), x(ixI^S,1:ndim)
2222  double precision, intent(out) :: v(ixI^S,ndir)
2223 
2224  double precision :: rho(ixI^S)
2225  integer :: idir
2226 
2227  call mhd_get_rho(w,x,ixi^l,ixo^l,rho)
2228 
2229  rho(ixo^s)=1.d0/rho(ixo^s)
2230  ! Convert momentum to velocity
2231  do idir = 1, ndir
2232  v(ixo^s, idir) = w(ixo^s, mom(idir))*rho(ixo^s)
2233  end do
2234 
2235  end subroutine mhd_get_v_origin
2236 
2237  !> Calculate v vector
2238  subroutine mhd_get_v_boris(w,x,ixI^L,ixO^L,v)
2240 
2241  integer, intent(in) :: ixI^L, ixO^L
2242  double precision, intent(in) :: w(ixI^S,nw), x(ixI^S,1:ndim)
2243  double precision, intent(out) :: v(ixI^S,ndir)
2244 
2245  double precision :: rho(ixI^S), gamma2(ixO^S)
2246  integer :: idir
2247 
2248  call mhd_get_rho(w,x,ixi^l,ixo^l,rho)
2249 
2250  rho(ixo^s)=1.d0/rho(ixo^s)
2251  gamma2=1.d0/(1.d0+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*rho(ixo^s)*inv_squared_c)
2252  ! Convert momentum to velocity
2253  do idir = 1, ndir
2254  v(ixo^s, idir) = w(ixo^s, mom(idir))*rho(ixo^s)*gamma2
2255  end do
2256 
2257  end subroutine mhd_get_v_boris
2258 
2259  !> Calculate v component
2260  subroutine mhd_get_v_idim(w,x,ixI^L,ixO^L,idim,v)
2262 
2263  integer, intent(in) :: ixi^l, ixo^l, idim
2264  double precision, intent(in) :: w(ixi^s,nw), x(ixi^s,1:ndim)
2265  double precision, intent(out) :: v(ixi^s)
2266 
2267  double precision :: rho(ixi^s)
2268 
2269  call mhd_get_rho(w,x,ixi^l,ixo^l,rho)
2270 
2271  if(mhd_boris_simplification) then
2272  v(ixo^s) = w(ixo^s, mom(idim)) / rho(ixo^s) &
2273  /(1.d0+sum(w(ixo^s,mag(:))**2,dim=ndim+1)/rho(ixo^s)*inv_squared_c)
2274  else
2275  v(ixo^s) = w(ixo^s, mom(idim)) / rho(ixo^s)
2276  end if
2277 
2278  end subroutine mhd_get_v_idim
2279 
2280  !> Calculate cmax_idim=csound+abs(v_idim) within ixO^L
2281  subroutine mhd_get_cmax_origin(w,x,ixI^L,ixO^L,idim,cmax)
2283 
2284  integer, intent(in) :: ixI^L, ixO^L, idim
2285  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2286  double precision, intent(inout) :: cmax(ixI^S)
2287  double precision :: vel(ixI^S)
2288 
2289  call mhd_get_csound(w,x,ixi^l,ixo^l,idim,cmax)
2290  call mhd_get_v_idim(w,x,ixi^l,ixo^l,idim,vel)
2291 
2292  cmax(ixo^s)=abs(vel(ixo^s))+cmax(ixo^s)
2293 
2294  end subroutine mhd_get_cmax_origin
2295 
2296  !> Calculate cmax_idim for semirelativistic MHD
2297  subroutine mhd_get_cmax_semirelati(w,x,ixI^L,ixO^L,idim,cmax)
2299 
2300  integer, intent(in) :: ixI^L, ixO^L, idim
2301  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2302  double precision, intent(inout):: cmax(ixI^S)
2303  double precision :: wprim(ixI^S,nw)
2304  double precision :: csound(ixO^S), AvMinCs2(ixO^S), idim_Alfven_speed2(ixO^S)
2305  double precision :: inv_rho(ixO^S), Alfven_speed2(ixO^S), gamma2(ixO^S), B(ixO^S,1:ndir)
2306 
2307  if(b0field) then
2308  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+block%B0(ixo^s,1:ndir,b0i)
2309  else
2310  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
2311  end if
2312  inv_rho = 1.d0/w(ixo^s,rho_)
2313 
2314  alfven_speed2=sum(b(ixo^s,:)**2,dim=ndim+1)*inv_rho
2315  gamma2 = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2316 
2317  wprim=w
2318  call mhd_to_primitive(ixi^l,ixo^l,wprim,x)
2319  cmax(ixo^s)=1.d0-gamma2*wprim(ixo^s,mom(idim))**2*inv_squared_c
2320  ! equation (69)
2321  alfven_speed2=alfven_speed2*cmax(ixo^s)
2322 
2323  ! squared sound speed
2324  if(mhd_energy) then
2325  csound=mhd_gamma*wprim(ixo^s,p_)*inv_rho
2326  else
2327  csound=mhd_gamma*mhd_adiab*w(ixo^s,rho_)**gamma_1
2328  end if
2329 
2330  idim_alfven_speed2=b(ixo^s,idim)**2*inv_rho
2331 
2332  ! Va_hat^2+a_hat^2 equation (57)
2333  alfven_speed2=alfven_speed2+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2334 
2335  avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ixo^s)
2336 
2337  where(avmincs2<zero)
2338  avmincs2=zero
2339  end where
2340 
2341  ! equation (68) fast magnetosonic wave speed
2342  csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2343  cmax(ixo^s)=gamma2*abs(wprim(ixo^s,mom(idim)))+csound
2344 
2345  end subroutine mhd_get_cmax_semirelati
2346 
2347  subroutine mhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
2349 
2350  integer, intent(in) :: ixI^L, ixO^L
2351  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2352  double precision, intent(inout) :: a2max(ndim)
2353  double precision :: a2(ixI^S,ndim,nw)
2354  integer :: gxO^L,hxO^L,jxO^L,kxO^L,i,j
2355 
2356  a2=zero
2357  do i = 1,ndim
2358  !> 4th order
2359  hxo^l=ixo^l-kr(i,^d);
2360  gxo^l=hxo^l-kr(i,^d);
2361  jxo^l=ixo^l+kr(i,^d);
2362  kxo^l=jxo^l+kr(i,^d);
2363  a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
2364  -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
2365  a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/dxlevel(i)**2
2366  end do
2367  end subroutine mhd_get_a2max
2368 
2369  !> get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
2370  subroutine mhd_get_tcutoff(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
2372  use mod_geometry
2373  integer, intent(in) :: ixI^L,ixO^L
2374  double precision, intent(in) :: x(ixI^S,1:ndim)
2375  double precision, intent(inout) :: w(ixI^S,1:nw)
2376  double precision, intent(out) :: Tco_local,Tmax_local
2377 
2378  double precision, parameter :: trac_delta=0.25d0
2379  double precision :: tmp1(ixI^S),Te(ixI^S),lts(ixI^S)
2380  double precision, dimension(ixI^S,1:ndir) :: bunitvec
2381  double precision, dimension(ixI^S,1:ndim) :: gradT
2382  double precision :: Bdir(ndim)
2383  double precision :: ltrc,ltrp,altr(ixI^S)
2384  integer :: idims,jxO^L,hxO^L,ixA^D,ixB^D
2385  integer :: jxP^L,hxP^L,ixP^L
2386  logical :: lrlt(ixI^S)
2387 
2388  ! reuse lts as rhoc
2389  call mhd_get_rho(w,x,ixi^l,ixi^l,lts)
2390  if(mhd_internal_e) then
2391  tmp1(ixi^s)=w(ixi^s,e_)*gamma_1
2392  else
2393  call phys_get_pthermal(w,x,ixi^l,ixi^l,tmp1)
2394  end if
2395  te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)
2396  tco_local=zero
2397  tmax_local=maxval(te(ixo^s))
2398 
2399  {^ifoned
2400  select case(mhd_trac_type)
2401  case(0)
2402  !> test case, fixed cutoff temperature
2403  block%wextra(ixi^s,tcoff_)=2.5d5/unit_temperature
2404  case(1)
2405  hxo^l=ixo^l-1;
2406  jxo^l=ixo^l+1;
2407  lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
2408  lrlt=.false.
2409  where(lts(ixo^s) > trac_delta)
2410  lrlt(ixo^s)=.true.
2411  end where
2412  if(any(lrlt(ixo^s))) then
2413  tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
2414  end if
2415  case(2)
2416  !> iijima et al. 2021, LTRAC method
2417  ltrc=1.5d0
2418  ltrp=4.d0
2419  ixp^l=ixo^l^ladd1;
2420  hxo^l=ixo^l-1;
2421  jxo^l=ixo^l+1;
2422  hxp^l=ixp^l-1;
2423  jxp^l=ixp^l+1;
2424  lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
2425  lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2426  lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
2427  block%wextra(ixo^s,tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
2428  case default
2429  call mpistop("mhd_trac_type not allowed for 1D simulation")
2430  end select
2431  }
2432  {^nooned
2433  select case(mhd_trac_type)
2434  case(0)
2435  !> test case, fixed cutoff temperature
2436  block%wextra(ixi^s,tcoff_)=2.5d5/unit_temperature
2437  case(1,4,6)
2438  ! temperature gradient at cell centers
2439  do idims=1,ndim
2440  call gradient(te,ixi^l,ixo^l,idims,tmp1)
2441  gradt(ixo^s,idims)=tmp1(ixo^s)
2442  end do
2443  ! B vector
2444  if(b0field) then
2445  bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+block%B0(ixo^s,:,0)
2446  else
2447  bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
2448  end if
2449  if(mhd_trac_type .gt. 1) then
2450  ! B direction at cell center
2451  bdir=zero
2452  {do ixa^d=0,1\}
2453  ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
2454  bdir(1:ndim)=bdir(1:ndim)+bunitvec(ixb^d,1:ndim)
2455  {end do\}
2456  if(sum(bdir(:)**2) .gt. zero) then
2457  bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
2458  end if
2459  block%special_values(3:ndim+2)=bdir(1:ndim)
2460  end if
2461  tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
2462  where(tmp1(ixo^s)/=0.d0)
2463  tmp1(ixo^s)=1.d0/tmp1(ixo^s)
2464  elsewhere
2465  tmp1(ixo^s)=bigdouble
2466  end where
2467  ! b unit vector: magnetic field direction vector
2468  do idims=1,ndim
2469  bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
2470  end do
2471  ! temperature length scale inversed
2472  lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
2473  ! fraction of cells size to temperature length scale
2474  if(slab_uniform) then
2475  lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
2476  else
2477  lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
2478  end if
2479  lrlt=.false.
2480  where(lts(ixo^s) > trac_delta)
2481  lrlt(ixo^s)=.true.
2482  end where
2483  if(any(lrlt(ixo^s))) then
2484  block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
2485  else
2486  block%special_values(1)=zero
2487  end if
2488  block%special_values(2)=tmax_local
2489  case(2)
2490  !> iijima et al. 2021, LTRAC method
2491  ltrc=1.5d0
2492  ltrp=4.d0
2493  ixp^l=ixo^l^ladd1;
2494  ! temperature gradient at cell centers
2495  do idims=1,ndim
2496  call gradient(te,ixi^l,ixp^l,idims,tmp1)
2497  gradt(ixp^s,idims)=tmp1(ixp^s)
2498  end do
2499  ! B vector
2500  if(b0field) then
2501  bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))+block%B0(ixp^s,:,0)
2502  else
2503  bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))
2504  end if
2505  tmp1(ixp^s)=dsqrt(sum(bunitvec(ixp^s,:)**2,dim=ndim+1))
2506  where(tmp1(ixp^s)/=0.d0)
2507  tmp1(ixp^s)=1.d0/tmp1(ixp^s)
2508  elsewhere
2509  tmp1(ixp^s)=bigdouble
2510  end where
2511  ! b unit vector: magnetic field direction vector
2512  do idims=1,ndim
2513  bunitvec(ixp^s,idims)=bunitvec(ixp^s,idims)*tmp1(ixp^s)
2514  end do
2515  ! temperature length scale inversed
2516  lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
2517  ! fraction of cells size to temperature length scale
2518  if(slab_uniform) then
2519  lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
2520  else
2521  lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
2522  end if
2523  lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2524 
2525  altr=zero
2526  do idims=1,ndim
2527  hxo^l=ixo^l-kr(idims,^d);
2528  jxo^l=ixo^l+kr(idims,^d);
2529  altr(ixo^s)=altr(ixo^s)+0.25d0*(lts(hxo^s)+two*lts(ixo^s)+lts(jxo^s))*bunitvec(ixo^s,idims)**2
2530  end do
2531  block%wextra(ixo^s,tcoff_)=te(ixo^s)*altr(ixo^s)**0.4d0
2532  ! need one ghost layer for thermal conductivity
2533  {block%wextra(ixomin^d-1^d%ixP^s,tcoff_)=block%wextra(ixomin^d^d%ixP^s,tcoff_) \}
2534  {block%wextra(ixomax^d+1^d%ixP^s,tcoff_)=block%wextra(ixomax^d^d%ixP^s,tcoff_) \}
2535  case(3,5)
2536  !> do nothing here
2537  case default
2538  call mpistop("unknown mhd_trac_type")
2539  end select
2540  }
2541  end subroutine mhd_get_tcutoff
2542 
2543  !> get H speed for H-correction to fix the carbuncle problem at grid-aligned shock front
2544  subroutine mhd_get_h_speed(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2546 
2547  integer, intent(in) :: ixI^L, ixO^L, idim
2548  double precision, intent(in) :: wprim(ixI^S, nw)
2549  double precision, intent(in) :: x(ixI^S,1:ndim)
2550  double precision, intent(out) :: Hspeed(ixI^S,1:number_species)
2551 
2552  double precision :: csound(ixI^S,ndim),tmp(ixI^S)
2553  integer :: jxC^L, ixC^L, ixA^L, id, ix^D
2554 
2555  hspeed=0.d0
2556  ixa^l=ixo^l^ladd1;
2557  do id=1,ndim
2558  call mhd_get_csound_prim(wprim,x,ixi^l,ixa^l,id,tmp)
2559  csound(ixa^s,id)=tmp(ixa^s)
2560  end do
2561  ixcmax^d=ixomax^d;
2562  ixcmin^d=ixomin^d+kr(idim,^d)-1;
2563  jxcmax^d=ixcmax^d+kr(idim,^d);
2564  jxcmin^d=ixcmin^d+kr(idim,^d);
2565  hspeed(ixc^s,1)=0.5d0*abs(wprim(jxc^s,mom(idim))+csound(jxc^s,idim)-wprim(ixc^s,mom(idim))+csound(ixc^s,idim))
2566 
2567  do id=1,ndim
2568  if(id==idim) cycle
2569  ixamax^d=ixcmax^d+kr(id,^d);
2570  ixamin^d=ixcmin^d+kr(id,^d);
2571  hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixa^s,mom(id))+csound(ixa^s,id)-wprim(ixc^s,mom(id))+csound(ixc^s,id)))
2572  ixamax^d=ixcmax^d-kr(id,^d);
2573  ixamin^d=ixcmin^d-kr(id,^d);
2574  hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixc^s,mom(id))+csound(ixc^s,id)-wprim(ixa^s,mom(id))+csound(ixa^s,id)))
2575  end do
2576 
2577  do id=1,ndim
2578  if(id==idim) cycle
2579  ixamax^d=jxcmax^d+kr(id,^d);
2580  ixamin^d=jxcmin^d+kr(id,^d);
2581  hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixa^s,mom(id))+csound(ixa^s,id)-wprim(jxc^s,mom(id))+csound(jxc^s,id)))
2582  ixamax^d=jxcmax^d-kr(id,^d);
2583  ixamin^d=jxcmin^d-kr(id,^d);
2584  hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(jxc^s,mom(id))+csound(jxc^s,id)-wprim(ixa^s,mom(id))+csound(ixa^s,id)))
2585  end do
2586 
2587  end subroutine mhd_get_h_speed
2588 
2589  !> Estimating bounds for the minimum and maximum signal velocities without split
2590  subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2592 
2593  integer, intent(in) :: ixI^L, ixO^L, idim
2594  double precision, intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2595  double precision, intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2596  double precision, intent(in) :: x(ixI^S,1:ndim)
2597  double precision, intent(inout) :: cmax(ixI^S,1:number_species)
2598  double precision, intent(inout), optional :: cmin(ixI^S,1:number_species)
2599  double precision, intent(in) :: Hspeed(ixI^S,1:number_species)
2600 
2601  double precision :: wmean(ixI^S,nw)
2602  double precision, dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
2603  integer :: ix^D
2604 
2605  select case (boundspeed)
2606  case (1)
2607  ! This implements formula (10.52) from "Riemann Solvers and Numerical
2608  ! Methods for Fluid Dynamics" by Toro.
2609  tmp1(ixo^s)=sqrt(wlp(ixo^s,rho_))
2610  tmp2(ixo^s)=sqrt(wrp(ixo^s,rho_))
2611  tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2612  umean(ixo^s)=(wlp(ixo^s,mom(idim))*tmp1(ixo^s)+wrp(ixo^s,mom(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2613  call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2614  call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2615  dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2616  0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2617  (wrp(ixo^s,mom(idim))-wlp(ixo^s,mom(idim)))**2
2618  dmean(ixo^s)=sqrt(dmean(ixo^s))
2619  if(present(cmin)) then
2620  cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2621  cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2622  if(h_correction) then
2623  {do ix^db=ixomin^db,ixomax^db\}
2624  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2625  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2626  {end do\}
2627  end if
2628  else
2629  cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2630  end if
2631  case (2)
2632  wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2633  tmp1(ixo^s)=wmean(ixo^s,mom(idim))/wmean(ixo^s,rho_)
2634  call mhd_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
2635  if(present(cmin)) then
2636  cmax(ixo^s,1)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
2637  cmin(ixo^s,1)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
2638  if(h_correction) then
2639  {do ix^db=ixomin^db,ixomax^db\}
2640  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2641  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2642  {end do\}
2643  end if
2644  else
2645  cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2646  end if
2647  case (3)
2648  ! Miyoshi 2005 JCP 208, 315 equation (67)
2649  call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2650  call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2651  csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2652  if(present(cmin)) then
2653  cmin(ixo^s,1)=min(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))-csoundl(ixo^s)
2654  cmax(ixo^s,1)=max(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
2655  if(h_correction) then
2656  {do ix^db=ixomin^db,ixomax^db\}
2657  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2658  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2659  {end do\}
2660  end if
2661  else
2662  cmax(ixo^s,1)=max(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
2663  end if
2664  end select
2665 
2666  end subroutine mhd_get_cbounds
2667 
2668  !> Estimating bounds for the minimum and maximum signal velocities without split
2669  subroutine mhd_get_cbounds_semirelati(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2671 
2672  integer, intent(in) :: ixI^L, ixO^L, idim
2673  double precision, intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2674  double precision, intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2675  double precision, intent(in) :: x(ixI^S,1:ndim)
2676  double precision, intent(inout) :: cmax(ixI^S,1:number_species)
2677  double precision, intent(inout), optional :: cmin(ixI^S,1:number_species)
2678  double precision, intent(in) :: Hspeed(ixI^S,1:number_species)
2679 
2680  double precision, dimension(ixO^S) :: csoundL, csoundR, gamma2L, gamma2R
2681 
2682  ! Miyoshi 2005 JCP 208, 315 equation (67)
2683  call mhd_get_csound_semirelati(wlp,x,ixi^l,ixo^l,idim,csoundl,gamma2l)
2684  call mhd_get_csound_semirelati(wrp,x,ixi^l,ixo^l,idim,csoundr,gamma2r)
2685  csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2686  if(present(cmin)) then
2687  cmin(ixo^s,1)=min(gamma2l*wlp(ixo^s,mom(idim)),gamma2r*wrp(ixo^s,mom(idim)))-csoundl(ixo^s)
2688  cmax(ixo^s,1)=max(gamma2l*wlp(ixo^s,mom(idim)),gamma2r*wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
2689  else
2690  cmax(ixo^s,1)=max(gamma2l*wlp(ixo^s,mom(idim)),gamma2r*wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
2691  end if
2692 
2693  end subroutine mhd_get_cbounds_semirelati
2694 
2695  !> Estimating bounds for the minimum and maximum signal velocities with rho split
2696  subroutine mhd_get_cbounds_split_rho(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2698 
2699  integer, intent(in) :: ixI^L, ixO^L, idim
2700  double precision, intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2701  double precision, intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2702  double precision, intent(in) :: x(ixI^S,1:ndim)
2703  double precision, intent(inout) :: cmax(ixI^S,1:number_species)
2704  double precision, intent(inout), optional :: cmin(ixI^S,1:number_species)
2705  double precision, intent(in) :: Hspeed(ixI^S,1:number_species)
2706 
2707  double precision :: wmean(ixI^S,nw)
2708  double precision, dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
2709  integer :: ix^D
2710  double precision :: rho(ixI^S)
2711 
2712  select case (boundspeed)
2713  case (1)
2714  ! This implements formula (10.52) from "Riemann Solvers and Numerical
2715  ! Methods for Fluid Dynamics" by Toro.
2716  tmp1(ixo^s)=sqrt(wlp(ixo^s,rho_)+block%equi_vars(ixo^s,equi_rho0_,b0i))
2717  tmp2(ixo^s)=sqrt(wrp(ixo^s,rho_)+block%equi_vars(ixo^s,equi_rho0_,b0i))
2718  tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2719  umean(ixo^s)=(wlp(ixo^s,mom(idim))*tmp1(ixo^s)+wrp(ixo^s,mom(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2720  call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2721  call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2722  dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2723  0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2724  (wrp(ixo^s,mom(idim))-wlp(ixo^s,mom(idim)))**2
2725  dmean(ixo^s)=sqrt(dmean(ixo^s))
2726  if(present(cmin)) then
2727  cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2728  cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2729  if(h_correction) then
2730  {do ix^db=ixomin^db,ixomax^db\}
2731  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2732  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2733  {end do\}
2734  end if
2735  else
2736  cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2737  end if
2738  case (2)
2739  wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2740  tmp1(ixo^s)=wmean(ixo^s,mom(idim))/(wmean(ixo^s,rho_)+block%equi_vars(ixo^s,equi_rho0_,b0i))
2741  call mhd_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
2742  if(present(cmin)) then
2743  cmax(ixo^s,1)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
2744  cmin(ixo^s,1)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
2745  if(h_correction) then
2746  {do ix^db=ixomin^db,ixomax^db\}
2747  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2748  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2749  {end do\}
2750  end if
2751  else
2752  cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2753  end if
2754  case (3)
2755  ! Miyoshi 2005 JCP 208, 315 equation (67)
2756  call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2757  call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2758  csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2759  if(present(cmin)) then
2760  cmin(ixo^s,1)=min(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))-csoundl(ixo^s)
2761  cmax(ixo^s,1)=max(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
2762  if(h_correction) then
2763  {do ix^db=ixomin^db,ixomax^db\}
2764  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2765  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2766  {end do\}
2767  end if
2768  else
2769  cmax(ixo^s,1)=max(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
2770  end if
2771  end select
2772 
2773  end subroutine mhd_get_cbounds_split_rho
2774 
2775  !> prepare velocities for ct methods
2776  subroutine mhd_get_ct_velocity(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2778 
2779  integer, intent(in) :: ixI^L, ixO^L, idim
2780  double precision, intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2781  double precision, intent(in) :: cmax(ixI^S)
2782  double precision, intent(in), optional :: cmin(ixI^S)
2783  type(ct_velocity), intent(inout):: vcts
2784 
2785  integer :: idimE,idimN
2786 
2787  ! calculate velocities related to different UCT schemes
2788  select case(type_ct)
2789  case('average')
2790  case('uct_contact')
2791  if(.not.allocated(vcts%vnorm)) allocate(vcts%vnorm(ixi^s,1:ndim))
2792  ! get average normal velocity at cell faces
2793  vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,mom(idim))+wrp(ixo^s,mom(idim)))
2794  case('uct_hll')
2795  if(.not.allocated(vcts%vbarC)) then
2796  allocate(vcts%vbarC(ixi^s,1:ndir,2),vcts%vbarLC(ixi^s,1:ndir,2),vcts%vbarRC(ixi^s,1:ndir,2))
2797  allocate(vcts%cbarmin(ixi^s,1:ndim),vcts%cbarmax(ixi^s,1:ndim))
2798  end if
2799  ! Store magnitude of characteristics
2800  if(present(cmin)) then
2801  vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
2802  vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2803  else
2804  vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2805  vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
2806  end if
2807 
2808  idimn=mod(idim,ndir)+1 ! 'Next' direction
2809  idime=mod(idim+1,ndir)+1 ! Electric field direction
2810  ! Store velocities
2811  vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,mom(idimn))
2812  vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,mom(idimn))
2813  vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
2814  +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2815  /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2816 
2817  vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,mom(idime))
2818  vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,mom(idime))
2819  vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
2820  +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2821  /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2822  case default
2823  call mpistop('choose average, uct_contact,or uct_hll for type_ct!')
2824  end select
2825 
2826  end subroutine mhd_get_ct_velocity
2827 
2828  !> Calculate fast magnetosonic wave speed
2829  subroutine mhd_get_csound(w,x,ixI^L,ixO^L,idim,csound)
2831 
2832  integer, intent(in) :: ixI^L, ixO^L, idim
2833  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2834  double precision, intent(out):: csound(ixI^S)
2835  double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2836  double precision :: inv_rho(ixO^S)
2837 
2838  if(has_equi_rho0) then
2839  inv_rho(ixo^s) = 1d0/(w(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,b0i))
2840  else
2841  inv_rho(ixo^s) = 1d0/w(ixo^s,rho_)
2842  endif
2843 
2844  call mhd_get_csound2(w,x,ixi^l,ixo^l,csound)
2845 
2846  ! store |B|^2 in v
2847  b2(ixo^s) = mhd_mag_en_all(w,ixi^l,ixo^l)
2848 
2849  cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
2850  avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2851  * mhd_mag_i_all(w,ixi^l,ixo^l,idim)**2 * inv_rho
2852 
2853  where(avmincs2(ixo^s)<zero)
2854  avmincs2(ixo^s)=zero
2855  end where
2856 
2857  avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2858 
2859  if (.not. mhd_hall) then
2860  csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2861  if (mhd_boris_simplification) then
2862  ! equation (88)
2863  csound(ixo^s) = mhd_gamma_alfven(w, ixi^l,ixo^l) * csound(ixo^s)
2864  end if
2865  else
2866  ! take the Hall velocity into account:
2867  ! most simple estimate, high k limit:
2868  ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
2869  kmax = dpi/min({dxlevel(^d)},bigdouble)*half
2870  csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2871  mhd_etah * sqrt(b2(ixo^s))*inv_rho*kmax)
2872  end if
2873 
2874  end subroutine mhd_get_csound
2875 
2876  !> Calculate fast magnetosonic wave speed
2877  subroutine mhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
2879 
2880  integer, intent(in) :: ixI^L, ixO^L, idim
2881  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2882  double precision, intent(out):: csound(ixI^S)
2883  double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2884  double precision :: inv_rho(ixO^S)
2885  double precision :: tmp(ixI^S)
2886 
2887  call mhd_get_rho(w,x,ixi^l,ixo^l,tmp)
2888  inv_rho(ixo^s) = 1d0/tmp(ixo^s)
2889 
2890 
2891  if(mhd_energy) then
2892  if(has_equi_pe0) then
2893  csound(ixo^s) = w(ixo^s,e_) + block%equi_vars(ixo^s,equi_pe0_,b0i)
2894  else
2895  csound(ixo^s) = w(ixo^s,e_)
2896  endif
2897  csound(ixo^s)=mhd_gamma*csound(ixo^s)*inv_rho
2898  else
2899  csound(ixo^s)=mhd_gamma*mhd_adiab*tmp(ixo^s)**gamma_1
2900  end if
2901 
2902  ! store |B|^2 in v
2903  b2(ixo^s) = mhd_mag_en_all(w,ixi^l,ixo^l)
2904  cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
2905  avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2906  * mhd_mag_i_all(w,ixi^l,ixo^l,idim)**2 * inv_rho
2907 
2908  where(avmincs2(ixo^s)<zero)
2909  avmincs2(ixo^s)=zero
2910  end where
2911 
2912  avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2913 
2914  if (.not. mhd_hall) then
2915  csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2916  if (mhd_boris_simplification) then
2917  ! equation (88)
2918  csound(ixo^s) = mhd_gamma_alfven(w, ixi^l,ixo^l) * csound(ixo^s)
2919  end if
2920  else
2921  ! take the Hall velocity into account:
2922  ! most simple estimate, high k limit:
2923  ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
2924  kmax = dpi/min({dxlevel(^d)},bigdouble)*half
2925  csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2926  mhd_etah * sqrt(b2(ixo^s))*inv_rho*kmax)
2927  end if
2928 
2929  end subroutine mhd_get_csound_prim
2930 
2931  !> Calculate cmax_idim for semirelativistic MHD
2932  subroutine mhd_get_csound_semirelati(w,x,ixI^L,ixO^L,idim,csound,gamma2)
2934 
2935  integer, intent(in) :: ixI^L, ixO^L, idim
2936  ! here w is primitive variables
2937  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2938  double precision, intent(out):: csound(ixO^S), gamma2(ixO^S)
2939  double precision :: AvMinCs2(ixO^S), kmax
2940  double precision :: inv_rho(ixO^S), Alfven_speed2(ixO^S), idim_Alfven_speed2(ixO^S),B(ixO^S,1:ndir)
2941 
2942  if(b0field) then
2943  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+block%B0(ixo^s,1:ndir,b0i)
2944  else
2945  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
2946  end if
2947 
2948  inv_rho = 1.d0/w(ixo^s,rho_)
2949 
2950  alfven_speed2=sum(b(ixo^s,:)**2,dim=ndim+1)*inv_rho
2951  gamma2 = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2952 
2953  avmincs2=1.d0-gamma2*w(ixo^s,mom(idim))**2*inv_squared_c
2954  ! equatoin (69)
2955  alfven_speed2=alfven_speed2*avmincs2
2956 
2957  ! squared sound speed
2958  if(mhd_energy) then
2959  csound(ixo^s)=mhd_gamma*w(ixo^s,p_)*inv_rho
2960  else
2961  csound(ixo^s)=mhd_gamma*mhd_adiab*w(ixo^s,rho_)**gamma_1
2962  end if
2963 
2964  idim_alfven_speed2=b(ixo^s,idim)**2*inv_rho
2965 
2966  ! Va_hat^2+a_hat^2 equation (57)
2967  alfven_speed2=alfven_speed2+csound(ixo^s)*(1.d0+idim_alfven_speed2*inv_squared_c)
2968 
2969  avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound(ixo^s)*idim_alfven_speed2*avmincs2
2970 
2971  where(avmincs2<zero)
2972  avmincs2=zero
2973  end where
2974 
2975  ! equation (68) fast magnetosonic speed
2976  csound(ixo^s) = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2977 
2978  end subroutine mhd_get_csound_semirelati
2979 
2980  !> Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho-b**2/2) within ixO^L
2981  subroutine mhd_get_pthermal_origin(w,x,ixI^L,ixO^L,pth)
2984 
2985  integer, intent(in) :: ixI^L, ixO^L
2986  double precision, intent(in) :: w(ixI^S,nw)
2987  double precision, intent(in) :: x(ixI^S,1:ndim)
2988  double precision, intent(out):: pth(ixI^S)
2989  integer :: iw, ix^D
2990 
2991  if(mhd_energy) then
2992  if(mhd_internal_e) then
2993  pth(ixo^s)=gamma_1*w(ixo^s,e_)
2994  else
2995  pth(ixo^s)=gamma_1*(w(ixo^s,e_)&
2996  - mhd_kin_en(w,ixi^l,ixo^l)&
2997  - mhd_mag_en(w,ixi^l,ixo^l))
2998  end if
2999  if(has_equi_pe0) then
3000  pth(ixo^s) = pth(ixo^s) + block%equi_vars(ixo^s,equi_pe0_,b0i)
3001  endif
3002  else
3003  call mhd_get_rho(w,x,ixi^l,ixo^l,pth)
3004  pth(ixo^s)=mhd_adiab*pth(ixo^s)**mhd_gamma
3005  end if
3006 
3007  if (fix_small_values) then
3008  {do ix^db= ixo^lim^db\}
3009  if(pth(ix^d)<small_pressure) then
3010  pth(ix^d)=small_pressure
3011  end if
3012  {enddo^d&\}
3013  end if
3014 
3015  if (check_small_values) then
3016  {do ix^db= ixo^lim^db\}
3017  if(pth(ix^d)<small_pressure) then
3018  write(*,*) "Error: small value of gas pressure",pth(ix^d),&
3019  " encountered when call mhd_get_pthermal"
3020  write(*,*) "Iteration: ", it, " Time: ", global_time
3021  write(*,*) "Location: ", x(ix^d,:)
3022  write(*,*) "Cell number: ", ix^d
3023  do iw=1,nw
3024  write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
3025  end do
3026  ! use erroneous arithmetic operation to crash the run
3027  if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
3028  write(*,*) "Saving status at the previous time step"
3029  crash=.true.
3030  end if
3031  {enddo^d&\}
3032  end if
3033 
3034  end subroutine mhd_get_pthermal_origin
3035 
3036  !> Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho-b**2/2) within ixO^L
3037  subroutine mhd_get_pthermal_semirelati(w,x,ixI^L,ixO^L,pth)
3040 
3041  integer, intent(in) :: ixI^L, ixO^L
3042  double precision, intent(in) :: w(ixI^S,nw)
3043  double precision, intent(in) :: x(ixI^S,1:ndim)
3044  double precision, intent(out):: pth(ixI^S)
3045  integer :: iw, ix^D
3046 
3047  double precision :: wprim(ixI^S,nw)
3048 
3049  if(mhd_energy) then
3050  wprim=w
3051  call mhd_to_primitive_semirelati(ixi^l,ixo^l,wprim,x)
3052  pth(ixo^s)=wprim(ixo^s,p_)
3053  else
3054  pth(ixo^s)=mhd_adiab*w(ixo^s,rho_)**mhd_gamma
3055  end if
3056 
3057  if (check_small_values) then
3058  {do ix^db= ixo^lim^db\}
3059  if(pth(ix^d)<small_pressure) then
3060  write(*,*) "Error: small value of gas pressure",pth(ix^d),&
3061  " encountered when call mhd_get_pthermal_semirelati"
3062  write(*,*) "Iteration: ", it, " Time: ", global_time
3063  write(*,*) "Location: ", x(ix^d,:)
3064  write(*,*) "Cell number: ", ix^d
3065  do iw=1,nw
3066  write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
3067  end do
3068  ! use erroneous arithmetic operation to crash the run
3069  if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
3070  write(*,*) "Saving status at the previous time step"
3071  crash=.true.
3072  end if
3073  {enddo^d&\}
3074  end if
3075 
3076  end subroutine mhd_get_pthermal_semirelati
3077 
3078  !> Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho-b**2/2) within ixO^L
3079  subroutine mhd_get_pthermal_hde(w,x,ixI^L,ixO^L,pth)
3082 
3083  integer, intent(in) :: ixI^L, ixO^L
3084  double precision, intent(in) :: w(ixI^S,nw)
3085  double precision, intent(in) :: x(ixI^S,1:ndim)
3086  double precision, intent(out):: pth(ixI^S)
3087  integer :: iw, ix^D
3088 
3089  if(mhd_energy) then
3090  pth(ixo^s)=gamma_1*(w(ixo^s,e_)-mhd_kin_en(w,ixi^l,ixo^l))
3091  else
3092  pth(ixo^s)=mhd_adiab*w(ixo^s,rho_)**mhd_gamma
3093  end if
3094 
3095  if (fix_small_values) then
3096  {do ix^db= ixo^lim^db\}
3097  if(pth(ix^d)<small_pressure) then
3098  pth(ix^d)=small_pressure
3099  end if
3100  {enddo^d&\}
3101  end if
3102 
3103  if (check_small_values) then
3104  {do ix^db= ixo^lim^db\}
3105  if(pth(ix^d)<small_pressure) then
3106  write(*,*) "Error: small value of gas pressure",pth(ix^d),&
3107  " encountered when call mhd_get_pthermal_hde"
3108  write(*,*) "Iteration: ", it, " Time: ", global_time
3109  write(*,*) "Location: ", x(ix^d,:)
3110  write(*,*) "Cell number: ", ix^d
3111  do iw=1,nw
3112  write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
3113  end do
3114  ! use erroneous arithmetic operation to crash the run
3115  if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
3116  write(*,*) "Saving status at the previous time step"
3117  crash=.true.
3118  end if
3119  {enddo^d&\}
3120  end if
3121 
3122  end subroutine mhd_get_pthermal_hde
3123 
3124  !> Calculate temperature=p/rho when in e_ the internal energy is stored
3125  subroutine mhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
3127  integer, intent(in) :: ixI^L, ixO^L
3128  double precision, intent(in) :: w(ixI^S, 1:nw)
3129  double precision, intent(in) :: x(ixI^S, 1:ndim)
3130  double precision, intent(out):: res(ixI^S)
3131  res(ixo^s) = gamma_1 * w(ixo^s, e_) /w(ixo^s,rho_)
3132  end subroutine mhd_get_temperature_from_eint
3133 
3134  !> Calculate temperature=p/rho when in e_ the total energy is stored
3135  !> this does not check the values of mhd_energy and mhd_internal_e,
3136  !> mhd_energy = .true. and mhd_internal_e = .false.
3137  !> also check small_values is avoided
3138  subroutine mhd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
3140  integer, intent(in) :: ixI^L, ixO^L
3141  double precision, intent(in) :: w(ixI^S, 1:nw)
3142  double precision, intent(in) :: x(ixI^S, 1:ndim)
3143  double precision, intent(out):: res(ixI^S)
3144  res(ixo^s)=(gamma_1*(w(ixo^s,e_)&
3145  - mhd_kin_en(w,ixi^l,ixo^l)&
3146  - mhd_mag_en(w,ixi^l,ixo^l)))/w(ixo^s,rho_)
3147  end subroutine mhd_get_temperature_from_etot
3148 
3149  !> Calculate temperature from hydrodynamic energy
3150  subroutine mhd_get_temperature_from_hde(w, x, ixI^L, ixO^L, res)
3152  integer, intent(in) :: ixI^L, ixO^L
3153  double precision, intent(in) :: w(ixI^S, 1:nw)
3154  double precision, intent(in) :: x(ixI^S, 1:ndim)
3155  double precision, intent(out):: res(ixI^S)
3156  res(ixo^s)=gamma_1*(w(ixo^s,e_)&
3157  - mhd_kin_en(w,ixi^l,ixo^l))/w(ixo^s,rho_)
3158  end subroutine mhd_get_temperature_from_hde
3159 
3160  subroutine mhd_get_temperature_from_eint_with_equi(w, x, ixI^L, ixO^L, res)
3162  integer, intent(in) :: ixI^L, ixO^L
3163  double precision, intent(in) :: w(ixI^S, 1:nw)
3164  double precision, intent(in) :: x(ixI^S, 1:ndim)
3165  double precision, intent(out):: res(ixI^S)
3166  res(ixo^s) = (gamma_1 * w(ixo^s, e_) + block%equi_vars(ixo^s,equi_pe0_,b0i)) /&
3167  (w(ixo^s,rho_) +block%equi_vars(ixo^s,equi_rho0_,b0i))
3169 
3170  subroutine mhd_get_temperature_equi(w,x, ixI^L, ixO^L, res)
3172  integer, intent(in) :: ixI^L, ixO^L
3173  double precision, intent(in) :: w(ixI^S, 1:nw)
3174  double precision, intent(in) :: x(ixI^S, 1:ndim)
3175  double precision, intent(out):: res(ixI^S)
3176  res(ixo^s)= block%equi_vars(ixo^s,equi_pe0_,b0i)/block%equi_vars(ixo^s,equi_rho0_,b0i)
3177  end subroutine mhd_get_temperature_equi
3178 
3179  subroutine mhd_get_rho_equi(w, x, ixI^L, ixO^L, res)
3181  integer, intent(in) :: ixI^L, ixO^L
3182  double precision, intent(in) :: w(ixI^S, 1:nw)
3183  double precision, intent(in) :: x(ixI^S, 1:ndim)
3184  double precision, intent(out):: res(ixI^S)
3185  res(ixo^s) = block%equi_vars(ixo^s,equi_rho0_,b0i)
3186  end subroutine mhd_get_rho_equi
3187 
3188  subroutine mhd_get_pe_equi(w,x, ixI^L, ixO^L, res)
3190  integer, intent(in) :: ixI^L, ixO^L
3191  double precision, intent(in) :: w(ixI^S, 1:nw)
3192  double precision, intent(in) :: x(ixI^S, 1:ndim)
3193  double precision, intent(out):: res(ixI^S)
3194  res(ixo^s) = block%equi_vars(ixo^s,equi_pe0_,b0i)
3195  end subroutine mhd_get_pe_equi
3196 
3197  subroutine mhd_get_temperature_from_etot_with_equi(w, x, ixI^L, ixO^L, res)
3199  integer, intent(in) :: ixI^L, ixO^L
3200  double precision, intent(in) :: w(ixI^S, 1:nw)
3201  double precision, intent(in) :: x(ixI^S, 1:ndim)
3202  double precision, intent(out):: res(ixI^S)
3203  res(ixo^s)=(gamma_1*(w(ixo^s,e_)&
3204  - mhd_kin_en(w,ixi^l,ixo^l)&
3205  - mhd_mag_en(w,ixi^l,ixo^l)) + block%equi_vars(ixo^s,equi_pe0_,b0i))&
3206  /(w(ixo^s,rho_) +block%equi_vars(ixo^s,equi_rho0_,b0i))
3207 
3209 
3210  !> Calculate the square of the thermal sound speed csound2 within ixO^L.
3211  !> csound2=gamma*p/rho
3212  subroutine mhd_get_csound2(w,x,ixI^L,ixO^L,csound2)
3214  integer, intent(in) :: ixi^l, ixo^l
3215  double precision, intent(in) :: w(ixi^s,nw)
3216  double precision, intent(in) :: x(ixi^s,1:ndim)
3217  double precision, intent(out) :: csound2(ixi^s)
3218  double precision :: rho(ixi^s)
3219 
3220  call mhd_get_rho(w,x,ixi^l,ixo^l,rho)
3221  if(mhd_energy) then
3222  call mhd_get_pthermal(w,x,ixi^l,ixo^l,csound2)
3223  csound2(ixo^s)=mhd_gamma*csound2(ixo^s)/rho(ixo^s)
3224  else
3225  csound2(ixo^s)=mhd_gamma*mhd_adiab*rho(ixo^s)**gamma_1
3226  end if
3227  end subroutine mhd_get_csound2
3228 
3229  !> Calculate total pressure within ixO^L including magnetic pressure
3230  subroutine mhd_get_p_total(w,x,ixI^L,ixO^L,p)
3232 
3233  integer, intent(in) :: ixI^L, ixO^L
3234  double precision, intent(in) :: w(ixI^S,nw)
3235  double precision, intent(in) :: x(ixI^S,1:ndim)
3236  double precision, intent(out) :: p(ixI^S)
3237 
3238  call mhd_get_pthermal(w,x,ixi^l,ixo^l,p)
3239 
3240  p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
3241 
3242  end subroutine mhd_get_p_total
3243 
3244  !> Calculate fluxes within ixO^L without any splitting
3245  subroutine mhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3247  use mod_geometry
3248 
3249  integer, intent(in) :: ixI^L, ixO^L, idim
3250  ! conservative w
3251  double precision, intent(in) :: wC(ixI^S,nw)
3252  ! primitive w
3253  double precision, intent(in) :: w(ixI^S,nw)
3254  double precision, intent(in) :: x(ixI^S,1:ndim)
3255  double precision,intent(out) :: f(ixI^S,nwflux)
3256 
3257  double precision :: ptotal(ixO^S)
3258  double precision :: tmp(ixI^S)
3259  double precision :: vHall(ixI^S,1:ndir)
3260  integer :: idirmin, iw, idir, jdir, kdir
3261  double precision, allocatable, dimension(:^D&,:) :: Jambi, btot
3262  double precision, allocatable, dimension(:^D&) :: tmp2, tmp3
3263 
3264  ! Get flux of density
3265  f(ixo^s,rho_)=w(ixo^s,mom(idim))*w(ixo^s,rho_)
3266 
3267  if(mhd_energy) then
3268  ptotal(ixo^s)=w(ixo^s,p_)+0.5d0*sum(w(ixo^s,mag(:))**2,dim=ndim+1)
3269  else
3270  ptotal(ixo^s)=mhd_adiab*w(ixo^s,rho_)**mhd_gamma+0.5d0*sum(w(ixo^s,mag(:))**2,dim=ndim+1)
3271  end if
3272 
3273  if (mhd_hall) then
3274  call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3275  end if
3276 
3277  ! Get flux of tracer
3278  do iw=1,mhd_n_tracer
3279  f(ixo^s,tracer(iw))=w(ixo^s,mom(idim))*w(ixo^s,tracer(iw))
3280  end do
3281 
3282  ! Get flux of momentum
3283  ! f_i[m_k]=v_i*m_k-b_k*b_i [+ptotal if i==k]
3284  do idir=1,ndir
3285  if(idim==idir) then
3286  f(ixo^s,mom(idir))=wc(ixo^s,mom(idir))*w(ixo^s,mom(idim))+ptotal(ixo^s)-&
3287  w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3288  else
3289  f(ixo^s,mom(idir))=wc(ixo^s,mom(idir))*w(ixo^s,mom(idim))-&
3290  w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3291  end if
3292  end do
3293 
3294  ! Get flux of energy
3295  ! f_i[e]=v_i*e+v_i*ptotal-b_i*(b_k*v_k)
3296  if(mhd_energy) then
3297  if (mhd_internal_e) then
3298  f(ixo^s,e_)=w(ixo^s,mom(idim))*wc(ixo^s,e_)
3299  else
3300  f(ixo^s,e_)=w(ixo^s,mom(idim))*(wc(ixo^s,e_)+ptotal(ixo^s))&
3301  -w(ixo^s,mag(idim))*sum(w(ixo^s,mag(:))*w(ixo^s,mom(:)),dim=ndim+1)
3302  if(mhd_solve_eaux) f(ixo^s,eaux_)=w(ixo^s,mom(idim))*wc(ixo^s,eaux_)
3303  if(mhd_hall) then
3304  ! f_i[e]= f_i[e] + vHall_i*(b_k*b_k) - b_i*(vHall_k*b_k)
3305  f(ixo^s,e_) = f(ixo^s,e_) + vhall(ixo^s,idim) * &
3306  sum(w(ixo^s, mag(:))**2,dim=ndim+1) &
3307  - w(ixo^s,mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1)
3308  end if
3309  end if
3310  end if
3311 
3312  ! compute flux of magnetic field
3313  ! f_i[b_k]=v_i*b_k-v_k*b_i
3314  do idir=1,ndir
3315  if (idim==idir) then
3316  ! f_i[b_i] should be exactly 0, so we do not use the transport flux
3317  if (mhd_glm) then
3318  f(ixo^s,mag(idir))=w(ixo^s,psi_)
3319  else
3320  f(ixo^s,mag(idir))=zero
3321  end if
3322  else
3323  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))
3324  if (mhd_hall) then
3325  ! f_i[b_k] = f_i[b_k] + vHall_i*b_k - vHall_k*b_i
3326  f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3327  - vhall(ixo^s,idir)*w(ixo^s,mag(idim)) &
3328  + vhall(ixo^s,idim)*w(ixo^s,mag(idir))
3329  end if
3330  end if
3331  end do
3332 
3333  if (mhd_glm) then
3334  !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645
3335  f(ixo^s,psi_) = cmax_global**2*w(ixo^s,mag(idim))
3336  end if
3337 
3338  ! Contributions of ambipolar term in explicit scheme
3339  if(mhd_ambipolar_exp.and. .not.stagger_grid) then
3340  ! ambipolar electric field
3341  ! E_ambi=-eta_ambi*JxBxB=-JaxBxB=B^2*Ja-(Ja dot B)*B
3342  !Ja=eta_ambi*J=J * mhd_eta_ambi/rho**2
3343  allocate(jambi(ixi^s,1:3))
3344  call mhd_get_jambi(w,x,ixi^l,ixo^l,jambi)
3345  allocate(btot(ixo^s,1:3))
3346  btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
3347  allocate(tmp2(ixo^s),tmp3(ixo^s))
3348  !tmp2 = Btot^2
3349  tmp2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1)
3350  !tmp3 = J_ambi dot Btot
3351  tmp3(ixo^s) = sum(jambi(ixo^s,:)*btot(ixo^s,:),dim=ndim+1)
3352 
3353  select case(idim)
3354  case(1)
3355  tmp(ixo^s)=w(ixo^s,mag(3)) *jambi(ixo^s,2) - w(ixo^s,mag(2)) * jambi(ixo^s,3)
3356  f(ixo^s,mag(2))= f(ixo^s,mag(2)) - tmp2(ixo^s) * jambi(ixo^s,3) + tmp3(ixo^s) * btot(ixo^s,3)
3357  f(ixo^s,mag(3))= f(ixo^s,mag(3)) + tmp2(ixo^s) * jambi(ixo^s,2) - tmp3(ixo^s) * btot(ixo^s,2)
3358  case(2)
3359  tmp(ixo^s)=w(ixo^s,mag(1)) *jambi(ixo^s,3) - w(ixo^s,mag(3)) * jambi(ixo^s,1)
3360  f(ixo^s,mag(1))= f(ixo^s,mag(1)) + tmp2(ixo^s) * jambi(ixo^s,3) - tmp3(ixo^s) * btot(ixo^s,3)
3361  f(ixo^s,mag(3))= f(ixo^s,mag(3)) - tmp2(ixo^s) * jambi(ixo^s,1) + tmp3(ixo^s) * btot(ixo^s,1)
3362  case(3)
3363  tmp(ixo^s)=w(ixo^s,mag(2)) *jambi(ixo^s,1) - w(ixo^s,mag(1)) * jambi(ixo^s,2)
3364  f(ixo^s,mag(1))= f(ixo^s,mag(1)) - tmp2(ixo^s) * jambi(ixo^s,2) + tmp3(ixo^s) * btot(ixo^s,2)
3365  f(ixo^s,mag(2))= f(ixo^s,mag(2)) + tmp2(ixo^s) * jambi(ixo^s,1) - tmp3(ixo^s) * btot(ixo^s,1)
3366  endselect
3367 
3368  if(mhd_energy .and. .not. mhd_internal_e) then
3369  f(ixo^s,e_) = f(ixo^s,e_) + tmp2(ixo^s) * tmp(ixo^s)
3370  endif
3371 
3372  deallocate(jambi,btot,tmp2,tmp3)
3373  endif
3374 
3375  end subroutine mhd_get_flux
3376 
3377  !> Calculate fluxes within ixO^L without any splitting
3378  subroutine mhd_get_flux_hde(wC,w,x,ixI^L,ixO^L,idim,f)
3380  use mod_geometry
3381 
3382  integer, intent(in) :: ixI^L, ixO^L, idim
3383  ! conservative w
3384  double precision, intent(in) :: wC(ixI^S,nw)
3385  ! primitive w
3386  double precision, intent(in) :: w(ixI^S,nw)
3387  double precision, intent(in) :: x(ixI^S,1:ndim)
3388  double precision,intent(out) :: f(ixI^S,nwflux)
3389 
3390  double precision :: pgas(ixO^S), ptotal(ixO^S)
3391  double precision :: tmp(ixI^S)
3392  integer :: idirmin, iw, idir, jdir, kdir
3393  double precision, allocatable, dimension(:^D&,:) :: Jambi, btot
3394  double precision, allocatable, dimension(:^D&) :: tmp2, tmp3
3395 
3396  ! Get flux of density
3397  f(ixo^s,rho_)=w(ixo^s,mom(idim))*w(ixo^s,rho_)
3398  ! pgas is time dependent only
3399  if(mhd_energy) then
3400  pgas(ixo^s)=w(ixo^s,p_)
3401  else
3402  pgas(ixo^s)=mhd_adiab*w(ixo^s,rho_)**mhd_gamma
3403  end if
3404 
3405  ptotal(ixo^s)=pgas(ixo^s)+0.5d0*sum(w(ixo^s,mag(:))**2,dim=ndim+1)
3406 
3407  ! Get flux of tracer
3408  do iw=1,mhd_n_tracer
3409  f(ixo^s,tracer(iw))=w(ixo^s,mom(idim))*w(ixo^s,tracer(iw))
3410  end do
3411 
3412  ! Get flux of momentum
3413  ! f_i[m_k]=v_i*m_k-b_k*b_i [+ptotal if i==k]
3414  do idir=1,ndir
3415  if(idim==idir) then
3416  f(ixo^s,mom(idir))=wc(ixo^s,mom(idir))*w(ixo^s,mom(idim))+ptotal(ixo^s)-&
3417  w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3418  else
3419  f(ixo^s,mom(idir))=wc(ixo^s,mom(idir))*w(ixo^s,mom(idim))-&
3420  w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3421  end if
3422  end do
3423 
3424  ! Get flux of energy
3425  if(mhd_energy) then
3426  f(ixo^s,e_)=w(ixo^s,mom(idim))*(wc(ixo^s,e_)+pgas(ixo^s))
3427  end if
3428 
3429  ! compute flux of magnetic field
3430  ! f_i[b_k]=v_i*b_k-v_k*b_i
3431  do idir=1,ndir
3432  if (idim==idir) then
3433  ! f_i[b_i] should be exactly 0, so we do not use the transport flux
3434  if (mhd_glm) then
3435  f(ixo^s,mag(idir))=w(ixo^s,psi_)
3436  else
3437  f(ixo^s,mag(idir))=zero
3438  end if
3439  else
3440  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))
3441  end if
3442  end do
3443 
3444  if (mhd_glm) then
3445  !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645
3446  f(ixo^s,psi_) = cmax_global**2*w(ixo^s,mag(idim))
3447  end if
3448 
3449  ! Contributions of ambipolar term in explicit scheme
3450  if(mhd_ambipolar_exp.and. .not.stagger_grid) then
3451  ! ambipolar electric field
3452  ! E_ambi=-eta_ambi*JxBxB=-JaxBxB=B^2*Ja-(Ja dot B)*B
3453  !Ja=eta_ambi*J=J * mhd_eta_ambi/rho**2
3454  allocate(jambi(ixi^s,1:3))
3455  call mhd_get_jambi(w,x,ixi^l,ixo^l,jambi)
3456  allocate(btot(ixo^s,1:3))
3457  btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
3458  allocate(tmp2(ixo^s),tmp3(ixo^s))
3459  !tmp2 = Btot^2
3460  tmp2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1)
3461  !tmp3 = J_ambi dot Btot
3462  tmp3(ixo^s) = sum(jambi(ixo^s,:)*btot(ixo^s,:),dim=ndim+1)
3463 
3464  select case(idim)
3465  case(1)
3466  tmp(ixo^s)=w(ixo^s,mag(3)) *jambi(ixo^s,2) - w(ixo^s,mag(2)) * jambi(ixo^s,3)
3467  f(ixo^s,mag(2))= f(ixo^s,mag(2)) - tmp2(ixo^s) * jambi(ixo^s,3) + tmp3(ixo^s) * btot(ixo^s,3)
3468  f(ixo^s,mag(3))= f(ixo^s,mag(3)) + tmp2(ixo^s) * jambi(ixo^s,2) - tmp3(ixo^s) * btot(ixo^s,2)
3469  case(2)
3470  tmp(ixo^s)=w(ixo^s,mag(1)) *jambi(ixo^s,3) - w(ixo^s,mag(3)) * jambi(ixo^s,1)
3471  f(ixo^s,mag(1))= f(ixo^s,mag(1)) + tmp2(ixo^s) * jambi(ixo^s,3) - tmp3(ixo^s) * btot(ixo^s,3)
3472  f(ixo^s,mag(3))= f(ixo^s,mag(3)) - tmp2(ixo^s) * jambi(ixo^s,1) + tmp3(ixo^s) * btot(ixo^s,1)
3473  case(3)
3474  tmp(ixo^s)=w(ixo^s,mag(2)) *jambi(ixo^s,1) - w(ixo^s,mag(1)) * jambi(ixo^s,2)
3475  f(ixo^s,mag(1))= f(ixo^s,mag(1)) - tmp2(ixo^s) * jambi(ixo^s,2) + tmp3(ixo^s) * btot(ixo^s,2)
3476  f(ixo^s,mag(2))= f(ixo^s,mag(2)) + tmp2(ixo^s) * jambi(ixo^s,1) - tmp3(ixo^s) * btot(ixo^s,1)
3477  endselect
3478 
3479  if(mhd_energy) then
3480  f(ixo^s,e_) = f(ixo^s,e_) + tmp2(ixo^s) * tmp(ixo^s)
3481  endif
3482 
3483  deallocate(jambi,btot,tmp2,tmp3)
3484  endif
3485 
3486  end subroutine mhd_get_flux_hde
3487 
3488  !> Calculate fluxes within ixO^L with possible splitting
3489  subroutine mhd_get_flux_split(wC,w,x,ixI^L,ixO^L,idim,f)
3491  use mod_geometry
3492 
3493  integer, intent(in) :: ixI^L, ixO^L, idim
3494  ! conservative w
3495  double precision, intent(in) :: wC(ixI^S,nw)
3496  ! primitive w
3497  double precision, intent(in) :: w(ixI^S,nw)
3498  double precision, intent(in) :: x(ixI^S,1:ndim)
3499  double precision,intent(out) :: f(ixI^S,nwflux)
3500 
3501  double precision :: pgas(ixO^S), ptotal(ixO^S), B(ixO^S,1:ndir)
3502  double precision :: tmp(ixI^S)
3503  double precision :: vHall(ixI^S,1:ndir)
3504  integer :: idirmin, iw, idir, jdir, kdir
3505  double precision, allocatable, dimension(:^D&,:) :: Jambi, btot
3506  double precision, allocatable, dimension(:^D&) :: tmp2, tmp3
3507  double precision :: tmp4(ixO^S)
3508 
3509 
3510  call mhd_get_rho(w,x,ixi^l,ixo^l,tmp)
3511  ! Get flux of density
3512  f(ixo^s,rho_)=w(ixo^s,mom(idim))*tmp(ixo^s)
3513  ! pgas is time dependent only
3514  if(mhd_energy) then
3515  pgas(ixo^s)=w(ixo^s,p_)
3516  else
3517  pgas(ixo^s)=mhd_adiab*tmp(ixo^s)**mhd_gamma
3518  if(has_equi_pe0) then
3519  pgas(ixo^s)=pgas(ixo^s)-block%equi_vars(ixo^s,equi_pe0_,b0i)
3520  endif
3521  end if
3522 
3523  if (mhd_hall) then
3524  call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3525  end if
3526 
3527  if(b0field) then
3528  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+block%B0(ixo^s,1:ndir,idim)
3529  pgas(ixo^s)=pgas(ixo^s)+sum(w(ixo^s,mag(:))*block%B0(ixo^s,:,idim),dim=ndim+1)
3530  else
3531  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
3532  end if
3533 
3534  ptotal(ixo^s)=pgas(ixo^s)+0.5d0*sum(w(ixo^s,mag(:))**2,dim=ndim+1)
3535 
3536  ! Get flux of tracer
3537  do iw=1,mhd_n_tracer
3538  f(ixo^s,tracer(iw))=w(ixo^s,mom(idim))*w(ixo^s,tracer(iw))
3539  end do
3540 
3541  ! Get flux of momentum
3542  ! f_i[m_k]=v_i*m_k-b_k*b_i [+ptotal if i==k]
3543  if(b0field) then
3544  do idir=1,ndir
3545  if(idim==idir) then
3546  f(ixo^s,mom(idir))=wc(ixo^s,mom(idir))*w(ixo^s,mom(idim))+ptotal(ixo^s)-&
3547  w(ixo^s,mag(idir))*b(ixo^s,idim)-&
3548  block%B0(ixo^s,idir,idim)*w(ixo^s,mag(idim))
3549  else
3550  f(ixo^s,mom(idir))=wc(ixo^s,mom(idir))*w(ixo^s,mom(idim))-&
3551  w(ixo^s,mag(idir))*b(ixo^s,idim)-&
3552  block%B0(ixo^s,idir,idim)*w(ixo^s,mag(idim))
3553  end if
3554  end do
3555  else
3556  do idir=1,ndir
3557  if(idim==idir) then
3558  f(ixo^s,mom(idir))=wc(ixo^s,mom(idir))*w(ixo^s,mom(idim))+ptotal(ixo^s)-&
3559  w(ixo^s,mag(idir))*b(ixo^s,idim)
3560  else
3561  f(ixo^s,mom(idir))=wc(ixo^s,mom(idir))*w(ixo^s,mom(idim))-&
3562  w(ixo^s,mag(idir))*b(ixo^s,idim)
3563  end if
3564  end do
3565  end if
3566 
3567  ! Get flux of energy
3568  ! f_i[e]=v_i*e+v_i*ptotal-b_i*(b_k*v_k)
3569  if(mhd_energy) then
3570  if (mhd_internal_e) then
3571  f(ixo^s,e_)=w(ixo^s,mom(idim))*wc(ixo^s,e_)
3572  else
3573  f(ixo^s,e_)=w(ixo^s,mom(idim))*(wc(ixo^s,e_)+ptotal(ixo^s))&
3574  -b(ixo^s,idim)*sum(w(ixo^s,mag(:))*w(ixo^s,mom(:)),dim=ndim+1)
3575  if(mhd_solve_eaux) f(ixo^s,eaux_)=w(ixo^s,mom(idim))*wc(ixo^s,eaux_)
3576 
3577  if (mhd_hall) then
3578  ! f_i[e]= f_i[e] + vHall_i*(b_k*b_k) - b_i*(vHall_k*b_k)
3579  if (mhd_etah>zero) then
3580  f(ixo^s,e_) = f(ixo^s,e_) + vhall(ixo^s,idim) * &
3581  sum(w(ixo^s, mag(:))*b(ixo^s,:),dim=ndim+1) &
3582  - b(ixo^s,idim) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1)
3583  end if
3584  end if
3585  end if
3586  if(has_equi_pe0) then
3587  f(ixo^s,e_)= f(ixo^s,e_) &
3588  + w(ixo^s,mom(idim)) * block%equi_vars(ixo^s,equi_pe0_,idim) * inv_gamma_1
3589  end if
3590  end if
3591 
3592  ! compute flux of magnetic field
3593  ! f_i[b_k]=v_i*b_k-v_k*b_i
3594  do idir=1,ndir
3595  if (idim==idir) then
3596  ! f_i[b_i] should be exactly 0, so we do not use the transport flux
3597  if (mhd_glm) then
3598  f(ixo^s,mag(idir))=w(ixo^s,psi_)
3599  else
3600  f(ixo^s,mag(idir))=zero
3601  end if
3602  else
3603  f(ixo^s,mag(idir))=w(ixo^s,mom(idim))*b(ixo^s,idir)-b(ixo^s,idim)*w(ixo^s,mom(idir))
3604 
3605  if (mhd_hall) then
3606  ! f_i[b_k] = f_i[b_k] + vHall_i*b_k - vHall_k*b_i
3607  if (mhd_etah>zero) then
3608  f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3609  - vhall(ixo^s,idir)*b(ixo^s,idim) &
3610  + vhall(ixo^s,idim)*b(ixo^s,idir)
3611  end if
3612  end if
3613 
3614  end if
3615  end do
3616 
3617  if (mhd_glm) then
3618  !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645
3619  f(ixo^s,psi_) = cmax_global**2*w(ixo^s,mag(idim))
3620  end if
3621 
3622  ! Contributions of ambipolar term in explicit scheme
3623  if(mhd_ambipolar_exp.and. .not.stagger_grid) then
3624  ! ambipolar electric field
3625  ! E_ambi=-eta_ambi*JxBxB=-JaxBxB=B^2*Ja-(Ja dot B)*B
3626  !Ja=eta_ambi*J=J * mhd_eta_ambi/rho**2
3627  allocate(jambi(ixi^s,1:3))
3628  call mhd_get_jambi(w,x,ixi^l,ixo^l,jambi)
3629  allocate(btot(ixo^s,1:3))
3630  if(b0field) then
3631  do idir=1,3
3632  btot(ixo^s, idir) = w(ixo^s,mag(idir)) + block%B0(ixo^s,idir,idim)
3633  enddo
3634  else
3635  btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
3636  endif
3637  allocate(tmp2(ixo^s),tmp3(ixo^s))
3638  !tmp2 = Btot^2
3639  tmp2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1)
3640  !tmp3 = J_ambi dot Btot
3641  tmp3(ixo^s) = sum(jambi(ixo^s,:)*btot(ixo^s,:),dim=ndim+1)
3642 
3643  select case(idim)
3644  case(1)
3645  tmp(ixo^s)=w(ixo^s,mag(3)) *jambi(ixo^s,2) - w(ixo^s,mag(2)) * jambi(ixo^s,3)
3646  if(b0field) tmp4(ixo^s) = w(ixo^s,mag(2)) * btot(ixo^s,3) - w(ixo^s,mag(3)) * btot(ixo^s,2)
3647  f(ixo^s,mag(2))= f(ixo^s,mag(2)) - tmp2(ixo^s) * jambi(ixo^s,3) + tmp3(ixo^s) * btot(ixo^s,3)
3648  f(ixo^s,mag(3))= f(ixo^s,mag(3)) + tmp2(ixo^s) * jambi(ixo^s,2) - tmp3(ixo^s) * btot(ixo^s,2)
3649  case(2)
3650  tmp(ixo^s)=w(ixo^s,mag(1)) *jambi(ixo^s,3) - w(ixo^s,mag(3)) * jambi(ixo^s,1)
3651  if(b0field) tmp4(ixo^s) = w(ixo^s,mag(3)) * btot(ixo^s,1) - w(ixo^s,mag(1)) * btot(ixo^s,3)
3652  f(ixo^s,mag(1))= f(ixo^s,mag(1)) + tmp2(ixo^s) * jambi(ixo^s,3) - tmp3(ixo^s) * btot(ixo^s,3)
3653  f(ixo^s,mag(3))= f(ixo^s,mag(3)) - tmp2(ixo^s) * jambi(ixo^s,1) + tmp3(ixo^s) * btot(ixo^s,1)
3654  case(3)
3655  tmp(ixo^s)=w(ixo^s,mag(2)) *jambi(ixo^s,1) - w(ixo^s,mag(1)) * jambi(ixo^s,2)
3656  if(b0field) tmp4(ixo^s) = w(ixo^s,mag(1)) * btot(ixo^s,2) - w(ixo^s,mag(2)) * btot(ixo^s,1)
3657  f(ixo^s,mag(1))= f(ixo^s,mag(1)) - tmp2(ixo^s) * jambi(ixo^s,2) + tmp3(ixo^s) * btot(ixo^s,2)
3658  f(ixo^s,mag(2))= f(ixo^s,mag(2)) + tmp2(ixo^s) * jambi(ixo^s,1) - tmp3(ixo^s) * btot(ixo^s,1)
3659  endselect
3660 
3661  if(mhd_energy .and. .not. mhd_internal_e) then
3662  f(ixo^s,e_) = f(ixo^s,e_) + tmp2(ixo^s) * tmp(ixo^s)
3663  if(b0field) f(ixo^s,e_) = f(ixo^s,e_) + tmp3(ixo^s) * tmp4(ixo^s)
3664  endif
3665 
3666  deallocate(jambi,btot,tmp2,tmp3)
3667  endif
3668 
3669  end subroutine mhd_get_flux_split
3670 
3671  !> Calculate semirelativistic fluxes within ixO^L without any splitting
3672  subroutine mhd_get_flux_semirelati(wC,w,x,ixI^L,ixO^L,idim,f)
3674  use mod_geometry
3675 
3676  integer, intent(in) :: ixI^L, ixO^L, idim
3677  ! conservative w
3678  double precision, intent(in) :: wC(ixI^S,nw)
3679  ! primitive w
3680  double precision, intent(in) :: w(ixI^S,nw)
3681  double precision, intent(in) :: x(ixI^S,1:ndim)
3682  double precision,intent(out) :: f(ixI^S,nwflux)
3683 
3684  double precision :: pgas(ixO^S)
3685  double precision :: SA(ixO^S), E(ixO^S,1:ndir), B(ixO^S,1:ndir)
3686  integer :: idirmin, iw, idir, jdir, kdir
3687 
3688  ! gas thermal pressure
3689  if(mhd_energy) then
3690  pgas(ixo^s)=w(ixo^s,p_)
3691  else
3692  pgas(ixo^s)=mhd_adiab*w(ixo^s,rho_)**mhd_gamma
3693  end if
3694 
3695  ! Get flux of density
3696  f(ixo^s,rho_)=w(ixo^s,mom(idim))*w(ixo^s,rho_)
3697 
3698  ! Get flux of tracer
3699  do iw=1,mhd_n_tracer
3700  f(ixo^s,tracer(iw))=w(ixo^s,mom(idim))*w(ixo^s,tracer(iw))
3701  end do
3702  ! E=-uxB=Bxu
3703  if(b0field) then
3704  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+block%B0(ixo^s,1:ndir,idim)
3705  pgas(ixo^s)=pgas(ixo^s)+sum(w(ixo^s,mag(:))*block%B0(ixo^s,:,idim),dim=ndim+1)
3706  else
3707  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
3708  end if
3709  e=0.d0
3710  do idir=1,ndir; do jdir=1,ndir; do kdir=1,ndir
3711  if(lvc(idir,jdir,kdir)==1)then
3712  e(ixo^s,idir)=e(ixo^s,idir)+b(ixo^s,jdir)*w(ixo^s,mom(kdir))
3713  else if(lvc(idir,jdir,kdir)==-1)then
3714  e(ixo^s,idir)=e(ixo^s,idir)-b(ixo^s,jdir)*w(ixo^s,mom(kdir))
3715  end if
3716  end do; end do; end do
3717 
3718  pgas(ixo^s)=pgas(ixo^s)+half*(sum(w(ixo^s,mag(:))**2,dim=ndim+1)+&
3719  sum(e(ixo^s,:)**2,dim=ndim+1)*inv_squared_c)
3720 
3721  ! Get flux of momentum
3722  if(b0field) then
3723  do idir=1,ndir
3724  if(idim==idir) then
3725  f(ixo^s,mom(idir))=w(ixo^s,rho_)*w(ixo^s,mom(idir))*w(ixo^s,mom(idim))+pgas&
3726  -w(ixo^s,mag(idir))*b(ixo^s,idim)-e(ixo^s,idir)*e(ixo^s,idim)*inv_squared_c&
3727  -block%B0(ixo^s,idir,idim)*w(ixo^s,mag(idim))
3728  else
3729  f(ixo^s,mom(idir))=w(ixo^s,rho_)*w(ixo^s,mom(idir))*w(ixo^s,mom(idim))&
3730  -w(ixo^s,mag(idir))*b(ixo^s,idim)-e(ixo^s,idir)*e(ixo^s,idim)*inv_squared_c&
3731  -block%B0(ixo^s,idir,idim)*w(ixo^s,mag(idim))
3732  end if
3733  end do
3734  else
3735  do idir=1,ndir
3736  if(idim==idir) then
3737  f(ixo^s,mom(idir))=w(ixo^s,rho_)*w(ixo^s,mom(idir))*w(ixo^s,mom(idim))+pgas&
3738  -w(ixo^s,mag(idir))*b(ixo^s,idim)-e(ixo^s,idir)*e(ixo^s,idim)*inv_squared_c
3739  else
3740  f(ixo^s,mom(idir))=w(ixo^s,rho_)*w(ixo^s,mom(idir))*w(ixo^s,mom(idim))&
3741  -w(ixo^s,mag(idir))*b(ixo^s,idim)-e(ixo^s,idir)*e(ixo^s,idim)*inv_squared_c
3742  end if
3743  end do
3744  end if
3745 
3746  ! Get flux of energy
3747  if(mhd_energy) then
3748  sa=0.d0
3749  do jdir=1,ndir; do kdir=1,ndir
3750  if(lvc(idim,jdir,kdir)==1)then
3751  sa(ixo^s)=sa(ixo^s)+e(ixo^s,jdir)*w(ixo^s,mag(kdir))
3752  else if(lvc(idim,jdir,kdir)==-1) then
3753  sa(ixo^s)=sa(ixo^s)-e(ixo^s,jdir)*w(ixo^s,mag(kdir))
3754  end if
3755  end do; end do
3756  f(ixo^s,e_)=w(ixo^s,mom(idim))*(half*w(ixo^s,rho_)*sum(w(ixo^s,mom(:))**2,dim=ndim+1)+&
3757  mhd_gamma*pgas*inv_gamma_1)+sa(ixo^s)
3758  end if
3759 
3760  ! compute flux of magnetic field
3761  ! f_i[b_k]=v_i*b_k-v_k*b_i
3762  do idir=1,ndir
3763  if (idim==idir) then
3764  ! f_i[b_i] should be exactly 0, so we do not use the transport flux
3765  if (mhd_glm) then
3766  f(ixo^s,mag(idir))=w(ixo^s,psi_)
3767  else
3768  f(ixo^s,mag(idir))=zero
3769  end if
3770  else
3771  f(ixo^s,mag(idir))=w(ixo^s,mom(idim))*b(ixo^s,idir)-b(ixo^s,idim)*w(ixo^s,mom(idir))
3772  end if
3773  end do
3774 
3775  if (mhd_glm) then
3776  !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645
3777  f(ixo^s,psi_) = cmax_global**2*w(ixo^s,mag(idim))
3778  end if
3779 
3780  end subroutine mhd_get_flux_semirelati
3781 
3782  !> Source terms J.E in internal energy.
3783  !> For the ambipolar term E = ambiCoef * JxBxB=ambiCoef * B^2(-J_perpB)
3784  !=> the source term J.E = ambiCoef * B^2 * J_perpB^2 = ambiCoef * JxBxB^2/B^2
3785  !> ambiCoef is calculated as mhd_ambi_coef/rho^2, see also the subroutine mhd_get_Jambi
3786  subroutine add_source_ambipolar_internal_energy(qdt,ixI^L,ixO^L,wCT,w,x,ie)
3788  integer, intent(in) :: ixI^L, ixO^L,ie
3789  double precision, intent(in) :: qdt
3790  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3791  double precision, intent(inout) :: w(ixI^S,1:nw)
3792  double precision :: tmp(ixI^S)
3793  double precision :: jxbxb(ixI^S,1:3)
3794 
3795  call mhd_get_jxbxb(wct,x,ixi^l,ixo^l,jxbxb)
3796  tmp(ixo^s) = sum(jxbxb(ixo^s,1:3)**2,dim=ndim+1) / mhd_mag_en_all(wct, ixi^l, ixo^l)
3797  call multiplyambicoef(ixi^l,ixo^l,tmp,wct,x)
3798  w(ixo^s,ie)=w(ixo^s,ie)+qdt * tmp
3799 
3801 
3802  subroutine mhd_get_jxbxb(w,x,ixI^L,ixO^L,res)
3804 
3805  integer, intent(in) :: ixI^L, ixO^L
3806  double precision, intent(in) :: w(ixI^S,nw)
3807  double precision, intent(in) :: x(ixI^S,1:ndim)
3808  double precision, intent(out) :: res(:^D&,:)
3809 
3810  double precision :: btot(ixI^S,1:3)
3811  integer :: idir, idirmin
3812  double precision :: current(ixI^S,7-2*ndir:3)
3813  double precision :: tmp(ixI^S),b2(ixI^S)
3814 
3815  res=0.d0
3816  ! Calculate current density and idirmin
3817  call get_current(w,ixi^l,ixo^l,idirmin,current)
3818  !!!here we know that current has nonzero values only for components in the range idirmin, 3
3819 
3820  if(b0field) then
3821  do idir=1,3
3822  btot(ixo^s, idir) = w(ixo^s,mag(idir)) + block%B0(ixo^s,idir,b0i)
3823  enddo
3824  else
3825  btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
3826  endif
3827 
3828  tmp(ixo^s) = sum(current(ixo^s,idirmin:3)*btot(ixo^s,idirmin:3),dim=ndim+1) !J.B
3829  b2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1) !B^2
3830  do idir=1,idirmin-1
3831  res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s)
3832  enddo
3833  do idir=idirmin,3
3834  res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s) - current(ixo^s,idir) * b2(ixo^s)
3835  enddo
3836  end subroutine mhd_get_jxbxb
3837 
3838  !> Sets the sources for the ambipolar
3839  !> this is used for the STS method
3840  ! The sources are added directly (instead of fluxes as in the explicit)
3841  !> at the corresponding indices
3842  !> store_flux_var is explicitly called for each of the fluxes one by one
3843  subroutine sts_set_source_ambipolar(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
3845  use mod_fix_conserve
3846 
3847  integer, intent(in) :: ixI^L, ixO^L,igrid,nflux
3848  double precision, intent(in) :: x(ixI^S,1:ndim)
3849  double precision, intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
3850  double precision, intent(in) :: my_dt
3851  logical, intent(in) :: fix_conserve_at_step
3852 
3853  double precision, dimension(ixI^S,1:3) :: tmp,ff
3854  double precision :: fluxall(ixI^S,1:nflux,1:ndim)
3855  double precision :: fE(ixI^S,7-2*ndim:3)
3856  double precision :: btot(ixI^S,1:3),tmp2(ixI^S)
3857  integer :: i, ixA^L, ie_
3858 
3859  ixa^l=ixo^l^ladd1;
3860 
3861  fluxall=zero
3862 
3863  call mhd_get_jxbxb(w,x,ixi^l,ixa^l,tmp)
3864 
3865  !set electric field in tmp: E=nuA * jxbxb, where nuA=-etaA/rho^2
3866  do i=1,3
3867  !tmp(ixA^S,i) = -(mhd_eta_ambi/w(ixA^S, rho_)**2) * tmp(ixA^S,i)
3868  call multiplyambicoef(ixi^l,ixa^l,tmp(ixi^s,i),w,x)
3869  enddo
3870 
3871  if(mhd_energy .and. .not.mhd_internal_e) then
3872  !btot should be only mag. pert.
3873  btot(ixa^s,1:3)=0.d0
3874  !if(B0field) then
3875  ! do i=1,ndir
3876  ! btot(ixA^S, i) = w(ixA^S,mag(i)) + block%B0(ixA^S,i,0)
3877  ! enddo
3878  !else
3879  btot(ixa^s,1:ndir) = w(ixa^s,mag(1:ndir))
3880  !endif
3881  call cross_product(ixi^l,ixa^l,tmp,btot,ff)
3882  call get_flux_on_cell_face(ixi^l,ixo^l,ff,tmp2)
3883  if(fix_conserve_at_step) fluxall(ixi^s,1,1:ndim)=ff(ixi^s,1:ndim)
3884  !- sign comes from the fact that the flux divergence is a source now
3885  wres(ixo^s,e_)=-tmp2(ixo^s)
3886  endif
3887 
3888  if(stagger_grid) then
3889  if(ndir>ndim) then
3890  !!!Bz
3891  ff(ixa^s,1) = tmp(ixa^s,2)
3892  ff(ixa^s,2) = -tmp(ixa^s,1)
3893  ff(ixa^s,3) = 0.d0
3894  call get_flux_on_cell_face(ixi^l,ixo^l,ff,tmp2)
3895  if(fix_conserve_at_step) fluxall(ixi^s,1+ndir,1:ndim)=ff(ixi^s,1:ndim)
3896  wres(ixo^s,mag(ndir))=-tmp2(ixo^s)
3897  end if
3898  fe=0.d0
3899  call update_faces_ambipolar(ixi^l,ixo^l,w,x,tmp,fe,btot)
3900  ixamax^d=ixomax^d;
3901  ixamin^d=ixomin^d-1;
3902  wres(ixa^s,mag(1:ndim))=-btot(ixa^s,1:ndim)
3903  else
3904  !write curl(ele) as the divergence
3905  !m1={0,ele[[3]],-ele[[2]]}
3906  !m2={-ele[[3]],0,ele[[1]]}
3907  !m3={ele[[2]],-ele[[1]],0}
3908 
3909  !!!Bx
3910  ff(ixa^s,1) = 0.d0
3911  ff(ixa^s,2) = tmp(ixa^s,3)
3912  ff(ixa^s,3) = -tmp(ixa^s,2)
3913  call get_flux_on_cell_face(ixi^l,ixo^l,ff,tmp2)
3914  if(fix_conserve_at_step) fluxall(ixi^s,2,1:ndim)=ff(ixi^s,1:ndim)
3915  !flux divergence is a source now
3916  wres(ixo^s,mag(1))=-tmp2(ixo^s)
3917  !!!By
3918  ff(ixa^s,1) = -tmp(ixa^s,3)
3919  ff(ixa^s,2) = 0.d0
3920  ff(ixa^s,3) = tmp(ixa^s,1)
3921  call get_flux_on_cell_face(ixi^l,ixo^l,ff,tmp2)
3922  if(fix_conserve_at_step) fluxall(ixi^s,3,1:ndim)=ff(ixi^s,1:ndim)
3923  wres(ixo^s,mag(2))=-tmp2(ixo^s)
3924 
3925  if(ndir==3) then
3926  !!!Bz
3927  ff(ixa^s,1) = tmp(ixa^s,2)
3928  ff(ixa^s,2) = -tmp(ixa^s,1)
3929  ff(ixa^s,3) = 0.d0
3930  call get_flux_on_cell_face(ixi^l,ixo^l,ff,tmp2)
3931  if(fix_conserve_at_step) fluxall(ixi^s,1+ndir,1:ndim)=ff(ixi^s,1:ndim)
3932  wres(ixo^s,mag(ndir))=-tmp2(ixo^s)
3933  end if
3934 
3935  end if
3936 
3937  if(fix_conserve_at_step) then
3938  fluxall=my_dt*fluxall
3939  call store_flux(igrid,fluxall,1,ndim,nflux)
3940  if(stagger_grid) then
3941  call store_edge(igrid,ixi^l,my_dt*fe,1,ndim)
3942  end if
3943  end if
3944 
3945  end subroutine sts_set_source_ambipolar
3946 
3947  !> get ambipolar electric field and the integrals around cell faces
3948  subroutine update_faces_ambipolar(ixI^L,ixO^L,w,x,ECC,fE,circ)
3950 
3951  integer, intent(in) :: ixI^L, ixO^L
3952  double precision, intent(in) :: w(ixI^S,1:nw)
3953  double precision, intent(in) :: x(ixI^S,1:ndim)
3954  ! amibipolar electric field at cell centers
3955  double precision, intent(in) :: ECC(ixI^S,1:3)
3956  double precision, intent(out) :: fE(ixI^S,7-2*ndim:3)
3957  double precision, intent(out) :: circ(ixI^S,1:ndim)
3958 
3959  integer :: hxC^L,ixC^L,ixA^L
3960  integer :: idim1,idim2,idir,ix^D
3961 
3962  fe=zero
3963  ! calcuate ambipolar electric field on cell edges from cell centers
3964  do idir=7-2*ndim,3
3965  ixcmax^d=ixomax^d;
3966  ixcmin^d=ixomin^d+kr(idir,^d)-1;
3967  {do ix^db=0,1\}
3968  if({ ix^d==1 .and. ^d==idir | .or.}) cycle
3969  ixamin^d=ixcmin^d+ix^d;
3970  ixamax^d=ixcmax^d+ix^d;
3971  fe(ixc^s,idir)=fe(ixc^s,idir)+ecc(ixa^s,idir)
3972  {end do\}
3973  fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0*block%dsC(ixc^s,idir)
3974  end do
3975 
3976  ! Calculate circulation on each face to get value of line integral of
3977  ! electric field in the positive idir direction.
3978  ixcmax^d=ixomax^d;
3979  ixcmin^d=ixomin^d-1;
3980 
3981  circ=zero
3982 
3983  do idim1=1,ndim ! Coordinate perpendicular to face
3984  do idim2=1,ndim
3985  do idir=7-2*ndim,3 ! Direction of line integral
3986  ! Assemble indices
3987  hxc^l=ixc^l-kr(idim2,^d);
3988  ! Add line integrals in direction idir
3989  circ(ixc^s,idim1)=circ(ixc^s,idim1)&
3990  +lvc(idim1,idim2,idir)&
3991  *(fe(ixc^s,idir)&
3992  -fe(hxc^s,idir))
3993  end do
3994  end do
3995  circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
3996  end do
3997 
3998  end subroutine update_faces_ambipolar
3999 
4000  !> use cell-center flux to get cell-face flux
4001  !> and get the source term as the divergence of the flux
4002  subroutine get_flux_on_cell_face(ixI^L,ixO^L,ff,src)
4004 
4005  integer, intent(in) :: ixI^L, ixO^L
4006  double precision, dimension(:^D&,:), intent(inout) :: ff
4007  double precision, intent(out) :: src(ixI^S)
4008 
4009  double precision :: ffc(ixI^S,1:ndim)
4010  double precision :: dxinv(ndim)
4011  integer :: idims, ix^D, ixA^L, ixB^L, ixC^L
4012 
4013  ixa^l=ixo^l^ladd1;
4014  dxinv=1.d0/dxlevel
4015  ! cell corner flux in ffc
4016  ffc=0.d0
4017  ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
4018  {do ix^db=0,1\}
4019  ixbmin^d=ixcmin^d+ix^d;
4020  ixbmax^d=ixcmax^d+ix^d;
4021  ffc(ixc^s,1:ndim)=ffc(ixc^s,1:ndim)+ff(ixb^s,1:ndim)
4022  {end do\}
4023  ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
4024  ! flux at cell face
4025  ff(ixi^s,1:ndim)=0.d0
4026  do idims=1,ndim
4027  ixb^l=ixo^l-kr(idims,^d);
4028  ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
4029  {do ix^db=0,1 \}
4030  if({ ix^d==0 .and. ^d==idims | .or.}) then
4031  ixbmin^d=ixcmin^d-ix^d;
4032  ixbmax^d=ixcmax^d-ix^d;
4033  ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
4034  end if
4035  {end do\}
4036  ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
4037  end do
4038  src=0.d0
4039  if(slab_uniform) then
4040  do idims=1,ndim
4041  ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
4042  ixb^l=ixo^l-kr(idims,^d);
4043  src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4044  end do
4045  else
4046  do idims=1,ndim
4047  ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
4048  ixb^l=ixo^l-kr(idims,^d);
4049  src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4050  end do
4051  src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
4052  end if
4053  end subroutine get_flux_on_cell_face
4054 
4055  !> Calculates the explicit dt for the ambipokar term
4056  !> This function is used by both explicit scheme and STS method
4057  function get_ambipolar_dt(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
4059 
4060  integer, intent(in) :: ixi^l, ixo^l
4061  double precision, intent(in) :: dx^d, x(ixi^s,1:ndim)
4062  double precision, intent(in) :: w(ixi^s,1:nw)
4063  double precision :: dtnew
4064 
4065  double precision :: coef
4066  double precision :: dxarr(ndim)
4067  double precision :: tmp(ixi^s)
4068 
4069  ^d&dxarr(^d)=dx^d;
4070  tmp(ixo^s) = mhd_mag_en_all(w, ixi^l, ixo^l)
4071  call multiplyambicoef(ixi^l,ixo^l,tmp,w,x)
4072  coef = maxval(abs(tmp(ixo^s)))
4073  if(coef/=0.d0) then
4074  coef=1.d0/coef
4075  else
4076  coef=bigdouble
4077  end if
4078  if(slab_uniform) then
4079  dtnew=minval(dxarr(1:ndim))**2.0d0*coef
4080  else
4081  dtnew=minval(block%ds(ixo^s,1:ndim))**2.0d0*coef
4082  end if
4083 
4084  end function get_ambipolar_dt
4085 
4086  !> multiply res by the ambipolar coefficient
4087  !> The ambipolar coefficient is calculated as -mhd_eta_ambi/rho^2
4088  !> The user may mask its value in the user file
4089  !> by implemneting usr_mask_ambipolar subroutine
4090  subroutine multiplyambicoef(ixI^L,ixO^L,res,w,x)
4092  integer, intent(in) :: ixi^l, ixo^l
4093  double precision, intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:ndim)
4094  double precision, intent(inout) :: res(ixi^s)
4095  double precision :: tmp(ixi^s)
4096  double precision :: rho(ixi^s)
4097 
4098  call mhd_get_rho(w,x,ixi^l,ixo^l,rho)
4099  tmp=0.d0
4100  tmp(ixo^s)=-mhd_eta_ambi/rho(ixo^s)**2
4101  if (associated(usr_mask_ambipolar)) then
4102  call usr_mask_ambipolar(ixi^l,ixo^l,w,x,tmp)
4103  end if
4104 
4105  res(ixo^s) = tmp(ixo^s) * res(ixo^s)
4106  end subroutine multiplyambicoef
4107 
4108  !> w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
4109  subroutine mhd_add_source(qdt,ixI^L,ixO^L,wCT,w,x,qsourcesplit,active,wCTprim)
4113  use mod_gravity, only: gravity_add_source
4114  use mod_cak_force, only: cak_add_source
4115 
4116  integer, intent(in) :: ixI^L, ixO^L
4117  double precision, intent(in) :: qdt
4118  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4119  double precision, intent(inout) :: w(ixI^S,1:nw)
4120  logical, intent(in) :: qsourcesplit
4121  logical, intent(inout) :: active
4122  double precision, intent(in), optional :: wCTprim(ixI^S,1:nw)
4123 
4124  if (.not. qsourcesplit) then
4125  if(mhd_internal_e) then
4126  ! Source for solving internal energy
4127  active = .true.
4128  call internal_energy_add_source(qdt,ixi^l,ixo^l,wct,w,x,e_)
4129  else
4130  if(mhd_solve_eaux) then
4131  ! Source for auxiliary internal energy equation
4132  active = .true.
4133  call internal_energy_add_source(qdt,ixi^l,ixo^l,wct,w,x,eaux_)
4134  endif
4135  if(has_equi_pe0) then
4136  active = .true.
4137  call add_pe0_divv(qdt,ixi^l,ixo^l,wct,w,x)
4138  endif
4139  endif
4140 
4141  ! Source for B0 splitting
4142  if (b0field.or.b0field) then
4143  active = .true.
4144  call add_source_b0split(qdt,ixi^l,ixo^l,wct,w,x)
4145  end if
4146 
4147  ! Sources for resistivity in eqs. for e, B1, B2 and B3
4148  if (abs(mhd_eta)>smalldouble)then
4149  active = .true.
4150  call add_source_res2(qdt,ixi^l,ixo^l,wct,w,x)
4151  end if
4152 
4153  if (mhd_eta_hyper>0.d0)then
4154  active = .true.
4155  call add_source_hyperres(qdt,ixi^l,ixo^l,wct,w,x)
4156  end if
4157  ! add sources for semirelativistic MHD
4158  if (unsplit_semirelativistic) then
4159  active = .true.
4160  call add_source_semirelativistic(qdt,ixi^l,ixo^l,wct,w,x,wctprim)
4161  end if
4162  ! add sources for hydrodynamic energy version of MHD
4163  if (mhd_hydrodynamic_e) then
4164  active = .true.
4165  call add_source_hydrodynamic_e(qdt,ixi^l,ixo^l,wct,w,x,wctprim)
4166  end if
4167  end if
4168 
4169  {^nooned
4170  if(.not.source_split_divb .and. .not.qsourcesplit .and. istep==nstep) then
4171  ! Sources related to div B
4172  select case (type_divb)
4173  case (divb_none)
4174  ! Do nothing
4175  case (divb_glm)
4176  active = .true.
4177  call add_source_glm(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
4178  case (divb_powel)
4179  active = .true.
4180  call add_source_powel(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
4181  case (divb_janhunen)
4182  active = .true.
4183  call add_source_janhunen(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
4184  case (divb_linde)
4185  active = .true.
4186  call add_source_linde(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
4187  case (divb_lindejanhunen)
4188  active = .true.
4189  call add_source_linde(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
4190  call add_source_janhunen(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
4191  case (divb_lindepowel)
4192  active = .true.
4193  call add_source_linde(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
4194  call add_source_powel(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
4195  case (divb_lindeglm)
4196  active = .true.
4197  call add_source_linde(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
4198  call add_source_glm(dt,ixi^l,ixo^l,pso(block%igrid)%w,w,x)
4199  case (divb_ct)
4200  continue ! Do nothing
4201  case (divb_multigrid)
4202  continue ! Do nothing
4203  case default
4204  call mpistop('Unknown divB fix')
4205  end select
4206  else if(source_split_divb .and. qsourcesplit) then
4207  ! Sources related to div B
4208  select case (type_divb)
4209  case (divb_none)
4210  ! Do nothing
4211  case (divb_glm)
4212  active = .true.
4213  call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
4214  case (divb_powel)
4215  active = .true.
4216  call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
4217  case (divb_janhunen)
4218  active = .true.
4219  call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
4220  case (divb_linde)
4221  active = .true.
4222  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
4223  case (divb_lindejanhunen)
4224  active = .true.
4225  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
4226  call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
4227  case (divb_lindepowel)
4228  active = .true.
4229  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
4230  call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
4231  case (divb_lindeglm)
4232  active = .true.
4233  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
4234  call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
4235  case (divb_ct)
4236  continue ! Do nothing
4237  case (divb_multigrid)
4238  continue ! Do nothing
4239  case default
4240  call mpistop('Unknown divB fix')
4241  end select
4242  end if
4243  }
4244 
4245  if(mhd_radiative_cooling) then
4246  call radiative_cooling_add_source(qdt,ixi^l,ixo^l,wct,&
4247  w,x,qsourcesplit,active, rc_fl)
4248  end if
4249 
4250  if(mhd_viscosity) then
4251  call viscosity_add_source(qdt,ixi^l,ixo^l,wct,&
4252  w,x,mhd_energy,qsourcesplit,active)
4253  end if
4254 
4255  if(mhd_gravity) then
4256  call gravity_add_source(qdt,ixi^l,ixo^l,wct,&
4257  w,x,gravity_energy,qsourcesplit,active)
4258  end if
4259 
4260  if (mhd_cak_force) then
4261  call cak_add_source(qdt,ixi^l,ixo^l,wct,w,x,mhd_energy,qsourcesplit,active)
4262  end if
4263 
4264  end subroutine mhd_add_source
4265 
4266  subroutine add_pe0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
4268  use mod_geometry
4269 
4270  integer, intent(in) :: ixI^L, ixO^L
4271  double precision, intent(in) :: qdt
4272  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4273  double precision, intent(inout) :: w(ixI^S,1:nw)
4274  double precision :: v(ixI^S,1:ndir)
4275  double precision :: divv(ixI^S)
4276 
4277  call mhd_get_v(wct,x,ixi^l,ixi^l,v)
4278 
4279  if(slab_uniform) then
4280  if(nghostcells .gt. 2) then
4281  call divvector(v,ixi^l,ixo^l,divv,sixthorder=.true.)
4282  else
4283  call divvector(v,ixi^l,ixo^l,divv,fourthorder=.true.)
4284  end if
4285  else
4286  call divvector(v,ixi^l,ixo^l,divv)
4287  end if
4288  w(ixo^s,e_)=w(ixo^s,e_)-qdt*block%equi_vars(ixo^s,equi_pe0_,0)*divv(ixo^s)
4289 
4290  end subroutine add_pe0_divv
4291 
4292  !> Compute the Lorentz force (JxB)
4293  subroutine get_lorentz_force(ixI^L,ixO^L,w,JxB)
4295  integer, intent(in) :: ixI^L, ixO^L
4296  double precision, intent(in) :: w(ixI^S,1:nw)
4297  double precision, intent(inout) :: JxB(ixI^S,3)
4298  double precision :: a(ixI^S,3), b(ixI^S,3)
4299  integer :: idir, idirmin
4300  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
4301  double precision :: current(ixI^S,7-2*ndir:3)
4302 
4303  b=0.0d0
4304  do idir = 1, ndir
4305  b(ixo^s, idir) = mhd_mag_i_all(w, ixi^l, ixo^l,idir)
4306  end do
4307 
4308  ! store J current in a
4309  call get_current(w,ixi^l,ixo^l,idirmin,current)
4310 
4311  a=0.0d0
4312  do idir=7-2*ndir,3
4313  a(ixo^s,idir)=current(ixo^s,idir)
4314  end do
4315 
4316  call cross_product(ixi^l,ixo^l,a,b,jxb)
4317  end subroutine get_lorentz_force
4318 
4319  !> Compute 1/(1+v_A^2/c^2) for semirelativistic MHD, where v_A is the Alfven
4320  !> velocity
4321  subroutine mhd_gamma2_alfven(ixI^L, ixO^L, w, gamma_A2)
4323  integer, intent(in) :: ixI^L, ixO^L
4324  double precision, intent(in) :: w(ixI^S, nw)
4325  double precision, intent(out) :: gamma_A2(ixO^S)
4326  double precision :: rho(ixI^S)
4327 
4328  ! mhd_get_rho cannot be used as x is not a param
4329  if(has_equi_rho0) then
4330  rho(ixo^s) = w(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,b0i)
4331  else
4332  rho(ixo^s) = w(ixo^s,rho_)
4333  endif
4334  ! Compute the inverse of 1 + B^2/(rho * c^2)
4335  gamma_a2(ixo^s) = 1.0d0/(1.0d0+mhd_mag_en_all(w, ixi^l, ixo^l)/rho(ixo^s)*inv_squared_c)
4336  end subroutine mhd_gamma2_alfven
4337 
4338  !> Compute 1/sqrt(1+v_A^2/c^2) for semirelativisitic MHD, where v_A is the
4339  !> Alfven velocity
4340  function mhd_gamma_alfven(w, ixI^L, ixO^L) result(gamma_A)
4342  integer, intent(in) :: ixi^l, ixo^l
4343  double precision, intent(in) :: w(ixi^s, nw)
4344  double precision :: gamma_a(ixo^s)
4345 
4346  call mhd_gamma2_alfven(ixi^l, ixo^l, w, gamma_a)
4347  gamma_a = sqrt(gamma_a)
4348  end function mhd_gamma_alfven
4349 
4350  subroutine internal_energy_add_source(qdt,ixI^L,ixO^L,wCT,w,x,ie)
4352  use mod_geometry
4353 
4354  integer, intent(in) :: ixI^L, ixO^L,ie
4355  double precision, intent(in) :: qdt
4356  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4357  double precision, intent(inout) :: w(ixI^S,1:nw)
4358 
4359  double precision :: v(ixI^S,1:ndir),divv(ixI^S)
4360 
4361  call mhd_get_v(wct,x,ixi^l,ixi^l,v)
4362  if(slab_uniform) then
4363  if(nghostcells .gt. 2) then
4364  call divvector(v,ixi^l,ixo^l,divv,sixthorder=.true.)
4365  else
4366  call divvector(v,ixi^l,ixo^l,divv,fourthorder=.true.)
4367  end if
4368  else
4369  call divvector(v,ixi^l,ixo^l,divv)
4370  end if
4371  w(ixo^s,ie)=w(ixo^s,ie)-qdt*wct(ixo^s,ie)*gamma_1*divv(ixo^s)
4372  if(mhd_ambipolar)then
4373  call add_source_ambipolar_internal_energy(qdt,ixi^l,ixo^l,wct,w,x,ie)
4374  end if
4375  if(fix_small_values) then
4376  call mhd_handle_small_ei(w,x,ixi^l,ixo^l,ie,'internal_energy_add_source')
4377  end if
4378  end subroutine internal_energy_add_source
4379 
4380  subroutine mhd_get_rho(w,x,ixI^L,ixO^L,rho)
4382  integer, intent(in) :: ixi^l, ixo^l
4383  double precision, intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:ndim)
4384  double precision, intent(out) :: rho(ixi^s)
4385 
4386  if(has_equi_rho0) then
4387  rho(ixo^s) = w(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,b0i)
4388  else
4389  rho(ixo^s) = w(ixo^s,rho_)
4390  endif
4391 
4392  end subroutine mhd_get_rho
4393 
4394  !> handle small or negative internal energy
4395  subroutine mhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
4397  use mod_small_values
4398  integer, intent(in) :: ixI^L,ixO^L, ie
4399  double precision, intent(inout) :: w(ixI^S,1:nw)
4400  double precision, intent(in) :: x(ixI^S,1:ndim)
4401  character(len=*), intent(in) :: subname
4402 
4403  integer :: idir
4404  logical :: flag(ixI^S,1:nw)
4405  double precision :: rho(ixI^S)
4406 
4407  flag=.false.
4408  if(has_equi_pe0) then
4409  where(w(ixo^s,ie)+block%equi_vars(ixo^s,equi_pe0_,0)*inv_gamma_1<small_e)&
4410  flag(ixo^s,ie)=.true.
4411  else
4412  where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
4413  endif
4414  if(any(flag(ixo^s,ie))) then
4415  select case (small_values_method)
4416  case ("replace")
4417  if(has_equi_pe0) then
4418  where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
4419  block%equi_vars(ixo^s,equi_pe0_,0)*inv_gamma_1
4420  else
4421  where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
4422  endif
4423  case ("average")
4424  call small_values_average(ixi^l, ixo^l, w, x, flag, ie)
4425  case default
4426  ! small values error shows primitive variables
4427  w(ixo^s,e_)=w(ixo^s,e_)*gamma_1
4428  call mhd_get_rho(w,x,ixi^l,ixo^l,rho)
4429  do idir = 1, ndir
4430  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))/rho(ixo^s)
4431  end do
4432  call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
4433  end select
4434  end if
4435 
4436  end subroutine mhd_handle_small_ei
4437 
4438  !> Source terms after split off time-independent magnetic field
4439  subroutine add_source_b0split(qdt,ixI^L,ixO^L,wCT,w,x)
4441 
4442  integer, intent(in) :: ixI^L, ixO^L
4443  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4444  double precision, intent(inout) :: w(ixI^S,1:nw)
4445 
4446  double precision :: a(ixI^S,3), b(ixI^S,3), axb(ixI^S,3)
4447  integer :: idir
4448 
4449  a=0.d0
4450  b=0.d0
4451  ! for force-free field J0xB0 =0
4452  if(.not.b0field_forcefree) then
4453  ! store B0 magnetic field in b
4454  b(ixo^s,1:ndir)=block%B0(ixo^s,1:ndir,0)
4455 
4456  ! store J0 current in a
4457  do idir=7-2*ndir,3
4458  a(ixo^s,idir)=block%J0(ixo^s,idir)
4459  end do
4460  call cross_product(ixi^l,ixo^l,a,b,axb)
4461  axb(ixo^s,:)=axb(ixo^s,:)*qdt
4462  ! add J0xB0 source term in momentum equations
4463  w(ixo^s,mom(1:ndir))=w(ixo^s,mom(1:ndir))+axb(ixo^s,1:ndir)
4464  end if
4465 
4466  if(total_energy) then
4467  a=0.d0
4468  ! for free-free field -(vxB0) dot J0 =0
4469  b(ixo^s,:)=wct(ixo^s,mag(:))
4470  ! store full magnetic field B0+B1 in b
4471  if(.not.b0field_forcefree) b(ixo^s,:)=b(ixo^s,:)+block%B0(ixo^s,:,0)
4472  ! store velocity in a
4473  call mhd_get_v(wct,x,ixi^l,ixo^l,a(ixi^s,1:ndir))
4474  call cross_product(ixi^l,ixo^l,a,b,axb)
4475  axb(ixo^s,:)=axb(ixo^s,:)*qdt
4476  ! add -(vxB) dot J0 source term in energy equation
4477  do idir=7-2*ndir,3
4478  w(ixo^s,e_)=w(ixo^s,e_)-axb(ixo^s,idir)*block%J0(ixo^s,idir)
4479  end do
4480  if(mhd_ambipolar) then
4481  !reuse axb
4482  call mhd_get_jxbxb(wct,x,ixi^l,ixo^l,axb)
4483  ! source J0 * E
4484  do idir=7-2*ndim,3
4485  !set electric field in jxbxb: E=nuA * jxbxb, where nuA=-etaA/rho^2
4486  call multiplyambicoef(ixi^l,ixo^l,axb(ixi^s,idir),wct,x)
4487  w(ixo^s,e_)=w(ixo^s,e_)+axb(ixo^s,idir)*block%J0(ixo^s,idir)
4488  enddo
4489  endif
4490  end if
4491 
4492  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_B0')
4493 
4494  end subroutine add_source_b0split
4495 
4496  !> Source terms for semirelativistic MHD
4497  subroutine add_source_semirelativistic(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4499  use mod_geometry
4500 
4501  integer, intent(in) :: ixI^L, ixO^L
4502  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4503  double precision, intent(inout) :: w(ixI^S,1:nw)
4504  double precision, intent(in), optional :: wCTprim(ixI^S,1:nw)
4505 
4506  double precision :: B(ixI^S,3), v(ixI^S,3), E(ixI^S,3), divE(ixI^S)
4507  integer :: idir
4508 
4509  ! store B0 magnetic field in b
4510  if(b0field) then
4511  b(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))+block%B0(ixi^s,1:ndir,0)
4512  else
4513  b(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
4514  end if
4515  v(ixi^s,1:ndir)=wctprim(ixi^s,mom(1:ndir))
4516 
4517  call cross_product(ixi^l,ixi^l,b,v,e)
4518  ! add source term in momentum equations (1/c0^2-1/c^2)E divE
4519  call divvector(e,ixi^l,ixo^l,dive)
4520  dive(ixo^s)=qdt*(inv_squared_c0-inv_squared_c)*dive(ixo^s)
4521  do idir=1,ndir
4522  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))+e(ixo^s,idir)*dive(ixo^s)
4523  end do
4524 
4525  end subroutine add_source_semirelativistic
4526 
4527  !> Source terms for hydrodynamic energy version of MHD
4528  subroutine add_source_hydrodynamic_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4530  use mod_geometry
4531  use mod_usr_methods, only: usr_gravity
4532 
4533  integer, intent(in) :: ixI^L, ixO^L
4534  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4535  double precision, intent(inout) :: w(ixI^S,1:nw)
4536  double precision, intent(in), optional :: wCTprim(ixI^S,1:nw)
4537 
4538  double precision :: B(ixI^S,3), J(ixI^S,3), JxB(ixI^S,3)
4539  integer :: idir, idirmin, idims, ix^D
4540  double precision :: current(ixI^S,7-2*ndir:3)
4541  double precision :: bu(ixO^S,1:ndir), tmp(ixO^S), b2(ixO^S)
4542  double precision :: gravity_field(ixI^S,1:ndir), Vaoc
4543 
4544  b=0.0d0
4545  do idir = 1, ndir
4546  b(ixo^s, idir) = wct(ixo^s,mag(idir))
4547  end do
4548 
4549  call get_current(wct,ixi^l,ixo^l,idirmin,current)
4550 
4551  j=0.0d0
4552  do idir=7-2*ndir,3
4553  j(ixo^s,idir)=current(ixo^s,idir)
4554  end do
4555 
4556  ! get Lorentz force JxB
4557  call cross_product(ixi^l,ixo^l,j,b,jxb)
4558 
4559  if(mhd_semirelativistic) then
4560  ! (v . nabla) v
4561  do idir=1,ndir
4562  do idims=1,ndim
4563  call gradient(wctprim(ixi^s,mom(idir)),ixi^l,ixo^l,idims,j(ixi^s,idims))
4564  end do
4565  b(ixo^s,idir)=sum(wctprim(ixo^s,mom(1:ndir))*j(ixo^s,1:ndir),dim=ndim+1)
4566  end do
4567  ! nabla p
4568  do idir=1,ndir
4569  call gradient(wctprim(ixi^s,p_),ixi^l,ixo^l,idir,j(ixi^s,idir))
4570  end do
4571 
4572  if(mhd_gravity) then
4573  if (.not. associated(usr_gravity)) then
4574  write(*,*) "mod_usr.t: please point usr_gravity to a subroutine"
4575  write(*,*) "like the phys_gravity in mod_usr_methods.t"
4576  call mpistop("add_source_hydrodynamic_e: usr_gravity not defined")
4577  else
4578  gravity_field=0.d0
4579  call usr_gravity(ixi^l,ixo^l,wct,x,gravity_field(ixi^s,1:ndim))
4580  end if
4581  do idir=1,ndir
4582  b(ixo^s,idir)=wct(ixo^s,rho_)*(b(ixo^s,idir)-gravity_field(ixo^s,idir))+j(ixo^s,idir)-jxb(ixo^s,idir)
4583  end do
4584  else
4585  do idir=1,ndir
4586  b(ixo^s,idir)=wct(ixo^s,rho_)*b(ixo^s,idir)+j(ixo^s,idir)-jxb(ixo^s,idir)
4587  end do
4588  end if
4589 
4590  b2(ixo^s)=sum(wct(ixo^s,mag(:))**2,dim=ndim+1)
4591  tmp(ixo^s)=sqrt(b2(ixo^s))
4592  where(tmp(ixo^s)>smalldouble)
4593  tmp(ixo^s)=1.d0/tmp(ixo^s)
4594  else where
4595  tmp(ixo^s)=0.d0
4596  end where
4597  ! unit vector of magnetic field
4598  do idir=1,ndir
4599  bu(ixo^s,idir)=wct(ixo^s,mag(idir))*tmp(ixo^s)
4600  end do
4601 
4602  !b2(ixO^S)=b2(ixO^S)/w(ixO^S,rho_)*inv_squared_c
4603  !b2(ixO^S)=b2(ixO^S)/(1.d0+b2(ixO^S))
4604  {do ix^db=ixomin^db,ixomax^db\}
4605  ! Va^2/c^2
4606  vaoc=b2(ix^d)/w(ix^d,rho_)*inv_squared_c
4607  ! Va^2/c^2 / (1+Va^2/c^2)
4608  b2(ix^d)=vaoc/(1.d0+vaoc)
4609  {end do\}
4610  ! bu . F
4611  tmp(ixo^s)=sum(bu(ixo^s,1:ndir)*b(ixo^s,1:ndir),dim=ndim+1)
4612  ! Rempel 2017 ApJ 834, 10 equation (54)
4613  do idir=1,ndir
4614  j(ixo^s,idir)=b2(ixo^s)*(b(ixo^s,idir)-bu(ixo^s,idir)*tmp(ixo^s))
4615  end do
4616  ! Rempel 2017 ApJ 834, 10 equation (29) add SR force at momentum equation
4617  do idir=1,ndir
4618  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))+qdt*j(ixo^s,idir)
4619  end do
4620  ! Rempel 2017 ApJ 834, 10 equation (30) add work of Lorentz force and SR force
4621  w(ixo^s,e_)=w(ixo^s,e_)+qdt*sum(wctprim(ixo^s,mom(1:ndir))*&
4622  (jxb(ixo^s,1:ndir)+j(ixo^s,1:ndir)),dim=ndim+1)
4623  else
4624  ! add work of Lorentz force
4625  w(ixo^s,e_)=w(ixo^s,e_)+qdt*sum(wctprim(ixo^s,mom(1:ndir))*jxb(ixo^s,1:ndir),dim=ndim+1)
4626  end if
4627 
4628  end subroutine add_source_hydrodynamic_e
4629 
4630  !> Add resistive source to w within ixO Uses 3 point stencil (1 neighbour) in
4631  !> each direction, non-conservative. If the fourthorder precompiler flag is
4632  !> set, uses fourth order central difference for the laplacian. Then the
4633  !> stencil is 5 (2 neighbours).
4634  subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
4636  use mod_usr_methods
4637  use mod_geometry
4638 
4639  integer, intent(in) :: ixI^L, ixO^L
4640  double precision, intent(in) :: qdt
4641  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4642  double precision, intent(inout) :: w(ixI^S,1:nw)
4643  integer :: ixA^L,idir,jdir,kdir,idirmin,idim,jxO^L,hxO^L,ix
4644  integer :: lxO^L, kxO^L
4645 
4646  double precision :: tmp(ixI^S),tmp2(ixI^S)
4647 
4648  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
4649  double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
4650  double precision :: gradeta(ixI^S,1:ndim), Bf(ixI^S,1:ndir)
4651 
4652  ! Calculating resistive sources involve one extra layer
4653  if (mhd_4th_order) then
4654  ixa^l=ixo^l^ladd2;
4655  else
4656  ixa^l=ixo^l^ladd1;
4657  end if
4658 
4659  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
4660  call mpistop("Error in add_source_res1: Non-conforming input limits")
4661 
4662  ! Calculate current density and idirmin
4663  call get_current(wct,ixi^l,ixo^l,idirmin,current)
4664 
4665  if (mhd_eta>zero)then
4666  eta(ixa^s)=mhd_eta
4667  gradeta(ixo^s,1:ndim)=zero
4668  else
4669  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
4670  ! assumes that eta is not function of current?
4671  do idim=1,ndim
4672  call gradient(eta,ixi^l,ixo^l,idim,tmp)
4673  gradeta(ixo^s,idim)=tmp(ixo^s)
4674  end do
4675  end if
4676 
4677  if(b0field) then
4678  bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))+block%B0(ixi^s,1:ndir,0)
4679  else
4680  bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
4681  end if
4682 
4683  do idir=1,ndir
4684  ! Put B_idir into tmp2 and eta*Laplace B_idir into tmp
4685  if (mhd_4th_order) then
4686  tmp(ixo^s)=zero
4687  tmp2(ixi^s)=bf(ixi^s,idir)
4688  do idim=1,ndim
4689  lxo^l=ixo^l+2*kr(idim,^d);
4690  jxo^l=ixo^l+kr(idim,^d);
4691  hxo^l=ixo^l-kr(idim,^d);
4692  kxo^l=ixo^l-2*kr(idim,^d);
4693  tmp(ixo^s)=tmp(ixo^s)+&
4694  (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
4695  /(12.0d0 * dxlevel(idim)**2)
4696  end do
4697  else
4698  tmp(ixo^s)=zero
4699  tmp2(ixi^s)=bf(ixi^s,idir)
4700  do idim=1,ndim
4701  jxo^l=ixo^l+kr(idim,^d);
4702  hxo^l=ixo^l-kr(idim,^d);
4703  tmp(ixo^s)=tmp(ixo^s)+&
4704  (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/dxlevel(idim)**2
4705  end do
4706  end if
4707 
4708  ! Multiply by eta
4709  tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
4710 
4711  ! Subtract grad(eta) x J = eps_ijk d_j eta J_k if eta is non-constant
4712  if (mhd_eta<zero)then
4713  do jdir=1,ndim; do kdir=idirmin,3
4714  if (lvc(idir,jdir,kdir)/=0)then
4715  if (lvc(idir,jdir,kdir)==1)then
4716  tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4717  else
4718  tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4719  end if
4720  end if
4721  end do; end do
4722  end if
4723 
4724  ! Add sources related to eta*laplB-grad(eta) x J to B and e
4725  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
4726  if(total_energy) then
4727  w(ixo^s,e_)=w(ixo^s,e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
4728  end if
4729  end do ! idir
4730 
4731  if(mhd_energy) then
4732  ! de/dt+=eta*J**2
4733  tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4734  w(ixo^s,e_)=w(ixo^s,e_)+tmp(ixo^s)
4735  if(mhd_solve_eaux) then
4736  ! add eta*J**2 source term in the internal energy equation
4737  w(ixo^s,eaux_)=w(ixo^s,eaux_)+tmp(ixo^s)
4738  end if
4739  end if
4740 
4741  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_res1')
4742 
4743  end subroutine add_source_res1
4744 
4745  !> Add resistive source to w within ixO
4746  !> Uses 5 point stencil (2 neighbours) in each direction, conservative
4747  subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
4749  use mod_usr_methods
4750  use mod_geometry
4751 
4752  integer, intent(in) :: ixI^L, ixO^L
4753  double precision, intent(in) :: qdt
4754  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4755  double precision, intent(inout) :: w(ixI^S,1:nw)
4756 
4757  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
4758  double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S),curlj(ixI^S,1:3)
4759  double precision :: tmpvec(ixI^S,1:3),tmp(ixO^S)
4760  integer :: ixA^L,idir,idirmin,idirmin1
4761 
4762  ixa^l=ixo^l^ladd2;
4763 
4764  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
4765  call mpistop("Error in add_source_res2: Non-conforming input limits")
4766 
4767  ixa^l=ixo^l^ladd1;
4768  ! Calculate current density within ixL: J=curl B, thus J_i=eps_ijk*d_j B_k
4769  ! Determine exact value of idirmin while doing the loop.
4770  call get_current(wct,ixi^l,ixa^l,idirmin,current)
4771 
4772  tmpvec=zero
4773  if(mhd_eta>zero)then
4774  do idir=idirmin,3
4775  tmpvec(ixa^s,idir)=current(ixa^s,idir)*mhd_eta
4776  end do
4777  else
4778  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
4779  do idir=idirmin,3
4780  tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
4781  end do
4782  end if
4783 
4784  ! dB/dt= -curl(J*eta), thus B_i=B_i-eps_ijk d_j Jeta_k
4785  call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
4786  if(stagger_grid) then
4787  if(ndim==2.and.ndir==3) then
4788  ! if 2.5D
4789  w(ixo^s,mag(ndir)) = w(ixo^s,mag(ndir))-qdt*curlj(ixo^s,ndir)
4790  end if
4791  else
4792  w(ixo^s,mag(1:ndir)) = w(ixo^s,mag(1:ndir))-qdt*curlj(ixo^s,1:ndir)
4793  end if
4794 
4795  if(mhd_energy) then
4796  if(mhd_eta>zero)then
4797  tmp(ixo^s)=qdt*mhd_eta*sum(current(ixo^s,:)**2,dim=ndim+1)
4798  else
4799  tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4800  end if
4801  if(total_energy) then
4802  ! de/dt= +div(B x Jeta) = eta J^2 - B dot curl(eta J)
4803  ! de1/dt= eta J^2 - B1 dot curl(eta J)
4804  w(ixo^s,e_)=w(ixo^s,e_)+tmp(ixo^s)-&
4805  qdt*sum(wct(ixo^s,mag(1:ndir))*curlj(ixo^s,1:ndir),dim=ndim+1)
4806  else
4807  ! add eta*J**2 source term in the internal energy equation
4808  w(ixo^s,e_)=w(ixo^s,e_)+tmp(ixo^s)
4809  end if
4810  if(mhd_solve_eaux) then
4811  ! add eta*J**2 source term in the internal energy equation
4812  w(ixo^s,eaux_)=w(ixo^s,eaux_)+tmp(ixo^s)
4813  end if
4814  end if
4815 
4816  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_res2')
4817  end subroutine add_source_res2
4818 
4819  !> Add Hyper-resistive source to w within ixO
4820  !> Uses 9 point stencil (4 neighbours) in each direction.
4821  subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
4823  use mod_geometry
4824 
4825  integer, intent(in) :: ixI^L, ixO^L
4826  double precision, intent(in) :: qdt
4827  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4828  double precision, intent(inout) :: w(ixI^S,1:nw)
4829  !.. local ..
4830  double precision :: current(ixI^S,7-2*ndir:3)
4831  double precision :: tmpvec(ixI^S,1:3),tmpvec2(ixI^S,1:3),tmp(ixI^S),ehyper(ixI^S,1:3)
4832  integer :: ixA^L,idir,jdir,kdir,idirmin,idirmin1
4833 
4834  ixa^l=ixo^l^ladd3;
4835  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
4836  call mpistop("Error in add_source_hyperres: Non-conforming input limits")
4837 
4838  call get_current(wct,ixi^l,ixa^l,idirmin,current)
4839  tmpvec(ixa^s,1:ndir)=zero
4840  do jdir=idirmin,3
4841  tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
4842  end do
4843 
4844  ixa^l=ixo^l^ladd2;
4845  call curlvector(tmpvec,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
4846 
4847  ixa^l=ixo^l^ladd1;
4848  tmpvec(ixa^s,1:ndir)=zero
4849  call curlvector(tmpvec2,ixi^l,ixa^l,tmpvec,idirmin1,1,3)
4850  ehyper(ixa^s,1:ndir) = - tmpvec(ixa^s,1:ndir)*mhd_eta_hyper
4851 
4852  ixa^l=ixo^l;
4853  tmpvec2(ixa^s,1:ndir)=zero
4854  call curlvector(ehyper,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
4855 
4856  do idir=1,ndir
4857  w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
4858  end do
4859 
4860  if(total_energy) then
4861  ! de/dt= +div(B x Ehyper)
4862  ixa^l=ixo^l^ladd1;
4863  tmpvec2(ixa^s,1:ndir)=zero
4864  do idir=1,ndir; do jdir=1,ndir; do kdir=idirmin,3
4865  tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
4866  + lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
4867  end do; end do; end do
4868  tmp(ixo^s)=zero
4869  call divvector(tmpvec2,ixi^l,ixo^l,tmp)
4870  w(ixo^s,e_)=w(ixo^s,e_)+tmp(ixo^s)*qdt
4871  end if
4872 
4873  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_hyperres')
4874 
4875  end subroutine add_source_hyperres
4876 
4877  subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
4878  ! Add divB related sources to w within ixO
4879  ! corresponding to Dedner JCP 2002, 175, 645 _equation 24_
4880  ! giving the EGLM-MHD scheme or GLM-MHD scheme
4882  use mod_geometry
4883 
4884  integer, intent(in) :: ixI^L, ixO^L
4885  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4886  double precision, intent(inout) :: w(ixI^S,1:nw)
4887  double precision:: divb(ixI^S)
4888  integer :: idim,idir
4889  double precision :: gradPsi(ixI^S)
4890 
4891 
4892  ! dPsi/dt = - Ch^2/Cp^2 Psi
4893  if (mhd_glm_alpha < zero) then
4894  w(ixo^s,psi_) = abs(mhd_glm_alpha)*wct(ixo^s,psi_)
4895  else
4896  ! implicit update of Psi variable
4897  ! equation (27) in Mignone 2010 J. Com. Phys. 229, 2117
4898  if(slab_uniform) then
4899  w(ixo^s,psi_) = dexp(-qdt*cmax_global*mhd_glm_alpha/minval(dxlevel(:)))*w(ixo^s,psi_)
4900  else
4901  w(ixo^s,psi_) = dexp(-qdt*cmax_global*mhd_glm_alpha/minval(block%ds(ixo^s,:),dim=ndim+1))*w(ixo^s,psi_)
4902  end if
4903  end if
4904 
4905  if(mhd_glm_extended) then
4906  ! gradient of Psi
4907  if(total_energy) then
4908  do idim=1,ndim
4909  select case(typegrad)
4910  case("central")
4911  call gradient(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
4912  case("limited")
4913  call gradients(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
4914  end select
4915  ! e = e -qdt (b . grad(Psi))
4916  w(ixo^s,e_) = w(ixo^s,e_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
4917  end do
4918  end if
4919 
4920  ! We calculate now div B
4921  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
4922 
4923  ! m = m - qdt b div b
4924  do idir=1,ndir
4925  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-qdt*mhd_mag_i_all(w,ixi^l,ixo^l,idir)*divb(ixo^s)
4926  end do
4927  end if
4928 
4929  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_glm')
4930 
4931  end subroutine add_source_glm
4932 
4933  !> Add divB related sources to w within ixO corresponding to Powel
4934  subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
4936 
4937  integer, intent(in) :: ixI^L, ixO^L
4938  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4939  double precision, intent(inout) :: w(ixI^S,1:nw)
4940  double precision :: divb(ixI^S),v(ixI^S,1:ndir)
4941  integer :: idir
4942 
4943  ! We calculate now div B
4944  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
4945 
4946  ! calculate velocity
4947  call mhd_get_v(wct,x,ixi^l,ixo^l,v)
4948 
4949  if (total_energy) then
4950  ! e = e - qdt (v . b) * div b
4951  w(ixo^s,e_)=w(ixo^s,e_)-&
4952  qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=ndim+1)*divb(ixo^s)
4953  end if
4954 
4955  ! b = b - qdt v * div b
4956  do idir=1,ndir
4957  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
4958  end do
4959 
4960  ! m = m - qdt b div b
4961  do idir=1,ndir
4962  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-qdt*mhd_mag_i_all(w,ixi^l,ixo^l,idir)*divb(ixo^s)
4963  end do
4964 
4965  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_powel')
4966 
4967  end subroutine add_source_powel
4968 
4969  subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
4970  ! Add divB related sources to w within ixO
4971  ! corresponding to Janhunen, just the term in the induction equation.
4973 
4974  integer, intent(in) :: ixI^L, ixO^L
4975  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4976  double precision, intent(inout) :: w(ixI^S,1:nw)
4977  double precision :: divb(ixI^S),vel(ixI^S,1:ndir)
4978  integer :: idir
4979 
4980  ! We calculate now div B
4981  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
4982 
4983  call mhd_get_v(wct,x,ixi^l,ixo^l,vel)
4984  ! b = b - qdt v * div b
4985  do idir=1,ndir
4986  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*vel(ixo^s,idir)*divb(ixo^s)
4987  end do
4988 
4989  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_janhunen')
4990 
4991  end subroutine add_source_janhunen
4992 
4993  subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
4994  ! Add Linde's divB related sources to wnew within ixO
4996  use mod_geometry
4997 
4998  integer, intent(in) :: ixI^L, ixO^L
4999  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5000  double precision, intent(inout) :: w(ixI^S,1:nw)
5001  integer :: idim, idir, ixp^L, i^D, iside
5002  double precision :: divb(ixI^S),graddivb(ixI^S)
5003  logical, dimension(-1:1^D&) :: leveljump
5004 
5005  ! Calculate div B
5006  ixp^l=ixo^l^ladd1;
5007  call get_divb(wct,ixi^l,ixp^l,divb, mhd_divb_4thorder)
5008 
5009  ! for AMR stability, retreat one cell layer from the boarders of level jump
5010  {do i^db=-1,1\}
5011  if(i^d==0|.and.) cycle
5012  if(neighbor_type(i^d,block%igrid)==2 .or. neighbor_type(i^d,block%igrid)==4) then
5013  leveljump(i^d)=.true.
5014  else
5015  leveljump(i^d)=.false.
5016  end if
5017  {end do\}
5018 
5019  ixp^l=ixo^l;
5020  do idim=1,ndim
5021  select case(idim)
5022  {case(^d)
5023  do iside=1,2
5024  i^dd=kr(^dd,^d)*(2*iside-3);
5025  if (leveljump(i^dd)) then
5026  if (iside==1) then
5027  ixpmin^d=ixomin^d-i^d
5028  else
5029  ixpmax^d=ixomax^d-i^d
5030  end if
5031  end if
5032  end do
5033  \}
5034  end select
5035  end do
5036 
5037  ! Add Linde's diffusive terms
5038  do idim=1,ndim
5039  ! Calculate grad_idim(divb)
5040  select case(typegrad)
5041  case("central")
5042  call gradient(divb,ixi^l,ixp^l,idim,graddivb)
5043  case("limited")
5044  call gradients(divb,ixi^l,ixp^l,idim,graddivb)
5045  end select
5046 
5047  ! Multiply by Linde's eta*dt = divbdiff*(c_max*dx)*dt = divbdiff*dx**2
5048  if (slab_uniform) then
5049  graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
5050  else
5051  graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
5052  /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
5053  end if
5054 
5055  w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
5056 
5057  if (typedivbdiff=='all' .and. total_energy) then
5058  ! e += B_idim*eta*grad_idim(divb)
5059  w(ixp^s,e_)=w(ixp^s,e_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
5060  end if
5061  end do
5062 
5063  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_linde')
5064 
5065  end subroutine add_source_linde
5066 
5067  !> Calculate div B within ixO
5068  subroutine get_divb(w,ixI^L,ixO^L,divb, fourthorder)
5070  use mod_geometry
5071 
5072  integer, intent(in) :: ixi^l, ixo^l
5073  double precision, intent(in) :: w(ixi^s,1:nw)
5074  double precision, intent(inout) :: divb(ixi^s)
5075  logical, intent(in), optional :: fourthorder
5076 
5077  integer :: ixc^l, idir
5078 
5079  if(stagger_grid) then
5080  divb(ixo^s)=0.d0
5081  do idir=1,ndim
5082  ixc^l=ixo^l-kr(idir,^d);
5083  divb(ixo^s)=divb(ixo^s)+block%ws(ixo^s,idir)*block%surfaceC(ixo^s,idir)-&
5084  block%ws(ixc^s,idir)*block%surfaceC(ixc^s,idir)
5085  end do
5086  divb(ixo^s)=divb(ixo^s)/block%dvolume(ixo^s)
5087  else
5088  select case(typediv)
5089  case("central")
5090  call divvector(w(ixi^s,mag(1:ndir)),ixi^l,ixo^l,divb,fourthorder)
5091  case("limited")
5092  call divvectors(w(ixi^s,mag(1:ndir)),ixi^l,ixo^l,divb)
5093  end select
5094  end if
5095 
5096  end subroutine get_divb
5097 
5098  !> get dimensionless div B = |divB| * volume / area / |B|
5099  subroutine get_normalized_divb(w,ixI^L,ixO^L,divb)
5100 
5102 
5103  integer, intent(in) :: ixi^l, ixo^l
5104  double precision, intent(in) :: w(ixi^s,1:nw)
5105  double precision :: divb(ixi^s), dsurface(ixi^s)
5106 
5107  double precision :: invb(ixo^s)
5108  integer :: ixa^l,idims
5109 
5110  call get_divb(w,ixi^l,ixo^l,divb)
5111  invb(ixo^s)=sqrt(mhd_mag_en_all(w,ixi^l,ixo^l))
5112  where(invb(ixo^s)/=0.d0)
5113  invb(ixo^s)=1.d0/invb(ixo^s)
5114  end where
5115  if(slab_uniform) then
5116  divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/dxlevel(:))
5117  else
5118  ixamin^d=ixomin^d-1;
5119  ixamax^d=ixomax^d-1;
5120  dsurface(ixo^s)= sum(block%surfaceC(ixo^s,:),dim=ndim+1)
5121  do idims=1,ndim
5122  ixa^l=ixo^l-kr(idims,^d);
5123  dsurface(ixo^s)=dsurface(ixo^s)+block%surfaceC(ixa^s,idims)
5124  end do
5125  divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
5126  block%dvolume(ixo^s)/dsurface(ixo^s)
5127  end if
5128 
5129  end subroutine get_normalized_divb
5130 
5131  !> Calculate idirmin and the idirmin:3 components of the common current array
5132  !> make sure that dxlevel(^D) is set correctly.
5133  subroutine get_current(w,ixI^L,ixO^L,idirmin,current)
5135  use mod_geometry
5136 
5137  integer, intent(in) :: ixo^l, ixi^l
5138  double precision, intent(in) :: w(ixi^s,1:nw)
5139  integer, intent(out) :: idirmin
5140  integer :: idir, idirmin0
5141 
5142  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
5143  double precision :: current(ixi^s,7-2*ndir:3)
5144 
5145  idirmin0 = 7-2*ndir
5146 
5147  call curlvector(w(ixi^s,mag(1:ndir)),ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
5148 
5149  if(b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
5150  block%J0(ixo^s,idirmin0:3)
5151  end subroutine get_current
5152 
5153  !> If resistivity is not zero, check diffusion time limit for dt
5154  subroutine mhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
5156  use mod_usr_methods
5158  use mod_viscosity, only: viscosity_get_dt
5159  use mod_gravity, only: gravity_get_dt
5160  use mod_cak_force, only: cak_get_dt
5161 
5162  integer, intent(in) :: ixI^L, ixO^L
5163  double precision, intent(inout) :: dtnew
5164  double precision, intent(in) :: dx^D
5165  double precision, intent(in) :: w(ixI^S,1:nw)
5166  double precision, intent(in) :: x(ixI^S,1:ndim)
5167 
5168  integer :: idirmin,idim
5169  double precision :: dxarr(ndim)
5170  double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
5171 
5172  dtnew = bigdouble
5173 
5174  ^d&dxarr(^d)=dx^d;
5175  if (mhd_eta>zero)then
5176  dtnew=dtdiffpar*minval(dxarr(1:ndim))**2/mhd_eta
5177  else if (mhd_eta<zero)then
5178  call get_current(w,ixi^l,ixo^l,idirmin,current)
5179  call usr_special_resistivity(w,ixi^l,ixo^l,idirmin,x,current,eta)
5180  dtnew=bigdouble
5181  do idim=1,ndim
5182  if(slab_uniform) then
5183  dtnew=min(dtnew,&
5184  dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
5185  else
5186  dtnew=min(dtnew,&
5187  dtdiffpar/(smalldouble+maxval(eta(ixo^s)/block%ds(ixo^s,idim)**2)))
5188  end if
5189  end do
5190  end if
5191 
5192  if(mhd_eta_hyper>zero) then
5193  if(slab_uniform) then
5194  dtnew=min(dtdiffpar*minval(dxarr(1:ndim))**4/mhd_eta_hyper,dtnew)
5195  else
5196  dtnew=min(dtdiffpar*minval(block%ds(ixo^s,1:ndim))**4/mhd_eta_hyper,dtnew)
5197  end if
5198  end if
5199 
5200  if(mhd_radiative_cooling) then
5201  call cooling_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x,rc_fl)
5202  end if
5203 
5204  if(mhd_viscosity) then
5205  call viscosity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
5206  end if
5207 
5208  if(mhd_gravity) then
5209  call gravity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
5210  end if
5211 
5212  if(mhd_ambipolar_exp) then