diff --git a/configure.ac b/configure.ac index a8facd73fc9132081c90a5ecb640c4e877e1bd6a..795c9870356e1f24d7676760864f67ebaf7b6b74 100644 --- a/configure.ac +++ b/configure.ac @@ -2087,8 +2087,8 @@ AC_MSG_RESULT([ Make gravity glass : $gravity_glass_making External potential : $with_potential - Entropy floor : $with_entropy_floor Pressure floor : $with_pressure_floor + Entropy floor : $with_entropy_floor Cooling function : $with_cooling Chemistry : $with_chemistry Tracers : $with_tracers diff --git a/src/Makefile.am b/src/Makefile.am index 3d748b135c93cf409464a199232b243c699fabe1..4710f3c99371f3bd585c328db51426363651e04e 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -50,8 +50,9 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \ chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h space_getsid.h utilities.h \ mesh_gravity.h cbrt.h exp10.h velociraptor_interface.h swift_velociraptor_part.h outputlist.h \ logger_io.h tracers_io.h tracers.h tracers_struct.h star_formation_io.h fof.h fof_struct.h fof_io.h \ - star_formation_struct.h star_formation.h star_formation_iact.h \ + star_formation_struct.h star_formation.h \ star_formation_logger.h star_formation_logger_struct.h \ + pressure_floor.h pressure_floor_struct.h pressure_floor_iact.h \ velociraptor_struct.h velociraptor_io.h random.h memuse.h black_holes.h black_holes_io.h \ black_holes_properties.h black_holes_struct.h feedback.h feedback_struct.h feedback_properties.h @@ -172,11 +173,11 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h potential/isothermal/potential.h potential/disc_patch/potential.h \ potential/sine_wave/potential.h \ star_formation/none/star_formation.h star_formation/none/star_formation_struct.h \ - star_formation/none/star_formation_io.h star_formation/none/star_formation_iact.h \ + star_formation/none/star_formation_io.h \ star_formation/EAGLE/star_formation.h star_formation/EAGLE/star_formation_struct.h \ - star_formation/EAGLE/star_formation_io.h star_formation/EAGLE/star_formation_iact.h \ + star_formation/EAGLE/star_formation_io.h \ star_formation/GEAR/star_formation.h star_formation/GEAR/star_formation_struct.h \ - star_formation/GEAR/star_formation_io.h star_formation/GEAR/star_formation_iact.h \ + star_formation/GEAR/star_formation_io.h \ star_formation/EAGLE/star_formation_logger.h star_formation/EAGLE/star_formation_logger_struct.h \ star_formation/GEAR/star_formation_logger.h star_formation/GEAR/star_formation_logger_struct.h \ star_formation/none/star_formation_logger.h star_formation/none/star_formation_logger_struct.h \ @@ -223,8 +224,9 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h black_holes/EAGLE/black_holes_part.h black_holes/EAGLE/black_holes_iact.h \ black_holes/EAGLE/black_holes_properties.h \ black_holes/EAGLE/black_holes_struct.h \ - pressure_floor.h \ - pressure_floor/GEAR/pressure_floor.h pressure_floor/none/pressure_floor.h + pressure_floor/GEAR/pressure_floor.h pressure_floor/none/pressure_floor.h \ + pressure_floor/GEAR/pressure_floor_iact.h pressure_floor/none/pressure_floor_iact.h \ + pressure_floor/GEAR/pressure_floor_struct.h pressure_floor/none/pressure_floor_struct.h # Sources and flags for regular library diff --git a/src/adiabatic_index.h b/src/adiabatic_index.h index de7c3871cfb6e42739edf45a7d5b4882547d3cc2..d64968d736df0f0539a568632e8cf9c50a85145e 100644 --- a/src/adiabatic_index.h +++ b/src/adiabatic_index.h @@ -545,4 +545,38 @@ pow_three_gamma_minus_five_over_two(float x) { #endif } +/** + * @brief Return the argument to the power three (adiabatic index - 1). + * + * Computes \f$x^{3(\gamma - 1)}\f$. + * + * @param x Argument + */ +__attribute__((always_inline, const)) INLINE static float +pow_three_gamma_minus_one(float x) { + +#if defined(HYDRO_GAMMA_5_3) + + return x * x; /* x^(2) */ + +#elif defined(HYDRO_GAMMA_7_5) + + return powf(x, 1.2f); /* x^(6/5) */ + +#elif defined(HYDRO_GAMMA_4_3) + + return x; /* x^(1) */ + +#elif defined(HYDRO_GAMMA_2_1) + + return x * x * x; /* x^(3) */ + +#else + + error("The adiabatic index is not defined !"); + return 0.f; + +#endif +} + #endif /* SWIFT_ADIABATIC_INDEX_H */ diff --git a/src/cell.c b/src/cell.c index 0649d31b38279f3a75ec2b439f71de105425ae27..32a47e9451076ef18f70b737c68d074d2593026e 100644 --- a/src/cell.c +++ b/src/cell.c @@ -61,6 +61,7 @@ #include "hydro_properties.h" #include "memswap.h" #include "minmax.h" +#include "pressure_floor.h" #include "scheduler.h" #include "space.h" #include "space_getsid.h" @@ -4413,7 +4414,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) { if (part_is_active(p, e)) { hydro_init_part(p, &e->s->hs); chemistry_init_part(p, e->chemistry); - star_formation_init_part(p, xp, e->star_formation); + pressure_floor_init_part(p, xp); tracers_after_init(p, xp, e->internal_units, e->physical_constants, with_cosmology, e->cosmology, e->hydro_properties, e->cooling_func, e->time); diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h index c06b29f4417673bbdac79ab76226770c9c327406..a470eef3fafc32c99fe3a853dcf46051a3086441 100644 --- a/src/chemistry/EAGLE/chemistry.h +++ b/src/chemistry/EAGLE/chemistry.h @@ -243,8 +243,9 @@ static INLINE void chemistry_print_backend( } /** - * @brief Updates the metal mass fractions after diffusion at the end of the - * force loop. + * @brief Updates to the chemistry data after the hydro force loop. + * + * Nothing to do here in EAGLE. * * @param p The particle to act upon. * @param cosmo The current cosmological model. @@ -255,10 +256,13 @@ __attribute__((always_inline)) INLINE static void chemistry_end_force( /** * @brief Computes the chemistry-related time-step constraint. * + * No constraints in the EAGLE model (no diffusion etc.) --> FLT_MAX + * * @param phys_const The physical constants in internal units. * @param cosmo The current cosmological model. * @param us The internal system of units. * @param hydro_props The properties of the hydro scheme. + * @param cd The global properties of the chemistry scheme. * @param p Pointer to the particle data. */ __attribute__((always_inline)) INLINE static float chemistry_timestep( diff --git a/src/chemistry/GEAR/chemistry.h b/src/chemistry/GEAR/chemistry.h index 951d565337eae39bec05c0e142a020e6647811fe..34c4d10e5fb2f491b029c4c8d76cf04ca7a5429a 100644 --- a/src/chemistry/GEAR/chemistry.h +++ b/src/chemistry/GEAR/chemistry.h @@ -135,6 +135,15 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density( } } +/** + * @brief Updates to the chemistry data after the hydro force loop. + * + * @param p The particle to act upon. + * @param cosmo The current cosmological model. + */ +__attribute__((always_inline)) INLINE static void chemistry_end_force( + struct part* restrict p, const struct cosmology* cosmo) {} + /** * @brief Sets all particle fields to sensible values when the #part has 0 ngbs. * @@ -151,6 +160,27 @@ chemistry_part_has_no_neighbours(struct part* restrict p, error("Needs implementing!"); } +/** + * @brief Computes the chemistry-related time-step constraint. + * + * No constraints in the GEAR model (no diffusion) --> FLT_MAX + * + * @param phys_const The physical constants in internal units. + * @param cosmo The current cosmological model. + * @param us The internal system of units. + * @param hydro_props The properties of the hydro scheme. + * @param cd The global properties of the chemistry scheme. + * @param p Pointer to the particle data. + */ +__attribute__((always_inline)) INLINE static float chemistry_timestep( + const struct phys_const* restrict phys_const, + const struct cosmology* restrict cosmo, + const struct unit_system* restrict us, + const struct hydro_props* hydro_props, + const struct chemistry_global_data* cd, const struct part* restrict p) { + return FLT_MAX; +} + /** * @brief Sets the chemistry properties of the (x-)particles to a valid start * state. diff --git a/src/chemistry/none/chemistry.h b/src/chemistry/none/chemistry.h index 543a5e77eea245da9ec18de210c781d5be07d7fb..4e7e790deb86d18a35f2628105da3dd072340747 100644 --- a/src/chemistry/none/chemistry.h +++ b/src/chemistry/none/chemistry.h @@ -87,6 +87,34 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density( struct part* restrict p, const struct chemistry_global_data* cd, const struct cosmology* cosmo) {} +/** + * @brief Updates to the chemistry data after the hydro force loop. + * + * @param p The particle to act upon. + * @param cosmo The current cosmological model. + */ +__attribute__((always_inline)) INLINE static void chemistry_end_force( + struct part* restrict p, const struct cosmology* cosmo) {} + +/** + * @brief Computes the chemistry-related time-step constraint. + * + * @param phys_const The physical constants in internal units. + * @param cosmo The current cosmological model. + * @param us The internal system of units. + * @param hydro_props The properties of the hydro scheme. + * @param cd The global properties of the chemistry scheme. + * @param p Pointer to the particle data. + */ +__attribute__((always_inline)) INLINE static float chemistry_timestep( + const struct phys_const* restrict phys_const, + const struct cosmology* restrict cosmo, + const struct unit_system* restrict us, + const struct hydro_props* hydro_props, + const struct chemistry_global_data* cd, const struct part* restrict p) { + return FLT_MAX; +} + /** * @brief Sets all particle fields to sensible values when the #part has 0 ngbs. * diff --git a/src/hydro/AnarchyDU/hydro.h b/src/hydro/AnarchyDU/hydro.h index ac4ed6c2f8176669aa7364109f95de08ddbf722c..a8b84c46aca9a117efe8173c492bbfd678106a5f 100644 --- a/src/hydro/AnarchyDU/hydro.h +++ b/src/hydro/AnarchyDU/hydro.h @@ -822,9 +822,11 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( * * @param p The particle. * @param xp The extended data of this particle. + * @param cosmo The cosmological model. */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( - struct part *restrict p, const struct xpart *restrict xp) { + struct part *restrict p, const struct xpart *restrict xp, + const struct cosmology *cosmo) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; diff --git a/src/hydro/AnarchyDU/hydro_io.h b/src/hydro/AnarchyDU/hydro_io.h index d78d16c35bf4ebb561116e1a504d24b248104206..4044a0c1f1cdbf4263de03071b046cd4c0884d2e 100644 --- a/src/hydro/AnarchyDU/hydro_io.h +++ b/src/hydro/AnarchyDU/hydro_io.h @@ -182,7 +182,7 @@ INLINE static void hydro_write_particles(const struct part* parts, list[4] = io_make_output_field( "InternalEnergies", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS, - 3. * hydro_gamma_minus_one, parts, u, + -3.f * hydro_gamma_minus_one, parts, u, "Co-moving thermal energies per unit mass of the particles"); list[5] = @@ -198,7 +198,7 @@ INLINE static void hydro_write_particles(const struct part* parts, xparts, convert_S, "Co-moving entropies per unit mass of the particles"); list[8] = io_make_output_field_convert_part( - "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, 3.f * hydro_gamma, parts, + "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, -3.f * hydro_gamma, parts, xparts, convert_P, "Co-moving pressures of the particles"); list[9] = io_make_output_field_convert_part( diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h index 9ba58828dd84df205539550c221746dfac37ed21..806d4aed204c71d3b9265d1bbea1098fdb83ebb1 100644 --- a/src/hydro/AnarchyPU/hydro.h +++ b/src/hydro/AnarchyPU/hydro.h @@ -828,9 +828,11 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( * * @param p The particle. * @param xp The extended data of this particle. + * @param cosmo The cosmological model. */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( - struct part *restrict p, const struct xpart *restrict xp) { + struct part *restrict p, const struct xpart *restrict xp, + const struct cosmology *cosmo) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; diff --git a/src/hydro/AnarchyPU/hydro_io.h b/src/hydro/AnarchyPU/hydro_io.h index 499ef5dc2d79f86f6f76c119b1b421260cf25cea..7dfe61f198163c458ea8ecc22c1cbaa7eb4b99dd 100644 --- a/src/hydro/AnarchyPU/hydro_io.h +++ b/src/hydro/AnarchyPU/hydro_io.h @@ -184,7 +184,7 @@ INLINE static void hydro_write_particles(const struct part* parts, list[4] = io_make_output_field( "InternalEnergies", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS, - 3. * hydro_gamma_minus_one, parts, u, + -3.f * hydro_gamma_minus_one, parts, u, "Co-moving thermal energies per unit mass of the particles"); list[5] = @@ -196,7 +196,7 @@ INLINE static void hydro_write_particles(const struct part* parts, "Co-moving mass densities of the particles"); list[7] = io_make_output_field( - "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, 3.f * hydro_gamma, parts, + "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, -3.f * hydro_gamma, parts, pressure_bar, "Co-moving smoothed pressures of the particles"); list[8] = io_make_output_field_convert_part( diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h index 558143f88d07225f5e6bfe5666d345a7891a7534..469fee137ee624ab649718f25c490f61e35cfb59 100644 --- a/src/hydro/Default/hydro.h +++ b/src/hydro/Default/hydro.h @@ -766,9 +766,11 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( * * @param p The particle. * @param xp The extended data of this particle. + * @param cosmo The cosmological model. */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( - struct part *restrict p, const struct xpart *restrict xp) { + struct part *restrict p, const struct xpart *restrict xp, + const struct cosmology *cosmo) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; diff --git a/src/hydro/Default/hydro_io.h b/src/hydro/Default/hydro_io.h index 5d923837765c3d4759db37cd67badb40d75763f3..1944f677ff16c5f002aa4ee2f830c84e808ffa63 100644 --- a/src/hydro/Default/hydro_io.h +++ b/src/hydro/Default/hydro_io.h @@ -183,7 +183,7 @@ INLINE static void hydro_write_particles(const struct part* parts, list[4] = io_make_output_field( "InternalEnergies", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS, - 3. * hydro_gamma_minus_one, parts, u, + -3.f * hydro_gamma_minus_one, parts, u, "Co-moving thermal energies per unit mass of the particles"); list[5] = @@ -195,7 +195,7 @@ INLINE static void hydro_write_particles(const struct part* parts, "Co-moving mass densities of the particles"); list[7] = io_make_output_field_convert_part( - "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, 3.f * hydro_gamma, parts, + "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, -3.f * hydro_gamma, parts, xparts, convert_P, "Co-moving pressures of the particles"); list[8] = io_make_output_field_convert_part( diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h index 93c0e58068341ab29a7ed0fb6b0d21f76980d952..3fa6f19cc145890feacbd7284368d5378654bf38 100644 --- a/src/hydro/Gadget2/hydro.h +++ b/src/hydro/Gadget2/hydro.h @@ -110,8 +110,7 @@ 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) { - const float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy); - return pressure_floor_get_pressure(p, p->rho, comoving_pressure); + return gas_pressure_from_entropy(p->rho, p->entropy); } /** @@ -123,10 +122,7 @@ __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) { - 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); + return gas_pressure_from_entropy(p->rho * cosmo->a3_inv, p->entropy); } /** @@ -394,7 +390,8 @@ hydro_set_drifted_physical_internal_energy(struct part *p, /* Compute the pressure */ float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy); - comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure); + comoving_pressure = + pressure_floor_get_comoving_pressure(p, comoving_pressure, cosmo); /* Compute the sound speed */ const float soundspeed = @@ -466,7 +463,7 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra( * @brief Prepares a particle for the density calculation. * * Zeroes all the relevant arrays in preparation for the sums taking place in - * the variaous density tasks + * the various density tasks * * @param p The particle to act upon * @param hs #hydro_space containing hydro specific space information. @@ -599,7 +596,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Compute the pressure */ float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy); - comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure); + comoving_pressure = + pressure_floor_get_comoving_pressure(p, comoving_pressure, cosmo); /* Compute the sound speed */ const float soundspeed = @@ -661,9 +659,11 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( * * @param p The particle. * @param xp The extended data of this particle. + * @param cosmo The cosmological model. */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( - struct part *restrict p, const struct xpart *restrict xp) { + struct part *restrict p, const struct xpart *restrict xp, + const struct cosmology *cosmo) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; @@ -675,7 +675,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( /* Re-compute the pressure */ float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy); - comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure); + comoving_pressure = + pressure_floor_get_comoving_pressure(p, comoving_pressure, cosmo); /* Compute the new sound speed */ const float soundspeed = @@ -742,7 +743,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( /* Re-compute the pressure */ float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy); - comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure); + comoving_pressure = + pressure_floor_get_comoving_pressure(p, comoving_pressure, cosmo); /* Compute the new sound speed */ const float soundspeed = @@ -854,7 +856,8 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( /* Compute the pressure */ float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy); - comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure); + comoving_pressure = + pressure_floor_get_comoving_pressure(p, comoving_pressure, cosmo); /* Compute the sound speed */ const float soundspeed = diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h index 9715558be39c853ec5a96262d95cdf4fe92309fa..20ad8e2d0c15094101e44999ff643f9d21619622 100644 --- a/src/hydro/Gadget2/hydro_io.h +++ b/src/hydro/Gadget2/hydro_io.h @@ -177,11 +177,11 @@ INLINE static void hydro_write_particles(const struct part* parts, list[7] = io_make_output_field_convert_part( "InternalEnergies", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS, - 3. * hydro_gamma_minus_one, parts, xparts, convert_part_u, + -3.f * hydro_gamma_minus_one, parts, xparts, convert_part_u, "Co-moving thermal energies per unit mass of the particles"); list[8] = io_make_output_field_convert_part( - "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, 3.f * hydro_gamma, parts, + "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, -3.f * hydro_gamma, parts, xparts, convert_part_P, "Co-moving pressures of the particles"); list[9] = io_make_output_field_convert_part( diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h index 853d2adf17bc069434562fa96ddb881f760f6830..5f4cca366b101c2341863f47d7a9adce51b860a2 100644 --- a/src/hydro/Gadget2/hydro_part.h +++ b/src/hydro/Gadget2/hydro_part.h @@ -35,6 +35,7 @@ #include "chemistry_struct.h" #include "cooling_struct.h" #include "logger.h" +#include "pressure_floor_struct.h" #include "star_formation_struct.h" #include "tracers_struct.h" @@ -155,8 +156,8 @@ struct part { /*! Black holes information (e.g. swallowing ID) */ struct black_holes_part_data black_holes_data; - /* Additional data used by the star formation */ - struct star_formation_part_data sf_data; + /* Additional data used by the pressure floor */ + struct pressure_floor_part_data pressure_floor_data; /* Time-step length */ timebin_t time_bin; diff --git a/src/hydro/GizmoMFM/hydro.h b/src/hydro/GizmoMFM/hydro.h index eb3d9044baecadeaba9165061fae33c816446042..890c7f4c113e289e167247bb4978c1a362b8ff5d 100644 --- a/src/hydro/GizmoMFM/hydro.h +++ b/src/hydro/GizmoMFM/hydro.h @@ -497,9 +497,11 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( * * @param p The particle. * @param xp The extended data of this particle. + * @param cosmo The cosmological model. */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( - struct part* restrict p, const struct xpart* restrict xp) { + struct part* restrict p, const struct xpart* restrict xp, + const struct cosmology* cosmo) { // MATTHIEU: Do we need something here? } diff --git a/src/hydro/GizmoMFM/hydro_io.h b/src/hydro/GizmoMFM/hydro_io.h index b1d51cff90de52285abcd63c519aac911dcba7c2..711eee9e3117c220e6ab5a27a6d8f3557ec7ce13 100644 --- a/src/hydro/GizmoMFM/hydro_io.h +++ b/src/hydro/GizmoMFM/hydro_io.h @@ -220,15 +220,15 @@ INLINE static void hydro_write_particles(const struct part* parts, "Co-moving mass densities of the particles"); list[7] = io_make_output_field_convert_part( - "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, 0.f, parts, xparts, convert_A, + "Entropies", FLOAT, 1, UNIT_CONV_ENTROPY, 0.f, parts, xparts, convert_A, "Co-moving entropies of the particles"); - list[8] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, - 3.f * hydro_gamma, parts, P, + list[8] = io_make_output_field("Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, + -3.f * hydro_gamma, parts, P, "Co-moving pressures of the particles"); list[9] = io_make_output_field_convert_part( - "TotalEnergies", FLOAT, 1, UNIT_CONV_ENERGY, 3.f * hydro_gamma_minus_one, + "TotalEnergies", FLOAT, 1, UNIT_CONV_ENERGY, -3.f * hydro_gamma_minus_one, parts, xparts, convert_Etot, "Total (co-moving) energy of the particles"); list[10] = io_make_output_field_convert_part( diff --git a/src/hydro/GizmoMFV/hydro.h b/src/hydro/GizmoMFV/hydro.h index e6a87b9a86fc7ec03a4bfe055e151ee8c7c580e4..58a8a19dccd2dd102beb2803ec04f1f555cdcec2 100644 --- a/src/hydro/GizmoMFV/hydro.h +++ b/src/hydro/GizmoMFV/hydro.h @@ -539,9 +539,11 @@ __attribute__((always_inline)) INLINE static void hydro_reset_gradient( * * @param p The particle. * @param xp The extended data of this particle. + * @param cosmo The cosmological model. */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( - struct part* restrict p, const struct xpart* restrict xp) { + struct part* restrict p, const struct xpart* restrict xp, + const struct cosmology* cosmo) { // MATTHIEU: Apply the entropy floor here. } diff --git a/src/hydro/GizmoMFV/hydro_io.h b/src/hydro/GizmoMFV/hydro_io.h index 7288df0c510ea48489353b5be4ef9d3f252d5f59..11b8e95867495ccbc9df6a67b3045a09052fb161 100644 --- a/src/hydro/GizmoMFV/hydro_io.h +++ b/src/hydro/GizmoMFV/hydro_io.h @@ -219,15 +219,15 @@ INLINE static void hydro_write_particles(const struct part* parts, "Co-moving mass densities of the particles"); list[7] = io_make_output_field_convert_part( - "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, 0.f, parts, xparts, convert_A, + "Entropies", FLOAT, 1, UNIT_CONV_ENTROPY, 0.f, parts, xparts, convert_A, "Co-moving entropies of the particles"); - list[8] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, - 3.f * hydro_gamma, parts, primitives.P, + list[8] = io_make_output_field("Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, + -3.f * hydro_gamma, parts, primitives.P, "Co-moving pressures of the particles"); list[9] = io_make_output_field_convert_part( - "TotalEnergies", FLOAT, 1, UNIT_CONV_ENERGY, 3.f * hydro_gamma_minus_one, + "TotalEnergies", FLOAT, 1, UNIT_CONV_ENERGY, -3.f * hydro_gamma_minus_one, parts, xparts, convert_Etot, "Total (co-moving) energy of the particles"); list[10] = io_make_output_field_convert_part( diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h index 0cac1de40f258ff9921ea98ddd59a848d4bb7a67..3d7f43579033afd8f5c29e765b31fee145d9c590 100644 --- a/src/hydro/Minimal/hydro.h +++ b/src/hydro/Minimal/hydro.h @@ -642,9 +642,11 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( * * @param p The particle. * @param xp The extended data of this particle. + * @param cosmo The cosmological model */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( - struct part *restrict p, const struct xpart *restrict xp) { + struct part *restrict p, const struct xpart *restrict xp, + const struct cosmology *cosmo) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h index a32e2c1a87c3219640ed94d56725632539121dd8..ba62489bcc8a7bff3b9d3cd45d9b9d0881772194 100644 --- a/src/hydro/Minimal/hydro_io.h +++ b/src/hydro/Minimal/hydro_io.h @@ -179,7 +179,7 @@ INLINE static void hydro_write_particles(const struct part* parts, list[4] = io_make_output_field( "InternalEnergies", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS, - 3. * hydro_gamma_minus_one, parts, u, + -3.f * hydro_gamma_minus_one, parts, u, "Co-moving thermal energies per unit mass of the particles"); list[5] = @@ -195,7 +195,7 @@ INLINE static void hydro_write_particles(const struct part* parts, xparts, convert_S, "Co-moving entropies per unit mass of the particles"); list[8] = io_make_output_field_convert_part( - "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, 3.f * hydro_gamma, parts, + "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, -3.f * hydro_gamma, parts, xparts, convert_P, "Co-moving pressures of the particles"); } diff --git a/src/hydro/Planetary/hydro.h b/src/hydro/Planetary/hydro.h index ee9aa95d5082e4bf21e8ac1ebf6710530638a974..3891bf3df8245b63ec40985910aa7d827b0c00b7 100644 --- a/src/hydro/Planetary/hydro.h +++ b/src/hydro/Planetary/hydro.h @@ -652,9 +652,11 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( * * @param p The particle. * @param xp The extended data of this particle. + * @param cosmo The cosmological model. */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( - struct part *restrict p, const struct xpart *restrict xp) { + struct part *restrict p, const struct xpart *restrict xp, + const struct cosmology *cosmo) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; diff --git a/src/hydro/Planetary/hydro_io.h b/src/hydro/Planetary/hydro_io.h index 6dd84b3e1b00beda160b4b51109b544ac0ad8b86..782e4de6dc05c8d5b485f9f4d2dadd685e44e3a7 100644 --- a/src/hydro/Planetary/hydro_io.h +++ b/src/hydro/Planetary/hydro_io.h @@ -176,7 +176,7 @@ INLINE static void hydro_write_particles(const struct part* parts, "Smoothing lengths (FWHM of the kernel) of the particles"); list[4] = io_make_output_field( "InternalEnergies", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS, - 3. * hydro_gamma_minus_one, parts, u, + -3.f * hydro_gamma_minus_one, parts, u, "Thermal energies per unit mass of the particles"); list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, @@ -190,7 +190,7 @@ INLINE static void hydro_write_particles(const struct part* parts, io_make_output_field("MaterialIDs", INT, 1, UNIT_CONV_NO_UNITS, 0.f, parts, mat_id, "Material IDs of the particles"); list[9] = io_make_output_field_convert_part( - "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, 3.f * hydro_gamma, parts, + "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, -3.f * hydro_gamma, parts, xparts, convert_P, "Pressures of the particles"); list[10] = io_make_output_field_convert_part( "Potentials", FLOAT, 1, UNIT_CONV_POTENTIAL, 0.f, parts, xparts, diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h index 57e0091a7d79fb5404a30c73bca31b2c555e01b1..e49445f59edd7a3a09860bfffe81a56d6a05abd4 100644 --- a/src/hydro/PressureEnergy/hydro.h +++ b/src/hydro/PressureEnergy/hydro.h @@ -213,20 +213,30 @@ hydro_get_drifted_physical_entropy(const struct part *restrict p, } /** - * @brief Returns the comoving sound speed of a particle + * @brief Update the sound speed of a particle * - * @param p The particle of interest + * @param p The particle of interest. + * @param cosmo The cosmological model. */ -__attribute__((always_inline)) INLINE static float -hydro_get_comoving_soundspeed(const struct part *restrict p) { +__attribute__((always_inline)) INLINE static void hydro_update_soundspeed( + struct part *restrict p, const struct cosmology *cosmo) { /* Compute the sound speed -- see theory section for justification */ /* IDEAL GAS ONLY -- P-U does not work with generic EoS. */ 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); + pressure_floor_get_comoving_pressure(p, p->pressure_bar, cosmo); + p->force.soundspeed = gas_soundspeed_from_pressure(p->rho, comoving_pressure); +} + +/** + * @brief Returns the comoving sound speed of a particle + * + * @param p The particle of interest. + */ +__attribute__((always_inline)) INLINE static float +hydro_get_comoving_soundspeed(const struct part *restrict p) { - return square_rooted; + return p->force.soundspeed; } /** @@ -239,10 +249,8 @@ __attribute__((always_inline)) INLINE static float hydro_get_physical_soundspeed(const struct part *restrict p, const struct cosmology *cosmo) { - 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); + /* The pressure floor is already included in p->force.soundspeed */ + return cosmo->a_factor_sound_speed * p->force.soundspeed; } /** @@ -420,11 +428,8 @@ hydro_set_drifted_physical_internal_energy(struct part *p, /* Now recompute the extra quantities */ - /* Compute the sound speed */ - const float soundspeed = hydro_get_comoving_soundspeed(p); - /* Update variables. */ - p->force.soundspeed = soundspeed; + hydro_update_soundspeed(p, cosmo); } /** @@ -629,6 +634,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( const float abs_div_v = fabsf(p->density.div_v); /* Compute the sound speed -- see theory section for justification */ + hydro_update_soundspeed(p, cosmo); const float soundspeed = hydro_get_comoving_soundspeed(p); /* Compute the Balsara switch */ @@ -644,11 +650,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Get the pressures */ const float comoving_pressure_with_floor = - pressure_floor_get_pressure(p, p->rho, p->pressure_bar); + pressure_floor_get_comoving_pressure(p, p->pressure_bar, cosmo); /* 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; } @@ -681,9 +686,11 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( * * @param p The particle. * @param xp The extended data of this particle. + * @param cosmo The cosmological model. */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( - struct part *restrict p, const struct xpart *restrict xp) { + struct part *restrict p, const struct xpart *restrict xp, + const struct cosmology *cosmo) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; @@ -694,9 +701,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( p->u = xp->u_full; /* Compute the sound speed */ - const float soundspeed = hydro_get_comoving_soundspeed(p); - - p->force.soundspeed = soundspeed; + hydro_update_soundspeed(p, cosmo); } /** @@ -759,13 +764,11 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( } /* Compute the new sound speed */ - const float soundspeed = hydro_get_comoving_soundspeed(p); - - p->force.soundspeed = soundspeed; + hydro_update_soundspeed(p, cosmo); /* update the required variables */ const float comoving_pressure_with_floor = - pressure_floor_get_pressure(p, p->rho, p->pressure_bar); + pressure_floor_get_comoving_pressure(p, p->pressure_bar, cosmo); p->force.pressure_bar_with_floor = comoving_pressure_with_floor; } diff --git a/src/hydro/PressureEnergy/hydro_io.h b/src/hydro/PressureEnergy/hydro_io.h index e093fe628e5dc8f546fe9be8133692e95b30af2c..896d0137ca7aae2d17159e4dbf335a9263dce3a7 100644 --- a/src/hydro/PressureEnergy/hydro_io.h +++ b/src/hydro/PressureEnergy/hydro_io.h @@ -175,7 +175,7 @@ INLINE static void hydro_write_particles(const struct part* parts, list[4] = io_make_output_field( "InternalEnergies", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS, - 3. * hydro_gamma_minus_one, parts, u, + -3.f * hydro_gamma_minus_one, parts, u, "Co-moving thermal energies per unit mass of the particles"); list[5] = @@ -187,7 +187,7 @@ INLINE static void hydro_write_particles(const struct part* parts, "Co-moving mass densities of the particles"); list[7] = io_make_output_field( - "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, 3.f * hydro_gamma, parts, + "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, -3.f * hydro_gamma, parts, pressure_bar, "Co-moving smoothed pressures of the particles"); list[8] = io_make_output_field_convert_part( diff --git a/src/hydro/PressureEnergy/hydro_part.h b/src/hydro/PressureEnergy/hydro_part.h index 500569533a5328f3cefb681f5dfbfaf8ab4460c0..6231ad4e1356b91c220dad5e992d68c77975e950 100644 --- a/src/hydro/PressureEnergy/hydro_part.h +++ b/src/hydro/PressureEnergy/hydro_part.h @@ -34,6 +34,7 @@ #include "black_holes_struct.h" #include "chemistry_struct.h" #include "cooling_struct.h" +#include "pressure_floor_struct.h" #include "star_formation_struct.h" #include "tracers_struct.h" @@ -177,6 +178,9 @@ struct part { /*! Black holes information (e.g. swallowing ID) */ struct black_holes_part_data black_holes_data; + /* Additional data used by the pressure floor */ + struct pressure_floor_part_data pressure_floor_data; + /*! Time-step length */ timebin_t time_bin; diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h index c23b3467fd0dee575aacf92c108508be7d44ad32..e21750bdf24f19b0f784b8652bcb5a18a6a89e8f 100644 --- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h +++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h @@ -689,9 +689,11 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( * * @param p The particle. * @param xp The extended data of this particle. + * @param cosmo The cosmological model. */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( - struct part *restrict p, const struct xpart *restrict xp) { + struct part *restrict p, const struct xpart *restrict xp, + const struct cosmology *cosmo) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h index d89dc36ad018f89114a8e80eb2abb663e211d762..7da369847a5d684ebe30dbf7428ea3da9a8c216f 100644 --- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h +++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h @@ -176,7 +176,7 @@ INLINE static void hydro_write_particles(const struct part* parts, list[4] = io_make_output_field( "InternalEnergies", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS, - 3. * hydro_gamma_minus_one, parts, u, + -3.f * hydro_gamma_minus_one, parts, u, "Co-moving thermal energies per unit mass of the particles"); list[5] = @@ -188,7 +188,7 @@ INLINE static void hydro_write_particles(const struct part* parts, "Co-moving mass densities of the particles"); list[7] = io_make_output_field( - "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, 3.f * hydro_gamma, parts, + "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, -3.f * hydro_gamma, parts, pressure_bar, "Co-moving smoothed pressures of the particles"); list[8] = io_make_output_field_convert_part( @@ -196,12 +196,12 @@ INLINE static void hydro_write_particles(const struct part* parts, xparts, convert_S, "Co-moving entropies per unit mass of the particles"); list[9] = io_make_output_field_convert_part( - "Potentials", FLOAT, 1, UNIT_CONV_POTENTIAL, parts, xparts, + "Potentials", FLOAT, 1, UNIT_CONV_POTENTIAL, -1.f, parts, xparts, convert_part_potential, "Gravitational potentials of the particles"); - list[10] = io_make_output_field_convert_part( - "ViscosityParameters", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, parts, xparts, - convert_viscosity, "Visosity coefficient (alpha_visc) of the particles"); + list[10] = io_make_output_field( + "ViscosityParameters", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, parts, alpha, + "Visosity coefficient (alpha_visc) of the particles"); } /** diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h index 66b8e186c82ab1646ffb862ff2d401f0362886c3..83234a7ea145c426767cf4c6d36dc1af8ae490e4 100644 --- a/src/hydro/PressureEntropy/hydro.h +++ b/src/hydro/PressureEntropy/hydro.h @@ -628,9 +628,11 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( * * @param p The particle. * @param xp The extended data of this particle. + * @param cosmo The cosmological model. */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( - struct part *restrict p, const struct xpart *restrict xp) { + struct part *restrict p, const struct xpart *restrict xp, + const struct cosmology *cosmo) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; diff --git a/src/hydro/PressureEntropy/hydro_io.h b/src/hydro/PressureEntropy/hydro_io.h index d7591fbdea1576cf6fec645c7a9842f2197b60a3..8ac77f31efde3b708cbf631c792404b135e4c112 100644 --- a/src/hydro/PressureEntropy/hydro_io.h +++ b/src/hydro/PressureEntropy/hydro_io.h @@ -188,11 +188,11 @@ INLINE static void hydro_write_particles(const struct part* parts, list[7] = io_make_output_field_convert_part( "InternalEnergies", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS, - 3. * hydro_gamma_minus_one, parts, xparts, convert_u, + -3.f * hydro_gamma_minus_one, parts, xparts, convert_u, "Co-moving thermal energies per unit mass of the particles"); list[8] = io_make_output_field_convert_part( - "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, 3.f * hydro_gamma, parts, + "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, -3.f * hydro_gamma, parts, xparts, convert_P, "Co-moving smoothed pressures of the particles"); list[9] = io_make_output_field( diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h index 955a3db1d4480c43dd2fd57c6d976689d36a5652..8c0454025180598baaeafb4032699cddb44b26d1 100644 --- a/src/hydro/Shadowswift/hydro.h +++ b/src/hydro/Shadowswift/hydro.h @@ -397,9 +397,11 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( * * @param p The particle. * @param xp The extended data of this particle. + * @param cosmo The cosmological model. */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( - struct part* restrict p, const struct xpart* restrict xp) {} + struct part* restrict p, const struct xpart* restrict xp, + const struct cosmology* cosmo) {} /** * @brief Converts the hydrodynamic variables from the initial condition file to @@ -797,6 +799,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_density( */ __attribute__((always_inline)) INLINE static float hydro_get_physical_internal_energy(const struct part* restrict p, + const struct xpart* restrict xp, const struct cosmology* cosmo) { return cosmo->a_factor_internal_energy * @@ -874,6 +877,38 @@ hydro_set_drifted_physical_internal_energy(struct part* p, error("Need implementing"); } +/** + * @brief Gets the drifted physical internal energy of a particle + * + * @param p The particle of interest. + * @param cosmo Cosmology data structure + * + * @return The physical internal energy + */ +__attribute__((always_inline)) INLINE static float +hydro_get_drifted_physical_internal_energy(const struct part* p, + const struct cosmology* cosmo) { + error("Need implementing"); + + return 0; +} + +/** + * @brief Gets the drifted physical entropy of a particle + * + * @param p The particle of interest. + * @param cosmo Cosmology data structure + * + * @return The physical entropy + */ +__attribute__((always_inline)) INLINE static float +hydro_get_drifted_physical_entropy(const struct part* p, + const struct cosmology* cosmo) { + error("Need implementing"); + + return 0; +} + /** * @brief Update the value of the viscosity alpha for the scheme. * diff --git a/src/hydro/Shadowswift/hydro_io.h b/src/hydro/Shadowswift/hydro_io.h index 410422331469f980b398139f18ef09d995e0f655..31840a05ee436d5b535ae022bd79118851ef4e5b 100644 --- a/src/hydro/Shadowswift/hydro_io.h +++ b/src/hydro/Shadowswift/hydro_io.h @@ -142,39 +142,42 @@ INLINE static void hydro_write_particles(const struct part* parts, "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f, parts, xparts, convert_part_pos, "Co-moving positions of the particles"); - list[1] = io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, - primitives.v); + list[1] = io_make_output_field( + "Velocities", FLOAT, 3, UNIT_CONV_SPEED, 0.f, parts, primitives.v, + "Peculiar velocities of the stars. This is (a * dx/dt) where x is the " + "co-moving positions of the particles"); - list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, - conserved.mass); + list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f, parts, + conserved.mass, "Masses of the particles"); list[3] = io_make_output_field( "SmoothingLengths", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, parts, h, "Co-moving smoothing lengths (FWHM of the kernel) of the particles"); - list[4] = io_make_output_field_convert_part("InternalEnergies", FLOAT, 1, - UNIT_CONV_ENERGY_PER_UNIT_MASS, - parts, xparts, convert_u); + list[4] = io_make_output_field_convert_part( + "InternalEnergies", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS, + -3.f * hydro_gamma_minus_one, parts, xparts, convert_u, + "Co-moving thermal energies per unit mass of the particles"); list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, parts, id, "Unique IDs of the particles"); list[6] = io_make_output_field("Accelerations", FLOAT, 3, - UNIT_CONV_ACCELERATION, parts, 1.f, a_hydro, + UNIT_CONV_ACCELERATION, 1.f, parts, a_hydro, "Accelerations of the particles(does not " "work in non-cosmological runs)."); - list[7] = io_make_output_field("Densities", FLOAT, 1, UNIT_CONV_DENSITY, - parts, 3.f * hydro_gamma, primitives.rho, + list[7] = io_make_output_field("Densities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f, + parts, primitives.rho, "Co-moving mass densities of the particles"); - list[8] = io_make_output_field("Volumes", FLOAT, 1, UNIT_CONV_VOLUME, parts, - 3.f * hydro_gamma, cell.volume, - "Co-moving volumes of the particles"); + list[8] = + io_make_output_field("Volumes", FLOAT, 1, UNIT_CONV_VOLUME, -3.f, parts, + cell.volume, "Co-moving volumes of the particles"); list[9] = io_make_output_field("GradDensities", FLOAT, 3, UNIT_CONV_DENSITY, - parts, 1.f, primitives.gradients.rho, + 1.f, parts, primitives.gradients.rho, "Gradient densities of the particles"); list[10] = io_make_output_field_convert_part( @@ -182,11 +185,11 @@ INLINE static void hydro_write_particles(const struct part* parts, "Co-moving entropies of the particles"); list[11] = io_make_output_field("Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, - 3.f * hydro_gamma, parts, primitives.P, + -3.f * hydro_gamma, parts, primitives.P, "Co-moving pressures of the particles"); list[12] = io_make_output_field_convert_part( - "TotalEnergies", FLOAT, 1, UNIT_CONV_ENERGY, 3.f * hydro_gamma_minus_one, + "TotalEnergies", FLOAT, 1, UNIT_CONV_ENERGY, -3.f * hydro_gamma_minus_one, parts, xparts, convert_Etot, "Total (co-moving) energy of the particles"); } diff --git a/src/pressure_floor/GEAR/pressure_floor.h b/src/pressure_floor/GEAR/pressure_floor.h index dd715c155a095f9ad5f97b1d14cabfc94d9b11b0..de9ad6cf4450005becade48022be1be481fc1cc9 100644 --- a/src/pressure_floor/GEAR/pressure_floor.h +++ b/src/pressure_floor/GEAR/pressure_floor.h @@ -1,6 +1,6 @@ /******************************************************************************* * This file is part of SWIFT. - * Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl) + * 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 @@ -19,9 +19,19 @@ #ifndef SWIFT_PRESSURE_FLOOR_GEAR_H #define SWIFT_PRESSURE_FLOOR_GEAR_H +/* Forward declaration */ +__attribute__((always_inline)) static INLINE float +pressure_floor_get_comoving_pressure(const struct part* p, const float pressure, + const struct cosmology* cosmo); +__attribute__((always_inline)) static INLINE float +pressure_floor_get_physical_pressure(const struct part* p, const float pressure, + const struct cosmology* cosmo); + #include "adiabatic_index.h" #include "cosmology.h" +#include "dimension.h" #include "equation_of_state.h" +#include "hydro.h" #include "hydro_properties.h" #include "parser.h" #include "part.h" @@ -49,25 +59,54 @@ struct pressure_floor_properties { * * 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 pressure_physical The physical pressure without any pressure floor. + * @param cosmo The #cosmology model. + * + * @return The physical pressure with the floor. + */ +__attribute__((always_inline)) static INLINE float +pressure_floor_get_physical_pressure(const struct part* p, + const float pressure_physical, + const struct cosmology* cosmo) { + + const float H_phys = p->h * cosmo->a_inv * kernel_gamma; + const float rho = hydro_get_physical_density(p, cosmo); + + /* Compute the pressure floor */ + float floor = H_phys * H_phys * rho * pressure_floor_props.constants - + p->pressure_floor_data.sigma2; + floor *= rho * hydro_one_over_gamma; + + return fmaxf(pressure_physical, floor); +} + +/** + * @brief Compute the comoving pressure floor of a given #part. + * + * Note that the particle is not updated!! * * @param p The #part. - * @param rho The physical or comoving density. - * @param pressure The physical or comoving pressure without any pressure floor. + * @param pressure_comoving The comoving pressure without any pressure floor. + * @param cosmo The #cosmology model. * * @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) { +__attribute__((always_inline)) static INLINE float +pressure_floor_get_comoving_pressure(const struct part* p, + const float pressure_comoving, + const struct cosmology* cosmo) { - /* 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; + const float a_coef = pow_three_gamma_minus_one(cosmo->a); + const float rho = hydro_get_comoving_density(p); + + /* Compute the pressure floor */ + float floor = kernel_gamma * kernel_gamma * p->h * p->h * rho * + pressure_floor_props.constants; + floor -= p->pressure_floor_data.sigma2 * cosmo->a * cosmo->a; + floor *= a_coef * rho * hydro_one_over_gamma; - return fmax(pressure, floor); + return fmaxf(pressure_comoving, floor); } /** @@ -83,11 +122,10 @@ static INLINE float pressure_floor_get_pressure(const struct part *p, * @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) { +__attribute__((always_inline)) 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 = @@ -103,8 +141,8 @@ static INLINE void pressure_floor_init(struct pressure_floor_properties *props, * * @param props The pressure floor properties. */ -static INLINE void pressure_floor_print( - const struct pressure_floor_properties *props) { +__attribute__((always_inline)) 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); @@ -116,9 +154,76 @@ static INLINE void pressure_floor_print( * @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) { +__attribute__((always_inline)) INLINE static void pressure_floor_print_snapshot( + hid_t h_grp) { io_write_attribute_s(h_grp, "Pressure floor", "GEAR"); } + +/** + * @brief Finishes the density calculation. + * + * @param p The particle to act upon + * @param cosmo The current cosmological model. + */ +__attribute__((always_inline)) INLINE static void pressure_floor_end_density( + struct part* restrict p, const struct cosmology* cosmo) { + + /* To finish the turbulence estimation we devide by the density */ + p->pressure_floor_data.sigma2 /= + pow_dimension(p->h) * hydro_get_comoving_density(p); + + /* Add the cosmological term */ + p->pressure_floor_data.sigma2 *= cosmo->a2_inv; +} + +/** + * @brief Sets all particle fields to sensible values when the #part has 0 ngbs. + * + * @param p The particle to act upon + * @param xp The extended particle data to act upon + * @param cosmo The current cosmological model. + */ +__attribute__((always_inline)) INLINE static void +pressure_floor_part_has_no_neighbours(struct part* restrict p, + struct xpart* restrict xp, + const struct cosmology* cosmo) { + + /* If part has 0 neighbours, the estimation of turbulence is 0 */ + p->pressure_floor_data.sigma2 = 0.f; +} + +/** + * @brief Sets the pressure_floor properties of the (x-)particles to a valid + * start state. + * + * @param p Pointer to the particle data. + * @param xp Pointer to the extended particle data. + */ +__attribute__((always_inline)) INLINE static void pressure_floor_init_part( + struct part* restrict p, struct xpart* restrict xp) { + p->pressure_floor_data.sigma2 = 0.f; +} + +/** + * @brief Sets the pressure_floor properties of the (x-)particles to a valid + * start state. + * + * @param phys_const The physical constant in internal units. + * @param us The unit system. + * @param cosmo The current cosmological model. + * @param p Pointer to the particle data. + * @param xp Pointer to the extended particle data. + */ +__attribute__((always_inline)) INLINE static void +pressure_floor_first_init_part(const struct phys_const* restrict phys_const, + const struct unit_system* restrict us, + const struct cosmology* restrict cosmo, + struct part* restrict p, + struct xpart* restrict xp) { + + pressure_floor_init_part(p, xp); +} + #endif #endif /* SWIFT_PRESSURE_FLOOR_GEAR_H */ diff --git a/src/star_formation/GEAR/star_formation_iact.h b/src/pressure_floor/GEAR/pressure_floor_iact.h similarity index 72% rename from src/star_formation/GEAR/star_formation_iact.h rename to src/pressure_floor/GEAR/pressure_floor_iact.h index 7325b92af2840b317cf1a924a1e509b34bdffba3..5ffb0b0097bb1d117e024f0e01e7babe1f838c40 100644 --- a/src/star_formation/GEAR/star_formation_iact.h +++ b/src/pressure_floor/GEAR/pressure_floor_iact.h @@ -17,16 +17,16 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. * ******************************************************************************/ -#ifndef SWIFT_GEAR_STAR_FORMATION_IACT_H -#define SWIFT_GEAR_STAR_FORMATION_IACT_H +#ifndef SWIFT_GEAR_PRESSURE_FLOOR_IACT_H +#define SWIFT_GEAR_PRESSURE_FLOOR_IACT_H /** - * @file GEAR/star_formation_iact.h + * @file GEAR/pressure_floor_iact.h * @brief Density computation */ /** - * @brief do star_formation computation after the runner_iact_density (symmetric + * @brief do pressure_floor computation after the runner_iact_density (symmetric * version) * * Compute the velocity dispersion follow eq. 2 in Revaz & Jablonka 2018. @@ -40,7 +40,7 @@ * @param a Current scale factor. * @param H Current Hubble parameter. */ -__attribute__((always_inline)) INLINE static void runner_iact_star_formation( +__attribute__((always_inline)) INLINE static void runner_iact_pressure_floor( float r2, const float *dx, float hi, float hj, struct part *restrict pi, struct part *restrict pj, float a, float H) { @@ -53,20 +53,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_star_formation( /* Delta v */ float dv[3] = {pi->v[0] - pj->v[0], pi->v[1] - pj->v[1], pi->v[2] - pj->v[2]}; - /* Norms at power 2 */ - const float norm_v2 = dv[0] * dv[0] + dv[1] * dv[1] + dv[2] * dv[2]; - const float norm_x2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; - /* Compute the velocity dispersion */ - const float sigma2 = norm_v2 + H * norm_x2; + const float a2H = a * a * H; + const float sigma[3] = {dv[0] + a2H * dx[0], dv[1] + a2H * dx[1], + dv[2] + a2H * dx[2]}; + const float sigma2 = + sigma[0] * sigma[0] + sigma[1] * sigma[1] + sigma[2] * sigma[2]; /* Compute the velocity dispersion */ - pi->sf_data.sigma2 += sigma2 * wi * hydro_get_mass(pj); - pj->sf_data.sigma2 += sigma2 * wj * hydro_get_mass(pi); + pi->pressure_floor_data.sigma2 += sigma2 * wi * hydro_get_mass(pj); + pj->pressure_floor_data.sigma2 += sigma2 * wj * hydro_get_mass(pi); } /** - * @brief do star_formation computation after the runner_iact_density (non + * @brief do pressure_floor computation after the runner_iact_density (non * symmetric version) * * @param r2 Comoving square distance between the two particles. @@ -79,7 +79,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_star_formation( * @param H Current Hubble parameter. */ __attribute__((always_inline)) INLINE static void -runner_iact_nonsym_star_formation(float r2, const float *dx, float hi, float hj, +runner_iact_nonsym_pressure_floor(float r2, const float *dx, float hi, float hj, struct part *restrict pi, const struct part *restrict pj, float a, float H) { @@ -90,15 +90,15 @@ runner_iact_nonsym_star_formation(float r2, const float *dx, float hi, float hj, /* Delta v */ float dv[3] = {pi->v[0] - pj->v[0], pi->v[1] - pj->v[1], pi->v[2] - pj->v[2]}; - /* Norms at power 2 */ - const float norm_v2 = dv[0] * dv[0] + dv[1] * dv[1] + dv[2] * dv[2]; - const float norm_x2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; - /* Compute the velocity dispersion */ - const float sigma2 = norm_v2 + H * norm_x2; + const float a2H = a * a * H; + const float sigma[3] = {dv[0] + a2H * dx[0], dv[1] + a2H * dx[1], + dv[2] + a2H * dx[2]}; + const float sigma2 = + sigma[0] * sigma[0] + sigma[1] * sigma[1] + sigma[2] * sigma[2]; /* Compute the velocity dispersion */ - pi->sf_data.sigma2 += sigma2 * wi * hydro_get_mass(pj); + pi->pressure_floor_data.sigma2 += sigma2 * wi * hydro_get_mass(pj); } -#endif /* SWIFT_GEAR_STAR_FORMATION_IACT_H */ +#endif /* SWIFT_GEAR_PRESSURE_FLOOR_IACT_H */ diff --git a/src/pressure_floor/GEAR/pressure_floor_struct.h b/src/pressure_floor/GEAR/pressure_floor_struct.h new file mode 100644 index 0000000000000000000000000000000000000000..1eb70b86dcfb79ce9818e1d7c598d49d5e654efd --- /dev/null +++ b/src/pressure_floor/GEAR/pressure_floor_struct.h @@ -0,0 +1,32 @@ +/******************************************************************************* + * 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_PART_GEAR_H +#define SWIFT_PRESSURE_FLOOR_PART_GEAR_H + +/** + * Structure containing the required variables for the pressure + * floor in the density loop. + */ +struct pressure_floor_part_data { + /*! Estimation of local turbulence (squared) + * Units: length^2 / time^2 (physical) */ + float sigma2; +}; + +#endif // SWIFT_PRESSURE_FLOOR_PART_GEAR_H diff --git a/src/pressure_floor/none/pressure_floor.h b/src/pressure_floor/none/pressure_floor.h index 2d0b95a71f9ccce6697e79760aa43b560933e7bd..84db9cbaf682b36fa68abc7d26273e15cd1da191 100644 --- a/src/pressure_floor/none/pressure_floor.h +++ b/src/pressure_floor/none/pressure_floor.h @@ -42,19 +42,33 @@ struct pressure_floor_properties {}; * * 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 physical_pressure The physical pressure without any pressure floor. + * @param cosmo The #cosmology model. + * + * @return The physical pressure with the floor. + */ +static INLINE float pressure_floor_get_physical_pressure( + const struct part* p, const float physical_pressure, + const struct cosmology* cosmo) { + return physical_pressure; +} + +/** + * @brief Compute the comoving pressure floor of a given #part. + * + * Note that the particle is not updated!! * * @param p The #part. - * @param rho The physical or comoving density. - * @param pressure The physical or comoving pressure without any pressure floor. + * @param comoving_pressure The comoving pressure without any pressure floor. + * @param cosmo The #cosmology model. * - * @return The physical or comoving pressure with the floor. + * @return The comoving pressure with the floor. */ -static INLINE float pressure_floor_get_pressure(const struct part *p, - const float rho, - const float pressure) { - return pressure; +static INLINE float pressure_floor_get_comoving_pressure( + const struct part* p, const float comoving_pressure, + const struct cosmology* cosmo) { + return comoving_pressure; } /** @@ -70,11 +84,11 @@ static INLINE float pressure_floor_get_pressure(const struct part *p, * @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) {} +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. @@ -82,7 +96,7 @@ static INLINE void pressure_floor_init(struct pressure_floor_properties *props, * @param props The pressure floor properties. */ static INLINE void pressure_floor_print( - const struct pressure_floor_properties *props) {} + const struct pressure_floor_properties* props) {} #ifdef HAVE_HDF5 @@ -96,4 +110,52 @@ INLINE static void pressure_floor_print_snapshot(hid_t h_grp) { } #endif +/** + * @brief Finishes the density calculation. + * + * @param p The particle to act upon + * @param cosmo The current cosmological model. + */ +__attribute__((always_inline)) INLINE static void pressure_floor_end_density( + struct part* restrict p, const struct cosmology* cosmo) {} + +/** + * @brief Sets all particle fields to sensible values when the #part has 0 ngbs. + * + * @param p The particle to act upon + * @param xp The extended particle data to act upon + * @param cosmo The current cosmological model. + */ +__attribute__((always_inline)) INLINE static void +pressure_floor_part_has_no_neighbours(struct part* restrict p, + struct xpart* restrict xp, + const struct cosmology* cosmo) {} + +/** + * @brief Sets the pressure_floor properties of the (x-)particles to a valid + * start state. + * + * @param p Pointer to the particle data. + * @param xp Pointer to the extended particle data. + */ +__attribute__((always_inline)) INLINE static void pressure_floor_init_part( + struct part* restrict p, struct xpart* restrict xp) {} + +/** + * @brief Sets the pressure_floor properties of the (x-)particles to a valid + * start state. + * + * @param phys_const The physical constant in internal units. + * @param us The unit system. + * @param cosmo The current cosmological model. + * @param p Pointer to the particle data. + * @param xp Pointer to the extended particle data. + */ +__attribute__((always_inline)) INLINE static void +pressure_floor_first_init_part(const struct phys_const* restrict phys_const, + const struct unit_system* restrict us, + const struct cosmology* restrict cosmo, + struct part* restrict p, + struct xpart* restrict xp) {} + #endif /* SWIFT_PRESSURE_FLOOR_NONE_H */ diff --git a/src/star_formation/none/star_formation_iact.h b/src/pressure_floor/none/pressure_floor_iact.h similarity index 75% rename from src/star_formation/none/star_formation_iact.h rename to src/pressure_floor/none/pressure_floor_iact.h index dd74115bec699029748806b512c9d6bd7fb829fe..fedacfe06be1a6de2a87f14ce2b5a6306f29ac09 100644 --- a/src/star_formation/none/star_formation_iact.h +++ b/src/pressure_floor/none/pressure_floor_iact.h @@ -1,4 +1,7 @@ /******************************************************************************* + * This file is part of SWIFT. + * Coypright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) + * 2019 Fabien Jeanquartier (fabien.jeanquartier@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 @@ -14,18 +17,20 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. * ******************************************************************************/ -#ifndef SWIFT_NONE_STAR_FORMATION_IACT_H -#define SWIFT_NONE_STAR_FORMATION_IACT_H +#ifndef SWIFT_NONE_PRESSURE_FLOOR_IACT_H +#define SWIFT_NONE_PRESSURE_FLOOR_IACT_H /** - * @file none/star_formation_iact.h + * @file NONE/pressure_floor_iact.h * @brief Density computation */ /** - * @brief do star_formation computation after the runner_iact_density (symmetric + * @brief do pressure_floor computation after the runner_iact_density (symmetric * version) * + * Compute the velocity dispersion follow eq. 2 in Revaz & Jablonka 2018. + * * @param r2 Comoving square distance between the two particles. * @param dx Comoving vector separating both particles (pi - pj). * @param hi Comoving smoothing-length of particle i. @@ -35,12 +40,12 @@ * @param a Current scale factor. * @param H Current Hubble parameter. */ -__attribute__((always_inline)) INLINE static void runner_iact_star_formation( +__attribute__((always_inline)) INLINE static void runner_iact_pressure_floor( float r2, const float *dx, float hi, float hj, struct part *restrict pi, struct part *restrict pj, float a, float H) {} /** - * @brief do star_formation computation after the runner_iact_density (non + * @brief do pressure_floor computation after the runner_iact_density (non * symmetric version) * * @param r2 Comoving square distance between the two particles. @@ -53,9 +58,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_star_formation( * @param H Current Hubble parameter. */ __attribute__((always_inline)) INLINE static void -runner_iact_nonsym_star_formation(float r2, const float *dx, float hi, float hj, +runner_iact_nonsym_pressure_floor(float r2, const float *dx, float hi, float hj, struct part *restrict pi, const struct part *restrict pj, float a, float H) {} -#endif /* SWIFT_NONE_STAR_FORMATION_IACT_H */ +#endif /* SWIFT_NONE_PRESSURE_FLOOR_IACT_H */ diff --git a/src/pressure_floor/none/pressure_floor_struct.h b/src/pressure_floor/none/pressure_floor_struct.h new file mode 100644 index 0000000000000000000000000000000000000000..d2eee232faff5c74f2f4af6597d919a7901c3473 --- /dev/null +++ b/src/pressure_floor/none/pressure_floor_struct.h @@ -0,0 +1,28 @@ +/******************************************************************************* + * 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_PART_NONE_H +#define SWIFT_PRESSURE_FLOOR_PART_NONE_H + +/** + * Structure containing the required variables for the pressure + * floor in the density loop. + */ +struct pressure_floor_part_data {}; + +#endif // SWIFT_PRESSURE_FLOOR_PART_NONE_H diff --git a/src/star_formation_iact.h b/src/pressure_floor_iact.h similarity index 59% rename from src/star_formation_iact.h rename to src/pressure_floor_iact.h index ef457214a23102bc33385705db41c89dc29d8b8f..258bf25a5a8429668c30b1d78532159ff32916b9 100644 --- a/src/star_formation_iact.h +++ b/src/pressure_floor_iact.h @@ -1,6 +1,6 @@ /******************************************************************************* * This file is part of SWIFT. - * Copyright (c) 2018 Folkert Nobels (nobels@strw.leidenuniv.nl) + * 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 @@ -16,26 +16,24 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. * ******************************************************************************/ -#ifndef SWIFT_STAR_FORMATION_IACT_H -#define SWIFT_STAR_FORMATION_IACT_H +#ifndef SWIFT_PRESSURE_FLOOR_IACT_H +#define SWIFT_PRESSURE_FLOOR_IACT_H /** - * @file src/star_formation_iact.h - * @brief Branches between the different star formation iact. + * @file src/pressure_floor_iact.h + * @brief Branches between the different pressure floor iact. */ /* Config parameters. */ #include "../config.h" -/* Import the right star formation law definition */ -#if defined(STAR_FORMATION_NONE) -#include "./star_formation/none/star_formation_iact.h" -#elif defined(STAR_FORMATION_EAGLE) -#include "./star_formation/EAGLE/star_formation_iact.h" -#elif defined(STAR_FORMATION_GEAR) -#include "./star_formation/GEAR/star_formation_iact.h" +/* Import the right pressure floor definition */ +#if defined(PRESSURE_FLOOR_NONE) +#include "./pressure_floor/none/pressure_floor_iact.h" +#elif defined(PRESSURE_FLOOR_GEAR) +#include "./pressure_floor/GEAR/pressure_floor_iact.h" #else -#error "Invalid choice of star formation law" +#error "Invalid choice of pressure floor" #endif -#endif /* SWIFT_STAR_FORMATION_IACT_H */ +#endif /* SWIFT_PRESSURE_FLOOR_IACT_H */ diff --git a/src/pressure_floor_struct.h b/src/pressure_floor_struct.h new file mode 100644 index 0000000000000000000000000000000000000000..e5538b37e03258fab45f9a2cb067f571ff7d16fd --- /dev/null +++ b/src/pressure_floor_struct.h @@ -0,0 +1,39 @@ +/******************************************************************************* + * 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/>. + * + ******************************************************************************/ +#ifndef SWIFT_PRESSURE_FLOOR_STRUCT_H +#define SWIFT_PRESSURE_FLOOR_STRUCT_H + +/** + * @file src/pressure_floor_struct.h + * @brief Branches between the different pressure floor data + */ + +/* Config parameters. */ +#include "../config.h" + +/* Import the right pressure floor definition */ +#if defined(PRESSURE_FLOOR_NONE) +#include "./pressure_floor/none/pressure_floor_struct.h" +#elif defined(PRESSURE_FLOOR_GEAR) +#include "./pressure_floor/GEAR/pressure_floor_struct.h" +#else +#error "Invalid choice of pressure floor structure." +#endif + +#endif /* SWIFT_PRESSURE_FLOOR_STRUCT_H */ diff --git a/src/proxy.c b/src/proxy.c index ee2ba818541fdb874fff9c865ac4bb2ee02b371c..4e7e979a68c311ecdde7d36f214a6d5dcded4f5e 100644 --- a/src/proxy.c +++ b/src/proxy.c @@ -69,7 +69,10 @@ void proxy_tags_exchange(struct proxy *proxies, int num_proxies, /* Run through the cells and get the size of the tags that will be sent off. */ int count_out = 0; - int offset_out[s->nr_cells]; + int *offset_out = + (int *)swift_malloc("tags_offsets_out", s->nr_cells * sizeof(int)); + if (offset_out == NULL) error("Error allocating memory for tag offsets"); + for (int k = 0; k < s->nr_cells; k++) { offset_out[k] = count_out; if (s->cells_top[k].mpi.sendto) { @@ -79,7 +82,10 @@ void proxy_tags_exchange(struct proxy *proxies, int num_proxies, /* Run through the proxies and get the count of incoming tags. */ int count_in = 0; - int offset_in[s->nr_cells]; + int *offset_in = + (int *)swift_malloc("tags_offsets_in", s->nr_cells * sizeof(int)); + if (offset_in == NULL) error("Error allocating memory for tag offsets"); + for (int k = 0; k < num_proxies; k++) { for (int j = 0; j < proxies[k].nr_cells_in; j++) { offset_in[proxies[k].cells_in[j] - s->cells_top] = count_in; @@ -170,6 +176,8 @@ void proxy_tags_exchange(struct proxy *proxies, int num_proxies, /* Clean up. */ swift_free("tags_in", tags_in); swift_free("tags_out", tags_out); + swift_free("tags_offsets_in", offset_in); + swift_free("tags_offsets_out", offset_out); free(reqs_in); free(cids_in); @@ -389,7 +397,10 @@ void proxy_cells_exchange(struct proxy *proxies, int num_proxies, s->nr_cells, sizeof(struct cell), /*chunk=*/0, /*extra_data=*/NULL); int count_out = 0; - int offset[s->nr_cells]; + int *offset = + (int *)swift_malloc("proxy_cell_offset", s->nr_cells * sizeof(int)); + if (offset == NULL) error("Error allocating memory for proxy cell offsets"); + for (int k = 0; k < s->nr_cells; k++) { offset[k] = count_out; if (s->cells_top[k].mpi.sendto) count_out += s->cells_top[k].mpi.pcell_size; @@ -472,6 +483,7 @@ void proxy_cells_exchange(struct proxy *proxies, int num_proxies, /* Clean up. */ free(reqs); swift_free("pcells", pcells); + swift_free("proxy_cell_offset", offset); #else error("SWIFT was not compiled with MPI support."); diff --git a/src/runner.c b/src/runner.c index 05dc33e97e74901bd87d1fe4da4259b6d8aa29e0..eb63cca223886d594c8bef9c2072762213a6f495 100644 --- a/src/runner.c +++ b/src/runner.c @@ -60,13 +60,14 @@ #include "logger.h" #include "memuse.h" #include "minmax.h" +#include "pressure_floor.h" +#include "pressure_floor_iact.h" #include "runner_doiact_vec.h" #include "scheduler.h" #include "sort_part.h" #include "space.h" #include "space_getsid.h" #include "star_formation.h" -#include "star_formation_iact.h" #include "star_formation_logger.h" #include "stars.h" #include "task.h" @@ -2032,7 +2033,6 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { const struct hydro_space *hs = &s->hs; const struct cosmology *cosmo = e->cosmology; const struct chemistry_global_data *chemistry = e->chemistry; - const struct star_formation *star_formation = e->star_formation; const int with_cosmology = (e->policy & engine_policy_cosmology); @@ -2129,7 +2129,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { /* Finish the density calculation */ hydro_end_density(p, cosmo); chemistry_end_density(p, chemistry, cosmo); - star_formation_end_density(p, star_formation, cosmo); + pressure_floor_end_density(p, cosmo); /* Compute one step of the Newton-Raphson scheme */ const float n_sum = p->density.wcount * h_old_dim; @@ -2276,7 +2276,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { /* Re-initialise everything */ hydro_init_part(p, hs); chemistry_init_part(p, chemistry); - star_formation_init_part(p, xp, star_formation); + pressure_floor_init_part(p, xp); tracers_after_init(p, xp, e->internal_units, e->physical_constants, with_cosmology, e->cosmology, e->hydro_properties, e->cooling_func, e->time); @@ -2298,8 +2298,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { if (has_no_neighbours) { hydro_part_has_no_neighbours(p, xp, cosmo); chemistry_part_has_no_neighbours(p, xp, chemistry, cosmo); - star_formation_part_has_no_neighbours(p, xp, star_formation, - cosmo); + pressure_floor_part_has_no_neighbours(p, xp, cosmo); } } else { @@ -2962,7 +2961,7 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) { #endif /* Prepare the values to be drifted */ - hydro_reset_predicted_values(p, xp); + hydro_reset_predicted_values(p, xp, cosmo); } } @@ -3618,6 +3617,7 @@ void runner_do_end_hydro_force(struct runner *r, struct cell *c, int timer) { /* Finish the force loop */ hydro_end_force(p, cosmo); + chemistry_end_force(p, cosmo); #ifdef SWIFT_BOUNDARY_PARTICLES diff --git a/src/runner_doiact.h b/src/runner_doiact.h index 1a39c8d49f2c4234c04982e255705f06ec1c5d38..8aabb05d177385c6bbee1a91eb2ea231ccbca3e4 100644 --- a/src/runner_doiact.h +++ b/src/runner_doiact.h @@ -213,7 +213,7 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci, IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); #endif } if (r2 < hjg2 && pj_active) { @@ -225,7 +225,7 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci, IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); #endif } @@ -325,14 +325,14 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci, IACT(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_star_formation(r2, dx, hi, hj, pi, pj, a, H); + runner_iact_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); #endif } else if (pi_active) { IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); #endif } else if (pj_active) { @@ -343,7 +343,7 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci, IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); #endif } } @@ -431,14 +431,14 @@ void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) { IACT(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_star_formation(r2, dx, hi, hj, pi, pj, a, H); + runner_iact_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); #endif } else if (doi) { IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); #endif } else if (doj) { @@ -449,7 +449,7 @@ void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) { IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); #endif } } /* loop over the parts in cj. */ @@ -536,14 +536,14 @@ void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) { IACT(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_star_formation(r2, dx, hi, hj, pi, pj, a, H); + runner_iact_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); #endif } else if (doi) { IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); #endif } else if (doj) { @@ -554,7 +554,7 @@ void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) { IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); #endif } } /* loop over the parts in cj. */ @@ -640,7 +640,7 @@ void DOPAIR_SUBSET_NAIVE(struct runner *r, struct cell *restrict ci, IACT_NONSYM(r2, dx, hi, pj->h, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hi, pj->h, pi, pj, a, H); - runner_iact_nonsym_star_formation(r2, dx, hi, pj->h, pi, pj, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hi, pj->h, pi, pj, a, H); #endif } } /* loop over the parts in cj. */ @@ -733,7 +733,7 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci, IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); #endif } } /* loop over the parts in cj. */ @@ -789,7 +789,7 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci, IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); #endif } } /* loop over the parts in cj. */ @@ -931,7 +931,7 @@ void DOSELF_SUBSET(struct runner *r, struct cell *restrict ci, IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); #endif } } /* loop over the parts in cj. */ @@ -1095,7 +1095,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid, IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); #endif } } /* loop over the parts in cj. */ @@ -1183,7 +1183,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid, IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); #endif } } /* loop over the parts in ci. */ @@ -1483,7 +1483,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid, IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); #endif } } /* loop over the active parts in cj. */ @@ -1554,13 +1554,13 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid, IACT(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_star_formation(r2, dx, hi, hj, pi, pj, a, H); + runner_iact_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); #endif } else { IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); #endif } } @@ -1663,7 +1663,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid, IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); #endif } } /* loop over the active parts in ci. */ @@ -1736,13 +1736,13 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid, IACT(r2, dx, hj, hi, pj, pi, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_chemistry(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_star_formation(r2, dx, hj, hi, pj, pi, a, H); + runner_iact_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); #endif } else { IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); #endif } } @@ -1931,7 +1931,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) { IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); #endif } } /* loop over all other particles. */ @@ -1983,14 +1983,14 @@ void DOSELF1(struct runner *r, struct cell *restrict c) { IACT(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_star_formation(r2, dx, hi, hj, pi, pj, a, H); + runner_iact_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); #endif } else if (doi) { IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); #endif } else if (doj) { @@ -2000,7 +2000,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) { IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); #endif } } @@ -2127,7 +2127,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) { IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); #endif } } /* loop over all other particles. */ @@ -2174,13 +2174,13 @@ void DOSELF2(struct runner *r, struct cell *restrict c) { IACT(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_star_formation(r2, dx, hi, hj, pi, pj, a, H); + runner_iact_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); #endif } else { IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); #endif } } diff --git a/src/space.c b/src/space.c index d8e37fd5a43a151daf8501b9011922983244d99e..322211ca6c313420a95a4d2e39937594b7ef3b96 100644 --- a/src/space.c +++ b/src/space.c @@ -56,9 +56,9 @@ #include "memuse.h" #include "minmax.h" #include "multipole.h" +#include "pressure_floor.h" #include "restart.h" #include "sort_part.h" -#include "star_formation.h" #include "star_formation_logger.h" #include "stars.h" #include "threadpool.h" @@ -4030,7 +4030,6 @@ void space_first_init_parts_mapper(void *restrict map_data, int count, const int with_gravity = e->policy & engine_policy_self_gravity; const struct chemistry_global_data *chemistry = e->chemistry; - const struct star_formation *star_formation = e->star_formation; const struct cooling_function_data *cool_func = e->cooling_func; /* Check that the smoothing lengths are non-zero */ @@ -4081,9 +4080,8 @@ void space_first_init_parts_mapper(void *restrict map_data, int count, /* Also initialise the chemistry */ chemistry_first_init_part(phys_const, us, cosmo, chemistry, &p[k], &xp[k]); - /* Also initialise the star formation */ - star_formation_first_init_part(phys_const, us, cosmo, star_formation, &p[k], - &xp[k]); + /* Also initialise the pressure floor */ + pressure_floor_first_init_part(phys_const, us, cosmo, &p[k], &xp[k]); /* And the cooling */ cooling_first_init_part(phys_const, us, cosmo, cool_func, &p[k], &xp[k]); @@ -4365,7 +4363,7 @@ void space_init_parts_mapper(void *restrict map_data, int count, for (int k = 0; k < count; k++) { hydro_init_part(&parts[k], hs); chemistry_init_part(&parts[k], e->chemistry); - star_formation_init_part(&parts[k], &xparts[k], e->star_formation); + pressure_floor_init_part(&parts[k], &xparts[k]); tracers_after_init(&parts[k], &xparts[k], e->internal_units, e->physical_constants, with_cosmology, e->cosmology, e->hydro_properties, e->cooling_func, e->time); diff --git a/src/star_formation/EAGLE/star_formation.h b/src/star_formation/EAGLE/star_formation.h index bdd6dcb414d2a559e113b54de290fb47b82e38f8..851f493801dc5cb0beee9cd07ea5415a5ad1ccf1 100644 --- a/src/star_formation/EAGLE/star_formation.h +++ b/src/star_formation/EAGLE/star_formation.h @@ -632,71 +632,4 @@ INLINE static void starformation_print_backend( starform->max_gas_density_HpCM3); } -/** - * @brief Finishes the density calculation. - * - * Nothing to do here. We do not need to compute any quantity in the hydro - * density loop for the EAGLE star formation model. - * - * @param p The particle to act upon - * @param cd The global star_formation information. - * @param cosmo The current cosmological model. - */ -__attribute__((always_inline)) INLINE static void star_formation_end_density( - struct part* restrict p, const struct star_formation* cd, - const struct cosmology* cosmo) {} - -/** - * @brief Sets all particle fields to sensible values when the #part has 0 ngbs. - * - * Nothing to do here. We do not need to compute any quantity in the hydro - * density loop for the EAGLE star formation model. - * - * @param p The particle to act upon - * @param xp The extended particle data to act upon - * @param cd #star_formation containing star_formation informations. - * @param cosmo The current cosmological model. - */ -__attribute__((always_inline)) INLINE static void -star_formation_part_has_no_neighbours(struct part* restrict p, - struct xpart* restrict xp, - const struct star_formation* cd, - const struct cosmology* cosmo) {} - -/** - * @brief Sets the star_formation properties of the (x-)particles to a valid - * start state. - * - * Nothing to do here. - * - * @param phys_const The physical constant in internal units. - * @param us The unit system. - * @param cosmo The current cosmological model. - * @param data The global star_formation information used for this run. - * @param p Pointer to the particle data. - * @param xp Pointer to the extended particle data. - */ -__attribute__((always_inline)) INLINE static void -star_formation_first_init_part(const struct phys_const* restrict phys_const, - const struct unit_system* restrict us, - const struct cosmology* restrict cosmo, - const struct star_formation* data, - const struct part* restrict p, - struct xpart* restrict xp) {} - -/** - * @brief Sets the star_formation properties of the (x-)particles to a valid - * start state. - * - * Nothing to do here. We do not need to compute any quantity in the hydro - * density loop for the EAGLE star formation model. - * - * @param p Pointer to the particle data. - * @param xp Pointer to the extended particle data. - * @param data The global star_formation information. - */ -__attribute__((always_inline)) INLINE static void star_formation_init_part( - struct part* restrict p, struct xpart* restrict xp, - const struct star_formation* data) {} - #endif /* SWIFT_EAGLE_STAR_FORMATION_H */ diff --git a/src/star_formation/EAGLE/star_formation_iact.h b/src/star_formation/EAGLE/star_formation_iact.h deleted file mode 100644 index ab917cbe7aa67cad93a92a4b24212c5f1dcf3aeb..0000000000000000000000000000000000000000 --- a/src/star_formation/EAGLE/star_formation_iact.h +++ /dev/null @@ -1,71 +0,0 @@ -/******************************************************************************* - * This file is part of SWIFT. - * Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk) - * - * 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_EAGLE_STAR_FORMATION_IACT_H -#define SWIFT_EAGLE_STAR_FORMATION_IACT_H - -/** - * @file EAGLE/star_formation_iact.h - * @brief Density computation - */ - -/** - * @brief do star_formation computation after the runner_iact_density (symmetric - * version) - * - * @param r2 Comoving square distance between the two particles. - * @param dx Comoving vector separating both particles (pi - pj). - * @param hi Comoving smoothing-length of particle i. - * @param hj Comoving smoothing-length of particle j. - * @param pi First particle. - * @param pj Second particle. - * @param a Current scale factor. - * @param H Current Hubble parameter. - */ -__attribute__((always_inline)) INLINE static void runner_iact_star_formation( - float r2, const float *dx, float hi, float hj, struct part *restrict pi, - struct part *restrict pj, float a, float H) { - - /* Nothing to do here. We do not need to compute any quantity in the hydro - density loop for the EAGLE star formation model. */ -} - -/** - * @brief do star_formation computation after the runner_iact_density (non - * symmetric version) - * - * @param r2 Comoving square distance between the two particles. - * @param dx Comoving vector separating both particles (pi - pj). - * @param hi Comoving smoothing-length of particle i. - * @param hj Comoving smoothing-length of particle j. - * @param pi First particle. - * @param pj Second particle (not updated). - * @param a Current scale factor. - * @param H Current Hubble parameter. - */ -__attribute__((always_inline)) INLINE static void -runner_iact_nonsym_star_formation(float r2, const float *dx, float hi, float hj, - struct part *restrict pi, - const struct part *restrict pj, float a, - float H) { - - /* Nothing to do here. We do not need to compute any quantity in the hydro - density loop for the EAGLE star formation model. */ -} - -#endif /* SWIFT_EAGLE_STAR_FORMATION_IACT_H */ diff --git a/src/star_formation/EAGLE/star_formation_struct.h b/src/star_formation/EAGLE/star_formation_struct.h index 8caac49d4b57652c5db9ae93e3789dc690e6d23f..41247e160a3eddbc9184c59b67cfa2a1d7259a05 100644 --- a/src/star_formation/EAGLE/star_formation_struct.h +++ b/src/star_formation/EAGLE/star_formation_struct.h @@ -29,6 +29,4 @@ struct star_formation_xpart_data { float SFR; }; -struct star_formation_part_data {}; - #endif /* SWIFT_EAGLE_STAR_FORMATION_STRUCT_H */ diff --git a/src/star_formation/GEAR/star_formation.h b/src/star_formation/GEAR/star_formation.h index ac423a51865609460e870f65b1eeeb266182e2ef..5fc3380fe6869bd5bcb9435fb0c129ac6fc0aad2 100644 --- a/src/star_formation/GEAR/star_formation.h +++ b/src/star_formation/GEAR/star_formation.h @@ -70,7 +70,7 @@ INLINE static int star_formation_is_star_forming( } /* Get the required variables */ - const float sigma2 = p->sf_data.sigma2; + const float sigma2 = p->pressure_floor_data.sigma2 * cosmo->a * cosmo->a; const float n_jeans_2_3 = starform->n_jeans_2_3; const float h = p->h; @@ -223,79 +223,4 @@ INLINE static void starformation_print_backend( message("Star formation law is 'GEAR'"); } -/** - * @brief Finishes the density calculation. - * - * @param p The particle to act upon - * @param xp The extended particle data to act upon - * @param sf The global star_formation information. - * @param cosmo The current cosmological model. - */ -__attribute__((always_inline)) INLINE static void star_formation_end_density( - struct part* restrict p, const struct star_formation* sf, - const struct cosmology* cosmo) { - - // TODO move into pressure floor - /* To finish the turbulence estimation we devide by the density */ - p->sf_data.sigma2 /= - pow_dimension(p->h) * hydro_get_physical_density(p, cosmo); - - /* Add the cosmological factor */ - p->sf_data.sigma2 *= cosmo->a * cosmo->a; -} - -/** - * @brief Sets all particle fields to sensible values when the #part has 0 ngbs. - * - * @param p The particle to act upon - * @param xp The extended particle data to act upon - * @param cd #star_formation containing star_formation informations. - * @param cosmo The current cosmological model. - */ -__attribute__((always_inline)) INLINE static void -star_formation_part_has_no_neighbours(struct part* restrict p, - struct xpart* restrict xp, - const struct star_formation* cd, - const struct cosmology* cosmo) { - - // TODO move into pressure floor - /* If part has 0 neighbours, the estimation of turbulence is 0 */ - p->sf_data.sigma2 = 0.f; -} - -/** - * @brief Sets the star_formation properties of the (x-)particles to a valid - * start state. - * - * @param p Pointer to the particle data. - * @param xp Pointer to extended particle data - * @param data The global star_formation information. - */ -__attribute__((always_inline)) INLINE static void star_formation_init_part( - struct part* restrict p, struct xpart* restrict xp, - const struct star_formation* data) { - p->sf_data.sigma2 = 0.f; -} - -/** - * @brief Sets the star_formation properties of the (x-)particles to a valid - * start state. - * @param phys_const The physical constant in internal units. - * @param us The unit system. - * @param cosmo The current cosmological model. - * @param data The global star_formation information used for this run. - * @param p Pointer to the particle data. - */ -__attribute__((always_inline)) INLINE static void -star_formation_first_init_part(const struct phys_const* restrict phys_const, - const struct unit_system* restrict us, - const struct cosmology* restrict cosmo, - const struct star_formation* data, - struct part* restrict p, - struct xpart* restrict xp) { - - /* Nothing special here */ - star_formation_init_part(p, xp, data); -} - #endif /* SWIFT_GEAR_STAR_FORMATION_H */ diff --git a/src/star_formation/GEAR/star_formation_io.h b/src/star_formation/GEAR/star_formation_io.h index 3e021f7844c1deaeca40d7144d6f7b69cb6c2bdb..7e6dcefce633fd55f421e3c8672eaf8e6f18579c 100644 --- a/src/star_formation/GEAR/star_formation_io.h +++ b/src/star_formation/GEAR/star_formation_io.h @@ -56,11 +56,6 @@ INLINE static void starformation_init_backend( const struct unit_system* us, const struct hydro_props* hydro_props, struct star_formation* starform) { - // TODO move into pressure floor - starform->n_jeans_2_3 = - parser_get_param_float(parameter_file, "GEARStarFormation:NJeans"); - starform->n_jeans_2_3 = pow(starform->n_jeans_2_3, 2. / 3.); - /* Star formation efficiency */ starform->star_formation_efficiency = parser_get_param_double( parameter_file, "GEARStarFormation:star_formation_efficiency"); @@ -69,6 +64,9 @@ INLINE static void starformation_init_backend( starform->maximal_temperature = parser_get_param_double( parameter_file, "GEARStarFormation:maximal_temperature"); + /* Get the jeans factor */ + starform->n_jeans_2_3 = pow(pressure_floor_props.n_jeans, 2. / 3.); + /* Apply unit change */ starform->maximal_temperature *= units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE); diff --git a/src/star_formation/GEAR/star_formation_struct.h b/src/star_formation/GEAR/star_formation_struct.h index 2e9a7548f83ca6ae9bb78ee7bcf4be69a0a31489..50a735ff45f27ccb197ae6089f13237f67a59e3f 100644 --- a/src/star_formation/GEAR/star_formation_struct.h +++ b/src/star_formation/GEAR/star_formation_struct.h @@ -25,20 +25,13 @@ */ struct star_formation_xpart_data {}; -struct star_formation_part_data { - // TODO move it to the pressure floor - /*! Estimation of local turbulence (squared) */ - float sigma2; -}; - /** * @brief Global star formation properties */ struct star_formation { - // TODO move it to pressure floor - /*! Number of particle required to resolved the Jeans criterion (at power 2/3) - */ + /*! Number of particle required to resolved the + * Jeans criterion (at power 2/3) */ float n_jeans_2_3; /*! Maximal temperature for forming a star */ diff --git a/src/star_formation/none/star_formation.h b/src/star_formation/none/star_formation.h index 6b83010ec5f8935778af8c9ad21010ff3452fc0e..0f53e951cb5842e5be3bb9bbe64eb6686f822b1e 100644 --- a/src/star_formation/none/star_formation.h +++ b/src/star_formation/none/star_formation.h @@ -163,64 +163,4 @@ INLINE static void starformation_print_backend( message("Star formation law is 'No Star Formation'"); } -/** - * @brief Finishes the density calculation. - * - * @param p The particle to act upon - * @param cd The global star_formation information. - * @param cosmo The current cosmological model. - */ -__attribute__((always_inline)) INLINE static void star_formation_end_density( - struct part* restrict p, const struct star_formation* cd, - const struct cosmology* cosmo) {} - -/** - * @brief Sets all particle fields to sensible values when the #part has 0 ngbs. - * - * @param p The particle to act upon - * @param xp The extended particle data to act upon - * @param cd #star_formation containing star_formation informations. - * @param cosmo The current cosmological model. - */ -__attribute__((always_inline)) INLINE static void -star_formation_part_has_no_neighbours(struct part* restrict p, - struct xpart* restrict xp, - const struct star_formation* cd, - const struct cosmology* cosmo) {} - -/** - * @brief Sets the star_formation properties of the (x-)particles to a valid - * start state. - * - * Nothing to do here. - * - * @param phys_const The physical constant in internal units. - * @param us The unit system. - * @param cosmo The current cosmological model. - * @param data The global star_formation information used for this run. - * @param p Pointer to the particle data. - * @param xp Pointer to the extended particle data. - */ -__attribute__((always_inline)) INLINE static void -star_formation_first_init_part(const struct phys_const* restrict phys_const, - const struct unit_system* restrict us, - const struct cosmology* restrict cosmo, - const struct star_formation* data, - const struct part* restrict p, - struct xpart* restrict xp) {} - -/** - * @brief Sets the star_formation properties of the (x-)particles to a valid - * start state. - * - * Nothing to do here. - * - * @param p Pointer to the particle data. - * @param xp Pointer to the extended particle data. - * @param data The global star_formation information. - */ -__attribute__((always_inline)) INLINE static void star_formation_init_part( - struct part* restrict p, struct xpart* restrict xp, - const struct star_formation* data) {} - #endif /* SWIFT_NONE_STAR_FORMATION_H */ diff --git a/src/star_formation/none/star_formation_struct.h b/src/star_formation/none/star_formation_struct.h index 2f5241a58caf1ca70fa98a40d467c8ff5a3237f7..27a2adaf83d0a02a0d08e7eef8b45bea630689e4 100644 --- a/src/star_formation/none/star_formation_struct.h +++ b/src/star_formation/none/star_formation_struct.h @@ -25,10 +25,4 @@ */ struct star_formation_xpart_data {}; -/** - * @brief Star-formation-related properties stored in the particle - * data. - */ -struct star_formation_part_data {}; - #endif /* SWIFT_NONE_STAR_FORMATION_STRUCT_H */ diff --git a/src/star_formation_struct.h b/src/star_formation_struct.h index 2a62d284b435c353525311979b343754856364e8..92386d532fb7e0ad445477bf9e3ec35fe597fe2f 100644 --- a/src/star_formation_struct.h +++ b/src/star_formation_struct.h @@ -27,7 +27,7 @@ /* Config parameters. */ #include "../config.h" -/* Import the right cooling definition */ +/* Import the right star formation definition */ #if defined(STAR_FORMATION_NONE) #include "./star_formation/none/star_formation_struct.h" #elif defined(STAR_FORMATION_EAGLE) diff --git a/src/stars/GEAR/stars_io.h b/src/stars/GEAR/stars_io.h index ebd72aa50a4194bf8c6f747e55d265ace0550c35..dc40e86e29b19a370d549576fd955464bf0e609d 100644 --- a/src/stars/GEAR/stars_io.h +++ b/src/stars/GEAR/stars_io.h @@ -108,7 +108,7 @@ INLINE static void stars_write_particles(const struct spart *sparts, "Temperatures at the time of birth of the gas " "particles that turned into stars"); - list[7] = io_make_output_field("BirthMasses", FLOAT, 1, UNIT_CONV_MASS, 0.f, + list[8] = io_make_output_field("BirthMasses", FLOAT, 1, UNIT_CONV_MASS, 0.f, sparts, birth.mass, "Masses of the star particles at birth time"); diff --git a/src/timestep.h b/src/timestep.h index c2b1a10fcb3b0426e7c34625d65c1fd5353d25e9..cd9faaea612c8a666ae9077ccc5e3d85fd4677f9 100644 --- a/src/timestep.h +++ b/src/timestep.h @@ -146,8 +146,14 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep( new_dt_grav = min(new_dt_self_grav, new_dt_ext_grav); } - /* Final time-step is minimum of hydro and gravity */ - float new_dt = min3(new_dt_hydro, new_dt_cooling, new_dt_grav); + /* Compute the next timestep (chemistry condition, e.g. diffusion) */ + const float new_dt_chemistry = + chemistry_timestep(e->physical_constants, e->cosmology, e->internal_units, + e->hydro_properties, e->chemistry, p); + + /* Final time-step is minimum of hydro, gravity and subgrid */ + float new_dt = + min4(new_dt_hydro, new_dt_cooling, new_dt_grav, new_dt_chemistry); /* Limit change in smoothing length */ const float dt_h_change = diff --git a/src/tools.c b/src/tools.c index bd467e1841b78d45a74b037c63e8ed58dd92183f..0643fb7922c0e3d56b70c6a0d1a30e3ca13154c6 100644 --- a/src/tools.c +++ b/src/tools.c @@ -45,8 +45,8 @@ #include "hydro.h" #include "part.h" #include "periodic.h" +#include "pressure_floor_iact.h" #include "runner.h" -#include "star_formation_iact.h" #include "stars.h" /** @@ -224,7 +224,7 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) { /* Interact */ runner_iact_nonsym_density(r2, dx, hi, pj->h, pi, pj, a, H); runner_iact_nonsym_chemistry(r2, dx, hi, pj->h, pi, pj, a, H); - runner_iact_nonsym_star_formation(r2, dx, hi, pj->h, pi, pj, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hi, pj->h, pi, pj, a, H); } } } @@ -257,7 +257,7 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) { /* Interact */ runner_iact_nonsym_density(r2, dx, hj, pi->h, pj, pi, a, H); runner_iact_nonsym_chemistry(r2, dx, hj, pi->h, pj, pi, a, H); - runner_iact_nonsym_star_formation(r2, dx, hj, pi->h, pj, pi, a, H); + runner_iact_nonsym_pressure_floor(r2, dx, hj, pi->h, pj, pi, a, H); } } } @@ -535,7 +535,7 @@ void self_all_density(struct runner *r, struct cell *ci) { /* Interact */ runner_iact_nonsym_density(r2, dxi, hi, hj, pi, pj, a, H); runner_iact_nonsym_chemistry(r2, dxi, hi, hj, pi, pj, a, H); - runner_iact_nonsym_star_formation(r2, dxi, hi, hj, pi, pj, a, H); + runner_iact_nonsym_pressure_floor(r2, dxi, hi, hj, pi, pj, a, H); } /* Hit or miss? */ @@ -548,7 +548,7 @@ void self_all_density(struct runner *r, struct cell *ci) { /* Interact */ runner_iact_nonsym_density(r2, dxi, hj, hi, pj, pi, a, H); runner_iact_nonsym_chemistry(r2, dxi, hj, hi, pj, pi, a, H); - runner_iact_nonsym_star_formation(r2, dxi, hj, hi, pj, pi, a, H); + runner_iact_nonsym_pressure_floor(r2, dxi, hj, hi, pj, pi, a, H); } } }