MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_hd_phys.t
Go to the documentation of this file.
1 !> Hydrodynamics physics module
2 module mod_hd_phys
6  use mod_physics
7  use mod_comm_lib, only: mpistop
8  implicit none
9  private
10 
11  !> Whether an energy equation is used
12  logical, public, protected :: hd_energy = .true.
13 
14  !> Whether thermal conduction is added
15  logical, public, protected :: hd_thermal_conduction = .false.
16  type(tc_fluid), allocatable, public :: tc_fl
17  type(te_fluid), allocatable, public :: te_fl_hd
18 
19  !> Whether radiative cooling is added
20  logical, public, protected :: hd_radiative_cooling = .false.
21  type(rc_fluid), allocatable, public :: rc_fl
22 
23  !> Whether dust is added
24  logical, public, protected :: hd_dust = .false.
25 
26  !> Whether viscosity is added
27  logical, public, protected :: hd_viscosity = .false.
28 
29  !> Whether gravity is added
30  logical, public, protected :: hd_gravity = .false.
31 
32  !> Whether particles module is added
33  logical, public, protected :: hd_particles = .false.
34 
35  !> Whether rotating frame is activated
36  logical, public, protected :: hd_rotating_frame = .false.
37 
38  !> Whether CAK radiation line force is activated
39  logical, public, protected :: hd_cak_force = .false.
40 
41  !> Number of tracer species
42  integer, public, protected :: hd_n_tracer = 0
43 
44  !> Whether plasma is partially ionized
45  logical, public, protected :: hd_partial_ionization = .false.
46 
47  !> Index of the density (in the w array)
48  integer, public, protected :: rho_
49 
50  !> Indices of the momentum density
51  integer, allocatable, public, protected :: mom(:)
52 
53  !> Indices of the tracers
54  integer, allocatable, public, protected :: tracer(:)
55 
56  !> Index of the energy density (-1 if not present)
57  integer, public, protected :: e_
58 
59  !> Index of the gas pressure (-1 if not present) should equal e_
60  integer, public, protected :: p_
61 
62  !> Indices of temperature
63  integer, public, protected :: te_
64 
65  !> Index of the cutoff temperature for the TRAC method
66  integer, public, protected :: tcoff_
67 
68  !> The adiabatic index
69  double precision, public :: hd_gamma = 5.d0/3.0d0
70 
71  !> The adiabatic constant
72  double precision, public :: hd_adiab = 1.0d0
73 
74  !> The small_est allowed energy
75  double precision, protected :: small_e
76 
77  !> Whether TRAC method is used
78  logical, public, protected :: hd_trac = .false.
79  integer, public, protected :: hd_trac_type = 1
80 
81  !> Allows overruling default corner filling (for debug mode, since otherwise corner primitives fail)
82  logical, public, protected :: hd_force_diagonal = .false.
83 
84  !> Helium abundance over Hydrogen
85  double precision, public, protected :: he_abundance=0.1d0
86  !> Ionization fraction of H
87  !> H_ion_fr = H+/(H+ + H)
88  double precision, public, protected :: h_ion_fr=1d0
89  !> Ionization fraction of He
90  !> He_ion_fr = (He2+ + He+)/(He2+ + He+ + He)
91  double precision, public, protected :: he_ion_fr=1d0
92  !> Ratio of number He2+ / number He+ + He2+
93  !> He_ion_fr2 = He2+/(He2+ + He+)
94  double precision, public, protected :: he_ion_fr2=1d0
95  ! used for eq of state when it is not defined by units,
96  ! the units do not contain terms related to ionization fraction
97  ! and it is p = RR * rho * T
98  double precision, public, protected :: rr=1d0
99  ! remove the below flag and assume default value = .false.
100  ! when eq state properly implemented everywhere
101  ! and not anymore through units
102  logical, public, protected :: eq_state_units = .true.
103 
104  procedure(sub_get_pthermal), pointer :: hd_get_rfactor => null()
105  ! Public methods
106  public :: hd_phys_init
107  public :: hd_kin_en
108  public :: hd_get_pthermal
109  public :: hd_get_csound2
110  public :: hd_to_conserved
111  public :: hd_to_primitive
112  public :: hd_check_params
113  public :: hd_check_w
114 
115 contains
116 
117  !> Read this module's parameters from a file
118  subroutine hd_read_params(files)
120  character(len=*), intent(in) :: files(:)
121  integer :: n
122 
123  namelist /hd_list/ hd_energy, hd_n_tracer, hd_gamma, hd_adiab, &
128 
129  do n = 1, size(files)
130  open(unitpar, file=trim(files(n)), status="old")
131  read(unitpar, hd_list, end=111)
132 111 close(unitpar)
133  end do
134 
135  end subroutine hd_read_params
136 
137  !> Write this module's parameters to a snapsoht
138  subroutine hd_write_info(fh)
140  integer, intent(in) :: fh
141  integer, parameter :: n_par = 1
142  double precision :: values(n_par)
143  character(len=name_len) :: names(n_par)
144  integer, dimension(MPI_STATUS_SIZE) :: st
145  integer :: er
146 
147  call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
148 
149  names(1) = "gamma"
150  values(1) = hd_gamma
151  call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
152  call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
153  end subroutine hd_write_info
154 
155  !> Initialize the module
156  subroutine hd_phys_init()
160  use mod_dust, only: dust_init
161  use mod_viscosity, only: viscosity_init
162  use mod_gravity, only: gravity_init
163  use mod_particles, only: particles_init
165  use mod_cak_force, only: cak_init
169  use mod_usr_methods, only: usr_rfactor
170 
171  integer :: itr, idir
172 
173  call hd_read_params(par_files)
174 
175  physics_type = "hd"
176  phys_energy = hd_energy
177  phys_total_energy = hd_energy
178  phys_internal_e = .false.
179  phys_gamma = hd_gamma
180  phys_partial_ionization=hd_partial_ionization
181 
183  if(phys_trac) then
184  if(ndim .eq. 1) then
185  if(hd_trac_type .gt. 2) then
186  hd_trac_type=1
187  if(mype==0) write(*,*) 'WARNING: set hd_trac_type=1'
188  end if
190  else
191  phys_trac=.false.
192  if(mype==0) write(*,*) 'WARNING: set hd_trac=F when ndim>=2'
193  end if
194  end if
195 
196  ! set default gamma for polytropic/isothermal process
197  if(.not.hd_energy) then
198  if(hd_thermal_conduction) then
199  hd_thermal_conduction=.false.
200  if(mype==0) write(*,*) 'WARNING: set hd_thermal_conduction=F when hd_energy=F'
201  end if
202  if(hd_radiative_cooling) then
203  hd_radiative_cooling=.false.
204  if(mype==0) write(*,*) 'WARNING: set hd_radiative_cooling=F when hd_energy=F'
205  end if
206  end if
207  if(.not.eq_state_units) then
208  if(hd_partial_ionization) then
209  hd_partial_ionization=.false.
210  if(mype==0) write(*,*) 'WARNING: set hd_partial_ionization=F when eq_state_units=F'
211  end if
212  end if
214 
215  allocate(start_indices(number_species),stop_indices(number_species))
216 
217  ! set the index of the first flux variable for species 1
218  start_indices(1)=1
219  ! Determine flux variables
220  rho_ = var_set_rho()
221 
222  allocate(mom(ndir))
223  mom(:) = var_set_momentum(ndir)
224 
225  ! Set index of energy variable
226  if (hd_energy) then
227  e_ = var_set_energy()
228  p_ = e_
229  else
230  e_ = -1
231  p_ = -1
232  end if
233 
234  phys_get_dt => hd_get_dt
235  phys_get_cmax => hd_get_cmax
236  phys_get_a2max => hd_get_a2max
237  phys_get_tcutoff => hd_get_tcutoff
238  phys_get_cbounds => hd_get_cbounds
239  phys_get_flux => hd_get_flux
240  phys_add_source_geom => hd_add_source_geom
241  phys_add_source => hd_add_source
242  phys_to_conserved => hd_to_conserved
243  phys_to_primitive => hd_to_primitive
244  phys_check_params => hd_check_params
245  phys_check_w => hd_check_w
246  phys_get_pthermal => hd_get_pthermal
247  phys_get_v => hd_get_v
248  phys_get_rho => hd_get_rho
249  phys_write_info => hd_write_info
250  phys_handle_small_values => hd_handle_small_values
251 
252  ! Whether diagonal ghost cells are required for the physics
253  phys_req_diagonal = .false.
254 
255  ! derive units from basic units
256  call hd_physical_units()
257 
258  if (hd_dust) then
259  call dust_init(rho_, mom(:), e_)
260  endif
261 
262  if (hd_force_diagonal) then
263  ! ensure corners are filled, otherwise divide by zero when getting primitives
264  ! --> only for debug purposes
265  phys_req_diagonal = .true.
266  endif
267 
268  allocate(tracer(hd_n_tracer))
269 
270  ! Set starting index of tracers
271  do itr = 1, hd_n_tracer
272  tracer(itr) = var_set_fluxvar("trc", "trp", itr, need_bc=.false.)
273  end do
274 
275  ! set number of variables which need update ghostcells
276  nwgc=nwflux+nwaux
277 
278  ! set the index of the last flux variable for species 1
279  stop_indices(1)=nwflux
280 
281  ! set temperature as an auxiliary variable to get ionization degree
282  if(hd_partial_ionization) then
283  te_ = var_set_auxvar('Te','Te')
284  else
285  te_ = -1
286  end if
287 
288  if(hd_trac) then
289  tcoff_ = var_set_wextra()
290  iw_tcoff=tcoff_
291  else
292  tcoff_ = -1
293  end if
294 
295  ! choose Rfactor in ideal gas law
296  if(hd_partial_ionization) then
298  phys_update_temperature => hd_update_temperature
299  else if(associated(usr_rfactor)) then
300  hd_get_rfactor=>usr_rfactor
301  else
302  hd_get_rfactor=>rfactor_from_constant_ionization
303  end if
304 
305  ! initialize thermal conduction module
306  if (hd_thermal_conduction) then
307  if (.not. hd_energy) &
308  call mpistop("thermal conduction needs hd_energy=T")
309  phys_req_diagonal = .true.
310 
311  call sts_init()
313 
314  allocate(tc_fl)
319  tc_fl%get_temperature_from_conserved => hd_get_temperature_from_etot
320  tc_fl%get_temperature_from_eint => hd_get_temperature_from_eint
321  tc_fl%get_rho => hd_get_rho
322  tc_fl%e_ = e_
323  end if
324 
325  ! Initialize radiative cooling module
326  if (hd_radiative_cooling) then
327  if (.not. hd_energy) &
328  call mpistop("radiative cooling needs hd_energy=T")
330  allocate(rc_fl)
332  rc_fl%get_rho => hd_get_rho
333  rc_fl%get_pthermal => hd_get_pthermal
334  rc_fl%get_var_Rfactor => hd_get_rfactor
335  rc_fl%e_ = e_
336  rc_fl%Tcoff_ = tcoff_
337  end if
338  allocate(te_fl_hd)
339  te_fl_hd%get_rho=> hd_get_rho
340  te_fl_hd%get_pthermal=> hd_get_pthermal
341  te_fl_hd%get_var_Rfactor => hd_get_rfactor
342 {^ifthreed
343  phys_te_images => hd_te_images
344 }
345  ! Initialize viscosity module
346  if (hd_viscosity) call viscosity_init(phys_wider_stencil,phys_req_diagonal)
347 
348  ! Initialize gravity module
349  if (hd_gravity) call gravity_init()
350 
351  ! Initialize rotating_frame module
353 
354  ! Initialize CAK radiation force module
355  if (hd_cak_force) call cak_init(hd_gamma)
356 
357  ! Initialize particles module
358  if (hd_particles) then
359  call particles_init()
360  phys_req_diagonal = .true.
361  end if
362 
363  ! Check whether custom flux types have been defined
364  if (.not. allocated(flux_type)) then
365  allocate(flux_type(ndir, nw))
366  flux_type = flux_default
367  else if (any(shape(flux_type) /= [ndir, nw])) then
368  call mpistop("phys_check error: flux_type has wrong shape")
369  end if
370 
371  nvector = 1 ! No. vector vars
372  allocate(iw_vector(nvector))
373  iw_vector(1) = mom(1) - 1
374  ! initialize ionization degree table
376 
377  end subroutine hd_phys_init
378 
379 {^ifthreed
380  subroutine hd_te_images
383  select case(convert_type)
384  case('EIvtiCCmpi','EIvtuCCmpi')
386  case('ESvtiCCmpi','ESvtuCCmpi')
388  case('SIvtiCCmpi','SIvtuCCmpi')
390  case('WIvtiCCmpi','WIvtuCCmpi')
392  case default
393  call mpistop("Error in synthesize emission: Unknown convert_type")
394  end select
395  end subroutine hd_te_images
396 }
397 !!start th cond
398  ! wrappers for STS functions in thermal_conductivity module
399  ! which take as argument the tc_fluid (defined in the physics module)
400  subroutine hd_sts_set_source_tc_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
402  use mod_fix_conserve
404  integer, intent(in) :: ixI^L, ixO^L, igrid, nflux
405  double precision, intent(in) :: x(ixI^S,1:ndim)
406  double precision, intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
407  double precision, intent(in) :: my_dt
408  logical, intent(in) :: fix_conserve_at_step
409  call sts_set_source_tc_hd(ixi^l,ixo^l,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux,tc_fl)
410  end subroutine hd_sts_set_source_tc_hd
411 
412  function hd_get_tc_dt_hd(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
413  !Check diffusion time limit dt < dx_i**2/((gamma-1)*tc_k_para_i/rho)
414  !where tc_k_para_i=tc_k_para*B_i**2/B**2
415  !and T=p/rho
418 
419  integer, intent(in) :: ixi^l, ixo^l
420  double precision, intent(in) :: dx^d, x(ixi^s,1:ndim)
421  double precision, intent(in) :: w(ixi^s,1:nw)
422  double precision :: dtnew
423 
424  dtnew=get_tc_dt_hd(w,ixi^l,ixo^l,dx^d,x,tc_fl)
425  end function hd_get_tc_dt_hd
426 
427  subroutine hd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
428  ! move this in a different routine as in mhd if needed in more places
430  use mod_small_values
431 
432  integer, intent(in) :: ixI^L,ixO^L
433  double precision, intent(inout) :: w(ixI^S,1:nw)
434  double precision, intent(in) :: x(ixI^S,1:ndim)
435  integer, intent(in) :: step
436 
437  integer :: idir
438  logical :: flag(ixI^S,1:nw)
439  character(len=140) :: error_msg
440 
441  flag=.false.
442  where(w(ixo^s,e_)<small_e) flag(ixo^s,e_)=.true.
443  if(any(flag(ixo^s,e_))) then
444  select case (small_values_method)
445  case ("replace")
446  where(flag(ixo^s,e_)) w(ixo^s,e_)=small_e
447  case ("average")
448  call small_values_average(ixi^l, ixo^l, w, x, flag, e_)
449  case default
450  ! small values error shows primitive variables
451  w(ixo^s,e_)=w(ixo^s,e_)*(hd_gamma - 1.0d0)
452  do idir = 1, ndir
453  w(ixo^s, iw_mom(idir)) = w(ixo^s, iw_mom(idir))/w(ixo^s,rho_)
454  end do
455  write(error_msg,"(a,i3)") "Thermal conduction step ", step
456  call small_values_error(w, x, ixi^l, ixo^l, flag, error_msg)
457  end select
458  end if
459  end subroutine hd_tc_handle_small_e
460 
461  ! fill in tc_fluid fields from namelist
462  subroutine tc_params_read_hd(fl)
464  type(tc_fluid), intent(inout) :: fl
465  integer :: n
466  logical :: tc_saturate=.false.
467  double precision :: tc_k_para=0d0
468 
469  namelist /tc_list/ tc_saturate, tc_k_para
470 
471  do n = 1, size(par_files)
472  open(unitpar, file=trim(par_files(n)), status="old")
473  read(unitpar, tc_list, end=111)
474 111 close(unitpar)
475  end do
476  fl%tc_saturate = tc_saturate
477  fl%tc_k_para = tc_k_para
478 
479  end subroutine tc_params_read_hd
480 
481  subroutine hd_get_rho(w,x,ixI^L,ixO^L,rho)
483  integer, intent(in) :: ixI^L, ixO^L
484  double precision, intent(in) :: w(ixI^S,1:nw),x(ixI^S,1:ndim)
485  double precision, intent(out) :: rho(ixI^S)
486 
487  rho(ixo^s) = w(ixo^s,rho_)
488 
489  end subroutine hd_get_rho
490 
491 !!end th cond
492 !!rad cool
493  subroutine rc_params_read(fl)
495  use mod_constants, only: bigdouble
496  use mod_basic_types, only: std_len
497  type(rc_fluid), intent(inout) :: fl
498  integer :: n
499  ! list parameters
500  integer :: ncool = 4000
501  double precision :: cfrac=0.1d0
502 
503  !> Name of cooling curve
504  character(len=std_len) :: coolcurve='JCcorona'
505 
506  !> Name of cooling method
507  character(len=std_len) :: coolmethod='exact'
508 
509  !> Fixed temperature not lower than tlow
510  logical :: Tfix=.false.
511 
512  !> Lower limit of temperature
513  double precision :: tlow=bigdouble
514 
515  !> Add cooling source in a split way (.true.) or un-split way (.false.)
516  logical :: rc_split=.false.
517 
518 
519  namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
520 
521  do n = 1, size(par_files)
522  open(unitpar, file=trim(par_files(n)), status="old")
523  read(unitpar, rc_list, end=111)
524 111 close(unitpar)
525  end do
526 
527  fl%ncool=ncool
528  fl%coolcurve=coolcurve
529  fl%coolmethod=coolmethod
530  fl%tlow=tlow
531  fl%Tfix=tfix
532  fl%rc_split=rc_split
533  fl%cfrac=cfrac
534  end subroutine rc_params_read
535 !! end rad cool
536 
537  subroutine hd_check_params
540 
541  if (.not. hd_energy) then
542  if (hd_gamma <= 0.0d0) call mpistop ("Error: hd_gamma <= 0")
543  if (hd_adiab < 0.0d0) call mpistop ("Error: hd_adiab < 0")
545  else
546  if (hd_gamma <= 0.0d0 .or. hd_gamma == 1.0d0) &
547  call mpistop ("Error: hd_gamma <= 0 or hd_gamma == 1.0")
548  small_e = small_pressure/(hd_gamma - 1.0d0)
549  end if
550 
551  if (hd_dust) call dust_check_params()
552  if(use_imex_scheme) then
553  ! implicit dust update
554  phys_implicit_update => dust_implicit_update
555  phys_evaluate_implicit => dust_evaluate_implicit
556  endif
557 
558  end subroutine hd_check_params
559 
560  subroutine hd_physical_units
562  double precision :: mp,kB
563  double precision :: a,b
564  ! Derive scaling units
565  if(si_unit) then
566  mp=mp_si
567  kb=kb_si
568  else
569  mp=mp_cgs
570  kb=kb_cgs
571  end if
572  if(eq_state_units) then
573  a = 1d0 + 4d0 * he_abundance
574  if(hd_partial_ionization) then
575  b = 2.d0+3.d0*he_abundance
576  else
577  b = 1d0 + h_ion_fr + he_abundance*(he_ion_fr*(he_ion_fr2 + 1d0)+1d0)
578  end if
579  rr = 1d0
580  else
581  a = 1d0
582  b = 1d0
583  rr = (1d0 + h_ion_fr + he_abundance*(he_ion_fr*(he_ion_fr2 + 1d0)+1d0))/(1d0 + 4d0 * he_abundance)
584  end if
585  if(unit_density/=1.d0) then
587  else
588  ! unit of numberdensity is independent by default
590  end if
591  if(unit_velocity/=1.d0) then
594  else if(unit_pressure/=1.d0) then
597  else
598  ! unit of temperature is independent by default
601  end if
602  if(unit_time/=1.d0) then
604  else
605  ! unit of length is independent by default
607  end if
609 
610  end subroutine hd_physical_units
611 
612  !> Returns logical argument flag where values are ok
613  subroutine hd_check_w(primitive, ixI^L, ixO^L, w, flag)
615  use mod_dust, only: dust_check_w
616 
617  logical, intent(in) :: primitive
618  integer, intent(in) :: ixi^l, ixo^l
619  double precision, intent(in) :: w(ixi^s, nw)
620  logical, intent(inout) :: flag(ixi^s,1:nw)
621  double precision :: tmp(ixi^s)
622 
623  flag=.false.
624 
625  if (hd_energy) then
626  if (primitive) then
627  where(w(ixo^s, e_) < small_pressure) flag(ixo^s,e_) = .true.
628  else
629  tmp(ixo^s) = (hd_gamma - 1.0d0)*(w(ixo^s, e_) - &
630  hd_kin_en(w, ixi^l, ixo^l))
631  where(tmp(ixo^s) < small_pressure) flag(ixo^s,e_) = .true.
632  endif
633  end if
634 
635  where(w(ixo^s, rho_) < small_density) flag(ixo^s,rho_) = .true.
636 
637  if(hd_dust) call dust_check_w(ixi^l,ixo^l,w,flag)
638 
639  end subroutine hd_check_w
640 
641  !> Transform primitive variables into conservative ones
642  subroutine hd_to_conserved(ixI^L, ixO^L, w, x)
644  use mod_dust, only: dust_to_conserved
645  integer, intent(in) :: ixi^l, ixo^l
646  double precision, intent(inout) :: w(ixi^s, nw)
647  double precision, intent(in) :: x(ixi^s, 1:ndim)
648  double precision :: invgam
649  integer :: idir, itr
650 
651  !!if (fix_small_values) then
652  !! call hd_handle_small_values(.true., w, x, ixI^L, ixO^L, 'hd_to_conserved')
653  !!end if
654 
655  if (hd_energy) then
656  invgam = 1.d0/(hd_gamma - 1.0d0)
657  ! Calculate total energy from pressure and kinetic energy
658  w(ixo^s, e_) = w(ixo^s, e_) * invgam + &
659  0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) * w(ixo^s, rho_)
660  end if
661 
662  ! Convert velocity to momentum
663  do idir = 1, ndir
664  w(ixo^s, mom(idir)) = w(ixo^s, rho_) * w(ixo^s, mom(idir))
665  end do
666 
667  if (hd_dust) then
668  call dust_to_conserved(ixi^l, ixo^l, w, x)
669  end if
670 
671  end subroutine hd_to_conserved
672 
673  !> Transform conservative variables into primitive ones
674  subroutine hd_to_primitive(ixI^L, ixO^L, w, x)
676  use mod_dust, only: dust_to_primitive
677  integer, intent(in) :: ixi^l, ixo^l
678  double precision, intent(inout) :: w(ixi^s, nw)
679  double precision, intent(in) :: x(ixi^s, 1:ndim)
680  integer :: itr, idir
681  double precision :: inv_rho(ixo^s)
682 
683  if (fix_small_values) then
684  call hd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'hd_to_primitive')
685  end if
686 
687  inv_rho = 1.0d0 / w(ixo^s, rho_)
688 
689  if (hd_energy) then
690  ! Compute pressure
691  w(ixo^s, e_) = (hd_gamma - 1.0d0) * (w(ixo^s, e_) - &
692  hd_kin_en(w, ixi^l, ixo^l, inv_rho))
693  end if
694 
695  ! Convert momentum to velocity
696  do idir = 1, ndir
697  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir)) * inv_rho
698  end do
699 
700  ! Convert dust momentum to dust velocity
701  if (hd_dust) then
702  call dust_to_primitive(ixi^l, ixo^l, w, x)
703  end if
704 
705  end subroutine hd_to_primitive
706 
707  !> Transform internal energy to total energy
708  subroutine hd_ei_to_e(ixI^L,ixO^L,w,x)
710  integer, intent(in) :: ixI^L, ixO^L
711  double precision, intent(inout) :: w(ixI^S, nw)
712  double precision, intent(in) :: x(ixI^S, 1:ndim)
713 
714  ! Calculate total energy from internal and kinetic energy
715  w(ixo^s,e_)=w(ixo^s,e_)&
716  +hd_kin_en(w,ixi^l,ixo^l)
717 
718  end subroutine hd_ei_to_e
719 
720  !> Transform total energy to internal energy
721  subroutine hd_e_to_ei(ixI^L,ixO^L,w,x)
723  integer, intent(in) :: ixI^L, ixO^L
724  double precision, intent(inout) :: w(ixI^S, nw)
725  double precision, intent(in) :: x(ixI^S, 1:ndim)
726 
727  ! Calculate ei = e - ek
728  w(ixo^s,e_)=w(ixo^s,e_)&
729  -hd_kin_en(w,ixi^l,ixo^l)
730 
731  end subroutine hd_e_to_ei
732 
733  subroutine e_to_rhos(ixI^L, ixO^L, w, x)
735 
736  integer, intent(in) :: ixI^L, ixO^L
737  double precision :: w(ixI^S, nw)
738  double precision, intent(in) :: x(ixI^S, 1:ndim)
739 
740  if (hd_energy) then
741  w(ixo^s, e_) = (hd_gamma - 1.0d0) * w(ixo^s, rho_)**(1.0d0 - hd_gamma) * &
742  (w(ixo^s, e_) - hd_kin_en(w, ixi^l, ixo^l))
743  else
744  call mpistop("energy from entropy can not be used with -eos = iso !")
745  end if
746  end subroutine e_to_rhos
747 
748  subroutine rhos_to_e(ixI^L, ixO^L, w, x)
750 
751  integer, intent(in) :: ixI^L, ixO^L
752  double precision :: w(ixI^S, nw)
753  double precision, intent(in) :: x(ixI^S, 1:ndim)
754 
755  if (hd_energy) then
756  w(ixo^s, e_) = w(ixo^s, rho_)**(hd_gamma - 1.0d0) * w(ixo^s, e_) &
757  / (hd_gamma - 1.0d0) + hd_kin_en(w, ixi^l, ixo^l)
758  else
759  call mpistop("entropy from energy can not be used with -eos = iso !")
760  end if
761  end subroutine rhos_to_e
762 
763  !> Calculate v_i = m_i / rho within ixO^L
764  subroutine hd_get_v_idim(w, x, ixI^L, ixO^L, idim, v)
766  integer, intent(in) :: ixI^L, ixO^L, idim
767  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S, 1:ndim)
768  double precision, intent(out) :: v(ixI^S)
769 
770  v(ixo^s) = w(ixo^s, mom(idim)) / w(ixo^s, rho_)
771  end subroutine hd_get_v_idim
772 
773  !> Calculate velocity vector v_i = m_i / rho within ixO^L
774  subroutine hd_get_v(w,x,ixI^L,ixO^L,v)
776 
777  integer, intent(in) :: ixI^L, ixO^L
778  double precision, intent(in) :: w(ixI^S,nw), x(ixI^S,1:^ND)
779  double precision, intent(out) :: v(ixI^S,1:ndir)
780 
781  integer :: idir
782 
783  do idir=1,ndir
784  v(ixo^s,idir) = w(ixo^s, mom(idir)) / w(ixo^s, rho_)
785  end do
786 
787  end subroutine hd_get_v
788 
789  !> Calculate cmax_idim = csound + abs(v_idim) within ixO^L
790  subroutine hd_get_cmax(w, x, ixI^L, ixO^L, idim, cmax)
792  use mod_dust, only: dust_get_cmax_prim
794 
795  integer, intent(in) :: ixI^L, ixO^L, idim
796  ! w in primitive form
797  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S, 1:ndim)
798  double precision, intent(inout) :: cmax(ixI^S)
799 
800  if(hd_energy) then
801  cmax(ixo^s)=dabs(w(ixo^s,mom(idim)))+dsqrt(hd_gamma*w(ixo^s,p_)/w(ixo^s,rho_))
802  else
803  if (.not. associated(usr_set_pthermal)) then
804  cmax(ixo^s) = hd_adiab * w(ixo^s, rho_)**hd_gamma
805  else
806  call usr_set_pthermal(w,x,ixi^l,ixo^l,cmax)
807  end if
808  cmax(ixo^s)=dabs(w(ixo^s,mom(idim)))+dsqrt(hd_gamma*cmax(ixo^s)/w(ixo^s,rho_))
809  end if
810 
811  if (hd_dust) then
812  call dust_get_cmax_prim(w, x, ixi^l, ixo^l, idim, cmax)
813  end if
814  end subroutine hd_get_cmax
815 
816  subroutine hd_get_a2max(w,x,ixI^L,ixO^L,a2max)
818 
819  integer, intent(in) :: ixI^L, ixO^L
820  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
821  double precision, intent(inout) :: a2max(ndim)
822  double precision :: a2(ixI^S,ndim,nw)
823  integer :: gxO^L,hxO^L,jxO^L,kxO^L,i,j
824 
825  a2=zero
826  do i = 1,ndim
827  !> 4th order
828  hxo^l=ixo^l-kr(i,^d);
829  gxo^l=hxo^l-kr(i,^d);
830  jxo^l=ixo^l+kr(i,^d);
831  kxo^l=jxo^l+kr(i,^d);
832  a2(ixo^s,i,1:nw)=dabs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
833  -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
834  a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/dxlevel(i)**2
835  end do
836  end subroutine hd_get_a2max
837 
838  !> get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
839  subroutine hd_get_tcutoff(ixI^L,ixO^L,w,x,tco_local,Tmax_local)
841  integer, intent(in) :: ixI^L,ixO^L
842  double precision, intent(in) :: x(ixI^S,1:ndim)
843  ! in primitive form
844  double precision, intent(inout) :: w(ixI^S,1:nw)
845  double precision, intent(out) :: tco_local, Tmax_local
846 
847  double precision, parameter :: trac_delta=0.25d0
848  double precision :: tmp1(ixI^S),Te(ixI^S),lts(ixI^S)
849  double precision :: ltr(ixI^S),ltrc,ltrp,tcoff(ixI^S)
850  integer :: jxO^L,hxO^L
851  integer :: jxP^L,hxP^L,ixP^L
852  logical :: lrlt(ixI^S)
853 
854  {^ifoned
855  call hd_get_rfactor(w,x,ixi^l,ixi^l,te)
856  te(ixi^s)=w(ixi^s,p_)/(te(ixi^s)*w(ixi^s,rho_))
857 
858  tco_local=zero
859  tmax_local=maxval(te(ixo^s))
860  select case(hd_trac_type)
861  case(0)
862  w(ixi^s,tcoff_)=3.d5/unit_temperature
863  case(1)
864  hxo^l=ixo^l-1;
865  jxo^l=ixo^l+1;
866  lts(ixo^s)=0.5d0*dabs(te(jxo^s)-te(hxo^s))/te(ixo^s)
867  lrlt=.false.
868  where(lts(ixo^s) > trac_delta)
869  lrlt(ixo^s)=.true.
870  end where
871  if(any(lrlt(ixo^s))) then
872  tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
873  end if
874  case(2)
875  !> iijima et al. 2021, LTRAC method
876  ltrc=1.5d0
877  ltrp=2.5d0
878  ixp^l=ixo^l^ladd1;
879  hxo^l=ixo^l-1;
880  jxo^l=ixo^l+1;
881  hxp^l=ixp^l-1;
882  jxp^l=ixp^l+1;
883  lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
884  ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
885  w(ixo^s,tcoff_)=te(ixo^s)*&
886  (0.25*(ltr(jxo^s)+two*ltr(ixo^s)+ltr(hxo^s)))**0.4d0
887  case default
888  call mpistop("mhd_trac_type not allowed for 1D simulation")
889  end select
890  }
891  end subroutine hd_get_tcutoff
892 
893  !> Calculate cmax_idim = csound + abs(v_idim) within ixO^L
894  subroutine hd_get_cbounds(wLC, wRC, wLp, wRp, x, ixI^L, ixO^L, idim,Hspeed,cmax, cmin)
896  use mod_dust, only: dust_get_cmax
897  use mod_variables
898 
899  integer, intent(in) :: ixI^L, ixO^L, idim
900  ! conservative left and right status
901  double precision, intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
902  ! primitive left and right status
903  double precision, intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
904  double precision, intent(in) :: x(ixI^S, 1:ndim)
905  double precision, intent(inout) :: cmax(ixI^S,1:number_species)
906  double precision, intent(inout), optional :: cmin(ixI^S,1:number_species)
907  double precision, intent(in) :: Hspeed(ixI^S,1:number_species)
908 
909  double precision :: wmean(ixI^S,nw)
910  double precision, dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
911  integer :: ix^D
912 
913  select case(boundspeed)
914  case (1)
915  ! This implements formula (10.52) from "Riemann Solvers and Numerical
916  ! Methods for Fluid Dynamics" by Toro.
917 
918  tmp1(ixo^s)=dsqrt(wlp(ixo^s,rho_))
919  tmp2(ixo^s)=dsqrt(wrp(ixo^s,rho_))
920  tmp3(ixo^s)=1.d0/(dsqrt(wlp(ixo^s,rho_))+dsqrt(wrp(ixo^s,rho_)))
921  umean(ixo^s)=(wlp(ixo^s,mom(idim))*tmp1(ixo^s)+wrp(ixo^s,mom(idim))*tmp2(ixo^s))*tmp3(ixo^s)
922 
923  if(hd_energy) then
924  csoundl(ixo^s)=hd_gamma*wlp(ixo^s,p_)/wlp(ixo^s,rho_)
925  csoundr(ixo^s)=hd_gamma*wrp(ixo^s,p_)/wrp(ixo^s,rho_)
926  else
927  call hd_get_csound2(wlc,x,ixi^l,ixo^l,csoundl)
928  call hd_get_csound2(wrc,x,ixi^l,ixo^l,csoundr)
929  end if
930 
931  dmean(ixo^s) = (tmp1(ixo^s)*csoundl(ixo^s)+tmp2(ixo^s)*csoundr(ixo^s)) * &
932  tmp3(ixo^s) + 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2 * &
933  (wrp(ixo^s,mom(idim))-wlp(ixo^s,mom(idim)))**2
934 
935  dmean(ixo^s)=dsqrt(dmean(ixo^s))
936  if(present(cmin)) then
937  cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
938  cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
939  if(h_correction) then
940  {do ix^db=ixomin^db,ixomax^db\}
941  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
942  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
943  {end do\}
944  end if
945  else
946  cmax(ixo^s,1)=dabs(umean(ixo^s))+dmean(ixo^s)
947  end if
948 
949  if (hd_dust) then
950  wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
951  call dust_get_cmax(wmean, x, ixi^l, ixo^l, idim, cmax, cmin)
952  end if
953 
954  case (2)
955  wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
956  tmp1(ixo^s)=wmean(ixo^s,mom(idim))/wmean(ixo^s,rho_)
957  call hd_get_csound2(wmean,x,ixi^l,ixo^l,csoundr)
958  csoundr(ixo^s) = dsqrt(csoundr(ixo^s))
959 
960  if(present(cmin)) then
961  cmax(ixo^s,1)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
962  cmin(ixo^s,1)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
963  if(h_correction) then
964  {do ix^db=ixomin^db,ixomax^db\}
965  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
966  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
967  {end do\}
968  end if
969  else
970  cmax(ixo^s,1)=dabs(tmp1(ixo^s))+csoundr(ixo^s)
971  end if
972 
973  if (hd_dust) then
974  call dust_get_cmax(wmean, x, ixi^l, ixo^l, idim, cmax, cmin)
975  end if
976  case (3)
977  ! Miyoshi 2005 JCP 208, 315 equation (67)
978  if(hd_energy) then
979  csoundl(ixo^s)=hd_gamma*wlp(ixo^s,p_)/wlp(ixo^s,rho_)
980  csoundr(ixo^s)=hd_gamma*wrp(ixo^s,p_)/wrp(ixo^s,rho_)
981  else
982  call hd_get_csound2(wlc,x,ixi^l,ixo^l,csoundl)
983  call hd_get_csound2(wrc,x,ixi^l,ixo^l,csoundr)
984  end if
985  csoundl(ixo^s)=max(dsqrt(csoundl(ixo^s)),dsqrt(csoundr(ixo^s)))
986  if(present(cmin)) then
987  cmin(ixo^s,1)=min(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))-csoundl(ixo^s)
988  cmax(ixo^s,1)=max(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
989  if(h_correction) then
990  {do ix^db=ixomin^db,ixomax^db\}
991  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
992  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
993  {end do\}
994  end if
995  else
996  cmax(ixo^s,1)=max(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
997  end if
998  if (hd_dust) then
999  wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
1000  call dust_get_cmax(wmean, x, ixi^l, ixo^l, idim, cmax, cmin)
1001  end if
1002  end select
1003 
1004  end subroutine hd_get_cbounds
1005 
1006  !> Calculate the square of the thermal sound speed csound2 within ixO^L.
1007  !> csound2=gamma*p/rho
1008  subroutine hd_get_csound2(w,x,ixI^L,ixO^L,csound2)
1010  integer, intent(in) :: ixi^l, ixo^l
1011  double precision, intent(in) :: w(ixi^s,nw)
1012  double precision, intent(in) :: x(ixi^s,1:ndim)
1013  double precision, intent(out) :: csound2(ixi^s)
1014 
1015  call hd_get_pthermal(w,x,ixi^l,ixo^l,csound2)
1016  csound2(ixo^s)=hd_gamma*csound2(ixo^s)/w(ixo^s,rho_)
1017 
1018  end subroutine hd_get_csound2
1019 
1020  !> Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho) within ixO^L
1021  subroutine hd_get_pthermal(w, x, ixI^L, ixO^L, pth)
1025 
1026  integer, intent(in) :: ixi^l, ixo^l
1027  double precision, intent(in) :: w(ixi^s, 1:nw)
1028  double precision, intent(in) :: x(ixi^s, 1:ndim)
1029  double precision, intent(out):: pth(ixi^s)
1030  integer :: iw, ix^d
1031 
1032  if (hd_energy) then
1033  pth(ixo^s) = (hd_gamma - 1.0d0) * (w(ixo^s, e_) - &
1034  hd_kin_en(w, ixi^l, ixo^l))
1035  else
1036  if (.not. associated(usr_set_pthermal)) then
1037  pth(ixo^s) = hd_adiab * w(ixo^s, rho_)**hd_gamma
1038  else
1039  call usr_set_pthermal(w,x,ixi^l,ixo^l,pth)
1040  end if
1041  end if
1042 
1043  if (fix_small_values) then
1044  {do ix^db= ixo^lim^db\}
1045  if(pth(ix^d)<small_pressure) then
1046  pth(ix^d)=small_pressure
1047  endif
1048  {enddo^d&\}
1049  else if (check_small_values) then
1050  {do ix^db= ixo^lim^db\}
1051  if(pth(ix^d)<small_pressure) then
1052  write(*,*) "Error: small value of gas pressure",pth(ix^d),&
1053  " encountered when call hd_get_pthermal"
1054  write(*,*) "Iteration: ", it, " Time: ", global_time
1055  write(*,*) "Location: ", x(ix^d,:)
1056  write(*,*) "Cell number: ", ix^d
1057  do iw=1,nw
1058  write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
1059  end do
1060  ! use erroneous arithmetic operation to crash the run
1061  if(trace_small_values) write(*,*) dsqrt(pth(ix^d)-bigdouble)
1062  write(*,*) "Saving status at the previous time step"
1063  crash=.true.
1064  end if
1065  {enddo^d&\}
1066  end if
1067 
1068  end subroutine hd_get_pthermal
1069 
1070  !> Calculate temperature=p/rho when in e_ the total energy is stored
1071  subroutine hd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
1073  integer, intent(in) :: ixI^L, ixO^L
1074  double precision, intent(in) :: w(ixI^S, 1:nw)
1075  double precision, intent(in) :: x(ixI^S, 1:ndim)
1076  double precision, intent(out):: res(ixI^S)
1077 
1078  double precision :: R(ixI^S)
1079 
1080  call hd_get_rfactor(w,x,ixi^l,ixo^l,r)
1081  call hd_get_pthermal(w, x, ixi^l, ixo^l, res)
1082  res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,rho_))
1083  end subroutine hd_get_temperature_from_etot
1084 
1085  !> Calculate temperature=p/rho when in e_ the internal energy is stored
1086  subroutine hd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
1088  integer, intent(in) :: ixI^L, ixO^L
1089  double precision, intent(in) :: w(ixI^S, 1:nw)
1090  double precision, intent(in) :: x(ixI^S, 1:ndim)
1091  double precision, intent(out):: res(ixI^S)
1092 
1093  double precision :: R(ixI^S)
1094 
1095  call hd_get_rfactor(w,x,ixi^l,ixo^l,r)
1096  res(ixo^s) = (hd_gamma - 1.0d0) * w(ixo^s, e_)/(w(ixo^s,rho_)*r(ixo^s))
1097  end subroutine hd_get_temperature_from_eint
1098 
1099  ! Calculate flux f_idim[iw]
1100  subroutine hd_get_flux_cons(w, x, ixI^L, ixO^L, idim, f)
1102  use mod_dust, only: dust_get_flux
1103 
1104  integer, intent(in) :: ixI^L, ixO^L, idim
1105  double precision, intent(in) :: w(ixI^S, 1:nw), x(ixI^S, 1:ndim)
1106  double precision, intent(out) :: f(ixI^S, nwflux)
1107  double precision :: pth(ixI^S), v(ixI^S),frame_vel(ixI^S)
1108  integer :: idir, itr
1109 
1110  call hd_get_pthermal(w, x, ixi^l, ixo^l, pth)
1111  call hd_get_v_idim(w, x, ixi^l, ixo^l, idim, v)
1112 
1113  f(ixo^s, rho_) = v(ixo^s) * w(ixo^s, rho_)
1114 
1115  ! Momentum flux is v_i*m_i, +p in direction idim
1116  do idir = 1, ndir
1117  f(ixo^s, mom(idir)) = v(ixo^s) * w(ixo^s, mom(idir))
1118  end do
1119 
1120  f(ixo^s, mom(idim)) = f(ixo^s, mom(idim)) + pth(ixo^s)
1121 
1122  if(hd_energy) then
1123  ! Energy flux is v_i*(e + p)
1124  f(ixo^s, e_) = v(ixo^s) * (w(ixo^s, e_) + pth(ixo^s))
1125  end if
1126 
1127  do itr = 1, hd_n_tracer
1128  f(ixo^s, tracer(itr)) = v(ixo^s) * w(ixo^s, tracer(itr))
1129  end do
1130 
1131  ! Dust fluxes
1132  if (hd_dust) then
1133  call dust_get_flux(w, x, ixi^l, ixo^l, idim, f)
1134  end if
1135 
1136  end subroutine hd_get_flux_cons
1137 
1138  ! Calculate flux f_idim[iw]
1139  subroutine hd_get_flux(wC, w, x, ixI^L, ixO^L, idim, f)
1141  use mod_dust, only: dust_get_flux_prim
1142  use mod_viscosity, only: visc_get_flux_prim ! viscInDiv
1143 
1144  integer, intent(in) :: ixI^L, ixO^L, idim
1145  ! conservative w
1146  double precision, intent(in) :: wC(ixI^S, 1:nw)
1147  ! primitive w
1148  double precision, intent(in) :: w(ixI^S, 1:nw)
1149  double precision, intent(in) :: x(ixI^S, 1:ndim)
1150  double precision, intent(out) :: f(ixI^S, nwflux)
1151  double precision :: pth(ixI^S),frame_vel(ixI^S)
1152  integer :: idir, itr
1153 
1154  if (hd_energy) then
1155  pth(ixo^s) = w(ixo^s,p_)
1156  else
1157  call hd_get_pthermal(wc, x, ixi^l, ixo^l, pth)
1158  end if
1159 
1160  f(ixo^s, rho_) = w(ixo^s,mom(idim)) * w(ixo^s, rho_)
1161 
1162  ! Momentum flux is v_i*m_i, +p in direction idim
1163  do idir = 1, ndir
1164  f(ixo^s, mom(idir)) = w(ixo^s,mom(idim)) * wc(ixo^s, mom(idir))
1165  end do
1166 
1167  f(ixo^s, mom(idim)) = f(ixo^s, mom(idim)) + pth(ixo^s)
1168 
1169  if(hd_energy) then
1170  ! Energy flux is v_i*(e + p)
1171  f(ixo^s, e_) = w(ixo^s,mom(idim)) * (wc(ixo^s, e_) + w(ixo^s,p_))
1172  end if
1173 
1174  do itr = 1, hd_n_tracer
1175  f(ixo^s, tracer(itr)) = w(ixo^s,mom(idim)) * w(ixo^s, tracer(itr))
1176  end do
1177 
1178  ! Dust fluxes
1179  if (hd_dust) then
1180  call dust_get_flux_prim(w, x, ixi^l, ixo^l, idim, f)
1181  end if
1182 
1183  ! Viscosity fluxes - viscInDiv
1184  if (hd_viscosity) then
1185  call visc_get_flux_prim(w, x, ixi^l, ixo^l, idim, f, hd_energy)
1186  endif
1187 
1188  end subroutine hd_get_flux
1189 
1190  !> Add geometrical source terms to w
1191  !>
1192  !> Notice that the expressions of the geometrical terms depend only on ndir,
1193  !> not ndim. Eg, they are the same in 2.5D and in 3D, for any geometry.
1194  !>
1195  !> Ileyk : to do :
1196  !> - address the source term for the dust in case (coordinate == spherical)
1197  subroutine hd_add_source_geom(qdt, dtfactor, ixI^L, ixO^L, wCT, wprim, w, x)
1200  use mod_viscosity, only: visc_add_source_geom ! viscInDiv
1203  use mod_geometry
1204  integer, intent(in) :: ixI^L, ixO^L
1205  double precision, intent(in) :: qdt, dtfactor, x(ixI^S, 1:ndim)
1206  double precision, intent(inout) :: wCT(ixI^S, 1:nw), wprim(ixI^S,1:nw),w(ixI^S, 1:nw)
1207  ! to change and to set as a parameter in the parfile once the possibility to
1208  ! solve the equations in an angular momentum conserving form has been
1209  ! implemented (change tvdlf.t eg)
1210  double precision :: pth(ixI^S), source(ixI^S), minrho
1211  integer :: iw,idir, h1x^L{^NOONED, h2x^L}
1212  integer :: mr_,mphi_ ! Polar var. names
1213  integer :: irho, ifluid, n_fluids
1214  double precision :: exp_factor(ixI^S), del_exp_factor(ixI^S), exp_factor_primitive(ixI^S)
1215 
1216  if (hd_dust) then
1217  n_fluids = 1 + dust_n_species
1218  else
1219  n_fluids = 1
1220  end if
1221 
1222  select case (coordinate)
1223 
1224  case(cartesian_expansion)
1225  !the user provides the functions of exp_factor and del_exp_factor
1226  if(associated(usr_set_surface)) call usr_set_surface(ixi^l,x,block%dx,exp_factor,del_exp_factor,exp_factor_primitive)
1227  if(hd_energy) then
1228  source(ixo^s)=wprim(ixo^s, p_)
1229  else
1230  if(.not. associated(usr_set_pthermal)) then
1231  source(ixo^s)=hd_adiab * wprim(ixo^s, rho_)**hd_gamma
1232  else
1233  call usr_set_pthermal(wct,x,ixi^l,ixo^l,source)
1234  end if
1235  end if
1236  source(ixo^s) = source(ixo^s)*del_exp_factor(ixo^s)/exp_factor(ixo^s)
1237  w(ixo^s,mom(1)) = w(ixo^s,mom(1)) + qdt*source(ixo^s)
1238 
1239  case (cylindrical)
1240  do ifluid = 0, n_fluids-1
1241  ! s[mr]=(pthermal+mphi**2/rho)/radius
1242  if (ifluid == 0) then
1243  ! gas
1244  irho = rho_
1245  mr_ = mom(r_)
1246  if(phi_>0) mphi_ = mom(phi_)
1247  if(hd_energy) then
1248  source(ixo^s)=wprim(ixo^s, p_)
1249  else
1250  if(.not. associated(usr_set_pthermal)) then
1251  source(ixo^s)=hd_adiab * wprim(ixo^s, rho_)**hd_gamma
1252  else
1253  call usr_set_pthermal(wct,x,ixi^l,ixo^l,source)
1254  end if
1255  end if
1256  minrho = 0.0d0
1257  else
1258  ! dust : no pressure
1259  irho = dust_rho(ifluid)
1260  mr_ = dust_mom(r_, ifluid)
1261  if(phi_>0) mphi_ = dust_mom(phi_, ifluid)
1262  source(ixi^s) = zero
1263  minrho = 0.0d0
1264  end if
1265  if(phi_ > 0) then
1266  where (wct(ixo^s, irho) > minrho)
1267  source(ixo^s) = source(ixo^s) + wct(ixo^s,mphi_)*wprim(ixo^s,mphi_)
1268  w(ixo^s, mr_) = w(ixo^s, mr_) + qdt*source(ixo^s)/x(ixo^s,r_)
1269  end where
1270  ! s[mphi]=(-mphi*vr)/radius
1271  where (wct(ixo^s, irho) > minrho)
1272  source(ixo^s) = -wct(ixo^s, mphi_) * wprim(ixo^s, mr_)
1273  w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt * source(ixo^s) / x(ixo^s, r_)
1274  end where
1275  else
1276  ! s[mr]=2pthermal/radius
1277  w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * source(ixo^s) / x(ixo^s, r_)
1278  end if
1279  end do
1280  case (spherical)
1281  if (hd_dust) then
1282  call mpistop("Dust geom source terms not implemented yet with spherical geometries")
1283  end if
1284  mr_ = mom(r_)
1285  if(phi_>0) mphi_ = mom(phi_)
1286  h1x^l=ixo^l-kr(1,^d); {^nooned h2x^l=ixo^l-kr(2,^d);}
1287  if(hd_energy) then
1288  pth(ixo^s)=wprim(ixo^s, p_)
1289  else
1290  if(.not. associated(usr_set_pthermal)) then
1291  pth(ixo^s)=hd_adiab * wprim(ixo^s, rho_)**hd_gamma
1292  else
1293  call usr_set_pthermal(wct,x,ixi^l,ixo^l,pth)
1294  end if
1295  end if
1296  ! s[mr]=((vtheta**2+vphi**2)*rho+2*p)/r
1297  source(ixo^s) = pth(ixo^s) * x(ixo^s, 1) &
1298  *(block%surfaceC(ixo^s, 1) - block%surfaceC(h1x^s, 1)) &
1299  /block%dvolume(ixo^s)
1300  do idir = 2, ndir
1301  source(ixo^s) = source(ixo^s) + wprim(ixo^s, mom(idir))**2 * wprim(ixo^s, rho_)
1302  end do
1303  w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * source(ixo^s) / x(ixo^s, 1)
1304 
1305  {^nooned
1306  ! s[mtheta]=-(vr*vtheta*rho)/r+cot(theta)*(vphi**2*rho+p)/r
1307  source(ixo^s) = pth(ixo^s) * x(ixo^s, 1) &
1308  * (block%surfaceC(ixo^s, 2) - block%surfaceC(h2x^s, 2)) &
1309  / block%dvolume(ixo^s)
1310  if (ndir == 3) then
1311  source(ixo^s) = source(ixo^s) + (wprim(ixo^s, mom(3))**2 * wprim(ixo^s, rho_)) / tan(x(ixo^s, 2))
1312  end if
1313  source(ixo^s) = source(ixo^s) - (wprim(ixo^s, mom(2)) * wprim(ixo^s, mr_)) * wprim(ixo^s, rho_)
1314  w(ixo^s, mom(2)) = w(ixo^s, mom(2)) + qdt * source(ixo^s) / x(ixo^s, 1)
1315 
1316  if (ndir == 3) then
1317  ! s[mphi]=-(vphi*vr/rho)/r-cot(theta)*(vtheta*vphi/rho)/r
1318  source(ixo^s) = -(wprim(ixo^s, mom(3)) * wprim(ixo^s, mr_)) * wprim(ixo^s, rho_)&
1319  - (wprim(ixo^s, mom(2)) * wprim(ixo^s, mom(3))) * wprim(ixo^s, rho_) / tan(x(ixo^s, 2))
1320  w(ixo^s, mom(3)) = w(ixo^s, mom(3)) + qdt * source(ixo^s) / x(ixo^s, 1)
1321  end if
1322  }
1323  end select
1324 
1325  if (hd_viscosity) call visc_add_source_geom(qdt,ixi^l,ixo^l,wprim,w,x)
1326 
1327  if (hd_rotating_frame) then
1328  if (hd_dust) then
1329  call mpistop("Rotating frame not implemented yet with dust")
1330  else
1331  call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
1332  end if
1333  end if
1334 
1335  end subroutine hd_add_source_geom
1336 
1337  ! w[iw]= w[iw]+qdt*S[wCT, qtC, x] where S is the source based on wCT within ixO
1338  subroutine hd_add_source(qdt,dtfactor, ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
1343  use mod_usr_methods, only: usr_gravity
1345  use mod_cak_force, only: cak_add_source
1346 
1347  integer, intent(in) :: ixI^L, ixO^L
1348  double precision, intent(in) :: qdt, dtfactor
1349  double precision, intent(in) :: wCT(ixI^S, 1:nw),wCTprim(ixI^S,1:nw), x(ixI^S, 1:ndim)
1350  double precision, intent(inout) :: w(ixI^S, 1:nw)
1351  logical, intent(in) :: qsourcesplit
1352  logical, intent(inout) :: active
1353 
1354  double precision :: gravity_field(ixI^S, 1:ndim)
1355  integer :: idust, idim
1356 
1357  if(hd_dust .and. .not. use_imex_scheme) then
1358  call dust_add_source(qdt,ixi^l,ixo^l,wct,w,x,qsourcesplit,active)
1359  end if
1360 
1361  if(hd_radiative_cooling) then
1362  call radiative_cooling_add_source(qdt,ixi^l,ixo^l,wct,wctprim,w,x,&
1363  qsourcesplit,active, rc_fl)
1364  end if
1365 
1366  if(hd_viscosity) then
1367  call viscosity_add_source(qdt,ixi^l,ixo^l,wct,w,x,&
1368  hd_energy,qsourcesplit,active)
1369  end if
1370 
1371  if (hd_gravity) then
1372  call gravity_add_source(qdt,ixi^l,ixo^l,wct,wctprim,w,x,&
1373  hd_energy,.false.,qsourcesplit,active)
1374 
1375  if (hd_dust .and. qsourcesplit .eqv. grav_split) then
1376  active = .true.
1377 
1378  call usr_gravity(ixi^l, ixo^l, wct, x, gravity_field)
1379  do idust = 1, dust_n_species
1380  do idim = 1, ndim
1381  w(ixo^s, dust_mom(idim, idust)) = w(ixo^s, dust_mom(idim, idust)) &
1382  + qdt * gravity_field(ixo^s, idim) * wct(ixo^s, dust_rho(idust))
1383  end do
1384  end do
1385  end if
1386  end if
1387 
1388  if (hd_cak_force) then
1389  call cak_add_source(qdt,ixi^l,ixo^l,wct,w,x,hd_energy,qsourcesplit,active)
1390  end if
1391 
1392  if(hd_partial_ionization) then
1393  if(.not.qsourcesplit) then
1394  active = .true.
1395  call hd_update_temperature(ixi^l,ixo^l,wct,w,x)
1396  end if
1397  end if
1398 
1399  end subroutine hd_add_source
1400 
1401  subroutine hd_get_dt(w, ixI^L, ixO^L, dtnew, dx^D, x)
1403  use mod_dust, only: dust_get_dt
1405  use mod_viscosity, only: viscosity_get_dt
1406  use mod_gravity, only: gravity_get_dt
1407  use mod_cak_force, only: cak_get_dt
1408 
1409  integer, intent(in) :: ixI^L, ixO^L
1410  double precision, intent(in) :: dx^D, x(ixI^S, 1:^ND)
1411  double precision, intent(in) :: w(ixI^S, 1:nw)
1412  double precision, intent(inout) :: dtnew
1413 
1414  dtnew = bigdouble
1415 
1416  if(hd_dust) then
1417  call dust_get_dt(w, ixi^l, ixo^l, dtnew, dx^d, x)
1418  end if
1419 
1420  if(hd_radiative_cooling) then
1421  call cooling_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x,rc_fl)
1422  end if
1423 
1424  if(hd_viscosity) then
1425  call viscosity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
1426  end if
1427 
1428  if(hd_gravity) then
1429  call gravity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
1430  end if
1431 
1432  if (hd_cak_force) then
1433  call cak_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
1434  end if
1435 
1436  end subroutine hd_get_dt
1437 
1438  function hd_kin_en(w, ixI^L, ixO^L, inv_rho) result(ke)
1439  use mod_global_parameters, only: nw, ndim
1440  integer, intent(in) :: ixi^l, ixo^l
1441  double precision, intent(in) :: w(ixi^s, nw)
1442  double precision :: ke(ixo^s)
1443  double precision, intent(in), optional :: inv_rho(ixo^s)
1444 
1445  if (present(inv_rho)) then
1446  ke = 0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) * inv_rho
1447  else
1448  ke = 0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) / w(ixo^s, rho_)
1449  end if
1450  end function hd_kin_en
1451 
1452  function hd_inv_rho(w, ixI^L, ixO^L) result(inv_rho)
1453  use mod_global_parameters, only: nw, ndim
1454  integer, intent(in) :: ixi^l, ixo^l
1455  double precision, intent(in) :: w(ixi^s, nw)
1456  double precision :: inv_rho(ixo^s)
1457 
1458  ! Can make this more robust
1459  inv_rho = 1.0d0 / w(ixo^s, rho_)
1460  end function hd_inv_rho
1461 
1462  subroutine hd_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
1463  ! handles hydro (density,pressure,velocity) bootstrapping
1464  ! any negative dust density is flagged as well (and throws an error)
1465  ! small_values_method=replace also for dust
1467  use mod_small_values
1469  logical, intent(in) :: primitive
1470  integer, intent(in) :: ixI^L,ixO^L
1471  double precision, intent(inout) :: w(ixI^S,1:nw)
1472  double precision, intent(in) :: x(ixI^S,1:ndim)
1473  character(len=*), intent(in) :: subname
1474 
1475  integer :: n,idir
1476  logical :: flag(ixI^S,1:nw)
1477 
1478  call hd_check_w(primitive, ixi^l, ixo^l, w, flag)
1479 
1480  if (any(flag)) then
1481  select case (small_values_method)
1482  case ("replace")
1483  where(flag(ixo^s,rho_)) w(ixo^s,rho_) = small_density
1484  do idir = 1, ndir
1485  if(small_values_fix_iw(mom(idir))) then
1486  where(flag(ixo^s,rho_)) w(ixo^s, mom(idir)) = 0.0d0
1487  end if
1488  end do
1489  if(hd_energy)then
1490  if(small_values_fix_iw(e_)) then
1491  if(primitive) then
1492  where(flag(ixo^s,rho_)) w(ixo^s, p_) = small_pressure
1493  else
1494  where(flag(ixo^s,rho_)) w(ixo^s, e_) = small_e + hd_kin_en(w,ixi^l,ixo^l)
1495  endif
1496  end if
1497  endif
1498 
1499  if(hd_energy) then
1500  if(primitive) then
1501  where(flag(ixo^s,e_)) w(ixo^s,p_) = small_pressure
1502  else
1503  where(flag(ixo^s,e_))
1504  ! Add kinetic energy
1505  w(ixo^s,e_) = small_e + hd_kin_en(w,ixi^l,ixo^l)
1506  end where
1507  end if
1508  end if
1509 
1510  if(hd_dust)then
1511  do n=1,dust_n_species
1512  where(flag(ixo^s,dust_rho(n))) w(ixo^s,dust_rho(n)) = 0.0d0
1513  do idir = 1, ndir
1514  where(flag(ixo^s,dust_rho(n))) w(ixo^s,dust_mom(idir,n)) = 0.0d0
1515  enddo
1516  enddo
1517  endif
1518  case ("average")
1519  if(primitive)then
1520  ! averaging for all primitive fields, including dust
1521  call small_values_average(ixi^l, ixo^l, w, x, flag)
1522  else
1523  ! do averaging of density
1524  call small_values_average(ixi^l, ixo^l, w, x, flag, rho_)
1525  if(hd_energy) then
1526  ! do averaging of pressure
1527  w(ixi^s,p_)=(hd_gamma-1.d0)*(w(ixi^s,e_) &
1528  -0.5d0*sum(w(ixi^s, mom(:))**2, dim=ndim+1)/w(ixi^s,rho_))
1529  call small_values_average(ixi^l, ixo^l, w, x, flag, p_)
1530  w(ixi^s,e_)=w(ixi^s,p_)/(hd_gamma-1.d0) &
1531  +0.5d0*sum(w(ixi^s, mom(:))**2, dim=ndim+1)/w(ixi^s,rho_)
1532  end if
1533  if(hd_dust)then
1534  do n=1,dust_n_species
1535  where(flag(ixo^s,dust_rho(n))) w(ixo^s,dust_rho(n)) = 0.0d0
1536  do idir = 1, ndir
1537  where(flag(ixo^s,dust_rho(n))) w(ixo^s,dust_mom(idir,n)) = 0.0d0
1538  enddo
1539  enddo
1540  endif
1541  endif
1542  case default
1543  if(.not.primitive) then
1544  !convert w to primitive
1545  ! Calculate pressure = (gamma-1) * (e-ek)
1546  if(hd_energy) then
1547  w(ixo^s,p_)=(hd_gamma-1.d0)*(w(ixo^s,e_)-hd_kin_en(w,ixi^l,ixo^l))
1548  end if
1549  ! Convert gas momentum to velocity
1550  do idir = 1, ndir
1551  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))/w(ixo^s,rho_)
1552  end do
1553  end if
1554  ! NOTE: dust entries may still have conserved values here
1555  call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
1556  end select
1557  end if
1558  end subroutine hd_handle_small_values
1559 
1560  subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
1563  integer, intent(in) :: ixI^L, ixO^L
1564  double precision, intent(in) :: w(ixI^S,1:nw)
1565  double precision, intent(in) :: x(ixI^S,1:ndim)
1566  double precision, intent(out):: Rfactor(ixI^S)
1567 
1568  double precision :: iz_H(ixO^S),iz_He(ixO^S)
1569 
1570  call ionization_degree_from_temperature(ixi^l,ixo^l,w(ixi^s,te_),iz_h,iz_he)
1571  ! assume the first and second ionization of Helium have the same degree
1572  rfactor(ixo^s)=(1.d0+iz_h(ixo^s)+0.1d0*(1.d0+iz_he(ixo^s)*(1.d0+iz_he(ixo^s))))/2.3d0
1573 
1575 
1576  subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
1578  integer, intent(in) :: ixI^L, ixO^L
1579  double precision, intent(in) :: w(ixI^S,1:nw)
1580  double precision, intent(in) :: x(ixI^S,1:ndim)
1581  double precision, intent(out):: Rfactor(ixI^S)
1582 
1583  rfactor(ixo^s)=rr
1584 
1585  end subroutine rfactor_from_constant_ionization
1586 
1587  subroutine hd_update_temperature(ixI^L,ixO^L,wCT,w,x)
1590 
1591  integer, intent(in) :: ixI^L, ixO^L
1592  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
1593  double precision, intent(inout) :: w(ixI^S,1:nw)
1594 
1595  double precision :: iz_H(ixO^S),iz_He(ixO^S), pth(ixI^S)
1596 
1597  call ionization_degree_from_temperature(ixi^l,ixo^l,wct(ixi^s,te_),iz_h,iz_he)
1598 
1599  call hd_get_pthermal(w,x,ixi^l,ixo^l,pth)
1600 
1601  w(ixo^s,te_)=(2.d0+3.d0*he_abundance)*pth(ixo^s)/(w(ixo^s,rho_)*(1.d0+iz_h(ixo^s)+&
1602  he_abundance*(iz_he(ixo^s)*(iz_he(ixo^s)+1.d0)+1.d0)))
1603 
1604  end subroutine hd_update_temperature
1605 
1606 end module mod_hd_phys
Calculate w(iw)=w(iw)+qdt*SOURCE[wCT,qtC,x] within ixO for all indices iw=iwmin......
Module with basic data types used in amrvac.
integer, parameter std_len
Default length for strings.
Module to include CAK radiation line force in (magneto)hydrodynamic models Computes both the force fr...
Definition: mod_cak_force.t:27
subroutine cak_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
Check time step for total radiation contribution.
subroutine cak_init(phys_gamma)
Initialize the module.
Definition: mod_cak_force.t:98
subroutine cak_add_source(qdt, ixIL, ixOL, wCT, w, x, energy, qsourcesplit, active)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
Module for physical and numeric constants.
Definition: mod_constants.t:2
double precision, parameter bigdouble
A very large real number.
Definition: mod_constants.t:11
Module for including dust species, which interact with the gas through a drag force.
Definition: mod_dust.t:3
subroutine, public dust_check_w(ixIL, ixOL, w, flag)
Definition: mod_dust.t:194
subroutine, public dust_get_flux(w, x, ixIL, ixOL, idim, f)
Definition: mod_dust.t:252
subroutine, public dust_implicit_update(dtfactor, qdt, qtC, psb, psa)
Implicit solve of psb=psa+dtfactor*dt*F_im(psb)
Definition: mod_dust.t:652
subroutine, public dust_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
Get dt related to dust and gas stopping time (Laibe 2011)
Definition: mod_dust.t:898
subroutine, public dust_add_source(qdt, ixIL, ixOL, wCT, w, x, qsourcesplit, active)
w[iw]= w[iw]+qdt*S[wCT, x] where S is the source based on wCT within ixO
Definition: mod_dust.t:530
subroutine, public dust_to_primitive(ixIL, ixOL, w, x)
Definition: mod_dust.t:228
integer, dimension(:, :), allocatable, public, protected dust_mom
Indices of the dust momentum densities.
Definition: mod_dust.t:47
subroutine, public dust_get_cmax(w, x, ixIL, ixOL, idim, cmax, cmin)
Definition: mod_dust.t:1011
integer, public, protected dust_n_species
The number of dust species.
Definition: mod_dust.t:37
integer, dimension(:), allocatable, public, protected dust_rho
Indices of the dust densities.
Definition: mod_dust.t:44
subroutine, public dust_get_flux_prim(w, x, ixIL, ixOL, idim, f)
Definition: mod_dust.t:275
subroutine, public dust_check_params()
Definition: mod_dust.t:154
subroutine, public dust_evaluate_implicit(qtC, psa)
inplace update of psa==>F_im(psa)
Definition: mod_dust.t:583
subroutine, public dust_init(g_rho, g_mom, g_energy)
Definition: mod_dust.t:95
subroutine, public dust_get_cmax_prim(w, x, ixIL, ixOL, idim, cmax, cmin)
Definition: mod_dust.t:1035
subroutine, public dust_to_conserved(ixIL, ixOL, w, x)
Definition: mod_dust.t:208
Module for flux conservation near refinement boundaries.
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
integer coordinate
Definition: mod_geometry.t:7
integer, parameter spherical
Definition: mod_geometry.t:11
integer, parameter cylindrical
Definition: mod_geometry.t:10
integer, parameter cartesian_expansion
Definition: mod_geometry.t:12
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
logical h_correction
If true, do H-correction to fix the carbuncle problem at grid-aligned shocks.
double precision small_pressure
double precision unit_time
Physical scaling factor for time.
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
double precision global_time
The global simulation time.
double precision unit_mass
Physical scaling factor for mass.
logical use_imex_scheme
whether IMEX in use or not
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer it
Number of time steps taken.
double precision unit_numberdensity
Physical scaling factor for number density.
character(len=std_len) convert_type
Which format to use when converting.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision unit_length
Physical scaling factor for length.
logical use_particles
Use particles module or not.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer mype
The rank of the current MPI task.
double precision, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
double precision unit_velocity
Physical scaling factor for velocity.
double precision unit_temperature
Physical scaling factor for temperature.
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
logical phys_trac
Use TRAC for MHD or 1D HD.
logical fix_small_values
fix small values with average or replace methods
logical crash
Save a snapshot before crash a run met unphysical values.
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
double precision small_density
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
integer, parameter unitconvert
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
Module for including gravity in (magneto)hydrodynamics simulations.
Definition: mod_gravity.t:2
logical grav_split
source split or not
Definition: mod_gravity.t:6
subroutine gravity_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
Definition: mod_gravity.t:87
subroutine gravity_add_source(qdt, ixIL, ixOL, wCT, wCTprim, w, x, energy, rhov, qsourcesplit, active)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
Definition: mod_gravity.t:43
subroutine gravity_init()
Initialize the module.
Definition: mod_gravity.t:26
Hydrodynamics physics module.
Definition: mod_hd_phys.t:2
subroutine, public hd_check_params
Definition: mod_hd_phys.t:538
logical, public, protected hd_energy
Whether an energy equation is used.
Definition: mod_hd_phys.t:12
logical, public, protected hd_dust
Whether dust is added.
Definition: mod_hd_phys.t:24
subroutine rfactor_from_temperature_ionization(w, x, ixIL, ixOL, Rfactor)
Definition: mod_hd_phys.t:1561
integer, public, protected e_
Index of the energy density (-1 if not present)
Definition: mod_hd_phys.t:57
logical, public, protected hd_radiative_cooling
Whether radiative cooling is added.
Definition: mod_hd_phys.t:20
double precision, public, protected rr
Definition: mod_hd_phys.t:98
subroutine hd_get_flux(wC, w, x, ixIL, ixOL, idim, f)
Definition: mod_hd_phys.t:1140
double precision, public hd_gamma
The adiabatic index.
Definition: mod_hd_phys.t:69
subroutine rhos_to_e(ixIL, ixOL, w, x)
Definition: mod_hd_phys.t:749
subroutine hd_get_a2max(w, x, ixIL, ixOL, a2max)
Definition: mod_hd_phys.t:817
integer, public, protected hd_trac_type
Definition: mod_hd_phys.t:79
subroutine hd_get_flux_cons(w, x, ixIL, ixOL, idim, f)
Definition: mod_hd_phys.t:1101
logical, public, protected hd_particles
Whether particles module is added.
Definition: mod_hd_phys.t:33
subroutine hd_te_images
Definition: mod_hd_phys.t:381
type(tc_fluid), allocatable, public tc_fl
Definition: mod_hd_phys.t:16
logical, public, protected hd_viscosity
Whether viscosity is added.
Definition: mod_hd_phys.t:27
subroutine rc_params_read(fl)
Definition: mod_hd_phys.t:494
subroutine hd_get_rho(w, x, ixIL, ixOL, rho)
Definition: mod_hd_phys.t:482
subroutine, public hd_check_w(primitive, ixIL, ixOL, w, flag)
Returns logical argument flag where values are ok.
Definition: mod_hd_phys.t:614
integer, public, protected tcoff_
Index of the cutoff temperature for the TRAC method.
Definition: mod_hd_phys.t:66
subroutine hd_e_to_ei(ixIL, ixOL, w, x)
Transform total energy to internal energy.
Definition: mod_hd_phys.t:722
double precision, public, protected he_ion_fr2
Ratio of number He2+ / number He+ + He2+ He_ion_fr2 = He2+/(He2+ + He+)
Definition: mod_hd_phys.t:94
subroutine hd_sts_set_source_tc_hd(ixIL, ixOL, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux)
Definition: mod_hd_phys.t:401
subroutine hd_ei_to_e(ixIL, ixOL, w, x)
Transform internal energy to total energy.
Definition: mod_hd_phys.t:709
subroutine hd_add_source(qdt, dtfactor, ixIL, ixOL, wCT, wCTprim, w, x, qsourcesplit, active)
Definition: mod_hd_phys.t:1339
subroutine hd_write_info(fh)
Write this module's parameters to a snapsoht.
Definition: mod_hd_phys.t:139
integer, public, protected te_
Indices of temperature.
Definition: mod_hd_phys.t:63
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
Definition: mod_hd_phys.t:51
subroutine hd_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
Definition: mod_hd_phys.t:1402
double precision function hd_get_tc_dt_hd(w, ixIL, ixOL, dxD, x)
Definition: mod_hd_phys.t:413
double precision, public, protected h_ion_fr
Ionization fraction of H H_ion_fr = H+/(H+ + H)
Definition: mod_hd_phys.t:88
logical, public, protected hd_cak_force
Whether CAK radiation line force is activated.
Definition: mod_hd_phys.t:39
subroutine hd_get_cbounds(wLC, wRC, wLp, wRp, x, ixIL, ixOL, idim, Hspeed, cmax, cmin)
Calculate cmax_idim = csound + abs(v_idim) within ixO^L.
Definition: mod_hd_phys.t:895
subroutine, public hd_phys_init()
Initialize the module.
Definition: mod_hd_phys.t:157
subroutine hd_get_temperature_from_eint(w, x, ixIL, ixOL, res)
Calculate temperature=p/rho when in e_ the internal energy is stored.
Definition: mod_hd_phys.t:1087
integer, dimension(:), allocatable, public, protected tracer
Indices of the tracers.
Definition: mod_hd_phys.t:54
logical, public, protected hd_thermal_conduction
Whether thermal conduction is added.
Definition: mod_hd_phys.t:15
logical, public, protected eq_state_units
Definition: mod_hd_phys.t:102
logical, public, protected hd_force_diagonal
Allows overruling default corner filling (for debug mode, since otherwise corner primitives fail)
Definition: mod_hd_phys.t:82
subroutine, public hd_get_pthermal(w, x, ixIL, ixOL, pth)
Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho) within ixO^L.
Definition: mod_hd_phys.t:1022
double precision, public hd_adiab
The adiabatic constant.
Definition: mod_hd_phys.t:72
subroutine hd_get_tcutoff(ixIL, ixOL, w, x, tco_local, Tmax_local)
get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
Definition: mod_hd_phys.t:840
integer, public, protected rho_
Index of the density (in the w array)
Definition: mod_hd_phys.t:48
subroutine tc_params_read_hd(fl)
Definition: mod_hd_phys.t:463
subroutine hd_get_v_idim(w, x, ixIL, ixOL, idim, v)
Calculate v_i = m_i / rho within ixO^L.
Definition: mod_hd_phys.t:765
logical, public, protected hd_partial_ionization
Whether plasma is partially ionized.
Definition: mod_hd_phys.t:45
double precision, public, protected he_ion_fr
Ionization fraction of He He_ion_fr = (He2+ + He+)/(He2+ + He+ + He)
Definition: mod_hd_phys.t:91
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
Definition: mod_hd_phys.t:85
logical, public, protected hd_gravity
Whether gravity is added.
Definition: mod_hd_phys.t:30
subroutine hd_physical_units
Definition: mod_hd_phys.t:561
subroutine rfactor_from_constant_ionization(w, x, ixIL, ixOL, Rfactor)
Definition: mod_hd_phys.t:1577
subroutine, public hd_to_conserved(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
Definition: mod_hd_phys.t:643
subroutine hd_add_source_geom(qdt, dtfactor, ixIL, ixOL, wCT, wprim, w, x)
Add geometrical source terms to w.
Definition: mod_hd_phys.t:1198
type(rc_fluid), allocatable, public rc_fl
Definition: mod_hd_phys.t:21
subroutine, public hd_to_primitive(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
Definition: mod_hd_phys.t:675
logical, public, protected hd_trac
Whether TRAC method is used.
Definition: mod_hd_phys.t:78
integer, public, protected hd_n_tracer
Number of tracer species.
Definition: mod_hd_phys.t:42
type(te_fluid), allocatable, public te_fl_hd
Definition: mod_hd_phys.t:17
subroutine e_to_rhos(ixIL, ixOL, w, x)
Definition: mod_hd_phys.t:734
subroutine hd_get_cmax(w, x, ixIL, ixOL, idim, cmax)
Calculate cmax_idim = csound + abs(v_idim) within ixO^L.
Definition: mod_hd_phys.t:791
logical, public, protected hd_rotating_frame
Whether rotating frame is activated.
Definition: mod_hd_phys.t:36
subroutine hd_get_v(w, x, ixIL, ixOL, v)
Calculate velocity vector v_i = m_i / rho within ixO^L.
Definition: mod_hd_phys.t:775
double precision function, dimension(ixo^s), public hd_kin_en(w, ixIL, ixOL, inv_rho)
Definition: mod_hd_phys.t:1439
subroutine hd_update_temperature(ixIL, ixOL, wCT, w, x)
Definition: mod_hd_phys.t:1588
integer, public, protected p_
Index of the gas pressure (-1 if not present) should equal e_.
Definition: mod_hd_phys.t:60
subroutine hd_tc_handle_small_e(w, x, ixIL, ixOL, step)
Definition: mod_hd_phys.t:428
subroutine hd_get_temperature_from_etot(w, x, ixIL, ixOL, res)
Calculate temperature=p/rho when in e_ the total energy is stored.
Definition: mod_hd_phys.t:1072
subroutine hd_handle_small_values(primitive, w, x, ixIL, ixOL, subname)
Definition: mod_hd_phys.t:1463
subroutine, public hd_get_csound2(w, x, ixIL, ixOL, csound2)
Calculate the square of the thermal sound speed csound2 within ixO^L. csound2=gamma*p/rho.
Definition: mod_hd_phys.t:1009
double precision function, dimension(ixo^s) hd_inv_rho(w, ixIL, ixOL)
Definition: mod_hd_phys.t:1453
module ionization degree - get ionization degree for given temperature
subroutine ionization_degree_from_temperature(ixIL, ixOL, Te, iz_H, iz_He)
Module containing all the particle routines.
Definition: mod_particles.t:2
subroutine particles_init()
Initialize particle data and parameters.
Definition: mod_particles.t:15
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
module radiative cooling – add optically thin radiative cooling for HD and MHD
subroutine radiative_cooling_init(fl, read_params)
subroutine cooling_get_dt(w, ixIL, ixOL, dtnew, dxD, x, fl)
subroutine radiative_cooling_init_params(phys_gamma, He_abund)
Radiative cooling initialization.
subroutine radiative_cooling_add_source(qdt, ixIL, ixOL, wCT, wCTprim, w, x, qsourcesplit, active, fl)
Module for including rotating frame in (magneto)hydrodynamics simulations The rotation vector is assu...
subroutine rotating_frame_add_source(qdt, dtfactor, ixIL, ixOL, wCT, w, x)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
subroutine rotating_frame_init()
Initialize the module.
Module for handling problematic values in simulations, such as negative pressures.
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
subroutine, public small_values_average(ixIL, ixOL, w, x, w_flag, windex)
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
subroutine, public small_values_error(wprim, x, ixIL, ixOL, w_flag, subname)
character(len=20), public small_values_method
How to handle small values.
Generic supertimestepping method 1) in amrvac.par in sts_list set the following parameters which have...
subroutine, public add_sts_method(sts_getdt, sts_set_sources, startVar, nflux, startwbc, nwbc, evolve_B)
subroutine which added programatically a term to be calculated using STS Params: sts_getdt function c...
subroutine, public set_conversion_methods_to_head(sts_before_first_cycle, sts_after_last_cycle)
Set the hooks called before the first cycle and after the last cycle in the STS update This method sh...
subroutine, public set_error_handling_to_head(sts_error_handling)
Set the hook of error handling in the STS update. This method is called before updating the BC....
subroutine, public sts_init()
Initialize sts module.
Thermal conduction for HD and MHD or RHD and RMHD or twofl (plasma-neutral) module Adaptation of mod_...
subroutine, public sts_set_source_tc_hd(ixIL, ixOL, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
subroutine, public tc_get_hd_params(fl, read_hd_params)
Init TC coefficients: HD case.
subroutine tc_init_params(phys_gamma)
double precision function, public get_tc_dt_hd(w, ixIL, ixOL, dxD, x, fl)
Get the explicit timestep for the TC (hd implementation)
subroutine get_euv_image(qunit, fl)
subroutine get_sxr_image(qunit, fl)
subroutine get_euv_spectrum(qunit, fl)
subroutine get_whitelight_image(qunit, fl)
Module with all the methods that users can customize in AMRVAC.
procedure(rfactor), pointer usr_rfactor
procedure(set_surface), pointer usr_set_surface
procedure(phys_gravity), pointer usr_gravity
procedure(hd_pthermal), pointer usr_set_pthermal
The module add viscous source terms and check time step.
Definition: mod_viscosity.t:10
subroutine viscosity_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
subroutine viscosity_add_source(qdt, ixIL, ixOL, wCT, w, x, energy, qsourcesplit, active)
Definition: mod_viscosity.t:90
subroutine viscosity_init(phys_wider_stencil, phys_req_diagonal)
Initialize the module.
Definition: mod_viscosity.t:59
subroutine visc_add_source_geom(qdt, ixIL, ixOL, wCT, w, x)
subroutine, public visc_get_flux_prim(w, x, ixIL, ixOL, idim, f, energy)