diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index 5175e7d47fb90560549aca43f4686aa399fb7fc6..0a81716a3e34bbb305916b1148a9d8199b84fa5a 100644 --- a/examples/parameter_example.yml +++ b/examples/parameter_example.yml @@ -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 diff --git a/src/hydro_properties.c b/src/hydro_properties.c index bc324ef27c4f07021bdbf3f47d8ee99f982e9b68..e4135a762472bcf734676a3f294cd401770f9f01 100644 --- a/src/hydro_properties.c +++ b/src/hydro_properties.c @@ -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); diff --git a/src/hydro_properties.h b/src/hydro_properties.h index 14500a3fa228b48b4ab7de86c218144e528035d2..2799f6c86eec7f0ebf140643c5ab7fa9b60e6273 100644 --- a/src/hydro_properties.h +++ b/src/hydro_properties.h @@ -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; diff --git a/src/space.c b/src/space.c index ea599a45d3957c2ea224c0e7c3278f6ddfbe163e..4dc1369882675f7ec4f4f1aa0d136bd85711345a 100644 --- a/src/space.c +++ b/src/space.c @@ -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);