diff --git a/configure.ac b/configure.ac
index 24b103480a7bbd2001d0ee473d6f2f7751a6c85f..0f8feb9907441e4ca6943f7786c70f7c2104e6d0 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1327,6 +1327,7 @@ with_subgrid_cooling=none
with_subgrid_chemistry=none
with_subgrid_tracers=none
with_subgrid_entropy_floor=none
+with_subgrid_pressure_floor=none
with_subgrid_stars=none
with_subgrid_star_formation=none
with_subgrid_feedback=none
@@ -1338,10 +1339,9 @@ case "$with_subgrid" in
none)
;;
GEAR)
- with_subgrid_cooling=grackle
+ with_subgrid_cooling=grackle_0
with_subgrid_chemistry=GEAR
- with_subgrid_tracers=none
- with_subgrid_entropy_floor=none
+ with_subgrid_pressure_floor=GEAR
with_subgrid_stars=GEAR
with_subgrid_star_formation=GEAR
with_subgrid_feedback=none
@@ -1400,6 +1400,12 @@ AC_ARG_WITH([gravity],
[with_gravity="default"]
)
+if test "$with_subgrid" = "EAGLE"; then
+ if test "$with_gravity" = "default"; then
+ with_gravity="with-potential"
+ fi
+fi
+
case "$with_gravity" in
with-potential)
AC_DEFINE([POTENTIAL_GRAVITY], [1], [Gravity scheme with potential calculation])
@@ -1643,7 +1649,8 @@ esac
# Cooling function
AC_ARG_WITH([cooling],
[AS_HELP_STRING([--with-cooling=<function>],
- [cooling function @<:@none, const-du, const-lambda, EAGLE, grackle, grackle1, grackle2, grackle3 default: none@:>@]
+ [cooling function @<:@none, const-du, const-lambda, EAGLE, grackle_* default: none@:>@.
+ For Grackle, you need to provide the primordial chemistry parameter (e.g. grackle_0)]
)],
[with_cooling="$withval"],
[with_cooling="none"]
@@ -1670,21 +1677,10 @@ case "$with_cooling" in
compton)
AC_DEFINE([COOLING_COMPTON], [1], [Compton cooling off the CMB])
;;
- grackle)
- AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
- AC_DEFINE([COOLING_GRACKLE_MODE], [0], [Grackle chemistry network, mode 0])
- ;;
- grackle1)
- AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
- AC_DEFINE([COOLING_GRACKLE_MODE], [1], [Grackle chemistry network, mode 1])
- ;;
- grackle2)
+ grackle_*)
AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
- AC_DEFINE([COOLING_GRACKLE_MODE], [2], [Grackle chemistry network, mode 2])
- ;;
- grackle3)
- AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
- AC_DEFINE([COOLING_GRACKLE_MODE], [3], [Grackle chemistry network, mode 3])
+ primordial_chemistry=${with_cooling:8}
+ AC_DEFINE_UNQUOTED([COOLING_GRACKLE_MODE], [$primordial_chemistry], [Grackle chemistry network])
;;
EAGLE)
AC_DEFINE([COOLING_EAGLE], [1], [Cooling following the EAGLE model])
@@ -1919,6 +1915,35 @@ case "$with_entropy_floor" in
;;
esac
+# Pressure floor
+AC_ARG_WITH([pressure-floor],
+ [AS_HELP_STRING([--with-pressure-floor=<floor>],
+ [pressure floor @<:@none, GEAR, default: none@:>@
+ The hydro model needs to be compatible.]
+ )],
+ [with_pressure_floor="$withval"],
+ [with_pressure_floor="none"]
+)
+if test "$with_subgrid" != "none"; then
+ if test "$with_pressure_floor" != "none"; then
+ AC_MSG_ERROR([Cannot provide with-subgrid and with-pressure-floor together])
+ else
+ with_pressure_floor="$with_subgrid_pressure_floor"
+ fi
+fi
+
+case "$with_pressure_floor" in
+ none)
+ AC_DEFINE([PRESSURE_FLOOR_NONE], [1], [No pressure floor])
+ ;;
+ GEAR)
+ AC_DEFINE([PRESSURE_FLOOR_GEAR], [1], [GEAR pressure floor])
+ ;;
+ *)
+ AC_MSG_ERROR([Unknown pressure floor model])
+ ;;
+esac
+
# Star formation
AC_ARG_WITH([star-formation],
[AS_HELP_STRING([--with-star-formation=<sfm>],
@@ -1964,6 +1989,13 @@ AC_DEFINE_UNQUOTED([SELF_GRAVITY_MULTIPOLE_ORDER], [$with_multipole_order], [Mul
AC_PATH_PROG([GIT_CMD], [git])
AC_SUBST([GIT_CMD])
+# Check that the gravity model is compatible with the subgrid
+if test $with_black_holes = "EAGLE"; then
+ if test $with_gravity != "with-potential"; then
+ AC_MSG_ERROR([The EAGLE BH model needs the gravity scheme to provide potentials. The code must be compile with --with-gravity=with-potential.])
+ fi
+fi
+
# Make the documentation. Add conditional to handle disable option.
DX_INIT_DOXYGEN(libswift,doc/Doxyfile,doc/)
AM_CONDITIONAL([HAVE_DOXYGEN], [test "$ac_cv_path_ac_pt_DX_DOXYGEN" != ""])
@@ -2046,6 +2078,7 @@ AC_MSG_RESULT([
External potential : $with_potential
Entropy floor : $with_entropy_floor
+ Pressure floor : $with_pressure_floor
Cooling function : $with_cooling
Chemistry : $with_chemistry
Tracers : $with_tracers
diff --git a/doc/RTD/source/SubgridModels/EAGLE/index.rst b/doc/RTD/source/SubgridModels/EAGLE/index.rst
index 56736cc4cca4bd9042a2c84d07dcf535f48c566c..cac6716c7568b7c7d2168a77586cc80123df7127 100644
--- a/doc/RTD/source/SubgridModels/EAGLE/index.rst
+++ b/doc/RTD/source/SubgridModels/EAGLE/index.rst
@@ -577,6 +577,11 @@ Black-hole creation
Black-hole accretion
~~~~~~~~~~~~~~~~~~~~
+.. _EAGLE_black_hole_reposition:
+
+Black-hole repositioning
+~~~~~~~~~~~~~~~~~~~~~~~~
+
.. _EAGLE_black_hole_feedback:
AGN feedback
diff --git a/doc/RTD/source/SubgridModels/GEAR/index.rst b/doc/RTD/source/SubgridModels/GEAR/index.rst
index 2a211759bfc4895fd07279b72f78200d6ea47546..152056ef2eb60da189bf18b3b3c9909f88ace6de 100644
--- a/doc/RTD/source/SubgridModels/GEAR/index.rst
+++ b/doc/RTD/source/SubgridModels/GEAR/index.rst
@@ -5,6 +5,29 @@
GEAR model
===========
+Pressure Floor
+~~~~~~~~~~~~~~
+
+In order to avoid the artificial collapse of unresolved clumps, a minimum in pressure is applied to the particles.
+This additional pressure can be seen as the pressure due to unresolved hydrodynamics turbulence and is given by:
+
+.. math::
+ P_\textrm{Jeans} = \frac{\rho}{\gamma} \left( \frac{4}{\pi} G h^2 \rho N_\textrm{Jeans}^{2/3} - \sigma^2 \right)
+
+where :math:`\rho` is the density, :math:`\gamma` the adiabatic index, :math:`G` is the gravitational constant,
+:math:`h` the smoothing length, :math:`N_\textrm{Jeans}` is the number of particle required in order to resolve a clump and
+:math:`\sigma` the velocity dispersion.
+
+
+This must be directly implemented into the hydro schemes, therefore only a subset of schemes (Gadget-2 and Pressure-Energy) are currently implemented.
+In order to implement it, you need equation 12 in `Hopkins 2013 <https://arxiv.org/abs/1206.5006>`_:
+
+.. math::
+ m_i \frac{\mathrm{d}v_i}{\mathrm{d}t} = - \sum_j x_i x_j \left[ \frac{P_i}{y_i^2} f_{ij} \nabla_i W_{ij}(h_i) + \frac{P_j}{y_j^2} f_{ji} \nabla_j W_{ji}(h_j) \right]
+
+and simply replace the :math:`P_i, P_j` by the pressure with the floor.
+Here the :math:`x, y` are simple weights that should never have the pressure floor included even if they are related to the pressure (e.g. pressure-entropy).
+
Cooling: Grackle
~~~~~~~~~~~~~~~~
diff --git a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
index 3270af86afffba13cfac6be7c7c21dd84980adc5..b5122afdba28ffaed1e52cdd8bd2db8a4e566c03 100644
--- a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
@@ -170,3 +170,4 @@ EAGLEAGN:
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
+ max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
\ No newline at end of file
diff --git a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
index 6461982f70c9635652b8253887a9998777f37c69..888009f27761feda91d3cb70df2366b9f04fd3eb 100644
--- a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
@@ -171,3 +171,4 @@ EAGLEAGN:
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
+ max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
\ No newline at end of file
diff --git a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
index c6f7f0e01bd292b03b0eb31c6c1b57d5163dafe7..f2e3b0656541ca698c20a85a308e538871d4fdc4 100644
--- a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
@@ -171,5 +171,4 @@ EAGLEAGN:
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
-
-
+ max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
diff --git a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
index 9205b9f6f560398b046625f6c65fcc887492991d..685d24d20f14f567175e1d7086a65455017199e0 100644
--- a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
@@ -162,3 +162,4 @@ EAGLEAGN:
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
+ max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
\ No newline at end of file
diff --git a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
index ccd89efbead0c511316bd664f6754ed0797c7c50..b782631d7e264e582455f70b4fb399f7db1f4731 100644
--- a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
@@ -169,3 +169,4 @@ EAGLEAGN:
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
+ max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
\ No newline at end of file
diff --git a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
index e2a98725a94c9fd34b45a95530f480562c0aaec8..53afda9445b4f8ee7c258b5c42af2a605eba4755 100644
--- a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
@@ -161,3 +161,4 @@ EAGLEAGN:
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
+ max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
\ No newline at end of file
diff --git a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
index 9136f85424beaf70b614cde0e70b67c134056831..694c1c111d8f5b2455d4cfcdc567b8854b12ff2c 100644
--- a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
@@ -171,3 +171,4 @@ EAGLEAGN:
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
+ max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
\ No newline at end of file
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml
index b108290fcd146461827d5858742bd0971ac66945..94a9082af95e96b161cf5fe469eb082774e5007f 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml
@@ -99,3 +99,6 @@ GrackleCooling:
GearChemistry:
InitialMetallicity: 0.01295
+
+GEARPressureFloor:
+ Jeans_factor: 10. # Number of particles required to suppose a resolved clump and avoid the pressure floor.
diff --git a/examples/main.c b/examples/main.c
index 34bd9ff7ec68941f36da9a4992cb8318972bdd98..cf8734532a9b8e48997729d5d1b086b15c6859fc 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -771,6 +771,13 @@ int main(int argc, char *argv[]) {
else
bzero(&entropy_floor, sizeof(struct entropy_floor_properties));
+ /* Initialise the pressure floor */
+ if (with_hydro)
+ pressure_floor_init(&pressure_floor_props, &prog_const, &us,
+ &hydro_properties, params);
+ else
+ bzero(&pressure_floor_props, sizeof(struct pressure_floor_properties));
+
/* Initialise the stars properties */
if (with_stars)
stars_props_init(&stars_properties, &prog_const, &us, params,
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 8afe5b7b39afacfa23aaf8240c1d70ad33278385..3607690f7cbf07c86eb2d84cd742c734b33bd266 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -288,6 +288,11 @@ EAGLEEntropyFloor:
Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
+# Parameters related to pressure floors ----------------------------------------------
+
+GEARPressureFloor:
+ Jeans_factor: 10. # Number of particles required to suppose a resolved clump and avoid the pressure floor.
+
# Parameters related to cooling function ----------------------------------------------
# Constant du/dt cooling function
@@ -341,6 +346,11 @@ EAGLEChemistry:
# Parameters related to star formation models -----------------------------------------------
+# GEAR star formation model (Revaz and Jablonka 2018)
+GEARStarFormation:
+ star_formation_efficiency: 0.01 # star formation efficiency (c_*)
+ maximal_temperature: 3e4 # Upper limit to the temperature of a star forming particle
+
# EAGLE star formation model (Schaye and Dalla Vecchia 2008)
EAGLEStarFormation:
EOS_density_norm_H_p_cm3: 0.1 # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
@@ -406,4 +416,4 @@ EAGLEAGN:
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
-
+ max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
diff --git a/src/Makefile.am b/src/Makefile.am
index 0f7cef5e3c6b861ebde639d90624a360f1be7047..3d748b135c93cf409464a199232b243c699fabe1 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -79,7 +79,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c engine_maketasks.c
collectgroup.c hydro_space.c equation_of_state.c \
chemistry.c cosmology.c restart.c mesh_gravity.c velociraptor_interface.c \
outputlist.c velociraptor_dummy.c logger_io.c memuse.c fof.c \
- hashmap.c \
+ hashmap.c pressure_floor.c \
$(EAGLE_COOLING_SOURCES) $(EAGLE_FEEDBACK_SOURCES)
# Include files for distribution, not installation.
@@ -98,28 +98,28 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
hydro.h hydro_io.h hydro_parameters.h \
hydro/Minimal/hydro.h hydro/Minimal/hydro_iact.h hydro/Minimal/hydro_io.h \
hydro/Minimal/hydro_debug.h hydro/Minimal/hydro_part.h \
- hydro/Minimal/hydro_parameters.h \
+ hydro/Minimal/hydro_parameters.h \
hydro/Default/hydro.h hydro/Default/hydro_iact.h hydro/Default/hydro_io.h \
hydro/Default/hydro_debug.h hydro/Default/hydro_part.h \
- hydro/Default/hydro_parameters.h \
+ hydro/Default/hydro_parameters.h \
hydro/Gadget2/hydro.h hydro/Gadget2/hydro_iact.h hydro/Gadget2/hydro_io.h \
hydro/Gadget2/hydro_debug.h hydro/Gadget2/hydro_part.h \
- hydro/Gadget2/hydro_parameters.h \
+ hydro/Gadget2/hydro_parameters.h \
hydro/PressureEntropy/hydro.h hydro/PressureEntropy/hydro_iact.h hydro/PressureEntropy/hydro_io.h \
hydro/PressureEntropy/hydro_debug.h hydro/PressureEntropy/hydro_part.h \
- hydro/PressureEntropy/hydro_parameters.h \
+ hydro/PressureEntropy/hydro_parameters.h \
hydro/PressureEnergy/hydro.h hydro/PressureEnergy/hydro_iact.h hydro/PressureEnergy/hydro_io.h \
hydro/PressureEnergy/hydro_debug.h hydro/PressureEnergy/hydro_part.h \
- hydro/PressureEnergy/hydro_parameters.h \
+ hydro/PressureEnergy/hydro_parameters.h \
hydro/PressureEnergyMorrisMonaghanAV/hydro.h hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h \
hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h \
- hydro/PressureEnergyMorrisMonaghanAV/hydro_parameters.h \
+ hydro/PressureEnergyMorrisMonaghanAV/hydro_parameters.h \
hydro/AnarchyPU/hydro.h hydro/PressureEnergy/hydro_iact.h hydro/PressureEnergy/hydro_io.h \
hydro/AnarchyPU/hydro_debug.h hydro/PressureEnergy/hydro_part.h \
- hydro/AnarchyPU/hydro_parameters.h \
+ hydro/AnarchyPU/hydro_parameters.h \
hydro/AnarchyDU/hydro.h hydro/PressureEnergy/hydro_iact.h hydro/PressureEnergy/hydro_io.h \
hydro/AnarchyDU/hydro_debug.h hydro/PressureEnergy/hydro_part.h \
- hydro/AnarchyDU/hydro_parameters.h \
+ hydro/AnarchyDU/hydro_parameters.h \
hydro/GizmoMFV/hydro.h hydro/GizmoMFV/hydro_iact.h \
hydro/GizmoMFV/hydro_io.h hydro/GizmoMFV/hydro_debug.h \
hydro/GizmoMFV/hydro_part.h \
@@ -131,7 +131,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
hydro/GizmoMFV/hydro_slope_limiters.h \
hydro/GizmoMFV/hydro_unphysical.h \
hydro/GizmoMFV/hydro_velocities.h \
- hydro/GizmoMFV/hydro_parameters.h \
+ hydro/GizmoMFV/hydro_parameters.h \
hydro/GizmoMFM/hydro.h hydro/GizmoMFM/hydro_iact.h \
hydro/GizmoMFM/hydro_io.h hydro/GizmoMFM/hydro_debug.h \
hydro/GizmoMFM/hydro_part.h \
@@ -142,7 +142,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
hydro/GizmoMFM/hydro_slope_limiters_face.h \
hydro/GizmoMFM/hydro_slope_limiters.h \
hydro/GizmoMFM/hydro_unphysical.h \
- hydro/GizmoMFM/hydro_parameters.h \
+ hydro/GizmoMFM/hydro_parameters.h \
hydro/Shadowswift/hydro_debug.h \
hydro/Shadowswift/hydro_gradients.h hydro/Shadowswift/hydro.h \
hydro/Shadowswift/hydro_iact.h \
@@ -190,7 +190,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
cooling/const_lambda/cooling_io.h \
cooling/grackle/cooling.h cooling/grackle/cooling_struct.h \
cooling/grackle/cooling_io.h \
- cooling/EAGLE/cooling.h cooling/EAGLE/cooling_struct.h \
+ cooling/EAGLE/cooling.h cooling/EAGLE/cooling_struct.h cooling/EAGLE/cooling_tables.h \
cooling/EAGLE/cooling_io.h cooling/EAGLE/interpolate.h cooling/EAGLE/cooling_rates.h \
chemistry/none/chemistry.h \
chemistry/none/chemistry_io.h \
@@ -221,7 +221,10 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
black_holes/Default/black_holes_struct.h \
black_holes/EAGLE/black_holes.h black_holes/EAGLE/black_holes_io.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_properties.h \
+ black_holes/EAGLE/black_holes_struct.h \
+ pressure_floor.h \
+ pressure_floor/GEAR/pressure_floor.h pressure_floor/none/pressure_floor.h
# Sources and flags for regular library
diff --git a/src/black_holes/Default/black_holes.h b/src/black_holes/Default/black_holes.h
index 56866764fe4b65603309856bfb2c5ae6b8d82f00..082134cc78f5a7f5e0aa9adeb231c92b9d66c233 100644
--- a/src/black_holes/Default/black_holes.h
+++ b/src/black_holes/Default/black_holes.h
@@ -187,6 +187,20 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
const struct phys_const* constants, const struct cosmology* cosmo,
const double dt) {}
+/**
+ * @brief Finish the calculation of the new BH position.
+ *
+ * Nothing to do here.
+ *
+ * @param bp The black hole particle.
+ * @param props The properties of the black hole scheme.
+ * @param constants The physical constants (in internal units).
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void black_holes_end_reposition(
+ struct bpart* restrict bp, const struct black_holes_props* props,
+ const struct phys_const* constants, const struct cosmology* cosmo) {}
+
/**
* @brief Reset acceleration fields of a particle
*
diff --git a/src/black_holes/Default/black_holes_iact.h b/src/black_holes/Default/black_holes_iact.h
index ccb1e8c90cd8fef42dd0bb5b6e2d9c6eb9382226..b26b8a50f40599402f4ea9886dbbd7495e8efd59 100644
--- a/src/black_holes/Default/black_holes_iact.h
+++ b/src/black_holes/Default/black_holes_iact.h
@@ -30,16 +30,15 @@
* @param pj Second particle (gas, not updated).
* @param xpj The extended data of the second particle (not updated).
* @param cosmo The cosmological model.
+ * @param grav_props The properties of the gravity scheme (softening, G, ...).
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_bh_gas_density(const float r2, const float *dx,
- const float hi, const float hj,
- struct bpart *restrict bi,
- const struct part *restrict pj,
- const struct xpart *restrict xpj,
- const struct cosmology *cosmo,
- const integertime_t ti_current) {
+runner_iact_nonsym_bh_gas_density(
+ const float r2, const float *dx, const float hi, const float hj,
+ struct bpart *restrict bi, const struct part *restrict pj,
+ const struct xpart *restrict xpj, const struct cosmology *cosmo,
+ const struct gravity_props *grav_props, const integertime_t ti_current) {
float wi, wi_dx;
@@ -80,16 +79,15 @@ runner_iact_nonsym_bh_gas_density(const float r2, const float *dx,
* @param pj Second particle (gas)
* @param xpj The extended data of the second particle.
* @param cosmo The cosmological model.
+ * @param grav_props The properties of the gravity scheme (softening, G, ...).
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
- const float hi, const float hj,
- const struct bpart *restrict bi,
- struct part *restrict pj,
- struct xpart *restrict xpj,
- const struct cosmology *cosmo,
- const integertime_t ti_current) {}
+runner_iact_nonsym_bh_gas_swallow(
+ const float r2, const float *dx, const float hi, const float hj,
+ const struct bpart *restrict bi, struct part *restrict pj,
+ struct xpart *restrict xpj, const struct cosmology *cosmo,
+ const struct gravity_props *grav_props, const integertime_t ti_current) {}
/**
* @brief Swallowing interaction between two BH particles (non-symmetric).
@@ -104,7 +102,7 @@ runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
* @param bi First particle (black hole).
* @param bj Second particle (black hole)
* @param cosmo The cosmological model.
- * @param grav_props The properties of the gravity scheme (softening, G, ...)
+ * @param grav_props The properties of the gravity scheme (softening, G, ...).
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
@@ -127,16 +125,15 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
* @param pj Second particle (gas)
* @param xpj The extended data of the second particle.
* @param cosmo The cosmological model.
+ * @param grav_props The properties of the gravity scheme (softening, G, ...).
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_bh_gas_feedback(const float r2, const float *dx,
- const float hi, const float hj,
- struct bpart *restrict bi,
- struct part *restrict pj,
- struct xpart *restrict xpj,
- const struct cosmology *cosmo,
- const integertime_t ti_current) {
+runner_iact_nonsym_bh_gas_feedback(
+ const float r2, const float *dx, const float hi, const float hj,
+ struct bpart *restrict bi, struct part *restrict pj,
+ struct xpart *restrict xpj, const struct cosmology *cosmo,
+ const struct gravity_props *grav_props, const integertime_t ti_current) {
#ifdef DEBUG_INTERACTIONS_BH
/* Update ngb counters */
if (si->num_ngb_force < MAX_NUM_OF_NEIGHBOURS_BH)
diff --git a/src/black_holes/Default/black_holes_io.h b/src/black_holes/Default/black_holes_io.h
index 5288ae595663871bb3187782712d52708fac7931..53322da1dd90a98fb0f1f2529774a20dc718b30c 100644
--- a/src/black_holes/Default/black_holes_io.h
+++ b/src/black_holes/Default/black_holes_io.h
@@ -102,6 +102,7 @@ INLINE static void convert_bpart_vel(const struct engine* e,
* @param bparts The b-particle array.
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
+ * @param with_cosmology Are we running a cosmological simulation?
*/
INLINE static void black_holes_write_particles(const struct bpart* bparts,
struct io_props* list,
diff --git a/src/black_holes/EAGLE/black_holes.h b/src/black_holes/EAGLE/black_holes.h
index 3bf4e2295e1e5f13993e94f76ef6b83882a5e2e7..a1b41f0954e1bec890329779d0edb7fcd31265f1 100644
--- a/src/black_holes/EAGLE/black_holes.h
+++ b/src/black_holes/EAGLE/black_holes.h
@@ -24,6 +24,7 @@
#include "black_holes_struct.h"
#include "cosmology.h"
#include "dimension.h"
+#include "gravity.h"
#include "kernel_hydro.h"
#include "minmax.h"
#include "physical_constants.h"
@@ -89,16 +90,43 @@ __attribute__((always_inline)) INLINE static void black_holes_init_bpart(
bp->circular_velocity_gas[2] = 0.f;
bp->ngb_mass = 0.f;
bp->num_ngbs = 0;
+ bp->reposition.x[0] = -FLT_MAX;
+ bp->reposition.x[1] = -FLT_MAX;
+ bp->reposition.x[2] = -FLT_MAX;
+ bp->reposition.min_potential = FLT_MAX;
}
/**
* @brief Predict additional particle fields forward in time when drifting
*
+ * The fields do not get predicted but we move the BH to its new position
+ * if a new one was calculated in the repositioning loop.
+ *
* @param bp The particle
* @param dt_drift The drift time-step for positions.
*/
__attribute__((always_inline)) INLINE static void black_holes_predict_extra(
- struct bpart* restrict bp, float dt_drift) {}
+ struct bpart* restrict bp, float dt_drift) {
+
+ /* Are we doing some repositioning? */
+ if (bp->reposition.min_potential != FLT_MAX) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (bp->reposition.x[0] == -FLT_MAX || bp->reposition.x[1] == -FLT_MAX ||
+ bp->reposition.x[2] == -FLT_MAX) {
+ error("Something went wrong with the new repositioning position");
+ }
+#endif
+
+ bp->x[0] = bp->reposition.x[0];
+ bp->x[1] = bp->reposition.x[1];
+ bp->x[2] = bp->reposition.x[2];
+
+ bp->gpart->x[0] = bp->reposition.x[0];
+ bp->gpart->x[1] = bp->reposition.x[1];
+ bp->gpart->x[2] = bp->reposition.x[2];
+ }
+}
/**
* @brief Sets the values to be predicted in the drifts to their values at a
@@ -430,6 +458,35 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
}
}
+/**
+ * @brief Finish the calculation of the new BH position.
+ *
+ * Here, we check that the BH should indeed be moved in the next drift.
+ *
+ * @param bp The black hole particle.
+ * @param props The properties of the black hole scheme.
+ * @param constants The physical constants (in internal units).
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void black_holes_end_reposition(
+ struct bpart* restrict bp, const struct black_holes_props* props,
+ const struct phys_const* constants, const struct cosmology* cosmo) {
+
+ const float potential = gravity_get_comoving_potential(bp->gpart);
+
+ /* Is the potential lower (i.e. the BH is at the bottom already)
+ * OR is the BH massive enough that we don't reposition? */
+ if (potential < bp->reposition.min_potential ||
+ bp->subgrid_mass > props->max_reposition_mass) {
+
+ /* No need to reposition */
+ bp->reposition.min_potential = FLT_MAX;
+ bp->reposition.x[0] = -FLT_MAX;
+ bp->reposition.x[1] = -FLT_MAX;
+ bp->reposition.x[2] = -FLT_MAX;
+ }
+}
+
/**
* @brief Reset acceleration fields of a particle
*
diff --git a/src/black_holes/EAGLE/black_holes_iact.h b/src/black_holes/EAGLE/black_holes_iact.h
index fb787df2e34fcfec7cb2192c45f5bf5b4040ebe5..8fa27f162dacf11f7127e7cd473b36c4e2e00236 100644
--- a/src/black_holes/EAGLE/black_holes_iact.h
+++ b/src/black_holes/EAGLE/black_holes_iact.h
@@ -20,6 +20,7 @@
#define SWIFT_EAGLE_BH_IACT_H
/* Local includes */
+#include "gravity.h"
#include "hydro.h"
#include "random.h"
#include "space.h"
@@ -35,16 +36,15 @@
* @param pj Second particle (gas, not updated).
* @param xpj The extended data of the second particle (not updated).
* @param cosmo The cosmological model.
+ * @param grav_props The properties of the gravity scheme (softening, G, ...).
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_bh_gas_density(const float r2, const float *dx,
- const float hi, const float hj,
- struct bpart *restrict bi,
- const struct part *restrict pj,
- const struct xpart *restrict xpj,
- const struct cosmology *cosmo,
- const integertime_t ti_current) {
+runner_iact_nonsym_bh_gas_density(
+ const float r2, const float *dx, const float hi, const float hj,
+ struct bpart *restrict bi, const struct part *restrict pj,
+ const struct xpart *restrict xpj, const struct cosmology *cosmo,
+ const struct gravity_props *grav_props, const integertime_t ti_current) {
float wi, wi_dx;
@@ -117,16 +117,15 @@ runner_iact_nonsym_bh_gas_density(const float r2, const float *dx,
* @param pj Second particle (gas)
* @param xpj The extended data of the second particle.
* @param cosmo The cosmological model.
+ * @param grav_props The properties of the gravity scheme (softening, G, ...).
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
- const float hi, const float hj,
- const struct bpart *restrict bi,
- struct part *restrict pj,
- struct xpart *restrict xpj,
- const struct cosmology *cosmo,
- const integertime_t ti_current) {
+runner_iact_nonsym_bh_gas_swallow(
+ const float r2, const float *dx, const float hi, const float hj,
+ struct bpart *restrict bi, struct part *restrict pj,
+ struct xpart *restrict xpj, const struct cosmology *cosmo,
+ const struct gravity_props *grav_props, const integertime_t ti_current) {
float wi;
@@ -140,6 +139,44 @@ runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
const float ui = r * hi_inv;
kernel_eval(ui, &wi);
+ /* Start by checking the repositioning criteria */
+
+ /* Note the factor 9 is taken from EAGLE. Will be turned into a parameter */
+ const float max_dist_repos2 =
+ kernel_gravity_softening_plummer_equivalent_inv *
+ kernel_gravity_softening_plummer_equivalent_inv * 9.f *
+ grav_props->epsilon_cur2;
+
+ /* This gas neighbour is close enough that we can consider it's potential
+ for repositioning */
+ if (r2 < max_dist_repos2) {
+
+ /* Compute relative peculiar velocity between the two BHs
+ * Recall that in SWIFT v is (v_pec * a) */
+ const float delta_v[3] = {bi->v[0] - pj->v[0], bi->v[1] - pj->v[1],
+ bi->v[2] - pj->v[2]};
+ const float v2 = delta_v[0] * delta_v[0] + delta_v[1] * delta_v[1] +
+ delta_v[2] * delta_v[2];
+
+ const float v2_pec = v2 * cosmo->a2_inv;
+
+ /* Check the velocity criterion */
+ if (v2_pec < 0.25f * bi->sound_speed_gas * bi->sound_speed_gas) {
+
+ const float potential = gravity_get_comoving_potential(pj->gpart);
+
+ /* Is the potential lower? */
+ if (potential < bi->reposition.min_potential) {
+
+ /* Store this as our new best */
+ bi->reposition.min_potential = potential;
+ bi->reposition.x[0] = pj->x[0];
+ bi->reposition.x[1] = pj->x[1];
+ bi->reposition.x[2] = pj->x[2];
+ }
+ }
+ }
+
/* Is the BH hungry? */
if (bi->subgrid_mass > bi->mass) {
@@ -188,13 +225,13 @@ runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
* @param bi First particle (black hole).
* @param bj Second particle (black hole)
* @param cosmo The cosmological model.
- * @param grav_props The properties of the gravity scheme (softening, G, ...)
+ * @param grav_props The properties of the gravity scheme (softening, G, ...).
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
const float hi, const float hj,
- const struct bpart *restrict bi,
+ struct bpart *restrict bi,
struct bpart *restrict bj,
const struct cosmology *cosmo,
const struct gravity_props *grav_props,
@@ -209,6 +246,33 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
const float v2_pec = v2 * cosmo->a2_inv;
+ /* Note the factor 9 is taken from EAGLE. Will be turned into a parameter */
+ const float max_dist_repos2 =
+ kernel_gravity_softening_plummer_equivalent_inv *
+ kernel_gravity_softening_plummer_equivalent_inv * 9.f *
+ grav_props->epsilon_cur2;
+
+ /* This gas neighbour is close enough that we can consider it's potential
+ for repositioning */
+ if (r2 < max_dist_repos2) {
+
+ /* Check the velocity criterion */
+ if (v2_pec < 0.25f * bi->sound_speed_gas * bi->sound_speed_gas) {
+
+ const float potential = gravity_get_comoving_potential(bj->gpart);
+
+ /* Is the potential lower? */
+ if (potential < bi->reposition.min_potential) {
+
+ /* Store this as our new best */
+ bi->reposition.min_potential = potential;
+ bi->reposition.x[0] = bj->x[0];
+ bi->reposition.x[1] = bj->x[1];
+ bi->reposition.x[2] = bj->x[2];
+ }
+ }
+ }
+
/* Find the most massive of the two BHs */
float M = bi->subgrid_mass;
float h = hi;
@@ -218,10 +282,11 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
}
/* Note the factor 9 is taken from EAGLE. Will be turned into a parameter */
- const float max_distance2 = kernel_gravity_softening_plummer_equivalent_inv *
- kernel_gravity_softening_plummer_equivalent_inv *
- grav_props->epsilon_baryon_cur *
- grav_props->epsilon_baryon_cur * 9.f;
+ const float max_dist_merge2 =
+ kernel_gravity_softening_plummer_equivalent_inv *
+ kernel_gravity_softening_plummer_equivalent_inv * 9.f *
+ grav_props->epsilon_baryon_cur *
+ grav_props->epsilon_baryon_cur;
const float G_Newton = grav_props->G_Newton;
@@ -235,7 +300,7 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
/* Merge if gravitationally bound AND if within max distance
* Note that we use the kernel support here as the size and not just the
* smoothing length */
- if (v2_pec < G_Newton * M / (kernel_gamma * h) && (r2 < max_distance2)) {
+ if (v2_pec < G_Newton * M / (kernel_gamma * h) && (r2 < max_dist_merge2)) {
/* This particle is swallowed by the BH with the largest ID of all the
* candidates wanting to swallow it */
@@ -270,16 +335,15 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
* @param pj Second particle (gas)
* @param xpj The extended data of the second particle.
* @param cosmo The cosmological model.
+ * @param grav_props The properties of the gravity scheme (softening, G, ...).
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_bh_gas_feedback(const float r2, const float *dx,
- const float hi, const float hj,
- const struct bpart *restrict bi,
- struct part *restrict pj,
- struct xpart *restrict xpj,
- const struct cosmology *cosmo,
- const integertime_t ti_current) {
+runner_iact_nonsym_bh_gas_feedback(
+ const float r2, const float *dx, const float hi, const float hj,
+ const struct bpart *restrict bi, struct part *restrict pj,
+ struct xpart *restrict xpj, const struct cosmology *cosmo,
+ const struct gravity_props *grav_props, const integertime_t ti_current) {
/* Get the heating probability */
const float prob = bi->to_distribute.AGN_heating_probability;
diff --git a/src/black_holes/EAGLE/black_holes_io.h b/src/black_holes/EAGLE/black_holes_io.h
index ea8b6f70d2847e65128cb23051ee059a1460b52d..cd6c4599f93b3a324ed2be16e6e3fb4ee06ec02c 100644
--- a/src/black_holes/EAGLE/black_holes_io.h
+++ b/src/black_holes/EAGLE/black_holes_io.h
@@ -104,6 +104,7 @@ INLINE static void convert_bpart_vel(const struct engine* e,
* @param bparts The b-particle array.
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
+ * @param with_cosmology Are we running a cosmological simulation?
*/
INLINE static void black_holes_write_particles(const struct bpart* bparts,
struct io_props* list,
diff --git a/src/black_holes/EAGLE/black_holes_part.h b/src/black_holes/EAGLE/black_holes_part.h
index 39ed9d12849c4bd8a4bff62b1ef4853d858ebaac..4d8a9c8e177cc64d7e3191c046ff95c374b2d4a4 100644
--- a/src/black_holes/EAGLE/black_holes_part.h
+++ b/src/black_holes/EAGLE/black_holes_part.h
@@ -122,6 +122,16 @@ struct bpart {
} to_distribute;
+ struct {
+
+ /*! Value of the minimum potential across all neighbours. */
+ float min_potential;
+
+ /*! New position of the BH following the reposition procedure */
+ double x[3];
+
+ } reposition;
+
/*! Chemistry information (e.g. metal content at birth, swallowed metal
* content, etc.) */
struct chemistry_bpart_data chemistry_data;
diff --git a/src/black_holes/EAGLE/black_holes_properties.h b/src/black_holes/EAGLE/black_holes_properties.h
index 50981b43e1ec23927d9958c95ec396fe4a4e4d3a..db50e5b3e08b1a66497aaf75da9d3acdea3287ee 100644
--- a/src/black_holes/EAGLE/black_holes_properties.h
+++ b/src/black_holes/EAGLE/black_holes_properties.h
@@ -74,6 +74,11 @@ struct black_holes_props {
/*! Number of gas neighbours to heat in a feedback event */
float num_ngbs_to_heat;
+ /* ---- Properties of the repositioning model --- */
+
+ /*! Maximal mass of BH to reposition */
+ float max_reposition_mass;
+
/* ---- Common conversion factors --------------- */
/*! Conversion factor from temperature to internal energy */
@@ -157,6 +162,14 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
bp->num_ngbs_to_heat =
parser_get_param_float(params, "EAGLEAGN:AGN_num_ngb_to_heat");
+ /* Reposition parameters --------------------------------- */
+
+ bp->max_reposition_mass =
+ parser_get_param_float(params, "EAGLEAGN:max_reposition_mass");
+
+ /* Convert to internal units */
+ bp->max_reposition_mass *= phys_const->const_solar_mass;
+
/* Common conversion factors ----------------------------- */
/* Calculate temperature to internal energy conversion factor (all internal
diff --git a/src/cell.c b/src/cell.c
index 0b949f3a27edf6e35e4c7679a9daf06f3829397b..3aa95a72c5e3c1af29b3f1180ea9f37d8e219020 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -4419,7 +4419,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, e->star_formation);
+ star_formation_init_part(p, xp, e->star_formation);
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/collectgroup.c b/src/collectgroup.c
index d17608d67a291ae412728b3fcd9ea80c192802bf..8dd9e974f41c979f4c50e84fb555f9b80db19e25 100644
--- a/src/collectgroup.c
+++ b/src/collectgroup.c
@@ -122,7 +122,8 @@ void collectgroup1_apply(const struct collectgroup1 *grp1, struct engine *e) {
e->total_nr_cells = grp1->total_nr_cells;
e->total_nr_tasks = grp1->total_nr_tasks;
e->tasks_per_cell_max = grp1->tasks_per_cell_max;
- e->sfh = grp1->sfh;
+
+ star_formation_logger_add_to_accumulator(&e->sfh, &grp1->sfh);
}
/**
diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c
index 3fda0c80040b6ad139b9100064b937e695e6030f..674631ae52329611eeddb490df1ed205333c0322 100644
--- a/src/cooling/EAGLE/cooling.c
+++ b/src/cooling/EAGLE/cooling.c
@@ -207,7 +207,7 @@ INLINE static double bisection_iter(
const float d_He, const double Lambda_He_reion_cgs,
const double ratefact_cgs,
const struct cooling_function_data *restrict cooling,
- const float abundance_ratio[chemistry_element_count + 2],
+ const float abundance_ratio[eagle_cooling_N_abundances],
const double dt_cgs, const long long ID) {
/* Bracketing */
@@ -419,7 +419,7 @@ void cooling_cool_part(const struct phys_const *phys_const,
* Note that we need to add S and Ca that are in the tables but not tracked
* by the particles themselves.
* The order is [H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe] */
- float abundance_ratio[chemistry_element_count + 2];
+ float abundance_ratio[eagle_cooling_N_abundances];
abundance_ratio_to_solar(p, cooling, abundance_ratio);
/* Get the Hydrogen and Helium mass fractions */
diff --git a/src/cooling/EAGLE/cooling_rates.h b/src/cooling/EAGLE/cooling_rates.h
index 6707a9a9d70b58eb218cfc8ebeca9b1f08fdcb12..bf908c3a7cfc6008995b0d60ab911ceedc2f60dc 100644
--- a/src/cooling/EAGLE/cooling_rates.h
+++ b/src/cooling/EAGLE/cooling_rates.h
@@ -43,16 +43,13 @@
* We also re-order the elements such that they match the order of the
* tables. This is [H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe].
*
- * The solar abundances table (from the cooling struct) is arranged as
- * [H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe].
- *
* @param p Pointer to #part struct.
* @param cooling #cooling_function_data struct.
* @param ratio_solar (return) Array of ratios to solar abundances.
*/
__attribute__((always_inline)) INLINE void abundance_ratio_to_solar(
const struct part *p, const struct cooling_function_data *cooling,
- float ratio_solar[chemistry_element_count + 2]) {
+ float ratio_solar[eagle_cooling_N_abundances]) {
ratio_solar[0] =
p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_H] *
@@ -313,7 +310,7 @@ __attribute__((always_inline)) INLINE double eagle_Compton_cooling_rate(
*/
INLINE static double eagle_metal_cooling_rate(
const double log10_u_cgs, const double redshift, const double n_H_cgs,
- const float solar_ratio[chemistry_element_count + 2], const int n_H_index,
+ const float solar_ratio[eagle_cooling_N_abundances], const int n_H_index,
const float d_n_H, const int He_index, const float d_He,
const struct cooling_function_data *cooling, double *element_lambda) {
@@ -537,7 +534,7 @@ INLINE static double eagle_metal_cooling_rate(
*/
INLINE static double eagle_cooling_rate(
const double log10_u_cgs, const double redshift, const double n_H_cgs,
- const float abundance_ratio[chemistry_element_count + 2],
+ const float abundance_ratio[eagle_cooling_N_abundances],
const int n_H_index, const float d_n_H, const int He_index,
const float d_He, const struct cooling_function_data *cooling) {
diff --git a/src/cooling/EAGLE/cooling_tables.c b/src/cooling/EAGLE/cooling_tables.c
index cddc8d50e0c2f00c30846ebaf65b2bf250e9cc8a..4261e9ac0a6fee9f77c03afe22b7a9b66ade487d 100644
--- a/src/cooling/EAGLE/cooling_tables.c
+++ b/src/cooling/EAGLE/cooling_tables.c
@@ -213,10 +213,6 @@ void read_cooling_header(const char *fname,
if (N_SolarAbundances != eagle_cooling_N_abundances)
error("Invalid solar abundances array length.");
- /* Check value */
- if (N_SolarAbundances != chemistry_element_count + 2)
- error("Number of abundances not compatible with the chemistry model.");
-
dataset = H5Dopen(tempfile_id, "/Header/Number_of_metals", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
&N_Elements);
diff --git a/src/cooling/EAGLE/newton_cooling.c b/src/cooling/EAGLE/newton_cooling.c
index c68a77614425f85da011ab45d2ec488710e068dc..e86ebb7a2f38c29b2d889fb96c4a7d5f791b699c 100644
--- a/src/cooling/EAGLE/newton_cooling.c
+++ b/src/cooling/EAGLE/newton_cooling.c
@@ -214,7 +214,7 @@ __attribute__((always_inline)) INLINE double eagle_convert_u_to_temp(
*/
INLINE static double eagle_metal_cooling_rate(
double log10_u_cgs, double redshift, double n_H_cgs,
- const float solar_ratio[chemistry_element_count + 2], int n_H_index,
+ const float solar_ratio[eagle_cooling_N_abundances], int n_H_index,
float d_n_H, int He_index, float d_He,
const struct cooling_function_data *restrict cooling, double *dlambda_du,
double *element_lambda) {
@@ -569,7 +569,7 @@ INLINE static double eagle_metal_cooling_rate(
*/
INLINE static double eagle_cooling_rate(
double log_u_cgs, double redshift, double n_H_cgs,
- const float abundance_ratio[chemistry_element_count + 2], int n_H_index,
+ const float abundance_ratio[eagle_cooling_N_abundances], int n_H_index,
float d_n_H, int He_index, float d_He,
const struct cooling_function_data *restrict cooling,
double *dLambdaNet_du) {
@@ -609,7 +609,7 @@ INLINE static double eagle_metal_cooling_rate(
const struct cosmology *restrict cosmo,
const struct cooling_function_data *restrict cooling,
const struct phys_const *restrict phys_const,
- const float abundance_ratio[chemistry_element_count + 2],
+ const float abundance_ratio[eagle_cooling_N_abundances],
float dt, int *bisection_flag) {
double logu, logu_old;
diff --git a/src/cooling/grackle/cooling.h b/src/cooling/grackle/cooling.h
index 1abbe5d3827726bde34502cd2b6a1d18cd309950..98d34dc630aa0c83198ddd60b6ee9e2775778d0f 100644
--- a/src/cooling/grackle/cooling.h
+++ b/src/cooling/grackle/cooling.h
@@ -52,7 +52,7 @@ static gr_float cooling_time(
const struct unit_system* restrict us,
const struct cosmology* restrict cosmo,
const struct cooling_function_data* restrict cooling,
- const struct part* restrict p, const struct xpart* restrict xp);
+ const struct part* restrict p, struct xpart* restrict xp);
static gr_float cooling_new_energy(
const struct phys_const* restrict phys_const,
const struct unit_system* restrict us,
@@ -573,7 +573,7 @@ __attribute__((always_inline)) INLINE static gr_float cooling_time(
const struct unit_system* restrict us,
const struct cosmology* restrict cosmo,
const struct cooling_function_data* restrict cooling,
- const struct part* restrict p, const struct xpart* restrict xp) {
+ const struct part* restrict p, struct xpart* restrict xp) {
/* set current time */
code_units units = cooling->units;
@@ -697,6 +697,17 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
xp->cooling_data.radiated_energy -= hydro_get_mass(p) * cooling_du_dt * dt;
}
+/**
+ * @brief Compute the temperature of a #part based on the cooling function.
+ *
+ * @param phys_const #phys_const data structure.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param us The internal system of units.
+ * @param cosmo #cosmology data structure.
+ * @param cooling #cooling_function_data struct.
+ * @param p #part data.
+ * @param xp Pointer to the #xpart data.
+ */
static INLINE float cooling_get_temperature(
const struct phys_const* restrict phys_const,
const struct hydro_props* restrict hydro_props,
@@ -704,9 +715,30 @@ static INLINE float cooling_get_temperature(
const struct cosmology* restrict cosmo,
const struct cooling_function_data* restrict cooling,
const struct part* restrict p, const struct xpart* restrict xp) {
+ // TODO use the grackle library
- error("This function needs implementing!!");
- return 0.;
+ /* Physical constants */
+ const double m_H = phys_const->const_proton_mass;
+ const double k_B = phys_const->const_boltzmann_k;
+
+ /* Gas properties */
+ const double T_transition = hydro_props->hydrogen_ionization_temperature;
+ const double mu_neutral = hydro_props->mu_neutral;
+ const double mu_ionised = hydro_props->mu_ionised;
+
+ /* Particle temperature */
+ const double u = hydro_get_physical_internal_energy(p, xp, cosmo);
+
+ /* Temperature over mean molecular weight */
+ const double T_over_mu = hydro_gamma_minus_one * u * m_H / k_B;
+
+ /* Are we above or below the HII -> HI transition? */
+ if (T_over_mu > (T_transition + 1.) / mu_ionised)
+ return T_over_mu * mu_ionised;
+ else if (T_over_mu < (T_transition - 1.) / mu_neutral)
+ return T_over_mu * mu_neutral;
+ else
+ return T_transition;
}
/**
diff --git a/src/engine.c b/src/engine.c
index e9c250939c20908392086cb811351343aab1c26c..75ff70b693c121770e7b1331ec3e7e1b87aecbbe 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -5069,6 +5069,11 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
parser_get_opt_param_double(params, "FOF:delta_time", -1.);
}
+ /* Initialize the star formation history structure */
+ if (e->policy & engine_policy_star_formation) {
+ star_formation_logger_accumulator_init(&e->sfh);
+ }
+
engine_init_output_lists(e, params);
}
diff --git a/src/engine.h b/src/engine.h
index 8296933561b5f7664904d439ec6ca20f1e37a4b8..967a430871648ac4cb91f390f8b85675a314d3a8 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -237,7 +237,7 @@ struct engine {
long long b_updates_since_rebuild;
/* Star formation logger information */
- struct star_formation_history sfh;
+ struct star_formation_history_accumulator sfh;
/* Properties of the previous step */
int step_props;
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index afc9c17d5b6ed691685bf42df2562ecb6f4409d8..93c0e58068341ab29a7ed0fb6b0d21f76980d952 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -41,6 +41,7 @@
#include "hydro_space.h"
#include "kernel_hydro.h"
#include "minmax.h"
+#include "pressure_floor.h"
#include "./hydro_parameters.h"
@@ -109,7 +110,8 @@ hydro_get_drifted_physical_internal_energy(const struct part *restrict p,
__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
const struct part *restrict p) {
- return gas_pressure_from_entropy(p->rho, p->entropy);
+ const float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+ return pressure_floor_get_pressure(p, p->rho, comoving_pressure);
}
/**
@@ -121,7 +123,10 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
__attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
const struct part *restrict p, const struct cosmology *cosmo) {
- return gas_pressure_from_entropy(p->rho * cosmo->a3_inv, p->entropy);
+ const float phys_pressure =
+ gas_pressure_from_entropy(p->rho * cosmo->a3_inv, p->entropy);
+ const float phys_rho = hydro_get_physical_density(p, cosmo);
+ return pressure_floor_get_pressure(p, phys_rho, phys_pressure);
}
/**
@@ -388,13 +393,15 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
const float rho_inv = 1.f / p->rho;
/* Compute the pressure */
- const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+ float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+ comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure);
/* Compute the sound speed */
- const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+ const float soundspeed =
+ gas_soundspeed_from_pressure(p->rho, comoving_pressure);
/* Divide the pressure by the density squared to get the SPH term */
- const float P_over_rho2 = pressure * rho_inv * rho_inv;
+ const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv;
/* Update variables. */
p->force.P_over_rho2 = P_over_rho2;
@@ -591,13 +598,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const float abs_div_physical_v = fabsf(div_physical_v);
/* Compute the pressure */
- const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+ float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+ comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure);
/* Compute the sound speed */
- const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+ const float soundspeed =
+ gas_soundspeed_from_pressure(p->rho, comoving_pressure);
/* Divide the pressure by the density squared to get the SPH term */
- const float P_over_rho2 = pressure * rho_inv * rho_inv;
+ const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv;
/* Compute the Balsara switch */
/* Pre-multiply in the AV factor; hydro_props are not passed to the iact
@@ -665,14 +674,16 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
p->entropy = xp->entropy_full;
/* Re-compute the pressure */
- const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+ float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+ comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure);
/* Compute the new sound speed */
- const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+ const float soundspeed =
+ gas_soundspeed_from_pressure(p->rho, comoving_pressure);
/* Divide the pressure by the density squared to get the SPH term */
const float rho_inv = 1.f / p->rho;
- const float P_over_rho2 = pressure * rho_inv * rho_inv;
+ const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv;
/* Update variables */
p->force.soundspeed = soundspeed;
@@ -730,14 +741,16 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->rho *= expf(w2);
/* Re-compute the pressure */
- const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+ float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+ comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure);
/* Compute the new sound speed */
- const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+ const float soundspeed =
+ gas_soundspeed_from_pressure(p->rho, comoving_pressure);
/* Divide the pressure by the density squared to get the SPH term */
const float rho_inv = 1.f / p->rho;
- const float P_over_rho2 = pressure * rho_inv * rho_inv;
+ const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv;
/* Update variables */
p->force.soundspeed = soundspeed;
@@ -840,14 +853,16 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
}
/* Compute the pressure */
- const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+ float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+ comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure);
/* Compute the sound speed */
- const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+ const float soundspeed =
+ gas_soundspeed_from_pressure(p->rho, comoving_pressure);
/* Divide the pressure by the density squared to get the SPH term */
const float rho_inv = 1.f / p->rho;
- const float P_over_rho2 = pressure * rho_inv * rho_inv;
+ const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv;
p->force.soundspeed = soundspeed;
p->force.P_over_rho2 = P_over_rho2;
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index 7a8844d560561dae80c676a0b5bb72b34416d080..853d2adf17bc069434562fa96ddb881f760f6830 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/src/hydro/Gadget2/hydro_part.h
@@ -155,6 +155,9 @@ 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;
+
/* Time-step length */
timebin_t time_bin;
diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h
index e3a4f98160150717e50651760d8af68af2a70662..15f9a958a6c5d4fc199a1c1d7d85d3f909520d23 100644
--- a/src/hydro/PressureEnergy/hydro.h
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -46,6 +46,7 @@
#include "hydro_space.h"
#include "kernel_hydro.h"
#include "minmax.h"
+#include "pressure_floor.h"
#include "./hydro_parameters.h"
@@ -221,7 +222,9 @@ hydro_get_comoving_soundspeed(const struct part *restrict p) {
/* Compute the sound speed -- see theory section for justification */
/* IDEAL GAS ONLY -- P-U does not work with generic EoS. */
- const float square_rooted = sqrtf(hydro_gamma * p->pressure_bar / p->rho);
+ const float comoving_pressure =
+ pressure_floor_get_pressure(p, p->rho, p->pressure_bar);
+ const float square_rooted = sqrtf(hydro_gamma * comoving_pressure / p->rho);
return square_rooted;
}
@@ -236,7 +239,10 @@ __attribute__((always_inline)) INLINE static float
hydro_get_physical_soundspeed(const struct part *restrict p,
const struct cosmology *cosmo) {
- return cosmo->a_factor_sound_speed * p->force.soundspeed;
+ const float phys_rho = hydro_get_physical_density(p, cosmo);
+
+ return pressure_floor_get_pressure(
+ p, phys_rho, cosmo->a_factor_sound_speed * p->force.soundspeed);
}
/**
@@ -640,10 +646,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
hydro_one_over_gamma_minus_one) /
(1.f + common_factor * p->density.wcount_dh);
+ /* Get the pressures */
+ const float comoving_pressure_with_floor =
+ pressure_floor_get_pressure(p, p->rho, p->pressure_bar);
+
/* Update variables. */
p->force.f = grad_h_term;
p->force.soundspeed = soundspeed;
p->force.balsara = balsara;
+ p->force.pressure_bar_with_floor = comoving_pressure_with_floor;
}
/**
@@ -755,6 +766,11 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
const float soundspeed = hydro_get_comoving_soundspeed(p);
p->force.soundspeed = soundspeed;
+
+ /* update the required variables */
+ const float comoving_pressure_with_floor =
+ pressure_floor_get_pressure(p, p->rho, p->pressure_bar);
+ p->force.pressure_bar_with_floor = comoving_pressure_with_floor;
}
/**
diff --git a/src/hydro/PressureEnergy/hydro_iact.h b/src/hydro/PressureEnergy/hydro_iact.h
index b59fdbbe418d0442fd290b3596fd08531fae88e6..aab8237c7b26c50d9e04610c6c1029f97e5d73ed 100644
--- a/src/hydro/PressureEnergy/hydro_iact.h
+++ b/src/hydro/PressureEnergy/hydro_iact.h
@@ -259,11 +259,18 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Convolve with the kernel */
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
+ /* Compute the ratio of pressures */
+ const float pressure_inverse_i =
+ pi->force.pressure_bar_with_floor / (pi->pressure_bar * pi->pressure_bar);
+ const float pressure_inverse_j =
+ pj->force.pressure_bar_with_floor / (pj->pressure_bar * pj->pressure_bar);
+
/* SPH acceleration term */
- const float sph_acc_term =
- pj->u * pi->u * hydro_gamma_minus_one * hydro_gamma_minus_one *
- ((f_ij / pi->pressure_bar) * wi_dr + (f_ji / pj->pressure_bar) * wj_dr) *
- r_inv;
+ const float sph_acc_term = pj->u * pi->u * hydro_gamma_minus_one *
+ hydro_gamma_minus_one *
+ ((f_ij * pressure_inverse_i) * wi_dr +
+ (f_ji * pressure_inverse_j) * wj_dr) *
+ r_inv;
/* Assemble the acceleration */
const float acc = sph_acc_term + visc_acc_term;
@@ -278,11 +285,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->a_hydro[2] += mi * acc * dx[2];
/* Get the time derivative for u. */
+
const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one *
- pj->u * pi->u * (f_ij / pi->pressure_bar) *
+ pj->u * pi->u * (f_ij * pressure_inverse_i) *
wi_dr * dvdr * r_inv;
+
const float sph_du_term_j = hydro_gamma_minus_one * hydro_gamma_minus_one *
- pi->u * pj->u * (f_ji / pj->pressure_bar) *
+ pi->u * pj->u * (f_ji * pressure_inverse_j) *
wj_dr * dvdr * r_inv;
/* Viscosity term */
@@ -386,11 +395,18 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Convolve with the kernel */
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
+ /* Compute the ratio of pressures */
+ const float pressure_inverse_i =
+ pi->force.pressure_bar_with_floor / (pi->pressure_bar * pi->pressure_bar);
+ const float pressure_inverse_j =
+ pj->force.pressure_bar_with_floor / (pj->pressure_bar * pj->pressure_bar);
+
/* SPH acceleration term */
- const float sph_acc_term =
- pj->u * pi->u * hydro_gamma_minus_one * hydro_gamma_minus_one *
- ((f_ij / pi->pressure_bar) * wi_dr + (f_ji / pj->pressure_bar) * wj_dr) *
- r_inv;
+ const float sph_acc_term = pj->u * pi->u * hydro_gamma_minus_one *
+ hydro_gamma_minus_one *
+ ((f_ij * pressure_inverse_i) * wi_dr +
+ (f_ji * pressure_inverse_j) * wj_dr) *
+ r_inv;
/* Assemble the acceleration */
const float acc = sph_acc_term + visc_acc_term;
@@ -402,7 +418,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Get the time derivative for u. */
const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one *
- pj->u * pi->u * (f_ij / pi->pressure_bar) *
+ pj->u * pi->u * (f_ij * pressure_inverse_i) *
wi_dr * dvdr * r_inv;
/* Viscosity term */
diff --git a/src/hydro/PressureEnergy/hydro_part.h b/src/hydro/PressureEnergy/hydro_part.h
index ded83d70176409332b669517ed203fdf95ecfd5e..c8f02f518cfde536965c0e6d000999a0e07e4aab 100644
--- a/src/hydro/PressureEnergy/hydro_part.h
+++ b/src/hydro/PressureEnergy/hydro_part.h
@@ -168,6 +168,9 @@ struct part {
/*! Balsara switch */
float balsara;
+
+ /*! Pressure term including the pressure floor */
+ float pressure_bar_with_floor;
} force;
};
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index b24aff140b4b9cb53af0c09a136a2aa23bf4938c..d28cc1dff363517e72ebb3607ec5aa8121f828b2 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -34,6 +34,7 @@
#include "hydro.h"
#include "kernel_hydro.h"
#include "parser.h"
+#include "pressure_floor.h"
#include "units.h"
#define hydro_props_default_max_iterations 30
@@ -186,6 +187,9 @@ void hydro_props_print(const struct hydro_props *p) {
/* Print equation of state first */
eos_print(&eos);
+ /* Then the pressure floor */
+ pressure_floor_print(&pressure_floor_props);
+
/* Now SPH */
message("Hydrodynamic scheme: %s in %dD.", SPH_IMPLEMENTATION,
(int)hydro_dimension);
@@ -238,6 +242,7 @@ void hydro_props_print(const struct hydro_props *p) {
void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
eos_print_snapshot(h_grpsph, &eos);
+ pressure_floor_print_snapshot(h_grpsph);
io_write_attribute_i(h_grpsph, "Dimension", (int)hydro_dimension);
io_write_attribute_s(h_grpsph, "Scheme", SPH_IMPLEMENTATION);
diff --git a/src/pressure_floor.c b/src/pressure_floor.c
new file mode 100644
index 0000000000000000000000000000000000000000..2901241f1a2ff42d526a0122008dcf6ab2d8ea6f
--- /dev/null
+++ b/src/pressure_floor.c
@@ -0,0 +1,24 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* This object's header. */
+#include "pressure_floor.h"
+
+/* Pressure floor for the physics model. */
+struct pressure_floor_properties pressure_floor_props;
diff --git a/src/pressure_floor.h b/src/pressure_floor.h
new file mode 100644
index 0000000000000000000000000000000000000000..4389dfe0891dd6bfbc1e6730cac1c6ddcc920547
--- /dev/null
+++ b/src/pressure_floor.h
@@ -0,0 +1,54 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_PRESSURE_FLOOR_H
+#define SWIFT_PRESSURE_FLOOR_H
+
+/**
+ * @file src/pressure_floor.h
+ * @brief Branches between the different pressure floor models
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local includes */
+#include "common_io.h"
+#include "cosmology.h"
+#include "error.h"
+#include "inline.h"
+
+extern struct pressure_floor_properties pressure_floor_props;
+
+/* Check if pressure floor is implemented in hydro */
+#ifndef PRESSURE_FLOOR_NONE
+#if defined(GADGET2_SPH) || defined(HOPKINS_PU_SPH)
+/* Implemented */
+#else
+#error Pressure floor not implemented with this hydro scheme
+#endif
+
+#endif
+/* Import the right pressure floor definition */
+#if defined(PRESSURE_FLOOR_NONE)
+#include "./pressure_floor/none/pressure_floor.h"
+#elif defined(PRESSURE_FLOOR_GEAR)
+#include "./pressure_floor/GEAR/pressure_floor.h"
+#endif
+
+#endif /* SWIFT_PRESSURE_FLOOR_H */
diff --git a/src/pressure_floor/GEAR/pressure_floor.h b/src/pressure_floor/GEAR/pressure_floor.h
new file mode 100644
index 0000000000000000000000000000000000000000..dd715c155a095f9ad5f97b1d14cabfc94d9b11b0
--- /dev/null
+++ b/src/pressure_floor/GEAR/pressure_floor.h
@@ -0,0 +1,124 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_PRESSURE_FLOOR_GEAR_H
+#define SWIFT_PRESSURE_FLOOR_GEAR_H
+
+#include "adiabatic_index.h"
+#include "cosmology.h"
+#include "equation_of_state.h"
+#include "hydro_properties.h"
+#include "parser.h"
+#include "part.h"
+#include "units.h"
+
+/**
+ * @file src/pressure_floor/GEAR/pressure_floor.h
+ * @brief Pressure floor used in the GEAR model
+ */
+
+/**
+ * @brief Properties of the pressure floor in the GEAR model.
+ */
+struct pressure_floor_properties {
+
+ /*! Jeans factor */
+ float n_jeans;
+
+ /*! The constants in internal units (4 G N_jeans^(2/3) / PI) */
+ float constants;
+};
+
+/**
+ * @brief Compute the physical pressure floor of a given #part.
+ *
+ * Note that the particle is not updated!!
+ *
+ * The inputs can be either in physical or comoving coordinates (the output is
+ * in the same coordinates).
+ *
+ * @param p The #part.
+ * @param rho The physical or comoving density.
+ * @param pressure The physical or comoving pressure without any pressure floor.
+ *
+ * @return The physical or comoving pressure with the floor.
+ */
+static INLINE float pressure_floor_get_pressure(const struct part *p,
+ const float rho,
+ const float pressure) {
+
+ /* Compute pressure floor */
+ float floor = p->h * p->h * rho * pressure_floor_props.constants;
+ // TODO add sigma (will be done once the SF is merged)
+ floor *= rho * hydro_one_over_gamma;
+
+ return fmax(pressure, floor);
+}
+
+/**
+ * @brief Initialise the pressure floor by reading the parameters and converting
+ * to internal units.
+ *
+ * The input temperatures and number densities are converted to pressure and
+ * density assuming a neutral gas of primoridal abundance.
+ *
+ * @param params The YAML parameter file.
+ * @param us The system of units used internally.
+ * @param phys_const The physical constants.
+ * @param hydro_props The propoerties of the hydro scheme.
+ * @param props The pressure floor properties to fill.
+ */
+static INLINE void pressure_floor_init(struct pressure_floor_properties *props,
+ const struct phys_const *phys_const,
+ const struct unit_system *us,
+ const struct hydro_props *hydro_props,
+ struct swift_params *params) {
+
+ /* Read the Jeans factor */
+ props->n_jeans =
+ parser_get_param_float(params, "GEARPressureFloor:Jeans_factor");
+
+ /* Compute the constants */
+ props->constants =
+ 4.0 * M_1_PI * phys_const->const_newton_G * pow(props->n_jeans, 2. / 3.);
+}
+
+/**
+ * @brief Print the properties of the pressure floor to stdout.
+ *
+ * @param props The pressure floor properties.
+ */
+static INLINE void pressure_floor_print(
+ const struct pressure_floor_properties *props) {
+
+ message("Pressure floor is 'GEAR' with:");
+ message("Jeans factor: %g", props->n_jeans);
+}
+
+#ifdef HAVE_HDF5
+
+/**
+ * @brief Writes the current model of pressure floor to the file
+ * @param h_grp The HDF5 group in which to write
+ */
+INLINE static void pressure_floor_print_snapshot(hid_t h_grp) {
+
+ io_write_attribute_s(h_grp, "Pressure floor", "GEAR");
+}
+#endif
+#endif /* SWIFT_PRESSURE_FLOOR_GEAR_H */
diff --git a/src/pressure_floor/none/pressure_floor.h b/src/pressure_floor/none/pressure_floor.h
new file mode 100644
index 0000000000000000000000000000000000000000..2d0b95a71f9ccce6697e79760aa43b560933e7bd
--- /dev/null
+++ b/src/pressure_floor/none/pressure_floor.h
@@ -0,0 +1,99 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_PRESSURE_FLOOR_NONE_H
+#define SWIFT_PRESSURE_FLOOR_NONE_H
+
+#include "adiabatic_index.h"
+#include "cosmology.h"
+#include "equation_of_state.h"
+#include "hydro_properties.h"
+#include "parser.h"
+#include "part.h"
+#include "units.h"
+
+/**
+ * @file src/pressure_floor/none/pressure_floor.h
+ * @brief Model without pressure floor
+ */
+
+/**
+ * @brief Properties of the pressure floor in the NONE model.
+ */
+struct pressure_floor_properties {};
+
+/**
+ * @brief Compute the physical pressure floor of a given #part.
+ *
+ * Note that the particle is not updated!!
+ *
+ * The inputs can be either in physical or comoving coordinates (the output is
+ * in the same coordinates).
+ *
+ * @param p The #part.
+ * @param rho The physical or comoving density.
+ * @param pressure The physical or comoving pressure without any pressure floor.
+ *
+ * @return The physical or comoving pressure with the floor.
+ */
+static INLINE float pressure_floor_get_pressure(const struct part *p,
+ const float rho,
+ const float pressure) {
+ return pressure;
+}
+
+/**
+ * @brief Initialise the pressure floor by reading the parameters and converting
+ * to internal units.
+ *
+ * The input temperatures and number densities are converted to pressure and
+ * density assuming a neutral gas of primoridal abundance.
+ *
+ * @param params The YAML parameter file.
+ * @param us The system of units used internally.
+ * @param phys_const The physical constants.
+ * @param hydro_props The propoerties of the hydro scheme.
+ * @param props The pressure floor properties to fill.
+ */
+static INLINE void pressure_floor_init(struct pressure_floor_properties *props,
+ const struct phys_const *phys_const,
+ const struct unit_system *us,
+ const struct hydro_props *hydro_props,
+ struct swift_params *params) {}
+
+/**
+ * @brief Print the properties of the pressure floor to stdout.
+ *
+ * @param props The pressure floor properties.
+ */
+static INLINE void pressure_floor_print(
+ const struct pressure_floor_properties *props) {}
+
+#ifdef HAVE_HDF5
+
+/**
+ * @brief Writes the current model of pressure floor to the file
+ * @param h_grp The HDF5 group in which to write
+ */
+INLINE static void pressure_floor_print_snapshot(hid_t h_grp) {
+
+ io_write_attribute_s(h_grp, "Pressure floor", "none");
+}
+#endif
+
+#endif /* SWIFT_PRESSURE_FLOOR_NONE_H */
diff --git a/src/runner.c b/src/runner.c
index 9249924a92dbcffba30baafa76df1a18e9962c16..7d41019ed3dc00364e7c170e3310d93f303e7414 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -905,6 +905,10 @@ void runner_do_black_holes_swallow_ghost(struct runner *r, struct cell *c,
if (bpart_is_active(bp, e)) {
+ /* Compute the final operations for repositioning of this BH */
+ black_holes_end_reposition(bp, e->black_holes_properties,
+ e->physical_constants, e->cosmology);
+
/* Get particle time-step */
double dt;
if (with_cosmology) {
@@ -2272,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, star_formation);
+ star_formation_init_part(p, xp, star_formation);
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/runner_doiact.h b/src/runner_doiact.h
index 6caa287cf726f85778ef5abdc184acaf759e8b0e..1a39c8d49f2c4234c04982e255705f06ec1c5d38 100644
--- a/src/runner_doiact.h
+++ b/src/runner_doiact.h
@@ -884,7 +884,6 @@ void DOSELF_SUBSET(struct runner *r, struct cell *restrict ci,
const int count_i = ci->hydro.count;
struct part *restrict parts_j = ci->hydro.parts;
-
/* Loop over the parts in ci. */
for (int pid = 0; pid < count; pid++) {
diff --git a/src/runner_doiact_black_holes.h b/src/runner_doiact_black_holes.h
index 71ce7aefcfc3689cc8476884729e4e7cd77bef66..ce159c7ac24a508bc625070ed50b3aad7dd9fa8d 100644
--- a/src/runner_doiact_black_holes.h
+++ b/src/runner_doiact_black_holes.h
@@ -159,7 +159,8 @@ void DOSELF1_BH(struct runner *r, struct cell *c, int timer) {
#endif
if (r2 < hig2) {
- IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, cosmo, ti_current);
+ IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, cosmo, e->gravity_properties,
+ ti_current);
}
} /* loop over the parts in ci. */
} /* loop over the bparts in ci. */
@@ -310,7 +311,8 @@ void DO_NONSYM_PAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
#endif
if (r2 < hig2) {
- IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, cosmo, ti_current);
+ IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, cosmo, e->gravity_properties,
+ ti_current);
}
} /* loop over the parts in cj. */
} /* loop over the bparts in ci. */
@@ -475,7 +477,8 @@ void DOPAIR1_SUBSET_BH_NAIVE(struct runner *r, struct cell *restrict ci,
#endif
/* Hit or miss? */
if (r2 < hig2) {
- IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, cosmo, ti_current);
+ IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, cosmo, e->gravity_properties,
+ ti_current);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
@@ -550,7 +553,8 @@ void DOSELF1_SUBSET_BH(struct runner *r, struct cell *restrict ci,
/* Hit or miss? */
if (r2 < hig2) {
- IACT_BH_GAS(r2, dx, hi, pj->h, bi, pj, xpj, cosmo, ti_current);
+ IACT_BH_GAS(r2, dx, hi, pj->h, bi, pj, xpj, cosmo,
+ e->gravity_properties, ti_current);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
diff --git a/src/space.c b/src/space.c
index cea579d09f438ff07611eb4c674decd7a3547de8..2fd83785451b634b97610bea34fb0bce032f9acc 100644
--- a/src/space.c
+++ b/src/space.c
@@ -4373,7 +4373,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], e->star_formation);
+ star_formation_init_part(&parts[k], &xparts[k], e->star_formation);
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/space.h b/src/space.h
index 94ba0654352c6c33f1115337073bb76d099a0354..0b332716645e733636b7ab0da57a0a31b28e3d31 100644
--- a/src/space.h
+++ b/src/space.h
@@ -115,9 +115,6 @@ struct space {
/*! The minimum top-level cell width allowed. */
double cell_min;
- /*! Current maximum displacement for particles. */
- float dx_max;
-
/*! Space dimensions in number of top-cells. */
int cdim[3];
diff --git a/src/star_formation.c b/src/star_formation.c
index 698a64cc636dd79f00feac3f6cc88bf519fe09c1..60cff1e2e68feaf7e71705b5079294ec478fad42 100644
--- a/src/star_formation.c
+++ b/src/star_formation.c
@@ -24,6 +24,7 @@
#include "part.h"
#include "restart.h"
#include "star_formation.h"
+#include "star_formation_io.h"
#include "units.h"
/**
diff --git a/src/star_formation/EAGLE/star_formation.h b/src/star_formation/EAGLE/star_formation.h
index 694cac1512d91ef3d640185cbf9b614773263a75..f2c77e036842b4ca040c58a6bcad1513b03a42bd 100644
--- a/src/star_formation/EAGLE/star_formation.h
+++ b/src/star_formation/EAGLE/star_formation.h
@@ -691,6 +691,7 @@ star_formation_first_init_part(const struct phys_const* restrict phys_const,
* @param data The global star_formation information.
*/
__attribute__((always_inline)) INLINE static void star_formation_init_part(
- struct part* restrict p, const struct star_formation* data) {}
+ 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_logger.h b/src/star_formation/EAGLE/star_formation_logger.h
index d634c876e52d45588ed0b93c0afc09731317c037..fea15242e9a83b81cf5c7801c0546a1bbb5a8eae 100644
--- a/src/star_formation/EAGLE/star_formation_logger.h
+++ b/src/star_formation/EAGLE/star_formation_logger.h
@@ -87,6 +87,28 @@ INLINE static void star_formation_logger_add(
sf_update->SFR_inactive += sf_add->SFR_inactive;
}
+/**
+ * @brief add a star formation history struct to the engine star formation
+ * history accumulator struct
+ *
+ * @param sf_add the star formation accumulator struct which we want to add to
+ * the star formation history
+ * @param sf_update the star formation structure which we want to update
+ */
+INLINE static void star_formation_logger_add_to_accumulator(
+ struct star_formation_history_accumulator *sf_update,
+ const struct star_formation_history *sf_add) {
+
+ /* Update the SFH structure */
+ sf_update->new_stellar_mass = sf_add->new_stellar_mass;
+
+ sf_update->SFR_active = sf_add->SFR_active;
+
+ sf_update->SFRdt_active = sf_add->SFRdt_active;
+
+ sf_update->SFR_inactive = sf_add->SFR_inactive;
+}
+
/**
* @brief Initialize the star formation history structure in the #engine
*
@@ -105,6 +127,24 @@ INLINE static void star_formation_logger_init(
sfh->SFR_inactive = 0.f;
}
+/**
+ * @brief Initialize the star formation history structure in the #engine
+ *
+ * @param sfh The pointer to the star formation history structure
+ */
+INLINE static void star_formation_logger_accumulator_init(
+ struct star_formation_history_accumulator *sfh) {
+
+ /* Initialize the collecting SFH structure to zero */
+ sfh->new_stellar_mass = 0.f;
+
+ sfh->SFR_active = 0.f;
+
+ sfh->SFRdt_active = 0.f;
+
+ sfh->SFR_inactive = 0.f;
+}
+
/**
* @brief Write the final SFH to a file
*
@@ -117,7 +157,7 @@ INLINE static void star_formation_logger_init(
*/
INLINE static void star_formation_logger_write_to_log_file(
FILE *fp, const double time, const double a, const double z,
- const struct star_formation_history sf, const int step) {
+ const struct star_formation_history_accumulator sf, const int step) {
/* Calculate the total SFR */
const float totalSFR = sf.SFR_active + sf.SFR_inactive;
diff --git a/src/star_formation/EAGLE/star_formation_logger_struct.h b/src/star_formation/EAGLE/star_formation_logger_struct.h
index 2a23659c4d931735d1b82a6143b3d9f871f7137a..c03a00c97ead46f552350a43574c5bbe7ac6df1b 100644
--- a/src/star_formation/EAGLE/star_formation_logger_struct.h
+++ b/src/star_formation/EAGLE/star_formation_logger_struct.h
@@ -34,4 +34,21 @@ struct star_formation_history {
float SFRdt_active;
};
+/* Starformation history struct for the engine.
+ Allows to integrate in time some values.
+ Nothing to do in EAGLE => copy of star_formation_history */
+struct star_formation_history_accumulator {
+ /*! Total new stellar mass */
+ float new_stellar_mass;
+
+ /*! SFR of all particles */
+ float SFR_inactive;
+
+ /*! SFR of active particles */
+ float SFR_active;
+
+ /*! SFR*dt of active particles */
+ float SFRdt_active;
+};
+
#endif /* SWIFT_EAGLE_STAR_FORMATION_LOGGER_STRUCT_H */
diff --git a/src/star_formation/EAGLE/star_formation_struct.h b/src/star_formation/EAGLE/star_formation_struct.h
index 41247e160a3eddbc9184c59b67cfa2a1d7259a05..8caac49d4b57652c5db9ae93e3789dc690e6d23f 100644
--- a/src/star_formation/EAGLE/star_formation_struct.h
+++ b/src/star_formation/EAGLE/star_formation_struct.h
@@ -29,4 +29,6 @@ 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 c479feb5c66328f9fab8bf62593ca66b6658b79e..ac423a51865609460e870f65b1eeeb266182e2ef 100644
--- a/src/star_formation/GEAR/star_formation.h
+++ b/src/star_formation/GEAR/star_formation.h
@@ -1,6 +1,7 @@
/*******************************************************************************
* This file is part of SWIFT.
- * Copyright (c) 2018 Folkert Nobels (nobels@strw.leidenuniv.nl)
+ * 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
@@ -20,20 +21,24 @@
#define SWIFT_GEAR_STAR_FORMATION_H
/* Local includes */
+#include "cooling.h"
#include "cosmology.h"
+#include "engine.h"
#include "entropy_floor.h"
#include "error.h"
#include "hydro_properties.h"
#include "parser.h"
#include "part.h"
#include "physical_constants.h"
+#include "random.h"
+#include "star_formation_struct.h"
#include "units.h"
/**
* @brief Calculate if the gas has the potential of becoming
* a star.
*
- * No star formation should occur, so return 0.
+ * Use the star formation criterion given by eq. 3 in Revaz & Jablonka 2018.
*
* @param starform the star formation law properties to use.
* @param p the gas particles.
@@ -46,7 +51,7 @@
*
*/
INLINE static int star_formation_is_star_forming(
- const struct part* restrict p, const struct xpart* restrict xp,
+ struct part* restrict p, struct xpart* restrict xp,
const struct star_formation* starform, const struct phys_const* phys_const,
const struct cosmology* cosmo,
const struct hydro_props* restrict hydro_props,
@@ -54,14 +59,43 @@ INLINE static int star_formation_is_star_forming(
const struct cooling_function_data* restrict cooling,
const struct entropy_floor_properties* restrict entropy_floor) {
- return 0;
+ const float temperature = cooling_get_temperature(phys_const, hydro_props, us,
+ cosmo, cooling, p, xp);
+
+ const float temperature_max = starform->maximal_temperature;
+
+ /* Check the temperature criterion */
+ if (temperature > temperature_max) {
+ return 0;
+ }
+
+ /* Get the required variables */
+ const float sigma2 = p->sf_data.sigma2;
+ const float n_jeans_2_3 = starform->n_jeans_2_3;
+
+ const float h = p->h;
+ const float density = hydro_get_physical_density(p, cosmo);
+
+ // TODO use GRACKLE */
+ const float mu = hydro_props->mu_neutral;
+
+ /* Compute the density criterion */
+ const float coef =
+ M_PI_4 / (phys_const->const_newton_G * n_jeans_2_3 * h * h);
+ const float density_criterion =
+ coef * (hydro_gamma * phys_const->const_boltzmann_k * temperature /
+ (mu * phys_const->const_proton_mass) +
+ sigma2);
+
+ /* Check the density criterion */
+ return density > density_criterion;
}
/**
- * @brief Compute the star-formation rate of a given particle and store
- * it into the #xpart.
+ * @brief Compute the star-formation rate of a given particle.
*
- * Nothing to do here.
+ * Nothing to do here. Everything is done in
+ * #star_formation_should_convert_to_star.
*
* @param p #part.
* @param xp the #xpart.
@@ -71,15 +105,15 @@ INLINE static int star_formation_is_star_forming(
* @param dt_star The time-step of this particle.
*/
INLINE static void star_formation_compute_SFR(
- const struct part* restrict p, struct xpart* restrict xp,
+ struct part* restrict p, struct xpart* restrict xp,
const struct star_formation* starform, const struct phys_const* phys_const,
const struct cosmology* cosmo, const double dt_star) {}
/**
* @brief Decides whether a particle should be converted into a
* star or not.
- *
- * No SF should occur, so return 0.
+
+ * Compute the star formation rate from eq. 4 in Revaz & Jablonka 2012.
*
* @param p The #part.
* @param xp The #xpart.
@@ -89,18 +123,38 @@ INLINE static void star_formation_compute_SFR(
* @return 1 if a conversion should be done, 0 otherwise.
*/
INLINE static int star_formation_should_convert_to_star(
- const struct part* p, const struct xpart* xp,
- const struct star_formation* starform, const struct engine* e,
- const double dt_star) {
+ struct part* p, struct xpart* xp, const struct star_formation* starform,
+ const struct engine* e, const double dt_star) {
+
+ const struct phys_const* phys_const = e->physical_constants;
+ const struct cosmology* cosmo = e->cosmology;
+
+ /* Check that we are running a full time step */
+ if (dt_star == 0.) {
+ return 0;
+ }
+
+ /* Get a few variables */
+ const float G = phys_const->const_newton_G;
+ const float density = hydro_get_physical_density(p, cosmo);
+
+ /* Compute the probability */
+ const float inv_free_fall_time =
+ sqrtf(density * 32.f * G * 0.33333333f * M_1_PI);
+ const float prob = 1.f - exp(-starform->star_formation_efficiency *
+ inv_free_fall_time * dt_star);
+
+ /* Roll the dice... */
+ const float random_number =
+ random_unit_interval(p->id, e->ti_current, random_number_star_formation);
- return 0;
+ /* Can we form a star? */
+ return random_number < prob;
}
/**
* @brief Update the SF properties of a particle that is not star forming.
*
- * Nothing to do here.
- *
* @param p The #part.
* @param xp The #xpart.
* @param e The #engine.
@@ -115,8 +169,6 @@ INLINE static void star_formation_update_part_not_SFR(
* @brief Copies the properties of the gas particle over to the
* star particle.
*
- * Nothing to do here.
- *
* @param e The #engine
* @param p the gas particles.
* @param xp the additional properties of the gas particles.
@@ -133,21 +185,33 @@ INLINE static void star_formation_copy_properties(
const struct phys_const* phys_const,
const struct hydro_props* restrict hydro_props,
const struct unit_system* restrict us,
- const struct cooling_function_data* restrict cooling) {}
+ const struct cooling_function_data* restrict cooling) {
-/**
- * @brief initialization of the star formation law
- *
- * @param parameter_file The parsed parameter file
- * @param phys_const Physical constants in internal units
- * @param us The current internal system of units
- * @param starform the star formation law properties to initialize
- *
- */
-INLINE static void starformation_init_backend(
- struct swift_params* parameter_file, const struct phys_const* phys_const,
- const struct unit_system* us, const struct hydro_props* hydro_props,
- const struct star_formation* starform) {}
+ /* Store the current mass */
+ sp->mass = hydro_get_mass(p);
+ sp->birth.mass = sp->mass;
+
+ /* Store either the birth_scale_factor or birth_time depending */
+ if (with_cosmology) {
+ sp->birth_scale_factor = cosmo->a;
+ } else {
+ sp->birth_time = e->time;
+ }
+
+ // TODO copy only metals
+ /* Store the chemistry struct in the star particle */
+ sp->chemistry_data = p->chemistry_data;
+
+ /* Store the tracers data */
+ sp->tracers_data = xp->tracers_data;
+
+ /* Store the birth density in the star particle */
+ sp->birth.density = hydro_get_physical_density(p, cosmo);
+
+ /* Store the birth temperature*/
+ sp->birth.temperature = cooling_get_temperature(phys_const, hydro_props, us,
+ cosmo, cooling, p, xp);
+}
/**
* @brief Prints the used parameters of the star formation law
@@ -156,7 +220,6 @@ INLINE static void starformation_init_backend(
*/
INLINE static void starformation_print_backend(
const struct star_formation* starform) {
-
message("Star formation law is 'GEAR'");
}
@@ -164,12 +227,22 @@ INLINE static void starformation_print_backend(
* @brief Finishes the density calculation.
*
* @param p The particle to act upon
- * @param cd The global star_formation information.
+ * @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* cd,
- const struct cosmology* cosmo) {}
+ 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.
@@ -183,39 +256,46 @@ __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) {}
+ 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.
*
- * Nothing to do here.
- *
+ * @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.
- * @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) {}
+ 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 data The global star_formation information.
- */
-__attribute__((always_inline)) INLINE static void star_formation_init_part(
- struct part* restrict p, const struct star_formation* data) {}
+ /* 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_iact.h b/src/star_formation/GEAR/star_formation_iact.h
index 749b608068650a27cbe4c9a0ca4126d2740337f3..7325b92af2840b317cf1a924a1e509b34bdffba3 100644
--- a/src/star_formation/GEAR/star_formation_iact.h
+++ b/src/star_formation/GEAR/star_formation_iact.h
@@ -1,6 +1,7 @@
/*******************************************************************************
* This file is part of SWIFT.
- * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
+ * 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
@@ -28,6 +29,8 @@
* @brief do star_formation 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.
@@ -39,7 +42,28 @@
*/
__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) {}
+ struct part *restrict pj, float a, float H) {
+
+ float wi;
+ float wj;
+ /* Evaluation of the SPH kernel */
+ kernel_eval(sqrt(r2) / hi, &wi);
+ kernel_eval(sqrt(r2) / hj, &wj);
+
+ /* 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;
+
+ /* Compute the velocity dispersion */
+ pi->sf_data.sigma2 += sigma2 * wi * hydro_get_mass(pj);
+ pj->sf_data.sigma2 += sigma2 * wj * hydro_get_mass(pi);
+}
/**
* @brief do star_formation computation after the runner_iact_density (non
@@ -58,6 +82,23 @@ __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) {}
+ float H) {
+ float wi;
+ /* Evaluation of the SPH kernel */
+ kernel_eval(sqrt(r2) / hi, &wi);
+
+ /* 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;
+
+ /* Compute the velocity dispersion */
+ pi->sf_data.sigma2 += sigma2 * wi * hydro_get_mass(pj);
+}
#endif /* SWIFT_GEAR_STAR_FORMATION_IACT_H */
diff --git a/src/star_formation/GEAR/star_formation_io.h b/src/star_formation/GEAR/star_formation_io.h
index 6ef04c49c4abcd00175aaa164271628a9ff89360..3e021f7844c1deaeca40d7144d6f7b69cb6c2bdb 100644
--- a/src/star_formation/GEAR/star_formation_io.h
+++ b/src/star_formation/GEAR/star_formation_io.h
@@ -1,6 +1,7 @@
/*******************************************************************************
* This file is part of SWIFT.
- * Copyright (c) 2018 Folkert Nobels (nobels@strw.leidenuniv.nl)
+ * 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
@@ -37,8 +38,40 @@
__attribute__((always_inline)) INLINE static int star_formation_write_particles(
const struct part* parts, const struct xpart* xparts,
struct io_props* list) {
-
+ /* Nothing to write here */
return 0;
}
+/**
+ * @brief initialization of the star formation law
+ *
+ * @param parameter_file The parsed parameter file
+ * @param phys_const Physical constants in internal units
+ * @param us The current internal system of units
+ * @param starform the star formation law properties to initialize
+ *
+ */
+INLINE static void starformation_init_backend(
+ struct swift_params* parameter_file, const struct phys_const* phys_const,
+ 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");
+
+ /* Maximum temperature for star formation */
+ starform->maximal_temperature = parser_get_param_double(
+ parameter_file, "GEARStarFormation:maximal_temperature");
+
+ /* Apply unit change */
+ starform->maximal_temperature *=
+ units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
+}
+
#endif /* SWIFT_STAR_FORMATION_GEAR_IO_H */
diff --git a/src/star_formation/GEAR/star_formation_logger.h b/src/star_formation/GEAR/star_formation_logger.h
index 5b9e033d21d3d202f3c289d1dd6a843ba17fa524..84475909c0f524a4f930a48dbcc3b0943719f8b0 100644
--- a/src/star_formation/GEAR/star_formation_logger.h
+++ b/src/star_formation/GEAR/star_formation_logger.h
@@ -1,6 +1,7 @@
/*******************************************************************************
* This file is part of SWIFT.
- * Copyright (c) 2019 Folkert Nobels (nobels@strw.leidenuniv.nl)
+ * 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
@@ -27,16 +28,28 @@
#include "hydro.h"
#include "part.h"
#include "star_formation_logger_struct.h"
+#include "units.h"
/**
- * @brief Update the stellar mass in the current cell after creating
- * the new star particle spart sp
+ * @brief Update the stellar quantities in the current cell after creating
+ * the new star particle spart sp.
*
+ * @param time_step, the current time step of the simulation
* @param sp new created star particle
* @param sf the star_formation_history struct of the current cell
*/
INLINE static void star_formation_logger_log_new_spart(
- struct spart *sp, struct star_formation_history *sf) {}
+ const struct spart *sp, struct star_formation_history *sf) {
+
+ /* Add mass of created sparticle to the total stellar mass in this cell */
+ sf->new_stellar_mass += sp->mass;
+
+ /* Increase the number of stars */
+ sf->number_new_stars += 1;
+
+ /* No need to deal with the integrated quantities, only the engine's one is
+ * updated */
+}
/**
* @brief Initialize the star formation history struct in the case the cell is
@@ -45,41 +58,42 @@ INLINE static void star_formation_logger_log_new_spart(
* @param sf the star_formation_history struct we want to initialize
*/
INLINE static void star_formation_logger_log_inactive_cell(
- struct star_formation_history *sf) {}
+ struct star_formation_history *sf) {
-/**
- * @brief add a star formation history struct to an other star formation history
- * struct
- *
- * @param sf_add the star formation struct which we want to add to the star
- * formation history
- * @param sf_update the star formation structure which we want to update
- */
-INLINE static void star_formation_logger_add(
- struct star_formation_history *sf_update,
- const struct star_formation_history *sf_add) {}
+ /* Initialize the stellar mass to zero */
+ sf->new_stellar_mass = 0.f;
+ /* initialize number of stars to zero*/
+ sf->number_new_stars = 0;
+}
/**
* @brief Initialize the star formation history structure in the #engine
*
* @param sfh The pointer to the star formation history structure
*/
INLINE static void star_formation_logger_init(
- struct star_formation_history *sfh) {}
+ struct star_formation_history *sfh) {
+ /* Initialize the collecting SFH structure to zero */
+ sfh->new_stellar_mass = 0.f;
+ sfh->number_new_stars = 0;
+}
/**
- * @brief Write the final SFH to a file
+ * @brief add a star formation history struct to an other star formation history
+ * struct
*
- * @param fp The file to write to.
- * @param time the simulation time (time since Big Bang) in internal units.
- * @param a the scale factor.
- * @param z the redshift.
- * @param sf the #star_formation_history struct.
- * @param step The time-step of the simulation.
+ * @param sf_add the star formation struct which we want to add to the star
+ * formation history
+ * @param sf_update the star formation structure which we want to update
*/
-INLINE static void star_formation_logger_write_to_log_file(
- FILE *fp, const double time, const double a, const double z,
- const struct star_formation_history sf, const int step) {}
+INLINE static void star_formation_logger_add(
+ struct star_formation_history *sf_update,
+ const struct star_formation_history *sf_add) {
+
+ /* Update the SFH structure */
+ sf_update->number_new_stars += sf_add->number_new_stars;
+ sf_update->new_stellar_mass += sf_add->new_stellar_mass;
+}
/**
* @brief Initialize the SFH logger file
@@ -90,13 +104,47 @@ INLINE static void star_formation_logger_write_to_log_file(
*/
INLINE static void star_formation_logger_init_log_file(
FILE *fp, const struct unit_system *restrict us,
- const struct phys_const *phys_const) {}
+ const struct phys_const *phys_const) {
+
+ /* Write some general text to the logger file */
+ fprintf(fp, "# Star Formation History Logger file\n");
+ fprintf(fp, "######################################################\n");
+ fprintf(fp, "# The quantities are all given in internal physical units!\n");
+ fprintf(fp, "#\n");
+ fprintf(fp, "# (0) Simulation step\n");
+ fprintf(fp,
+ "# (1) Time since Big Bang (cosmological run), Time since start of "
+ "the simulation (non-cosmological run).\n");
+ fprintf(fp, "# Unit = %e seconds\n", us->UnitTime_in_cgs);
+ fprintf(fp, "# Unit = %e yr or %e Myr\n", 1.f / phys_const->const_year,
+ 1.f / phys_const->const_year / 1e6);
+ fprintf(fp, "# (2) Scale factor (no unit)\n");
+ fprintf(fp, "# (3) Redshift (no unit)\n");
+ fprintf(fp,
+ "# (4) Total number of stars formed in the simulation (no unit)\n");
+ fprintf(fp, "# (5) Total stellar mass formed in the simulation.\n");
+ fprintf(fp, "# Unit = %e gram\n", us->UnitMass_in_cgs);
+ fprintf(fp, "# Unit = %e solar mass\n",
+ 1.f / phys_const->const_solar_mass);
+ fprintf(fp,
+ "# (6) Number of stars formed in the current time step (no unit).\n");
+ fprintf(fp, "# (7) Mass of stars formed in the current time step.\n");
+ fprintf(fp, "# Unit = %e gram\n", us->UnitMass_in_cgs);
+ fprintf(fp, "# Unit = %e solar mass\n",
+ 1.f / phys_const->const_solar_mass);
+ fprintf(fp,
+ "# (0) (1) (2) (3) (4) "
+ " (5) (6) (7)\n");
+}
/**
* @brief Add the SFR tracer to the total active SFR of this cell
*
+ * Nothing to do here
+ *
* @param p the #part
* @param xp the #xpart
+ *
* @param sf the SFH logger struct
* @param dt_star The length of the time-step in physical internal units.
*/
@@ -108,6 +156,8 @@ INLINE static void star_formation_logger_log_active_part(
* @brief Add the SFR tracer to the total inactive SFR of this cell as long as
* the SFR tracer is larger than 0
*
+ * Nothing to do here
+ *
* @param p the #part
* @param xp the #xpart
* @param sf the SFH logger struct
@@ -116,4 +166,57 @@ INLINE static void star_formation_logger_log_inactive_part(
const struct part *p, const struct xpart *xp,
struct star_formation_history *sf) {}
+/**
+ * @brief add a star formation history struct to an other star formation history
+ * struct
+ *
+ * @param sf_add the star formation struct which we want to add to the star
+ * formation history
+ * @param sf_update the star formation structure which we want to update
+ */
+INLINE static void star_formation_logger_add_to_accumulator(
+ struct star_formation_history_accumulator *sf_update,
+ const struct star_formation_history *sf_add) {
+
+ /* Update the SFH structure */
+ sf_update->number_new_stars = sf_add->number_new_stars;
+ sf_update->new_stellar_mass = sf_add->new_stellar_mass;
+ sf_update->total_number_stars += sf_add->number_new_stars;
+ sf_update->total_stellar_mass += sf_add->new_stellar_mass;
+}
+
+/**
+ * @brief Write the final SFH to a file
+ *
+ * @param fp The file to write to.
+ * @param time the simulation time (time since Big Bang) in internal units.
+ * @param a the scale factor.
+ * @param z the redshift.
+ * @param sf the #star_formation_history struct.
+ * @param step The time-step of the simulation.
+ */
+INLINE static void star_formation_logger_write_to_log_file(
+ FILE *fp, const double time, const double a, const double z,
+ struct star_formation_history_accumulator sf, const int step) {
+
+ fprintf(fp, "%6d %16e %12.7f %14e %14ld %14e %14ld %14e\n", step, time, a, z,
+ sf.total_number_stars, sf.total_stellar_mass, sf.number_new_stars,
+ sf.new_stellar_mass);
+}
+
+/**
+ * @brief Initialize the star formation history struct in the #engine when
+ * starting the simulation.
+ *
+ * @param sfh the star_formation_history struct we want to initialize
+ */
+INLINE static void star_formation_logger_accumulator_init(
+ struct star_formation_history_accumulator *sfh) {
+ /* Initialize all values to zero */
+ sfh->new_stellar_mass = 0.f;
+ sfh->number_new_stars = 0;
+ sfh->total_number_stars = 0;
+ sfh->total_stellar_mass = 0.f;
+}
+
#endif /* SWIFT_GEAR_STARFORMATION_LOGGER_H */
diff --git a/src/star_formation/GEAR/star_formation_logger_struct.h b/src/star_formation/GEAR/star_formation_logger_struct.h
index 04b5dfdc038f7b684cfb1f2079d13eb312624b3f..63e3af06a7cd5af375662a1ba28dc5000a69dc3f 100644
--- a/src/star_formation/GEAR/star_formation_logger_struct.h
+++ b/src/star_formation/GEAR/star_formation_logger_struct.h
@@ -1,6 +1,7 @@
/*******************************************************************************
* This file is part of SWIFT.
- * Copyright (c) 2019 Folkert Nobels (nobels@strw.leidenuniv.nl)
+ * 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
@@ -19,7 +20,33 @@
#ifndef SWIFT_GEAR_STAR_FORMATION_LOGGER_STRUCT_H
#define SWIFT_GEAR_STAR_FORMATION_LOGGER_STRUCT_H
-/* Starformation history struct */
-struct star_formation_history {};
+/**
+ * Structure containing the star formation information from the cells.
+ */
+struct star_formation_history {
+ /*! Stellar mass created in the current timestep */
+ float new_stellar_mass;
-#endif /* SWIFT_NONE_STAR_FORMATION_STRUCT_H */
+ /*! Number of stars created in this timestep */
+ long int number_new_stars;
+};
+
+/**
+ * Structure containing the global star formation information (including time
+ * integrated variables).
+ */
+struct star_formation_history_accumulator {
+ /*! Total stellar mass from the begining of the simulation */
+ float total_stellar_mass;
+
+ /*! Total number of stars */
+ long int total_number_stars;
+
+ /*! Stellar mass created in the current timestep */
+ float new_stellar_mass;
+
+ /*! Number of stars created in this timestep */
+ long int number_new_stars;
+};
+
+#endif /* SWIFT_GEAR_STAR_FORMATION_LOGGER_STRUCT_H */
diff --git a/src/star_formation/GEAR/star_formation_struct.h b/src/star_formation/GEAR/star_formation_struct.h
index 9b4e216fd2955f29d89dade6ee46c0e1af715cdb..2e9a7548f83ca6ae9bb78ee7bcf4be69a0a31489 100644
--- a/src/star_formation/GEAR/star_formation_struct.h
+++ b/src/star_formation/GEAR/star_formation_struct.h
@@ -25,7 +25,27 @@
*/
struct star_formation_xpart_data {};
-/* Starformation struct */
-struct star_formation {};
+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)
+ */
+ float n_jeans_2_3;
+
+ /*! Maximal temperature for forming a star */
+ float maximal_temperature;
+
+ /*! Star formation efficiency */
+ float star_formation_efficiency;
+};
#endif /* SWIFT_GEAR_STAR_FORMATION_STRUCT_H */
diff --git a/src/star_formation/none/star_formation.h b/src/star_formation/none/star_formation.h
index 412bda2cc356729bb9f3e3657689fb5bb7b36d8f..6b83010ec5f8935778af8c9ad21010ff3452fc0e 100644
--- a/src/star_formation/none/star_formation.h
+++ b/src/star_formation/none/star_formation.h
@@ -220,6 +220,7 @@ star_formation_first_init_part(const struct phys_const* restrict phys_const,
* @param data The global star_formation information.
*/
__attribute__((always_inline)) INLINE static void star_formation_init_part(
- struct part* restrict p, const struct star_formation* data) {}
+ 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_logger.h b/src/star_formation/none/star_formation_logger.h
index b4e6987c03d295348fc8c22d66cb20d10e54378c..552df0c6cae533d1eb2678cb23dd675e0a058715 100644
--- a/src/star_formation/none/star_formation_logger.h
+++ b/src/star_formation/none/star_formation_logger.h
@@ -59,6 +59,18 @@ INLINE static void star_formation_logger_add(
struct star_formation_history *sf_update,
const struct star_formation_history *sf_add) {}
+/**
+ * @brief add a star formation history accumulator struct to an other star
+ * formation history struct
+ *
+ * @param sf_add the star formation accumulator struct which we want to add to
+ * the star formation history
+ * @param sf_update the star formation structure which we want to update
+ */
+INLINE static void star_formation_logger_add_to_accumulator(
+ struct star_formation_history_accumulator *sf_update,
+ const struct star_formation_history *sf_add) {}
+
/**
* @brief Initialize the star formation history structure in the #engine
*
@@ -67,6 +79,14 @@ INLINE static void star_formation_logger_add(
INLINE static void star_formation_logger_init(
struct star_formation_history *sfh) {}
+/**
+ * @brief Initialize the star formation history structure in the #engine
+ *
+ * @param sfh The pointer to the star formation history structure
+ */
+INLINE static void star_formation_logger_accumulator_init(
+ struct star_formation_history_accumulator *sfh) {}
+
/**
* @brief Write the final SFH to a file
*
@@ -74,12 +94,12 @@ INLINE static void star_formation_logger_init(
* @param time the simulation time (time since Big Bang) in internal units.
* @param a the scale factor.
* @param z the redshift.
- * @param sf the #star_formation_history struct.
+ * @param sf the #star_formation_history_accumulator struct.
* @param step The time-step of the simulation.
*/
INLINE static void star_formation_logger_write_to_log_file(
FILE *fp, const double time, const double a, const double z,
- const struct star_formation_history sf, const int step) {}
+ const struct star_formation_history_accumulator sf, const int step) {}
/**
* @brief Initialize the SFH logger file
diff --git a/src/star_formation/none/star_formation_logger_struct.h b/src/star_formation/none/star_formation_logger_struct.h
index 9efda271da96faf2088169fd75d0e3c01247a429..b60c64f2eb47894db828a1bde0ca20803892c7db 100644
--- a/src/star_formation/none/star_formation_logger_struct.h
+++ b/src/star_formation/none/star_formation_logger_struct.h
@@ -22,4 +22,10 @@
/* Starformation history struct */
struct star_formation_history {};
+/* Starformation history accumulator struct.
+ This structure is only defined in the engine and
+ allows the user to integrate some quantities over time.
+*/
+struct star_formation_history_accumulator {};
+
#endif /* SWIFT_NONE_STAR_FORMATION_STRUCT_H */
diff --git a/src/star_formation/none/star_formation_struct.h b/src/star_formation/none/star_formation_struct.h
index 27a2adaf83d0a02a0d08e7eef8b45bea630689e4..2f5241a58caf1ca70fa98a40d467c8ff5a3237f7 100644
--- a/src/star_formation/none/star_formation_struct.h
+++ b/src/star_formation/none/star_formation_struct.h
@@ -25,4 +25,10 @@
*/
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/stars.h b/src/stars.h
index dd8390e0206580fc2a07a08e51bb69c6ee5ab5ed..dea6e07a87cabd7d1778ec2a850be4f3b16b04b0 100644
--- a/src/stars.h
+++ b/src/stars.h
@@ -29,6 +29,9 @@
#elif defined(STARS_EAGLE)
#include "./stars/EAGLE/stars.h"
#include "./stars/EAGLE/stars_iact.h"
+#elif defined(STARS_GEAR)
+#include "./stars/GEAR/stars.h"
+#include "./stars/GEAR/stars_iact.h"
#else
#error "Invalid choice of star model"
#endif
diff --git a/src/stars/GEAR/stars.h b/src/stars/GEAR/stars.h
new file mode 100644
index 0000000000000000000000000000000000000000..467aaa164ba9762d85f8dfae85db86ff76ae779e
--- /dev/null
+++ b/src/stars/GEAR/stars.h
@@ -0,0 +1,172 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (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_GEAR_STARS_H
+#define SWIFT_GEAR_STARS_H
+
+#include <float.h>
+#include "minmax.h"
+
+/**
+ * @brief Computes the gravity time-step of a given star particle.
+ *
+ * @param sp Pointer to the s-particle data.
+ */
+__attribute__((always_inline)) INLINE static float stars_compute_timestep(
+ const struct spart* const sp) {
+
+ return FLT_MAX;
+}
+
+/**
+ * @brief Initialises the s-particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param sp The particle to act upon
+ * @param stars_properties The properties of the stellar model.
+ */
+__attribute__((always_inline)) INLINE static void stars_first_init_spart(
+ struct spart* sp, const struct stars_props* stars_properties) {
+
+ sp->time_bin = 0;
+}
+
+/**
+ * @brief Prepares a s-particle for its interactions
+ *
+ * @param sp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void stars_init_spart(
+ struct spart* sp) {
+
+#ifdef DEBUG_INTERACTIONS_STARS
+ for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
+ sp->ids_ngbs_density[i] = -1;
+ sp->num_ngb_density = 0;
+#endif
+
+ sp->density.wcount = 0.f;
+ sp->density.wcount_dh = 0.f;
+}
+
+/**
+ * @brief Predict additional particle fields forward in time when drifting
+ *
+ * @param sp The particle
+ * @param dt_drift The drift time-step for positions.
+ */
+__attribute__((always_inline)) INLINE static void stars_predict_extra(
+ struct spart* restrict sp, float dt_drift) {}
+
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param sp The particle.
+ */
+__attribute__((always_inline)) INLINE static void stars_reset_predicted_values(
+ struct spart* restrict sp) {}
+
+/**
+ * @brief Finishes the calculation of (non-gravity) forces acting on stars
+ *
+ * @param sp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void stars_end_feedback(
+ struct spart* sp) {}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * @param sp The particle to act upon
+ * @param dt The time-step for this kick
+ */
+__attribute__((always_inline)) INLINE static void stars_kick_extra(
+ struct spart* sp, float dt) {}
+
+/**
+ * @brief Finishes the calculation of density on stars
+ *
+ * @param sp The particle to act upon
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void stars_end_density(
+ struct spart* sp, const struct cosmology* cosmo) {
+
+ /* Some smoothing length multiples. */
+ const float h = sp->h;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+ const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
+
+ /* Finish the calculation by inserting the missing h-factors */
+ sp->density.wcount *= h_inv_dim;
+ sp->density.wcount_dh *= h_inv_dim_plus_one;
+}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #spart has 0
+ * ngbs.
+ *
+ * @param sp The particle to act upon
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void stars_spart_has_no_neighbours(
+ struct spart* restrict sp, const struct cosmology* cosmo) {
+
+ /* Some smoothing length multiples. */
+ const float h = sp->h;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+
+ /* Re-set problematic values */
+ sp->density.wcount = kernel_root * h_inv_dim;
+ sp->density.wcount_dh = 0.f;
+}
+
+/**
+ * @brief Reset acceleration fields of a particle
+ *
+ * This is the equivalent of hydro_reset_acceleration.
+ * We do not compute the acceleration on star, therefore no need to use it.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void stars_reset_feedback(
+ struct spart* restrict p) {
+
+#ifdef DEBUG_INTERACTIONS_STARS
+ for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
+ p->ids_ngbs_force[i] = -1;
+ p->num_ngb_force = 0;
+#endif
+}
+
+/**
+ * @brief Initializes constants related to stellar evolution, initializes imf,
+ * reads and processes yield tables
+ *
+ * @param params swift_params parameters structure
+ * @param stars stars_props data structure
+ */
+inline static void stars_evolve_init(struct swift_params* params,
+ struct stars_props* restrict stars) {}
+
+#endif /* SWIFT_GEAR_STARS_H */
diff --git a/src/stars/GEAR/stars_debug.h b/src/stars/GEAR/stars_debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..41953367f6c8771ffd7460cee493c8330ecd874b
--- /dev/null
+++ b/src/stars/GEAR/stars_debug.h
@@ -0,0 +1,31 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (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_GEAR_STARS_DEBUG_H
+#define SWIFT_GEAR_STARS_DEBUG_H
+
+__attribute__((always_inline)) INLINE static void stars_debug_particle(
+ const struct spart* p) {
+ printf(
+ "x=[%.3e,%.3e,%.3e], "
+ "v_full=[%.3e,%.3e,%.3e] p->mass=%.3e \n t_begin=%d, t_end=%d\n",
+ p->x[0], p->x[1], p->x[2], p->v_full[0], p->v_full[1], p->v_full[2],
+ p->mass, p->ti_begin, p->ti_end);
+}
+
+#endif /* SWIFT_GEAR_STARS_DEBUG_H */
diff --git a/src/stars/GEAR/stars_iact.h b/src/stars/GEAR/stars_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..c7bda43fc0c66fcc890e6b05bacf871edfd10b5a
--- /dev/null
+++ b/src/stars/GEAR/stars_iact.h
@@ -0,0 +1,94 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (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_GEAR_STARS_IACT_H
+#define SWIFT_GEAR_STARS_IACT_H
+
+/**
+ * @brief Density interaction between two particles (non-symmetric).
+ *
+ * @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 si First sparticle.
+ * @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_stars_density(const float r2, const float *dx,
+ const float hi, const float hj,
+ struct spart *restrict si,
+ const struct part *restrict pj, const float a,
+ const float H) {
+
+ float wi, wi_dx;
+
+ /* Get r and 1/r. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
+
+ /* Compute the kernel function */
+ const float hi_inv = 1.0f / hi;
+ const float ui = r * hi_inv;
+ kernel_deval(ui, &wi, &wi_dx);
+
+ /* Compute contribution to the number of neighbours */
+ si->density.wcount += wi;
+ si->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+#ifdef DEBUG_INTERACTIONS_STARS
+ /* Update ngb counters */
+ if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS_STARS)
+ si->ids_ngbs_density[si->num_ngb_density] = pj->id;
+
+ /* Update ngb counters */
+ ++si->num_ngb_density;
+#endif
+}
+
+/**
+ * @brief Feedback interaction between two particles (non-symmetric).
+ *
+ * @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 si First sparticle.
+ * @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_stars_feedback(const float r2, const float *dx,
+ const float hi, const float hj,
+ struct spart *restrict si,
+ struct part *restrict pj, const float a,
+ const float H) {
+#ifdef DEBUG_INTERACTIONS_STARS
+ /* Update ngb counters */
+ if (si->num_ngb_force < MAX_NUM_OF_NEIGHBOURS_STARS)
+ si->ids_ngbs_force[si->num_ngb_force] = pj->id;
+
+ /* Update ngb counters */
+ ++si->num_ngb_force;
+#endif
+}
+
+#endif
diff --git a/src/stars/GEAR/stars_io.h b/src/stars/GEAR/stars_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..ebd72aa50a4194bf8c6f747e55d265ace0550c35
--- /dev/null
+++ b/src/stars/GEAR/stars_io.h
@@ -0,0 +1,248 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (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_GEAR_STARS_IO_H
+#define SWIFT_GEAR_STARS_IO_H
+
+#include "io_properties.h"
+#include "stars_part.h"
+
+/**
+ * @brief Specifies which s-particle fields to read from a dataset
+ *
+ * @param sparts The s-particle array.
+ * @param list The list of i/o properties to read.
+ * @param num_fields The number of i/o fields to read.
+ */
+INLINE static void stars_read_particles(struct spart *sparts,
+ struct io_props *list,
+ int *num_fields) {
+
+ /* Say how much we want to read */
+ *num_fields = 5;
+
+ /* List what we want to read */
+ list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
+ UNIT_CONV_LENGTH, sparts, x);
+ list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
+ UNIT_CONV_SPEED, sparts, v);
+ list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
+ sparts, mass);
+ list[3] = io_make_input_field("ParticleIDs", LONGLONG, 1, COMPULSORY,
+ UNIT_CONV_NO_UNITS, sparts, id);
+ list[4] = io_make_input_field("SmoothingLength", FLOAT, 1, OPTIONAL,
+ UNIT_CONV_LENGTH, sparts, h);
+}
+
+/**
+ * @brief Specifies which s-particle fields to write to a dataset
+ *
+ * @param sparts The s-particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
+ * @param with_cosmology Is it a cosmological run?
+ */
+INLINE static void stars_write_particles(const struct spart *sparts,
+ struct io_props *list, int *num_fields,
+ const int with_cosmology) {
+
+ /* Say how much we want to write */
+ *num_fields = 9;
+
+ /* List what we want to write */
+ list[0] =
+ io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, 1.,
+ sparts, x, "Co-moving positions of the particles");
+
+ list[1] = io_make_output_field(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, 0.f, sparts, 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, 0.f,
+ sparts, mass, "Masses of the particles");
+
+ list[3] =
+ io_make_output_field("ParticleIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
+ sparts, id, "Unique IDs of the particles");
+
+ list[4] = io_make_output_field(
+ "SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, sparts, h,
+ "Co-moving smoothing lengths (FWHM of the kernel) of the particles");
+
+ if (with_cosmology) {
+ list[5] = io_make_output_field(
+ "BirthScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, sparts,
+ birth_scale_factor, "Scale-factors at which the stars were born");
+ } else {
+ list[5] = io_make_output_field("BirthTimes", FLOAT, 1, UNIT_CONV_TIME, 0.f,
+ sparts, birth_time,
+ "Times at which the stars were born");
+ }
+
+ list[6] = io_make_output_field(
+ "BirthDensities", FLOAT, 1, UNIT_CONV_DENSITY, 0.f, sparts, birth.density,
+ "Physical densities at the time of birth of the gas particles that "
+ "turned into stars (note that "
+ "we store the physical density at the birth redshift, no conversion is "
+ "needed)");
+
+ list[7] =
+ io_make_output_field("BirthTemperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE,
+ 0.f, sparts, birth.temperature,
+ "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,
+ sparts, birth.mass,
+ "Masses of the star particles at birth time");
+
+#ifdef DEBUG_INTERACTIONS_STARS
+
+ list += *num_fields;
+ *num_fields += 4;
+
+ list[0] = io_make_output_field("Num_ngb_density", INT, 1, UNIT_CONV_NO_UNITS,
+ sparts, num_ngb_density);
+ list[1] = io_make_output_field("Num_ngb_force", INT, 1, UNIT_CONV_NO_UNITS,
+ sparts, num_ngb_force);
+ list[2] = io_make_output_field("Ids_ngb_density", LONGLONG,
+ MAX_NUM_OF_NEIGHBOURS_STARS,
+ UNIT_CONV_NO_UNITS, sparts, ids_ngbs_density);
+ list[3] = io_make_output_field("Ids_ngb_force", LONGLONG,
+ MAX_NUM_OF_NEIGHBOURS_STARS,
+ UNIT_CONV_NO_UNITS, sparts, ids_ngbs_force);
+#endif
+}
+
+/**
+ * @brief Initialize the global properties of the stars scheme.
+ *
+ * By default, takes the values provided by the hydro.
+ *
+ * @param sp The #stars_props.
+ * @param phys_const The physical constants in the internal unit system.
+ * @param us The internal unit system.
+ * @param params The parsed parameters.
+ * @param p The already read-in properties of the hydro scheme.
+ * @param cosmo The cosmological model.
+ */
+INLINE static void stars_props_init(struct stars_props *sp,
+ const struct phys_const *phys_const,
+ const struct unit_system *us,
+ struct swift_params *params,
+ const struct hydro_props *p,
+ const struct cosmology *cosmo) {
+
+ /* Kernel properties */
+ sp->eta_neighbours = parser_get_opt_param_float(
+ params, "Stars:resolution_eta", p->eta_neighbours);
+
+ /* Tolerance for the smoothing length Newton-Raphson scheme */
+ sp->h_tolerance =
+ parser_get_opt_param_float(params, "Stars:h_tolerance", p->h_tolerance);
+
+ /* Get derived properties */
+ sp->target_neighbours = pow_dimension(sp->eta_neighbours) * kernel_norm;
+ const float delta_eta = sp->eta_neighbours * (1.f + sp->h_tolerance);
+ sp->delta_neighbours =
+ (pow_dimension(delta_eta) - pow_dimension(sp->eta_neighbours)) *
+ kernel_norm;
+
+ /* Number of iterations to converge h */
+ sp->max_smoothing_iterations = parser_get_opt_param_int(
+ params, "Stars:max_ghost_iterations", p->max_smoothing_iterations);
+
+ /* Time integration properties */
+ const float max_volume_change =
+ parser_get_opt_param_float(params, "Stars:max_volume_change", -1);
+ if (max_volume_change == -1)
+ sp->log_max_h_change = p->log_max_h_change;
+ else
+ sp->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv));
+}
+
+/**
+ * @brief Print the global properties of the stars scheme.
+ *
+ * @param sp The #stars_props.
+ */
+INLINE static void stars_props_print(const struct stars_props *sp) {
+
+ /* Now stars */
+ message("Stars kernel: %s with eta=%f (%.2f neighbours).", kernel_name,
+ sp->eta_neighbours, sp->target_neighbours);
+
+ message("Stars relative tolerance in h: %.5f (+/- %.4f neighbours).",
+ sp->h_tolerance, sp->delta_neighbours);
+
+ message(
+ "Stars integration: Max change of volume: %.2f "
+ "(max|dlog(h)/dt|=%f).",
+ pow_dimension(expf(sp->log_max_h_change)), sp->log_max_h_change);
+
+ message("Maximal iterations in ghost task set to %d",
+ sp->max_smoothing_iterations);
+}
+
+#if defined(HAVE_HDF5)
+INLINE static void stars_props_print_snapshot(hid_t h_grpstars,
+ const struct stars_props *sp) {
+
+ io_write_attribute_s(h_grpstars, "Kernel function", kernel_name);
+ io_write_attribute_f(h_grpstars, "Kernel target N_ngb",
+ sp->target_neighbours);
+ io_write_attribute_f(h_grpstars, "Kernel delta N_ngb", sp->delta_neighbours);
+ io_write_attribute_f(h_grpstars, "Kernel eta", sp->eta_neighbours);
+ io_write_attribute_f(h_grpstars, "Smoothing length tolerance",
+ sp->h_tolerance);
+ io_write_attribute_f(h_grpstars, "Volume log(max(delta h))",
+ sp->log_max_h_change);
+ io_write_attribute_f(h_grpstars, "Volume max change time-step",
+ pow_dimension(expf(sp->log_max_h_change)));
+ io_write_attribute_i(h_grpstars, "Max ghost iterations",
+ sp->max_smoothing_iterations);
+}
+#endif
+
+/**
+ * @brief Write a #stars_props struct to the given FILE as a stream of bytes.
+ *
+ * @param p the struct
+ * @param stream the file stream
+ */
+INLINE static void stars_props_struct_dump(const struct stars_props *p,
+ FILE *stream) {
+ restart_write_blocks((void *)p, sizeof(struct stars_props), 1, stream,
+ "starsprops", "stars props");
+}
+
+/**
+ * @brief Restore a stars_props struct from the given FILE as a stream of
+ * bytes.
+ *
+ * @param p the struct
+ * @param stream the file stream
+ */
+INLINE static void stars_props_struct_restore(const struct stars_props *p,
+ FILE *stream) {
+ restart_read_blocks((void *)p, sizeof(struct stars_props), 1, stream, NULL,
+ "stars props");
+}
+
+#endif /* SWIFT_GEAR_STAR_IO_H */
diff --git a/src/stars/GEAR/stars_part.h b/src/stars/GEAR/stars_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..bf68a580ef009f5a814fa17123301b5585d6084c
--- /dev/null
+++ b/src/stars/GEAR/stars_part.h
@@ -0,0 +1,154 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (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_GEAR_STAR_PART_H
+#define SWIFT_GEAR_STAR_PART_H
+
+/* Some standard headers. */
+#include <stdlib.h>
+
+/* Read additional subgrid models */
+#include "chemistry_struct.h"
+#include "feedback_struct.h"
+#include "tracers_struct.h"
+
+/**
+ * @brief Particle fields for the star particles.
+ *
+ * All quantities related to gravity are stored in the associate #gpart.
+ */
+struct spart {
+
+ /*! Particle ID. */
+ long long id;
+
+ /*! Pointer to corresponding gravity part. */
+ struct gpart* gpart;
+
+ /*! Particle position. */
+ double x[3];
+
+ /* Offset between current position and position at last tree rebuild. */
+ float x_diff[3];
+
+ /* Offset between current position and position at last tree rebuild. */
+ float x_diff_sort[3];
+
+ /*! Particle velocity. */
+ float v[3];
+
+ /*! Star mass */
+ float mass;
+
+ /* Particle cutoff radius. */
+ float h;
+
+ /*! Union for the birth time and birth scale factor */
+ union {
+
+ /*! Birth time */
+ float birth_time;
+
+ /*! Birth scale factor */
+ float birth_scale_factor;
+ };
+
+ /*! Particle time bin */
+ timebin_t time_bin;
+
+ struct {
+
+ /* Number of neighbours. */
+ float wcount;
+
+ /* Number of neighbours spatial derivative. */
+ float wcount_dh;
+
+ } density;
+
+ struct {
+ /*! birth density*/
+ float density;
+
+ /*! birth temperature*/
+ float temperature;
+
+ /*! birth mass */
+ float mass;
+ } birth;
+
+ /*! Feedback structure */
+ struct feedback_spart_data feedback_data;
+
+ /*! Tracer structure */
+ struct tracers_xpart_data tracers_data;
+
+ /*! Chemistry structure */
+ struct chemistry_part_data chemistry_data;
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+ /* Time of the last drift */
+ integertime_t ti_drift;
+
+ /* Time of the last kick */
+ integertime_t ti_kick;
+
+#endif
+
+#ifdef DEBUG_INTERACTIONS_STARS
+ /*! Number of interactions in the density SELF and PAIR */
+ int num_ngb_density;
+
+ /*! List of interacting particles in the density SELF and PAIR */
+ long long ids_ngbs_density[MAX_NUM_OF_NEIGHBOURS_STARS];
+
+ /*! Number of interactions in the force SELF and PAIR */
+ int num_ngb_force;
+
+ /*! List of interacting particles in the force SELF and PAIR */
+ long long ids_ngbs_force[MAX_NUM_OF_NEIGHBOURS_STARS];
+#endif
+
+} SWIFT_STRUCT_ALIGN;
+
+/**
+ * @brief Contains all the constants and parameters of the stars scheme
+ */
+struct stars_props {
+
+ /*! Resolution parameter */
+ float eta_neighbours;
+
+ /*! Target weightd number of neighbours (for info only)*/
+ float target_neighbours;
+
+ /*! Smoothing length tolerance */
+ float h_tolerance;
+
+ /*! Tolerance on neighbour number (for info only)*/
+ float delta_neighbours;
+
+ /*! Maximal number of iterations to converge h */
+ int max_smoothing_iterations;
+
+ /*! Maximal change of h over one time-step */
+ float log_max_h_change;
+};
+
+#endif /* SWIFT_GEAR_STAR_PART_H */
diff --git a/src/swift.h b/src/swift.h
index b9f8818d8b833231971abb1afb36ee4507648488..7790ec6f9203a6534e2f569017d90090352d8f28 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -64,6 +64,7 @@
#include "periodic.h"
#include "physical_constants.h"
#include "potential.h"
+#include "pressure_floor.h"
#include "profiler.h"
#include "queue.h"
#include "random.h"