diff --git a/configure.ac b/configure.ac index 864ab1acc759efcc6ac7247810e5d17dbca60260..219891509d63522d3e7338de0dd65de045fd915d 100644 --- a/configure.ac +++ b/configure.ac @@ -976,7 +976,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, tillotson, hubbard-macfarlane, planetary default: ideal-gas@:>@] + [equation of state @<:@ideal-gas, isothermal-gas, tillotson, hubbard-macfarlane, aneos, planetary default: ideal-gas@:>@] )], [with_eos="$withval"], [with_eos="ideal-gas"] @@ -989,13 +989,16 @@ case "$with_eos" in AC_DEFINE([EOS_ISOTHERMAL_GAS], [1], [Isothermal gas equation of state]) ;; tillotson) - AC_DEFINE([EOS_TILLOTSON], [1], [Tillotson equation of state]) + AC_DEFINE([EOS_TILLOTSON], [1], [Tillotson equations of state]) ;; hubbard-macfarlane) - AC_DEFINE([EOS_HUBBARD_MACFARLANE], [1], [Hubbard & MacFarlane (1980) equation of state]) + AC_DEFINE([EOS_HUBBARD_MACFARLANE], [1], [Hubbard & MacFarlane (1980) equations of state]) + ;; + aneos) + AC_DEFINE([EOS_ANEOS], [1], [ANEOS equations of state]) ;; planetary) - AC_DEFINE([EOS_PLANETARY], [1], [Planetary equations of state]) + AC_DEFINE([EOS_PLANETARY], [1], [All planetary equations 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 cf3a2a8e3f5fb87c7ed9d27014fbee8afbd0ea49..e6767d235df83002f5bcacd5c4717466a82fdf94 100644 --- a/src/equation_of_state.h +++ b/src/equation_of_state.h @@ -38,6 +38,8 @@ #include "./equation_of_state/tillotson/equation_of_state.h" #elif defined(EOS_HUBBARD_MACFARLANE) #include "./equation_of_state/hubbard_macfarlane/equation_of_state.h" +#elif defined(EOS_ANEOS) +#include "./equation_of_state/aneos/equation_of_state.h" #elif defined(EOS_PLANETARY) #include "./equation_of_state/planetary/equation_of_state.h" #else diff --git a/src/equation_of_state/aneos/equation_of_state.h b/src/equation_of_state/aneos/equation_of_state.h new file mode 100644 index 0000000000000000000000000000000000000000..4d2472f7a7427db7c32a141255108dd698f3c61a --- /dev/null +++ b/src/equation_of_state/aneos/equation_of_state.h @@ -0,0 +1,246 @@ +/******************************************************************************* + * 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_ANEOS_EQUATION_OF_STATE_H +#define SWIFT_ANEOS_EQUATION_OF_STATE_H + +/** + * @file equation_of_state/aneos/equation_of_state.h + * + * WIP PLACEHOLDER ONLY! + */ + +/* Some standard headers. */ +#include <math.h> + +/* Local headers. */ +#include "adiabatic_index.h" +#include "common_io.h" +#include "inline.h" +#include "units.h" +#include "physical_constants.h" + +extern struct eos_parameters eos; +/* ------------------------------------------------------------------------- */ + +// ANEOS parameters +struct ANEOS_params { + int mat_id; + int num_rho, num_u; +}; + +struct eos_parameters { + struct ANEOS_params ANEOS_iron, MANEOS_forsterite; +}; + +// Material identifier flags (material_ID = type_ID * type_factor + unit_ID) +#define type_factor 10 +enum type_id { + type_ANEOS = 3 +}; +enum material_id { + // ANEOS + ANEOS_iron = type_ANEOS*type_factor, + MANEOS_forsterite = type_ANEOS*type_factor + 1 +}; + +// Parameter values for each material (cgs units) +INLINE static void set_ANEOS_iron(struct ANEOS_params *mat) { + mat->mat_id = ANEOS_iron; + mat->num_rho = 100; + mat->num_u = 100; +} +INLINE static void set_ANEOS_forsterite(struct ANEOS_params *mat) { + mat->mat_id = MANEOS_forsterite; + mat->num_rho = 100; + mat->num_u = 100; +} + +// Convert from cgs to internal units +INLINE static void convert_units_ANEOS( + struct ANEOS_params *mat, const struct unit_system* us) { + +} + +/** + * @brief Returns the internal energy given density and entropy + * + * NOT IMPLEMENTED! + * + * @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, int mat_id) { + + return 0; +} + +/** + * @brief Returns the pressure given density and entropy + * + * NOT IMPLEMENTED! + * + * @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, int mat_id) { + + return 0; +} + +/** + * @brief Returns the entropy given density and pressure. + * + * NOT IMPLEMENTED! + * + * @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, int mat_id) { + + return 0; +} + +/** + * @brief Returns the sound speed given density and entropy + * + * NOT IMPLEMENTED! + * + * @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, int mat_id) { + + return 0; +} + +/** + * @brief Returns the entropy given density and internal energy + * + * NOT IMPLEMENTED! + * + * @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, int mat_id) { + + return 0; +} + +/** + * @brief Returns the pressure given density and internal energy + * + * @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, int mat_id) { + + return 0; +} + +/** + * @brief Returns the internal energy given density and pressure. + * + * NOT IMPLEMENTED! + * + * @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, int mat_id) { + + return 0; +} + +/** + * @brief Returns the sound speed given density and internal energy + * + * @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, int mat_id) { + + return 0; +} + +/** + * @brief Returns the sound speed given density and pressure + * + * @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, int mat_id) { + + return 0; +} + +/** + * @brief Initialize the eos parameters + * + * @param e The #eos_parameters + * @param params The parsed parameters + */ +__attribute__((always_inline)) INLINE static void eos_init( + struct eos_parameters *e, const struct phys_const *phys_const, + const struct unit_system *us, const struct swift_params *params) { + + // Set the parameters and load tables etc. for each material + set_ANEOS_iron(&e->ANEOS_iron); + set_ANEOS_forsterite(&e->MANEOS_forsterite); + + // Convert to internal units + convert_units_ANEOS(&e->ANEOS_iron, us); + convert_units_ANEOS(&e->MANEOS_forsterite, us); +} + +/** + * @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: ANEOS."); +} + +#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_s(h_grpsph, "Equation of state", "ANEOS"); +} +#endif + +#endif /* SWIFT_ANEOS_EQUATION_OF_STATE_H */ diff --git a/src/equation_of_state/planetary/equation_of_state.h b/src/equation_of_state/planetary/equation_of_state.h index 4f9199dc4737b3aa271eb61236d593bc9bff1584..5d3378ae004fbc4b141044180a68b65c3e380f04 100644 --- a/src/equation_of_state/planetary/equation_of_state.h +++ b/src/equation_of_state/planetary/equation_of_state.h @@ -61,7 +61,8 @@ struct HM80_params { log_u_max, log_u_step, inv_log_u_step, bulk_mod; float **table_P_rho_u; }; -// Table file names /// to be read in from the parameter file once tested... +// Table file names +/// to be read in from the parameter file instead once tested... #define HM80_HHe_table_file "/gpfs/data/dc-kege1/gihr_data/P_rho_u_HHe.txt" #define HM80_ice_table_file "/gpfs/data/dc-kege1/gihr_data/P_rho_u_ice.txt" #define HM80_rock_table_file "/gpfs/data/dc-kege1/gihr_data/P_rho_u_roc.txt" @@ -218,6 +219,18 @@ INLINE static void load_HM80_table(struct HM80_params *mat, char *table_file) { fclose(f); } +// ANEOS +INLINE static void set_ANEOS_iron(struct ANEOS_params *mat) { + mat->mat_id = ANEOS_iron; + mat->num_rho = 100; + mat->num_u = 100; +} +INLINE static void set_ANEOS_forsterite(struct ANEOS_params *mat) { + mat->mat_id = ANEOS_forsterite; + mat->num_rho = 100; + mat->num_u = 100; +} + // Convert from cgs to internal units // Tillotson INLINE static void convert_units_Til( @@ -255,6 +268,11 @@ INLINE static void convert_units_HM80( mat->bulk_mod /= units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE); } +// ANEOS +INLINE static void convert_units_ANEOS( + struct ANEOS_params *mat, const struct unit_system* us) { + +} /** * @brief Returns the internal energy given density and entropy @@ -376,6 +394,11 @@ INLINE static float Til_pressure_from_internal_energy(float density, float u, INLINE static float HM80_pressure_from_internal_energy(float density, float u, struct HM80_params *mat) { float P; + + if (u <= 0) { + return 0; + } + int rho_idx, u_idx; float intp_rho, intp_u; const float log_rho = log(density); @@ -498,11 +521,6 @@ gas_pressure_from_internal_energy(float density, float u, int mat_id) { mat_HM80 = &eos.HM80_HHe; // Ignored, just here to keep the compiler happy }; - if (u <= 0) { - P = 0; - break; - } - P = HM80_pressure_from_internal_energy(density, u, mat_HM80); break; @@ -906,9 +924,11 @@ __attribute__((always_inline)) INLINE static void eos_init( load_HM80_table(&e->HM80_ice, HM80_ice_table_file); load_HM80_table(&e->HM80_rock, HM80_rock_table_file); - // ANEOS (WIP) + // ANEOS + set_ANEOS_iron(&e->ANEOS_iron); + set_ANEOS_forsterite(&e->ANEOS_forsterite); - // Convert from cgs units to internal units + // Convert to internal units // Tillotson convert_units_Til(&e->Til_iron, us); convert_units_Til(&e->Til_granite, us); @@ -919,7 +939,9 @@ __attribute__((always_inline)) INLINE static void eos_init( convert_units_HM80(&e->HM80_ice, us); convert_units_HM80(&e->HM80_rock, us); - // ANEOS (WIP) + // ANEOS + convert_units_ANEOS(&e->ANEOS_iron, us); + convert_units_ANEOS(&e->ANEOS_forsterite, us); } /** diff --git a/src/equation_of_state/tillotson/equation_of_state.h b/src/equation_of_state/tillotson/equation_of_state.h index 06d9fcb3dcbab9b48fc0deac525a2f33ec007588..c82fead13bf1c0bb61e24f2f1ee929ab43fc13af 100644 --- a/src/equation_of_state/tillotson/equation_of_state.h +++ b/src/equation_of_state/tillotson/equation_of_state.h @@ -22,7 +22,7 @@ /** * @file equation_of_state/tillotson/equation_of_state.h * - * Only P(rho, u), c_s(rho, u), and c_s(rho, P) are implemented for now! + * Only P(rho, u), c(rho, u), and c(rho, P) are implemented for now! * So, must be used with the MinimalMultiMat SPH formulation. */