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

Merge branch 'add_barotropic_eos' into 'master'

Add barotropic EoS

See merge request !1773
parents 847f25e6 128f366b
Branches
Tags
1 merge request!1773Add barotropic EoS
...@@ -25,4 +25,5 @@ Sylvia Ploeckinger ploeckinger@lorentz.leidenuniv.nl ...@@ -25,4 +25,5 @@ Sylvia Ploeckinger ploeckinger@lorentz.leidenuniv.nl
Willem Elbers willem.h.elbers@durham.ac.uk Willem Elbers willem.h.elbers@durham.ac.uk
TK Chan chantsangkeung@gmail.com TK Chan chantsangkeung@gmail.com
Marcel van Daalen daalen@strw.leidenuniv.nl Marcel van Daalen daalen@strw.leidenuniv.nl
Filip Husko filip.husko@durham.ac.uk Filip Husko filip.husko@durham.ac.uk
\ No newline at end of file Orestis Karapiperis karapiperis@lorentz.leidenuniv.nl
...@@ -2188,7 +2188,7 @@ esac ...@@ -2188,7 +2188,7 @@ esac
# Equation of state # Equation of state
AC_ARG_WITH([equation-of-state], AC_ARG_WITH([equation-of-state],
[AS_HELP_STRING([--with-equation-of-state=<EoS>], [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="$withval"],
[with_eos="ideal-gas"] [with_eos="ideal-gas"]
...@@ -2199,6 +2199,9 @@ case "$with_eos" in ...@@ -2199,6 +2199,9 @@ case "$with_eos" in
;; ;;
isothermal-gas) isothermal-gas)
AC_DEFINE([EOS_ISOTHERMAL_GAS], [1], [Isothermal gas equation of state]) 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) planetary)
AC_DEFINE([EOS_PLANETARY], [1], [All planetary equations of state]) AC_DEFINE([EOS_PLANETARY], [1], [All planetary equations of state])
......
...@@ -7,20 +7,20 @@ ...@@ -7,20 +7,20 @@
Equations of State Equations of State
================== ==================
Currently, SWIFT offers two different gas equations of state (EoS) Currently, SWIFT offers three different gas equations of state (EoS)
implemented: ``ideal`` and ``isothermal``; as well as a variety of EoS for implemented: ``ideal``, ``isothermal``, and ``barotropic``; as well as a variety
"planetary" materials. The EoS describe the relations between our of EoS for "planetary" materials. The EoS describe the relations between our
main thermodynamical variables: the internal energy per unit mass main thermodynamical variables: the internal energy per unit mass :math:`u`, the
(\\(u\\)), the mass density (\\(\\rho\\)), the entropy (\\(A\\)) and mass density :math:`\rho`, the entropy :math:`A` and the pressure :math:`P`.
the pressure (\\(P\\)). It is selected af configure time via the option ``--with-equation-of-state``.
Gas EoS 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 the speed of sound. The adiabatic index can be changed at configure
time by choosing one of the allowed values of the option 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 The tables below give the expression for the thermodynamic quantities
on each row entry as a function of the gas density and the 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. ...@@ -29,27 +29,38 @@ thermodynamical quantity given in the header of each column.
.. csv-table:: Ideal Gas .. csv-table:: Ideal Gas
:header: "Variable", "A", "u", "P" :header: "Variable", "A", "u", "P"
"A", "", "\\( \\left( \\gamma - 1 \\right) u \\rho^{1-\\gamma} \\)", "\\(P \\rho^{-\\gamma} \\)" "A", "", :math:`\left( \gamma - 1 \right) u \rho^{1-\gamma}`, :math:`P \rho^{-\gamma}`
"u", "\\( A \\frac{ \\rho^{ \\gamma - 1 } }{\\gamma - 1 } \\)", "", "\\(\\frac{1}{\\gamma - 1} \\frac{P}{\\rho}\\)" "u", :math:`A \frac{ \rho^{ \gamma - 1 } }{\gamma - 1 }`, "", :math:`\frac{1}{\gamma - 1} \frac{P}{\rho}`
"P", "\\( A \\rho^\\gamma \\)", "\\( \\left( \\gamma - 1\\right) u \\rho \\)", "" "P", :math:`A \rho^\gamma`, :math:`\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} }\\)" :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 .. 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", "" "u", "", "const", ""
"P", "", "\\(\\left( \\gamma - 1\\right) u \\rho \\)", "" "P", "", :math:`\left( \gamma - 1\right) u \rho`, ""
"\\( c_s\\)", "", "\\(\\sqrt{ u \\gamma \\left( \\gamma - 1 \\right) } \\)", "" :math:`c_s`, "", :math:`\sqrt{ u \gamma \left( \gamma - 1 \right) }`, ""
Note that when running with an isothermal equation of state, the value .. csv-table:: Barotropic Gas
of the tracked thermodynamic variable (e.g. the entropy in a :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 density-entropy scheme or the internal enegy in a density-energy SPH
formulation) written to the snapshots is meaningless. The pressure, formulation) written to the snapshots is meaningless. The pressure, however, is
however, is always correct in all scheme. 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 Planetary EoS
...@@ -66,7 +77,7 @@ See :ref:`new_option` for a full list of required changes. ...@@ -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 You will need to provide an ``equation_of_state.h`` file containing: the
definition of ``eos_parameters``, IO functions and transformations between the definition of ``eos_parameters``, IO functions and transformations between the
different variables: \\(u(\\rho, A)\\), \\(u(\\rho, P)\\), \\(P(\\rho,A)\\), different variables: :math:`u(\rho, A)`, :math:`u(\rho, P)`, :math:`P(\rho,A)`,
\\(P(\\rho, u)\\), \\(A(\\rho, P)\\), \\(A(\\rho, u)\\), \\(c_s(\\rho, A)\\), :math:`P(\rho, u)`, :math:`A(\rho, P)`, :math:`A(\rho, u)`, :math:`c_s(\rho, A)`,
\\(c_s(\\rho, u)\\) and \\(c_s(\\rho, P)\\). See other equation of state files :math:`c_s(\rho, u)` and :math:`c_s(\rho, P)`. See other equation of state files
to have implementation details. to have implementation details.
...@@ -1230,7 +1230,9 @@ parameter sets the thermal energy per unit mass. ...@@ -1230,7 +1230,9 @@ parameter sets the thermal energy per unit mass.
.. code:: YAML .. code:: YAML
EoS: 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. # Select which planetary EoS material(s) to enable for use.
planetary_use_idg_def: 0 # Default ideal gas, material ID 0 planetary_use_idg_def: 0 # Default ideal gas, material ID 0
planetary_use_Til_iron: 1 # Tillotson iron, material ID 100 planetary_use_Til_iron: 1 # Tillotson iron, material ID 100
......
...@@ -295,6 +295,8 @@ LineOfSight: ...@@ -295,6 +295,8 @@ LineOfSight:
EoS: 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 internal units)
barotropic_core_density: 1e-13 # Core density (in internal units)
# Select which planetary EoS material(s) to enable for use. # Select which planetary EoS material(s) to enable for use.
planetary_use_idg_def: 0 # Default ideal gas, material ID 0 planetary_use_idg_def: 0 # Default ideal gas, material ID 0
planetary_use_Til_iron: 1 # Tillotson iron, material ID 100 planetary_use_Til_iron: 1 # Tillotson iron, material ID 100
......
...@@ -214,7 +214,7 @@ nobase_noinst_HEADERS += gravity/MultiSoftening/gravity.h gravity/MultiSoftening ...@@ -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_debug.h gravity/MultiSoftening/gravity_part.h
nobase_noinst_HEADERS += gravity/MultiSoftening/gravity_csds.h nobase_noinst_HEADERS += gravity/MultiSoftening/gravity_csds.h
nobase_noinst_HEADERS += equation_of_state.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 += signal_velocity.h
nobase_noinst_HEADERS += hydro.h hydro_io.h hydro_csds.h hydro_parameters.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 nobase_noinst_HEADERS += hydro/None/hydro.h hydro/None/hydro_iact.h hydro/None/hydro_io.h
......
...@@ -35,6 +35,9 @@ ...@@ -35,6 +35,9 @@
struct eos_parameters eos = {.Til_iron.rho_0 = -1.f}; struct eos_parameters eos = {.Til_iron.rho_0 = -1.f};
#elif defined(EOS_ISOTHERMAL_GAS) #elif defined(EOS_ISOTHERMAL_GAS)
struct eos_parameters eos = {.isothermal_internal_energy = -1.}; 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 #else
struct eos_parameters eos; struct eos_parameters eos;
#endif #endif
......
...@@ -36,6 +36,8 @@ ...@@ -36,6 +36,8 @@
#include "./equation_of_state/isothermal/equation_of_state.h" #include "./equation_of_state/isothermal/equation_of_state.h"
#elif defined(EOS_PLANETARY) #elif defined(EOS_PLANETARY)
#include "./equation_of_state/planetary/equation_of_state.h" #include "./equation_of_state/planetary/equation_of_state.h"
#elif defined(EOS_BAROTROPIC_GAS)
#include "./equation_of_state/barotropic/equation_of_state.h"
#else #else
#error "Invalid choice of equation of state" #error "Invalid choice of equation of state"
#endif #endif
......
/*******************************************************************************
* 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 */
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment