diff --git a/AUTHORS b/AUTHORS index 5520d77633752e7ea863766984e516a67cef5a87..aaa5886c77432a291e0ae5e710256168a30a3f56 100644 --- a/AUTHORS +++ b/AUTHORS @@ -25,4 +25,5 @@ Sylvia Ploeckinger ploeckinger@lorentz.leidenuniv.nl Willem Elbers willem.h.elbers@durham.ac.uk TK Chan chantsangkeung@gmail.com Marcel van Daalen daalen@strw.leidenuniv.nl -Filip Husko filip.husko@durham.ac.uk \ No newline at end of file +Filip Husko filip.husko@durham.ac.uk +Orestis Karapiperis karapiperis@lorentz.leidenuniv.nl diff --git a/configure.ac b/configure.ac index 0eac6e9da29d0d4a8a909bdcabd59256c041b651..ef684d5d6ef1c51d8c5f4ae38191630c3a379087 100644 --- a/configure.ac +++ b/configure.ac @@ -2188,7 +2188,7 @@ esac # Equation of state AC_ARG_WITH([equation-of-state], [AS_HELP_STRING([--with-equation-of-state=<EoS>], - [equation of state @<:@ideal-gas, isothermal-gas, planetary default: ideal-gas@:>@] + [equation of state @<:@ideal-gas, isothermal-gas, barotropic-gas, planetary default: ideal-gas@:>@] )], [with_eos="$withval"], [with_eos="ideal-gas"] @@ -2199,6 +2199,9 @@ case "$with_eos" in ;; isothermal-gas) AC_DEFINE([EOS_ISOTHERMAL_GAS], [1], [Isothermal gas equation of state]) + ;; + barotropic-gas) + AC_DEFINE([EOS_BAROTROPIC_GAS], [1], [Barotropic gas equation of state]) ;; planetary) AC_DEFINE([EOS_PLANETARY], [1], [All planetary equations of state]) diff --git a/doc/RTD/source/EquationOfState/index.rst b/doc/RTD/source/EquationOfState/index.rst index 509497a4891fa79fa9b3762e7134487fe1a80b55..8b378136c8620546cc19285cfa0bda9c81d97429 100644 --- a/doc/RTD/source/EquationOfState/index.rst +++ b/doc/RTD/source/EquationOfState/index.rst @@ -7,20 +7,20 @@ Equations of State ================== -Currently, SWIFT offers two different gas equations of state (EoS) -implemented: ``ideal`` and ``isothermal``; as well as a variety of EoS for -"planetary" materials. The EoS describe the relations between our -main thermodynamical variables: the internal energy per unit mass -(\\(u\\)), the mass density (\\(\\rho\\)), the entropy (\\(A\\)) and -the pressure (\\(P\\)). +Currently, SWIFT offers three different gas equations of state (EoS) +implemented: ``ideal``, ``isothermal``, and ``barotropic``; as well as a variety +of EoS for "planetary" materials. The EoS describe the relations between our +main thermodynamical variables: the internal energy per unit mass :math:`u`, the +mass density :math:`\rho`, the entropy :math:`A` and the pressure :math:`P`. +It is selected af configure time via the option ``--with-equation-of-state``. Gas EoS ------- -We write the adiabatic index as \\(\\gamma \\) and \\( c_s \\) denotes +We write the adiabatic index as :math:`\gamma` and :math:`c_s` denotes the speed of sound. The adiabatic index can be changed at configure time by choosing one of the allowed values of the option -``--with-adiabatic-index``. The default value is \\(\\gamma = 5/3 \\). +``--with-adiabatic-index``. The default value is :math:`\gamma = 5/3`. The tables below give the expression for the thermodynamic quantities on each row entry as a function of the gas density and the @@ -29,27 +29,38 @@ thermodynamical quantity given in the header of each column. .. csv-table:: Ideal Gas :header: "Variable", "A", "u", "P" - "A", "", "\\( \\left( \\gamma - 1 \\right) u \\rho^{1-\\gamma} \\)", "\\(P \\rho^{-\\gamma} \\)" - "u", "\\( A \\frac{ \\rho^{ \\gamma - 1 } }{\\gamma - 1 } \\)", "", "\\(\\frac{1}{\\gamma - 1} \\frac{P}{\\rho}\\)" - "P", "\\( A \\rho^\\gamma \\)", "\\( \\left( \\gamma - 1\\right) u \\rho \\)", "" - "\\(c_s\\)", "\\(\\sqrt{ \\gamma \\rho^{\\gamma - 1} A}\\)", "\\(\\sqrt{ u \\gamma \\left( \\gamma - 1 \\right) } \\)", "\\(\\sqrt{ \\frac{\\gamma P}{\\rho} }\\)" + "A", "", :math:`\left( \gamma - 1 \right) u \rho^{1-\gamma}`, :math:`P \rho^{-\gamma}` + "u", :math:`A \frac{ \rho^{ \gamma - 1 } }{\gamma - 1 }`, "", :math:`\frac{1}{\gamma - 1} \frac{P}{\rho}` + "P", :math:`A \rho^\gamma`, :math:`\left( \gamma - 1\right) u \rho`, "" + :math:`c_s`, :math:`\sqrt{ \gamma \rho^{\gamma - 1} A}`, :math:`\sqrt{ u \gamma \left( \gamma - 1 \right) }`, :math:`\sqrt{ \frac{\gamma P}{\rho} }` .. csv-table:: Isothermal Gas - :header: "Variable", "A", "u", "P" + :header: "Variable", "-", "-", "-" - "A", "", "\\(\\left( \\gamma - 1 \\right) u \\rho^{1-\\gamma}\\)", "" + "A", "", :math:`\left( \gamma - 1 \right) u \rho^{1-\gamma}`, "" "u", "", "const", "" - "P", "", "\\(\\left( \\gamma - 1\\right) u \\rho \\)", "" - "\\( c_s\\)", "", "\\(\\sqrt{ u \\gamma \\left( \\gamma - 1 \\right) } \\)", "" - -Note that when running with an isothermal equation of state, the value -of the tracked thermodynamic variable (e.g. the entropy in a + "P", "", :math:`\left( \gamma - 1\right) u \rho`, "" + :math:`c_s`, "", :math:`\sqrt{ u \gamma \left( \gamma - 1 \right) }`, "" + +.. csv-table:: Barotropic Gas + :header: "Variable", "-", "-", "-" + + "A", "", :math:`\rho^{1-\gamma} c_0^2 \sqrt{1 + \left( \frac{\rho}{\rho_c} \right) }`, "" + "u", "", :math:`\frac{1}(\gamma -1)c_0^2 \sqrt{1 + \left( \frac{\rho}{\rho_c} \right) }`, "" + "P", "", :math:`\rho c_0^2 \sqrt{1 + \left( \frac{\rho}{\rho_c} \right) }`, "" + :math:`c_s`, "", :math:`\sqrt{ c_0^2 \sqrt{1 + \left( \frac{\rho}{\rho_c} \right) }}`, "" + +Note that when running with an isothermal or barotropic equation of state, the +value of the tracked thermodynamic variable (e.g. the entropy in a density-entropy scheme or the internal enegy in a density-energy SPH -formulation) written to the snapshots is meaningless. The pressure, -however, is always correct in all scheme. +formulation) written to the snapshots is meaningless. The pressure, however, is +always correct in all scheme. +For the isothermal equation of state, the internal energy is specified at +runtime via the parameter file. In the case of the barotropic gas, the vacuum +sound speed :math:`c_0` and core density :math:`\rho_c` are similarly specified. Planetary EoS @@ -66,7 +77,7 @@ See :ref:`new_option` for a full list of required changes. You will need to provide an ``equation_of_state.h`` file containing: the definition of ``eos_parameters``, IO functions and transformations between the -different variables: \\(u(\\rho, A)\\), \\(u(\\rho, P)\\), \\(P(\\rho,A)\\), -\\(P(\\rho, u)\\), \\(A(\\rho, P)\\), \\(A(\\rho, u)\\), \\(c_s(\\rho, A)\\), -\\(c_s(\\rho, u)\\) and \\(c_s(\\rho, P)\\). See other equation of state files +different variables: :math:`u(\rho, A)`, :math:`u(\rho, P)`, :math:`P(\rho,A)`, +:math:`P(\rho, u)`, :math:`A(\rho, P)`, :math:`A(\rho, u)`, :math:`c_s(\rho, A)`, +:math:`c_s(\rho, u)` and :math:`c_s(\rho, P)`. See other equation of state files to have implementation details. diff --git a/doc/RTD/source/ParameterFiles/parameter_description.rst b/doc/RTD/source/ParameterFiles/parameter_description.rst index d8a070462bccabcf5c71ebe1f24cdaab843c323d..2086d52cd46df3d8ee41daaafe58b6e1fc0cd50e 100644 --- a/doc/RTD/source/ParameterFiles/parameter_description.rst +++ b/doc/RTD/source/ParameterFiles/parameter_description.rst @@ -1230,7 +1230,9 @@ parameter sets the thermal energy per unit mass. .. code:: YAML EoS: - isothermal_internal_energy: 20.26784 # Thermal energy per unit mass for the case of isothermal equation of state (in internal units). + isothermal_internal_energy: 20.26784 # Thermal energy per unit mass for the case of isothermal equation of state (in internal units). + barotropic_vacuum_sound_speed: 2e4 # Vacuum sound speed in the case of the barotropic equation of state (in internal units). + barotropic_core_density: 1e-13 # Core density in the case of the barotropic equation of state (in internal units). # Select which planetary EoS material(s) to enable for use. planetary_use_idg_def: 0 # Default ideal gas, material ID 0 planetary_use_Til_iron: 1 # Tillotson iron, material ID 100 diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index 5e1cdcd036aee31427f2e94771f201bfa3d2b7e7..c94a480fa1658b3b1c541698171dc686379b1d5c 100644 --- a/examples/parameter_example.yml +++ b/examples/parameter_example.yml @@ -295,6 +295,8 @@ LineOfSight: EoS: isothermal_internal_energy: 20.26784 # Thermal energy per unit mass for the case of isothermal equation of state (in internal units). + barotropic_vacuum_sound_speed: 2e4 # Vacuum sound speed (in internal units) + barotropic_core_density: 1e-13 # Core density (in internal units) # Select which planetary EoS material(s) to enable for use. planetary_use_idg_def: 0 # Default ideal gas, material ID 0 planetary_use_Til_iron: 1 # Tillotson iron, material ID 100 diff --git a/src/Makefile.am b/src/Makefile.am index ac7ea7da7f880a0e7643aa318cdf2c8259f1a550..92640cac8c04a1ca4950ad05719c9f3140a6b592 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -214,7 +214,7 @@ nobase_noinst_HEADERS += gravity/MultiSoftening/gravity.h gravity/MultiSoftening nobase_noinst_HEADERS += gravity/MultiSoftening/gravity_debug.h gravity/MultiSoftening/gravity_part.h nobase_noinst_HEADERS += gravity/MultiSoftening/gravity_csds.h nobase_noinst_HEADERS += equation_of_state.h -nobase_noinst_HEADERS += equation_of_state/ideal_gas/equation_of_state.h equation_of_state/isothermal/equation_of_state.h +nobase_noinst_HEADERS += equation_of_state/ideal_gas/equation_of_state.h equation_of_state/isothermal/equation_of_state.h equation_of_state/barotropic/equation_of_state.h nobase_noinst_HEADERS += signal_velocity.h nobase_noinst_HEADERS += hydro.h hydro_io.h hydro_csds.h hydro_parameters.h nobase_noinst_HEADERS += hydro/None/hydro.h hydro/None/hydro_iact.h hydro/None/hydro_io.h diff --git a/src/equation_of_state.c b/src/equation_of_state.c index 5b1c0dc50b18f321e57a90e00a0a0363d90a2f1e..d6fb9cb64a904d4de02bf87a49acde954fad0a95 100644 --- a/src/equation_of_state.c +++ b/src/equation_of_state.c @@ -35,6 +35,9 @@ struct eos_parameters eos = {.Til_iron.rho_0 = -1.f}; #elif defined(EOS_ISOTHERMAL_GAS) struct eos_parameters eos = {.isothermal_internal_energy = -1.}; +#elif defined(EOS_BAROTROPIC_GAS) +struct eos_parameters eos = {.vacuum_sound_speed2 = -1., + .inverse_core_density = -1.}; #else struct eos_parameters eos; #endif diff --git a/src/equation_of_state.h b/src/equation_of_state.h index 769e3dbb6888f52efc16d52858b26ed0d8150515..a002118020acadbc8d50f6f18013f8e150a6cc55 100644 --- a/src/equation_of_state.h +++ b/src/equation_of_state.h @@ -36,6 +36,8 @@ #include "./equation_of_state/isothermal/equation_of_state.h" #elif defined(EOS_PLANETARY) #include "./equation_of_state/planetary/equation_of_state.h" +#elif defined(EOS_BAROTROPIC_GAS) +#include "./equation_of_state/barotropic/equation_of_state.h" #else #error "Invalid choice of equation of state" #endif diff --git a/src/equation_of_state/barotropic/equation_of_state.h b/src/equation_of_state/barotropic/equation_of_state.h new file mode 100644 index 0000000000000000000000000000000000000000..8528e48664918af181307d5b3bdf5fdbf0575dbb --- /dev/null +++ b/src/equation_of_state/barotropic/equation_of_state.h @@ -0,0 +1,280 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2023 Orestis Karapiperis (karapiperis@lorentz.leidenuniv.nl). + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_BAROTROPIC_GAS_EQUATION_OF_STATE_H +#define SWIFT_BAROTROPIC_GAS_EQUATION_OF_STATE_H + +/* Some standard headers. */ +#include <math.h> + +/* Local headers. */ +#include "adiabatic_index.h" +#include "common_io.h" +#include "inline.h" +#include "physical_constants.h" + +extern struct eos_parameters eos; + +/** + * @brief The parameters of the equation of state for the gas. + * + * Barotropic equation of state from Hennebelle et al., 2008, A&A, 477, 9 + * reimplemented following Pakmor et al., 2011, MNRAS, 418, 1392 + */ +struct eos_parameters { + + /*! Square of barotropic sound speed in vacuum */ + float vacuum_sound_speed2; + + /*! Inverse of the core density */ + float inverse_core_density; +}; + +/** + * @brief Returns the internal energy given density and entropy + * + * Since we are using a barotropic EoS, the entropy value is ignored. + * Computes \f$u = c_0^2 \frac{1 + \left(\frac{\rho}{\rho_c}\right)^\gamma + * }{\gamma - 1}\f$. + * + * @param density The density \f$\rho\f$. + * @param entropy The entropy \f$A\f$. + */ +__attribute__((always_inline, const)) INLINE static float +gas_internal_energy_from_entropy(float density, float entropy) { + + const float density_frac = density * eos.inverse_core_density; + const float density_factor = pow_gamma(density_frac); + + return eos.vacuum_sound_speed2 * sqrtf(1.0f + density_factor) * + hydro_one_over_gamma_minus_one; +} + +/** + * @brief Returns the pressure given density and entropy + * + * Since we are using a barotropic EoS, the entropy value is ignored. + * Computes \f$P = c_0^2 \left(1 + + * \left(\frac{\rho}{\rho_c}\right)^\gamma\right)\f$. + * + * @param density The density \f$\rho\f$. + * @param entropy The entropy \f$A\f$. + */ +__attribute__((always_inline, const)) INLINE static float +gas_pressure_from_entropy(float density, float entropy) { + + const float density_frac = density * eos.inverse_core_density; + const float density_factor = pow_gamma(density_frac); + + return eos.vacuum_sound_speed2 * density * sqrtf(1.0f + density_factor); +} + +/** + * @brief Returns the entropy given density and pressure. + * + * Since we are using a barotropic EoS, the pressure value is ignored. + * Computes \f$A = \rho^{1-\gamma}c_0^2 \left(1 + + * \left(\frac{\rho}{\rho_c}\right)^\gamma\right)\f$. + * + * @param density The density \f$\rho\f$. + * @param pressure The pressure \f$P\f$. + * @return The entropy \f$A\f$. + */ +__attribute__((always_inline, const)) INLINE static float +gas_entropy_from_pressure(float density, float pressure) { + + const float density_frac = density * eos.inverse_core_density; + const float density_factor = pow_gamma(density_frac); + + return eos.vacuum_sound_speed2 * pow_minus_gamma_minus_one(density) * + sqrtf(1.0f + density_factor); +} + +/** + * @brief Returns the sound speed given density and entropy + * + * Since we are using a barotropic EoS, the entropy is ignored. + * Computes \f$c = \sqrt{c_0^2 \left(1 + + * \left(\frac{\rho}{\rho_c}\right)^\gamma\right)}\f$. + * + * @param density The density \f$\rho\f$. + * @param entropy The entropy \f$A\f$. + */ +__attribute__((always_inline, const)) INLINE static float +gas_soundspeed_from_entropy(float density, float entropy) { + + const float density_frac = density * eos.inverse_core_density; + const float density_factor = pow_gamma(density_frac); + + return sqrtf(eos.vacuum_sound_speed2 * density * + sqrtf(1.0f + density_factor)); +} + +/** + * @brief Returns the entropy given density and internal energy + * + * Since we are using a barotropic EoS, the internal energy value is ignored. + * Computes \f$A = \rho^{1-\gamma}c_0^2 \left(1 + + * \left(\frac{\rho}{\rho_c}\right)^\gamma\right)\f$. + * + * @param density The density \f$\rho\f$ + * @param u The internal energy \f$u\f$ + */ +__attribute__((always_inline, const)) INLINE static float +gas_entropy_from_internal_energy(float density, float u) { + + const float density_frac = density * eos.inverse_core_density; + const float density_factor = pow_gamma(density_frac); + + return eos.vacuum_sound_speed2 * pow_minus_gamma_minus_one(density) * + sqrtf(1.0f + density_factor); +} + +/** + * @brief Returns the pressure given density and internal energy + * + * Since we are using a barotropic EoS, the internal energy value is ignored. + * Computes \f$P = c_0^2 \left(1 + + * \left(\frac{\rho}{\rho_c}\right)^\gamma\right)\f$. + * + * @param density The density \f$\rho\f$ + * @param u The internal energy \f$u\f$ + */ +__attribute__((always_inline, const)) INLINE static float +gas_pressure_from_internal_energy(float density, float u) { + + const float density_frac = density * eos.inverse_core_density; + const float density_factor = pow_gamma(density_frac); + + return eos.vacuum_sound_speed2 * density * sqrtf(1.0f + density_factor); +} + +/** + * @brief Returns the internal energy given density and pressure. + * + * Since we are using a barotropic EoS, the pressure value is ignored. + * Computes \f$u = c_0^2 \frac{1 + \left(\frac{\rho}{\rho_c}\right)^\gamma + * }{\gamma - 1}\f$. + * + * @param density The density \f$\rho\f$. + * @param pressure The pressure \f$P\f$. + * @return The internal energy \f$u\f$. + */ +__attribute__((always_inline, const)) INLINE static float +gas_internal_energy_from_pressure(float density, float pressure) { + + const float density_frac = density * eos.inverse_core_density; + const float density_factor = pow_gamma(density_frac); + + return eos.vacuum_sound_speed2 * sqrtf(1.0f + density_factor) * + hydro_one_over_gamma_minus_one; +} + +/** + * @brief Returns the sound speed given density and internal energy + * + * Since we are using a barotropic EoS, the internal energy value is ignored. + * Computes \f$c = \sqrt{c_0^2 \left(1 + + * \left(\frac{\rho}{\rho_c}\right)^\gamma\right)}\f$. + * + * @param density The density \f$\rho\f$ + * @param u The internal energy \f$u\f$ + */ +__attribute__((always_inline, const)) INLINE static float +gas_soundspeed_from_internal_energy(float density, float u) { + + const float density_frac = density * eos.inverse_core_density; + const float density_factor = pow_gamma(density_frac); + + return sqrtf(eos.vacuum_sound_speed2 * density * + sqrtf(1.0f + density_factor)); +} + +/** + * @brief Returns the sound speed given density and pressure + * + * Since we are using a barotropic EoS, the pressure value is ignored. + * Computes \f$c = \sqrt{c_0^2 \left(1 + + * \left(\frac{\rho}{\rho_c}\right)^\gamma\right)}\f$. + * + * @param density The density \f$\rho\f$ + * @param P The pressure \f$P\f$ + */ +__attribute__((always_inline, const)) INLINE static float +gas_soundspeed_from_pressure(float density, float P) { + + const float density_frac = density * eos.inverse_core_density; + const float density_factor = pow_gamma(density_frac); + + return sqrtf(eos.vacuum_sound_speed2 * density * + sqrtf(1.0f + density_factor)); +} + +/** + * @brief Initialize the eos parameters + * + * Read the vacuum sound speed and core density from the parameter file. + * + * @param e The #eos_parameters. + * @param phys_const The physical constants in the internal unit system. + * @param us The internal unit system. + * @param params The parsed parameters. + */ +INLINE static void eos_init(struct eos_parameters *e, + const struct phys_const *phys_const, + const struct unit_system *us, + struct swift_params *params) { + + const float vacuum_sound_speed = + 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"); +} +/** + * @brief Print the equation of state + * + * @param e The #eos_parameters + */ +INLINE static void eos_print(const struct eos_parameters *e) { + + message( + "Equation of state: Barotropic gas with vacuum sound speed set to %f and " + "core density set to %f.", + sqrtf(e->vacuum_sound_speed2), 1. / e->inverse_core_density); + + message("Adiabatic index gamma: %f.", hydro_gamma); +} + +#if defined(HAVE_HDF5) +/** + * @brief Write equation of state information to the snapshot + * + * @param h_grpsph The HDF5 group in which to write + * @param e The #eos_parameters + */ +INLINE static void eos_print_snapshot(hid_t h_grpsph, + const struct eos_parameters *e) { + + io_write_attribute_f(h_grpsph, "Adiabatic index", hydro_gamma); + + io_write_attribute_s(h_grpsph, "Equation of state", "Barotropic gas"); +} +#endif + +#endif /* SWIFT_BAROTROPIC_GAS_EQUATION_OF_STATE_H */