MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_rhd_phys.t
Go to the documentation of this file.
1!> Radiation-Hydrodynamics physics module
2!> Module aims at solving the Hydrodynamic equations toghether with
3!> the zeroth moment of the radiative transfer equation. A closure is
4!> provided by the flux limited diffusion (FLD)-approximation in the mod_fld.t module. See
5!> [1]Moens, N., Sundqvist, J. O., El Mellah, I., Poniatowski, L., Teunissen, J., and Keppens, R.,
6!> “Radiation-hydrodynamics with MPI-AMRVAC . Flux-limited diffusion”,
7!> <i>Astronomy and Astrophysics</i>, vol. 657, 2022. doi:10.1051/0004-6361/202141023.
8!> For more information.
9!> Another possible closure in the works is the anisotropic flux limited diffusion approximation (AFLD) described in mod_afld.t.
10
15 use mod_physics
16 use mod_comm_lib, only: mpistop
17 implicit none
18 private
19
20 !> Whether an energy equation is used
21 logical, public, protected :: rhd_energy = .true.
22
23 !> Whether thermal conduction is added
24 logical, public, protected :: rhd_thermal_conduction = .false.
25 type(tc_fluid), allocatable, public :: tc_fl
26 type(te_fluid), allocatable, public :: te_fl_rhd
27
28 !> Whether radiative cooling is added
29 logical, public, protected :: rhd_radiative_cooling = .false.
30 type(rc_fluid), allocatable, public :: rc_fl
31
32 !> Whether dust is added
33 logical, public, protected :: rhd_dust = .false.
34
35 !> Whether viscosity is added
36 logical, public, protected :: rhd_viscosity = .false.
37
38 !> Whether gravity is added
39 logical, public, protected :: rhd_gravity = .false.
40
41 !> Whether particles module is added
42 logical, public, protected :: rhd_particles = .false.
43
44 !> Whether rotating frame is activated
45 logical, public, protected :: rhd_rotating_frame = .false.
46
47 !> Number of tracer species
48 integer, public, protected :: rhd_n_tracer = 0
49
50 !> Index of the density (in the w array)
51 integer, public, protected :: rho_
52
53 !> Indices of the momentum density
54 integer, allocatable, public, protected :: mom(:)
55
56 !> Indices of the tracers
57 integer, allocatable, public, protected :: tracer(:)
58
59 !> Index of the energy density (-1 if not present)
60 integer, public, protected :: e_
61
62 !> Index of the gas pressure (-1 if not present) should equal e_
63 integer, public, protected :: p_
64
65 !> Index of the radiation energy
66 integer, public, protected :: r_e
67
68 !> Indices of temperature
69 integer, public, protected :: te_
70
71 !> Index of the cutoff temperature for the TRAC method
72 integer, public, protected :: tcoff_
73
74 !> The adiabatic index
75 double precision, public :: rhd_gamma = 5.d0/3.0d0
76
77 !> The adiabatic constant
78 double precision, public :: rhd_adiab = 1.0d0
79
80 !> The small_est allowed energy
81 double precision, protected :: small_e
82
83 !> The smallest allowed radiation energy
84 double precision, public, protected :: small_r_e = 0.d0
85
86 !> Whether TRAC method is used
87 logical, public, protected :: rhd_trac = .false.
88 integer, public, protected :: rhd_trac_type = 1
89
90 !> Helium abundance over Hydrogen
91 double precision, public, protected :: he_abundance=0.1d0
92
93 !> Formalism to treat radiation
94 character(len=8), public :: rhd_radiation_formalism = 'fld'
95
96 !> In the case of no rhd_energy, how to compute pressure
97 character(len=8), public :: rhd_pressure = 'Trad'
98
99 !> Treat radiation fld_Rad_force
100 logical, public, protected :: rhd_radiation_force = .true.
101
102 !> Treat radiation-gas energy interaction
103 logical, public, protected :: rhd_energy_interact = .true.
104
105 !> Treat radiation energy diffusion
106 logical, public, protected :: rhd_radiation_diffusion = .true.
107
108 !> Treat radiation advection
109 logical, public, protected :: rhd_radiation_advection = .true.
110
111 !> Whether plasma is partially ionized
112 logical, public, protected :: rhd_partial_ionization = .false.
113
114 !> Do a running mean over the radiation pressure when determining dt
115 logical, protected :: radio_acoustic_filter = .false.
116 integer, protected :: size_ra_filter = 1
117
118 !> kb/(m_p mu)* 1/a_rad**4,
119 double precision, public :: kbmpmua4
120
121 !> Use the speed of light to calculate the timestep, usefull for debugging
122 logical :: dt_c = .false.
123
124 !> Ionization fraction of H
125 !> H_ion_fr = H+/(H+ + H)
126 double precision, public, protected :: h_ion_fr=1d0
127 !> Ionization fraction of He
128 !> He_ion_fr = (He2+ + He+)/(He2+ + He+ + He)
129 double precision, public, protected :: he_ion_fr=1d0
130 !> Ratio of number He2+ / number He+ + He2+
131 !> He_ion_fr2 = He2+/(He2+ + He+)
132 double precision, public, protected :: he_ion_fr2=1d0
133 ! used for eq of state when it is not defined by units,
134 ! the units do not contain terms related to ionization fraction
135 ! and it is p = RR * rho * T
136 double precision, public, protected :: rr=1d0
137 ! remove the below flag and assume default value = .false.
138 ! when eq state properly implemented everywhere
139 ! and not anymore through units
140 logical, public, protected :: eq_state_units = .true.
141
142 procedure(sub_get_pthermal), pointer :: rhd_get_rfactor => null()
143 ! Public methods
144 public :: rhd_phys_init
145 public :: rhd_kin_en
146 public :: rhd_get_pthermal
147 public :: rhd_get_pradiation
148 public :: rhd_get_ptot
149 public :: rhd_get_csound2
150 public :: rhd_to_conserved
151 public :: rhd_to_primitive
152 public :: rhd_check_params
153 public :: rhd_check_w
154 public :: rhd_get_tgas
155 public :: rhd_get_trad
156 public :: rhd_set_mg_bounds
157
158contains
159
160 !> Read this module's parameters from a file
161 subroutine rhd_read_params(files)
163 character(len=*), intent(in) :: files(:)
164 integer :: n
165
166 namelist /rhd_list/ rhd_energy, rhd_pressure, rhd_n_tracer, rhd_gamma, rhd_adiab, &
172 rhd_radiation_advection, radio_acoustic_filter, size_ra_filter, dt_c, rhd_partial_ionization
173
174 do n = 1, size(files)
175 open(unitpar, file=trim(files(n)), status="old")
176 read(unitpar, rhd_list, end=111)
177111 close(unitpar)
178 end do
179
180 end subroutine rhd_read_params
181
182 !> Write this module's parameters to a snapsoht
183 subroutine rhd_write_info(fh)
185 integer, intent(in) :: fh
186 integer, parameter :: n_par = 1
187 double precision :: values(n_par)
188 character(len=name_len) :: names(n_par)
189 integer, dimension(MPI_STATUS_SIZE) :: st
190 integer :: er
191
192 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
193
194 names(1) = "gamma"
195 values(1) = rhd_gamma
196 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
197 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
198 end subroutine rhd_write_info
199
200 !> Initialize the module
201 subroutine rhd_phys_init()
205 use mod_dust, only: dust_init
207 use mod_gravity, only: gravity_init
210 use mod_fld
211 use mod_afld
216
217 integer :: itr, idir
218
219 call rhd_read_params(par_files)
220
221 physics_type = "rhd"
222 phys_energy = rhd_energy
223 phys_total_energy = rhd_energy
224 phys_gamma = rhd_gamma
225
227 if(phys_trac) then
228 if(ndim .eq. 1) then
229 if(rhd_trac_type .gt. 2) then
231 if(mype==0) write(*,*) 'WARNING: set rhd_trac_type=1'
232 end if
234 else
235 phys_trac=.false.
236 if(mype==0) write(*,*) 'WARNING: set rhd_trac=F when ndim>=2'
237 end if
238 end if
239
240 ! set default gamma for polytropic/isothermal process
241 if(.not.rhd_energy) then
244 if(mype==0) write(*,*) 'WARNING: set rhd_thermal_conduction=F when rhd_energy=F'
245 end if
246 if(rhd_radiative_cooling) then
248 if(mype==0) write(*,*) 'WARNING: set rhd_radiative_cooling=F when rhd_energy=F'
249 end if
250 end if
251 if(.not.eq_state_units) then
254 if(mype==0) write(*,*) 'WARNING: set rhd_partial_ionization=F when eq_state_units=F'
255 end if
256 end if
258
259 allocate(start_indices(number_species),stop_indices(number_species))
260
261 ! set the index of the first flux variable for species 1
262 start_indices(1)=1
263 ! Determine flux variables
264 rho_ = var_set_rho()
265
266 allocate(mom(ndir))
267 mom(:) = var_set_momentum(ndir)
268
269 ! Set index of energy variable
270 if (rhd_energy) then
271 e_ = var_set_energy()
272 p_ = e_
273 else
274 e_ = -1
275 p_ = -1
276 end if
277
278 !> set radiation energy
279 r_e = var_set_radiation_energy()
280
281 ! set temperature as an auxiliary variable to get ionization degree
283 te_ = var_set_auxvar('Te','Te')
284 else
285 te_ = -1
286 end if
287
288 phys_get_dt => rhd_get_dt
289 phys_get_cmax => rhd_get_cmax
290 phys_get_a2max => rhd_get_a2max
291 phys_get_tcutoff => rhd_get_tcutoff
292 phys_get_cbounds => rhd_get_cbounds
293 phys_get_flux => rhd_get_flux
294 phys_add_source_geom => rhd_add_source_geom
295 phys_add_source => rhd_add_source
296 phys_to_conserved => rhd_to_conserved
297 phys_to_primitive => rhd_to_primitive
298 ! phys_ei_to_e => rhd_ei_to_e
299 ! phys_e_to_ei => rhd_e_to_ei
300 phys_check_params => rhd_check_params
301 phys_check_w => rhd_check_w
302 phys_get_pthermal => rhd_get_pthermal
303 phys_write_info => rhd_write_info
304 phys_handle_small_values => rhd_handle_small_values
305 phys_set_mg_bounds => rhd_set_mg_bounds
306 phys_get_trad => rhd_get_trad
307 phys_get_tgas => rhd_get_tgas
308
309 ! derive units from basic units
310 call rhd_physical_units()
311
312 if (rhd_dust) then
313 call dust_init(rho_, mom(:), e_)
314 endif
315
316 !> Initiate radiation-closure module
317 select case (rhd_radiation_formalism)
318 case('fld')
320 case('afld')
321 call afld_init(he_abundance, rhd_radiation_diffusion, rhd_gamma)
322 case default
323 call mpistop('Radiation formalism unknown')
324 end select
325
326 allocate(tracer(rhd_n_tracer))
327
328 ! Set starting index of tracers
329 do itr = 1, rhd_n_tracer
330 tracer(itr) = var_set_fluxvar("trc", "trp", itr, need_bc=.false.)
331 end do
332
333 ! set number of variables which need update ghostcells
334 nwgc=nwflux+nwaux
335
336 ! set the index of the last flux variable for species 1
337 stop_indices(1)=nwflux
338
339 if(rhd_trac) then
340 tcoff_ = var_set_wextra()
341 iw_tcoff=tcoff_
342 else
343 tcoff_ = -1
344 end if
345
346 ! choose Rfactor in ideal gas law
348 rhd_get_rfactor=>rfactor_from_temperature_ionization
349 phys_update_temperature => rhd_update_temperature
350 else if(associated(usr_rfactor)) then
351 rhd_get_rfactor=>usr_rfactor
352 else
353 rhd_get_rfactor=>rfactor_from_constant_ionization
354 end if
355
356 ! initialize thermal conduction module
357 if (rhd_thermal_conduction) then
358 if (.not. rhd_energy) &
359 call mpistop("thermal conduction needs rhd_energy=T")
360
361 call sts_init()
363
364 allocate(tc_fl)
365 call tc_get_hd_params(tc_fl,tc_params_read_rhd)
366 tc_fl%get_temperature_from_conserved => rhd_get_temperature_from_etot
367 call add_sts_method(rhd_get_tc_dt_rhd,rhd_sts_set_source_tc_rhd,e_,1,e_,1,.false.)
368 call set_conversion_methods_to_head(rhd_e_to_ei, rhd_ei_to_e)
369 call set_error_handling_to_head(rhd_tc_handle_small_e)
370 tc_fl%get_temperature_from_eint => rhd_get_temperature_from_eint
371 tc_fl%get_rho => rhd_get_rho
372 tc_fl%e_ = e_
373 end if
374
375 ! Initialize radiative cooling module
376 if (rhd_radiative_cooling) then
377 if (.not. rhd_energy) &
378 call mpistop("radiative cooling needs rhd_energy=T")
380 allocate(rc_fl)
381 call radiative_cooling_init(rc_fl,rc_params_read)
382 rc_fl%get_rho => rhd_get_rho
383 rc_fl%get_pthermal => rhd_get_pthermal
384 rc_fl%get_var_Rfactor => rhd_get_rfactor
385 rc_fl%e_ = e_
386 rc_fl%Tcoff_ = tcoff_
387 end if
388 allocate(te_fl_rhd)
389 te_fl_rhd%get_rho=> rhd_get_rho
390 te_fl_rhd%get_pthermal=> rhd_get_pthermal
391 te_fl_rhd%get_var_Rfactor => rhd_get_rfactor
392{^ifthreed
393 phys_te_images => rhd_te_images
394}
395 ! Initialize viscosity module
396 if (rhd_viscosity) call viscosity_init(phys_wider_stencil)
397
398 ! Initialize gravity module
399 if (rhd_gravity) call gravity_init()
400
401 ! Initialize rotating_frame module
403
404 ! Initialize particles module
405 if (rhd_particles) then
406 call particles_init()
407 end if
408
409 ! Check whether custom flux types have been defined
410 if (.not. allocated(flux_type)) then
411 allocate(flux_type(ndir, nw))
412 flux_type = flux_default
413 else if (any(shape(flux_type) /= [ndir, nw])) then
414 call mpistop("phys_check error: flux_type has wrong shape")
415 end if
416
417 nvector = 1 ! No. vector vars
418 allocate(iw_vector(nvector))
419 iw_vector(1) = mom(1) - 1
420
421 !> Usefull constante
422 kbmpmua4 = unit_pressure**(-3.d0/4.d0)*unit_density*const_kb/(const_mp*fld_mu)*const_rad_a**(-1.d0/4.d0)
423 ! initialize ionization degree table
425 end subroutine rhd_phys_init
426
427{^ifthreed
428 subroutine rhd_te_images()
431 select case(convert_type)
432 case('EIvtiCCmpi','EIvtuCCmpi')
434 case('ESvtiCCmpi','ESvtuCCmpi')
436 case('SIvtiCCmpi','SIvtuCCmpi')
438 case default
439 call mpistop("Error in synthesize emission: Unknown convert_type")
440 end select
441 end subroutine rhd_te_images
442}
443!!start th cond
444 ! wrappers for STS functions in thermal_conductivity module
445 ! which take as argument the tc_fluid (defined in the physics module)
446 subroutine rhd_sts_set_source_tc_rhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
450 integer, intent(in) :: ixi^l, ixo^l, igrid, nflux
451 double precision, intent(in) :: x(ixi^s,1:ndim)
452 double precision, intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
453 double precision, intent(in) :: my_dt
454 logical, intent(in) :: fix_conserve_at_step
455 call sts_set_source_tc_hd(ixi^l,ixo^l,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux,tc_fl)
456 end subroutine rhd_sts_set_source_tc_rhd
457
458
459 function rhd_get_tc_dt_rhd(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
460 !Check diffusion time limit dt < dx_i**2/((gamma-1)*tc_k_para_i/rho)
461 !where tc_k_para_i=tc_k_para*B_i**2/B**2
462 !and T=p/rho
465
466 integer, intent(in) :: ixi^l, ixo^l
467 double precision, intent(in) :: dx^d, x(ixi^s,1:ndim)
468 double precision, intent(in) :: w(ixi^s,1:nw)
469 double precision :: dtnew
470
471 dtnew=get_tc_dt_hd(w,ixi^l,ixo^l,dx^d,x,tc_fl)
472 end function rhd_get_tc_dt_rhd
473
474
475 subroutine rhd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
476 ! move this in a different routine as in mhd if needed in more places
479
480 integer, intent(in) :: ixi^l,ixo^l
481 double precision, intent(inout) :: w(ixi^s,1:nw)
482 double precision, intent(in) :: x(ixi^s,1:ndim)
483 integer, intent(in) :: step
484
485 integer :: idir
486 logical :: flag(ixi^s,1:nw)
487 character(len=140) :: error_msg
488
489 flag=.false.
490 where(w(ixo^s,e_)<small_e) flag(ixo^s,e_)=.true.
491 if(any(flag(ixo^s,e_))) then
492 select case (small_values_method)
493 case ("replace")
494 where(flag(ixo^s,e_)) w(ixo^s,e_)=small_e
495 case ("average")
496 call small_values_average(ixi^l, ixo^l, w, x, flag, e_)
497 case default
498 ! small values error shows primitive variables
499 w(ixo^s,e_)=w(ixo^s,e_)*(rhd_gamma - 1.0d0)
500 do idir = 1, ndir
501 w(ixo^s, iw_mom(idir)) = w(ixo^s, iw_mom(idir))/w(ixo^s,rho_)
502 end do
503 write(error_msg,"(a,i3)") "Thermal conduction step ", step
504 call small_values_error(w, x, ixi^l, ixo^l, flag, error_msg)
505 end select
506 end if
507 end subroutine rhd_tc_handle_small_e
508
509 ! fill in tc_fluid fields from namelist
510 subroutine tc_params_read_rhd(fl)
513 type(tc_fluid), intent(inout) :: fl
514 integer :: n
515 logical :: tc_saturate=.false.
516 double precision :: tc_k_para=0d0
517
518 namelist /tc_list/ tc_saturate, tc_k_para
519
520 do n = 1, size(par_files)
521 open(unitpar, file=trim(par_files(n)), status="old")
522 read(unitpar, tc_list, end=111)
523111 close(unitpar)
524 end do
525 fl%tc_saturate = tc_saturate
526 fl%tc_k_para = tc_k_para
527
528 end subroutine tc_params_read_rhd
529
530 subroutine rhd_get_rho(w,x,ixI^L,ixO^L,rho)
532 integer, intent(in) :: ixi^l, ixo^l
533 double precision, intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:ndim)
534 double precision, intent(out) :: rho(ixi^s)
535
536 rho(ixo^s) = w(ixo^s,rho_)
537
538 end subroutine rhd_get_rho
539
540!!end th cond
541!!rad cool
542 subroutine rc_params_read(fl)
544 use mod_constants, only: bigdouble
545 use mod_basic_types, only: std_len
546 type(rc_fluid), intent(inout) :: fl
547 integer :: n
548 ! list parameters
549 integer :: ncool = 4000
550 double precision :: cfrac=0.1d0
551
552 !> Name of cooling curve
553 character(len=std_len) :: coolcurve='JCcorona'
554
555 !> Name of cooling method
556 character(len=std_len) :: coolmethod='exact'
557
558 !> Fixed temperature not lower than tlow
559 logical :: tfix=.false.
560
561 !> Lower limit of temperature
562 double precision :: tlow=bigdouble
563
564 !> Add cooling source in a split way (.true.) or un-split way (.false.)
565 logical :: rc_split=.false.
566
567
568 namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
569
570 do n = 1, size(par_files)
571 open(unitpar, file=trim(par_files(n)), status="old")
572 read(unitpar, rc_list, end=111)
573111 close(unitpar)
574 end do
575
576 fl%ncool=ncool
577 fl%coolcurve=coolcurve
578 fl%coolmethod=coolmethod
579 fl%tlow=tlow
580 fl%Tfix=tfix
581 fl%rc_split=rc_split
582 fl%cfrac=cfrac
583 end subroutine rc_params_read
584!! end rad cool
585
588 use mod_dust, only: dust_check_params
589
590 if (.not. rhd_energy .and. rhd_pressure == 'adiabatic') then
591 if (rhd_gamma <= 0.0d0) call mpistop ("Error: rhd_gamma <= 0")
592 if (rhd_adiab < 0.0d0) call mpistop ("Error: rhd_adiab < 0")
594 elseif (rhd_pressure == 'Tcond') then
595 small_pressure = smalldouble
596 else
597 if (rhd_gamma <= 0.0d0 .or. rhd_gamma == 1.0d0) &
598 call mpistop ("Error: rhd_gamma <= 0 or rhd_gamma == 1.0")
599 small_e = small_pressure/(rhd_gamma - 1.0d0)
600 end if
601
603
604 if (rhd_dust) call dust_check_params()
605
606 ! if (rhd_radiation_diffusion .and. (.not. use_imex_scheme) ) &
607 if (rhd_radiation_diffusion .and. (.not. use_imex_scheme) .and. ((rhd_radiation_formalism .eq. 'fld') .or. (rhd_radiation_formalism .eq. 'afld')) ) &
608 call mpistop("Use an IMEX scheme when doing FLD")
609
611
612 end subroutine rhd_check_params
613
614 !> Set the boundaries for the diffusion of E
619
620 integer :: ib
621
622 ! Set boundary conditions for the multigrid solver
623 do ib = 1, 2*ndim
624 select case (typeboundary(r_e, ib))
625 case (bc_symm)
626 ! d/dx u = 0
627 mg%bc(ib, mg_iphi)%bc_type = mg_bc_neumann
628 mg%bc(ib, mg_iphi)%bc_value = 0.0_dp
629 case (bc_asymm)
630 ! u = 0
631 mg%bc(ib, mg_iphi)%bc_type = mg_bc_dirichlet
632 mg%bc(ib, mg_iphi)%bc_value = 0.0_dp
633 case (bc_cont)
634 ! d/dx u = 0
635 ! mg%bc(iB, mg_iphi)%bc_type = mg_bc_continuous
636 mg%bc(ib, mg_iphi)%bc_type = mg_bc_neumann
637 mg%bc(ib, mg_iphi)%bc_value = 0.0_dp
638 case (bc_periodic)
639 ! Nothing to do here
640 case (bc_noinflow)
641 call usr_special_mg_bc(ib)
642 case (bc_special)
643 call usr_special_mg_bc(ib)
644 case default
645 call mpistop("divE_multigrid warning: unknown b.c. ")
646 end select
647 end do
648 end subroutine rhd_set_mg_bounds
649
650 subroutine rhd_physical_units
652 double precision :: mp,kb
653 double precision :: a,b
654 ! Derive scaling units
655 if(si_unit) then
656 mp=mp_si
657 kb=kb_si
658 else
659 mp=mp_cgs
660 kb=kb_cgs
661 end if
662 if(eq_state_units) then
663 a=1d0+4d0*he_abundance
665 b=1d0+h_ion_fr+he_abundance*(he_ion_fr*(he_ion_fr2+1d0)+1d0)
666 else
667 b=2d0+3d0*he_abundance
668 end if
669 rr=1d0
670 else
671 a=1d0
672 b=1d0
673 rr=(1d0+h_ion_fr+he_abundance*(he_ion_fr*(he_ion_fr2+1d0)+1d0))/(1d0+4d0*he_abundance)
674 end if
675 if(unit_density/=1.d0 .or. unit_numberdensity/=1.d0) then
676 if(unit_density/=1.d0) then
678 else if(unit_numberdensity/=1.d0) then
680 end if
681 if(unit_temperature/=1.d0) then
684 if(unit_length/=1.d0) then
686 else if(unit_time/=1.d0) then
688 end if
689 else if(unit_pressure/=1.d0) then
692 if(unit_length/=1.d0) then
694 else if(unit_time/=1.d0) then
696 end if
697 else if(unit_velocity/=1.d0) then
700 if(unit_length/=1.d0) then
702 else if(unit_time/=1.d0) then
704 end if
705 else if(unit_time/=1.d0) then
709 end if
710 else if(unit_temperature/=1.d0) then
711 ! units of temperature and velocity are dependent
712 if(unit_pressure/=1.d0) then
716 if(unit_length/=1.d0) then
718 else if(unit_time/=1.d0) then
720 end if
721 end if
722 else if(unit_pressure/=1.d0) then
723 if(unit_velocity/=1.d0) then
727 if(unit_length/=1.d0) then
729 else if(unit_time/=1.d0) then
731 end if
732 else if(unit_time/=0.d0) then
737 end if
738 end if
739
740 !> Units for radiative flux and opacity
743
744 end subroutine rhd_physical_units
745
746 !> Returns logical argument flag where values are ok
747 subroutine rhd_check_w(primitive, ixI^L, ixO^L, w, flag)
749 use mod_dust, only: dust_check_w
750
751 logical, intent(in) :: primitive
752 integer, intent(in) :: ixi^l, ixo^l
753 double precision, intent(in) :: w(ixi^s, nw)
754 logical, intent(inout) :: flag(ixi^s,1:nw)
755 double precision :: tmp(ixi^s)
756
757 flag=.false.
758
759 if (rhd_energy) then
760 if (primitive) then
761 where(w(ixo^s, e_) < small_pressure) flag(ixo^s,e_) = .true.
762 else
763 tmp(ixo^s) = (rhd_gamma - 1.0d0)*(w(ixo^s, e_) - &
764 rhd_kin_en(w, ixi^l, ixo^l))
765 where(tmp(ixo^s) < small_pressure) flag(ixo^s,e_) = .true.
766 endif
767 end if
768
769 where(w(ixo^s, rho_) < small_density) flag(ixo^s,rho_) = .true.
770
771 where(w(ixo^s, r_e) < small_r_e) flag(ixo^s,r_e) = .true.
772
773 if(rhd_dust) call dust_check_w(ixi^l,ixo^l,w,flag)
774
775 end subroutine rhd_check_w
776
777 !> Transform primitive variables into conservative ones
778 subroutine rhd_to_conserved(ixI^L, ixO^L, w, x)
780 use mod_dust, only: dust_to_conserved
781 integer, intent(in) :: ixi^l, ixo^l
782 double precision, intent(inout) :: w(ixi^s, nw)
783 double precision, intent(in) :: x(ixi^s, 1:ndim)
784 double precision :: invgam
785 integer :: idir, itr
786
787 !!if (fix_small_values) then
788 !! call rhd_handle_small_values(.true., w, x, ixI^L, ixO^L, 'rhd_to_conserved')
789 !!end if
790
791 if (rhd_energy) then
792 invgam = 1.d0/(rhd_gamma - 1.0d0)
793 ! Calculate total energy from pressure and kinetic energy
794 w(ixo^s, e_) = w(ixo^s, e_) * invgam + &
795 0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) * w(ixo^s, rho_)
796 end if
797
798 ! Convert velocity to momentum
799 do idir = 1, ndir
800 w(ixo^s, mom(idir)) = w(ixo^s, rho_) * w(ixo^s, mom(idir))
801 end do
802
803 if (rhd_dust) then
804 call dust_to_conserved(ixi^l, ixo^l, w, x)
805 end if
806
807 end subroutine rhd_to_conserved
808
809 !> Transform conservative variables into primitive ones
810 subroutine rhd_to_primitive(ixI^L, ixO^L, w, x)
812 use mod_dust, only: dust_to_primitive
813 integer, intent(in) :: ixi^l, ixo^l
814 double precision, intent(inout) :: w(ixi^s, nw)
815 double precision, intent(in) :: x(ixi^s, 1:ndim)
816 integer :: itr, idir
817 double precision :: inv_rho(ixo^s)
818
819 if (fix_small_values) then
820 call rhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'rhd_to_primitive')
821 end if
822
823 inv_rho = 1.0d0 / w(ixo^s, rho_)
824
825 if (rhd_energy) then
826 ! Compute pressure
827 w(ixo^s, e_) = (rhd_gamma - 1.0d0) * (w(ixo^s, e_) - &
828 rhd_kin_en(w, ixi^l, ixo^l, inv_rho))
829 end if
830
831 ! Convert momentum to velocity
832 do idir = 1, ndir
833 w(ixo^s, mom(idir)) = w(ixo^s, mom(idir)) * inv_rho
834 end do
835
836 ! Convert dust momentum to dust velocity
837 if (rhd_dust) then
838 call dust_to_primitive(ixi^l, ixo^l, w, x)
839 end if
840
841 end subroutine rhd_to_primitive
842
843 !> Transform internal energy to total energy
844 subroutine rhd_ei_to_e(ixI^L,ixO^L,w,x)
846 integer, intent(in) :: ixi^l, ixo^l
847 double precision, intent(inout) :: w(ixi^s, nw)
848 double precision, intent(in) :: x(ixi^s, 1:ndim)
849
850 ! Calculate total energy from internal and kinetic energy
851 w(ixo^s,e_)=w(ixo^s,e_)&
852 +rhd_kin_en(w,ixi^l,ixo^l)
853
854 end subroutine rhd_ei_to_e
855
856 !> Transform total energy to internal energy
857 subroutine rhd_e_to_ei(ixI^L,ixO^L,w,x)
859 integer, intent(in) :: ixi^l, ixo^l
860 double precision, intent(inout) :: w(ixi^s, nw)
861 double precision, intent(in) :: x(ixi^s, 1:ndim)
862
863 ! Calculate ei = e - ek
864 w(ixo^s,e_)=w(ixo^s,e_)&
865 -rhd_kin_en(w,ixi^l,ixo^l)
866
867 end subroutine rhd_e_to_ei
868
869 subroutine e_to_rhos(ixI^L, ixO^L, w, x)
871
872 integer, intent(in) :: ixi^l, ixo^l
873 double precision :: w(ixi^s, nw)
874 double precision, intent(in) :: x(ixi^s, 1:ndim)
875
876 if (rhd_energy) then
877 w(ixo^s, e_) = (rhd_gamma - 1.0d0) * w(ixo^s, rho_)**(1.0d0 - rhd_gamma) * &
878 (w(ixo^s, e_) - rhd_kin_en(w, ixi^l, ixo^l))
879 else
880 call mpistop("energy from entropy can not be used with -eos = iso !")
881 end if
882 end subroutine e_to_rhos
883
884 subroutine rhos_to_e(ixI^L, ixO^L, w, x)
886
887 integer, intent(in) :: ixi^l, ixo^l
888 double precision :: w(ixi^s, nw)
889 double precision, intent(in) :: x(ixi^s, 1:ndim)
890
891 if (rhd_energy) then
892 w(ixo^s, e_) = w(ixo^s, rho_)**(rhd_gamma - 1.0d0) * w(ixo^s, e_) &
893 / (rhd_gamma - 1.0d0) + rhd_kin_en(w, ixi^l, ixo^l)
894 else
895 call mpistop("entropy from energy can not be used with -eos = iso !")
896 end if
897 end subroutine rhos_to_e
898
899 !> Calculate v_i = m_i / rho within ixO^L
900 subroutine rhd_get_v(w, x, ixI^L, ixO^L, idim, v)
902 integer, intent(in) :: ixi^l, ixo^l, idim
903 double precision, intent(in) :: w(ixi^s, nw), x(ixi^s, 1:ndim)
904 double precision, intent(out) :: v(ixi^s)
905
906 v(ixo^s) = w(ixo^s, mom(idim)) / w(ixo^s, rho_)
907 end subroutine rhd_get_v
908
909 !> Calculate cmax_idim = csound + abs(v_idim) within ixO^L
910 subroutine rhd_get_cmax(w, x, ixI^L, ixO^L, idim, cmax)
914
915 integer, intent(in) :: ixi^l, ixo^l, idim
916 ! w in primitive form
917 double precision, intent(in) :: w(ixi^s, nw), x(ixi^s, 1:ndim)
918 double precision, intent(inout) :: cmax(ixi^s)
919
920 if(rhd_energy) then
921 cmax(ixo^s)=dabs(w(ixo^s,mom(idim)))+dsqrt(rhd_gamma*w(ixo^s,p_)/w(ixo^s,rho_))
922 else
923 if (.not. associated(usr_set_pthermal)) then
924 select case (rhd_pressure)
925 case ('Trad')
926 cmax(ixo^s) = (w(ixo^s,r_e)*unit_pressure/const_rad_a)**0.25d0&
927 /unit_temperature*w(ixo^s, rho_)
928 case ('adiabatic')
929 cmax(ixo^s) = rhd_adiab * w(ixo^s, rho_)**rhd_gamma
930 case ('Tcond') !> Thermal conduction?!
931 cmax(ixo^s) = (rhd_gamma-1.d0)*w(ixo^s,r_e)
932 case default
933 call mpistop('rhd_pressure unknown, use Trad or adiabatic')
934 end select
935 else
936 call usr_set_pthermal(w,x,ixi^l,ixo^l,cmax)
937 end if
938 cmax(ixo^s)=dabs(w(ixo^s,mom(idim)))+dsqrt(rhd_gamma*cmax(ixo^s)/w(ixo^s,rho_))
939 end if
940
941 if (rhd_dust) then
942 call dust_get_cmax_prim(w, x, ixi^l, ixo^l, idim, cmax)
943 end if
944 end subroutine rhd_get_cmax
945
946 subroutine rhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
948
949 integer, intent(in) :: ixi^l, ixo^l
950 double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
951 double precision, intent(inout) :: a2max(ndim)
952 double precision :: a2(ixi^s,ndim,nw)
953 integer :: gxo^l,hxo^l,jxo^l,kxo^l,i,j
954
955 a2=zero
956 do i = 1,ndim
957 !> 4th order
958 hxo^l=ixo^l-kr(i,^d);
959 gxo^l=hxo^l-kr(i,^d);
960 jxo^l=ixo^l+kr(i,^d);
961 kxo^l=jxo^l+kr(i,^d);
962 a2(ixo^s,i,1:nw)=dabs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
963 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
964 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/dxlevel(i)**2
965 end do
966 end subroutine rhd_get_a2max
967
968 !> get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
969 subroutine rhd_get_tcutoff(ixI^L,ixO^L,w,x,tco_local,Tmax_local)
971 integer, intent(in) :: ixi^l,ixo^l
972 double precision, intent(in) :: x(ixi^s,1:ndim)
973 double precision, intent(inout) :: w(ixi^s,1:nw)
974 double precision, intent(out) :: tco_local, tmax_local
975
976 double precision, parameter :: trac_delta=0.25d0
977 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
978 double precision :: ltr(ixi^s),ltrc,ltrp,tcoff(ixi^s)
979 integer :: jxo^l,hxo^l
980 integer :: jxp^l,hxp^l,ixp^l
981 logical :: lrlt(ixi^s)
982
983 {^ifoned
984 call rhd_get_temperature_from_etot(w,x,ixi^l,ixi^l,te)
985
986 tco_local=zero
987 tmax_local=maxval(te(ixo^s))
988 select case(rhd_trac_type)
989 case(0)
990 w(ixi^s,tcoff_)=3.d5/unit_temperature
991 case(1)
992 hxo^l=ixo^l-1;
993 jxo^l=ixo^l+1;
994 lts(ixo^s)=0.5d0*dabs(te(jxo^s)-te(hxo^s))/te(ixo^s)
995 lrlt=.false.
996 where(lts(ixo^s) > trac_delta)
997 lrlt(ixo^s)=.true.
998 end where
999 if(any(lrlt(ixo^s))) then
1000 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1001 end if
1002 case(2)
1003 !> iijima et al. 2021, LTRAC method
1004 ltrc=1.5d0
1005 ltrp=2.5d0
1006 ixp^l=ixo^l^ladd1;
1007 hxo^l=ixo^l-1;
1008 jxo^l=ixo^l+1;
1009 hxp^l=ixp^l-1;
1010 jxp^l=ixp^l+1;
1011 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
1012 ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1013 w(ixo^s,tcoff_)=te(ixo^s)*&
1014 (0.25*(ltr(jxo^s)+two*ltr(ixo^s)+ltr(hxo^s)))**0.4d0
1015 case default
1016 call mpistop("mrhd_trac_type not allowed for 1D simulation")
1017 end select
1018 }
1019 end subroutine rhd_get_tcutoff
1020
1021 !> Calculate cmax_idim = csound + abs(v_idim) within ixO^L
1022 subroutine rhd_get_cbounds(wLC, wRC, wLp, wRp, x, ixI^L, ixO^L, idim,Hspeed,cmax, cmin)
1024 use mod_dust, only: dust_get_cmax
1025 use mod_variables
1026
1027 integer, intent(in) :: ixi^l, ixo^l, idim
1028 ! conservative left and right status
1029 double precision, intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
1030 ! primitive left and right status
1031 double precision, intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
1032 double precision, intent(in) :: x(ixi^s, 1:ndim)
1033 double precision, intent(inout) :: cmax(ixi^s,1:number_species)
1034 double precision, intent(inout), optional :: cmin(ixi^s,1:number_species)
1035 double precision, intent(in) :: hspeed(ixi^s,1:number_species)
1036
1037 double precision :: wmean(ixi^s,nw)
1038 double precision, dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
1039 integer :: ix^d
1040
1041 select case(boundspeed)
1042 case (1)
1043 ! This implements formula (10.52) from "Riemann Solvers and Numerical
1044 ! Methods for Fluid Dynamics" by Toro.
1045
1046 tmp1(ixo^s)=dsqrt(wlp(ixo^s,rho_))
1047 tmp2(ixo^s)=dsqrt(wrp(ixo^s,rho_))
1048 tmp3(ixo^s)=1.d0/(dsqrt(wlp(ixo^s,rho_))+dsqrt(wrp(ixo^s,rho_)))
1049 umean(ixo^s)=(wlp(ixo^s,mom(idim))*tmp1(ixo^s)+wrp(ixo^s,mom(idim))*tmp2(ixo^s))*tmp3(ixo^s)
1050
1051 if(rhd_energy) then
1052 csoundl(ixo^s)=rhd_gamma*wlp(ixo^s,p_)/wlp(ixo^s,rho_)
1053 csoundr(ixo^s)=rhd_gamma*wrp(ixo^s,p_)/wrp(ixo^s,rho_)
1054 else
1055 select case (rhd_pressure)
1056 case ('Trad')
1057 csoundl(ixo^s)=rhd_gamma*kbmpmua4*wlp(ixo^s,r_e)**(1.d0/4)
1058 csoundr(ixo^s)=rhd_gamma*kbmpmua4*wrp(ixo^s,r_e)**(1.d0/4)
1059 case ('adiabatic')
1060 csoundl(ixo^s)=rhd_gamma*rhd_adiab*wlp(ixo^s,rho_)**(rhd_gamma-one)
1061 csoundr(ixo^s)=rhd_gamma*rhd_adiab*wrp(ixo^s,rho_)**(rhd_gamma-one)
1062 end select
1063 end if
1064
1065 dmean(ixo^s) = (tmp1(ixo^s)*csoundl(ixo^s)+tmp2(ixo^s)*csoundr(ixo^s)) * &
1066 tmp3(ixo^s) + 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2 * &
1067 (wrp(ixo^s,mom(idim))-wlp(ixo^s,mom(idim)))**2
1068
1069 dmean(ixo^s)=dsqrt(dmean(ixo^s))
1070 if(present(cmin)) then
1071 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
1072 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
1073 if(h_correction) then
1074 {do ix^db=ixomin^db,ixomax^db\}
1075 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
1076 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
1077 {end do\}
1078 end if
1079 else
1080 cmax(ixo^s,1)=dabs(umean(ixo^s))+dmean(ixo^s)
1081 end if
1082
1083 if (rhd_dust) then
1084 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
1085 call dust_get_cmax(wmean, x, ixi^l, ixo^l, idim, cmax, cmin)
1086 end if
1087
1088 case (2)
1089 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
1090 tmp1(ixo^s)=wmean(ixo^s,mom(idim))/wmean(ixo^s,rho_)
1091 call rhd_get_csound2(wmean,x,ixi^l,ixo^l,csoundr)
1092 csoundr(ixo^s) = dsqrt(csoundr(ixo^s))
1093
1094 if(present(cmin)) then
1095 cmax(ixo^s,1)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
1096 cmin(ixo^s,1)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
1097 if(h_correction) then
1098 {do ix^db=ixomin^db,ixomax^db\}
1099 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
1100 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
1101 {end do\}
1102 end if
1103 else
1104 cmax(ixo^s,1)=dabs(tmp1(ixo^s))+csoundr(ixo^s)
1105 end if
1106
1107 if (rhd_dust) then
1108 call dust_get_cmax(wmean, x, ixi^l, ixo^l, idim, cmax, cmin)
1109 end if
1110 case (3)
1111 ! Miyoshi 2005 JCP 208, 315 equation (67)
1112 if(rhd_energy) then
1113 csoundl(ixo^s)=rhd_gamma*wlp(ixo^s,p_)/wlp(ixo^s,rho_)
1114 csoundr(ixo^s)=rhd_gamma*wrp(ixo^s,p_)/wrp(ixo^s,rho_)
1115 else
1116 select case (rhd_pressure)
1117 case ('Trad')
1118 csoundl(ixo^s)=rhd_gamma*kbmpmua4*wlp(ixo^s,r_e)**(1.d0/4)
1119 csoundr(ixo^s)=rhd_gamma*kbmpmua4*wrp(ixo^s,r_e)**(1.d0/4)
1120 case ('adiabatic')
1121 csoundl(ixo^s)=rhd_gamma*rhd_adiab*wlp(ixo^s,rho_)**(rhd_gamma-one)
1122 csoundr(ixo^s)=rhd_gamma*rhd_adiab*wrp(ixo^s,rho_)**(rhd_gamma-one)
1123 end select
1124 end if
1125 csoundl(ixo^s)=max(dsqrt(csoundl(ixo^s)),dsqrt(csoundr(ixo^s)))
1126 if(present(cmin)) then
1127 cmin(ixo^s,1)=min(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))-csoundl(ixo^s)
1128 cmax(ixo^s,1)=max(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
1129 if(h_correction) then
1130 {do ix^db=ixomin^db,ixomax^db\}
1131 cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
1132 cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
1133 {end do\}
1134 end if
1135 else
1136 cmax(ixo^s,1)=max(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
1137 end if
1138 if (rhd_dust) then
1139 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
1140 call dust_get_cmax(wmean, x, ixi^l, ixo^l, idim, cmax, cmin)
1141 end if
1142 end select
1143
1144 end subroutine rhd_get_cbounds
1145
1146 !> Calculate the square of the thermal sound speed csound2 within ixO^L.
1147 !> csound2=gamma*p/rho
1148 subroutine rhd_get_csound2(w,x,ixI^L,ixO^L,csound2)
1150 integer, intent(in) :: ixi^l, ixo^l
1151 double precision, intent(in) :: w(ixi^s,nw)
1152 double precision, intent(in) :: x(ixi^s,1:ndim)
1153 double precision, intent(out) :: csound2(ixi^s)
1154
1155 call rhd_get_ptot(w,x,ixi^l,ixo^l,csound2)
1156 csound2(ixo^s)=max(rhd_gamma,4.d0/3.d0)*csound2(ixo^s)/w(ixo^s,rho_)
1157 !> Turner & Stone 2001
1158
1159 end subroutine rhd_get_csound2
1160
1161 !> Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho) within ixO^L
1162 subroutine rhd_get_pthermal(w, x, ixI^L, ixO^L, pth)
1166
1167 integer, intent(in) :: ixi^l, ixo^l
1168 double precision, intent(in) :: w(ixi^s, 1:nw)
1169 double precision, intent(in) :: x(ixi^s, 1:ndim)
1170 double precision, intent(out):: pth(ixi^s)
1171 integer :: iw, ix^d
1172
1173 if (rhd_energy) then
1174 pth(ixi^s) = (rhd_gamma - 1.d0) * (w(ixi^s, e_) - &
1175 0.5d0 * sum(w(ixi^s, mom(:))**2, dim=ndim+1) / w(ixi^s, rho_))
1176 else
1177 if (.not. associated(usr_set_pthermal)) then
1178 select case (rhd_pressure)
1179 case ('Trad')
1180 pth(ixi^s) = (w(ixi^s,r_e)*unit_pressure/const_rad_a)**0.25d0&
1181 /unit_temperature*w(ixi^s, rho_)
1182 case ('adiabatic')
1183 pth(ixi^s) = rhd_adiab * w(ixi^s, rho_)**rhd_gamma
1184 case ('Tcond') !> Thermal conduction?!
1185 pth(ixi^s) = (rhd_gamma-1.d0)*w(ixi^s,r_e)
1186 case default
1187 call mpistop('rhd_pressure unknown, use Trad or adiabatic')
1188 end select
1189 else
1190 call usr_set_pthermal(w,x,ixi^l,ixo^l,pth)
1191 end if
1192 end if
1193
1194 if (fix_small_values) then
1195 {do ix^db= ixo^lim^db\}
1196 if(pth(ix^d)<small_pressure) then
1197 pth(ix^d)=small_pressure
1198 endif
1199 {enddo^d&\}
1200 else if (check_small_values) then
1201 {do ix^db= ixo^lim^db\}
1202 if(pth(ix^d)<small_pressure) then
1203 write(*,*) "Error: small value of gas pressure",pth(ix^d),&
1204 " encountered when call rhd_get_pthermal"
1205 write(*,*) "Iteration: ", it, " Time: ", global_time
1206 write(*,*) "Location: ", x(ix^d,:)
1207 write(*,*) "Cell number: ", ix^d
1208 do iw=1,nw
1209 write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
1210 end do
1211 ! use erroneous arithmetic operation to crash the run
1212 if(trace_small_values) write(*,*) dsqrt(pth(ix^d)-bigdouble)
1213 write(*,*) "Saving status at the previous time step"
1214 crash=.true.
1215 end if
1216 {enddo^d&\}
1217 end if
1218
1219 end subroutine rhd_get_pthermal
1220
1221 !> Calculate radiation pressure within ixO^L
1222 subroutine rhd_get_pradiation(w, x, ixI^L, ixO^L, prad)
1224 use mod_fld
1225 use mod_afld
1226
1227 integer, intent(in) :: ixi^l, ixo^l
1228 double precision, intent(in) :: w(ixi^s, 1:nw)
1229 double precision, intent(in) :: x(ixi^s, 1:ndim)
1230 double precision, intent(out):: prad(ixo^s, 1:ndim, 1:ndim)
1231
1232 select case (rhd_radiation_formalism)
1233 case('fld')
1234 call fld_get_radpress(w, x, ixi^l, ixo^l, prad, nghostcells)
1235 case('afld')
1236 call afld_get_radpress(w, x, ixi^l, ixo^l, prad, nghostcells)
1237 case default
1238 call mpistop('Radiation formalism unknown')
1239 end select
1240 end subroutine rhd_get_pradiation
1241
1242 !> calculates the sum of the gas pressure and max Prad tensor element
1243 subroutine rhd_get_ptot(w, x, ixI^L, ixO^L, ptot)
1245
1246 integer, intent(in) :: ixi^l, ixo^l
1247 double precision, intent(in) :: w(ixi^s, 1:nw)
1248 double precision, intent(in) :: x(ixi^s, 1:ndim)
1249 double precision :: pth(ixi^s)
1250 double precision :: prad_tensor(ixo^s, 1:ndim, 1:ndim)
1251 double precision :: prad_max(ixo^s)
1252 double precision, intent(out):: ptot(ixi^s)
1253 integer :: ix^d
1254
1255 call rhd_get_pthermal(w, x, ixi^l, ixo^l, pth)
1256 call rhd_get_pradiation(w, x, ixi^l, ixo^l, prad_tensor)
1257
1258 {do ix^d = ixomin^d,ixomax^d\}
1259 prad_max(ix^d) = maxval(prad_tensor(ix^d,:,:))
1260 {enddo\}
1261
1262 !> filter cmax
1263 if (radio_acoustic_filter) then
1264 call rhd_radio_acoustic_filter(x, ixi^l, ixo^l, prad_max)
1265 endif
1266
1267 ptot(ixo^s) = pth(ixo^s) + prad_max(ixo^s)
1268
1269 end subroutine rhd_get_ptot
1270
1271 !> Filter peaks in cmax due to radiation energy density, used for debugging
1272 subroutine rhd_radio_acoustic_filter(x, ixI^L, ixO^L, prad_max)
1274
1275 integer, intent(in) :: ixi^l, ixo^l
1276 double precision, intent(in) :: x(ixi^s, 1:ndim)
1277 double precision, intent(inout) :: prad_max(ixo^s)
1278
1279 double precision :: tmp_prad(ixi^s)
1280 integer :: ix^d, filter, idim
1281
1282 if (size_ra_filter .lt. 1) call mpistop("ra filter of size < 1 makes no sense")
1283 if (size_ra_filter .gt. nghostcells) call mpistop("ra filter of size < nghostcells makes no sense")
1284
1285 tmp_prad(ixi^s) = zero
1286 tmp_prad(ixo^s) = prad_max(ixo^s)
1287
1288 do filter = 1,size_ra_filter
1289 do idim = 1,ndim
1290 ! {do ix^D = ixOmin^D+filter,ixOmax^D-filter\}
1291 {do ix^d = ixomin^d,ixomax^d\}
1292 prad_max(ix^d) = min(tmp_prad(ix^d),tmp_prad(ix^d+filter*kr(idim,^d)))
1293 prad_max(ix^d) = min(tmp_prad(ix^d),tmp_prad(ix^d-filter*kr(idim,^d)))
1294 {enddo\}
1295 enddo
1296 enddo
1297 end subroutine rhd_radio_acoustic_filter
1298
1299 !> Calculate temperature=p/rho when in e_ the total energy is stored
1300 subroutine rhd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
1302 integer, intent(in) :: ixi^l, ixo^l
1303 double precision, intent(in) :: w(ixi^s, 1:nw)
1304 double precision, intent(in) :: x(ixi^s, 1:ndim)
1305 double precision, intent(out):: res(ixi^s)
1306
1307 double precision :: r(ixi^s)
1308
1309 call rhd_get_rfactor(w,x,ixi^l,ixo^l,r)
1310 call rhd_get_pthermal(w, x, ixi^l, ixo^l, res)
1311 res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,rho_))
1312 end subroutine rhd_get_temperature_from_etot
1313
1314
1315 !> Calculate temperature=p/rho when in e_ the internal energy is stored
1316 subroutine rhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
1318 integer, intent(in) :: ixi^l, ixo^l
1319 double precision, intent(in) :: w(ixi^s, 1:nw)
1320 double precision, intent(in) :: x(ixi^s, 1:ndim)
1321 double precision, intent(out):: res(ixi^s)
1322 double precision :: r(ixi^s)
1323
1324 call rhd_get_rfactor(w,x,ixi^l,ixo^l,r)
1325 res(ixo^s) = (rhd_gamma - 1.0d0) * w(ixo^s, e_)/(w(ixo^s,rho_)*r(ixo^s))
1326 end subroutine rhd_get_temperature_from_eint
1327
1328 !> Calculates gas temperature
1329 subroutine rhd_get_tgas(w, x, ixI^L, ixO^L, tgas)
1331
1332 integer, intent(in) :: ixi^l, ixo^l
1333 double precision, intent(in) :: w(ixi^s, 1:nw)
1334 double precision, intent(in) :: x(ixi^s, 1:ndim)
1335 double precision :: pth(ixi^s)
1336 double precision, intent(out):: tgas(ixi^s)
1337
1338 call rhd_get_pthermal(w, x, ixi^l, ixo^l, pth)
1339 tgas(ixi^s) = pth(ixi^s)/w(ixi^s,rho_)
1340
1341 end subroutine rhd_get_tgas
1342
1343 !> Calculates radiation temperature
1344 subroutine rhd_get_trad(w, x, ixI^L, ixO^L, trad)
1346 use mod_constants
1347
1348 integer, intent(in) :: ixi^l, ixo^l
1349 double precision, intent(in) :: w(ixi^s, 1:nw)
1350 double precision, intent(in) :: x(ixi^s, 1:ndim)
1351 double precision, intent(out):: trad(ixi^s)
1352
1353 trad(ixi^s) = (w(ixi^s,r_e)*unit_pressure&
1354 /const_rad_a)**(1.d0/4.d0)/unit_temperature
1355
1356 end subroutine rhd_get_trad
1357
1358
1359 !these are very similar to the subroutines without 1, used in mod_thermal_conductivity
1360 !but no check on whether energy variable is present
1361 subroutine rhd_ei_to_e1(ixI^L,ixO^L,w,x)
1363 integer, intent(in) :: ixi^l, ixo^l
1364 double precision, intent(inout) :: w(ixi^s, nw)
1365 double precision, intent(in) :: x(ixi^s, 1:ndim)
1366
1367 ! Calculate total energy from internal and kinetic energy
1368 w(ixo^s,e_)=w(ixo^s,e_)&
1369 +rhd_kin_en(w,ixi^l,ixo^l)
1370
1371 end subroutine rhd_ei_to_e1
1372
1373 !> Transform total energy to internal energy
1374 !but no check on whether energy variable is present
1375 subroutine rhd_e_to_ei1(ixI^L,ixO^L,w,x)
1377 integer, intent(in) :: ixi^l, ixo^l
1378 double precision, intent(inout) :: w(ixi^s, nw)
1379 double precision, intent(in) :: x(ixi^s, 1:ndim)
1380
1381 ! Calculate ei = e - ek
1382 w(ixo^s,e_)=w(ixo^s,e_)&
1383 -rhd_kin_en(w,ixi^l,ixo^l)
1384
1385 end subroutine rhd_e_to_ei1
1386
1387 ! Calculate flux f_idim[iw]
1388 subroutine rhd_get_flux_cons(w, x, ixI^L, ixO^L, idim, f)
1390 use mod_dust, only: dust_get_flux
1391
1392 integer, intent(in) :: ixi^l, ixo^l, idim
1393 double precision, intent(in) :: w(ixi^s, 1:nw), x(ixi^s, 1:ndim)
1394 double precision, intent(out) :: f(ixi^s, nwflux)
1395 double precision :: pth(ixi^s), v(ixi^s),frame_vel(ixi^s)
1396 integer :: idir, itr
1397
1398 call rhd_get_pthermal(w, x, ixi^l, ixo^l, pth)
1399 call rhd_get_v(w, x, ixi^l, ixo^l, idim, v)
1400
1401 f(ixo^s, rho_) = v(ixo^s) * w(ixo^s, rho_)
1402
1403 ! Momentum flux is v_i*m_i, +p in direction idim
1404 do idir = 1, ndir
1405 f(ixo^s, mom(idir)) = v(ixo^s) * w(ixo^s, mom(idir))
1406 end do
1407
1408 f(ixo^s, mom(idim)) = f(ixo^s, mom(idim)) + pth(ixo^s)
1409
1410 if(rhd_energy) then
1411 ! Energy flux is v_i*(e + p)
1412 f(ixo^s, e_) = v(ixo^s) * (w(ixo^s, e_) + pth(ixo^s))
1413 end if
1414
1415 if (rhd_radiation_advection) then
1416 f(ixo^s, r_e) = v(ixo^s) * w(ixo^s,r_e)
1417 else
1418 f(ixo^s, r_e) = zero
1419 endif
1420
1421 do itr = 1, rhd_n_tracer
1422 f(ixo^s, tracer(itr)) = v(ixo^s) * w(ixo^s, tracer(itr))
1423 end do
1424
1425 ! Dust fluxes
1426 if (rhd_dust) then
1427 call dust_get_flux(w, x, ixi^l, ixo^l, idim, f)
1428 end if
1429
1430 end subroutine rhd_get_flux_cons
1431
1432 ! Calculate flux f_idim[iw]
1433 subroutine rhd_get_flux(wC, w, x, ixI^L, ixO^L, idim, f)
1435 use mod_dust, only: dust_get_flux_prim
1436 use mod_viscosity, only: visc_get_flux_prim ! viscInDiv
1437
1438 integer, intent(in) :: ixi^l, ixo^l, idim
1439 ! conservative w
1440 double precision, intent(in) :: wc(ixi^s, 1:nw)
1441 ! primitive w
1442 double precision, intent(in) :: w(ixi^s, 1:nw)
1443 double precision, intent(in) :: x(ixi^s, 1:ndim)
1444 double precision, intent(out) :: f(ixi^s, nwflux)
1445 double precision :: pth(ixi^s),frame_vel(ixi^s)
1446 integer :: idir, itr
1447
1448 if (rhd_energy) then
1449 pth(ixo^s) = w(ixo^s,p_)
1450 else
1451 call rhd_get_pthermal(wc, x, ixi^l, ixo^l, pth)
1452 end if
1453
1454 f(ixo^s, rho_) = w(ixo^s,mom(idim)) * w(ixo^s, rho_)
1455
1456 ! Momentum flux is v_i*m_i, +p in direction idim
1457 do idir = 1, ndir
1458 f(ixo^s, mom(idir)) = w(ixo^s,mom(idim)) * wc(ixo^s, mom(idir))
1459 end do
1460
1461 f(ixo^s, mom(idim)) = f(ixo^s, mom(idim)) + pth(ixo^s)
1462
1463 if(rhd_energy) then
1464 ! Energy flux is v_i*(e + p)
1465 f(ixo^s, e_) = w(ixo^s,mom(idim)) * (wc(ixo^s, e_) + w(ixo^s,p_))
1466 end if
1467
1468 if (rhd_radiation_advection) then
1469 f(ixo^s, r_e) = w(ixo^s,mom(idim)) * wc(ixo^s,r_e)
1470 else
1471 f(ixo^s, r_e) = zero
1472 endif
1473
1474 do itr = 1, rhd_n_tracer
1475 f(ixo^s, tracer(itr)) = w(ixo^s,mom(idim)) * w(ixo^s, tracer(itr))
1476 end do
1477
1478 ! Dust fluxes
1479 if (rhd_dust) then
1480 call dust_get_flux_prim(w, x, ixi^l, ixo^l, idim, f)
1481 end if
1482
1483 ! Viscosity fluxes - viscInDiv
1484 if (rhd_viscosity) then
1485 call visc_get_flux_prim(w, x, ixi^l, ixo^l, idim, f, rhd_energy)
1486 endif
1487
1488 end subroutine rhd_get_flux
1489
1490 !> Add geometrical source terms to w
1491 !>
1492 !> Notice that the expressions of the geometrical terms depend only on ndir,
1493 !> not ndim. Eg, they are the same in 2.5D and in 3D, for any geometry.
1494 !>
1495 subroutine rhd_add_source_geom(qdt, dtfactor, ixI^L, ixO^L, wCT, wprim, w, x)
1498 use mod_viscosity, only: visc_add_source_geom ! viscInDiv
1501 use mod_geometry
1502 integer, intent(in) :: ixi^l, ixo^l
1503 double precision, intent(in) :: qdt, dtfactor, x(ixi^s, 1:ndim)
1504 double precision, intent(inout) :: wct(ixi^s, 1:nw), wprim(ixi^s,1:nw), w(ixi^s, 1:nw)
1505 ! to change and to set as a parameter in the parfile once the possibility to
1506 ! solve the equations in an angular momentum conserving form has been
1507 ! implemented (change tvdlf.t eg)
1508 double precision :: pth(ixi^s), source(ixi^s), minrho
1509 integer :: iw,idir, h1x^l{^nooned, h2x^l}
1510 integer :: mr_,mphi_ ! Polar var. names
1511 integer :: irho, ifluid, n_fluids
1512 double precision :: exp_factor(ixi^s), del_exp_factor(ixi^s), exp_factor_primitive(ixi^s)
1513
1514 if (rhd_dust) then
1515 n_fluids = 1 + dust_n_species
1516 else
1517 n_fluids = 1
1518 end if
1519
1520 select case (coordinate)
1521
1523 !the user provides the functions of exp_factor and del_exp_factor
1524 if(associated(usr_set_surface)) call usr_set_surface(ixi^l,x,block%dx,exp_factor,del_exp_factor,exp_factor_primitive)
1525 if(rhd_energy) then
1526 source(ixo^s) =wprim(ixo^s, p_)
1527 else
1528 if(.not. associated(usr_set_pthermal)) then
1529 source(ixo^s)=rhd_adiab * wprim(ixo^s, rho_)**rhd_gamma
1530 else
1531 call usr_set_pthermal(wct,x,ixi^l,ixo^l,source)
1532 end if
1533 end if
1534 source(ixo^s) = source(ixo^s)*del_exp_factor(ixo^s)/exp_factor(ixo^s)
1535 w(ixo^s,mom(1)) = w(ixo^s,mom(1)) + qdt*source(ixo^s)
1536
1537 case (cylindrical)
1538 if ((rhd_radiation_formalism .eq. 'fld') .or. (rhd_radiation_formalism .eq. 'afld')) then
1539 call mpistop("Diffusion term not implemented yet with cylkindrical geometries")
1540 end if
1541
1542 do ifluid = 0, n_fluids-1
1543 ! s[mr]=(pthermal+mphi**2/rho)/radius
1544 if (ifluid == 0) then
1545 ! gas
1546 irho = rho_
1547 mr_ = mom(r_)
1548 if(phi_>0) mphi_ = mom(phi_)
1549 if(rhd_energy) then
1550 source(ixo^s) =wprim(ixo^s, p_)
1551 else
1552 if(.not. associated(usr_set_pthermal)) then
1553 source(ixo^s)=rhd_adiab * wprim(ixo^s, rho_)**rhd_gamma
1554 else
1555 call usr_set_pthermal(wct,x,ixi^l,ixo^l,source)
1556 end if
1557 end if
1558 minrho = 0.0d0
1559 else
1560 ! dust : no pressure
1561 irho = dust_rho(ifluid)
1562 mr_ = dust_mom(r_, ifluid)
1563 if(phi_>0) mphi_ = dust_mom(phi_, ifluid)
1564 source(ixi^s) = zero
1565 minrho = 0.0d0
1566 end if
1567 if (phi_ > 0) then
1568 source(ixo^s) = source(ixo^s) + wprim(ixo^s, mphi_)**2 * wprim(ixo^s, irho)
1569 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * source(ixo^s) / x(ixo^s, r_)
1570 ! s[mphi]=(-vphi*vr*rho)/radius
1571 source(ixo^s) = -wprim(ixo^s, mphi_) * wprim(ixo^s, mr_) * wprim(ixo^s, irho)
1572 w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt * source(ixo^s) / x(ixo^s, r_)
1573 else
1574 ! s[mr]=2pthermal/radius
1575 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * source(ixo^s) / x(ixo^s, r_)
1576 end if
1577 end do
1578 case (spherical)
1579
1580 if ((rhd_radiation_formalism .eq. 'fld') .or. (rhd_radiation_formalism .eq. 'afld')) then
1581 call mpistop("Diffusion term not implemented yet with spherical geometries")
1582 end if
1583
1584 if (rhd_dust) then
1585 call mpistop("Dust geom source terms not implemented yet with spherical geometries")
1586 end if
1587 mr_ = mom(r_)
1588 if(phi_>0) mphi_ = mom(phi_)
1589 h1x^l=ixo^l-kr(1,^d); {^nooned h2x^l=ixo^l-kr(2,^d);}
1590 ! s[mr]=((vtheta**2+vphi**2)*rho+2*p)/r
1591 if(rhd_energy) then
1592 pth(ixo^s) =wprim(ixo^s, p_)
1593 else
1594 if(.not. associated(usr_set_pthermal)) then
1595 pth(ixo^s)=rhd_adiab * wprim(ixo^s, rho_)**rhd_gamma
1596 else
1597 call usr_set_pthermal(wct,x,ixi^l,ixo^l,pth)
1598 end if
1599 end if
1600 source(ixo^s) = pth(ixo^s) * x(ixo^s, 1) &
1601 *(block%surfaceC(ixo^s, 1) - block%surfaceC(h1x^s, 1)) &
1602 /block%dvolume(ixo^s)
1603 do idir = 2, ndir
1604 source(ixo^s) = source(ixo^s) + wprim(ixo^s, mom(idir))**2 * wprim(ixo^s, rho_)
1605 end do
1606 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * source(ixo^s) / x(ixo^s, 1)
1607
1608 {^nooned
1609 ! s[mtheta]=-(vr*vtheta*rho)/r+cot(theta)*(vphi**2*rho+p)/r
1610 source(ixo^s) = pth(ixo^s) * x(ixo^s, 1) &
1611 * (block%surfaceC(ixo^s, 2) - block%surfaceC(h2x^s, 2)) &
1612 / block%dvolume(ixo^s)
1613 if (ndir == 3) then
1614 source(ixo^s) = source(ixo^s) + (wprim(ixo^s, mom(3))**2 * wprim(ixo^s, rho_)) / tan(x(ixo^s, 2))
1615 end if
1616 source(ixo^s) = source(ixo^s) - (wprim(ixo^s, mom(2)) * wprim(ixo^s, mr_)) * wprim(ixo^s, rho_)
1617 w(ixo^s, mom(2)) = w(ixo^s, mom(2)) + qdt * source(ixo^s) / x(ixo^s, 1)
1618
1619 if (ndir == 3) then
1620 ! s[mphi]=-(vphi*vr/rho)/r-cot(theta)*(vtheta*vphi/rho)/r
1621 source(ixo^s) = -(wprim(ixo^s, mom(3)) * wprim(ixo^s, mr_)) * wprim(ixo^s, rho_)&
1622 - (wprim(ixo^s, mom(2)) * wprim(ixo^s, mom(3))) * wprim(ixo^s, rho_) / tan(x(ixo^s, 2))
1623 w(ixo^s, mom(3)) = w(ixo^s, mom(3)) + qdt * source(ixo^s) / x(ixo^s, 1)
1624 end if
1625 }
1626 end select
1627
1628 if (rhd_viscosity) call visc_add_source_geom(qdt,ixi^l,ixo^l,wprim,w,x)
1629
1630 if (rhd_rotating_frame) then
1631 if (rhd_dust) then
1632 call mpistop("Rotating frame not implemented yet with dust")
1633 else
1634 call rotating_frame_add_source(qdt,dtfactor,ixi^l,ixo^l,wprim,w,x)
1635 end if
1636 end if
1637
1638 end subroutine rhd_add_source_geom
1639
1640 ! w[iw]= w[iw]+qdt*S[wCT, qtC, x] where S is the source based on wCT within ixO
1641 subroutine rhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
1646 use mod_usr_methods, only: usr_gravity
1648
1649 integer, intent(in) :: ixi^l, ixo^l
1650 double precision, intent(in) :: qdt,dtfactor
1651 double precision, intent(in) :: wct(ixi^s, 1:nw),wctprim(ixi^s,1:nw),x(ixi^s, 1:ndim)
1652 double precision, intent(inout) :: w(ixi^s, 1:nw)
1653 logical, intent(in) :: qsourcesplit
1654 logical, intent(inout) :: active
1655
1656 double precision :: gravity_field(ixi^s, 1:ndim)
1657 integer :: idust, idim
1658
1659 if(rhd_dust) then
1660 call dust_add_source(qdt,ixi^l,ixo^l,wct,w,x,qsourcesplit,active)
1661 end if
1662
1663 if(rhd_radiative_cooling) then
1664 call radiative_cooling_add_source(qdt,ixi^l,ixo^l,wct,wctprim,w,x,&
1665 qsourcesplit,active, rc_fl)
1666 end if
1667
1668 if(rhd_viscosity) then
1669 call viscosity_add_source(qdt,ixi^l,ixo^l,wct,w,x,&
1670 rhd_energy,qsourcesplit,active)
1671 end if
1672
1673 if (rhd_gravity) then
1674 call gravity_add_source(qdt,ixi^l,ixo^l,wct,wctprim,w,x,&
1675 rhd_energy,.false.,qsourcesplit,active)
1676
1677 if (rhd_dust .and. qsourcesplit .eqv. grav_split) then
1678 active = .true.
1679
1680 call usr_gravity(ixi^l, ixo^l, wct, x, gravity_field)
1681 do idust = 1, dust_n_species
1682 do idim = 1, ndim
1683 w(ixo^s, dust_mom(idim, idust)) = w(ixo^s, dust_mom(idim, idust)) &
1684 + qdt * gravity_field(ixo^s, idim) * wct(ixo^s, dust_rho(idust))
1685 end do
1686 end do
1687 end if
1688 end if
1689
1690 !> This is where the radiation force and heating/cooling are added/
1691 call rhd_add_radiation_source(qdt,ixi^l,ixo^l,wct,w,x,qsourcesplit,active)
1692
1693 if(rhd_partial_ionization) then
1694 if(.not.qsourcesplit) then
1695 active = .true.
1696 call rhd_update_temperature(ixi^l,ixo^l,wct,w,x)
1697 end if
1698 end if
1699
1700 end subroutine rhd_add_source
1701
1702 subroutine rhd_add_radiation_source(qdt,ixI^L,ixO^L,wCT,w,x,qsourcesplit,active)
1703 use mod_constants
1705 use mod_usr_methods
1706 use mod_fld!, only: fld_get_diffcoef_central, get_fld_rad_force, get_fld_energy_interact
1707 use mod_afld!, only: afld_get_diffcoef_central, get_afld_rad_force, get_afld_energy_interact
1708
1709 integer, intent(in) :: ixi^l, ixo^l
1710 double precision, intent(in) :: qdt, x(ixi^s,1:ndim)
1711 double precision, intent(in) :: wct(ixi^s,1:nw)
1712 double precision, intent(inout) :: w(ixi^s,1:nw)
1713 logical, intent(in) :: qsourcesplit
1714 logical, intent(inout) :: active
1715 double precision :: cmax(ixi^s)
1716
1717 select case(rhd_radiation_formalism)
1718 case('fld')
1719 if(fld_diff_scheme .eq. 'mg') call fld_get_diffcoef_central(w, wct, x, ixi^l, ixo^l)
1720 !> radiation force
1721 if(rhd_radiation_force) call get_fld_rad_force(qdt,ixi^l,ixo^l,wct,w,x,&
1722 rhd_energy,qsourcesplit,active)
1723 call rhd_handle_small_values(.true., w, x, ixi^l, ixo^l, 'fld_e_interact')
1724 case('afld')
1725 if (fld_diff_scheme .eq. 'mg') call afld_get_diffcoef_central(w, wct, x, ixi^l, ixo^l)
1726 !> radiation force
1727 if (rhd_radiation_force) call get_afld_rad_force(qdt,ixi^l,ixo^l,wct,w,x,&
1728 rhd_energy,qsourcesplit,active)
1729 call rhd_handle_small_values(.true., w, x, ixi^l, ixo^l, 'fld_e_interact')
1730 !> photon tiring, heating and cooling
1731 if (rhd_energy) then
1732 if (rhd_energy_interact) call get_afld_energy_interact(qdt,ixi^l,ixo^l,wct,w,x,&
1733 rhd_energy,qsourcesplit,active)
1734 endif
1735 case default
1736 call mpistop('Radiation formalism unknown')
1737 end select
1738
1739 ! ! !> NOT necessary for calculation, just want to know the grid-dependent-timestep
1740 ! call rhd_get_cmax(w, x, ixI^L, ixO^L, 2, cmax)
1741 ! w(ixI^S,i_test) = cmax(ixI^S)
1742 end subroutine rhd_add_radiation_source
1743
1744 subroutine rhd_get_dt(w, ixI^L, ixO^L, dtnew, dx^D, x)
1746 use mod_dust, only: dust_get_dt
1749 use mod_gravity, only: gravity_get_dt
1750 use mod_fld, only: fld_radforce_get_dt
1752
1753 integer, intent(in) :: ixi^l, ixo^l
1754 double precision, intent(in) :: dx^d, x(ixi^s, 1:^nd)
1755 double precision, intent(in) :: w(ixi^s, 1:nw)
1756 double precision, intent(inout) :: dtnew
1757
1758 dtnew = bigdouble
1759
1760 if (.not. dt_c) then
1761
1762 if(rhd_dust) then
1763 call dust_get_dt(w, ixi^l, ixo^l, dtnew, dx^d, x)
1764 end if
1765
1766 if(rhd_radiation_force) then
1767 select case(rhd_radiation_formalism)
1768 case('fld')
1769 call fld_radforce_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
1770 case('afld')
1771 call afld_radforce_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
1772 case default
1773 call mpistop('Radiation formalism unknown')
1774 end select
1775 endif
1776
1777 if(rhd_radiative_cooling) then
1778 call cooling_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x,rc_fl)
1779 end if
1780
1781 if(rhd_viscosity) then
1782 call viscosity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
1783 end if
1784
1785 if(rhd_gravity) then
1786 call gravity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
1787 end if
1788 else
1789 {^ifoned dtnew = dx1*unit_velocity/const_c}
1790 {^nooned dtnew = min(dx^d*unit_velocity/const_c)}
1791 endif
1792
1793 end subroutine rhd_get_dt
1794
1795 function rhd_kin_en(w, ixI^L, ixO^L, inv_rho) result(ke)
1796 use mod_global_parameters, only: nw, ndim
1797 integer, intent(in) :: ixi^l, ixo^l
1798 double precision, intent(in) :: w(ixi^s, nw)
1799 double precision :: ke(ixo^s)
1800 double precision, intent(in), optional :: inv_rho(ixo^s)
1801
1802 if (present(inv_rho)) then
1803 ke = 0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) * inv_rho
1804 else
1805 ke = 0.5d0 * sum(w(ixo^s, mom(:))**2, dim=ndim+1) / w(ixo^s, rho_)
1806 end if
1807 end function rhd_kin_en
1808
1809 function rhd_inv_rho(w, ixI^L, ixO^L) result(inv_rho)
1810 use mod_global_parameters, only: nw, ndim
1811 integer, intent(in) :: ixi^l, ixo^l
1812 double precision, intent(in) :: w(ixi^s, nw)
1813 double precision :: inv_rho(ixo^s)
1814
1815 ! Can make this more robust
1816 inv_rho = 1.0d0 / w(ixo^s, rho_)
1817 end function rhd_inv_rho
1818
1819 subroutine rhd_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
1820 ! handles hydro (density,pressure,velocity) bootstrapping
1821 ! any negative dust density is flagged as well (and throws an error)
1822 ! small_values_method=replace also for dust
1826 logical, intent(in) :: primitive
1827 integer, intent(in) :: ixi^l,ixo^l
1828 double precision, intent(inout) :: w(ixi^s,1:nw)
1829 double precision, intent(in) :: x(ixi^s,1:ndim)
1830 character(len=*), intent(in) :: subname
1831
1832 integer :: n,idir
1833 logical :: flag(ixi^s,1:nw)
1834
1835 call rhd_check_w(primitive, ixi^l, ixo^l, w, flag)
1836
1837 if (any(flag)) then
1838 select case (small_values_method)
1839 case ("replace")
1840 where(flag(ixo^s,rho_)) w(ixo^s,rho_) = small_density
1841 do idir = 1, ndir
1842 if(small_values_fix_iw(mom(idir))) then
1843 where(flag(ixo^s,rho_)) w(ixo^s, mom(idir)) = 0.0d0
1844 end if
1845 end do
1846
1847 if (small_values_fix_iw(r_e)) then
1848 where(flag(ixo^s,r_e)) w(ixo^s,r_e) = small_r_e
1849 end if
1850
1851 if(rhd_energy)then
1852 if(small_values_fix_iw(e_)) then
1853 if(primitive) then
1854 where(flag(ixo^s,rho_)) w(ixo^s, p_) = small_pressure
1855 else
1856 where(flag(ixo^s,rho_)) w(ixo^s, e_) = small_e + rhd_kin_en(w,ixi^l,ixo^l)
1857 endif
1858 end if
1859 endif
1860
1861 if(rhd_energy) then
1862 if(primitive) then
1863 where(flag(ixo^s,e_)) w(ixo^s,p_) = small_pressure
1864 else
1865 where(flag(ixo^s,e_))
1866 ! Add kinetic energy
1867 w(ixo^s,e_) = small_e + rhd_kin_en(w,ixi^l,ixo^l)
1868 end where
1869 end if
1870 end if
1871
1872 if(rhd_dust)then
1873 do n=1,dust_n_species
1874 where(flag(ixo^s,dust_rho(n))) w(ixo^s,dust_rho(n)) = 0.0d0
1875 do idir = 1, ndir
1876 where(flag(ixo^s,dust_rho(n))) w(ixo^s,dust_mom(idir,n)) = 0.0d0
1877 enddo
1878 enddo
1879 endif
1880 case ("average")
1881 if(primitive)then
1882 ! averaging for all primitive fields, including dust
1883 call small_values_average(ixi^l, ixo^l, w, x, flag)
1884 else
1885 ! do averaging of density
1886 call small_values_average(ixi^l, ixo^l, w, x, flag, rho_)
1887 if(rhd_energy) then
1888 ! do averaging of pressure
1889 w(ixi^s,p_)=(rhd_gamma-1.d0)*(w(ixi^s,e_) &
1890 -0.5d0*sum(w(ixi^s, mom(:))**2, dim=ndim+1)/w(ixi^s,rho_))
1891 call small_values_average(ixi^l, ixo^l, w, x, flag, p_)
1892 w(ixi^s,e_)=w(ixi^s,p_)/(rhd_gamma-1.d0) &
1893 +0.5d0*sum(w(ixi^s, mom(:))**2, dim=ndim+1)/w(ixi^s,rho_)
1894 end if
1895 ! do averaging of radiation energy
1896 call small_values_average(ixi^l, ixo^l, w, x, flag, r_e)
1897 if(rhd_dust)then
1898 do n=1,dust_n_species
1899 where(flag(ixo^s,dust_rho(n))) w(ixo^s,dust_rho(n)) = 0.0d0
1900 do idir = 1, ndir
1901 where(flag(ixo^s,dust_rho(n))) w(ixo^s,dust_mom(idir,n)) = 0.0d0
1902 enddo
1903 enddo
1904 endif
1905 endif
1906 case default
1907 if(.not.primitive) then
1908 !convert w to primitive
1909 ! Calculate pressure = (gamma-1) * (e-ek)
1910 if(rhd_energy) then
1911 w(ixo^s,p_)=(rhd_gamma-1.d0)*(w(ixo^s,e_)-rhd_kin_en(w,ixi^l,ixo^l))
1912 end if
1913 ! Convert gas momentum to velocity
1914 do idir = 1, ndir
1915 w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))/w(ixo^s,rho_)
1916 end do
1917 end if
1918 ! NOTE: dust entries may still have conserved values here
1919 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
1920 end select
1921 end if
1922 end subroutine rhd_handle_small_values
1923
1924 subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
1927 integer, intent(in) :: ixi^l, ixo^l
1928 double precision, intent(in) :: w(ixi^s,1:nw)
1929 double precision, intent(in) :: x(ixi^s,1:ndim)
1930 double precision, intent(out):: rfactor(ixi^s)
1931
1932 double precision :: iz_h(ixo^s),iz_he(ixo^s)
1933
1934 call ionization_degree_from_temperature(ixi^l,ixo^l,w(ixi^s,te_),iz_h,iz_he)
1935 ! assume the first and second ionization of Helium have the same degree
1936 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
1937
1938 end subroutine rfactor_from_temperature_ionization
1939
1940 subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
1942 integer, intent(in) :: ixi^l, ixo^l
1943 double precision, intent(in) :: w(ixi^s,1:nw)
1944 double precision, intent(in) :: x(ixi^s,1:ndim)
1945 double precision, intent(out):: rfactor(ixi^s)
1946
1947 rfactor(ixo^s)=rr
1948
1949 end subroutine rfactor_from_constant_ionization
1950
1951 subroutine rhd_update_temperature(ixI^L,ixO^L,wCT,w,x)
1954
1955 integer, intent(in) :: ixi^l, ixo^l
1956 double precision, intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:ndim)
1957 double precision, intent(inout) :: w(ixi^s,1:nw)
1958
1959 double precision :: iz_h(ixo^s),iz_he(ixo^s), pth(ixi^s)
1960
1961 call ionization_degree_from_temperature(ixi^l,ixo^l,wct(ixi^s,te_),iz_h,iz_he)
1962
1963 call rhd_get_pthermal(w,x,ixi^l,ixo^l,pth)
1964
1965 w(ixo^s,te_)=(2.d0+3.d0*he_abundance)*pth(ixo^s)/(w(ixo^s,rho_)*(1.d0+iz_h(ixo^s)+&
1966 he_abundance*(iz_he(ixo^s)*(iz_he(ixo^s)+1.d0)+1.d0)))
1967
1968 end subroutine rhd_update_temperature
1969
1970end module mod_rhd_phys
Calculate w(iw)=w(iw)+qdt*SOURCE[wCT,qtC,x] within ixO for all indices iw=iwmin......
Module for including anisotropic flux limited diffusion (AFLD)-approximation in Radiation-hydrodynami...
Definition mod_afld.t:8
subroutine afld_get_diffcoef_central(w, wct, x, ixil, ixol)
Calculates cell-centered diffusion coefficient to be used in multigrid.
Definition mod_afld.t:684
subroutine, public get_afld_rad_force(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 This subroutine handles th...
Definition mod_afld.t:141
subroutine, public afld_init(he_abundance, rhd_radiation_diffusion, afld_gamma)
Initialising FLD-module: Read opacities Initialise Multigrid adimensionalise kappa Add extra variable...
Definition mod_afld.t:93
subroutine, public afld_radforce_get_dt(w, ixil, ixol, dtnew, dxd, x)
Definition mod_afld.t:217
subroutine, public afld_get_radpress(w, x, ixil, ixol, rad_pressure, nth)
Calculate Radiation Pressure Returns Radiation Pressure as tensor.
Definition mod_afld.t:518
subroutine, public get_afld_energy_interact(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 This subroutine handles th...
Definition mod_afld.t:244
Module with basic data types used in amrvac.
integer, parameter std_len
Default length for strings.
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.
double precision, parameter const_rad_a
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_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
Module for flux conservation near refinement boundaries.
Nicolas Moens Module for including flux limited diffusion (FLD)-approximation in Radiation-hydrodynam...
Definition mod_fld.t:9
double precision, public fld_mu
mean particle mass
Definition mod_fld.t:19
subroutine, public get_fld_rad_force(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 This subroutine handles th...
Definition mod_fld.t:146
subroutine, public fld_get_radpress(w, x, ixil, ixol, rad_pressure, nth)
Calculate Radiation Pressure Returns Radiation Pressure as tensor.
Definition mod_fld.t:456
subroutine, public fld_init(he_abundance, radiation_diffusion, energy_interact, r_gamma)
Initialising FLD-module: Read opacities Initialise Multigrid adimensionalise kappa Add extra variable...
Definition mod_fld.t:94
subroutine fld_get_diffcoef_central(w, wct, x, ixil, ixol)
Calculates cell-centered diffusion coefficient to be used in multigrid.
Definition mod_fld.t:719
character(len=8) fld_diff_scheme
Which method to solve diffusion part.
Definition mod_fld.t:35
subroutine, public fld_radforce_get_dt(w, ixil, ixol, dtnew, dxd, x)
Definition mod_fld.t:181
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.
integer, parameter bc_noinflow
double precision unit_time
Physical scaling factor for time.
double precision unit_density
Physical scaling factor for density.
double precision unit_opacity
Physical scaling factor for Opacity.
integer, parameter unitpar
file handle for IO
double precision global_time
The global simulation time.
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.
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
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.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision unit_velocity
Physical scaling factor for velocity.
double precision unit_temperature
Physical scaling factor for temperature.
double precision unit_radflux
Physical scaling factor for radiation flux.
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
integer nghostcells
Number of ghost cells surrounding a grid.
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
logical use_multigrid
Use multigrid (only available in 2D and 3D)
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
module ionization degree - get ionization degree for given temperature
subroutine ionization_degree_from_temperature(ixil, ixol, te, iz_h, iz_he)
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
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)
Radiation-Hydrodynamics physics module Module aims at solving the Hydrodynamic equations toghether wi...
logical, public, protected rhd_radiative_cooling
Whether radiative cooling is added.
subroutine, public rhd_to_primitive(ixil, ixol, w, x)
Transform conservative variables into primitive ones.
integer, public, protected rhd_n_tracer
Number of tracer species.
integer, public, protected p_
Index of the gas pressure (-1 if not present) should equal e_.
type(tc_fluid), allocatable, public tc_fl
logical, public, protected rhd_dust
Whether dust is added.
integer, public, protected rhd_trac_type
integer, public, protected te_
Indices of temperature.
logical, public, protected rhd_rotating_frame
Whether rotating frame is activated.
logical, public, protected rhd_viscosity
Whether viscosity is added.
double precision, public kbmpmua4
kb/(m_p mu)* 1/a_rad**4,
subroutine, public rhd_check_params
logical, public, protected rhd_trac
Whether TRAC method is used.
subroutine, public rhd_phys_init()
Initialize the module.
character(len=8), public rhd_radiation_formalism
Formalism to treat radiation.
double precision function, dimension(ixo^s), public rhd_kin_en(w, ixil, ixol, inv_rho)
integer, public, protected rho_
Index of the density (in the w array)
integer, public, protected r_e
Index of the radiation energy.
logical, public, protected rhd_radiation_diffusion
Treat radiation energy diffusion.
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
character(len=8), public rhd_pressure
In the case of no rhd_energy, how to compute pressure.
subroutine, public rhd_get_trad(w, x, ixil, ixol, trad)
Calculates radiation temperature.
subroutine, public rhd_check_w(primitive, ixil, ixol, w, flag)
Returns logical argument flag where values are ok.
logical, public, protected rhd_partial_ionization
Whether plasma is partially ionized.
logical, public, protected rhd_radiation_force
Treat radiation fld_Rad_force.
logical, public, protected rhd_energy
Whether an energy equation is used.
double precision, public rhd_adiab
The adiabatic constant.
subroutine, public rhd_get_tgas(w, x, ixil, ixol, tgas)
Calculates gas temperature.
double precision, public, protected small_r_e
The smallest allowed radiation energy.
subroutine, public rhd_get_ptot(w, x, ixil, ixol, ptot)
calculates the sum of the gas pressure and max Prad tensor element
double precision, public, protected h_ion_fr
Ionization fraction of H H_ion_fr = H+/(H+ + H)
integer, dimension(:), allocatable, public, protected tracer
Indices of the tracers.
subroutine, public rhd_to_conserved(ixil, ixol, w, x)
Transform primitive variables into conservative ones.
integer, public, protected e_
Index of the energy density (-1 if not present)
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
subroutine, public rhd_set_mg_bounds
Set the boundaries for the diffusion of E.
subroutine, public rhd_get_csound2(w, x, ixil, ixol, csound2)
Calculate the square of the thermal sound speed csound2 within ixO^L. csound2=gamma*p/rho.
double precision, public, protected he_ion_fr2
Ratio of number He2+ / number He+ + He2+ He_ion_fr2 = He2+/(He2+ + He+)
logical, public, protected rhd_thermal_conduction
Whether thermal conduction is added.
double precision, public, protected rr
subroutine, public rhd_get_pthermal(w, x, ixil, ixol, pth)
Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho) within ixO^L.
logical, public, protected rhd_radiation_advection
Treat radiation advection.
logical, public, protected rhd_particles
Whether particles module is added.
double precision, public, protected he_ion_fr
Ionization fraction of He He_ion_fr = (He2+ + He+)/(He2+ + He+ + He)
type(te_fluid), allocatable, public te_fl_rhd
logical, public, protected rhd_energy_interact
Treat radiation-gas energy interaction.
logical, public, protected eq_state_units
type(rc_fluid), allocatable, public rc_fl
logical, public, protected rhd_gravity
Whether gravity is added.
integer, public, protected tcoff_
Index of the cutoff temperature for the TRAC method.
subroutine, public rhd_get_pradiation(w, x, ixil, ixol, prad)
Calculate radiation pressure within ixO^L.
double precision, public rhd_gamma
The adiabatic index.
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)
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(special_mg_bc), pointer usr_special_mg_bc
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)