diff --git a/README_planetary b/README_planetary index 4e352ea388f2f26e95c7915556666166cf65622c..ca809a7e5a3668b357e8e0fcc0bd5b2dc45ad86f 100644 --- a/README_planetary +++ b/README_planetary @@ -1,8 +1,17 @@ -SWIFT branch to add equations of state for planetary materials such as iron, -silicate rocks, ices, and atmospheres; as well as allowing different particles -to have different materials for the purpose of e.g. giant impact simulations. +SWIFT 'planetary' branch Jacob Kegerreis + jacob.kegerreis@durham.ac.uk +# Aims: + ++ Implement various equations of state, e.g. Tillotson, ANEOS + ++ Allow multiple SPH materials with different EOS for different particles + ++ Add examples, e.g. canonical Moon-forming giant impact + +# Status: + diff --git a/configure.ac b/configure.ac index e730619ea17ce87df7b5c62225094e1d7d0e34db..bb7c0eb10c37c2af7d899913c3a154ea18666dfe 100644 --- a/configure.ac +++ b/configure.ac @@ -890,7 +890,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 default: ideal-gas@:>@] + [equation of state @<:@ideal-gas, isothermal-gas, tillotson-iron default: ideal-gas@:>@] )], [with_eos="$withval"], [with_eos="ideal-gas"] @@ -902,6 +902,9 @@ case "$with_eos" in isothermal-gas) AC_DEFINE([EOS_ISOTHERMAL_GAS], [1], [Isothermal gas equation of state]) ;; + tillotson-iron) + AC_DEFINE([EOS_TILLOTSON_IRON], [1], [Tillotson iron equation of state]) + ;; *) AC_MSG_ERROR([Unknown equation of state: $with_eos]) ;; diff --git a/src/equation_of_state.h b/src/equation_of_state.h index 195b52514f2acc0c40959e09c088a06f0a411869..1f7bf5d77f6e448ad2d2c5617144166684320c56 100644 --- a/src/equation_of_state.h +++ b/src/equation_of_state.h @@ -34,6 +34,8 @@ #include "./equation_of_state/ideal_gas/equation_of_state.h" #elif defined(EOS_ISOTHERMAL_GAS) #include "./equation_of_state/isothermal/equation_of_state.h" +#elif defined(EOS_TILLOTSON_IRON) +#include "./equation_of_state/tillotson_iron/equation_of_state.h" #else #error "Invalid choice of equation of state" #endif diff --git a/src/equation_of_state/tillotson_iron/equation_of_state.h b/src/equation_of_state/tillotson_iron/equation_of_state.h new file mode 100644 index 0000000000000000000000000000000000000000..54a1e15b6600bbbac8a7b681e0421d1cb5914654 --- /dev/null +++ b/src/equation_of_state/tillotson_iron/equation_of_state.h @@ -0,0 +1,203 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk). + * + * 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_TILLOTSON_IRON_EQUATION_OF_STATE_H +#define SWIFT_TILLOTSON_IRON_EQUATION_OF_STATE_H + +/* Some standard headers. */ +#include <math.h> + +/* Local headers. */ +#include "adiabatic_index.h" +#include "common_io.h" +#include "debug.h" +#include "inline.h" + +extern struct eos_parameters eos; +/* ------------------------------------------------------------------------- */ + +struct eos_parameters {}; + +/** + * @brief Returns the internal energy given density and entropy + * + * Computes \f$u = \frac{S\rho^{\gamma-1} }{\gamma - 1}\f$. + * + * @param density The density \f$\rho\f$. + * @param entropy The entropy \f$S\f$. + */ +__attribute__((always_inline)) INLINE static float +gas_internal_energy_from_entropy(float density, float entropy) { + + return entropy * pow_gamma_minus_one(density) * + hydro_one_over_gamma_minus_one; +} + +/** + * @brief Returns the pressure given density and entropy + * + * Computes \f$P = S\rho^\gamma\f$. + * + * @param density The density \f$\rho\f$. + * @param entropy The entropy \f$S\f$. + */ +__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy( + float density, float entropy) { + + return entropy * pow_gamma(density); +} + +/** + * @brief Returns the entropy given density and pressure. + * + * Computes \f$A = \frac{P}{\rho^\gamma}\f$. + * + * @param density The density \f$\rho\f$. + * @param pressure The pressure \f$P\f$. + * @return The entropy \f$A\f$. + */ +__attribute__((always_inline)) INLINE static float gas_entropy_from_pressure( + float density, float pressure) { + + return pressure * pow_minus_gamma(density); +} + +/** + * @brief Returns the sound speed given density and entropy + * + * Computes \f$c = \sqrt{\gamma S \rho^{\gamma-1}}\f$. + * + * @param density The density \f$\rho\f$. + * @param entropy The entropy \f$S\f$. + */ +__attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy( + float density, float entropy) { + + return sqrtf(hydro_gamma * pow_gamma_minus_one(density) * entropy); +} + +/** + * @brief Returns the entropy given density and internal energy + * + * Computes \f$S = \frac{(\gamma - 1)u}{\rho^{\gamma-1}}\f$. + * + * @param density The density \f$\rho\f$ + * @param u The internal energy \f$u\f$ + */ +__attribute__((always_inline)) INLINE static float +gas_entropy_from_internal_energy(float density, float u) { + + return hydro_gamma_minus_one * u * pow_minus_gamma_minus_one(density); +} + +/** + * @brief Returns the pressure given density and internal energy + * + * Computes \f$P = (\gamma - 1)u\rho\f$. + * + * @param density The density \f$\rho\f$ + * @param u The internal energy \f$u\f$ + */ +__attribute__((always_inline)) INLINE static float +gas_pressure_from_internal_energy(float density, float u) { + + return hydro_gamma_minus_one * u * density; +} + +/** + * @brief Returns the internal energy given density and pressure. + * + * Computes \f$u = \frac{1}{\gamma - 1}\frac{P}{\rho}\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)) INLINE static float +gas_internal_energy_from_pressure(float density, float pressure) { + return hydro_one_over_gamma_minus_one * pressure / density; +} + +/** + * @brief Returns the sound speed given density and internal energy + * + * Computes \f$c = \sqrt{\gamma (\gamma - 1) u }\f$. + * + * @param density The density \f$\rho\f$ + * @param u The internal energy \f$u\f$ + */ +__attribute__((always_inline)) INLINE static float +gas_soundspeed_from_internal_energy(float density, float u) { + + return sqrtf(u * hydro_gamma * hydro_gamma_minus_one); +} + +/** + * @brief Returns the sound speed given density and pressure + * + * Computes \f$c = \sqrt{\frac{\gamma P}{\rho} }\f$. + * + * @param density The density \f$\rho\f$ + * @param P The pressure \f$P\f$ + */ +__attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure( + float density, float P) { + + const float density_inv = 1.f / density; + return sqrtf(hydro_gamma * P * density_inv); +} + +/** + * @brief Initialize the eos parameters + * + * @param e The #eos_paramters + * @param params The parsed parameters + */ +__attribute__((always_inline)) INLINE static void eos_init( + struct eos_parameters *e, const struct swift_params *params) {} + +/** + * @brief Print the equation of state + * + * @param e The #eos_parameters + */ +__attribute__((always_inline)) INLINE static void eos_print( + const struct eos_parameters *e) { + + message("Equation of state: Ideal gas."); + + 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 + */ +__attribute__((always_inline)) 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", "Ideal gas"); +} +#endif + +#endif /* SWIFT_TILLOTSON_IRON_EQUATION_OF_STATE_H */