diff --git a/configure.ac b/configure.ac index f551aaee9b82cce57644b38a2c46b7f8ae77cc32..45c0ab52248f05f45115ddd12e9a4120bbd38fcd 100644 --- a/configure.ac +++ b/configure.ac @@ -1327,6 +1327,7 @@ with_subgrid_cooling=none with_subgrid_chemistry=none with_subgrid_tracers=none with_subgrid_entropy_floor=none +with_subgrid_pressure_floor=none with_subgrid_stars=none with_subgrid_star_formation=none with_subgrid_feedback=none @@ -1338,10 +1339,9 @@ case "$with_subgrid" in none) ;; GEAR) - with_subgrid_cooling=grackle + with_subgrid_cooling=grackle_0 with_subgrid_chemistry=GEAR - with_subgrid_tracers=none - with_subgrid_entropy_floor=none + with_subgrid_pressure_floor=GEAR with_subgrid_stars=GEAR with_subgrid_star_formation=GEAR with_subgrid_feedback=none @@ -1646,7 +1646,8 @@ esac # Cooling function AC_ARG_WITH([cooling], [AS_HELP_STRING([--with-cooling=<function>], - [cooling function @<:@none, const-du, const-lambda, EAGLE, grackle, grackle1, grackle2, grackle3 default: none@:>@] + [cooling function @<:@none, const-du, const-lambda, EAGLE, grackle_* default: none@:>@. + For Grackle, you need to provide the primordial chemistry parameter (e.g. grackle_0)] )], [with_cooling="$withval"], [with_cooling="none"] @@ -1673,21 +1674,10 @@ case "$with_cooling" in compton) AC_DEFINE([COOLING_COMPTON], [1], [Compton cooling off the CMB]) ;; - grackle) + grackle_*) AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library]) - AC_DEFINE([COOLING_GRACKLE_MODE], [0], [Grackle chemistry network, mode 0]) - ;; - grackle1) - AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library]) - AC_DEFINE([COOLING_GRACKLE_MODE], [1], [Grackle chemistry network, mode 1]) - ;; - grackle2) - AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library]) - AC_DEFINE([COOLING_GRACKLE_MODE], [2], [Grackle chemistry network, mode 2]) - ;; - grackle3) - AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library]) - AC_DEFINE([COOLING_GRACKLE_MODE], [3], [Grackle chemistry network, mode 3]) + primordial_chemistry=${with_cooling:8} + AC_DEFINE_UNQUOTED([COOLING_GRACKLE_MODE], [$primordial_chemistry], [Grackle chemistry network]) ;; EAGLE) AC_DEFINE([COOLING_EAGLE], [1], [Cooling following the EAGLE model]) @@ -1922,6 +1912,35 @@ case "$with_entropy_floor" in ;; esac +# Pressure floor +AC_ARG_WITH([pressure-floor], + [AS_HELP_STRING([--with-pressure-floor=<floor>], + [pressure floor @<:@none, GEAR, default: none@:>@ + The hydro model needs to be compatible.] + )], + [with_pressure_floor="$withval"], + [with_pressure_floor="none"] +) +if test "$with_subgrid" != "none"; then + if test "$with_pressure_floor" != "none"; then + AC_MSG_ERROR([Cannot provide with-subgrid and with-pressure-floor together]) + else + with_pressure_floor="$with_subgrid_pressure_floor" + fi +fi + +case "$with_pressure_floor" in + none) + AC_DEFINE([PRESSURE_FLOOR_NONE], [1], [No pressure floor]) + ;; + GEAR) + AC_DEFINE([PRESSURE_FLOOR_GEAR], [1], [GEAR pressure floor]) + ;; + *) + AC_MSG_ERROR([Unknown pressure floor model]) + ;; +esac + # Star formation AC_ARG_WITH([star-formation], [AS_HELP_STRING([--with-star-formation=<sfm>], @@ -2056,6 +2075,7 @@ AC_MSG_RESULT([ External potential : $with_potential Entropy floor : $with_entropy_floor + Pressure floor : $with_pressure_floor Cooling function : $with_cooling Chemistry : $with_chemistry Tracers : $with_tracers diff --git a/doc/RTD/source/SubgridModels/GEAR/index.rst b/doc/RTD/source/SubgridModels/GEAR/index.rst index 2a211759bfc4895fd07279b72f78200d6ea47546..152056ef2eb60da189bf18b3b3c9909f88ace6de 100644 --- a/doc/RTD/source/SubgridModels/GEAR/index.rst +++ b/doc/RTD/source/SubgridModels/GEAR/index.rst @@ -5,6 +5,29 @@ GEAR model =========== +Pressure Floor +~~~~~~~~~~~~~~ + +In order to avoid the artificial collapse of unresolved clumps, a minimum in pressure is applied to the particles. +This additional pressure can be seen as the pressure due to unresolved hydrodynamics turbulence and is given by: + +.. math:: + P_\textrm{Jeans} = \frac{\rho}{\gamma} \left( \frac{4}{\pi} G h^2 \rho N_\textrm{Jeans}^{2/3} - \sigma^2 \right) + +where :math:`\rho` is the density, :math:`\gamma` the adiabatic index, :math:`G` is the gravitational constant, +:math:`h` the smoothing length, :math:`N_\textrm{Jeans}` is the number of particle required in order to resolve a clump and +:math:`\sigma` the velocity dispersion. + + +This must be directly implemented into the hydro schemes, therefore only a subset of schemes (Gadget-2 and Pressure-Energy) are currently implemented. +In order to implement it, you need equation 12 in `Hopkins 2013 <https://arxiv.org/abs/1206.5006>`_: + +.. math:: + m_i \frac{\mathrm{d}v_i}{\mathrm{d}t} = - \sum_j x_i x_j \left[ \frac{P_i}{y_i^2} f_{ij} \nabla_i W_{ij}(h_i) + \frac{P_j}{y_j^2} f_{ji} \nabla_j W_{ji}(h_j) \right] + +and simply replace the :math:`P_i, P_j` by the pressure with the floor. +Here the :math:`x, y` are simple weights that should never have the pressure floor included even if they are related to the pressure (e.g. pressure-entropy). + Cooling: Grackle ~~~~~~~~~~~~~~~~ diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml index b108290fcd146461827d5858742bd0971ac66945..94a9082af95e96b161cf5fe469eb082774e5007f 100644 --- a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml +++ b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml @@ -99,3 +99,6 @@ GrackleCooling: GearChemistry: InitialMetallicity: 0.01295 + +GEARPressureFloor: + Jeans_factor: 10. # Number of particles required to suppose a resolved clump and avoid the pressure floor. diff --git a/examples/main.c b/examples/main.c index b3577b0e2adec7650bb6acd73b128552071f47d8..491b1ac1e72895ba7f9028d7c0b1bf9181f70360 100644 --- a/examples/main.c +++ b/examples/main.c @@ -777,6 +777,13 @@ int main(int argc, char *argv[]) { else bzero(&entropy_floor, sizeof(struct entropy_floor_properties)); + /* Initialise the pressure floor */ + if (with_hydro) + pressure_floor_init(&pressure_floor_props, &prog_const, &us, + &hydro_properties, params); + else + bzero(&pressure_floor_props, sizeof(struct pressure_floor_properties)); + /* Initialise the stars properties */ if (with_stars) stars_props_init(&stars_properties, &prog_const, &us, params, diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index 246c6be2a8dd12d321c9c85a4da40cad8b8d03cf..5ef0679d072816b7ee0e042d5e473ec6d02a7599 100644 --- a/examples/parameter_example.yml +++ b/examples/parameter_example.yml @@ -285,6 +285,11 @@ EAGLEEntropyFloor: Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin. Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor +# Parameters related to pressure floors ---------------------------------------------- + +GEARPressureFloor: + Jeans_factor: 10. # Number of particles required to suppose a resolved clump and avoid the pressure floor. + # Parameters related to cooling function ---------------------------------------------- # Constant du/dt cooling function diff --git a/src/Makefile.am b/src/Makefile.am index 0f7cef5e3c6b861ebde639d90624a360f1be7047..ecf8399e5293689a3e00262b91bd2ebf6a828e2d 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -79,7 +79,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c engine_maketasks.c collectgroup.c hydro_space.c equation_of_state.c \ chemistry.c cosmology.c restart.c mesh_gravity.c velociraptor_interface.c \ outputlist.c velociraptor_dummy.c logger_io.c memuse.c fof.c \ - hashmap.c \ + hashmap.c pressure_floor.c \ $(EAGLE_COOLING_SOURCES) $(EAGLE_FEEDBACK_SOURCES) # Include files for distribution, not installation. diff --git a/src/cooling/grackle/cooling.h b/src/cooling/grackle/cooling.h index 57f4b146c451557365e08d260ddd45276844af9c..98d34dc630aa0c83198ddd60b6ee9e2775778d0f 100644 --- a/src/cooling/grackle/cooling.h +++ b/src/cooling/grackle/cooling.h @@ -52,7 +52,7 @@ static gr_float cooling_time( const struct unit_system* restrict us, const struct cosmology* restrict cosmo, const struct cooling_function_data* restrict cooling, - const struct part* restrict p, const struct xpart* restrict xp); + const struct part* restrict p, struct xpart* restrict xp); static gr_float cooling_new_energy( const struct phys_const* restrict phys_const, const struct unit_system* restrict us, @@ -573,7 +573,7 @@ __attribute__((always_inline)) INLINE static gr_float cooling_time( const struct unit_system* restrict us, const struct cosmology* restrict cosmo, const struct cooling_function_data* restrict cooling, - const struct part* restrict p, const struct xpart* restrict xp) { + const struct part* restrict p, struct xpart* restrict xp) { /* set current time */ code_units units = cooling->units; diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h index afc9c17d5b6ed691685bf42df2562ecb6f4409d8..93c0e58068341ab29a7ed0fb6b0d21f76980d952 100644 --- a/src/hydro/Gadget2/hydro.h +++ b/src/hydro/Gadget2/hydro.h @@ -41,6 +41,7 @@ #include "hydro_space.h" #include "kernel_hydro.h" #include "minmax.h" +#include "pressure_floor.h" #include "./hydro_parameters.h" @@ -109,7 +110,8 @@ hydro_get_drifted_physical_internal_energy(const struct part *restrict p, __attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure( const struct part *restrict p) { - return gas_pressure_from_entropy(p->rho, p->entropy); + const float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy); + return pressure_floor_get_pressure(p, p->rho, comoving_pressure); } /** @@ -121,7 +123,10 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure( __attribute__((always_inline)) INLINE static float hydro_get_physical_pressure( const struct part *restrict p, const struct cosmology *cosmo) { - return gas_pressure_from_entropy(p->rho * cosmo->a3_inv, p->entropy); + const float phys_pressure = + gas_pressure_from_entropy(p->rho * cosmo->a3_inv, p->entropy); + const float phys_rho = hydro_get_physical_density(p, cosmo); + return pressure_floor_get_pressure(p, phys_rho, phys_pressure); } /** @@ -388,13 +393,15 @@ hydro_set_drifted_physical_internal_energy(struct part *p, const float rho_inv = 1.f / p->rho; /* Compute the pressure */ - const float pressure = gas_pressure_from_entropy(p->rho, p->entropy); + float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy); + comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure); /* Compute the sound speed */ - const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); + const float soundspeed = + gas_soundspeed_from_pressure(p->rho, comoving_pressure); /* Divide the pressure by the density squared to get the SPH term */ - const float P_over_rho2 = pressure * rho_inv * rho_inv; + const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv; /* Update variables. */ p->force.P_over_rho2 = P_over_rho2; @@ -591,13 +598,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( const float abs_div_physical_v = fabsf(div_physical_v); /* Compute the pressure */ - const float pressure = gas_pressure_from_entropy(p->rho, p->entropy); + float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy); + comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure); /* Compute the sound speed */ - const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); + const float soundspeed = + gas_soundspeed_from_pressure(p->rho, comoving_pressure); /* Divide the pressure by the density squared to get the SPH term */ - const float P_over_rho2 = pressure * rho_inv * rho_inv; + const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv; /* Compute the Balsara switch */ /* Pre-multiply in the AV factor; hydro_props are not passed to the iact @@ -665,14 +674,16 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( p->entropy = xp->entropy_full; /* Re-compute the pressure */ - const float pressure = gas_pressure_from_entropy(p->rho, p->entropy); + float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy); + comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure); /* Compute the new sound speed */ - const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); + const float soundspeed = + gas_soundspeed_from_pressure(p->rho, comoving_pressure); /* Divide the pressure by the density squared to get the SPH term */ const float rho_inv = 1.f / p->rho; - const float P_over_rho2 = pressure * rho_inv * rho_inv; + const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv; /* Update variables */ p->force.soundspeed = soundspeed; @@ -730,14 +741,16 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( p->rho *= expf(w2); /* Re-compute the pressure */ - const float pressure = gas_pressure_from_entropy(p->rho, p->entropy); + float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy); + comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure); /* Compute the new sound speed */ - const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); + const float soundspeed = + gas_soundspeed_from_pressure(p->rho, comoving_pressure); /* Divide the pressure by the density squared to get the SPH term */ const float rho_inv = 1.f / p->rho; - const float P_over_rho2 = pressure * rho_inv * rho_inv; + const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv; /* Update variables */ p->force.soundspeed = soundspeed; @@ -840,14 +853,16 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( } /* Compute the pressure */ - const float pressure = gas_pressure_from_entropy(p->rho, p->entropy); + float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy); + comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure); /* Compute the sound speed */ - const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); + const float soundspeed = + gas_soundspeed_from_pressure(p->rho, comoving_pressure); /* Divide the pressure by the density squared to get the SPH term */ const float rho_inv = 1.f / p->rho; - const float P_over_rho2 = pressure * rho_inv * rho_inv; + const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv; p->force.soundspeed = soundspeed; p->force.P_over_rho2 = P_over_rho2; diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h index e3a4f98160150717e50651760d8af68af2a70662..15f9a958a6c5d4fc199a1c1d7d85d3f909520d23 100644 --- a/src/hydro/PressureEnergy/hydro.h +++ b/src/hydro/PressureEnergy/hydro.h @@ -46,6 +46,7 @@ #include "hydro_space.h" #include "kernel_hydro.h" #include "minmax.h" +#include "pressure_floor.h" #include "./hydro_parameters.h" @@ -221,7 +222,9 @@ hydro_get_comoving_soundspeed(const struct part *restrict p) { /* Compute the sound speed -- see theory section for justification */ /* IDEAL GAS ONLY -- P-U does not work with generic EoS. */ - const float square_rooted = sqrtf(hydro_gamma * p->pressure_bar / p->rho); + const float comoving_pressure = + pressure_floor_get_pressure(p, p->rho, p->pressure_bar); + const float square_rooted = sqrtf(hydro_gamma * comoving_pressure / p->rho); return square_rooted; } @@ -236,7 +239,10 @@ __attribute__((always_inline)) INLINE static float hydro_get_physical_soundspeed(const struct part *restrict p, const struct cosmology *cosmo) { - return cosmo->a_factor_sound_speed * p->force.soundspeed; + const float phys_rho = hydro_get_physical_density(p, cosmo); + + return pressure_floor_get_pressure( + p, phys_rho, cosmo->a_factor_sound_speed * p->force.soundspeed); } /** @@ -640,10 +646,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( hydro_one_over_gamma_minus_one) / (1.f + common_factor * p->density.wcount_dh); + /* Get the pressures */ + const float comoving_pressure_with_floor = + pressure_floor_get_pressure(p, p->rho, p->pressure_bar); + /* Update variables. */ p->force.f = grad_h_term; p->force.soundspeed = soundspeed; p->force.balsara = balsara; + p->force.pressure_bar_with_floor = comoving_pressure_with_floor; } /** @@ -755,6 +766,11 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( const float soundspeed = hydro_get_comoving_soundspeed(p); p->force.soundspeed = soundspeed; + + /* update the required variables */ + const float comoving_pressure_with_floor = + pressure_floor_get_pressure(p, p->rho, p->pressure_bar); + p->force.pressure_bar_with_floor = comoving_pressure_with_floor; } /** diff --git a/src/hydro/PressureEnergy/hydro_iact.h b/src/hydro/PressureEnergy/hydro_iact.h index b59fdbbe418d0442fd290b3596fd08531fae88e6..aab8237c7b26c50d9e04610c6c1029f97e5d73ed 100644 --- a/src/hydro/PressureEnergy/hydro_iact.h +++ b/src/hydro/PressureEnergy/hydro_iact.h @@ -259,11 +259,18 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( /* Convolve with the kernel */ const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv; + /* Compute the ratio of pressures */ + const float pressure_inverse_i = + pi->force.pressure_bar_with_floor / (pi->pressure_bar * pi->pressure_bar); + const float pressure_inverse_j = + pj->force.pressure_bar_with_floor / (pj->pressure_bar * pj->pressure_bar); + /* SPH acceleration term */ - const float sph_acc_term = - pj->u * pi->u * hydro_gamma_minus_one * hydro_gamma_minus_one * - ((f_ij / pi->pressure_bar) * wi_dr + (f_ji / pj->pressure_bar) * wj_dr) * - r_inv; + const float sph_acc_term = pj->u * pi->u * hydro_gamma_minus_one * + hydro_gamma_minus_one * + ((f_ij * pressure_inverse_i) * wi_dr + + (f_ji * pressure_inverse_j) * wj_dr) * + r_inv; /* Assemble the acceleration */ const float acc = sph_acc_term + visc_acc_term; @@ -278,11 +285,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( pj->a_hydro[2] += mi * acc * dx[2]; /* Get the time derivative for u. */ + const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one * - pj->u * pi->u * (f_ij / pi->pressure_bar) * + pj->u * pi->u * (f_ij * pressure_inverse_i) * wi_dr * dvdr * r_inv; + const float sph_du_term_j = hydro_gamma_minus_one * hydro_gamma_minus_one * - pi->u * pj->u * (f_ji / pj->pressure_bar) * + pi->u * pj->u * (f_ji * pressure_inverse_j) * wj_dr * dvdr * r_inv; /* Viscosity term */ @@ -386,11 +395,18 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( /* Convolve with the kernel */ const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv; + /* Compute the ratio of pressures */ + const float pressure_inverse_i = + pi->force.pressure_bar_with_floor / (pi->pressure_bar * pi->pressure_bar); + const float pressure_inverse_j = + pj->force.pressure_bar_with_floor / (pj->pressure_bar * pj->pressure_bar); + /* SPH acceleration term */ - const float sph_acc_term = - pj->u * pi->u * hydro_gamma_minus_one * hydro_gamma_minus_one * - ((f_ij / pi->pressure_bar) * wi_dr + (f_ji / pj->pressure_bar) * wj_dr) * - r_inv; + const float sph_acc_term = pj->u * pi->u * hydro_gamma_minus_one * + hydro_gamma_minus_one * + ((f_ij * pressure_inverse_i) * wi_dr + + (f_ji * pressure_inverse_j) * wj_dr) * + r_inv; /* Assemble the acceleration */ const float acc = sph_acc_term + visc_acc_term; @@ -402,7 +418,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( /* Get the time derivative for u. */ const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one * - pj->u * pi->u * (f_ij / pi->pressure_bar) * + pj->u * pi->u * (f_ij * pressure_inverse_i) * wi_dr * dvdr * r_inv; /* Viscosity term */ diff --git a/src/hydro/PressureEnergy/hydro_part.h b/src/hydro/PressureEnergy/hydro_part.h index ded83d70176409332b669517ed203fdf95ecfd5e..c8f02f518cfde536965c0e6d000999a0e07e4aab 100644 --- a/src/hydro/PressureEnergy/hydro_part.h +++ b/src/hydro/PressureEnergy/hydro_part.h @@ -168,6 +168,9 @@ struct part { /*! Balsara switch */ float balsara; + + /*! Pressure term including the pressure floor */ + float pressure_bar_with_floor; } force; }; diff --git a/src/hydro_properties.c b/src/hydro_properties.c index 0082546ce14030a37fa63d87ea88bc99153b8213..d272fb7381664d91f763119ae68043ce59a77e1f 100644 --- a/src/hydro_properties.c +++ b/src/hydro_properties.c @@ -34,6 +34,7 @@ #include "hydro.h" #include "kernel_hydro.h" #include "parser.h" +#include "pressure_floor.h" #include "units.h" #define hydro_props_default_max_iterations 30 @@ -186,6 +187,9 @@ void hydro_props_print(const struct hydro_props *p) { /* Print equation of state first */ eos_print(&eos); + /* Then the pressure floor */ + pressure_floor_print(&pressure_floor_props); + /* Now SPH */ message("Hydrodynamic scheme: %s in %dD.", SPH_IMPLEMENTATION, (int)hydro_dimension); @@ -238,6 +242,7 @@ void hydro_props_print(const struct hydro_props *p) { void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) { eos_print_snapshot(h_grpsph, &eos); + pressure_floor_print_snapshot(h_grpsph); io_write_attribute_i(h_grpsph, "Dimension", (int)hydro_dimension); io_write_attribute_s(h_grpsph, "Scheme", SPH_IMPLEMENTATION); diff --git a/src/pressure_floor.c b/src/pressure_floor.c new file mode 100644 index 0000000000000000000000000000000000000000..2901241f1a2ff42d526a0122008dcf6ab2d8ea6f --- /dev/null +++ b/src/pressure_floor.c @@ -0,0 +1,24 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ + +/* This object's header. */ +#include "pressure_floor.h" + +/* Pressure floor for the physics model. */ +struct pressure_floor_properties pressure_floor_props; diff --git a/src/pressure_floor.h b/src/pressure_floor.h new file mode 100644 index 0000000000000000000000000000000000000000..4389dfe0891dd6bfbc1e6730cac1c6ddcc920547 --- /dev/null +++ b/src/pressure_floor.h @@ -0,0 +1,54 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_PRESSURE_FLOOR_H +#define SWIFT_PRESSURE_FLOOR_H + +/** + * @file src/pressure_floor.h + * @brief Branches between the different pressure floor models + */ + +/* Config parameters. */ +#include "../config.h" + +/* Local includes */ +#include "common_io.h" +#include "cosmology.h" +#include "error.h" +#include "inline.h" + +extern struct pressure_floor_properties pressure_floor_props; + +/* Check if pressure floor is implemented in hydro */ +#ifndef PRESSURE_FLOOR_NONE +#if defined(GADGET2_SPH) || defined(HOPKINS_PU_SPH) +/* Implemented */ +#else +#error Pressure floor not implemented with this hydro scheme +#endif + +#endif +/* Import the right pressure floor definition */ +#if defined(PRESSURE_FLOOR_NONE) +#include "./pressure_floor/none/pressure_floor.h" +#elif defined(PRESSURE_FLOOR_GEAR) +#include "./pressure_floor/GEAR/pressure_floor.h" +#endif + +#endif /* SWIFT_PRESSURE_FLOOR_H */ diff --git a/src/pressure_floor/GEAR/pressure_floor.h b/src/pressure_floor/GEAR/pressure_floor.h new file mode 100644 index 0000000000000000000000000000000000000000..dd715c155a095f9ad5f97b1d14cabfc94d9b11b0 --- /dev/null +++ b/src/pressure_floor/GEAR/pressure_floor.h @@ -0,0 +1,124 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_PRESSURE_FLOOR_GEAR_H +#define SWIFT_PRESSURE_FLOOR_GEAR_H + +#include "adiabatic_index.h" +#include "cosmology.h" +#include "equation_of_state.h" +#include "hydro_properties.h" +#include "parser.h" +#include "part.h" +#include "units.h" + +/** + * @file src/pressure_floor/GEAR/pressure_floor.h + * @brief Pressure floor used in the GEAR model + */ + +/** + * @brief Properties of the pressure floor in the GEAR model. + */ +struct pressure_floor_properties { + + /*! Jeans factor */ + float n_jeans; + + /*! The constants in internal units (4 G N_jeans^(2/3) / PI) */ + float constants; +}; + +/** + * @brief Compute the physical pressure floor of a given #part. + * + * Note that the particle is not updated!! + * + * The inputs can be either in physical or comoving coordinates (the output is + * in the same coordinates). + * + * @param p The #part. + * @param rho The physical or comoving density. + * @param pressure The physical or comoving pressure without any pressure floor. + * + * @return The physical or comoving pressure with the floor. + */ +static INLINE float pressure_floor_get_pressure(const struct part *p, + const float rho, + const float pressure) { + + /* Compute pressure floor */ + float floor = p->h * p->h * rho * pressure_floor_props.constants; + // TODO add sigma (will be done once the SF is merged) + floor *= rho * hydro_one_over_gamma; + + return fmax(pressure, floor); +} + +/** + * @brief Initialise the pressure floor by reading the parameters and converting + * to internal units. + * + * The input temperatures and number densities are converted to pressure and + * density assuming a neutral gas of primoridal abundance. + * + * @param params The YAML parameter file. + * @param us The system of units used internally. + * @param phys_const The physical constants. + * @param hydro_props The propoerties of the hydro scheme. + * @param props The pressure floor properties to fill. + */ +static INLINE void pressure_floor_init(struct pressure_floor_properties *props, + const struct phys_const *phys_const, + const struct unit_system *us, + const struct hydro_props *hydro_props, + struct swift_params *params) { + + /* Read the Jeans factor */ + props->n_jeans = + parser_get_param_float(params, "GEARPressureFloor:Jeans_factor"); + + /* Compute the constants */ + props->constants = + 4.0 * M_1_PI * phys_const->const_newton_G * pow(props->n_jeans, 2. / 3.); +} + +/** + * @brief Print the properties of the pressure floor to stdout. + * + * @param props The pressure floor properties. + */ +static INLINE void pressure_floor_print( + const struct pressure_floor_properties *props) { + + message("Pressure floor is 'GEAR' with:"); + message("Jeans factor: %g", props->n_jeans); +} + +#ifdef HAVE_HDF5 + +/** + * @brief Writes the current model of pressure floor to the file + * @param h_grp The HDF5 group in which to write + */ +INLINE static void pressure_floor_print_snapshot(hid_t h_grp) { + + io_write_attribute_s(h_grp, "Pressure floor", "GEAR"); +} +#endif +#endif /* SWIFT_PRESSURE_FLOOR_GEAR_H */ diff --git a/src/pressure_floor/none/pressure_floor.h b/src/pressure_floor/none/pressure_floor.h new file mode 100644 index 0000000000000000000000000000000000000000..2d0b95a71f9ccce6697e79760aa43b560933e7bd --- /dev/null +++ b/src/pressure_floor/none/pressure_floor.h @@ -0,0 +1,99 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_PRESSURE_FLOOR_NONE_H +#define SWIFT_PRESSURE_FLOOR_NONE_H + +#include "adiabatic_index.h" +#include "cosmology.h" +#include "equation_of_state.h" +#include "hydro_properties.h" +#include "parser.h" +#include "part.h" +#include "units.h" + +/** + * @file src/pressure_floor/none/pressure_floor.h + * @brief Model without pressure floor + */ + +/** + * @brief Properties of the pressure floor in the NONE model. + */ +struct pressure_floor_properties {}; + +/** + * @brief Compute the physical pressure floor of a given #part. + * + * Note that the particle is not updated!! + * + * The inputs can be either in physical or comoving coordinates (the output is + * in the same coordinates). + * + * @param p The #part. + * @param rho The physical or comoving density. + * @param pressure The physical or comoving pressure without any pressure floor. + * + * @return The physical or comoving pressure with the floor. + */ +static INLINE float pressure_floor_get_pressure(const struct part *p, + const float rho, + const float pressure) { + return pressure; +} + +/** + * @brief Initialise the pressure floor by reading the parameters and converting + * to internal units. + * + * The input temperatures and number densities are converted to pressure and + * density assuming a neutral gas of primoridal abundance. + * + * @param params The YAML parameter file. + * @param us The system of units used internally. + * @param phys_const The physical constants. + * @param hydro_props The propoerties of the hydro scheme. + * @param props The pressure floor properties to fill. + */ +static INLINE void pressure_floor_init(struct pressure_floor_properties *props, + const struct phys_const *phys_const, + const struct unit_system *us, + const struct hydro_props *hydro_props, + struct swift_params *params) {} + +/** + * @brief Print the properties of the pressure floor to stdout. + * + * @param props The pressure floor properties. + */ +static INLINE void pressure_floor_print( + const struct pressure_floor_properties *props) {} + +#ifdef HAVE_HDF5 + +/** + * @brief Writes the current model of pressure floor to the file + * @param h_grp The HDF5 group in which to write + */ +INLINE static void pressure_floor_print_snapshot(hid_t h_grp) { + + io_write_attribute_s(h_grp, "Pressure floor", "none"); +} +#endif + +#endif /* SWIFT_PRESSURE_FLOOR_NONE_H */ diff --git a/src/swift.h b/src/swift.h index b9f8818d8b833231971abb1afb36ee4507648488..7790ec6f9203a6534e2f569017d90090352d8f28 100644 --- a/src/swift.h +++ b/src/swift.h @@ -64,6 +64,7 @@ #include "periodic.h" #include "physical_constants.h" #include "potential.h" +#include "pressure_floor.h" #include "profiler.h" #include "queue.h" #include "random.h"