diff --git a/doc/RTD/source/Planetary/equations_of_state.rst b/doc/RTD/source/Planetary/equations_of_state.rst
index d20cbd83ff0efdbc1f73c3a59a310d71155abda1..37b175c53737288374fe1a4afec07a280c3bbd53 100644
--- a/doc/RTD/source/Planetary/equations_of_state.rst
+++ b/doc/RTD/source/Planetary/equations_of_state.rst
@@ -36,6 +36,8 @@ planetary initial conditions, `WoMa <https://github.com/srbonilla/WoMa>`_:
+ Granite: ``101``
+ Water: ``102``
+ Basalt: ``103``
+ + Ice: ``104``
+ + Custom user-provided parameters: ``190``, ``191``, ..., ``199``
+ Hubbard \& MacFarlane (1980): ``2``
+ Hydrogen-helium atmosphere: ``200``
+ Ice H20-CH4-NH3 mix: ``201``
@@ -45,6 +47,10 @@ planetary initial conditions, `WoMa <https://github.com/srbonilla/WoMa>`_:
+ Basalt (7530): ``301``
+ Water (7154): ``302``
+ Senft \& Stewart (2008) water: ``303``
+ + AQUA, Haldemann, J. et al. (2020) water: ``304``
+ + Chabrier, G. et al. (2019) Hydrogen: ``305``
+ + Chabrier, G. et al. (2019) Helium: ``306``
+ + Chabrier & Debras (2021) H/He mixture Y=0.245 (Jupiter): ``307``
+ ANEOS (in SESAME-style tables): ``4``
+ Forsterite (Stewart et al. 2019): ``400``
+ Iron (Stewart, zenodo.org/record/3866507): ``401``
@@ -56,7 +62,7 @@ The data files for the tabulated EoS can be downloaded using
the ``examples/Planetary/EoSTables/get_eos_tables.sh`` script.
To enable one or multiple EoS, the corresponding ``planetary_use_*:``
-flag(s) must be set to ``1`` in the parameter file for a simulation,
+flag(s) must be set from ``0`` to ``1`` in the parameter file for a simulation,
along with the path to any table files, which are set by the
``planetary_*_table_file:`` parameters,
as detailed in :ref:`Parameters_eos` and ``examples/parameter_example.yml``.
diff --git a/examples/Planetary/EoSTables/get_eos_tables.sh b/examples/Planetary/EoSTables/get_eos_tables.sh
index aaae425fd5dd9e22490291f966d1f9c76e76f30e..e47f08816dac41099cde2b980772c7320c9567bb 100755
--- a/examples/Planetary/EoSTables/get_eos_tables.sh
+++ b/examples/Planetary/EoSTables/get_eos_tables.sh
@@ -7,6 +7,10 @@ wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/SESAME_basalt_7530.txt
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/SESAME_iron_2140.txt
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/SESAME_water_7154.txt
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/SS08_water.txt
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/AQUA_H20.txt
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/CMS19_H.txt
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/CMS19_He.txt
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/CD21_HHe.txt
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/ANEOS_forsterite_S19.txt
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/ANEOS_iron_S20.txt
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index ca3258a13db07cd11957cddcdb85d8004e56f705..e980b021603b6632c42f95b30336c2b7b21cc75a 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -314,6 +314,17 @@ EoS:
planetary_use_Til_granite: 1 # Tillotson granite, material ID 101
planetary_use_Til_water: 0 # Tillotson water, material ID 102
planetary_use_Til_basalt: 0 # Tillotson basalt, material ID 103
+ planetary_use_Til_ice: 0 # Tillotson ice, material ID 104
+ planetary_use_Til_custom_0: 0 # Tillotson generic user-provided parameters, material IDs 19[0-9]
+ planetary_use_Til_custom_1: 0
+ planetary_use_Til_custom_2: 0
+ planetary_use_Til_custom_3: 0
+ planetary_use_Til_custom_4: 0
+ planetary_use_Til_custom_5: 0
+ planetary_use_Til_custom_6: 0
+ planetary_use_Til_custom_7: 0
+ planetary_use_Til_custom_8: 0
+ planetary_use_Til_custom_9: 0
planetary_use_HM80_HHe: 0 # Hubbard & MacFarlane (1980) hydrogen-helium atmosphere, material ID 200
planetary_use_HM80_ice: 0 # Hubbard & MacFarlane (1980) H20-CH4-NH3 ice mix, material ID 201
planetary_use_HM80_rock: 0 # Hubbard & MacFarlane (1980) SiO2-MgO-FeS-FeO rock mix, material ID 202
@@ -321,6 +332,10 @@ EoS:
planetary_use_SESAME_basalt: 0 # SESAME basalt 7530, material ID 301
planetary_use_SESAME_water: 0 # SESAME water 7154, material ID 302
planetary_use_SS08_water: 0 # Senft & Stewart (2008) SESAME-like water, material ID 303
+ planetary_use_AQUA: 0 # Haldemann, J. et al. (2020) water, material ID 304
+ planetary_use_CMS19_H: 0 # Chabrier, G. et al. (2019) Hydrogen, material ID 305
+ planetary_use_CMS19_He: 0 # Chabrier, G. et al. (2019) Hydrogen, material ID 306
+ planetary_use_CD21_HHe: 0 # Chabrier & Debras (2021) H/He mixture Y=0.245 (Jupiter), material ID 307
planetary_use_ANEOS_forsterite: 0 # ANEOS forsterite (Stewart et al. 2019), material ID 400
planetary_use_ANEOS_iron: 0 # ANEOS iron (Stewart 2020), material ID 401
planetary_use_ANEOS_Fe85Si15: 0 # ANEOS Fe85Si15 (Stewart 2020), material ID 402
@@ -335,6 +350,16 @@ EoS:
planetary_use_custom_8: 0
planetary_use_custom_9: 0
# Tablulated EoS file paths.
+ planetary_Til_custom_0_param_file: ./EoSTables/Til_custom_0.txt
+ planetary_Til_custom_1_param_file: ./EoSTables/Til_custom_1.txt
+ planetary_Til_custom_2_param_file: ./EoSTables/Til_custom_2.txt
+ planetary_Til_custom_3_param_file: ./EoSTables/Til_custom_3.txt
+ planetary_Til_custom_4_param_file: ./EoSTables/Til_custom_4.txt
+ planetary_Til_custom_5_param_file: ./EoSTables/Til_custom_5.txt
+ planetary_Til_custom_6_param_file: ./EoSTables/Til_custom_6.txt
+ planetary_Til_custom_7_param_file: ./EoSTables/Til_custom_7.txt
+ planetary_Til_custom_8_param_file: ./EoSTables/Til_custom_8.txt
+ planetary_Til_custom_9_param_file: ./EoSTables/Til_custom_9.txt
planetary_HM80_HHe_table_file: ./EoSTables/HM80_HHe.txt
planetary_HM80_ice_table_file: ./EoSTables/HM80_ice.txt
planetary_HM80_rock_table_file: ./EoSTables/HM80_rock.txt
diff --git a/src/equation_of_state/barotropic/equation_of_state.h b/src/equation_of_state/barotropic/equation_of_state.h
index d91b0c0ce4cfe8b864f773c877d755f8dfcf77a7..71be2944437879bb0fa65142d7674fa830ed824c 100644
--- a/src/equation_of_state/barotropic/equation_of_state.h
+++ b/src/equation_of_state/barotropic/equation_of_state.h
@@ -56,7 +56,7 @@ struct eos_parameters {
* @param entropy The entropy \f$A\f$.
*/
__attribute__((always_inline, const)) INLINE static float
-gas_internal_energy_from_entropy(float density, float entropy) {
+gas_internal_energy_from_entropy(const float density, const float entropy) {
const float density_frac = density * eos.inverse_core_density;
const float density_factor = pow_gamma(density_frac);
@@ -76,7 +76,7 @@ gas_internal_energy_from_entropy(float density, float entropy) {
* @param entropy The entropy \f$A\f$.
*/
__attribute__((always_inline, const)) INLINE static float
-gas_pressure_from_entropy(float density, float entropy) {
+gas_pressure_from_entropy(const float density, const float entropy) {
const float density_frac = density * eos.inverse_core_density;
const float density_factor = pow_gamma(density_frac);
@@ -96,7 +96,7 @@ gas_pressure_from_entropy(float density, float entropy) {
* @return The entropy \f$A\f$.
*/
__attribute__((always_inline, const)) INLINE static float
-gas_entropy_from_pressure(float density, float pressure) {
+gas_entropy_from_pressure(const float density, const float pressure) {
const float density_frac = density * eos.inverse_core_density;
const float density_factor = pow_gamma(density_frac);
@@ -116,7 +116,7 @@ gas_entropy_from_pressure(float density, float pressure) {
* @param entropy The entropy \f$A\f$.
*/
__attribute__((always_inline, const)) INLINE static float
-gas_soundspeed_from_entropy(float density, float entropy) {
+gas_soundspeed_from_entropy(const float density, const float entropy) {
const float density_frac = density * eos.inverse_core_density;
const float density_factor = pow_gamma(density_frac);
@@ -135,7 +135,7 @@ gas_soundspeed_from_entropy(float density, float entropy) {
* @param u The internal energy \f$u\f$
*/
__attribute__((always_inline, const)) INLINE static float
-gas_entropy_from_internal_energy(float density, float u) {
+gas_entropy_from_internal_energy(const float density, const float u) {
const float density_frac = density * eos.inverse_core_density;
const float density_factor = pow_gamma(density_frac);
@@ -155,7 +155,7 @@ gas_entropy_from_internal_energy(float density, float u) {
* @param u The internal energy \f$u\f$
*/
__attribute__((always_inline, const)) INLINE static float
-gas_pressure_from_internal_energy(float density, float u) {
+gas_pressure_from_internal_energy(const float density, const float u) {
const float density_frac = density * eos.inverse_core_density;
const float density_factor = pow_gamma(density_frac);
@@ -175,7 +175,7 @@ gas_pressure_from_internal_energy(float density, float u) {
* @return The internal energy \f$u\f$.
*/
__attribute__((always_inline, const)) INLINE static float
-gas_internal_energy_from_pressure(float density, float pressure) {
+gas_internal_energy_from_pressure(const float density, const float pressure) {
const float density_frac = density * eos.inverse_core_density;
const float density_factor = pow_gamma(density_frac);
@@ -195,7 +195,7 @@ gas_internal_energy_from_pressure(float density, float pressure) {
* @param u The internal energy \f$u\f$
*/
__attribute__((always_inline, const)) INLINE static float
-gas_soundspeed_from_internal_energy(float density, float u) {
+gas_soundspeed_from_internal_energy(const float density, const float u) {
const float density_frac = density * eos.inverse_core_density;
const float density_factor = pow_gamma(density_frac);
@@ -214,7 +214,7 @@ gas_soundspeed_from_internal_energy(float density, float u) {
* @param P The pressure \f$P\f$
*/
__attribute__((always_inline, const)) INLINE static float
-gas_soundspeed_from_pressure(float density, float P) {
+gas_soundspeed_from_pressure(const float density, const float P) {
const float density_frac = density * eos.inverse_core_density;
const float density_factor = pow_gamma(density_frac);
@@ -241,7 +241,7 @@ INLINE static void eos_init(struct eos_parameters *e,
parser_get_param_float(params, "EoS:barotropic_vacuum_sound_speed");
e->vacuum_sound_speed2 = vacuum_sound_speed * vacuum_sound_speed;
e->inverse_core_density =
- 1. / parser_get_param_float(params, "EoS:barotropic_core_density");
+ 1.f / parser_get_param_float(params, "EoS:barotropic_core_density");
}
/**
* @brief Print the equation of state
diff --git a/src/equation_of_state/ideal_gas/equation_of_state.h b/src/equation_of_state/ideal_gas/equation_of_state.h
index cfdff5c7fd4af8052cf4a36b8bace66718e235e5..61d76e4fb36d204f9c40ef3c76a7d875ff59a5be 100644
--- a/src/equation_of_state/ideal_gas/equation_of_state.h
+++ b/src/equation_of_state/ideal_gas/equation_of_state.h
@@ -46,7 +46,7 @@ struct eos_parameters {};
* @param entropy The entropy \f$A\f$.
*/
__attribute__((always_inline, const)) INLINE static float
-gas_internal_energy_from_entropy(float density, float entropy) {
+gas_internal_energy_from_entropy(const float density, const float entropy) {
return entropy * pow_gamma_minus_one(density) *
hydro_one_over_gamma_minus_one;
@@ -61,7 +61,7 @@ gas_internal_energy_from_entropy(float density, float entropy) {
* @param entropy The entropy \f$A\f$.
*/
__attribute__((always_inline, const)) INLINE static float
-gas_pressure_from_entropy(float density, float entropy) {
+gas_pressure_from_entropy(const float density, const float entropy) {
return entropy * pow_gamma(density);
}
@@ -76,7 +76,7 @@ gas_pressure_from_entropy(float density, float entropy) {
* @return The entropy \f$A\f$.
*/
__attribute__((always_inline, const)) INLINE static float
-gas_entropy_from_pressure(float density, float pressure) {
+gas_entropy_from_pressure(const float density, const float pressure) {
return pressure * pow_minus_gamma(density);
}
@@ -90,7 +90,7 @@ gas_entropy_from_pressure(float density, float pressure) {
* @param entropy The entropy \f$A\f$.
*/
__attribute__((always_inline, const)) INLINE static float
-gas_soundspeed_from_entropy(float density, float entropy) {
+gas_soundspeed_from_entropy(const float density, const float entropy) {
return sqrtf(hydro_gamma * pow_gamma_minus_one(density) * entropy);
}
@@ -104,7 +104,7 @@ gas_soundspeed_from_entropy(float density, float entropy) {
* @param u The internal energy \f$u\f$
*/
__attribute__((always_inline, const)) INLINE static float
-gas_entropy_from_internal_energy(float density, float u) {
+gas_entropy_from_internal_energy(const float density, const float u) {
return hydro_gamma_minus_one * u * pow_minus_gamma_minus_one(density);
}
@@ -118,7 +118,7 @@ gas_entropy_from_internal_energy(float density, float u) {
* @param u The internal energy \f$u\f$
*/
__attribute__((always_inline, const)) INLINE static float
-gas_pressure_from_internal_energy(float density, float u) {
+gas_pressure_from_internal_energy(const float density, const float u) {
return hydro_gamma_minus_one * u * density;
}
@@ -133,7 +133,7 @@ gas_pressure_from_internal_energy(float density, float u) {
* @return The internal energy \f$u\f$.
*/
__attribute__((always_inline, const)) INLINE static float
-gas_internal_energy_from_pressure(float density, float pressure) {
+gas_internal_energy_from_pressure(const float density, const float pressure) {
return hydro_one_over_gamma_minus_one * pressure / density;
}
@@ -147,7 +147,7 @@ gas_internal_energy_from_pressure(float density, float pressure) {
* @param u The internal energy \f$u\f$
*/
__attribute__((always_inline, const)) INLINE static float
-gas_soundspeed_from_internal_energy(float density, float u) {
+gas_soundspeed_from_internal_energy(const float density, const float u) {
return sqrtf(u * hydro_gamma * hydro_gamma_minus_one);
}
@@ -161,7 +161,7 @@ gas_soundspeed_from_internal_energy(float density, float u) {
* @param P The pressure \f$P\f$
*/
__attribute__((always_inline, const)) INLINE static float
-gas_soundspeed_from_pressure(float density, float P) {
+gas_soundspeed_from_pressure(const float density, const float P) {
return sqrtf(hydro_gamma * P / density);
}
diff --git a/src/equation_of_state/isothermal/equation_of_state.h b/src/equation_of_state/isothermal/equation_of_state.h
index 3ee1e18ebb0b8dcc6be8e02d4ee640edadfc35d6..0ff1b1a5063146e9630c588683d57c0993edff60 100644
--- a/src/equation_of_state/isothermal/equation_of_state.h
+++ b/src/equation_of_state/isothermal/equation_of_state.h
@@ -50,7 +50,7 @@ struct eos_parameters {
* @param entropy The entropy \f$S\f$.
*/
__attribute__((always_inline)) INLINE static float
-gas_internal_energy_from_entropy(float density, float entropy) {
+gas_internal_energy_from_entropy(const float density, const float entropy) {
return eos.isothermal_internal_energy;
}
@@ -65,7 +65,7 @@ gas_internal_energy_from_entropy(float density, float entropy) {
* @param entropy The entropy \f$S\f$.
*/
__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
- float density, float entropy) {
+ const float density, const float entropy) {
return hydro_gamma_minus_one * eos.isothermal_internal_energy * density;
}
@@ -81,7 +81,7 @@ __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
* @return The entropy \f$A\f$.
*/
__attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
- float density, float pressure) {
+ const float density, const float pressure) {
return hydro_gamma_minus_one * eos.isothermal_internal_energy *
pow_minus_gamma_minus_one(density);
@@ -98,7 +98,7 @@ __attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
* @param entropy The entropy \f$S\f$.
*/
__attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
- float density, float entropy) {
+ const float density, const float entropy) {
return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
hydro_gamma_minus_one);
@@ -114,7 +114,7 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
* @param u The internal energy \f$u\f$
*/
__attribute__((always_inline)) INLINE static float
-gas_entropy_from_internal_energy(float density, float u) {
+gas_entropy_from_internal_energy(const float density, const float u) {
return hydro_gamma_minus_one * eos.isothermal_internal_energy *
pow_minus_gamma_minus_one(density);
@@ -130,7 +130,7 @@ gas_entropy_from_internal_energy(float density, float u) {
* @param u The internal energy \f$u\f$
*/
__attribute__((always_inline)) INLINE static float
-gas_pressure_from_internal_energy(float density, float u) {
+gas_pressure_from_internal_energy(const float density, const float u) {
return hydro_gamma_minus_one * eos.isothermal_internal_energy * density;
}
@@ -145,7 +145,7 @@ gas_pressure_from_internal_energy(float density, float u) {
* @return The internal energy \f$u\f$ (which is constant).
*/
__attribute__((always_inline)) INLINE static float
-gas_internal_energy_from_pressure(float density, float pressure) {
+gas_internal_energy_from_pressure(const float density, const float pressure) {
return eos.isothermal_internal_energy;
}
@@ -160,7 +160,7 @@ gas_internal_energy_from_pressure(float density, float pressure) {
* @param u The internal energy \f$u\f$
*/
__attribute__((always_inline)) INLINE static float
-gas_soundspeed_from_internal_energy(float density, float u) {
+gas_soundspeed_from_internal_energy(const float density, const float u) {
return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
hydro_gamma_minus_one);
@@ -177,7 +177,7 @@ gas_soundspeed_from_internal_energy(float density, float u) {
* @param P The pressure \f$P\f$
*/
__attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
- float density, float P) {
+ const float density, const float P) {
return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
hydro_gamma_minus_one);
diff --git a/src/equation_of_state/planetary/equation_of_state.h b/src/equation_of_state/planetary/equation_of_state.h
index a6d1a4e336a03de108dd6787edb9e66346c1bb43..4a46cadb80cb76fd72aba393026914b94252dbb5 100755
--- a/src/equation_of_state/planetary/equation_of_state.h
+++ b/src/equation_of_state/planetary/equation_of_state.h
@@ -102,6 +102,10 @@ enum eos_planetary_material_id {
eos_planetary_id_Til_basalt =
eos_planetary_type_Til * eos_planetary_type_factor + 3,
+ /*! Tillotson ice */
+ eos_planetary_id_Til_ice =
+ eos_planetary_type_Til * eos_planetary_type_factor + 4,
+
/* Hubbard & MacFarlane (1980) Uranus/Neptune */
/*! Hydrogen-helium atmosphere */
@@ -134,6 +138,23 @@ enum eos_planetary_material_id {
eos_planetary_id_SS08_water =
eos_planetary_type_SESAME * eos_planetary_type_factor + 3,
+ /*! AQUA (Haldemann et al. 2020) SESAME-like water */
+ eos_planetary_id_AQUA =
+ eos_planetary_type_SESAME * eos_planetary_type_factor + 4,
+
+ /*! CMS19 hydrogen (Chabrier et al. 2019) SESAME-like hydrogen */
+ eos_planetary_id_CMS19_H =
+ eos_planetary_type_SESAME * eos_planetary_type_factor + 5,
+
+ /*! CMS19 helium (Chabrier et al. 2019) SESAME-like helium */
+ eos_planetary_id_CMS19_He =
+ eos_planetary_type_SESAME * eos_planetary_type_factor + 6,
+
+ /*! CD21 hydrogen-helium (Chabrier & Debras 2021) SESAME-like H-He mixture
+ (Y=0.245) */
+ eos_planetary_id_CD21_HHe =
+ eos_planetary_type_SESAME * eos_planetary_type_factor + 7,
+
/* ANEOS */
/*! ANEOS forsterite (Stewart et al. 2019) -- in SESAME-style tables */
@@ -149,6 +170,10 @@ enum eos_planetary_material_id {
eos_planetary_type_ANEOS * eos_planetary_type_factor + 2,
};
+/* Base material ID for custom Tillotson EoS */
+#define eos_planetary_Til_custom_base_id \
+ (eos_planetary_type_Til * eos_planetary_type_factor + 90)
+
/* Individual EOS function headers. */
#include "hm80.h"
#include "ideal_gas.h"
@@ -160,9 +185,11 @@ enum eos_planetary_material_id {
*/
struct eos_parameters {
struct idg_params idg_def;
- struct Til_params Til_iron, Til_granite, Til_water, Til_basalt;
+ struct Til_params Til_iron, Til_granite, Til_water, Til_basalt, Til_ice,
+ Til_custom[10];
struct HM80_params HM80_HHe, HM80_ice, HM80_rock;
- struct SESAME_params SESAME_iron, SESAME_basalt, SESAME_water, SS08_water;
+ struct SESAME_params SESAME_iron, SESAME_basalt, SESAME_water, SS08_water,
+ AQUA, CMS19_H, CMS19_He, CD21_HHe;
struct SESAME_params ANEOS_forsterite, ANEOS_iron, ANEOS_Fe85Si15;
struct SESAME_params custom[10];
};
@@ -223,8 +250,20 @@ gas_internal_energy_from_entropy(float density, float entropy,
&eos.Til_basalt);
break;
+ case eos_planetary_id_Til_ice:
+ return Til_internal_energy_from_entropy(density, entropy,
+ &eos.Til_ice);
+ break;
+
default:
- return -1.f;
+ // Custom user-provided Tillotson
+ if (mat_id >= eos_planetary_Til_custom_base_id) {
+ const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+ return Til_internal_energy_from_entropy(density, entropy,
+ &eos.Til_custom[i_custom]);
+ } else {
+ return -1.f;
+ }
};
break;
@@ -278,6 +317,42 @@ gas_internal_energy_from_entropy(float density, float entropy,
&eos.SS08_water);
break;
+ case eos_planetary_id_AQUA:
+#ifdef SWIFT_DEBUG_CHECKS
+ if (eos.AQUA.mat_id != mat_id)
+ error("EoS not enabled. Please set EoS:planetary_use_AQUA: 1");
+#endif
+ return SESAME_internal_energy_from_entropy(density, entropy,
+ &eos.AQUA);
+ break;
+
+ case eos_planetary_id_CMS19_H:
+#ifdef SWIFT_DEBUG_CHECKS
+ if (eos.CMS19_H.mat_id != mat_id)
+ error("EoS not enabled. Please set EoS:planetary_use_CMS19_H: 1");
+#endif
+ return SESAME_internal_energy_from_entropy(density, entropy,
+ &eos.CMS19_H);
+ break;
+
+ case eos_planetary_id_CMS19_He:
+#ifdef SWIFT_DEBUG_CHECKS
+ if (eos.CMS19_He.mat_id != mat_id)
+ error("EoS not enabled. Please set EoS:planetary_use_CMS19_He: 1");
+#endif
+ return SESAME_internal_energy_from_entropy(density, entropy,
+ &eos.CMS19_He);
+ break;
+
+ case eos_planetary_id_CD21_HHe:
+#ifdef SWIFT_DEBUG_CHECKS
+ if (eos.CD21_HHe.mat_id != mat_id)
+ error("EoS not enabled. Please set EoS:planetary_use_CD21_HHe: 1");
+#endif
+ return SESAME_internal_energy_from_entropy(density, entropy,
+ &eos.CD21_HHe);
+ break;
+
default:
return -1.f;
};
@@ -372,8 +447,19 @@ __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
return Til_pressure_from_entropy(density, entropy, &eos.Til_basalt);
break;
+ case eos_planetary_id_Til_ice:
+ return Til_pressure_from_entropy(density, entropy, &eos.Til_ice);
+ break;
+
default:
- return -1.f;
+ // Custom user-provided Tillotson
+ if (mat_id >= eos_planetary_Til_custom_base_id) {
+ const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+ return Til_pressure_from_entropy(density, entropy,
+ &eos.Til_custom[i_custom]);
+ } else {
+ return -1.f;
+ }
};
break;
@@ -424,6 +510,22 @@ __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
&eos.SS08_water);
break;
+ case eos_planetary_id_AQUA:
+ return SESAME_pressure_from_entropy(density, entropy, &eos.AQUA);
+ break;
+
+ case eos_planetary_id_CMS19_H:
+ return SESAME_pressure_from_entropy(density, entropy, &eos.CMS19_H);
+ break;
+
+ case eos_planetary_id_CMS19_He:
+ return SESAME_pressure_from_entropy(density, entropy, &eos.CMS19_He);
+ break;
+
+ case eos_planetary_id_CD21_HHe:
+ return SESAME_pressure_from_entropy(density, entropy, &eos.CD21_HHe);
+ break;
+
default:
return -1.f;
};
@@ -519,8 +621,19 @@ __attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
return Til_entropy_from_pressure(density, P, &eos.Til_basalt);
break;
+ case eos_planetary_id_Til_ice:
+ return Til_entropy_from_pressure(density, P, &eos.Til_ice);
+ break;
+
default:
- return -1.f;
+ // Custom user-provided Tillotson
+ if (mat_id >= eos_planetary_Til_custom_base_id) {
+ const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+ return Til_entropy_from_pressure(density, P,
+ &eos.Til_custom[i_custom]);
+ } else {
+ return -1.f;
+ }
};
break;
@@ -567,6 +680,22 @@ __attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
return SESAME_entropy_from_pressure(density, P, &eos.SS08_water);
break;
+ case eos_planetary_id_AQUA:
+ return SESAME_entropy_from_pressure(density, P, &eos.AQUA);
+ break;
+
+ case eos_planetary_id_CMS19_H:
+ return SESAME_entropy_from_pressure(density, P, &eos.CMS19_H);
+ break;
+
+ case eos_planetary_id_CMS19_He:
+ return SESAME_entropy_from_pressure(density, P, &eos.CMS19_He);
+ break;
+
+ case eos_planetary_id_CD21_HHe:
+ return SESAME_entropy_from_pressure(density, P, &eos.CD21_HHe);
+ break;
+
default:
return -1.f;
};
@@ -659,8 +788,19 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
return Til_soundspeed_from_entropy(density, entropy, &eos.Til_basalt);
break;
+ case eos_planetary_id_Til_ice:
+ return Til_soundspeed_from_entropy(density, entropy, &eos.Til_ice);
+ break;
+
default:
- return -1.f;
+ // Custom user-provided Tillotson
+ if (mat_id >= eos_planetary_Til_custom_base_id) {
+ const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+ return Til_soundspeed_from_entropy(density, entropy,
+ &eos.Til_custom[i_custom]);
+ } else {
+ return -1.f;
+ }
};
break;
@@ -711,6 +851,24 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
&eos.SS08_water);
break;
+ case eos_planetary_id_AQUA:
+ return SESAME_soundspeed_from_entropy(density, entropy, &eos.AQUA);
+ break;
+
+ case eos_planetary_id_CMS19_H:
+ return SESAME_soundspeed_from_entropy(density, entropy, &eos.CMS19_H);
+ break;
+
+ case eos_planetary_id_CMS19_He:
+ return SESAME_soundspeed_from_entropy(density, entropy,
+ &eos.CMS19_He);
+ break;
+
+ case eos_planetary_id_CD21_HHe:
+ return SESAME_soundspeed_from_entropy(density, entropy,
+ &eos.CD21_HHe);
+ break;
+
default:
return -1.f;
};
@@ -805,8 +963,19 @@ gas_entropy_from_internal_energy(float density, float u,
return Til_entropy_from_internal_energy(density, u, &eos.Til_basalt);
break;
+ case eos_planetary_id_Til_ice:
+ return Til_entropy_from_internal_energy(density, u, &eos.Til_ice);
+ break;
+
default:
- return -1.f;
+ // Custom user-provided Tillotson
+ if (mat_id >= eos_planetary_Til_custom_base_id) {
+ const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+ return Til_entropy_from_internal_energy(density, u,
+ &eos.Til_custom[i_custom]);
+ } else {
+ return -1.f;
+ }
};
break;
@@ -857,6 +1026,22 @@ gas_entropy_from_internal_energy(float density, float u,
&eos.SS08_water);
break;
+ case eos_planetary_id_AQUA:
+ return SESAME_entropy_from_internal_energy(density, u, &eos.AQUA);
+ break;
+
+ case eos_planetary_id_CMS19_H:
+ return SESAME_entropy_from_internal_energy(density, u, &eos.CMS19_H);
+ break;
+
+ case eos_planetary_id_CMS19_He:
+ return SESAME_entropy_from_internal_energy(density, u, &eos.CMS19_He);
+ break;
+
+ case eos_planetary_id_CD21_HHe:
+ return SESAME_entropy_from_internal_energy(density, u, &eos.CD21_HHe);
+ break;
+
default:
return -1.f;
};
@@ -978,11 +1163,26 @@ gas_pressure_from_internal_energy(float density, float u,
return Til_pressure_from_internal_energy(density, u, &eos.Til_basalt);
break;
+ case eos_planetary_id_Til_ice:
+#ifdef SWIFT_DEBUG_CHECKS
+ if (eos.Til_ice.mat_id != mat_id)
+ error("EoS not enabled. Please set EoS:planetary_use_Til_ice: 1");
+#endif
+ return Til_pressure_from_internal_energy(density, u, &eos.Til_ice);
+ break;
+
default:
+ // Custom user-provided Tillotson
+ if (mat_id >= eos_planetary_Til_custom_base_id) {
+ const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+ return Til_pressure_from_internal_energy(density, u,
+ &eos.Til_custom[i_custom]);
+ } else {
#ifdef SWIFT_DEBUG_CHECKS
- error("Unknown material ID! mat_id = %d", mat_id);
+ error("Unknown material ID! mat_id = %d", mat_id);
#endif
- return -1.f;
+ return -1.f;
+ }
};
break;
@@ -1070,6 +1270,40 @@ gas_pressure_from_internal_energy(float density, float u,
&eos.SS08_water);
break;
+ case eos_planetary_id_AQUA:
+#ifdef SWIFT_DEBUG_CHECKS
+ if (eos.AQUA.mat_id != mat_id)
+ error("EoS not enabled. Please set EoS:planetary_use_AQUA: 1");
+#endif
+ return SESAME_pressure_from_internal_energy(density, u, &eos.AQUA);
+ break;
+
+ case eos_planetary_id_CMS19_H:
+#ifdef SWIFT_DEBUG_CHECKS
+ if (eos.CMS19_H.mat_id != mat_id)
+ error("EoS not enabled. Please set EoS:planetary_use_CMS19_H: 1");
+#endif
+ return SESAME_pressure_from_internal_energy(density, u, &eos.CMS19_H);
+ break;
+
+ case eos_planetary_id_CMS19_He:
+#ifdef SWIFT_DEBUG_CHECKS
+ if (eos.CMS19_He.mat_id != mat_id)
+ error("EoS not enabled. Please set EoS:planetary_use_CMS19_He: 1");
+#endif
+ return SESAME_pressure_from_internal_energy(density, u,
+ &eos.CMS19_He);
+ break;
+
+ case eos_planetary_id_CD21_HHe:
+#ifdef SWIFT_DEBUG_CHECKS
+ if (eos.CD21_HHe.mat_id != mat_id)
+ error("EoS not enabled. Please set EoS:planetary_use_CD21_HHe: 1");
+#endif
+ return SESAME_pressure_from_internal_energy(density, u,
+ &eos.CD21_HHe);
+ break;
+
default:
#ifdef SWIFT_DEBUG_CHECKS
error("Unknown material ID! mat_id = %d", mat_id);
@@ -1200,8 +1434,19 @@ gas_internal_energy_from_pressure(float density, float P,
return Til_internal_energy_from_pressure(density, P, &eos.Til_basalt);
break;
+ case eos_planetary_id_Til_ice:
+ return Til_internal_energy_from_pressure(density, P, &eos.Til_ice);
+ break;
+
default:
- return -1.f;
+ // Custom user-provided Tillotson
+ if (mat_id >= eos_planetary_Til_custom_base_id) {
+ const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+ return Til_internal_energy_from_pressure(density, P,
+ &eos.Til_custom[i_custom]);
+ } else {
+ return -1.f;
+ }
};
break;
@@ -1252,6 +1497,24 @@ gas_internal_energy_from_pressure(float density, float P,
&eos.SS08_water);
break;
+ case eos_planetary_id_AQUA:
+ return SESAME_internal_energy_from_pressure(density, P, &eos.AQUA);
+ break;
+
+ case eos_planetary_id_CMS19_H:
+ return SESAME_internal_energy_from_pressure(density, P, &eos.CMS19_H);
+ break;
+
+ case eos_planetary_id_CMS19_He:
+ return SESAME_internal_energy_from_pressure(density, P,
+ &eos.CMS19_He);
+ break;
+
+ case eos_planetary_id_CD21_HHe:
+ return SESAME_internal_energy_from_pressure(density, P,
+ &eos.CD21_HHe);
+ break;
+
default:
return -1.f;
};
@@ -1350,8 +1613,19 @@ gas_soundspeed_from_internal_energy(float density, float u,
&eos.Til_basalt);
break;
+ case eos_planetary_id_Til_ice:
+ return Til_soundspeed_from_internal_energy(density, u, &eos.Til_ice);
+ break;
+
default:
- return -1.f;
+ // Custom user-provided Tillotson
+ if (mat_id >= eos_planetary_Til_custom_base_id) {
+ const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+ return Til_soundspeed_from_internal_energy(
+ density, u, &eos.Til_custom[i_custom]);
+ } else {
+ return -1.f;
+ }
};
break;
@@ -1405,6 +1679,25 @@ gas_soundspeed_from_internal_energy(float density, float u,
&eos.SS08_water);
break;
+ case eos_planetary_id_AQUA:
+ return SESAME_soundspeed_from_internal_energy(density, u, &eos.AQUA);
+ break;
+
+ case eos_planetary_id_CMS19_H:
+ return SESAME_soundspeed_from_internal_energy(density, u,
+ &eos.CMS19_H);
+ break;
+
+ case eos_planetary_id_CMS19_He:
+ return SESAME_soundspeed_from_internal_energy(density, u,
+ &eos.CMS19_He);
+ break;
+
+ case eos_planetary_id_CD21_HHe:
+ return SESAME_soundspeed_from_internal_energy(density, u,
+ &eos.CD21_HHe);
+ break;
+
default:
return -1.f;
};
@@ -1499,8 +1792,19 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
return Til_soundspeed_from_pressure(density, P, &eos.Til_basalt);
break;
+ case eos_planetary_id_Til_ice:
+ return Til_soundspeed_from_pressure(density, P, &eos.Til_ice);
+ break;
+
default:
- return -1.f;
+ // Custom user-provided Tillotson
+ if (mat_id >= eos_planetary_Til_custom_base_id) {
+ const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+ return Til_soundspeed_from_pressure(density, P,
+ &eos.Til_custom[i_custom]);
+ } else {
+ return -1.f;
+ }
};
break;
@@ -1548,6 +1852,22 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
return SESAME_soundspeed_from_pressure(density, P, &eos.SS08_water);
break;
+ case eos_planetary_id_AQUA:
+ return SESAME_soundspeed_from_pressure(density, P, &eos.AQUA);
+ break;
+
+ case eos_planetary_id_CMS19_H:
+ return SESAME_soundspeed_from_pressure(density, P, &eos.CMS19_H);
+ break;
+
+ case eos_planetary_id_CMS19_He:
+ return SESAME_soundspeed_from_pressure(density, P, &eos.CMS19_He);
+ break;
+
+ case eos_planetary_id_CD21_HHe:
+ return SESAME_soundspeed_from_pressure(density, P, &eos.CD21_HHe);
+ break;
+
default:
return -1.f;
};
@@ -1591,106 +1911,740 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
}
/**
- * @brief Initialize the eos parameters
+ * @brief Returns the temperature given density and internal energy
*
- * @param e The #eos_parameters
- * @param params The parsed parameters
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
*/
-__attribute__((always_inline)) INLINE static void eos_init(
- struct eos_parameters *e, const struct phys_const *phys_const,
- const struct unit_system *us, struct swift_params *params) {
+__attribute__((always_inline)) INLINE static float
+gas_temperature_from_internal_energy(float density, float u,
+ enum eos_planetary_material_id mat_id) {
+ const enum eos_planetary_type_id type =
+ (enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
- // Prepare any/all requested EoS: Set the parameters and material IDs, load
- // tables etc., and convert to internal units
+ /* Select the material base type */
+ switch (type) {
- // Ideal gas
- if (parser_get_opt_param_int(params, "EoS:planetary_use_idg_def", 0)) {
- set_idg_def(&e->idg_def, eos_planetary_id_idg_def);
- }
+ /* Ideal gas EoS */
+ case eos_planetary_type_idg:
- // Tillotson
- if (parser_get_opt_param_int(params, "EoS:planetary_use_Til_iron", 0)) {
- set_Til_iron(&e->Til_iron, eos_planetary_id_Til_iron);
- convert_units_Til(&e->Til_iron, us);
- }
- if (parser_get_opt_param_int(params, "EoS:planetary_use_Til_granite", 0)) {
- set_Til_granite(&e->Til_granite, eos_planetary_id_Til_granite);
- convert_units_Til(&e->Til_granite, us);
- }
- if (parser_get_opt_param_int(params, "EoS:planetary_use_Til_water", 0)) {
- set_Til_water(&e->Til_water, eos_planetary_id_Til_water);
- convert_units_Til(&e->Til_water, us);
- }
- if (parser_get_opt_param_int(params, "EoS:planetary_use_Til_basalt", 0)) {
- set_Til_basalt(&e->Til_basalt, eos_planetary_id_Til_basalt);
- convert_units_Til(&e->Til_basalt, us);
- }
+ /* Select the material of this type */
+ switch (mat_id) {
+ case eos_planetary_id_idg_def:
+ return idg_temperature_from_internal_energy(density, u, &eos.idg_def);
+ break;
- // Hubbard & MacFarlane (1980)
- if (parser_get_opt_param_int(params, "EoS:planetary_use_HM80_HHe", 0)) {
- char HM80_HHe_table_file[PARSER_MAX_LINE_SIZE];
- set_HM80_HHe(&e->HM80_HHe, eos_planetary_id_HM80_HHe);
- parser_get_param_string(params, "EoS:planetary_HM80_HHe_table_file",
- HM80_HHe_table_file);
- load_table_HM80(&e->HM80_HHe, HM80_HHe_table_file);
- prepare_table_HM80(&e->HM80_HHe);
- convert_units_HM80(&e->HM80_HHe, us);
- }
- if (parser_get_opt_param_int(params, "EoS:planetary_use_HM80_ice", 0)) {
- char HM80_ice_table_file[PARSER_MAX_LINE_SIZE];
- set_HM80_ice(&e->HM80_ice, eos_planetary_id_HM80_ice);
- parser_get_param_string(params, "EoS:planetary_HM80_ice_table_file",
- HM80_ice_table_file);
- load_table_HM80(&e->HM80_ice, HM80_ice_table_file);
- prepare_table_HM80(&e->HM80_ice);
- convert_units_HM80(&e->HM80_ice, us);
- }
- if (parser_get_opt_param_int(params, "EoS:planetary_use_HM80_rock", 0)) {
- char HM80_rock_table_file[PARSER_MAX_LINE_SIZE];
- set_HM80_rock(&e->HM80_rock, eos_planetary_id_HM80_rock);
- parser_get_param_string(params, "EoS:planetary_HM80_rock_table_file",
- HM80_rock_table_file);
- load_table_HM80(&e->HM80_rock, HM80_rock_table_file);
- prepare_table_HM80(&e->HM80_rock);
- convert_units_HM80(&e->HM80_rock, us);
- }
+ default:
+#ifdef SWIFT_DEBUG_CHECKS
+ error("Unknown material ID! mat_id = %d", mat_id);
+#endif
+ return -1.f;
+ };
+ break;
- // SESAME
- if (parser_get_opt_param_int(params, "EoS:planetary_use_SESAME_iron", 0)) {
- char SESAME_iron_table_file[PARSER_MAX_LINE_SIZE];
- set_SESAME_iron(&e->SESAME_iron, eos_planetary_id_SESAME_iron);
- parser_get_param_string(params, "EoS:planetary_SESAME_iron_table_file",
- SESAME_iron_table_file);
- load_table_SESAME(&e->SESAME_iron, SESAME_iron_table_file);
- prepare_table_SESAME(&e->SESAME_iron);
- convert_units_SESAME(&e->SESAME_iron, us);
- }
- if (parser_get_opt_param_int(params, "EoS:planetary_use_SESAME_basalt", 0)) {
- char SESAME_basalt_table_file[PARSER_MAX_LINE_SIZE];
- set_SESAME_basalt(&e->SESAME_basalt, eos_planetary_id_SESAME_basalt);
- parser_get_param_string(params, "EoS:planetary_SESAME_basalt_table_file",
- SESAME_basalt_table_file);
- load_table_SESAME(&e->SESAME_basalt, SESAME_basalt_table_file);
- prepare_table_SESAME(&e->SESAME_basalt);
- convert_units_SESAME(&e->SESAME_basalt, us);
- }
- if (parser_get_opt_param_int(params, "EoS:planetary_use_SESAME_water", 0)) {
- char SESAME_water_table_file[PARSER_MAX_LINE_SIZE];
- set_SESAME_water(&e->SESAME_water, eos_planetary_id_SESAME_water);
- parser_get_param_string(params, "EoS:planetary_SESAME_water_table_file",
- SESAME_water_table_file);
- load_table_SESAME(&e->SESAME_water, SESAME_water_table_file);
- prepare_table_SESAME(&e->SESAME_water);
- convert_units_SESAME(&e->SESAME_water, us);
- }
- if (parser_get_opt_param_int(params, "EoS:planetary_use_SS08_water", 0)) {
- char SS08_water_table_file[PARSER_MAX_LINE_SIZE];
- set_SS08_water(&e->SS08_water, eos_planetary_id_SS08_water);
- parser_get_param_string(params, "EoS:planetary_SS08_water_table_file",
- SS08_water_table_file);
- load_table_SESAME(&e->SS08_water, SS08_water_table_file);
- prepare_table_SESAME(&e->SS08_water);
- convert_units_SESAME(&e->SS08_water, us);
+ /* Tillotson EoS */
+ case eos_planetary_type_Til:
+
+ /* Select the material of this type */
+ switch (mat_id) {
+ case eos_planetary_id_Til_iron:
+ return Til_temperature_from_internal_energy(density, u,
+ &eos.Til_iron);
+ break;
+
+ case eos_planetary_id_Til_granite:
+ return Til_temperature_from_internal_energy(density, u,
+ &eos.Til_granite);
+ break;
+
+ case eos_planetary_id_Til_water:
+ return Til_temperature_from_internal_energy(density, u,
+ &eos.Til_water);
+ break;
+
+ case eos_planetary_id_Til_basalt:
+ return Til_temperature_from_internal_energy(density, u,
+ &eos.Til_basalt);
+ break;
+
+ case eos_planetary_id_Til_ice:
+ return Til_temperature_from_internal_energy(density, u, &eos.Til_ice);
+ break;
+
+ default:
+ // Custom user-provided Tillotson
+ if (mat_id >= eos_planetary_Til_custom_base_id) {
+ const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+ return Til_temperature_from_internal_energy(
+ density, u, &eos.Til_custom[i_custom]);
+ } else {
+ return -1.f;
+ }
+ };
+ break;
+
+ /* Hubbard & MacFarlane (1980) EoS */
+ case eos_planetary_type_HM80:
+
+ /* Select the material of this type */
+ switch (mat_id) {
+ case eos_planetary_id_HM80_HHe:
+ return HM80_temperature_from_internal_energy(density, u,
+ &eos.HM80_HHe);
+ break;
+
+ case eos_planetary_id_HM80_ice:
+ return HM80_temperature_from_internal_energy(density, u,
+ &eos.HM80_ice);
+ break;
+
+ case eos_planetary_id_HM80_rock:
+ return HM80_temperature_from_internal_energy(density, u,
+ &eos.HM80_rock);
+ break;
+
+ default:
+#ifdef SWIFT_DEBUG_CHECKS
+ error("Unknown material ID! mat_id = %d", mat_id);
+#endif
+ return -1.f;
+ };
+ break;
+
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
+
+ /* Select the material of this type */
+ switch (mat_id) {
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_temperature_from_internal_energy(density, u,
+ &eos.SESAME_iron);
+ break;
+
+ case eos_planetary_id_SESAME_basalt:
+ return SESAME_temperature_from_internal_energy(density, u,
+ &eos.SESAME_basalt);
+ break;
+
+ case eos_planetary_id_SESAME_water:
+ return SESAME_temperature_from_internal_energy(density, u,
+ &eos.SESAME_water);
+ break;
+
+ case eos_planetary_id_SS08_water:
+ return SESAME_temperature_from_internal_energy(density, u,
+ &eos.SS08_water);
+ break;
+
+ case eos_planetary_id_AQUA:
+ return SESAME_temperature_from_internal_energy(density, u, &eos.AQUA);
+ break;
+
+ case eos_planetary_id_CMS19_H:
+ return SESAME_temperature_from_internal_energy(density, u,
+ &eos.CMS19_H);
+ break;
+
+ case eos_planetary_id_CMS19_He:
+ return SESAME_temperature_from_internal_energy(density, u,
+ &eos.CMS19_He);
+ break;
+
+ case eos_planetary_id_CD21_HHe:
+ return SESAME_temperature_from_internal_energy(density, u,
+ &eos.CD21_HHe);
+ break;
+
+ default:
+ return -1.f;
+ };
+ break;
+
+ /* ANEOS -- using SESAME-style tables */
+ case eos_planetary_type_ANEOS:;
+
+ /* Select the material of this type */
+ switch (mat_id) {
+ case eos_planetary_id_ANEOS_forsterite:
+ return SESAME_temperature_from_internal_energy(density, u,
+ &eos.ANEOS_forsterite);
+ break;
+
+ case eos_planetary_id_ANEOS_iron:
+ return SESAME_temperature_from_internal_energy(density, u,
+ &eos.ANEOS_iron);
+ break;
+
+ case eos_planetary_id_ANEOS_Fe85Si15:
+ return SESAME_temperature_from_internal_energy(density, u,
+ &eos.ANEOS_Fe85Si15);
+ break;
+
+ default:
+ return -1.f;
+ };
+ break;
+
+ /*! Generic user-provided custom tables */
+ case eos_planetary_type_custom: {
+ const int i_custom =
+ mat_id - eos_planetary_type_custom * eos_planetary_type_factor;
+ return SESAME_temperature_from_internal_energy(density, u,
+ &eos.custom[i_custom]);
+ break;
+ }
+
+ default:
+ return -1.f;
+ }
+}
+
+/**
+ * @brief Returns the density given pressure and temperature
+ *
+ * @param P The pressure \f$P\f$
+ * @param T The temperature \f$T\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_density_from_pressure_and_temperature(
+ float P, float T, enum eos_planetary_material_id mat_id) {
+ const enum eos_planetary_type_id type =
+ (enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
+
+ /* Select the material base type */
+ switch (type) {
+
+ /* Ideal gas EoS */
+ case eos_planetary_type_idg:
+
+ /* Select the material of this type */
+ switch (mat_id) {
+ case eos_planetary_id_idg_def:
+ return idg_density_from_pressure_and_temperature(P, T, &eos.idg_def);
+ break;
+
+ default:
+#ifdef SWIFT_DEBUG_CHECKS
+ error("Unknown material ID! mat_id = %d", mat_id);
+#endif
+ return -1.f;
+ };
+ break;
+
+ /* Tillotson EoS */
+ case eos_planetary_type_Til:
+
+ /* Select the material of this type */
+ switch (mat_id) {
+ case eos_planetary_id_Til_iron:
+ return Til_density_from_pressure_and_temperature(P, T, &eos.Til_iron);
+ break;
+
+ case eos_planetary_id_Til_granite:
+ return Til_density_from_pressure_and_temperature(P, T,
+ &eos.Til_granite);
+ break;
+
+ case eos_planetary_id_Til_water:
+ return Til_density_from_pressure_and_temperature(P, T,
+ &eos.Til_water);
+ break;
+
+ case eos_planetary_id_Til_basalt:
+ return Til_density_from_pressure_and_temperature(P, T,
+ &eos.Til_basalt);
+ break;
+
+ case eos_planetary_id_Til_ice:
+ return Til_density_from_pressure_and_temperature(P, T, &eos.Til_ice);
+ break;
+
+ default:
+ // Custom user-provided Tillotson
+ if (mat_id >= eos_planetary_Til_custom_base_id) {
+ const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+ return Til_density_from_pressure_and_temperature(
+ P, T, &eos.Til_custom[i_custom]);
+ } else {
+ return -1.f;
+ }
+ };
+ break;
+
+ /* Hubbard & MacFarlane (1980) EoS */
+ case eos_planetary_type_HM80:
+
+ /* Select the material of this type */
+ switch (mat_id) {
+ case eos_planetary_id_HM80_HHe:
+ return HM80_density_from_pressure_and_temperature(P, T,
+ &eos.HM80_HHe);
+ break;
+
+ case eos_planetary_id_HM80_ice:
+ return HM80_density_from_pressure_and_temperature(P, T,
+ &eos.HM80_ice);
+ break;
+
+ case eos_planetary_id_HM80_rock:
+ return HM80_density_from_pressure_and_temperature(P, T,
+ &eos.HM80_rock);
+ break;
+
+ default:
+#ifdef SWIFT_DEBUG_CHECKS
+ error("Unknown material ID! mat_id = %d", mat_id);
+#endif
+ return -1.f;
+ };
+ break;
+
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
+
+ /* Select the material of this type */
+ switch (mat_id) {
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_density_from_pressure_and_temperature(P, T,
+ &eos.SESAME_iron);
+ break;
+
+ case eos_planetary_id_SESAME_basalt:
+ return SESAME_density_from_pressure_and_temperature(
+ P, T, &eos.SESAME_basalt);
+ break;
+
+ case eos_planetary_id_SESAME_water:
+ return SESAME_density_from_pressure_and_temperature(
+ P, T, &eos.SESAME_water);
+ break;
+
+ case eos_planetary_id_SS08_water:
+ return SESAME_density_from_pressure_and_temperature(P, T,
+ &eos.SS08_water);
+ break;
+
+ case eos_planetary_id_AQUA:
+ return SESAME_density_from_pressure_and_temperature(P, T, &eos.AQUA);
+ break;
+
+ case eos_planetary_id_CMS19_H:
+ return SESAME_density_from_pressure_and_temperature(P, T,
+ &eos.CMS19_H);
+ break;
+
+ case eos_planetary_id_CMS19_He:
+ return SESAME_density_from_pressure_and_temperature(P, T,
+ &eos.CMS19_He);
+ break;
+
+ case eos_planetary_id_CD21_HHe:
+ return SESAME_density_from_pressure_and_temperature(P, T,
+ &eos.CD21_HHe);
+ break;
+
+ default:
+ return -1.f;
+ };
+ break;
+
+ /* ANEOS -- using SESAME-style tables */
+ case eos_planetary_type_ANEOS:;
+
+ /* Select the material of this type */
+ switch (mat_id) {
+ case eos_planetary_id_ANEOS_forsterite:
+ return SESAME_density_from_pressure_and_temperature(
+ P, T, &eos.ANEOS_forsterite);
+ break;
+
+ case eos_planetary_id_ANEOS_iron:
+ return SESAME_density_from_pressure_and_temperature(P, T,
+ &eos.ANEOS_iron);
+ break;
+
+ case eos_planetary_id_ANEOS_Fe85Si15:
+ return SESAME_density_from_pressure_and_temperature(
+ P, T, &eos.ANEOS_Fe85Si15);
+ break;
+
+ default:
+ return -1.f;
+ };
+ break;
+
+ /*! Generic user-provided custom tables */
+ case eos_planetary_type_custom: {
+ const int i_custom =
+ mat_id - eos_planetary_type_custom * eos_planetary_type_factor;
+ return SESAME_density_from_pressure_and_temperature(
+ P, T, &eos.custom[i_custom]);
+ break;
+ }
+
+ default:
+ return -1.f;
+ }
+}
+
+/**
+ * @brief Returns the density given pressure and internal energy
+ *
+ * @param P The pressure \f$P\f$
+ * @param T The temperature \f$T\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_density_from_pressure_and_internal_energy(
+ float P, float u, float rho_ref, float rho_sph,
+ enum eos_planetary_material_id mat_id) {
+ const enum eos_planetary_type_id type =
+ (enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
+
+ /* Select the material base type */
+ switch (type) {
+
+ /* Ideal gas EoS */
+ case eos_planetary_type_idg:
+
+ /* Select the material of this type */
+ switch (mat_id) {
+ case eos_planetary_id_idg_def:
+ return idg_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.idg_def);
+ break;
+
+ default:
+#ifdef SWIFT_DEBUG_CHECKS
+ error("Unknown material ID! mat_id = %d", mat_id);
+#endif
+ return -1.f;
+ };
+ break;
+
+ /* Tillotson EoS */
+ case eos_planetary_type_Til:
+
+ /* Select the material of this type */
+ switch (mat_id) {
+ case eos_planetary_id_Til_iron:
+ return Til_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.Til_iron);
+ break;
+
+ case eos_planetary_id_Til_granite:
+ return Til_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.Til_granite);
+ break;
+
+ case eos_planetary_id_Til_water:
+ return Til_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.Til_water);
+ break;
+
+ case eos_planetary_id_Til_basalt:
+ return Til_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.Til_basalt);
+ break;
+
+ case eos_planetary_id_Til_ice:
+ return Til_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.Til_ice);
+ break;
+
+ default:
+ // Custom user-provided Tillotson
+ if (mat_id >= eos_planetary_Til_custom_base_id) {
+ const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+ return Til_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.Til_custom[i_custom]);
+ } else {
+ return -1.f;
+ }
+ };
+ break;
+
+ /* Hubbard & MacFarlane (1980) EoS */
+ case eos_planetary_type_HM80:
+
+ /* Select the material of this type */
+ switch (mat_id) {
+ case eos_planetary_id_HM80_HHe:
+ return HM80_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.HM80_HHe);
+ break;
+
+ case eos_planetary_id_HM80_ice:
+ return HM80_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.HM80_ice);
+ break;
+
+ case eos_planetary_id_HM80_rock:
+ return HM80_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.HM80_rock);
+ break;
+
+ default:
+#ifdef SWIFT_DEBUG_CHECKS
+ error("Unknown material ID! mat_id = %d", mat_id);
+#endif
+ return -1.f;
+ };
+ break;
+
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
+
+ /* Select the material of this type */
+ switch (mat_id) {
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.SESAME_iron);
+ break;
+
+ case eos_planetary_id_SESAME_basalt:
+ return SESAME_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.SESAME_basalt);
+ break;
+
+ case eos_planetary_id_SESAME_water:
+ return SESAME_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.SESAME_water);
+ break;
+
+ case eos_planetary_id_SS08_water:
+ return SESAME_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.SS08_water);
+ break;
+
+ case eos_planetary_id_AQUA:
+ return SESAME_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.AQUA);
+ break;
+
+ case eos_planetary_id_CMS19_H:
+ return SESAME_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.CMS19_H);
+ break;
+
+ case eos_planetary_id_CMS19_He:
+ return SESAME_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.CMS19_He);
+ break;
+
+ case eos_planetary_id_CD21_HHe:
+ return SESAME_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.CD21_HHe);
+ break;
+
+ default:
+ return -1.f;
+ };
+ break;
+
+ /* ANEOS -- using SESAME-style tables */
+ case eos_planetary_type_ANEOS:;
+
+ /* Select the material of this type */
+ switch (mat_id) {
+ case eos_planetary_id_ANEOS_forsterite:
+ return SESAME_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.ANEOS_forsterite);
+ break;
+
+ case eos_planetary_id_ANEOS_iron:
+ return SESAME_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.ANEOS_iron);
+ break;
+
+ case eos_planetary_id_ANEOS_Fe85Si15:
+ return SESAME_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.ANEOS_Fe85Si15);
+ break;
+
+ default:
+ return -1.f;
+ };
+ break;
+
+ /*! Generic user-provided custom tables */
+ case eos_planetary_type_custom: {
+ const int i_custom =
+ mat_id - eos_planetary_type_custom * eos_planetary_type_factor;
+ return SESAME_density_from_pressure_and_internal_energy(
+ P, u, rho_ref, rho_sph, &eos.custom[i_custom]);
+ break;
+ }
+
+ default:
+ return -1.f;
+ }
+}
+
+/**
+ * @brief Initialize the eos parameters
+ *
+ * @param e The #eos_parameters
+ * @param params The parsed parameters
+ */
+__attribute__((always_inline)) INLINE static void eos_init(
+ struct eos_parameters *e, const struct phys_const *phys_const,
+ const struct unit_system *us, struct swift_params *params) {
+
+ // Prepare any/all requested EoS: Set the parameters and material IDs, load
+ // tables etc., and convert to internal units
+
+ // Ideal gas
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_idg_def", 0)) {
+ set_idg_def(&e->idg_def, eos_planetary_id_idg_def);
+ }
+
+ // Tillotson
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_Til_iron", 0)) {
+ set_Til_iron(&e->Til_iron, eos_planetary_id_Til_iron);
+ set_Til_u_cold(&e->Til_iron, eos_planetary_id_Til_iron);
+ convert_units_Til(&e->Til_iron, us);
+ }
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_Til_granite", 0)) {
+ set_Til_granite(&e->Til_granite, eos_planetary_id_Til_granite);
+ set_Til_u_cold(&e->Til_granite, eos_planetary_id_Til_granite);
+ convert_units_Til(&e->Til_granite, us);
+ }
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_Til_water", 0)) {
+ set_Til_water(&e->Til_water, eos_planetary_id_Til_water);
+ set_Til_u_cold(&e->Til_water, eos_planetary_id_Til_water);
+ convert_units_Til(&e->Til_water, us);
+ }
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_Til_basalt", 0)) {
+ set_Til_basalt(&e->Til_basalt, eos_planetary_id_Til_basalt);
+ set_Til_u_cold(&e->Til_basalt, eos_planetary_id_Til_basalt);
+ convert_units_Til(&e->Til_basalt, us);
+ }
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_Til_ice", 0)) {
+ set_Til_ice(&e->Til_ice, eos_planetary_id_Til_ice);
+ set_Til_u_cold(&e->Til_ice, eos_planetary_id_Til_ice);
+ convert_units_Til(&e->Til_ice, us);
+ }
+
+ // Custom user-provided Tillotson
+ for (int i_custom = 0; i_custom <= 9; i_custom++) {
+ char param_name[PARSER_MAX_LINE_SIZE];
+ sprintf(param_name, "EoS:planetary_use_Til_custom_%d", i_custom);
+ if (parser_get_opt_param_int(params, param_name, 0)) {
+ char Til_custom_file[PARSER_MAX_LINE_SIZE];
+ int mat_id = eos_planetary_Til_custom_base_id + i_custom;
+
+ sprintf(param_name, "EoS:planetary_Til_custom_%d_param_file", i_custom);
+ parser_get_param_string(params, param_name, Til_custom_file);
+
+ set_Til_custom(&e->Til_custom[i_custom],
+ (enum eos_planetary_material_id)mat_id, Til_custom_file);
+ set_Til_u_cold(&e->Til_custom[i_custom],
+ (enum eos_planetary_material_id)mat_id);
+ convert_units_Til(&e->Til_custom[i_custom], us);
+ }
+ }
+
+ // Hubbard & MacFarlane (1980)
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_HM80_HHe", 0)) {
+ char HM80_HHe_table_file[PARSER_MAX_LINE_SIZE];
+ set_HM80_HHe(&e->HM80_HHe, eos_planetary_id_HM80_HHe);
+ parser_get_param_string(params, "EoS:planetary_HM80_HHe_table_file",
+ HM80_HHe_table_file);
+ load_table_HM80(&e->HM80_HHe, HM80_HHe_table_file);
+ prepare_table_HM80(&e->HM80_HHe);
+ convert_units_HM80(&e->HM80_HHe, us);
+ }
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_HM80_ice", 0)) {
+ char HM80_ice_table_file[PARSER_MAX_LINE_SIZE];
+ set_HM80_ice(&e->HM80_ice, eos_planetary_id_HM80_ice);
+ parser_get_param_string(params, "EoS:planetary_HM80_ice_table_file",
+ HM80_ice_table_file);
+ load_table_HM80(&e->HM80_ice, HM80_ice_table_file);
+ prepare_table_HM80(&e->HM80_ice);
+ convert_units_HM80(&e->HM80_ice, us);
+ }
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_HM80_rock", 0)) {
+ char HM80_rock_table_file[PARSER_MAX_LINE_SIZE];
+ set_HM80_rock(&e->HM80_rock, eos_planetary_id_HM80_rock);
+ parser_get_param_string(params, "EoS:planetary_HM80_rock_table_file",
+ HM80_rock_table_file);
+ load_table_HM80(&e->HM80_rock, HM80_rock_table_file);
+ prepare_table_HM80(&e->HM80_rock);
+ convert_units_HM80(&e->HM80_rock, us);
+ }
+
+ // SESAME
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_SESAME_iron", 0)) {
+ char SESAME_iron_table_file[PARSER_MAX_LINE_SIZE];
+ set_SESAME_iron(&e->SESAME_iron, eos_planetary_id_SESAME_iron);
+ parser_get_param_string(params, "EoS:planetary_SESAME_iron_table_file",
+ SESAME_iron_table_file);
+ load_table_SESAME(&e->SESAME_iron, SESAME_iron_table_file);
+ prepare_table_SESAME(&e->SESAME_iron);
+ convert_units_SESAME(&e->SESAME_iron, us);
+ }
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_SESAME_basalt", 0)) {
+ char SESAME_basalt_table_file[PARSER_MAX_LINE_SIZE];
+ set_SESAME_basalt(&e->SESAME_basalt, eos_planetary_id_SESAME_basalt);
+ parser_get_param_string(params, "EoS:planetary_SESAME_basalt_table_file",
+ SESAME_basalt_table_file);
+ load_table_SESAME(&e->SESAME_basalt, SESAME_basalt_table_file);
+ prepare_table_SESAME(&e->SESAME_basalt);
+ convert_units_SESAME(&e->SESAME_basalt, us);
+ }
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_SESAME_water", 0)) {
+ char SESAME_water_table_file[PARSER_MAX_LINE_SIZE];
+ set_SESAME_water(&e->SESAME_water, eos_planetary_id_SESAME_water);
+ parser_get_param_string(params, "EoS:planetary_SESAME_water_table_file",
+ SESAME_water_table_file);
+ load_table_SESAME(&e->SESAME_water, SESAME_water_table_file);
+ prepare_table_SESAME(&e->SESAME_water);
+ convert_units_SESAME(&e->SESAME_water, us);
+ }
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_SS08_water", 0)) {
+ char SS08_water_table_file[PARSER_MAX_LINE_SIZE];
+ set_SS08_water(&e->SS08_water, eos_planetary_id_SS08_water);
+ parser_get_param_string(params, "EoS:planetary_SS08_water_table_file",
+ SS08_water_table_file);
+ load_table_SESAME(&e->SS08_water, SS08_water_table_file);
+ prepare_table_SESAME(&e->SS08_water);
+ convert_units_SESAME(&e->SS08_water, us);
+ }
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_AQUA", 0)) {
+ char AQUA_table_file[PARSER_MAX_LINE_SIZE];
+ set_AQUA(&e->AQUA, eos_planetary_id_AQUA);
+ parser_get_param_string(params, "EoS:planetary_AQUA_table_file",
+ AQUA_table_file);
+ load_table_SESAME(&e->AQUA, AQUA_table_file);
+ prepare_table_SESAME(&e->AQUA);
+ convert_units_SESAME(&e->AQUA, us);
+ }
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_CMS19_H", 0)) {
+ char CMS19_H_table_file[PARSER_MAX_LINE_SIZE];
+ set_CMS19_H(&e->CMS19_H, eos_planetary_id_CMS19_H);
+ parser_get_param_string(params, "EoS:planetary_CMS19_H_table_file",
+ CMS19_H_table_file);
+ load_table_SESAME(&e->CMS19_H, CMS19_H_table_file);
+ prepare_table_SESAME(&e->CMS19_H);
+ convert_units_SESAME(&e->CMS19_H, us);
+ }
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_CMS19_He", 0)) {
+ char CMS19_He_table_file[PARSER_MAX_LINE_SIZE];
+ set_CMS19_He(&e->CMS19_He, eos_planetary_id_CMS19_He);
+ parser_get_param_string(params, "EoS:planetary_CMS19_He_table_file",
+ CMS19_He_table_file);
+ load_table_SESAME(&e->CMS19_He, CMS19_He_table_file);
+ prepare_table_SESAME(&e->CMS19_He);
+ convert_units_SESAME(&e->CMS19_He, us);
+ }
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_CD21_HHe", 0)) {
+ char CD21_HHe_table_file[PARSER_MAX_LINE_SIZE];
+ set_CD21_HHe(&e->CD21_HHe, eos_planetary_id_CD21_HHe);
+ parser_get_param_string(params, "EoS:planetary_CD21_HHe_table_file",
+ CD21_HHe_table_file);
+ load_table_SESAME(&e->CD21_HHe, CD21_HHe_table_file);
+ prepare_table_SESAME(&e->CD21_HHe);
+ convert_units_SESAME(&e->CD21_HHe, us);
}
// ANEOS -- using SESAME-style tables
diff --git a/src/equation_of_state/planetary/hm80.h b/src/equation_of_state/planetary/hm80.h
index d969e86e88d5807e514fad8606d4a63188e99a5d..893ef60afa9b07f37e94367824d25b8f3fddc06b 100644
--- a/src/equation_of_state/planetary/hm80.h
+++ b/src/equation_of_state/planetary/hm80.h
@@ -179,7 +179,7 @@ INLINE static void convert_units_HM80(struct HM80_params *mat,
// gas_internal_energy_from_entropy
INLINE static float HM80_internal_energy_from_entropy(
- float density, float entropy, const struct HM80_params *mat) {
+ const float density, const float entropy, const struct HM80_params *mat) {
error("This EOS function is not yet implemented!");
@@ -187,7 +187,8 @@ INLINE static float HM80_internal_energy_from_entropy(
}
// gas_pressure_from_entropy
-INLINE static float HM80_pressure_from_entropy(float density, float entropy,
+INLINE static float HM80_pressure_from_entropy(const float density,
+ const float entropy,
const struct HM80_params *mat) {
error("This EOS function is not yet implemented!");
@@ -196,7 +197,8 @@ INLINE static float HM80_pressure_from_entropy(float density, float entropy,
}
// gas_entropy_from_pressure
-INLINE static float HM80_entropy_from_pressure(float density, float pressure,
+INLINE static float HM80_entropy_from_pressure(const float density,
+ const float pressure,
const struct HM80_params *mat) {
error("This EOS function is not yet implemented!");
@@ -206,7 +208,7 @@ INLINE static float HM80_entropy_from_pressure(float density, float pressure,
// gas_soundspeed_from_entropy
INLINE static float HM80_soundspeed_from_entropy(
- float density, float entropy, const struct HM80_params *mat) {
+ const float density, const float entropy, const struct HM80_params *mat) {
error("This EOS function is not yet implemented!");
@@ -215,29 +217,25 @@ INLINE static float HM80_soundspeed_from_entropy(
// gas_entropy_from_internal_energy
INLINE static float HM80_entropy_from_internal_energy(
- float density, float u, const struct HM80_params *mat) {
+ const float density, const float u, const struct HM80_params *mat) {
return 0.f;
}
// gas_pressure_from_internal_energy
INLINE static float HM80_pressure_from_internal_energy(
- float density, float u, const struct HM80_params *mat) {
-
- float log_P, log_P_1, log_P_2, log_P_3, log_P_4;
+ const float density, const float u, const struct HM80_params *mat) {
if (u <= 0.f) {
return 0.f;
}
- int idx_rho, idx_u;
- float intp_rho, intp_u;
const float log_rho = logf(density);
const float log_u = logf(u);
// 2D interpolation (bilinear with log(rho), log(u)) to find P(rho, u)
- idx_rho = floor((log_rho - mat->log_rho_min) * mat->inv_log_rho_step);
- idx_u = floor((log_u - mat->log_u_min) * mat->inv_log_u_step);
+ int idx_rho = floor((log_rho - mat->log_rho_min) * mat->inv_log_rho_step);
+ int idx_u = floor((log_u - mat->log_u_min) * mat->inv_log_u_step);
// If outside the table then extrapolate from the edge and edge-but-one values
if (idx_rho <= -1) {
@@ -251,26 +249,31 @@ INLINE static float HM80_pressure_from_internal_energy(
idx_u = mat->num_u - 2;
}
- intp_rho = (log_rho - mat->log_rho_min - idx_rho * mat->log_rho_step) *
- mat->inv_log_rho_step;
- intp_u =
+ const float intp_rho =
+ (log_rho - mat->log_rho_min - idx_rho * mat->log_rho_step) *
+ mat->inv_log_rho_step;
+ const float intp_u =
(log_u - mat->log_u_min - idx_u * mat->log_u_step) * mat->inv_log_u_step;
// Table values
- log_P_1 = mat->table_log_P_rho_u[idx_rho * mat->num_u + idx_u];
- log_P_2 = mat->table_log_P_rho_u[idx_rho * mat->num_u + idx_u + 1];
- log_P_3 = mat->table_log_P_rho_u[(idx_rho + 1) * mat->num_u + idx_u];
- log_P_4 = mat->table_log_P_rho_u[(idx_rho + 1) * mat->num_u + idx_u + 1];
-
- log_P = (1.f - intp_rho) * ((1.f - intp_u) * log_P_1 + intp_u * log_P_2) +
- intp_rho * ((1.f - intp_u) * log_P_3 + intp_u * log_P_4);
+ const float log_P_1 = mat->table_log_P_rho_u[idx_rho * mat->num_u + idx_u];
+ const float log_P_2 =
+ mat->table_log_P_rho_u[idx_rho * mat->num_u + idx_u + 1];
+ const float log_P_3 =
+ mat->table_log_P_rho_u[(idx_rho + 1) * mat->num_u + idx_u];
+ const float log_P_4 =
+ mat->table_log_P_rho_u[(idx_rho + 1) * mat->num_u + idx_u + 1];
+
+ const float log_P =
+ (1.f - intp_rho) * ((1.f - intp_u) * log_P_1 + intp_u * log_P_2) +
+ intp_rho * ((1.f - intp_u) * log_P_3 + intp_u * log_P_4);
return expf(log_P);
}
// gas_internal_energy_from_pressure
INLINE static float HM80_internal_energy_from_pressure(
- float density, float P, const struct HM80_params *mat) {
+ const float density, const float P, const struct HM80_params *mat) {
error("This EOS function is not yet implemented!");
@@ -279,9 +282,9 @@ INLINE static float HM80_internal_energy_from_pressure(
// gas_soundspeed_from_internal_energy
INLINE static float HM80_soundspeed_from_internal_energy(
- float density, float u, const struct HM80_params *mat) {
+ const float density, const float u, const struct HM80_params *mat) {
- float c, P;
+ float c;
// Bulk modulus
if (mat->bulk_mod != 0) {
@@ -289,10 +292,10 @@ INLINE static float HM80_soundspeed_from_internal_energy(
}
// Ideal gas
else {
- P = HM80_pressure_from_internal_energy(density, u, mat);
+ const float P = HM80_pressure_from_internal_energy(density, u, mat);
c = sqrtf(hydro_gamma * P / density);
- if (c <= 0) {
+ if (c <= 0.f) {
c = sqrtf(hydro_gamma * mat->P_min_for_c_min / density);
}
}
@@ -302,7 +305,35 @@ INLINE static float HM80_soundspeed_from_internal_energy(
// gas_soundspeed_from_pressure
INLINE static float HM80_soundspeed_from_pressure(
- float density, float P, const struct HM80_params *mat) {
+ const float density, const float P, const struct HM80_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0.f;
+}
+
+// gas_entropy_from_internal_energy
+INLINE static float HM80_temperature_from_internal_energy(
+ const float density, const float u, const struct HM80_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0.f;
+}
+
+// gas_density_from_pressure_and_temperature
+INLINE static float HM80_density_from_pressure_and_temperature(
+ const float P, const float T, const struct HM80_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0.f;
+}
+
+// gas_density_from_pressure_and_internal_energy
+INLINE static float HM80_density_from_pressure_and_internal_energy(
+ const float P, const float u, const float rho_ref, const float rho_sph,
+ const struct HM80_params *mat) {
error("This EOS function is not yet implemented!");
diff --git a/src/equation_of_state/planetary/ideal_gas.h b/src/equation_of_state/planetary/ideal_gas.h
index c6da2445dad98d37cd03caa444efb9744f80d238..a2fe001e5f92791854606d2c347ee6fbac1fd88e 100644
--- a/src/equation_of_state/planetary/ideal_gas.h
+++ b/src/equation_of_state/planetary/ideal_gas.h
@@ -58,7 +58,7 @@ INLINE static void set_idg_def(struct idg_params *mat,
* @param entropy The entropy \f$A\f$.
*/
INLINE static float idg_internal_energy_from_entropy(
- float density, float entropy, const struct idg_params *mat) {
+ const float density, const float entropy, const struct idg_params *mat) {
return entropy * powf(density, mat->gamma - 1.f) *
mat->one_over_gamma_minus_one;
@@ -72,7 +72,8 @@ INLINE static float idg_internal_energy_from_entropy(
* @param density The density \f$\rho\f$.
* @param entropy The entropy \f$A\f$.
*/
-INLINE static float idg_pressure_from_entropy(float density, float entropy,
+INLINE static float idg_pressure_from_entropy(const float density,
+ const float entropy,
const struct idg_params *mat) {
return entropy * powf(density, mat->gamma);
@@ -87,7 +88,8 @@ INLINE static float idg_pressure_from_entropy(float density, float entropy,
* @param pressure The pressure \f$P\f$.
* @return The entropy \f$A\f$.
*/
-INLINE static float idg_entropy_from_pressure(float density, float pressure,
+INLINE static float idg_entropy_from_pressure(const float density,
+ const float pressure,
const struct idg_params *mat) {
return pressure * powf(density, -mat->gamma);
@@ -101,7 +103,8 @@ INLINE static float idg_entropy_from_pressure(float density, float pressure,
* @param density The density \f$\rho\f$.
* @param entropy The entropy \f$A\f$.
*/
-INLINE static float idg_soundspeed_from_entropy(float density, float entropy,
+INLINE static float idg_soundspeed_from_entropy(const float density,
+ const float entropy,
const struct idg_params *mat) {
return sqrtf(mat->gamma * powf(density, mat->gamma - 1.f) * entropy);
@@ -116,7 +119,7 @@ INLINE static float idg_soundspeed_from_entropy(float density, float entropy,
* @param u The internal energy \f$u\f$
*/
INLINE static float idg_entropy_from_internal_energy(
- float density, float u, const struct idg_params *mat) {
+ const float density, const float u, const struct idg_params *mat) {
return (mat->gamma - 1.f) * u * powf(density, 1.f - mat->gamma);
}
@@ -130,7 +133,7 @@ INLINE static float idg_entropy_from_internal_energy(
* @param u The internal energy \f$u\f$
*/
INLINE static float idg_pressure_from_internal_energy(
- float density, float u, const struct idg_params *mat) {
+ const float density, const float u, const struct idg_params *mat) {
return (mat->gamma - 1.f) * u * density;
}
@@ -145,7 +148,7 @@ INLINE static float idg_pressure_from_internal_energy(
* @return The internal energy \f$u\f$.
*/
INLINE static float idg_internal_energy_from_pressure(
- float density, float pressure, const struct idg_params *mat) {
+ const float density, const float pressure, const struct idg_params *mat) {
return mat->one_over_gamma_minus_one * pressure / density;
}
@@ -159,7 +162,7 @@ INLINE static float idg_internal_energy_from_pressure(
* @param u The internal energy \f$u\f$
*/
INLINE static float idg_soundspeed_from_internal_energy(
- float density, float u, const struct idg_params *mat) {
+ const float density, const float u, const struct idg_params *mat) {
return sqrtf(u * mat->gamma * (mat->gamma - 1.f));
}
@@ -172,10 +175,37 @@ INLINE static float idg_soundspeed_from_internal_energy(
* @param density The density \f$\rho\f$
* @param P The pressure \f$P\f$
*/
-INLINE static float idg_soundspeed_from_pressure(float density, float P,
+INLINE static float idg_soundspeed_from_pressure(const float density,
+ const float P,
const struct idg_params *mat) {
return sqrtf(mat->gamma * P / density);
}
+// gas_temperature_from_internal_energy
+INLINE static float idg_temperature_from_internal_energy(
+ const float density, const float u, const struct idg_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0.f;
+}
+
+// gas_density_from_pressure_and_temperature
+INLINE static float idg_density_from_pressure_and_temperature(
+ const float P, const float T, const struct idg_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0.f;
+}
+
+// gas_density_from_pressure_and_internal_energy
+INLINE static float idg_density_from_pressure_and_internal_energy(
+ const float P, const float u, const float rho_ref, const float rho_sph,
+ const struct idg_params *mat) {
+
+ return mat->one_over_gamma_minus_one * P / u;
+}
+
#endif /* SWIFT_IDEAL_GAS_EQUATION_OF_STATE_H */
diff --git a/src/equation_of_state/planetary/sesame.h b/src/equation_of_state/planetary/sesame.h
old mode 100755
new mode 100644
index e86111664cce4e98625629a24d2729ddf5e67f8d..cb4b8cccbcaaa4507727f1a4a3a61caaf66501ac
--- a/src/equation_of_state/planetary/sesame.h
+++ b/src/equation_of_state/planetary/sesame.h
@@ -43,6 +43,7 @@
// SESAME parameters
struct SESAME_params {
float *table_log_rho;
+ float *table_log_T;
float *table_log_u_rho_T;
float *table_P_rho_T;
float *table_c_rho_T;
@@ -77,6 +78,30 @@ INLINE static void set_SS08_water(struct SESAME_params *mat,
mat->mat_id = mat_id;
mat->version_date = 20220714;
}
+INLINE static void set_AQUA(struct SESAME_params *mat,
+ enum eos_planetary_material_id mat_id) {
+ // Haldemann et al. (2020)
+ mat->mat_id = mat_id;
+ mat->version_date = 20220714;
+}
+INLINE static void set_CMS19_H(struct SESAME_params *mat,
+ enum eos_planetary_material_id mat_id) {
+ // Chabrier et al. (2019)
+ mat->mat_id = mat_id;
+ mat->version_date = 20220905;
+}
+INLINE static void set_CMS19_He(struct SESAME_params *mat,
+ enum eos_planetary_material_id mat_id) {
+ // Chabrier et al. (2019)
+ mat->mat_id = mat_id;
+ mat->version_date = 20220905;
+}
+INLINE static void set_CD21_HHe(struct SESAME_params *mat,
+ enum eos_planetary_material_id mat_id) {
+ // Chabrier & Debras (2021)
+ mat->mat_id = mat_id;
+ mat->version_date = 20220905;
+}
INLINE static void set_ANEOS_forsterite(struct SESAME_params *mat,
enum eos_planetary_material_id mat_id) {
// Stewart et al. (2019)
@@ -176,6 +201,7 @@ INLINE static void load_table_SESAME(struct SESAME_params *mat,
// Allocate table memory
mat->table_log_rho = (float *)malloc(mat->num_rho * sizeof(float));
+ mat->table_log_T = (float *)malloc(mat->num_T * sizeof(float));
mat->table_log_u_rho_T =
(float *)malloc(mat->num_rho * mat->num_T * sizeof(float));
mat->table_P_rho_T =
@@ -197,10 +223,16 @@ INLINE static void load_table_SESAME(struct SESAME_params *mat,
}
}
- // Temperatures (ignored)
+ // Temperatures (not log yet)
for (int i_T = -1; i_T < mat->num_T; i_T++) {
- c = fscanf(f, "%f", &ignore);
- if (c != 1) error("Failed to read the SESAME EoS table %s", table_file);
+ // Ignore the first elements of rho = 0, T = 0
+ if (i_T == -1) {
+ c = fscanf(f, "%f", &ignore);
+ if (c != 1) error("Failed to read the SESAME EoS table %s", table_file);
+ } else {
+ c = fscanf(f, "%f", &mat->table_log_T[i_T]);
+ if (c != 1) error("Failed to read the SESAME EoS table %s", table_file);
+ }
}
// Sp. int. energies (not log yet), pressures, sound speeds, and sp.
@@ -232,6 +264,10 @@ INLINE static void prepare_table_SESAME(struct SESAME_params *mat) {
for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) {
mat->table_log_rho[i_rho] = logf(mat->table_log_rho[i_rho]);
}
+ // Convert temperatures to log(temperature)
+ for (int i_T = 0; i_T < mat->num_T; i_T++) {
+ mat->table_log_T[i_T] = logf(mat->table_log_T[i_T]);
+ }
// Initialise tiny values
mat->u_tiny = FLT_MAX;
@@ -297,6 +333,11 @@ INLINE static void prepare_table_SESAME(struct SESAME_params *mat) {
mat->table_log_s_rho_T[i_rho * mat->num_T + i_T] =
logf(mat->table_log_s_rho_T[i_rho * mat->num_T + i_T]);
+
+ // Ensure P > 0
+ if (mat->table_P_rho_T[i_rho * mat->num_T + i_T] <= 0) {
+ mat->table_P_rho_T[i_rho * mat->num_T + i_T] = mat->P_tiny;
+ }
}
}
}
@@ -316,6 +357,13 @@ INLINE static void convert_units_SESAME(struct SESAME_params *mat,
units_cgs_conversion_factor(us, UNIT_CONV_DENSITY));
}
+ // Temperatures (log)
+ for (int i_T = 0; i_T < mat->num_T; i_T++) {
+ mat->table_log_T[i_T] +=
+ logf(units_cgs_conversion_factor(&si, UNIT_CONV_TEMPERATURE) /
+ units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE));
+ }
+
// Sp. int. energies (log), pressures, sound speeds, and sp. entropies
for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) {
for (int i_T = 0; i_T < mat->num_T; i_T++) {
@@ -352,28 +400,26 @@ INLINE static void convert_units_SESAME(struct SESAME_params *mat,
// gas_internal_energy_from_entropy
INLINE static float SESAME_internal_energy_from_entropy(
- float density, float entropy, const struct SESAME_params *mat) {
-
- float u, log_u_1, log_u_2, log_u_3, log_u_4;
+ const float density, const float entropy, const struct SESAME_params *mat) {
+ // Return zero if entropy is zero
if (entropy <= 0.f) {
return 0.f;
}
- int idx_rho, idx_s_1, idx_s_2;
- float intp_rho, intp_s_1, intp_s_2;
const float log_rho = logf(density);
const float log_s = logf(entropy);
- // 2D interpolation (bilinear with log(rho), log(s)) to find u(rho, s))
+ // 2D interpolation (bilinear with log(rho), log(s) to find u(rho, s))
+
// Density index
- idx_rho =
+ int idx_rho =
find_value_in_monot_incr_array(log_rho, mat->table_log_rho, mat->num_rho);
// Sp. entropy at this and the next density (in relevant slice of s array)
- idx_s_1 = find_value_in_monot_incr_array(
+ int idx_s_1 = find_value_in_monot_incr_array(
log_s, mat->table_log_s_rho_T + idx_rho * mat->num_T, mat->num_T);
- idx_s_2 = find_value_in_monot_incr_array(
+ int idx_s_2 = find_value_in_monot_incr_array(
log_s, mat->table_log_s_rho_T + (idx_rho + 1) * mat->num_T, mat->num_T);
// If outside the table then extrapolate from the edge and edge-but-one values
@@ -394,6 +440,7 @@ INLINE static float SESAME_internal_energy_from_entropy(
}
// Check for duplicates in SESAME tables before interpolation
+ float intp_rho, intp_s_1, intp_s_2;
if (mat->table_log_rho[idx_rho + 1] != mat->table_log_rho[idx_rho]) {
intp_rho = (log_rho - mat->table_log_rho[idx_rho]) /
(mat->table_log_rho[idx_rho + 1] - mat->table_log_rho[idx_rho]);
@@ -420,10 +467,13 @@ INLINE static float SESAME_internal_energy_from_entropy(
}
// Table values
- log_u_1 = mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_s_1];
- log_u_2 = mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_s_1 + 1];
- log_u_3 = mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_s_2];
- log_u_4 = mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_s_2 + 1];
+ const float log_u_1 = mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_s_1];
+ const float log_u_2 =
+ mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_s_1 + 1];
+ const float log_u_3 =
+ mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_s_2];
+ const float log_u_4 =
+ mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_s_2 + 1];
// If below the minimum s at this rho then just use the lowest table values
if ((idx_rho > 0.f) && ((intp_s_1 < 0.f) || (intp_s_2 < 0.f) ||
@@ -433,18 +483,17 @@ INLINE static float SESAME_internal_energy_from_entropy(
}
// Interpolate with the log values
- u = (1.f - intp_rho) * ((1.f - intp_s_1) * log_u_1 + intp_s_1 * log_u_2) +
+ const float u =
+ (1.f - intp_rho) * ((1.f - intp_s_1) * log_u_1 + intp_s_1 * log_u_2) +
intp_rho * ((1.f - intp_s_2) * log_u_3 + intp_s_2 * log_u_4);
// Convert back from log
- u = expf(u);
-
- return u;
+ return expf(u);
}
// gas_pressure_from_entropy
INLINE static float SESAME_pressure_from_entropy(
- float density, float entropy, const struct SESAME_params *mat) {
+ const float density, const float entropy, const struct SESAME_params *mat) {
error("This EOS function is not yet implemented!");
@@ -453,7 +502,8 @@ INLINE static float SESAME_pressure_from_entropy(
// gas_entropy_from_pressure
INLINE static float SESAME_entropy_from_pressure(
- float density, float pressure, const struct SESAME_params *mat) {
+ const float density, const float pressure,
+ const struct SESAME_params *mat) {
error("This EOS function is not yet implemented!");
@@ -462,7 +512,7 @@ INLINE static float SESAME_entropy_from_pressure(
// gas_soundspeed_from_entropy
INLINE static float SESAME_soundspeed_from_entropy(
- float density, float entropy, const struct SESAME_params *mat) {
+ const float density, const float entropy, const struct SESAME_params *mat) {
error("This EOS function is not yet implemented!");
@@ -471,35 +521,32 @@ INLINE static float SESAME_soundspeed_from_entropy(
// gas_entropy_from_internal_energy
INLINE static float SESAME_entropy_from_internal_energy(
- float density, float u, const struct SESAME_params *mat) {
+ const float density, const float u, const struct SESAME_params *mat) {
return 0.f;
}
// gas_pressure_from_internal_energy
INLINE static float SESAME_pressure_from_internal_energy(
- float density, float u, const struct SESAME_params *mat) {
-
- float P, P_1, P_2, P_3, P_4;
+ const float density, const float u, const struct SESAME_params *mat) {
if (u <= 0.f) {
return 0.f;
}
- int idx_rho, idx_u_1, idx_u_2;
- float intp_rho, intp_u_1, intp_u_2;
const float log_rho = logf(density);
const float log_u = logf(u);
- // 2D interpolation (bilinear with log(rho), log(u)) to find P(rho, u))
+ // 2D interpolation (bilinear with log(rho), log(u) to find P(rho, u))
+
// Density index
- idx_rho =
+ int idx_rho =
find_value_in_monot_incr_array(log_rho, mat->table_log_rho, mat->num_rho);
// Sp. int. energy at this and the next density (in relevant slice of u array)
- idx_u_1 = find_value_in_monot_incr_array(
+ int idx_u_1 = find_value_in_monot_incr_array(
log_u, mat->table_log_u_rho_T + idx_rho * mat->num_T, mat->num_T);
- idx_u_2 = find_value_in_monot_incr_array(
+ int idx_u_2 = find_value_in_monot_incr_array(
log_u, mat->table_log_u_rho_T + (idx_rho + 1) * mat->num_T, mat->num_T);
// If outside the table then extrapolate from the edge and edge-but-one values
@@ -520,18 +567,19 @@ INLINE static float SESAME_pressure_from_internal_energy(
}
// Check for duplicates in SESAME tables before interpolation
+ float intp_rho, intp_u_1, intp_u_2;
+ const int idx_u_rho_T = idx_rho * mat->num_T + idx_u_1;
if (mat->table_log_rho[idx_rho + 1] != mat->table_log_rho[idx_rho]) {
intp_rho = (log_rho - mat->table_log_rho[idx_rho]) /
(mat->table_log_rho[idx_rho + 1] - mat->table_log_rho[idx_rho]);
} else {
intp_rho = 1.f;
}
- if (mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] !=
- mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) {
- intp_u_1 =
- (log_u - mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) /
- (mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] -
- mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]);
+ if (mat->table_log_u_rho_T[idx_u_rho_T + 1] !=
+ mat->table_log_u_rho_T[idx_u_rho_T]) {
+ intp_u_1 = (log_u - mat->table_log_u_rho_T[idx_u_rho_T]) /
+ (mat->table_log_u_rho_T[idx_u_rho_T + 1] -
+ mat->table_log_u_rho_T[idx_u_rho_T]);
} else {
intp_u_1 = 1.f;
}
@@ -546,10 +594,10 @@ INLINE static float SESAME_pressure_from_internal_energy(
}
// Table values
- P_1 = mat->table_P_rho_T[idx_rho * mat->num_T + idx_u_1];
- P_2 = mat->table_P_rho_T[idx_rho * mat->num_T + idx_u_1 + 1];
- P_3 = mat->table_P_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2];
- P_4 = mat->table_P_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2 + 1];
+ float P_1 = mat->table_P_rho_T[idx_u_rho_T];
+ float P_2 = mat->table_P_rho_T[idx_u_rho_T + 1];
+ float P_3 = mat->table_P_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2];
+ float P_4 = mat->table_P_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2 + 1];
// If below the minimum u at this rho then just use the lowest table values
if ((idx_rho > 0.f) &&
@@ -582,19 +630,18 @@ INLINE static float SESAME_pressure_from_internal_energy(
P_2 = logf(P_2);
P_3 = logf(P_3);
P_4 = logf(P_4);
+ const float P_1_2 = (1.f - intp_u_1) * P_1 + intp_u_1 * P_2;
+ const float P_3_4 = (1.f - intp_u_2) * P_3 + intp_u_2 * P_4;
- P = (1.f - intp_rho) * ((1.f - intp_u_1) * P_1 + intp_u_1 * P_2) +
- intp_rho * ((1.f - intp_u_2) * P_3 + intp_u_2 * P_4);
+ const float P = (1.f - intp_rho) * P_1_2 + intp_rho * P_3_4;
// Convert back from log
- P = expf(P);
-
- return P;
+ return expf(P);
}
// gas_internal_energy_from_pressure
INLINE static float SESAME_internal_energy_from_pressure(
- float density, float P, const struct SESAME_params *mat) {
+ const float density, const float P, const struct SESAME_params *mat) {
error("This EOS function is not yet implemented!");
@@ -603,28 +650,26 @@ INLINE static float SESAME_internal_energy_from_pressure(
// gas_soundspeed_from_internal_energy
INLINE static float SESAME_soundspeed_from_internal_energy(
- float density, float u, const struct SESAME_params *mat) {
-
- float c, c_1, c_2, c_3, c_4;
+ const float density, const float u, const struct SESAME_params *mat) {
+ // Return zero if internal energy is non-positive
if (u <= 0.f) {
return 0.f;
}
- int idx_rho, idx_u_1, idx_u_2;
- float intp_rho, intp_u_1, intp_u_2;
const float log_rho = logf(density);
const float log_u = logf(u);
- // 2D interpolation (bilinear with log(rho), log(u)) to find c(rho, u))
+ // 2D interpolation (bilinear with log(rho), log(u) to find c(rho, u))
+
// Density index
- idx_rho =
+ int idx_rho =
find_value_in_monot_incr_array(log_rho, mat->table_log_rho, mat->num_rho);
// Sp. int. energy at this and the next density (in relevant slice of u array)
- idx_u_1 = find_value_in_monot_incr_array(
+ int idx_u_1 = find_value_in_monot_incr_array(
log_u, mat->table_log_u_rho_T + idx_rho * mat->num_T, mat->num_T);
- idx_u_2 = find_value_in_monot_incr_array(
+ int idx_u_2 = find_value_in_monot_incr_array(
log_u, mat->table_log_u_rho_T + (idx_rho + 1) * mat->num_T, mat->num_T);
// If outside the table then extrapolate from the edge and edge-but-one values
@@ -645,18 +690,19 @@ INLINE static float SESAME_soundspeed_from_internal_energy(
}
// Check for duplicates in SESAME tables before interpolation
+ float intp_rho, intp_u_1, intp_u_2;
if (mat->table_log_rho[idx_rho + 1] != mat->table_log_rho[idx_rho]) {
intp_rho = (log_rho - mat->table_log_rho[idx_rho]) /
(mat->table_log_rho[idx_rho + 1] - mat->table_log_rho[idx_rho]);
} else {
intp_rho = 1.f;
}
- if (mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] !=
- mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) {
- intp_u_1 =
- (log_u - mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) /
- (mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] -
- mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]);
+ const int idx_u_rho_T = idx_rho * mat->num_T + idx_u_1;
+ if (mat->table_log_u_rho_T[idx_u_rho_T + 1] !=
+ mat->table_log_u_rho_T[idx_u_rho_T]) {
+ intp_u_1 = (log_u - mat->table_log_u_rho_T[idx_u_rho_T]) /
+ (mat->table_log_u_rho_T[idx_u_rho_T + 1] -
+ mat->table_log_u_rho_T[idx_u_rho_T]);
} else {
intp_u_1 = 1.f;
}
@@ -671,10 +717,10 @@ INLINE static float SESAME_soundspeed_from_internal_energy(
}
// Table values
- c_1 = mat->table_c_rho_T[idx_rho * mat->num_T + idx_u_1];
- c_2 = mat->table_c_rho_T[idx_rho * mat->num_T + idx_u_1 + 1];
- c_3 = mat->table_c_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2];
- c_4 = mat->table_c_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2 + 1];
+ float c_1 = mat->table_c_rho_T[idx_u_rho_T];
+ float c_2 = mat->table_c_rho_T[idx_u_rho_T + 1];
+ float c_3 = mat->table_c_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2];
+ float c_4 = mat->table_c_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2 + 1];
// If more than two table values are non-positive then return zero
int num_non_pos = 0;
@@ -710,22 +756,536 @@ INLINE static float SESAME_soundspeed_from_internal_energy(
intp_u_2 = 0;
}
- c = (1.f - intp_rho) * ((1.f - intp_u_1) * c_1 + intp_u_1 * c_2) +
- intp_rho * ((1.f - intp_u_2) * c_3 + intp_u_2 * c_4);
+ const float c = (1.f - intp_rho) * ((1.f - intp_u_1) * c_1 + intp_u_1 * c_2) +
+ intp_rho * ((1.f - intp_u_2) * c_3 + intp_u_2 * c_4);
// Convert back from log
- c = expf(c);
-
- return c;
+ return expf(c);
}
// gas_soundspeed_from_pressure
INLINE static float SESAME_soundspeed_from_pressure(
- float density, float P, const struct SESAME_params *mat) {
+ const float density, const float P, const struct SESAME_params *mat) {
error("This EOS function is not yet implemented!");
return 0.f;
}
+// gas_temperature_from_internal_energy
+INLINE static float SESAME_temperature_from_internal_energy(
+ const float density, const float u, const struct SESAME_params *mat) {
+
+ // Return zero if zero internal energy
+ if (u <= 0.f) {
+ return 0.f;
+ }
+
+ const float log_rho = logf(density);
+ const float log_u = logf(u);
+
+ // 2D interpolation (bilinear with log(rho), log(u) to find T(rho, u)
+
+ // Density index
+ int idx_rho =
+ find_value_in_monot_incr_array(log_rho, mat->table_log_rho, mat->num_rho);
+
+ // Sp. int. energy at this and the next density (in relevant slice of u array)
+ int idx_u_1 = find_value_in_monot_incr_array(
+ log_u, mat->table_log_u_rho_T + idx_rho * mat->num_T, mat->num_T);
+ int idx_u_2 = find_value_in_monot_incr_array(
+ log_u, mat->table_log_u_rho_T + (idx_rho + 1) * mat->num_T, mat->num_T);
+
+ // If outside the table then extrapolate from the edge and edge-but-one values
+ if (idx_rho <= -1) {
+ idx_rho = 0;
+ } else if (idx_rho >= mat->num_rho) {
+ idx_rho = mat->num_rho - 2;
+ }
+ if (idx_u_1 <= -1) {
+ idx_u_1 = 0;
+ } else if (idx_u_1 >= mat->num_T) {
+ idx_u_1 = mat->num_T - 2;
+ }
+ if (idx_u_2 <= -1) {
+ idx_u_2 = 0;
+ } else if (idx_u_2 >= mat->num_T) {
+ idx_u_2 = mat->num_T - 2;
+ }
+
+ // Check for duplicates in SESAME tables before interpolation
+ float intp_u_1, intp_u_2;
+ const int idx_u_rho_T = idx_rho * mat->num_T + idx_u_1;
+ if (mat->table_log_u_rho_T[idx_u_rho_T + 1] !=
+ mat->table_log_u_rho_T[idx_u_rho_T]) {
+ intp_u_1 = (log_u - mat->table_log_u_rho_T[idx_u_rho_T]) /
+ (mat->table_log_u_rho_T[idx_u_rho_T + 1] -
+ mat->table_log_u_rho_T[idx_u_rho_T]);
+ } else {
+ intp_u_1 = 1.f;
+ }
+ if (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] !=
+ mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) {
+ intp_u_2 =
+ (log_u - mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) /
+ (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] -
+ mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]);
+ } else {
+ intp_u_2 = 1.f;
+ }
+
+ // Compute line points
+ float log_rho_1 = mat->table_log_rho[idx_rho];
+ float log_rho_2 = mat->table_log_rho[idx_rho + 1];
+ float log_T_1 = mat->table_log_T[idx_u_1];
+ log_T_1 +=
+ intp_u_1 * (mat->table_log_T[idx_u_1 + 1] - mat->table_log_T[idx_u_1]);
+ float log_T_2 = mat->table_log_T[idx_u_2];
+ log_T_2 +=
+ intp_u_2 * (mat->table_log_T[idx_u_2 + 1] - mat->table_log_T[idx_u_2]);
+
+ // Intersect line passing through (log_rho_1, log_T_1), (log_rho_2, log_T_2)
+ // with line density = log_rho
+
+ // Check for log_T_1 == log_T_2
+ float log_T;
+ if (log_T_1 == log_T_2) {
+ log_T = log_T_1;
+ } else {
+ // log_rho = slope*log_T + intercept
+ const float slope = (log_rho_1 - log_rho_2) / (log_T_1 - log_T_2);
+ const float intercept = log_rho_1 - slope * log_T_1;
+ log_T = (log_rho - intercept) / slope;
+ }
+
+ // Convert back from log
+ return expf(log_T);
+}
+
+// gas_density_from_pressure_and_temperature
+INLINE static float SESAME_density_from_pressure_and_temperature(
+ float P, float T, const struct SESAME_params *mat) {
+
+ // Return zero if pressure is non-positive
+ if (P <= 0.f) {
+ return 0.f;
+ }
+
+ const float log_T = logf(T);
+
+ // 2D interpolation (bilinear with log(T), P to find rho(T, P))
+
+ // Temperature index
+ int idx_T =
+ find_value_in_monot_incr_array(log_T, mat->table_log_T, mat->num_T);
+
+ // If outside the table then extrapolate from the edge and edge-but-one values
+ if (idx_T <= -1) {
+ idx_T = 0;
+ } else if (idx_T >= mat->num_T) {
+ idx_T = mat->num_T - 2;
+ }
+
+ // Pressure at this and the next temperature (in relevant vertical slice of P
+ // array)
+ int idx_P_1 = vertical_find_value_in_monot_incr_array(
+ P, mat->table_P_rho_T, mat->num_rho, mat->num_T, idx_T);
+ int idx_P_2 = vertical_find_value_in_monot_incr_array(
+ P, mat->table_P_rho_T, mat->num_rho, mat->num_T, idx_T + 1);
+
+ // If outside the table then extrapolate from the edge and edge-but-one values
+ if (idx_P_1 <= -1) {
+ idx_P_1 = 0;
+ } else if (idx_P_1 >= mat->num_rho) {
+ idx_P_1 = mat->num_rho - 2;
+ }
+ if (idx_P_2 <= -1) {
+ idx_P_2 = 0;
+ } else if (idx_P_2 >= mat->num_rho) {
+ idx_P_2 = mat->num_rho - 2;
+ }
+
+ // Check for duplicates in SESAME tables before interpolation
+ float intp_P_1, intp_P_2;
+ if (mat->table_P_rho_T[(idx_P_1 + 1) * mat->num_T + idx_T] !=
+ mat->table_P_rho_T[idx_P_1 * mat->num_T + idx_T]) {
+ intp_P_1 = (P - mat->table_P_rho_T[idx_P_1 * mat->num_T + idx_T]) /
+ (mat->table_P_rho_T[(idx_P_1 + 1) * mat->num_T + idx_T] -
+ mat->table_P_rho_T[idx_P_1 * mat->num_T + idx_T]);
+ } else {
+ intp_P_1 = 1.f;
+ }
+ if (mat->table_P_rho_T[(idx_P_2 + 1) * mat->num_T + (idx_T + 1)] !=
+ mat->table_P_rho_T[idx_P_2 * mat->num_T + (idx_T + 1)]) {
+ intp_P_2 = (P - mat->table_P_rho_T[idx_P_2 * mat->num_T + (idx_T + 1)]) /
+ (mat->table_P_rho_T[(idx_P_2 + 1) * mat->num_T + (idx_T + 1)] -
+ mat->table_P_rho_T[idx_P_2 * mat->num_T + (idx_T + 1)]);
+ } else {
+ intp_P_2 = 1.f;
+ }
+
+ // Compute line points
+ const float log_T_1 = mat->table_log_T[idx_T];
+ const float log_T_2 = mat->table_log_T[idx_T + 1];
+ const float log_rho_1 = (1.f - intp_P_1) * mat->table_log_rho[idx_P_1] +
+ intp_P_1 * mat->table_log_rho[idx_P_1 + 1];
+ const float log_rho_2 = (1.f - intp_P_2) * mat->table_log_rho[idx_P_2] +
+ intp_P_2 * mat->table_log_rho[idx_P_2 + 1];
+
+ // Intersect line passing through (log_rho_1, log_T_1), (log_rho_2, log_T_2)
+ // with line temperature = log_T
+
+ // Check for log_rho_1 == log_rho_2
+ float log_rho;
+ if (log_rho_1 == log_rho_2) {
+ log_rho = log_rho_1;
+ } else {
+ // log_T = slope*log_rho + intercept
+ const float slope = (log_T_1 - log_T_2) / (log_rho_1 - log_rho_2);
+ const float intercept = log_T_1 - slope * log_rho_1;
+ log_rho = (log_T - intercept) / slope;
+ }
+
+ // Convert back from log
+ return expf(log_rho);
+}
+
+// gas_density_from_pressure_and_internal_energy
+INLINE static float SESAME_density_from_pressure_and_internal_energy(
+ float P, float u, float rho_ref, float rho_sph,
+ const struct SESAME_params *mat) {
+
+ // Return the unchanged density if u or P is non-positive
+ if (u <= 0.f || P <= 0.f) {
+ return rho_sph;
+ }
+
+ // Convert inputs to log
+ const float log_u = logf(u);
+ const float log_P = logf(P);
+ const float log_rho_ref = logf(rho_ref);
+
+ // Find rounded down index of reference density. This is where we start our
+ // search
+ int idx_rho_ref = find_value_in_monot_incr_array(
+ log_rho_ref, mat->table_log_rho, mat->num_rho);
+
+ // If no roots are found in the current search range, we increase search range
+ // by search_factor_log_rho above and below the reference density each
+ // iteration.
+ const float search_factor_log_rho = logf(10.f);
+
+ // Initialise the minimum and maximum densities we're searching to at the
+ // reference density. These will change before the first iteration.
+ float log_rho_min = log_rho_ref;
+ float log_rho_max = log_rho_ref;
+
+ // Initialise search indices around rho_ref
+ int idx_rho_below_min, idx_rho_above_max;
+
+ // If we find a root, it will get stored as closest_root
+ float closest_root = -1.f;
+ float root_below;
+
+ // Initialise pressures
+ float P_above_lower, P_above_upper;
+ float P_below_lower, P_below_upper;
+ P_above_upper = 0.f;
+ P_below_lower = 0.f;
+
+ // Increase search range by search_factor_log_rho
+ log_rho_max += search_factor_log_rho;
+ idx_rho_above_max = find_value_in_monot_incr_array(
+ log_rho_max, mat->table_log_rho, mat->num_rho);
+ log_rho_min -= search_factor_log_rho;
+ idx_rho_below_min = find_value_in_monot_incr_array(
+ log_rho_min, mat->table_log_rho, mat->num_rho);
+
+ // When searching above/below, we are looking for where the pressure P(rho, u)
+ // of the table densities changes from being less than to more than, or vice
+ // versa, the desired pressure. If this is the case, there is a root between
+ // these table values of rho.
+
+ // First look for roots above rho_ref
+ int idx_u_rho_T;
+ float intp_rho;
+ for (int idx_rho = idx_rho_ref; idx_rho <= idx_rho_above_max; idx_rho++) {
+
+ // This is similar to P_u_rho, but we're not interested in intp_rho,
+ // and instead calculate the pressure for both intp_rho=0 and intp_rho=1
+
+ // Sp. int. energy at this and the next density (in relevant slice of u
+ // array)
+ int idx_u_1 = find_value_in_monot_incr_array(
+ log_u, mat->table_log_u_rho_T + idx_rho * mat->num_T, mat->num_T);
+ int idx_u_2 = find_value_in_monot_incr_array(
+ log_u, mat->table_log_u_rho_T + (idx_rho + 1) * mat->num_T, mat->num_T);
+
+ // If outside the table then extrapolate from the edge and edge-but-one
+ // values
+ if (idx_rho <= -1) {
+ idx_rho = 0;
+ } else if (idx_rho >= mat->num_rho) {
+ idx_rho = mat->num_rho - 2;
+ }
+ if (idx_u_1 <= -1) {
+ idx_u_1 = 0;
+ } else if (idx_u_1 >= mat->num_T) {
+ idx_u_1 = mat->num_T - 2;
+ }
+ if (idx_u_2 <= -1) {
+ idx_u_2 = 0;
+ } else if (idx_u_2 >= mat->num_T) {
+ idx_u_2 = mat->num_T - 2;
+ }
+
+ float intp_u_1, intp_u_2;
+ idx_u_rho_T = idx_rho * mat->num_T + idx_u_1;
+ if (mat->table_log_u_rho_T[idx_u_rho_T + 1] !=
+ mat->table_log_u_rho_T[idx_u_rho_T]) {
+ intp_u_1 = (log_u - mat->table_log_u_rho_T[idx_u_rho_T]) /
+ (mat->table_log_u_rho_T[idx_u_rho_T + 1] -
+ mat->table_log_u_rho_T[idx_u_rho_T]);
+ } else {
+ intp_u_1 = 1.f;
+ }
+ if (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] !=
+ mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) {
+ intp_u_2 =
+ (log_u -
+ mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) /
+ (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] -
+ mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]);
+ } else {
+ intp_u_2 = 1.f;
+ }
+
+ // Table values
+ float P_1 = mat->table_P_rho_T[idx_u_rho_T];
+ float P_2 = mat->table_P_rho_T[idx_u_rho_T + 1];
+ float P_3 = mat->table_P_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2];
+ float P_4 = mat->table_P_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2 + 1];
+
+ // If below the minimum u at this rho then just use the lowest table values
+ if ((idx_rho > 0.f) &&
+ ((intp_u_1 < 0.f) || (intp_u_2 < 0.f) || (P_1 > P_2) || (P_3 > P_4))) {
+ intp_u_1 = 0;
+ intp_u_2 = 0;
+ }
+
+ // If more than two table values are non-positive then return zero
+ int num_non_pos = 0;
+ if (P_1 <= 0.f) num_non_pos++;
+ if (P_2 <= 0.f) num_non_pos++;
+ if (P_3 <= 0.f) num_non_pos++;
+ if (P_4 <= 0.f) num_non_pos++;
+ if (num_non_pos > 0) {
+ // If just one or two are non-positive then replace them with a tiny value
+ // Unless already trying to extrapolate in which case return zero
+ if ((num_non_pos > 2) || (mat->P_tiny == 0.f) || (intp_u_1 < 0.f) ||
+ (intp_u_2 < 0.f)) {
+ break; // return rho_sph;
+ }
+ if (P_1 <= 0.f) P_1 = mat->P_tiny;
+ if (P_2 <= 0.f) P_2 = mat->P_tiny;
+ if (P_3 <= 0.f) P_3 = mat->P_tiny;
+ if (P_4 <= 0.f) P_4 = mat->P_tiny;
+ }
+
+ // Interpolate with the log values
+ P_1 = logf(P_1);
+ P_2 = logf(P_2);
+ P_3 = logf(P_3);
+ P_4 = logf(P_4);
+ float P_1_2 = (1.f - intp_u_1) * P_1 + intp_u_1 * P_2;
+ float P_3_4 = (1.f - intp_u_2) * P_3 + intp_u_2 * P_4;
+
+ // Pressure for intp_rho = 0
+ P_above_lower = expf(P_1_2);
+
+ // Because of linear interpolation, pressures are not exactly continuous
+ // as we go from one side of a grid point to another. See if there is
+ // a root between the last P_above_upper and the new P_above_lower,
+ // which are approx the same.
+ if (idx_rho != idx_rho_ref) {
+ if ((P_above_lower - P) * (P_above_upper - P) <= 0) {
+ closest_root = expf(mat->table_log_rho[idx_rho]);
+ break;
+ }
+ }
+
+ // Pressure for intp_rho = 1
+ P_above_upper = expf(P_3_4);
+
+ // Does the pressure of the adjacent table densities switch from being
+ // above to below the desired pressure, or vice versa? If so, there is a
+ // root.
+ if ((P_above_lower - P) * (P_above_upper - P) <= 0.f) {
+
+ // If there is a root, interpolate between the table values:
+ intp_rho = (log_P - P_1_2) / (P_3_4 - P_1_2);
+
+ closest_root = expf(mat->table_log_rho[idx_rho] +
+ intp_rho * (mat->table_log_rho[idx_rho + 1] -
+ mat->table_log_rho[idx_rho]));
+
+ // If the root is between the same table values as the reference value,
+ // then this is the closest root, so we can return it without further
+ // searching
+ if (idx_rho == idx_rho_ref) {
+ return closest_root;
+ }
+
+ // Found a root, so no need to search higher densities
+ break;
+ }
+ }
+
+ // If we found a root above, change search range below so that we're only
+ // looking for closer (in log) roots than the one we found
+ if (closest_root > 0.f) {
+ log_rho_min = log_rho_ref - (logf(closest_root) - log_rho_ref);
+ idx_rho_below_min = find_value_in_monot_incr_array(
+ log_rho_min, mat->table_log_rho, mat->num_rho);
+ }
+
+ // Now look for roots below rho_ref
+ for (int idx_rho = idx_rho_ref; idx_rho >= idx_rho_below_min; idx_rho--) {
+
+ // Sp. int. energy at this and the next density (in relevant slice of u
+ // array)
+ int idx_u_1 = find_value_in_monot_incr_array(
+ log_u, mat->table_log_u_rho_T + idx_rho * mat->num_T, mat->num_T);
+ int idx_u_2 = find_value_in_monot_incr_array(
+ log_u, mat->table_log_u_rho_T + (idx_rho + 1) * mat->num_T, mat->num_T);
+
+ // If outside the table then extrapolate from the edge and edge-but-one
+ // values
+ if (idx_rho <= -1) {
+ idx_rho = 0;
+ } else if (idx_rho >= mat->num_rho) {
+ idx_rho = mat->num_rho - 2;
+ }
+ if (idx_u_1 <= -1) {
+ idx_u_1 = 0;
+ } else if (idx_u_1 >= mat->num_T) {
+ idx_u_1 = mat->num_T - 2;
+ }
+ if (idx_u_2 <= -1) {
+ idx_u_2 = 0;
+ } else if (idx_u_2 >= mat->num_T) {
+ idx_u_2 = mat->num_T - 2;
+ }
+
+ idx_u_rho_T = idx_rho * mat->num_T + idx_u_1;
+ float intp_u_1, intp_u_2;
+ if (mat->table_log_u_rho_T[idx_u_rho_T + 1] !=
+ mat->table_log_u_rho_T[idx_u_rho_T]) {
+ intp_u_1 = (log_u - mat->table_log_u_rho_T[idx_u_rho_T]) /
+ (mat->table_log_u_rho_T[idx_u_rho_T + 1] -
+ mat->table_log_u_rho_T[idx_u_rho_T]);
+ } else {
+ intp_u_1 = 1.f;
+ }
+ if (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] !=
+ mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) {
+ intp_u_2 =
+ (log_u -
+ mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) /
+ (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] -
+ mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]);
+ } else {
+ intp_u_2 = 1.f;
+ }
+
+ // Table values
+ float P_1 = mat->table_P_rho_T[idx_u_rho_T];
+ float P_2 = mat->table_P_rho_T[idx_u_rho_T + 1];
+ float P_3 = mat->table_P_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2];
+ float P_4 = mat->table_P_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2 + 1];
+
+ // If below the minimum u at this rho then just use the lowest table values
+ if ((idx_rho > 0.f) &&
+ ((intp_u_1 < 0.f) || (intp_u_2 < 0.f) || (P_1 > P_2) || (P_3 > P_4))) {
+ intp_u_1 = 0;
+ intp_u_2 = 0;
+ }
+
+ // If more than two table values are non-positive then return zero
+ int num_non_pos = 0;
+ if (P_1 <= 0.f) num_non_pos++;
+ if (P_2 <= 0.f) num_non_pos++;
+ if (P_3 <= 0.f) num_non_pos++;
+ if (P_4 <= 0.f) num_non_pos++;
+ if (num_non_pos > 0) {
+ // If just one or two are non-positive then replace them with a tiny value
+ // Unless already trying to extrapolate in which case return zero
+ if ((num_non_pos > 2) || (mat->P_tiny == 0.f) || (intp_u_1 < 0.f) ||
+ (intp_u_2 < 0.f)) {
+ break;
+ }
+ if (P_1 <= 0.f) P_1 = mat->P_tiny;
+ if (P_2 <= 0.f) P_2 = mat->P_tiny;
+ if (P_3 <= 0.f) P_3 = mat->P_tiny;
+ if (P_4 <= 0.f) P_4 = mat->P_tiny;
+ }
+
+ // Interpolate with the log values
+ P_1 = logf(P_1);
+ P_2 = logf(P_2);
+ P_3 = logf(P_3);
+ P_4 = logf(P_4);
+ const float P_1_2 = (1.f - intp_u_1) * P_1 + intp_u_1 * P_2;
+ const float P_3_4 = (1.f - intp_u_2) * P_3 + intp_u_2 * P_4;
+
+ // Pressure for intp_rho = 1
+ P_below_upper = expf(P_3_4);
+ // Because of linear interpolation, pressures are not exactly continuous
+ // as we go from one side of a grid point to another. See if there is
+ // a root between the last P_below_lower and the new P_below_upper,
+ // which are approx the same.
+ if (idx_rho != idx_rho_ref) {
+ if ((P_below_lower - P) * (P_below_upper - P) <= 0) {
+ closest_root = expf(mat->table_log_rho[idx_rho + 1]);
+ break;
+ }
+ }
+ // Pressure for intp_rho = 0
+ P_below_lower = expf(P_1_2);
+
+ // Does the pressure of the adjacent table densities switch from being
+ // above to below the desired pressure, or vice versa? If so, there is a
+ // root.
+ if ((P_below_lower - P) * (P_below_upper - P) <= 0.f) {
+
+ // If there is a root, interpolate between the table values:
+ intp_rho = (log_P - P_1_2) / (P_3_4 - P_1_2);
+
+ root_below = expf(mat->table_log_rho[idx_rho] +
+ intp_rho * (mat->table_log_rho[idx_rho + 1] -
+ mat->table_log_rho[idx_rho]));
+
+ // If we found a root above, which one is closer to the reference rho?
+ if (closest_root > 0.f) {
+ if (fabsf(logf(root_below) - logf(rho_ref)) <
+ fabsf(logf(closest_root) - logf(rho_ref))) {
+ closest_root = root_below;
+ }
+ } else {
+ closest_root = root_below;
+ }
+ // Found a root, so no need to search higher densities
+ break;
+ }
+ }
+
+ // Return the root if we found one
+ if (closest_root > 0.f) {
+ return closest_root;
+ }
+
+ // If we don't find a root before we reach max_counter, return rho_ref. Maybe
+ // we should give an error here?
+ return rho_sph;
+}
#endif /* SWIFT_SESAME_EQUATION_OF_STATE_H */
diff --git a/src/equation_of_state/planetary/tillotson.h b/src/equation_of_state/planetary/tillotson.h
index 5831b1785485d5dc2563718c94fc9cc04a08b822..dc5b5576c33616446795451fd0229e6e58fdfec4 100644
--- a/src/equation_of_state/planetary/tillotson.h
+++ b/src/equation_of_state/planetary/tillotson.h
@@ -37,12 +37,16 @@
#include "equation_of_state.h"
#include "inline.h"
#include "physical_constants.h"
+#include "sesame.h"
#include "units.h"
// Tillotson parameters
struct Til_params {
- float rho_0, a, b, A, B, u_0, u_iv, u_cv, alpha, beta, eta_min, eta_zero,
- P_min;
+ float rho_0, a, b, A, B, u_0, u_iv, u_cv, alpha, beta;
+ float eta_min, eta_zero, P_min, C_V;
+ float *A1_u_cold;
+ float rho_cold_min, rho_cold_max;
+ float rho_min, rho_max;
enum eos_planetary_material_id mat_id;
};
@@ -62,7 +66,12 @@ INLINE static void set_Til_iron(struct Til_params *mat,
mat->beta = 5.0f;
mat->eta_min = 0.0f;
mat->eta_zero = 0.0f;
- mat->P_min = 0.0f;
+ mat->P_min = 0.01f;
+ mat->C_V = 449.0f;
+ mat->rho_cold_min = 100.0f;
+ mat->rho_cold_max = 1.0e5f;
+ mat->rho_min = 1.0f;
+ mat->rho_max = 1.0e5f;
}
INLINE static void set_Til_granite(struct Til_params *mat,
enum eos_planetary_material_id mat_id) {
@@ -79,7 +88,12 @@ INLINE static void set_Til_granite(struct Til_params *mat,
mat->beta = 5.0f;
mat->eta_min = 0.0f;
mat->eta_zero = 0.0f;
- mat->P_min = 0.0f;
+ mat->P_min = 0.01f;
+ mat->C_V = 790.0f;
+ mat->rho_cold_min = 100.0f;
+ mat->rho_cold_max = 1.0e5f;
+ mat->rho_min = 1.0f;
+ mat->rho_max = 1.0e5f;
}
INLINE static void set_Til_basalt(struct Til_params *mat,
enum eos_planetary_material_id mat_id) {
@@ -96,7 +110,12 @@ INLINE static void set_Til_basalt(struct Til_params *mat,
mat->beta = 5.0f;
mat->eta_min = 0.0f;
mat->eta_zero = 0.0f;
- mat->P_min = 0.0f;
+ mat->P_min = 0.01f;
+ mat->C_V = 790.0f;
+ mat->rho_cold_min = 100.0f;
+ mat->rho_cold_max = 1.0e5f;
+ mat->rho_min = 1.0f;
+ mat->rho_max = 1.0e5f;
}
INLINE static void set_Til_water(struct Til_params *mat,
enum eos_planetary_material_id mat_id) {
@@ -113,7 +132,77 @@ INLINE static void set_Til_water(struct Til_params *mat,
mat->beta = 5.0f;
mat->eta_min = 0.925f;
mat->eta_zero = 0.875f;
+ mat->P_min = 0.01f;
+ mat->C_V = 4186.0f;
+ mat->rho_cold_min = 100.0f;
+ mat->rho_cold_max = 1.0e5f;
+ mat->rho_min = 1.0f;
+ mat->rho_max = 1.0e5f;
+}
+INLINE static void set_Til_ice(struct Til_params *mat,
+ enum eos_planetary_material_id mat_id) {
+ mat->mat_id = mat_id;
+ mat->rho_0 = 1293.0f;
+ mat->a = 0.3f;
+ mat->b = 0.1f;
+ mat->A = 1.07e10f;
+ mat->B = 6.5e10f;
+ mat->u_0 = 1.0e7f;
+ mat->u_iv = 7.73e5f;
+ mat->u_cv = 3.04e6f;
+ mat->alpha = 10.0f;
+ mat->beta = 5.0f;
+ mat->eta_min = 0.925f;
+ mat->eta_zero = 0.875f;
mat->P_min = 0.0f;
+ mat->C_V = 2093.0f;
+ mat->rho_cold_min = 100.0f;
+ mat->rho_cold_max = 1.0e5f;
+ mat->rho_min = 1.0f;
+ mat->rho_max = 1.0e5f;
+}
+
+/*
+ Read the parameters from a file.
+
+ File contents
+ -------------
+ # header (5 lines)
+ rho_0 (kg/m3) a (-) b (-) A (Pa) B (Pa)
+ u_0 (J) u_iv (J) u_cv (J) alpha (-) beta (-)
+ eta_min (-) eta_zero (-) P_min (Pa) C_V (J kg^-1 K^-1)
+ rho_cold_min (kg/m3) rho_cold_max (kg/m3) rho_min (kg/m3) rho_max (kg/m3)
+*/
+INLINE static void set_Til_custom(struct Til_params *mat,
+ enum eos_planetary_material_id mat_id,
+ char *param_file) {
+ mat->mat_id = mat_id;
+
+ // Load table contents from file
+ FILE *f = fopen(param_file, "r");
+ if (f == NULL)
+ error("Failed to open the Tillotson EoS file '%s'", param_file);
+
+ // Skip header lines
+ skip_lines(f, 5);
+
+ // Read parameters (SI)
+ int c;
+ c = fscanf(f, "%f %f %f %f %f", &mat->rho_0, &mat->a, &mat->b, &mat->A,
+ &mat->B);
+ if (c != 5) error("Failed to read the Tillotson EoS file %s", param_file);
+
+ c = fscanf(f, "%f %f %f %f %f", &mat->u_0, &mat->u_iv, &mat->u_cv,
+ &mat->alpha, &mat->beta);
+ if (c != 5) error("Failed to read the Tillotson EoS file %s", param_file);
+
+ c = fscanf(f, "%f %f %f %f", &mat->eta_min, &mat->eta_zero, &mat->P_min,
+ &mat->C_V);
+ if (c != 4) error("Failed to read the Tillotson EoS file %s", param_file);
+
+ c = fscanf(f, "%f %f %f %f", &mat->rho_cold_min, &mat->rho_cold_max,
+ &mat->rho_min, &mat->rho_max);
+ if (c != 4) error("Failed to read the Tillotson EoS file %s", param_file);
}
// Convert to internal units
@@ -123,6 +212,8 @@ INLINE static void convert_units_Til(struct Til_params *mat,
struct unit_system si;
units_init_si(&si);
+ const int N = 10000;
+
// SI to cgs
mat->rho_0 *= units_cgs_conversion_factor(&si, UNIT_CONV_DENSITY);
mat->A *= units_cgs_conversion_factor(&si, UNIT_CONV_PRESSURE);
@@ -132,6 +223,19 @@ INLINE static void convert_units_Til(struct Til_params *mat,
mat->u_cv *= units_cgs_conversion_factor(&si, UNIT_CONV_ENERGY_PER_UNIT_MASS);
mat->P_min *= units_cgs_conversion_factor(&si, UNIT_CONV_PRESSURE);
+ for (int i = 0; i < N; i++) {
+ mat->A1_u_cold[i] *=
+ units_cgs_conversion_factor(&si, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+ }
+
+ mat->C_V *= units_cgs_conversion_factor(&si, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+ // Entropy units don't work? using internal kelvin
+
+ mat->rho_cold_min *= units_cgs_conversion_factor(&si, UNIT_CONV_DENSITY);
+ mat->rho_cold_max *= units_cgs_conversion_factor(&si, UNIT_CONV_DENSITY);
+ mat->rho_min *= units_cgs_conversion_factor(&si, UNIT_CONV_DENSITY);
+ mat->rho_max *= units_cgs_conversion_factor(&si, UNIT_CONV_DENSITY);
+
// cgs to internal
mat->rho_0 /= units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
mat->A /= units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
@@ -140,11 +244,24 @@ INLINE static void convert_units_Til(struct Til_params *mat,
mat->u_iv /= units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
mat->u_cv /= units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
mat->P_min /= units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
+
+ for (int i = 0; i < N; i++) {
+ mat->A1_u_cold[i] /=
+ units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+ }
+
+ mat->C_V /= units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+ // Entropy units don't work? using internal kelvin
+
+ mat->rho_cold_min /= units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
+ mat->rho_cold_max /= units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
+ mat->rho_min /= units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
+ mat->rho_max /= units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
}
// gas_internal_energy_from_entropy
INLINE static float Til_internal_energy_from_entropy(
- float density, float entropy, const struct Til_params *mat) {
+ const float density, const float entropy, const struct Til_params *mat) {
error("This EOS function is not yet implemented!");
@@ -152,7 +269,8 @@ INLINE static float Til_internal_energy_from_entropy(
}
// gas_pressure_from_entropy
-INLINE static float Til_pressure_from_entropy(float density, float entropy,
+INLINE static float Til_pressure_from_entropy(const float density,
+ const float entropy,
const struct Til_params *mat) {
error("This EOS function is not yet implemented!");
@@ -161,7 +279,8 @@ INLINE static float Til_pressure_from_entropy(float density, float entropy,
}
// gas_entropy_from_pressure
-INLINE static float Til_entropy_from_pressure(float density, float pressure,
+INLINE static float Til_entropy_from_pressure(const float density,
+ const float pressure,
const struct Til_params *mat) {
error("This EOS function is not yet implemented!");
@@ -170,7 +289,8 @@ INLINE static float Til_entropy_from_pressure(float density, float pressure,
}
// gas_soundspeed_from_entropy
-INLINE static float Til_soundspeed_from_entropy(float density, float entropy,
+INLINE static float Til_soundspeed_from_entropy(const float density,
+ const float entropy,
const struct Til_params *mat) {
error("This EOS function is not yet implemented!");
@@ -180,14 +300,14 @@ INLINE static float Til_soundspeed_from_entropy(float density, float entropy,
// gas_entropy_from_internal_energy
INLINE static float Til_entropy_from_internal_energy(
- float density, float u, const struct Til_params *mat) {
+ const float density, const float u, const struct Til_params *mat) {
return 0.f;
}
// gas_pressure_from_internal_energy
INLINE static float Til_pressure_from_internal_energy(
- float density, float u, const struct Til_params *mat) {
+ const float density, const float u, const struct Til_params *mat) {
const float eta = density / mat->rho_0;
const float eta_sq = eta * eta;
@@ -195,10 +315,9 @@ INLINE static float Til_pressure_from_internal_energy(
const float nu = 1.f / eta - 1.f;
const float w = u / (mat->u_0 * eta_sq) + 1.f;
const float w_inv = 1.f / w;
- float P_c, P_e, P;
// Condensed or cold
- P_c =
+ float P_c =
(mat->a + mat->b * w_inv) * density * u + mat->A * mu + mat->B * mu * mu;
if (eta < mat->eta_zero) {
P_c = 0.f;
@@ -206,11 +325,12 @@ INLINE static float Til_pressure_from_internal_energy(
P_c *= (eta - mat->eta_zero) / (mat->eta_min - mat->eta_zero);
}
// Expanded and hot
- P_e = mat->a * density * u +
- (mat->b * density * u * w_inv + mat->A * mu * expf(-mat->beta * nu)) *
- expf(-mat->alpha * nu * nu);
+ float P_e = mat->a * density * u + (mat->b * density * u * w_inv +
+ mat->A * mu * expf(-mat->beta * nu)) *
+ expf(-mat->alpha * nu * nu);
// Condensed or cold state
+ float P;
if ((1.f < eta) || (u < mat->u_iv)) {
P = P_c;
}
@@ -234,7 +354,7 @@ INLINE static float Til_pressure_from_internal_energy(
// gas_internal_energy_from_pressure
INLINE static float Til_internal_energy_from_pressure(
- float density, float P, const struct Til_params *mat) {
+ const float density, const float P, const struct Til_params *mat) {
error("This EOS function is not yet implemented!");
@@ -243,7 +363,7 @@ INLINE static float Til_internal_energy_from_pressure(
// gas_soundspeed_from_internal_energy
INLINE static float Til_soundspeed_from_internal_energy(
- float density, float u, const struct Til_params *mat) {
+ const float density, const float u, const struct Til_params *mat) {
const float rho_0_inv = 1.f / mat->rho_0;
const float eta = density * rho_0_inv;
@@ -256,25 +376,25 @@ INLINE static float Til_soundspeed_from_internal_energy(
const float w_inv_sq = w_inv * w_inv;
const float exp_beta = expf(-mat->beta * nu);
const float exp_alpha = expf(-mat->alpha * nu * nu);
- float P_c, P_e, c_sq_c, c_sq_e, c_sq;
// Condensed or cold
- P_c =
+ float P_c =
(mat->a + mat->b * w_inv) * density * u + mat->A * mu + mat->B * mu * mu;
if (eta < mat->eta_zero) {
P_c = 0.f;
} else if (eta < mat->eta_min) {
P_c *= (eta - mat->eta_zero) / (mat->eta_min - mat->eta_zero);
}
- c_sq_c = P_c * rho_inv * (1.f + mat->a + mat->b * w_inv) +
- mat->b * (w - 1.f) * w_inv_sq * (2.f * u - P_c * rho_inv) +
- rho_inv * (mat->A + mat->B * (eta_sq - 1.f));
+ float c_sq_c = P_c * rho_inv * (1.f + mat->a + mat->b * w_inv) +
+ mat->b * (w - 1.f) * w_inv_sq * (2.f * u - P_c * rho_inv) +
+ rho_inv * (mat->A + mat->B * (eta_sq - 1.f));
// Expanded and hot
- P_e = mat->a * density * u +
- (mat->b * density * u * w_inv + mat->A * mu * exp_beta) * exp_alpha;
+ float P_e =
+ mat->a * density * u +
+ (mat->b * density * u * w_inv + mat->A * mu * exp_beta) * exp_alpha;
- c_sq_e =
+ float c_sq_e =
P_e * rho_inv * (1.f + mat->a + mat->b * w_inv * exp_alpha) +
(mat->b * density * u * w_inv_sq / eta_sq *
(rho_inv / mat->u_0 * (2.f * u - P_e * rho_inv) +
@@ -285,6 +405,7 @@ INLINE static float Til_soundspeed_from_internal_energy(
exp_alpha;
// Condensed or cold state
+ float c_sq;
if ((1.f < eta) || (u < mat->u_iv)) {
c_sq = c_sq_c;
}
@@ -304,7 +425,8 @@ INLINE static float Til_soundspeed_from_internal_energy(
}
// gas_soundspeed_from_pressure
-INLINE static float Til_soundspeed_from_pressure(float density, float P,
+INLINE static float Til_soundspeed_from_pressure(const float density,
+ const float P,
const struct Til_params *mat) {
error("This EOS function is not yet implemented!");
@@ -312,4 +434,494 @@ INLINE static float Til_soundspeed_from_pressure(float density, float P,
return 0.f;
}
+// Compute u cold
+INLINE static float compute_u_cold(
+ const float density, const struct Til_params *mat,
+ const enum eos_planetary_material_id mat_id) {
+ const int N = 10000;
+ const float rho_0 = mat->rho_0;
+ const float drho = (density - rho_0) / N;
+
+ float x = rho_0;
+ float u_cold = 1e-9;
+ for (int i = 0; i < N; i++) {
+ x += drho;
+ u_cold +=
+ Til_pressure_from_internal_energy(x, u_cold, mat) * drho / (x * x);
+ }
+
+ return u_cold;
+}
+
+// Compute A1_u_cold
+INLINE static void set_Til_u_cold(struct Til_params *mat,
+ enum eos_planetary_material_id mat_id) {
+
+ const int N = 10000;
+ const float rho_min = 100.f;
+ const float rho_max = 100000.f;
+
+ // Allocate table memory
+ mat->A1_u_cold = (float *)malloc(N * sizeof(float));
+
+ float rho = rho_min;
+ const float drho = (rho_max - rho_min) / (N - 1);
+
+ for (int i = 0; i < N; i++) {
+ mat->A1_u_cold[i] = compute_u_cold(rho, mat, mat_id);
+ rho += drho;
+ }
+}
+
+// Compute u cold fast from precomputed values
+INLINE static float compute_fast_u_cold(const float density,
+ const struct Til_params *mat) {
+
+ const int N = 10000;
+ const float rho_min = mat->rho_cold_min;
+ const float rho_max = mat->rho_cold_max;
+
+ const float drho = (rho_max - rho_min) / (N - 1);
+
+ const int a = (int)((density - rho_min) / drho);
+ const int b = a + 1;
+
+ float u_cold;
+ if (a >= 0 && a < (N - 1)) {
+ u_cold = mat->A1_u_cold[a];
+ u_cold += ((mat->A1_u_cold[b] - mat->A1_u_cold[a]) / drho) *
+ (density - rho_min - a * drho);
+ } else if (density < rho_min) {
+ u_cold = mat->A1_u_cold[0];
+ } else {
+ u_cold = mat->A1_u_cold[N - 1];
+ u_cold += ((mat->A1_u_cold[N - 1] - mat->A1_u_cold[N - 2]) / drho) *
+ (density - rho_max);
+ }
+ return u_cold;
+}
+
+// gas_temperature_from_internal_energy
+INLINE static float Til_temperature_from_internal_energy(
+ const float density, const float u, const struct Til_params *mat) {
+
+ const float u_cold = compute_fast_u_cold(density, mat);
+
+ float T = (u - u_cold) / (mat->C_V);
+ if (T < 0.f) {
+ T = 0.f;
+ }
+
+ return T;
+}
+
+// gas_pressure_from_density_and_temperature
+INLINE static float Til_pressure_from_temperature(
+ const float density, const float T, const struct Til_params *mat) {
+
+ const float u = compute_fast_u_cold(density, mat) + mat->C_V * T;
+ const float P = Til_pressure_from_internal_energy(density, u, mat);
+
+ return P;
+}
+
+// gas_density_from_pressure_and_temperature
+INLINE static float Til_density_from_pressure_and_temperature(
+ const float P, const float T, const struct Til_params *mat) {
+
+ float rho_min = mat->rho_min;
+ float rho_max = mat->rho_max;
+ float rho_mid = (rho_min + rho_max) / 2.f;
+
+ // Check for P == 0 or T == 0
+ float P_des;
+ if (P <= mat->P_min) {
+ P_des = mat->P_min;
+ } else {
+ P_des = P;
+ }
+
+ float P_min = Til_pressure_from_temperature(rho_min, T, mat);
+ float P_mid = Til_pressure_from_temperature(rho_mid, T, mat);
+ float P_max = Til_pressure_from_temperature(rho_max, T, mat);
+
+ // quick fix?
+ if (P_des < P_min) {
+ P_des = P_min;
+ }
+
+ const float tolerance = 0.001 * rho_min;
+ int counter = 0;
+ const int max_counter = 200;
+ if (P_des >= P_min && P_des <= P_max) {
+ while ((rho_max - rho_min) > tolerance && counter < max_counter) {
+
+ P_min = Til_pressure_from_temperature(rho_min, T, mat);
+ P_mid = Til_pressure_from_temperature(rho_mid, T, mat);
+ P_max = Til_pressure_from_temperature(rho_max, T, mat);
+
+ const float f0 = P_des - P_min;
+ const float f2 = P_des - P_mid;
+
+ if ((f0 * f2) > 0) {
+ rho_min = rho_mid;
+ } else {
+ rho_max = rho_mid;
+ }
+
+ rho_mid = (rho_min + rho_max) / 2.f;
+ counter += 1;
+ }
+ } else {
+ error("Error in Til_density_from_pressure_and_temperature");
+ return 0.f;
+ }
+ return rho_mid;
+}
+
+// gas_density_from_pressure_and_internal_energy
+INLINE static float Til_density_from_pressure_and_internal_energy(
+ const float P, const float u, const float rho_ref, const float rho_sph,
+ const struct Til_params *mat) {
+
+ if (P <= mat->P_min || u == 0) {
+ return rho_sph;
+ }
+
+ // We start search on the same curve as rho_ref, since this is most likely
+ // curve to find rho on
+ float eta_ref = rho_ref / mat->rho_0;
+
+ /*
+ There are 5 possible curves:
+ 1: cold_min
+ 2: cold
+ 3: hybrid_min
+ 4: hybrid
+ 5: hot
+
+ These curves cover different eta ranges within three different regions of u:
+ u REGION 1 (u < u_iv):
+ eta < eta_min: cold_min
+ eta_min < eta: cold
+
+ u REGION 2 (u_iv < u < u_cv):
+ eta < eta_min: hybrid_min
+ eta_min < eta < 1: hybrid
+ 1 < eta: cold
+
+ u REGION 3 (u_cv < u):
+ eta < 1: hot
+ 1 < eta: cold
+
+ NOTE: for a lot of EoS, eta_min = 0, so search this region last if given the
+ option to save time for most EoS
+ */
+ enum Til_region {
+ Til_region_none,
+ Til_region_cold_min,
+ Til_region_cold,
+ Til_region_hybrid_min,
+ Til_region_hybrid,
+ Til_region_hot
+ };
+
+ // Based on our u region, what possible curves can rho be on? Ordered based on
+ // order we search for roots. Numbers correspond to curves in order given
+ // above. 0 is a dummy which breaks the loop.
+ enum Til_region possible_curves[3];
+ // u REGION 1
+ if (u <= mat->u_iv) {
+ if (eta_ref <= mat->eta_min) {
+ possible_curves[0] = Til_region_cold_min;
+ possible_curves[1] = Til_region_cold;
+ possible_curves[2] = Til_region_none;
+ } else {
+ possible_curves[0] = Til_region_cold;
+ possible_curves[1] = Til_region_cold_min;
+ possible_curves[2] = Til_region_none;
+ }
+ // u REGION 2
+ } else if (u <= mat->u_cv) {
+ if (eta_ref <= mat->eta_min) {
+ possible_curves[0] = Til_region_hybrid_min;
+ possible_curves[1] = Til_region_hybrid;
+ possible_curves[2] = Til_region_cold;
+ } else if (eta_ref <= 1) {
+ possible_curves[0] = Til_region_hybrid;
+ possible_curves[1] = Til_region_cold;
+ possible_curves[2] = Til_region_hybrid_min;
+ } else {
+ possible_curves[0] = Til_region_cold;
+ possible_curves[1] = Til_region_hybrid;
+ possible_curves[2] = Til_region_hybrid_min;
+ }
+ // u REGION 3
+ } else {
+ if (eta_ref <= 1) {
+ possible_curves[0] = Til_region_hot;
+ possible_curves[1] = Til_region_cold;
+ possible_curves[2] = Til_region_none;
+ } else {
+ possible_curves[0] = Til_region_cold;
+ possible_curves[1] = Til_region_hot;
+ possible_curves[2] = Til_region_none;
+ }
+ }
+ // Newton-Raphson
+ const int max_iter = 10;
+ const float tol = 1e-5;
+
+ // loops over possible curves
+ for (int i = 0; i < 3; i++) {
+
+ enum Til_region curve = possible_curves[i];
+
+ // if there are only two possible curves, break when we get to three and
+ // haven't found a root
+ if (curve == Til_region_none) {
+ break;
+ }
+
+ // Start iteration at reference value
+ float rho_iter = rho_ref;
+
+ // Constrain our initial guess to be on the curve we're currently looking at
+ // in the first loop, this is already satisfied.
+ if (i > 0) {
+
+ // u REGION 1
+ if (u <= mat->u_iv) {
+ if (curve == Til_region_cold_min) {
+ if (rho_iter > mat->eta_min * mat->rho_0) {
+ rho_iter = mat->eta_min * mat->rho_0;
+ }
+ } else if (curve == Til_region_cold) {
+ if (rho_iter < mat->eta_min * mat->rho_0) {
+ rho_iter = mat->eta_min * mat->rho_0;
+ }
+ } else {
+ error("Error in Til_density_from_pressure_and_internal_energy");
+ }
+ // u REGION 2
+ } else if (u <= mat->u_cv) {
+ if (curve == Til_region_hybrid_min) {
+ if (rho_iter > mat->eta_min * mat->rho_0) {
+ rho_iter = mat->eta_min * mat->rho_0;
+ }
+ } else if (curve == Til_region_hybrid) {
+ if (rho_iter < mat->eta_min * mat->rho_0) {
+ rho_iter = mat->eta_min * mat->rho_0;
+ }
+ if (rho_iter > mat->rho_0) {
+ rho_iter = mat->rho_0;
+ }
+ } else if (curve == Til_region_cold) {
+ if (rho_iter < mat->rho_0) {
+ rho_iter = mat->rho_0;
+ }
+ } else {
+ error("Error in Til_density_from_pressure_and_internal_energy");
+ }
+
+ // u REGION 3
+ } else {
+ if (curve == Til_region_hot) {
+ if (rho_iter > mat->rho_0) {
+ rho_iter = mat->rho_0;
+ }
+ } else if (curve == Til_region_cold) {
+ if (rho_iter < mat->rho_0) {
+ rho_iter = mat->rho_0;
+ }
+ } else {
+ error("Error in Til_density_from_pressure_and_internal_energy");
+ }
+ }
+ }
+
+ // Set this to an arbitrary number so we definitely dont't think we converge
+ // straigt away
+ float last_rho_iter = -1e5;
+
+ // Now iterate
+ for (int j = 0; j < max_iter; j++) {
+
+ const float eta_iter = rho_iter / mat->rho_0;
+ const float eta_iter_sq = eta_iter * eta_iter;
+ const float mu_iter = eta_iter - 1.0f;
+ const float nu_iter = 1.0f / eta_iter - 1.0f;
+ const float w_iter = u / (mat->u_0 * eta_iter_sq) + 1.0f;
+ const float w_iter_inv = 1.0f / w_iter;
+ const float exp1 = expf(-mat->beta * nu_iter);
+ const float exp2 = expf(-mat->alpha * nu_iter * nu_iter);
+
+ // Derivatives
+ const float dw_inv_drho_iter =
+ (2.f * mat->u_0 * u * eta_iter / mat->rho_0) /
+ ((u + mat->u_0 * eta_iter_sq) * (u + mat->u_0 * eta_iter_sq));
+ const float dmu_drho_iter = 1.f / mat->rho_0;
+ const float dmu_sq_drho_iter =
+ 2.f * rho_iter / (mat->rho_0 * mat->rho_0) - 2.f / mat->rho_0;
+ const float dexp1_drho_iter =
+ mat->beta * mat->rho_0 * exp1 / (rho_iter * rho_iter);
+ const float dexp2_drho_iter = 2.f * mat->alpha * mat->rho_0 *
+ (mat->rho_0 - rho_iter) * exp2 /
+ (rho_iter * rho_iter * rho_iter);
+
+ // Use P_fraction to determine whether we've converged on a root
+ float P_fraction;
+
+ // Newton-Raphson
+ float P_c_iter, dP_c_drho_iter, P_h_iter, dP_h_drho_iter;
+
+ // if "cold" or "hybrid"
+ if (curve == Til_region_cold || curve == Til_region_hybrid) {
+ P_c_iter = (mat->a + mat->b * w_iter_inv) * rho_iter * u +
+ mat->A * mu_iter + mat->B * mu_iter * mu_iter - P;
+ dP_c_drho_iter = (mat->a + mat->b * w_iter_inv) * u +
+ mat->b * u * rho_iter * dw_inv_drho_iter +
+ mat->A * dmu_drho_iter + mat->B * dmu_sq_drho_iter;
+ P_fraction = P_c_iter / P;
+
+ // if curve is cold then we've got everything we need
+ if (curve == Til_region_cold) {
+ rho_iter -= P_c_iter / dP_c_drho_iter;
+ // Don't use these:
+ P_h_iter = 0.f;
+ dP_h_drho_iter = 0.f;
+ }
+ // if "cold_min" or "hybrid_min"
+ // Can only have one version of the cold curve, therefore either use the
+ // min version or the normal version hence "else if"
+ } else if (curve == Til_region_cold_min ||
+ curve == Til_region_hybrid_min) {
+ P_c_iter = ((mat->a + mat->b * w_iter_inv) * rho_iter * u +
+ mat->A * mu_iter + mat->B * mu_iter * mu_iter) *
+ (eta_iter - mat->eta_zero) /
+ (mat->eta_min - mat->eta_zero) -
+ P;
+ dP_c_drho_iter =
+ ((mat->a + mat->b * w_iter_inv) * u +
+ mat->b * u * rho_iter * dw_inv_drho_iter + mat->A * dmu_drho_iter +
+ mat->B * dmu_sq_drho_iter) *
+ (eta_iter - mat->eta_zero) / (mat->eta_min - mat->eta_zero) +
+ ((mat->a + mat->b * w_iter_inv) * rho_iter * u + mat->A * mu_iter +
+ mat->B * mu_iter * mu_iter) *
+ (1 / (mat->rho_0 * (mat->eta_min - mat->eta_zero)));
+ P_fraction = P_c_iter / P;
+
+ // if curve is cold_min then we've got everything we need
+ if (curve == Til_region_cold_min) {
+ rho_iter -= P_c_iter / dP_c_drho_iter;
+ // Don't use these:
+ P_c_iter = 0.f;
+ dP_c_drho_iter = 0.f;
+ }
+ }
+
+ // if "hybrid_min" or "hybrid" or "hot"
+ if (curve == Til_region_hybrid_min || curve == Til_region_hybrid ||
+ curve == Til_region_hot) {
+ P_h_iter =
+ mat->a * rho_iter * u +
+ (mat->b * rho_iter * u * w_iter_inv + mat->A * mu_iter * exp1) *
+ exp2 -
+ P;
+ dP_h_drho_iter =
+ mat->a * u +
+ (mat->b * u * w_iter_inv +
+ mat->b * u * rho_iter * dw_inv_drho_iter +
+ mat->A * mu_iter * dexp1_drho_iter +
+ mat->A * exp1 * dmu_drho_iter) *
+ exp2 +
+ (mat->b * rho_iter * u * w_iter_inv + mat->A * mu_iter * exp1) *
+ dexp2_drho_iter;
+ P_fraction = P_h_iter / P;
+
+ // if curve is hot then we've got everything we need
+ if (curve == Til_region_hot) {
+ rho_iter -= P_h_iter / dP_h_drho_iter;
+ }
+ }
+
+ // If we are on a hybrid or hybrid_min curve, we combie hot and cold
+ // curves
+ if (curve == Til_region_hybrid_min || curve == Til_region_hybrid) {
+ const float P_hybrid_iter =
+ ((u - mat->u_iv) * P_h_iter + (mat->u_cv - u) * P_c_iter) /
+ (mat->u_cv - mat->u_iv);
+ const float dP_hybrid_drho_iter = ((u - mat->u_iv) * dP_h_drho_iter +
+ (mat->u_cv - u) * dP_c_drho_iter) /
+ (mat->u_cv - mat->u_iv);
+ rho_iter -= P_hybrid_iter / dP_hybrid_drho_iter;
+ P_fraction = P_hybrid_iter / P;
+ }
+
+ // Now we have to constrain the new rho_iter to the curve we're on
+ // u REGION 1
+ if (u <= mat->u_iv) {
+ if (curve == Til_region_cold_min) {
+ if (rho_iter > mat->eta_min * mat->rho_0) {
+ rho_iter = mat->eta_min * mat->rho_0;
+ }
+ } else if (curve == Til_region_cold) {
+ if (rho_iter < mat->eta_min * mat->rho_0) {
+ rho_iter = mat->eta_min * mat->rho_0;
+ }
+ } else {
+ error("Error in Til_density_from_pressure_and_internal_energy");
+ }
+ // u REGION 2
+ } else if (u <= mat->u_cv) {
+ if (curve == Til_region_hybrid_min) {
+ if (rho_iter > mat->eta_min * mat->rho_0) {
+ rho_iter = mat->eta_min * mat->rho_0;
+ }
+ } else if (curve == Til_region_hybrid) {
+ if (rho_iter < mat->eta_min * mat->rho_0) {
+ rho_iter = mat->eta_min * mat->rho_0;
+ }
+ if (rho_iter > mat->rho_0) {
+ rho_iter = mat->rho_0;
+ }
+ } else if (curve == Til_region_cold) {
+ if (rho_iter < mat->rho_0) {
+ rho_iter = mat->rho_0;
+ }
+ } else {
+ error("Error in Til_density_from_pressure_and_internal_energy");
+ }
+
+ // u REGION 3
+ } else {
+ if (curve == Til_region_hot) {
+ if (rho_iter > mat->rho_0) {
+ rho_iter = mat->rho_0;
+ }
+ } else if (curve == Til_region_cold) {
+ if (rho_iter < mat->rho_0) {
+ rho_iter = mat->rho_0;
+ }
+ } else {
+ error("Error in Til_density_from_pressure_and_internal_energy");
+ }
+ }
+
+ // Either we've converged ...
+ if (fabsf(P_fraction) < tol) {
+ return rho_iter;
+
+ // ... or we're stuck at the boundary ...
+ } else if (rho_iter == last_rho_iter) {
+ break;
+ }
+
+ // ... or we loop again
+ last_rho_iter = rho_iter;
+ }
+ }
+ return rho_sph;
+}
+
#endif /* SWIFT_TILLOTSON_EQUATION_OF_STATE_H */
diff --git a/src/utilities.h b/src/utilities.h
index 28faf6b1e59159af4ee4a7b87054ab6a36478c9c..747ef2154ed2d0b09c22208d3c88aa8a15fcbd3c 100644
--- a/src/utilities.h
+++ b/src/utilities.h
@@ -56,4 +56,47 @@ INLINE static int find_value_in_monot_incr_array(const float x,
return index_low;
}
+/**
+ * @brief Search for a value in a monotonically increasing array to find the
+ * index such that table[i,j] = array[i*n_col + j] < value <
+ * table[i + 1, j] = array[(i + 1)*n_col + j]
+ *
+ * @param x The value to find
+ * @param array The array to search
+ * @param n_row The number of rows of the table
+ * @param n_col The number of columns of the table
+ * @param j The column index to perform the search
+ *
+ * Return -1 and n for x below and above the array edge values respectively.
+ */
+INLINE static int vertical_find_value_in_monot_incr_array(const float x,
+ const float *array,
+ const int n_row,
+ const int n_col,
+ const int j) {
+
+ int i_low = 0;
+ int i_high = n_row - 1;
+ int i_mid;
+
+ // Until table[i_low,j] < x < table[i_high=i_low + 1, j]
+ while (i_high - i_low > 1) {
+ i_mid = (i_high + i_low) / 2; // Middle index
+
+ // Replace the low or high i with the middle
+ if (array[i_mid * n_col + j] <= x)
+ i_low = i_mid;
+ else
+ i_high = i_mid;
+ }
+
+ // Set index with the found i_low or an error value if outside the array
+ if (x < array[j])
+ return -1;
+ else if (array[(n_row - 1) * n_col + j] <= x)
+ return n_row;
+ else
+ return i_low;
+}
+
#endif /* SWIFT_UTILITIES_H */