Commit df9ca55b authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Read the minimally allowed temperature from the parameter file and convert it...

Read the minimally allowed temperature from the parameter file and convert it to a minimal internal energy per unit mass.
parent f0db61e8
......@@ -30,7 +30,8 @@ SPH:
h_max: 10. # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
max_volume_change: 1.4 # (Optional) Maximal allowed change of kernel volume over one time-step.
max_ghost_iterations: 30 # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
initial_temperature: 100 # (Optional) Initial temperature (in internal units) to set the gass particles to. Value is ignored if set to 0.
initial_temperature: 0 # (Optional) Initial temperature (in internal units) to set the gas particles at start-up. Value is ignored if set to 0.
minimal_temperature: 0 # (Optional) Minimal temperature (in internal units) allowed for the gas particles. Value is ignored if set to 0.
H_mass_fraction: 0.76 # (Optional) Hydrogen mass fraction used for initial conversion from temp to internal energy.
# Parameters for the self-gravity scheme
......
......@@ -37,6 +37,7 @@
#define hydro_props_default_h_max FLT_MAX
#define hydro_props_default_h_tolerance 1e-4
#define hydro_props_default_init_temp 0.f
#define hydro_props_default_min_temp 0.f
#define hydro_props_default_H_fraction 0.76
void hydro_props_init(struct hydro_props *p,
......@@ -82,6 +83,13 @@ void hydro_props_init(struct hydro_props *p,
p->initial_temperature = parser_get_opt_param_float(
params, "SPH:initial_temperature", hydro_props_default_init_temp);
/* Initial temperature */
p->minimal_temperature = parser_get_opt_param_float(
params, "SPH:minimal_temperature", hydro_props_default_min_temp);
if (p->initial_temperature < p->minimal_temperature)
error("Initial temperature lower than minimal allowed temperature!");
/* Hydrogen mass fraction */
p->hydrogen_mass_fraction = parser_get_opt_param_double(
params, "SPH:H_mass_fraction", hydro_props_default_H_fraction);
......@@ -101,6 +109,21 @@ void hydro_props_init(struct hydro_props *p,
mean_molecular_weight = 4. / (1. + 3. * p->hydrogen_mass_fraction);
p->initial_internal_energy = u_init / mean_molecular_weight;
/* Compute the minimal energy (Note the temp. read is in internal units) */
double u_min = phys_const->const_boltzmann_k / phys_const->const_proton_mass;
u_min *= p->initial_temperature;
u_min *= hydro_one_over_gamma_minus_one;
/* Correct for hydrogen mass fraction */
if (p->minimal_temperature *
units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE) >
1e4)
mean_molecular_weight = 4. / (8. - 5. * (1. - p->hydrogen_mass_fraction));
else
mean_molecular_weight = 4. / (1. + 3. * p->hydrogen_mass_fraction);
p->minimal_internal_energy = u_min / mean_molecular_weight;
}
void hydro_props_print(const struct hydro_props *p) {
......@@ -134,6 +157,9 @@ void hydro_props_print(const struct hydro_props *p) {
if (p->initial_temperature != hydro_props_default_init_temp)
message("Initial gas temperature set to %f", p->initial_temperature);
if (p->minimal_temperature != hydro_props_default_min_temp)
message("Minimal gas temperature set to %f", p->initial_temperature);
}
#if defined(HAVE_HDF5)
......@@ -156,6 +182,7 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
pow_dimension(expf(p->log_max_h_change)));
io_write_attribute_i(h_grpsph, "Max ghost iterations",
p->max_smoothing_iterations);
io_write_attribute_f(h_grpsph, "Minimal temperature", p->minimal_temperature);
io_write_attribute_f(h_grpsph, "Initial temperature", p->initial_temperature);
io_write_attribute_f(h_grpsph, "Initial energy per unit mass",
p->initial_internal_energy);
......
......@@ -66,6 +66,12 @@ struct hydro_props {
/*! Maximal change of h over one time-step */
float log_max_h_change;
/*! Minimal temperature allowed */
float minimal_temperature;
/*! Minimal internal energy per unit mass */
float minimal_internal_energy;
/*! Initial temperature */
float initial_temperature;
......
......@@ -2652,6 +2652,7 @@ void space_first_init_parts(struct space *s,
const struct hydro_props *hydro_props = s->e->hydro_properties;
const float u_init = hydro_props->initial_internal_energy;
const float u_min = hydro_props->minimal_internal_energy;
for (size_t i = 0; i < nr_parts; ++i) {
......@@ -2674,6 +2675,7 @@ void space_first_init_parts(struct space *s,
/* Overwrite the internal energy? */
if (u_init > 0.f) hydro_set_init_internal_energy(&p[i], u_init);
if (u_min > 0.f) hydro_set_init_internal_energy(&p[i], u_min);
/* Also initialise the chemistry */
chemistry_first_init_part(&p[i], &xp[i], chemistry);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment