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