Skip to content
Snippets Groups Projects
Commit 73166e84 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'planetary_remix_merge_eos' into 'master'

Add new planetary equations of state and EoS functions

See merge request !2115
parents 4e3e41be 8f8a50c8
Branches
Tags
1 merge request!2115Add new planetary equations of state and EoS functions
Showing with 2537 additions and 272 deletions
...@@ -36,6 +36,8 @@ planetary initial conditions, `WoMa <https://github.com/srbonilla/WoMa>`_: ...@@ -36,6 +36,8 @@ planetary initial conditions, `WoMa <https://github.com/srbonilla/WoMa>`_:
+ Granite: ``101`` + Granite: ``101``
+ Water: ``102`` + Water: ``102``
+ Basalt: ``103`` + Basalt: ``103``
+ Ice: ``104``
+ Custom user-provided parameters: ``190``, ``191``, ..., ``199``
+ Hubbard \& MacFarlane (1980): ``2`` + Hubbard \& MacFarlane (1980): ``2``
+ Hydrogen-helium atmosphere: ``200`` + Hydrogen-helium atmosphere: ``200``
+ Ice H20-CH4-NH3 mix: ``201`` + Ice H20-CH4-NH3 mix: ``201``
...@@ -45,6 +47,10 @@ planetary initial conditions, `WoMa <https://github.com/srbonilla/WoMa>`_: ...@@ -45,6 +47,10 @@ planetary initial conditions, `WoMa <https://github.com/srbonilla/WoMa>`_:
+ Basalt (7530): ``301`` + Basalt (7530): ``301``
+ Water (7154): ``302`` + Water (7154): ``302``
+ Senft \& Stewart (2008) water: ``303`` + 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`` + ANEOS (in SESAME-style tables): ``4``
+ Forsterite (Stewart et al. 2019): ``400`` + Forsterite (Stewart et al. 2019): ``400``
+ Iron (Stewart, zenodo.org/record/3866507): ``401`` + Iron (Stewart, zenodo.org/record/3866507): ``401``
...@@ -56,7 +62,7 @@ The data files for the tabulated EoS can be downloaded using ...@@ -56,7 +62,7 @@ The data files for the tabulated EoS can be downloaded using
the ``examples/Planetary/EoSTables/get_eos_tables.sh`` script. the ``examples/Planetary/EoSTables/get_eos_tables.sh`` script.
To enable one or multiple EoS, the corresponding ``planetary_use_*:`` 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 along with the path to any table files, which are set by the
``planetary_*_table_file:`` parameters, ``planetary_*_table_file:`` parameters,
as detailed in :ref:`Parameters_eos` and ``examples/parameter_example.yml``. as detailed in :ref:`Parameters_eos` and ``examples/parameter_example.yml``.
......
...@@ -7,6 +7,10 @@ wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/SESAME_basalt_7530.txt ...@@ -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_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/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/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_forsterite_S19.txt
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/ANEOS_iron_S20.txt wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/ANEOS_iron_S20.txt
......
...@@ -314,6 +314,17 @@ EoS: ...@@ -314,6 +314,17 @@ EoS:
planetary_use_Til_granite: 1 # Tillotson granite, material ID 101 planetary_use_Til_granite: 1 # Tillotson granite, material ID 101
planetary_use_Til_water: 0 # Tillotson water, material ID 102 planetary_use_Til_water: 0 # Tillotson water, material ID 102
planetary_use_Til_basalt: 0 # Tillotson basalt, material ID 103 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_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_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 planetary_use_HM80_rock: 0 # Hubbard & MacFarlane (1980) SiO2-MgO-FeS-FeO rock mix, material ID 202
...@@ -321,6 +332,10 @@ EoS: ...@@ -321,6 +332,10 @@ EoS:
planetary_use_SESAME_basalt: 0 # SESAME basalt 7530, material ID 301 planetary_use_SESAME_basalt: 0 # SESAME basalt 7530, material ID 301
planetary_use_SESAME_water: 0 # SESAME water 7154, material ID 302 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_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_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_iron: 0 # ANEOS iron (Stewart 2020), material ID 401
planetary_use_ANEOS_Fe85Si15: 0 # ANEOS Fe85Si15 (Stewart 2020), material ID 402 planetary_use_ANEOS_Fe85Si15: 0 # ANEOS Fe85Si15 (Stewart 2020), material ID 402
...@@ -335,6 +350,16 @@ EoS: ...@@ -335,6 +350,16 @@ EoS:
planetary_use_custom_8: 0 planetary_use_custom_8: 0
planetary_use_custom_9: 0 planetary_use_custom_9: 0
# Tablulated EoS file paths. # 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_HHe_table_file: ./EoSTables/HM80_HHe.txt
planetary_HM80_ice_table_file: ./EoSTables/HM80_ice.txt planetary_HM80_ice_table_file: ./EoSTables/HM80_ice.txt
planetary_HM80_rock_table_file: ./EoSTables/HM80_rock.txt planetary_HM80_rock_table_file: ./EoSTables/HM80_rock.txt
......
...@@ -56,7 +56,7 @@ struct eos_parameters { ...@@ -56,7 +56,7 @@ struct eos_parameters {
* @param entropy The entropy \f$A\f$. * @param entropy The entropy \f$A\f$.
*/ */
__attribute__((always_inline, const)) INLINE static float __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_frac = density * eos.inverse_core_density;
const float density_factor = pow_gamma(density_frac); const float density_factor = pow_gamma(density_frac);
...@@ -76,7 +76,7 @@ gas_internal_energy_from_entropy(float density, float entropy) { ...@@ -76,7 +76,7 @@ gas_internal_energy_from_entropy(float density, float entropy) {
* @param entropy The entropy \f$A\f$. * @param entropy The entropy \f$A\f$.
*/ */
__attribute__((always_inline, const)) INLINE static float __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_frac = density * eos.inverse_core_density;
const float density_factor = pow_gamma(density_frac); const float density_factor = pow_gamma(density_frac);
...@@ -96,7 +96,7 @@ gas_pressure_from_entropy(float density, float entropy) { ...@@ -96,7 +96,7 @@ gas_pressure_from_entropy(float density, float entropy) {
* @return The entropy \f$A\f$. * @return The entropy \f$A\f$.
*/ */
__attribute__((always_inline, const)) INLINE static float __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_frac = density * eos.inverse_core_density;
const float density_factor = pow_gamma(density_frac); const float density_factor = pow_gamma(density_frac);
...@@ -116,7 +116,7 @@ gas_entropy_from_pressure(float density, float pressure) { ...@@ -116,7 +116,7 @@ gas_entropy_from_pressure(float density, float pressure) {
* @param entropy The entropy \f$A\f$. * @param entropy The entropy \f$A\f$.
*/ */
__attribute__((always_inline, const)) INLINE static float __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_frac = density * eos.inverse_core_density;
const float density_factor = pow_gamma(density_frac); const float density_factor = pow_gamma(density_frac);
...@@ -135,7 +135,7 @@ gas_soundspeed_from_entropy(float density, float entropy) { ...@@ -135,7 +135,7 @@ gas_soundspeed_from_entropy(float density, float entropy) {
* @param u The internal energy \f$u\f$ * @param u The internal energy \f$u\f$
*/ */
__attribute__((always_inline, const)) INLINE static float __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_frac = density * eos.inverse_core_density;
const float density_factor = pow_gamma(density_frac); const float density_factor = pow_gamma(density_frac);
...@@ -155,7 +155,7 @@ gas_entropy_from_internal_energy(float density, float u) { ...@@ -155,7 +155,7 @@ gas_entropy_from_internal_energy(float density, float u) {
* @param u The internal energy \f$u\f$ * @param u The internal energy \f$u\f$
*/ */
__attribute__((always_inline, const)) INLINE static float __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_frac = density * eos.inverse_core_density;
const float density_factor = pow_gamma(density_frac); const float density_factor = pow_gamma(density_frac);
...@@ -175,7 +175,7 @@ gas_pressure_from_internal_energy(float density, float u) { ...@@ -175,7 +175,7 @@ gas_pressure_from_internal_energy(float density, float u) {
* @return The internal energy \f$u\f$. * @return The internal energy \f$u\f$.
*/ */
__attribute__((always_inline, const)) INLINE static float __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_frac = density * eos.inverse_core_density;
const float density_factor = pow_gamma(density_frac); const float density_factor = pow_gamma(density_frac);
...@@ -195,7 +195,7 @@ gas_internal_energy_from_pressure(float density, float pressure) { ...@@ -195,7 +195,7 @@ gas_internal_energy_from_pressure(float density, float pressure) {
* @param u The internal energy \f$u\f$ * @param u The internal energy \f$u\f$
*/ */
__attribute__((always_inline, const)) INLINE static float __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_frac = density * eos.inverse_core_density;
const float density_factor = pow_gamma(density_frac); const float density_factor = pow_gamma(density_frac);
...@@ -214,7 +214,7 @@ gas_soundspeed_from_internal_energy(float density, float u) { ...@@ -214,7 +214,7 @@ gas_soundspeed_from_internal_energy(float density, float u) {
* @param P The pressure \f$P\f$ * @param P The pressure \f$P\f$
*/ */
__attribute__((always_inline, const)) INLINE static float __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_frac = density * eos.inverse_core_density;
const float density_factor = pow_gamma(density_frac); const float density_factor = pow_gamma(density_frac);
...@@ -241,7 +241,7 @@ INLINE static void eos_init(struct eos_parameters *e, ...@@ -241,7 +241,7 @@ INLINE static void eos_init(struct eos_parameters *e,
parser_get_param_float(params, "EoS:barotropic_vacuum_sound_speed"); parser_get_param_float(params, "EoS:barotropic_vacuum_sound_speed");
e->vacuum_sound_speed2 = vacuum_sound_speed * vacuum_sound_speed; e->vacuum_sound_speed2 = vacuum_sound_speed * vacuum_sound_speed;
e->inverse_core_density = 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 * @brief Print the equation of state
......
...@@ -46,7 +46,7 @@ struct eos_parameters {}; ...@@ -46,7 +46,7 @@ struct eos_parameters {};
* @param entropy The entropy \f$A\f$. * @param entropy The entropy \f$A\f$.
*/ */
__attribute__((always_inline, const)) INLINE static float __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) * return entropy * pow_gamma_minus_one(density) *
hydro_one_over_gamma_minus_one; hydro_one_over_gamma_minus_one;
...@@ -61,7 +61,7 @@ gas_internal_energy_from_entropy(float density, float entropy) { ...@@ -61,7 +61,7 @@ gas_internal_energy_from_entropy(float density, float entropy) {
* @param entropy The entropy \f$A\f$. * @param entropy The entropy \f$A\f$.
*/ */
__attribute__((always_inline, const)) INLINE static float __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); return entropy * pow_gamma(density);
} }
...@@ -76,7 +76,7 @@ gas_pressure_from_entropy(float density, float entropy) { ...@@ -76,7 +76,7 @@ gas_pressure_from_entropy(float density, float entropy) {
* @return The entropy \f$A\f$. * @return The entropy \f$A\f$.
*/ */
__attribute__((always_inline, const)) INLINE static float __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); return pressure * pow_minus_gamma(density);
} }
...@@ -90,7 +90,7 @@ gas_entropy_from_pressure(float density, float pressure) { ...@@ -90,7 +90,7 @@ gas_entropy_from_pressure(float density, float pressure) {
* @param entropy The entropy \f$A\f$. * @param entropy The entropy \f$A\f$.
*/ */
__attribute__((always_inline, const)) INLINE static float __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); return sqrtf(hydro_gamma * pow_gamma_minus_one(density) * entropy);
} }
...@@ -104,7 +104,7 @@ gas_soundspeed_from_entropy(float density, float entropy) { ...@@ -104,7 +104,7 @@ gas_soundspeed_from_entropy(float density, float entropy) {
* @param u The internal energy \f$u\f$ * @param u The internal energy \f$u\f$
*/ */
__attribute__((always_inline, const)) INLINE static float __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); 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) { ...@@ -118,7 +118,7 @@ gas_entropy_from_internal_energy(float density, float u) {
* @param u The internal energy \f$u\f$ * @param u The internal energy \f$u\f$
*/ */
__attribute__((always_inline, const)) INLINE static float __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; return hydro_gamma_minus_one * u * density;
} }
...@@ -133,7 +133,7 @@ gas_pressure_from_internal_energy(float density, float u) { ...@@ -133,7 +133,7 @@ gas_pressure_from_internal_energy(float density, float u) {
* @return The internal energy \f$u\f$. * @return The internal energy \f$u\f$.
*/ */
__attribute__((always_inline, const)) INLINE static float __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; return hydro_one_over_gamma_minus_one * pressure / density;
} }
...@@ -147,7 +147,7 @@ gas_internal_energy_from_pressure(float density, float pressure) { ...@@ -147,7 +147,7 @@ gas_internal_energy_from_pressure(float density, float pressure) {
* @param u The internal energy \f$u\f$ * @param u The internal energy \f$u\f$
*/ */
__attribute__((always_inline, const)) INLINE static float __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); return sqrtf(u * hydro_gamma * hydro_gamma_minus_one);
} }
...@@ -161,7 +161,7 @@ gas_soundspeed_from_internal_energy(float density, float u) { ...@@ -161,7 +161,7 @@ gas_soundspeed_from_internal_energy(float density, float u) {
* @param P The pressure \f$P\f$ * @param P The pressure \f$P\f$
*/ */
__attribute__((always_inline, const)) INLINE static float __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); return sqrtf(hydro_gamma * P / density);
} }
......
...@@ -50,7 +50,7 @@ struct eos_parameters { ...@@ -50,7 +50,7 @@ struct eos_parameters {
* @param entropy The entropy \f$S\f$. * @param entropy The entropy \f$S\f$.
*/ */
__attribute__((always_inline)) INLINE static float __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; return eos.isothermal_internal_energy;
} }
...@@ -65,7 +65,7 @@ gas_internal_energy_from_entropy(float density, float entropy) { ...@@ -65,7 +65,7 @@ gas_internal_energy_from_entropy(float density, float entropy) {
* @param entropy The entropy \f$S\f$. * @param entropy The entropy \f$S\f$.
*/ */
__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy( __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; return hydro_gamma_minus_one * eos.isothermal_internal_energy * density;
} }
...@@ -81,7 +81,7 @@ __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy( ...@@ -81,7 +81,7 @@ __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
* @return The entropy \f$A\f$. * @return The entropy \f$A\f$.
*/ */
__attribute__((always_inline)) INLINE static float gas_entropy_from_pressure( __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 * return hydro_gamma_minus_one * eos.isothermal_internal_energy *
pow_minus_gamma_minus_one(density); pow_minus_gamma_minus_one(density);
...@@ -98,7 +98,7 @@ __attribute__((always_inline)) INLINE static float gas_entropy_from_pressure( ...@@ -98,7 +98,7 @@ __attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
* @param entropy The entropy \f$S\f$. * @param entropy The entropy \f$S\f$.
*/ */
__attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy( __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 * return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
hydro_gamma_minus_one); hydro_gamma_minus_one);
...@@ -114,7 +114,7 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy( ...@@ -114,7 +114,7 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
* @param u The internal energy \f$u\f$ * @param u The internal energy \f$u\f$
*/ */
__attribute__((always_inline)) INLINE static float __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 * return hydro_gamma_minus_one * eos.isothermal_internal_energy *
pow_minus_gamma_minus_one(density); pow_minus_gamma_minus_one(density);
...@@ -130,7 +130,7 @@ gas_entropy_from_internal_energy(float density, float u) { ...@@ -130,7 +130,7 @@ gas_entropy_from_internal_energy(float density, float u) {
* @param u The internal energy \f$u\f$ * @param u The internal energy \f$u\f$
*/ */
__attribute__((always_inline)) INLINE static float __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; return hydro_gamma_minus_one * eos.isothermal_internal_energy * density;
} }
...@@ -145,7 +145,7 @@ gas_pressure_from_internal_energy(float density, float u) { ...@@ -145,7 +145,7 @@ gas_pressure_from_internal_energy(float density, float u) {
* @return The internal energy \f$u\f$ (which is constant). * @return The internal energy \f$u\f$ (which is constant).
*/ */
__attribute__((always_inline)) INLINE static float __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; return eos.isothermal_internal_energy;
} }
...@@ -160,7 +160,7 @@ gas_internal_energy_from_pressure(float density, float pressure) { ...@@ -160,7 +160,7 @@ gas_internal_energy_from_pressure(float density, float pressure) {
* @param u The internal energy \f$u\f$ * @param u The internal energy \f$u\f$
*/ */
__attribute__((always_inline)) INLINE static float __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 * return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
hydro_gamma_minus_one); hydro_gamma_minus_one);
...@@ -177,7 +177,7 @@ gas_soundspeed_from_internal_energy(float density, float u) { ...@@ -177,7 +177,7 @@ gas_soundspeed_from_internal_energy(float density, float u) {
* @param P The pressure \f$P\f$ * @param P The pressure \f$P\f$
*/ */
__attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure( __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 * return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
hydro_gamma_minus_one); hydro_gamma_minus_one);
......
This diff is collapsed.
...@@ -179,7 +179,7 @@ INLINE static void convert_units_HM80(struct HM80_params *mat, ...@@ -179,7 +179,7 @@ INLINE static void convert_units_HM80(struct HM80_params *mat,
// gas_internal_energy_from_entropy // gas_internal_energy_from_entropy
INLINE static float HM80_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!"); error("This EOS function is not yet implemented!");
...@@ -187,7 +187,8 @@ INLINE static float HM80_internal_energy_from_entropy( ...@@ -187,7 +187,8 @@ INLINE static float HM80_internal_energy_from_entropy(
} }
// gas_pressure_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) { const struct HM80_params *mat) {
error("This EOS function is not yet implemented!"); error("This EOS function is not yet implemented!");
...@@ -196,7 +197,8 @@ INLINE static float HM80_pressure_from_entropy(float density, float entropy, ...@@ -196,7 +197,8 @@ INLINE static float HM80_pressure_from_entropy(float density, float entropy,
} }
// gas_entropy_from_pressure // 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) { const struct HM80_params *mat) {
error("This EOS function is not yet implemented!"); error("This EOS function is not yet implemented!");
...@@ -206,7 +208,7 @@ INLINE static float HM80_entropy_from_pressure(float density, float pressure, ...@@ -206,7 +208,7 @@ INLINE static float HM80_entropy_from_pressure(float density, float pressure,
// gas_soundspeed_from_entropy // gas_soundspeed_from_entropy
INLINE static float HM80_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!"); error("This EOS function is not yet implemented!");
...@@ -215,29 +217,25 @@ INLINE static float HM80_soundspeed_from_entropy( ...@@ -215,29 +217,25 @@ INLINE static float HM80_soundspeed_from_entropy(
// gas_entropy_from_internal_energy // gas_entropy_from_internal_energy
INLINE static float HM80_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; return 0.f;
} }
// gas_pressure_from_internal_energy // gas_pressure_from_internal_energy
INLINE static float HM80_pressure_from_internal_energy( INLINE static float HM80_pressure_from_internal_energy(
float density, float u, const struct HM80_params *mat) { const float density, const float u, const struct HM80_params *mat) {
float log_P, log_P_1, log_P_2, log_P_3, log_P_4;
if (u <= 0.f) { if (u <= 0.f) {
return 0.f; return 0.f;
} }
int idx_rho, idx_u;
float intp_rho, intp_u;
const float log_rho = logf(density); const float log_rho = logf(density);
const float log_u = logf(u); 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)
idx_rho = floor((log_rho - mat->log_rho_min) * mat->inv_log_rho_step); int 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_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 outside the table then extrapolate from the edge and edge-but-one values
if (idx_rho <= -1) { if (idx_rho <= -1) {
...@@ -251,26 +249,31 @@ INLINE static float HM80_pressure_from_internal_energy( ...@@ -251,26 +249,31 @@ INLINE static float HM80_pressure_from_internal_energy(
idx_u = mat->num_u - 2; idx_u = mat->num_u - 2;
} }
intp_rho = (log_rho - mat->log_rho_min - idx_rho * mat->log_rho_step) * const float intp_rho =
mat->inv_log_rho_step; (log_rho - mat->log_rho_min - idx_rho * mat->log_rho_step) *
intp_u = 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; (log_u - mat->log_u_min - idx_u * mat->log_u_step) * mat->inv_log_u_step;
// Table values // Table values
log_P_1 = mat->table_log_P_rho_u[idx_rho * mat->num_u + idx_u]; const float 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]; const float log_P_2 =
log_P_3 = mat->table_log_P_rho_u[(idx_rho + 1) * mat->num_u + idx_u]; mat->table_log_P_rho_u[idx_rho * mat->num_u + idx_u + 1];
log_P_4 = mat->table_log_P_rho_u[(idx_rho + 1) * 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];
log_P = (1.f - intp_rho) * ((1.f - intp_u) * log_P_1 + intp_u * log_P_2) + const float log_P_4 =
intp_rho * ((1.f - intp_u) * log_P_3 + intp_u * 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); return expf(log_P);
} }
// gas_internal_energy_from_pressure // gas_internal_energy_from_pressure
INLINE static float HM80_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!"); error("This EOS function is not yet implemented!");
...@@ -279,9 +282,9 @@ INLINE static float HM80_internal_energy_from_pressure( ...@@ -279,9 +282,9 @@ INLINE static float HM80_internal_energy_from_pressure(
// gas_soundspeed_from_internal_energy // gas_soundspeed_from_internal_energy
INLINE static float HM80_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 // Bulk modulus
if (mat->bulk_mod != 0) { if (mat->bulk_mod != 0) {
...@@ -289,10 +292,10 @@ INLINE static float HM80_soundspeed_from_internal_energy( ...@@ -289,10 +292,10 @@ INLINE static float HM80_soundspeed_from_internal_energy(
} }
// Ideal gas // Ideal gas
else { 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); c = sqrtf(hydro_gamma * P / density);
if (c <= 0) { if (c <= 0.f) {
c = sqrtf(hydro_gamma * mat->P_min_for_c_min / density); c = sqrtf(hydro_gamma * mat->P_min_for_c_min / density);
} }
} }
...@@ -302,7 +305,35 @@ INLINE static float HM80_soundspeed_from_internal_energy( ...@@ -302,7 +305,35 @@ INLINE static float HM80_soundspeed_from_internal_energy(
// gas_soundspeed_from_pressure // gas_soundspeed_from_pressure
INLINE static float HM80_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!"); error("This EOS function is not yet implemented!");
......
...@@ -58,7 +58,7 @@ INLINE static void set_idg_def(struct idg_params *mat, ...@@ -58,7 +58,7 @@ INLINE static void set_idg_def(struct idg_params *mat,
* @param entropy The entropy \f$A\f$. * @param entropy The entropy \f$A\f$.
*/ */
INLINE static float idg_internal_energy_from_entropy( 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) * return entropy * powf(density, mat->gamma - 1.f) *
mat->one_over_gamma_minus_one; mat->one_over_gamma_minus_one;
...@@ -72,7 +72,8 @@ INLINE static float idg_internal_energy_from_entropy( ...@@ -72,7 +72,8 @@ INLINE static float idg_internal_energy_from_entropy(
* @param density The density \f$\rho\f$. * @param density The density \f$\rho\f$.
* @param entropy The entropy \f$A\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) { const struct idg_params *mat) {
return entropy * powf(density, mat->gamma); return entropy * powf(density, mat->gamma);
...@@ -87,7 +88,8 @@ INLINE static float idg_pressure_from_entropy(float density, float entropy, ...@@ -87,7 +88,8 @@ INLINE static float idg_pressure_from_entropy(float density, float entropy,
* @param pressure The pressure \f$P\f$. * @param pressure The pressure \f$P\f$.
* @return The entropy \f$A\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) { const struct idg_params *mat) {
return pressure * powf(density, -mat->gamma); return pressure * powf(density, -mat->gamma);
...@@ -101,7 +103,8 @@ INLINE static float idg_entropy_from_pressure(float density, float pressure, ...@@ -101,7 +103,8 @@ INLINE static float idg_entropy_from_pressure(float density, float pressure,
* @param density The density \f$\rho\f$. * @param density The density \f$\rho\f$.
* @param entropy The entropy \f$A\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) { const struct idg_params *mat) {
return sqrtf(mat->gamma * powf(density, mat->gamma - 1.f) * entropy); 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, ...@@ -116,7 +119,7 @@ INLINE static float idg_soundspeed_from_entropy(float density, float entropy,
* @param u The internal energy \f$u\f$ * @param u The internal energy \f$u\f$
*/ */
INLINE static float idg_entropy_from_internal_energy( 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); return (mat->gamma - 1.f) * u * powf(density, 1.f - mat->gamma);
} }
...@@ -130,7 +133,7 @@ INLINE static float idg_entropy_from_internal_energy( ...@@ -130,7 +133,7 @@ INLINE static float idg_entropy_from_internal_energy(
* @param u The internal energy \f$u\f$ * @param u The internal energy \f$u\f$
*/ */
INLINE static float idg_pressure_from_internal_energy( 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; return (mat->gamma - 1.f) * u * density;
} }
...@@ -145,7 +148,7 @@ INLINE static float idg_pressure_from_internal_energy( ...@@ -145,7 +148,7 @@ INLINE static float idg_pressure_from_internal_energy(
* @return The internal energy \f$u\f$. * @return The internal energy \f$u\f$.
*/ */
INLINE static float idg_internal_energy_from_pressure( 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; return mat->one_over_gamma_minus_one * pressure / density;
} }
...@@ -159,7 +162,7 @@ INLINE static float idg_internal_energy_from_pressure( ...@@ -159,7 +162,7 @@ INLINE static float idg_internal_energy_from_pressure(
* @param u The internal energy \f$u\f$ * @param u The internal energy \f$u\f$
*/ */
INLINE static float idg_soundspeed_from_internal_energy( 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)); return sqrtf(u * mat->gamma * (mat->gamma - 1.f));
} }
...@@ -172,10 +175,37 @@ INLINE static float idg_soundspeed_from_internal_energy( ...@@ -172,10 +175,37 @@ INLINE static float idg_soundspeed_from_internal_energy(
* @param density The density \f$\rho\f$ * @param density The density \f$\rho\f$
* @param P The pressure \f$P\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) { const struct idg_params *mat) {
return sqrtf(mat->gamma * P / density); 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 */ #endif /* SWIFT_IDEAL_GAS_EQUATION_OF_STATE_H */
This diff is collapsed.
This diff is collapsed.
...@@ -56,4 +56,47 @@ INLINE static int find_value_in_monot_incr_array(const float x, ...@@ -56,4 +56,47 @@ INLINE static int find_value_in_monot_incr_array(const float x,
return index_low; 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 */ #endif /* SWIFT_UTILITIES_H */
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment