Commit 52b7687f authored by Jacob Kegerreis 's avatar Jacob Kegerreis
Browse files

Move planetary EOS functions into their own headers

parent 12d10953
......@@ -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, aneos, sesame, planetary default: ideal-gas@:>@]
[equation of state @<:@ideal-gas, isothermal-gas, planetary default: ideal-gas@:>@]
)],
[with_eos="$withval"],
[with_eos="ideal-gas"]
......@@ -988,18 +988,6 @@ case "$with_eos" in
isothermal-gas)
AC_DEFINE([EOS_ISOTHERMAL_GAS], [1], [Isothermal gas equation of state])
;;
tillotson)
AC_DEFINE([EOS_TILLOTSON], [1], [Tillotson equations of state])
;;
hubbard-macfarlane)
AC_DEFINE([EOS_HUBBARD_MACFARLANE], [1], [Hubbard & MacFarlane (1980) equations of state])
;;
aneos)
AC_DEFINE([EOS_ANEOS], [1], [ANEOS equations of state])
;;
sesame)
AC_DEFINE([EOS_SESAME], [1], [SESAME equations of state])
;;
planetary)
AC_DEFINE([EOS_PLANETARY], [1], [All planetary equations of state])
;;
......
......@@ -34,14 +34,6 @@
#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)
#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_SESAME)
#include "./equation_of_state/sesame/equation_of_state.h"
#elif defined(EOS_PLANETARY)
#include "./equation_of_state/planetary/equation_of_state.h"
#else
......
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
* 2018 Jacob Kegerreis (jacob.kegerreis@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/planetary/aneos.h
*
* Contains the (M)ANEOS EOS functions for
* equation_of_state/planetary/equation_of_state.h
*
* WORK IN PROGRESS!
*
*/
/* 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"
#include "equation_of_state.h"
/* ------------------------------------------------------------------------- */
// ANEOS
struct ANEOS_params {
int mat_id;
int num_rho, num_u;
};
// Parameter values for each material (cgs units)
INLINE static void set_ANEOS_iron(struct ANEOS_params *mat, int mat_id) {
mat->mat_id = mat_id;
mat->num_rho = 100;
mat->num_u = 100;
}
INLINE static void set_MANEOS_forsterite(struct ANEOS_params *mat, int mat_id) {
mat->mat_id = mat_id;
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) {
}
// gas_pressure_from_internal_energy
INLINE static float ANEOS_pressure_from_internal_energy(float density, float u,
struct ANEOS_params *mat) {
float P;
/// Placeholder
P = mat->num_rho;
return P;
}
// gas_soundspeed_from_internal_energy
INLINE static float ANEOS_soundspeed_from_internal_energy(float density, float u,
struct ANEOS_params *mat) {
float c;
/// Placeholder
c = mat->num_rho;
return c;
}
// gas_soundspeed_from_pressure
INLINE static float ANEOS_soundspeed_from_pressure(float density, float P,
struct ANEOS_params *mat) {
float c;
/// Placeholder
c = mat->num_rho;
return c;
}
#endif /* SWIFT_ANEOS_EQUATION_OF_STATE_H */
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
* 2018 Jacob Kegerreis (jacob.kegerreis@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_HUBBARD_MACFARLANE_EQUATION_OF_STATE_H
#define SWIFT_HUBBARD_MACFARLANE_EQUATION_OF_STATE_H
/**
* @file equation_of_state/planetary/hm80.h
*
* Contains the Hubbard & MacFarlane (1980) Uranus/Neptune EOS functions for
* equation_of_state/planetary/equation_of_state.h
*
*/
/* 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"
#include "equation_of_state.h"
/* ------------------------------------------------------------------------- */
// Hubbard & MacFarlane (1980) table parameters
struct HM80_params {
int mat_id;
int num_rho, num_u;
float log_rho_min, log_rho_max, log_rho_step, inv_log_rho_step, log_u_min,
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 instead once finished testing...
#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"
// Parameter values for each material (cgs units)
INLINE static void set_HM80_HHe(struct HM80_params *mat, int mat_id) {
mat->mat_id = mat_id;
mat->num_rho = 100;
mat->num_u = 100;
mat->log_rho_min = -9.2103404;
mat->log_rho_max = 1.6094379;
mat->log_rho_step = 0.1092907;
mat->log_u_min = 9.2103404;
mat->log_u_max = 22.3327037;
mat->log_u_step = 0.1325491;
mat->bulk_mod = 0;
mat->inv_log_rho_step = 1.f / mat->log_rho_step;
mat->inv_log_u_step = 1.f / mat->log_u_step;
}
INLINE static void set_HM80_ice(struct HM80_params *mat, int mat_id) {
mat->mat_id = mat_id;
mat->num_rho = 200;
mat->num_u = 200;
mat->log_rho_min = -6.9077553;
mat->log_rho_max = 2.7080502;
mat->log_rho_step = 0.0483206;
mat->log_u_min = 6.9077553;
mat->log_u_max = 22.3327037;
mat->log_u_step = 0.0775123;
mat->bulk_mod = 2.0e10;
mat->inv_log_rho_step = 1.f / mat->log_rho_step;
mat->inv_log_u_step = 1.f / mat->log_u_step;
}
INLINE static void set_HM80_rock(struct HM80_params *mat, int mat_id) {
mat->mat_id = mat_id;
mat->num_rho = 100;
mat->num_u = 100;
mat->log_rho_min = -6.9077553;
mat->log_rho_max = 2.9957323;
mat->log_rho_step = 0.1000352;
mat->log_u_min = 9.2103404;
mat->log_u_max = 20.7232658;
mat->log_u_step = 0.1162922;
mat->bulk_mod = 3.49e11;
mat->inv_log_rho_step = 1.f / mat->log_rho_step;
mat->inv_log_u_step = 1.f / mat->log_u_step;
}
// Read the table from file
INLINE static void load_HM80_table(struct HM80_params *mat, char *table_file) {
// Allocate table memory
mat->table_P_rho_u = (float **) malloc(mat->num_rho*sizeof(float *));
for (int i=0; i<mat->num_rho; i++) {
mat->table_P_rho_u[i] = (float *) malloc(mat->num_u*sizeof(float));
}
// Load table contents from file
FILE *f = fopen(table_file, "r");
for (int i=0; i<mat->num_rho; i++) {
for (int j=0; j<mat->num_u; j++) {
fscanf(f, "%f", &mat->table_P_rho_u[i][j]);
}
}
fclose(f);
}
// Convert from cgs to internal units
#define Mbar_to_Ba 1e12 // Convert Megabar to Barye
#define J_kg_to_erg_g 1e4 // Convert J/kg to erg/g
INLINE static void convert_units_HM80(
struct HM80_params *mat, const struct unit_system* us) {
// Table densities in cgs
mat->log_rho_min -= log(units_cgs_conversion_factor(us, UNIT_CONV_DENSITY));
mat->log_rho_max -= log(units_cgs_conversion_factor(us, UNIT_CONV_DENSITY));
// Table energies in SI
mat->log_u_min += log(J_kg_to_erg_g /
units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
mat->log_u_max += log(J_kg_to_erg_g /
units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
// Table Pressures in Mbar
for (int i=0; i<mat->num_rho; i++) {
for (int j=0; j<mat->num_u; j++) {
mat->table_P_rho_u[i][j] *= Mbar_to_Ba /
units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
}
}
mat->bulk_mod /= units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
}
// gas_pressure_from_internal_energy
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);
const float log_u = log(u);
// 2D interpolation (linear in log(rho), log(u)) to find P(rho, u)
rho_idx = floor((log_rho - mat->log_rho_min) * mat->inv_log_rho_step);
u_idx = floor((log_u - mat->log_u_min) * mat->inv_log_u_step);
intp_rho = (log_rho - mat->log_rho_min - rho_idx*mat->log_rho_step) *
mat->inv_log_rho_step;
intp_u = (log_u - mat->log_u_min - u_idx*mat->log_u_step) *
mat->inv_log_u_step;
// Return zero pressure if below the table minimum/a
// Extrapolate the pressure for low densities
if (rho_idx < 0) { // Too-low rho
P = exp(log((1-intp_u)*mat->table_P_rho_u[0][u_idx]
+ intp_u*mat->table_P_rho_u[0][u_idx+1])
+ log_rho - mat->log_rho_min);
if (u_idx < 0) { // and too-low u
P = 0;
}
}
else if (u_idx < 0) { // Too-low u
P = 0;
}
// Return an edge value if above the table maximum/a
else if (rho_idx >= mat->num_rho-1) { // Too-high rho
if (u_idx >= mat->num_u-1) { // and too-high u
P = mat->table_P_rho_u[mat->num_rho-1][mat->num_u-1];
}
else {
P = mat->table_P_rho_u[mat->num_rho-1][u_idx];
}
}
else if (u_idx >= mat->num_u-1) { // Too-high u
P = mat->table_P_rho_u[rho_idx][mat->num_u-1];
}
// Normal interpolation within the table
else {
P = (1-intp_rho) * ((1-intp_u)*mat->table_P_rho_u[rho_idx][u_idx] +
intp_u*mat->table_P_rho_u[rho_idx][u_idx+1]) +
intp_rho * ((1-intp_u)*mat->table_P_rho_u[rho_idx+1][u_idx] +
intp_u*mat->table_P_rho_u[rho_idx+1][u_idx+1]);
}
return P;
}
// gas_soundspeed_from_internal_energy
INLINE static float HM80_soundspeed_from_internal_energy(float density, float u,
struct HM80_params *mat) {
float c, P;
// Bulk modulus
if (mat->bulk_mod != 0) {
c = sqrt(mat->bulk_mod / density);
}
// Ideal gas
else {
P = HM80_pressure_from_internal_energy(density, u, mat);
c = sqrt(5.f/3.f * P / density);
}
return c;
}
// gas_soundspeed_from_pressure
INLINE static float HM80_soundspeed_from_pressure(float density, float P,
struct HM80_params *mat) {
float c;
// Bulk modulus
if (mat->bulk_mod != 0) {
c = sqrt(mat->bulk_mod / density);
}
// Ideal gas
else {
c = sqrt(5.f/3.f * P / density);
}
return c;
}
#endif /* SWIFT_HUBBARD_MACFARLANE_EQUATION_OF_STATE_H */
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
* 2018 Jacob Kegerreis (jacob.kegerreis@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_SESAME_EQUATION_OF_STATE_H
#define SWIFT_SESAME_EQUATION_OF_STATE_H
/**
* @file equation_of_state/planetary/sesame.h
*
* Contains the SESAME EOS functions for
* equation_of_state/planetary/equation_of_state.h
*
* WORK IN PROGRESS!
*
*/
/* 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"
#include "equation_of_state.h"
/* ------------------------------------------------------------------------- */
// SESAME (WIP)
struct SESAME_params {
int mat_id;
};
// Parameter values for each material (cgs units)
INLINE static void set_SESAME_iron(struct SESAME_params *mat, int mat_id) {
mat->mat_id = mat_id;
}
// Convert from cgs to internal units
INLINE static void convert_units_SESAME(
struct SESAME_params *mat, const struct unit_system* us) {
}
// gas_pressure_from_internal_energy
INLINE static float SESAME_pressure_from_internal_energy(float density, float u,
struct SESAME_params *mat) {
float P;
/// Placeholder
P = mat->mat_id;
return P;
}
// gas_soundspeed_from_internal_energy
INLINE static float SESAME_soundspeed_from_internal_energy(float density, float u,
struct SESAME_params *mat) {
float c;
/// Placeholder
c = mat->mat_id;
return c;
}
// gas_soundspeed_from_pressure
INLINE static float SESAME_soundspeed_from_pressure(float density, float P,
struct SESAME_params *mat) {
float c;
/// Placeholder
c = mat->mat_id;
return c;
}
#endif /* SWIFT_SESAME_EQUATION_OF_STATE_H */
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
* 2018 Jacob Kegerreis (jacob.kegerreis@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_EQUATION_OF_STATE_H
#define SWIFT_TILLOTSON_EQUATION_OF_STATE_H
/**
* @file equation_of_state/planetary/tillotson.h
*
* Contains the Tillotson EOS functions for
* equation_of_state/planetary/equation_of_state.h
*
*/
/* 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"
#include "equation_of_state.h"
/* ------------------------------------------------------------------------- */
// Tillotson parameters
struct Til_params {
int mat_id;
float rho_0, a, b, A, B, E_0, E_iv, E_cv, alpha, beta, eta_min, P_min;
float c_TEMPORARY;
};
// Parameter values for each material (cgs units)
INLINE static void set_Til_iron(struct Til_params *mat, int mat_id) {
mat->mat_id = mat_id;
mat->rho_0 = 7.800;
mat->a = 0.5;
mat->b = 1.5;
mat->A = 1.28e12;
mat->B = 1.05e12;
mat->E_0 = 9.5e10;
mat->E_iv = 2.4e10;
mat->E_cv = 8.67e10;
mat->alpha = 5.0;
mat->beta = 5.0;
mat->eta_min = 0.0;
mat->P_min = 0.0;
mat->c_TEMPORARY = 9.4e-4;
}
INLINE static void set_Til_granite(struct Til_params *mat, int mat_id) {
mat->mat_id = mat_id;
mat->rho_0 = 2.680;
mat->a = 0.5;