tboudreaux/SERiF
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

feat(eos): added helmholtz eos as module

Browse Source 
Aaron Dotter implimented a C++ version of Frank Timmes' fortran code helmholtz.f90. I have taken that and refactored it to work in the 4DSEE code style. This has mostly involved some light moving of stuff around. The biggest change is removing all globals, and reorienting memory to be heap allocated and contiguous. This is because there was too much memory being stack allocated.
This commit is contained in:
Emily Boudreaux
2025-03-05 16:59:04 -05:00
parent e43caf3027
commit 154004c8ca
4 changed files with 1205 additions and 865 deletions
Download Patch File Download Diff File Expand all files Collapse all files

23
src/eos/meson.build Normal file
View File

@@ -0,0 +1,23 @@
# Define the library
eos_sources = files(
'private/helm.cpp',
)
eos_headers = files(
'public/helm.h'
)
# Define the libconst library so it can be linked against by other parts of the build system
libeos = static_library('eos',
eos_sources,
include_directories: include_directories('public'),
cpp_args: ['-fvisibility=default'],
dependencies: [const_dep, quill_dep, probe_dep, config_dep, mfem_dep],
install : true)
eos_dep = declare_dependency(
include_directories: include_directories('public'),
link_with: libeos,
)
# Make headers accessible
install_headers(eos_headers, subdir : '4DSSE/eos')

823
src/eos/private/helm.cpp Normal file
View File

@@ -0,0 +1,823 @@
// this file contains a C++ conversion of Frank Timmes' fortran code
// helmholtz.f90, which implements the Helmholtz Free Energy EOS described
// by Timmes & Swesty (2000) doi:10.1086/313304 -- Aaron Dotter 2025
#include <iostream>
#include <fstream>
#include <sstream>
#include <cmath>
#include <string>
#include <stdexcept>
#include <array>
#include <numbers>
#include "helm.h"
#include "const.h"
#include "probe.h"
#include "config.h"
#include "quill/LogMacros.h"
using namespace std;
double** heap_allocate_contiguous_2D_memory(int rows, int cols) {
double **array = new double*[rows];
array[0] = new double[rows * cols];
for (int i = 1; i < rows; i++) {
array[i] = array[0] + i * cols;
}
return array;
}
void heap_deallocate_contiguous_2D_memory(double **array) {
delete[] array[0];
delete[] array;
}
// interpolating polynomila function definitions
namespace helmholtz {
double psi0(double z) { return z*z*z * (z * (-6.0 * z + 15.0) - 10.0) + 1.0; }
double dpsi0(double z) { return z*z * ( z * (-30.0 * z + 60.0) - 30.0); }
double ddpsi0(double z) { return z * ( z * (-120.0 * z + 180.0) - 60.0); }
double psi1(double z) { return z * ( z*z * ( z * (-3.0*z + 8.0) - 6.0) + 1.0); }
double dpsi1(double z) { return z*z * ( z * (-15.0*z + 32.0) - 18.0) +1.0; }
double ddpsi1(double z) { return z * (z * (-60.0*z + 96.0) - 36.0); }
double psi2(double z) { return 0.5 * z*z * ( z* ( z * (-z + 3.0) - 3.0) + 1.0); }
double dpsi2(double z) { return 0.5 * z * ( z * (z * (-5.0*z + 12.0) - 9.0) + 2.0); }
double ddpsi2(double z) { return 0.5*(z*( z * (-20.*z + 36.) - 18.) + 2.); }
double xpsi0(double z) { return z * z * (2.0 * z - 3.0) + 1.0; }
double xdpsi0(double z) { return z * (6.0*z - 6.0); }
double xpsi1(double z) { return z * ( z * (z - 2.0) + 1.0); }
double xdpsi1(double z) { return z * (3.0 * z - 4.0) + 1.0; }
double h3(std::array<double, 36> fi, double w0t, double w1t, double w0mt, double w1mt,
double w0d, double w1d, double w0md, double w1md) {
return fi[0] * w0d * w0t + fi[1] * w0md * w0t
+ fi[2] * w0d * w0mt + fi[3] * w0md * w0mt
+ fi[4] * w0d * w1t + fi[5] * w0md * w1t
+ fi[6] * w0d * w1mt + fi[7] * w0md * w1mt
+ fi[8] * w1d * w0t + fi[9] * w1md * w0t
+ fi[10] * w1d * w0mt + fi[11] * w1md * w0mt
+ fi[12] * w1d * w1t + fi[13] * w1md * w1t
+ fi[14] * w1d * w1mt + fi[15] * w1md * w1mt;
}
double h5(std::array<double, 36> fi, double w0t, double w1t, double w2t,double w0mt,
double w1mt, double w2mt, double w0d, double w1d, double w2d,
double w0md, double w1md, double w2md) {
return fi[0] * w0d * w0t + fi[1] * w0md * w0t
+ fi[2] * w0d * w0mt + fi[3] * w0md * w0mt
+ fi[4] * w0d * w1t + fi[5] * w0md * w1t
+ fi[6] * w0d * w1mt + fi[7] * w0md * w1mt
+ fi[8] * w0d * w2t + fi[ 9] * w0md * w2t
+ fi[10] * w0d * w2mt + fi[11] * w0md * w2mt
+ fi[12] * w1d * w0t + fi[13] * w1md * w0t
+ fi[14] * w1d * w0mt + fi[15] * w1md * w0mt
+ fi[16] * w2d * w0t + fi[17] * w2md * w0t
+ fi[18] * w2d * w0mt + fi[19] * w2md * w0mt
+ fi[20] * w1d * w1t + fi[21] * w1md * w1t
+ fi[22] * w1d * w1mt + fi[23] * w1md * w1mt
+ fi[24] * w2d * w1t + fi[25] * w2md * w1t
+ fi[26] * w2d * w1mt + fi[27] * w2md * w1mt
+ fi[28] * w1d * w2t + fi[29] * w1md * w2t
+ fi[30] * w1d * w2mt + fi[31] * w1md * w2mt
+ fi[32] * w2d * w2t + fi[33] * w2md * w2t
+ fi[34] * w2d * w2mt + fi[35] * w2md * w2mt;
}
// this function reads in the HELM table and stores in the above arrays
HELMTable read_helm_table(const std::string filename) {
Config& config = Config::getInstance();
std::string logFile = config.get<std::string>("EOS:Helm:LogFile", "log");
Probe::LogManager& logManager = Probe::LogManager::getInstance();
quill::Logger* logger = logManager.getLogger(logFile);
LOG_INFO(logger, "read_helm_table : Reading HELM table from file {}", filename);
HELMTable table;
string data;
int i, j;
//set T and Rho (d) arrays
for (j=0; j<table.jmax; j++) { table.t[j] = exp10(tlo + tstp*j); }
for (i=0; i<table.imax; i++) { table.d[i] = exp10(dlo + dstp*i); }
ifstream helm_table(filename);
if (!helm_table) {
LOG_ERROR(logger, "read_helm_table : Error opening file {}", filename);
throw std::runtime_error("Error opening file " + filename);
}
//read the Helmholtz free energy and its derivatives
for (j=0; j<table.jmax; j++) {
for (i=0; i<table.imax; i++){
getline(helm_table, data);
stringstream id(data);
id >> table.f[i][j];
id >> table.fd[i][j];
id >> table.ft[i][j];
id >> table.fdd[i][j];
id >> table.ftt[i][j];
id >> table.fdt[i][j];
id >> table.fddt[i][j];
id >> table.fdtt[i][j];
id >> table.fddtt[i][j];
}
}
//read the pressure derivative with density
for (j=0; j<table.jmax; j++) {
for (i=0; i<table.imax; i++){
getline(helm_table, data);
stringstream id(data);
id >> table.dpdf[i][j];
id >> table.dpdfd[i][j];
id >> table.dpdft[i][j];
id >> table.dpdfdt[i][j];
}
}
//read the electron chemical potential
for (j=0; j<table.jmax; j++) {
for (i=0; i<table.imax; i++){
getline(helm_table, data);
stringstream id(data);
id >> table.ef[i][j];
id >> table.efd[i][j];
id >> table.eft[i][j];
id >> table.efdt[i][j];
}
}
//read the number density
for (j=0; j<table.jmax; j++) {
for (i=0; i<table.imax; i++){
getline(helm_table, data);
stringstream id(data);
id >> table.xf[i][j];
id >> table.xfd[i][j];
id >> table.xft[i][j];
id >> table.xfdt[i][j];
}
}
helm_table.close(); //done reading
// construct the temperature and density deltas and their inverses
for (j=0; j<table.jmax; j++) {
double dth = table.t[j+1] - table.t[j];
double dt2 = dth * dth;
double dti = 1.0/dth;
double dt2i = 1.0/dt2;
table.dt_sav[j] = dth;
table.dt2_sav[j] = dt2;
table.dti_sav[j] = dti;
table.dt2i_sav[j] = dt2i;
}
for (i=0; i<table.imax; i++) {
double dd = table.d[i+1] - table.d[i];
double dd2 = dd * dd;
double ddi = 1.0/dd;
table.dd_sav[i] = dd;
table.dd2_sav[i] = dd2;
table.ddi_sav[i] = ddi;
}
table.loaded = true;
return table;
}
/***
this function evaluates the EOS components:
ion, radiation, electron-positron and Coulomb interaction
and returns the calculated quantities in the input
***/
EOS get_helm_EOS(EOSInput &q, const HELMTable &table) {
Config& config = Config::getInstance();
std::string logFile = config.get<std::string>("EOS:Helm:LogFile", "log");
Probe::LogManager& logManager = Probe::LogManager::getInstance();
quill::Logger* logger = logManager.getLogger(logFile);
Constants& constants = Constants::getInstance();
const double pi = std::numbers::pi;
const double amu = constants.get("u").value;
const double h = constants.get("h").value;
const double avo = constants.get("N_a").value;
const double kerg = constants.get("kB").value;
const double clight = constants.get("c").value;
const double kergavo = kerg*avo;
const double clight2 = clight * clight;
const double ssol = 5.6704e-5;
const double asol = 4 * ssol / clight;
const double asoli3 = asol / 3;
const double sioncon = (2.0 * pi * amu * kerg)/(h*h);
const double qe = 4.8032042712e-10;
const double esqu = qe * qe;
const double third = 1./3.;
std::array<double, 36> fi;
double T, rho, s, free, abar, zbar, ytot1, ye, ptot;
double prad, srad, erad, etot, stot;
double dpraddd, dpraddt, dpradda, dpraddz;
double deraddd, deraddt, deradda, deraddz;
double dsraddd, dsraddt, dsradda, dsraddz;
double dpiondd, dpiondt, dpionda, dpiondz;
double deiondd, deiondt, deionda, deiondz;
double dsiondd, dsiondt, dsionda, dsiondz;
double dpcouldd, dpcouldt, dpcoulda, dpcouldz;
double decouldd, decouldt, decoulda, decouldz;
double dscouldd, dscouldt, dscoulda, dscouldz;
double dpgasdd, dpgasdt, dpgasda, dpgasdz;
double degasdd, degasdt, degasda, degasdz;
double dsgasdd, dsgasdt, dsgasda, dsgasdz;
double dpepdd, dpepdt, dpepda, dpepdz;
double dsepdd, dsepdt, dsepda, dsepdz;
double deepdd, deepdt, deepda, deepdz;
double dxnidd, dxnida;
double dpresdd, dpresdt, dpresda, dpresdz;
double denerdd, denerdt, denerda, denerdz;
double dentrdd, dentrdt, dentrda, dentrdz;
double xni, sion, pion, eion, x, y, z;
double pgas, egas, sgas, pele, sele, eele, pcoul, scoul, ecoul;
int jat, iat;
T = q.T; rho = q.rho; abar = q.abar; zbar = q.zbar;
ytot1 = 1/abar;
ye = zbar * ytot1;
double rhoi = 1/rho;
double Ti = 1/T;
double kT= kerg*T;
double kTinv = 1/kT;
/*** radiation ***/
prad = asoli3 * pow(T,4);
dpraddd = 0.;
dpraddt = 4 * prad * Ti;
dpradda = 0.;
dpraddz = 0.;
erad = 3 * prad * rhoi;
deraddd = -erad * rhoi;
deraddt = 3 * dpraddt * rhoi;
deradda = 0.;
deraddz = 0.;
srad = (prad*rhoi + erad)*Ti;
dsraddd = (dpraddd*rhoi - prad*rhoi*rhoi + deraddd)*Ti;
dsraddt = (dpraddt*rhoi + deraddt - srad)*Ti;
dsradda = 0.;
dsraddz = 0.;
/*** ions ***/
xni = avo * ytot1 * rho;
dxnidd = avo * ytot1;
dxnida = -xni * ytot1;
pion = xni * kT;
dpiondd = dxnidd * kT;
dpiondt = xni * kerg;
dpionda = dxnida * kT;
dpiondz = 0.;
eion = 1.5 * pion * rhoi;
deiondd = (1.5 * dpiondd - eion)*rhoi;
deiondt = 1.5 * dpiondt*rhoi;
deionda = 1.5 * dpionda*rhoi;
deiondz = 0.;
x = abar * abar * sqrt(abar) * rhoi/avo;
s = sioncon * T;
z = x * s * sqrt(s);
y = log(z);
sion = (pion * rhoi + eion) * Ti + kergavo * ytot1 * y;
dsiondd = (dpiondd*rhoi - pion*rhoi*rhoi + deiondd)*Ti
- kergavo * rhoi * ytot1;
dsiondt = (dpiondt*rhoi + deiondt)*Ti - (pion*rhoi + eion) * Ti*Ti
+ 1.5 * kergavo * Ti*ytot1;
x = avo*kerg/abar;
dsionda = (dpionda*rhoi + deionda)*Ti + kergavo*ytot1*ytot1* (2.5 - y);
dsiondz = 0.;
// electron-positron
// double xnem = xni * zbar;
double din = ye*rho;
// hash the table location in t, rho:
jat = (log10(T) - tlo)*tstpi;
jat = max(0, min(jat, table.jmax-1)); // bounds
iat = max(0, min(iat, table.imax-1)); // bounds
iat = (log10(din)-dlo)*dstpi;
//using the indices, access the table:
fi[0] = table.f[iat][jat];
fi[1] = table.f[iat+1][jat];
fi[2] = table.f[iat][jat+1];
fi[3] = table.f[iat+1][jat+1];
fi[4] = table.ft[iat][jat];
fi[5] = table.ft[iat+1][jat];
fi[6] = table.ft[iat][jat+1];
fi[7] = table.ft[iat+1][jat+1];
fi[8] = table.ftt[iat][jat];
fi[9] = table.ftt[iat+1][jat];
fi[10] = table.ftt[iat][jat+1];
fi[11] = table.ftt[iat+1][jat+1];
fi[12] = table.fd[iat][jat];
fi[13] = table.fd[iat+1][jat];
fi[14] = table.fd[iat][jat+1];
fi[15] = table.fd[iat+1][jat+1];
fi[16] = table.fdd[iat][jat];
fi[17] = table.fdd[iat+1][jat];
fi[18] = table.fdd[iat][jat+1];
fi[19] = table.fdd[iat+1][jat+1];
fi[20] = table.fdt[iat][jat];
fi[21] = table.fdt[iat+1][jat];
fi[22] = table.fdt[iat][jat+1];
fi[23] = table.fdt[iat+1][jat+1];
fi[24] = table.fddt[iat][jat];
fi[25] = table.fddt[iat+1][jat];
fi[26] = table.fddt[iat][jat+1];
fi[27] = table.fddt[iat+1][jat+1];
fi[28] = table.fdtt[iat][jat];
fi[29] = table.fdtt[iat+1][jat];
fi[30] = table.fdtt[iat][jat+1];
fi[31] = table.fdtt[iat+1][jat+1];
fi[32] = table.fddtt[iat][jat];
fi[33] = table.fddtt[iat+1][jat];
fi[34] = table.fddtt[iat][jat+1];
fi[35] = table.fddtt[iat+1][jat+1];
double xt = (T - table.t[jat]) * table.dti_sav[jat];
double xd = (din - table.d[iat]) * table.ddi_sav[iat];
double mxt = 1-xt;
double mxd = 1-xd;
//density and temperature basis functions
double si0t = psi0(xt);
double si1t = psi1(xt)*table.dt_sav[jat];
double si2t = psi2(xt)*table.dt2_sav[jat];
double si0mt = psi0(mxt);
double si1mt = -psi1(mxt)*table.dt_sav[jat];
double si2mt = psi2(mxt)*table.dt2_sav[jat];
double si0d = psi0(xd);
double si1d = psi1(xd)*table.dd_sav[iat];
double si2d = psi2(xd)*table.dd2_sav[iat];
double si0md = psi0(mxd);
double si1md = -psi1(mxd)*table.dd_sav[iat];
double si2md = psi2(mxd)*table.dd2_sav[iat];
// derivatives of the weight functions
double dsi0t = dpsi0(xt)*table.dti_sav[jat];
double dsi1t = dpsi1(xt);
double dsi2t = dpsi2(xt)*table.dt_sav[jat];
double dsi0mt = -dpsi0(mxt)*table.dti_sav[jat];
double dsi1mt = dpsi1(mxt);
double dsi2mt = -dpsi2(mxt)*table.dt_sav[jat];
double dsi0d = dpsi0(xd)*table.ddi_sav[iat];
double dsi1d = dpsi1(xd);
double dsi2d = dpsi2(xd)*table.dd_sav[iat];
double dsi0md = -dpsi0(mxd)*table.ddi_sav[iat];
double dsi1md = dpsi1(mxd);
double dsi2md = -dpsi2(mxd)*table.dd_sav[iat];
// second derivatives of the weight functions
double ddsi0t = ddpsi0(xt)*table.dt2i_sav[jat];
double ddsi1t = ddpsi1(xt)*table.dti_sav[jat];
double ddsi2t = ddpsi2(xt);
double ddsi0mt = ddpsi0(mxt)*table.dt2i_sav[jat];
double ddsi1mt = -ddpsi1(mxt)*table.dti_sav[jat];
double ddsi2mt = ddpsi2(mxt);
// Refactor these to use LOG_DEBUG(logger, message)
LOG_DEBUG(logger, " xt, mxt = {} {}", xt, mxt);
LOG_DEBUG(logger, " dti_sav[jat] = {}", table.dti_sav[jat]);
LOG_DEBUG(logger, " dt2i_sav[jat] = {}", table.dt2i_sav[jat]);
LOG_DEBUG(logger, " ddpsi0(xt) = {}", ddpsi0(xt));
LOG_DEBUG(logger, " ddpsi1(xt) = {}", ddpsi1(xt));
LOG_DEBUG(logger, " ddpsi2(xt) = {}", ddpsi2(xt));
LOG_DEBUG(logger, " ddpsi0(mxt) = {}", ddpsi0(mxt));
LOG_DEBUG(logger, " ddpsi1(mxt) = {}", ddpsi1(mxt));
LOG_DEBUG(logger, " ddpsi2(mxt) = {}", ddpsi2(mxt));
LOG_DEBUG(logger, " ddsi0t = {}", ddsi0t);
LOG_DEBUG(logger, " ddsi1t = {}", ddsi1t);
LOG_DEBUG(logger, " ddsi2t = {}", ddsi2t);
LOG_DEBUG(logger, " ddsi0mt = {}", ddsi0mt);
LOG_DEBUG(logger, " ddsi1mt = {}", ddsi1mt);
LOG_DEBUG(logger, " ddsi2mt = {}", ddsi2mt);
ddsi0t = 0.0; ddsi1t = 0.0; ddsi2t = 1.0;
ddsi0mt = 0.0; ddsi1mt = 0.0; ddsi2mt = 0.0;
// free energy
free = h5(fi,si0t, si1t, si2t, si0mt, si1mt, si2mt,
si0d, si1d, si2d, si0md, si1md, si2md);
// derivative with respect to density
double df_d = h5(fi, si0t, si1t, si2t, si0mt, si1mt, si2mt,
dsi0d, dsi1d, dsi2d, dsi0md, dsi1md, dsi2md);
// derivative with respect to temperature
double df_t = h5(fi, dsi0t, dsi1t, dsi2t, dsi0mt, dsi1mt, dsi2mt,
si0d, si1d, si2d, si0md, si1md, si2md);
// second derivative with respect to temperature
double df_tt = h5(fi, ddsi0t, ddsi1t, ddsi2t, ddsi0mt, ddsi1mt, ddsi2mt,
si0d, si1d, si2d, si0md, si1md, si2md);
LOG_DEBUG(logger, "df_tt = {}", df_tt);
// derivative with respect to temperature and density
double df_dt = h5(fi, dsi0t, dsi1t, dsi2t, dsi0mt, dsi1mt, dsi2mt,
dsi0d, dsi1d, dsi2d, dsi0md, dsi1md, dsi2md);
// pressure derivatives
si0t = xpsi0(xt);
si1t = xpsi1(xt)*table.dt_sav[jat];
si0mt = xpsi0(mxt);
si1mt = -xpsi1(mxt)*table.dt_sav[jat];
si0d = xpsi0(xd);
si1d = xpsi1(xd)*table.dd_sav[iat];
si0md = xpsi0(mxd);
si1md = -xpsi1(mxd)*table.dd_sav[iat];
// derivatives of weight functions
dsi0t = xdpsi0(xt)*table.dti_sav[jat];
dsi1t = xdpsi1(xt);
dsi0mt = -xdpsi0(mxt)*table.dti_sav[jat];
dsi1mt = xdpsi1(mxt);
dsi0d = xdpsi0(xd)*table.ddi_sav[iat];
dsi1d = xdpsi1(xd);
dsi0md = -xdpsi0(mxd)*table.ddi_sav[iat];
dsi1md = xdpsi1(mxd);
// pressure derivative
fi[0] = table.dpdf[iat][jat];
fi[1] = table.dpdf[iat+1][jat];
fi[2] = table.dpdf[iat][jat+1];
fi[3] = table.dpdf[iat+1][jat+1];
fi[4] = table.dpdft[iat][jat];
fi[5] = table.dpdft[iat+1][jat];
fi[6] = table.dpdft[iat][jat+1];
fi[7] = table.dpdft[iat+1][jat+1];
fi[8] = table.dpdfd[iat][jat];
fi[9] = table.dpdfd[iat+1][jat];
fi[10] = table.dpdfd[iat][jat+1];
fi[11] = table.dpdfd[iat+1][jat+1];
fi[12] = table.dpdfdt[iat][jat];
fi[13] = table.dpdfdt[iat+1][jat];
fi[14] = table.dpdfdt[iat][jat+1];
fi[15] = table.dpdfdt[iat+1][jat+1];
// pressure derivative with density
dpepdd = h3(fi, si0t, si1t, si0mt, si1mt, si0d, si1d, si0md, si1md);
dpepdd = max(ye * dpepdd,1e-30);
// look in the electron chemical potential table only once
fi[0] = table.ef[iat][jat];
fi[1] = table.ef[iat+1][jat];
fi[2] = table.ef[iat][jat+1];
fi[3] = table.ef[iat+1][jat+1];
fi[4] = table.eft[iat][jat];
fi[5] = table.eft[iat+1][jat];
fi[6] = table.eft[iat][jat+1];
fi[7] = table.eft[iat+1][jat+1];
fi[8] = table.efd[iat][jat];
fi[9] = table.efd[iat+1][jat];
fi[10] = table.efd[iat][jat+1];
fi[11] = table.efd[iat+1][jat+1];
fi[12] = table.efdt[iat][jat];
fi[13] = table.efdt[iat+1][jat];
fi[14] = table.efdt[iat][jat+1];
fi[15] = table.efdt[iat+1][jat+1];
// electron chemical potential etaele
double etaele = h3(fi, si0t, si1t, si0mt, si1mt, si0d, si1d, si0md, si1md);
// derivative with respect to density
x = h3(fi, si0t, si1t, si0mt, si1mt, dsi0d, dsi1d, dsi0md, dsi1md);
// double detadd = ye * x;
// derivative with respect to temperature
// double detadt = h3(fi, dsi0t, dsi1t, dsi0mt, dsi1mt, si0d,
// si1d, si0md, si1md);
// derivative with respect to abar and zbar
// double detada = -x * din * ytot1;
// double detadz = x * rho * ytot1;
// look in the number density table only once
fi[0] = table.xf[iat][jat];
fi[1] = table.xf[iat+1][jat];
fi[2] = table.xf[iat][jat+1];
fi[3] = table.xf[iat+1][jat+1];
fi[4] = table.xft[iat][jat];
fi[5] = table.xft[iat+1][jat];
fi[6] = table.xft[iat][jat+1];
fi[7] = table.xft[iat+1][jat+1];
fi[8] = table.xfd[iat][jat];
fi[9] = table.xfd[iat+1][jat];
fi[10] = table.xfd[iat][jat+1];
fi[11] = table.xfd[iat+1][jat+1];
fi[12] = table.xfdt[iat][jat];
fi[13] = table.xfdt[iat+1][jat];
fi[14] = table.xfdt[iat][jat+1];
fi[15] = table.xfdt[iat+1][jat+1];
// electron chemical potential etaele
double xnefer = h3(fi, si0t, si1t, si0mt, si1mt, si0d, si1d, si0md, si1md);
// derivative with respect to density
x = h3(fi, si0t, si1t, si0mt, si1mt, dsi0d, dsi1d, dsi0md, dsi1md);
x = max(x,1e-30);
// double dxnedd = ye * x;
// derivative with respect to temperature
// double dxnedt = h3(fi, dsi0t, dsi1t, dsi0mt, dsi1mt, si0d, si1d, si0md, si1md);
// derivative with respect to abar and zbar
// double dxneda = -x * din * ytot1;
// double dxnedz = x * rho * ytot1;
// the desired electron-positron thermodynamic quantities
x = din * din;
pele = x * df_d;
dpepdt = x * df_dt;
//dpepdd is set above
s = dpepdd/ye - 2 * din * df_d;
dpepda = -ytot1 * (2 * pele + s * din);
dpepdz = rho * ytot1 * (2 * din * df_d + s);
x = ye * ye;
sele = -df_t * ye;
dsepdt = -df_tt * ye;
dsepdd = -df_dt * x;
dsepda = ytot1 * (ye * df_dt * din - sele);
dsepdz = -ytot1 * (ye * df_dt * rho + df_t);
eele = ye*free + T * sele;
deepdt = T * dsepdt;
deepdd = x * df_d + T * dsepdd;
deepda = -ye * ytot1 * (free + df_d * din) + T * dsepda;
deepdz = ytot1 * (free + ye * df_d * rho) + T * dsepdz;
// coulomb
const double a1 = -0.898004;
const double b1 = 0.96786;
const double c1 = 0.220703;
const double d1 = -0.86097;
const double e1 = 2.5269;
const double a2 = 0.29561;
const double b2 = 1.9885;
const double c2 = 0.288675;
z = 4 * pi * third;
s = z * xni;
double dsdd = z * dxnidd;
double dsda = z * dxnida;
double lami = pow(s,-1./3.);
double inv_lami = 1.0/lami;
z = -lami/3.;
double lamidd = z * dsdd/s;
double lamida = z * dsda/s;
double plasg = zbar * zbar * esqu * kTinv * inv_lami;
z = -plasg * inv_lami;
double plasgdd = z * lamidd;
double plasgda = z * lamida;
double plasgdt = -plasg*kTinv * kerg;
double plasgdz = 2 * plasg/zbar;
// yakovlev & shalybkov 1989 equations 82, 85, 86, 87
if (plasg >= 1.0) {
x = pow(plasg,0.25);
y = avo * ytot1 * kerg;
ecoul = y * T * (a1*plasg + b1*x + c1/x + d1);
pcoul = third * rho * ecoul;
scoul = -y * (3*b1*x - 5*c1/x + d1 * (log(plasg) - 1) - e1);
y = avo * ytot1 * kT * (a1 + 0.25/plasg*(b1*x - c1/x));
decouldd = y * plasgdd;
decouldt = y * plasgdt + ecoul/T;
decoulda = y * plasgda - ecoul/abar;
decouldz = y * plasgdz;
y = third * rho;
dpcouldd = third * ecoul + y * decouldd;
dpcouldt = y * decouldt;
dpcoulda = y * decoulda;
dpcouldz = y * decouldz;
y = -avo * kerg / (abar*plasg) * (0.75*b1*x+1.25*c1/x+d1);
dscouldd = y * plasgdd;
dscouldt = y * plasgdt;
dscoulda = y * plasgda - scoul/abar;
dscouldz = y * plasgdz;
}
// yakovlev & shalybkov 1989 equations 102, 103, 104
else {
x = plasg*sqrt(plasg);
y = pow(plasg,b2);
z = c2 * x - third * a2 * y;
pcoul = -pion * z;
ecoul = 3 * pcoul/rho;
scoul = -(avo/abar) * kerg * (c2*x -a2*(b2-1)/b2*y);
s = 1.5*c2*x/plasg - third*a2*b2*y/plasg;
dpcouldd = -dpiondd*z - pion*s*plasgdd;
dpcouldt = -dpiondt*z - pion*s*plasgdt;
dpcoulda = -dpionda*z - pion*s*plasgda;
dpcouldz = -dpiondz*z - pion*s*plasgdz;
s = 3/rho;
decouldd = s * dpcouldd - ecoul/rho;
decouldt = s * dpcouldt;
decoulda = s * dpcoulda;
decouldz = s * dpcouldz;
s = -avo*kerg/(abar*plasg)*(1.5*c2*x - a2*(b2-1)*y);
dscouldd = s * plasgdd;
dscouldt = s * plasgdt;
dscoulda = s * plasgda - scoul/abar;
dscouldz = s * plasgdz;
}
// "bomb proof"
x = prad + pion + pele + pcoul;
y = erad + eion + eele + ecoul;
z = srad + sion + sele + scoul;
if(x <= 0 || y <= 0) {
// zero out coulomb terms
pcoul = 0; dpcouldd=0; dpcouldt=0; dpcoulda=0; dpcouldz=0;
ecoul = 0; decouldd=0; decouldt=0; decoulda=0; decouldz=0;
scoul = 0; dscouldd=0; dscouldt=0; dscoulda=0; dscouldz=0;
}
// gas subtotals
pgas = pion + pele + pcoul;
egas = eion + eele + ecoul;
sgas = sion + sele + scoul;
LOG_DEBUG(logger, "pgas = {}", pgas);
LOG_DEBUG(logger, "prad = {}", prad);
LOG_DEBUG(logger, "pion = {}", pion);
LOG_DEBUG(logger, "pele = {}", pele);
LOG_DEBUG(logger, "pcoul = {}", pcoul);
dpgasdd = dpiondd + dpepdd + dpcouldd;
dpgasdt = dpiondt + dpepdt + dpcouldt;
dpgasda = dpionda + dpepda + dpcoulda;
dpgasdz = dpiondz + dpepdz + dpcouldz;
LOG_DEBUG(logger, "dpgasdd = {}", dpgasdd);
LOG_DEBUG(logger, "dpraddd = {}", dpraddd);
LOG_DEBUG(logger, "dpiondd = {}", dpiondd);
LOG_DEBUG(logger, "dpeledd = {}", dpepdd);
LOG_DEBUG(logger, "dpcouldd = {}", dpcouldd);
degasdd = deiondd + deepdd + decouldd;
degasdt = deiondt + deepdt + decouldt;
degasda = deionda + deepda + decoulda;
degasdz = deiondz + deepdz + decouldz;
LOG_DEBUG(logger, "degasdd = {}", degasdd);
LOG_DEBUG(logger, "deraddd = {}", deraddd);
LOG_DEBUG(logger, "deiondd = {}", deiondd);
LOG_DEBUG(logger, "deeledd = {}", deepdd);
LOG_DEBUG(logger, "decouldd = {}", decouldd);
LOG_DEBUG(logger, "degasdt = {}", degasdt);
LOG_DEBUG(logger, "deraddt = {}", deraddt);
LOG_DEBUG(logger, "deiondt = {}", deiondt);
LOG_DEBUG(logger, "deeledt = {}", deepdt);
LOG_DEBUG(logger, "decouldt = {}", decouldt);
dsgasdd = dsiondd + dsepdd + dscouldd;
dsgasdt = dsiondt + dsepdt + dscouldt;
dsgasda = dsionda + dsepda + dscoulda;
dsgasdz = dsiondz + dsepdz + dscouldz;
LOG_DEBUG(logger, "sgas = {}", sgas);
LOG_DEBUG(logger, "srad = {}", srad);
LOG_DEBUG(logger, "sion = {}", sion);
LOG_DEBUG(logger, "sele = {}", sele);
LOG_DEBUG(logger, "scoul = {}", scoul);
LOG_DEBUG(logger, "dsgasdd = {}", dsgasdd);
LOG_DEBUG(logger, "dsraddd = {}", dsraddd);
LOG_DEBUG(logger, "dsiondd = {}", dsiondd);
LOG_DEBUG(logger, "dseledd = {}", dsepdd);
LOG_DEBUG(logger, "dscouldd = {}", dscouldd);
LOG_DEBUG(logger, "dsgasdt = {}", dsgasdt);
LOG_DEBUG(logger, "dsraddt = {}", dsraddt);
LOG_DEBUG(logger, "dsiondt = {}", dsiondt);
LOG_DEBUG(logger, "dseledt = {}", dsepdt);
LOG_DEBUG(logger, "dscouldt = {}", dscouldt);
//total
ptot = prad + pgas;
stot = srad + sgas;
etot = erad + egas;
dpresdd = dpraddd + dpgasdd;
dpresdt = dpraddt + dpgasdt;
dpresda = dpradda + dpgasda;
dpresdz = dpraddz + dpgasdz;
denerdd = deraddd + degasdd;
denerdt = deraddt + degasdt;
denerda = deradda + degasda;
denerdz = deraddz + degasdz;
dentrdd = dsraddd + dsgasdd;
dentrdt = dsraddt + dsgasdt;
dentrda = dsradda + dsgasda;
dentrdz = dsraddz + dsgasdz;
// derived quantities
double zz = ptot * rhoi;
double zzi = rho/ptot;
double chit = (T/ptot)*dpresdt;
double chirho = zzi*dpresdd;
double cv = denerdt;
x = zz * chit/(T*cv);
double gamma3 = x + 1.0;
double gamma1 = chit * x + chirho;
double grad_ad = x/gamma1;
double gamma2 = 1.0/(1.0 - grad_ad);
double cp = cv * gamma1/chirho;
z = 1.0 + (etot + clight2)*zzi;
double csound = clight * sqrt(gamma1/z);
// package in q:
EOS eos;
eos.ptot = ptot; eos.etot = etot; eos.stot = stot;
eos.pgas = pgas; eos.egas = egas; eos.sgas = sgas;
eos.prad = prad; eos.erad = erad; eos.srad = srad;
eos.chiT = chit, eos.chiRho = chirho, eos.csound = csound;
eos.gamma1 = gamma1, eos.gamma2 = gamma2; eos.gamma3 = gamma3;
eos.cV = cv, eos.cP = cp, eos.grad_ad = grad_ad;
eos.etaele = etaele, eos.xnefer = xnefer, eos.ye = ye;
eos.dpresdd = dpresdd; eos.dentrdd = dentrdd; eos.denerdd = denerdd;
eos.dpresdt = dpresdt; eos.dentrdt = dentrdt; eos.denerdt = denerdt;
eos.dpresda = dpresda; eos.dentrda = dentrda; eos.denerda = denerda;
eos.dpresdz = dpresdz; eos.dentrdz = dentrdz; eos.denerdz = denerdz;
// maxwell relations
x = rho * rho;
eos.dse = T * dentrdt/denerdt - 1.0;
eos.dpe = (denerdd*x + T*dpresdt)/ptot - 1.0;
eos.dsp = -dentrdd*x/dpresdt - 1.0;
return eos;
// TODO: In future this might be made more preformant by avoiding the copy.
}
}

865
src/eos/private/helmholtz.cpp
View File

@@ -1,865 +0,0 @@
// this file contains a C++ conversion of Frank Timmes' fortran code
// helmholtz.f90, which implements the Helmholtz Free Energy EOS described
// by Timmes & Swesty (2000) doi:10.1086/313304 -- Aaron Dotter 2025
#include <iostream>
#include <fstream>
#include <sstream>
#include <cmath>
#include <format>
using namespace std;
// interpolating polynomila function definitions
double psi0(double z) {return z*z*z * (z * (-6.0 * z + 15.0) - 10.0) + 1.0;}
double dpsi0(double z) {return z*z * ( z * (-30.0 * z + 60.0) - 30.0);}
double ddpsi0(double z) {return z * ( z * (-120.0 * z + 180.0) - 60.0);}
double psi1(double z) {return z * ( z*z * ( z * (-3.0*z + 8.0) - 6.0) + 1.0);}
double dpsi1(double z) {return z*z * ( z * (-15.0*z + 32.0) - 18.0) +1.0;}
double ddpsi1(double z) {return z * (z * (-60.0*z + 96.0) - 36.0);}
double psi2(double z) {return 0.5 * z*z * ( z* ( z * (-z + 3.0) - 3.0) + 1.0);}
double dpsi2(double z) {return 0.5 * z * ( z * (z * (-5.0*z + 12.0) - 9.0) + 2.0);}
double ddpsi2(double z) {return 0.5*(z*( z * (-20.*z + 36.) - 18.) + 2.);}
double xpsi0(double z) {return z * z * (2.0 * z - 3.0) + 1.0;}
double xdpsi0(double z) {return z * (6.0*z - 6.0);}
double xpsi1(double z) {return z * ( z * (z - 2.0) + 1.0);}
double xdpsi1(double z) {return z * (3.0 * z - 4.0) + 1.0;}
double h3(double fi[36], double w0t, double w1t, double w0mt, double w1mt,
double w0d, double w1d, double w0md, double w1md) {
return fi[0] * w0d * w0t + fi[1] * w0md * w0t
+ fi[2] * w0d * w0mt + fi[3] * w0md * w0mt
+ fi[4] * w0d * w1t + fi[5] * w0md * w1t
+ fi[6] * w0d * w1mt + fi[7] * w0md * w1mt
+ fi[8] * w1d * w0t + fi[9] * w1md * w0t
+ fi[10] * w1d * w0mt + fi[11] * w1md * w0mt
+ fi[12] * w1d * w1t + fi[13] * w1md * w1t
+ fi[14] * w1d * w1mt + fi[15] * w1md * w1mt;
}
double h5(double fi[36], double w0t, double w1t, double w2t,double w0mt,
double w1mt, double w2mt, double w0d, double w1d, double w2d,
double w0md, double w1md, double w2md) {
return fi[0] * w0d * w0t + fi[1] * w0md * w0t
+ fi[2] * w0d * w0mt + fi[3] * w0md * w0mt
+ fi[4] * w0d * w1t + fi[5] * w0md * w1t
+ fi[6] * w0d * w1mt + fi[7] * w0md * w1mt
+ fi[8] * w0d * w2t + fi[ 9] * w0md * w2t
+ fi[10] * w0d * w2mt + fi[11] * w0md * w2mt
+ fi[12] * w1d * w0t + fi[13] * w1md * w0t
+ fi[14] * w1d * w0mt + fi[15] * w1md * w0mt
+ fi[16] * w2d * w0t + fi[17] * w2md * w0t
+ fi[18] * w2d * w0mt + fi[19] * w2md * w0mt
+ fi[20] * w1d * w1t + fi[21] * w1md * w1t
+ fi[22] * w1d * w1mt + fi[23] * w1md * w1mt
+ fi[24] * w2d * w1t + fi[25] * w2md * w1t
+ fi[26] * w2d * w1mt + fi[27] * w2md * w1mt
+ fi[28] * w1d * w2t + fi[29] * w1md * w2t
+ fi[30] * w1d * w2mt + fi[31] * w1md * w2mt
+ fi[32] * w2d * w2t + fi[33] * w2md * w2t
+ fi[34] * w2d * w2mt + fi[35] * w2md * w2mt;
}
// some global quantities including table size, limits, and deltas
const int IMAX=541, JMAX=201;
const double tlo = 3.0, thi = 13.0;
const double dlo = -12.0, dhi = 15.0;
const double tstp = (thi - tlo)/(JMAX-1), tstpi=1/tstp;
const double dstp = (dhi-dlo)/(IMAX-1), dstpi = 1/dstp;
double t[JMAX], d[IMAX];
double f[IMAX][JMAX], fd[IMAX][JMAX], ft[IMAX][JMAX];
double fdd[IMAX][JMAX], ftt[IMAX][JMAX], fdt[IMAX][JMAX];
double fddt[IMAX][JMAX], fdtt[IMAX][JMAX], fddtt[IMAX][JMAX];
double dpdf[IMAX][JMAX], dpdfd[IMAX][JMAX], dpdft[IMAX][JMAX];
double dpdfdt[IMAX][JMAX], ef[IMAX][JMAX], efd[IMAX][JMAX];
double eft[IMAX][JMAX], efdt[IMAX][JMAX], xf[IMAX][JMAX];
double xfd[IMAX][JMAX], xft[IMAX][JMAX], xfdt[IMAX][JMAX];
double dt_sav[JMAX], dt2_sav[JMAX], dti_sav[JMAX], dt2i_sav[JMAX];
double dd_sav[IMAX], dd2_sav[IMAX], ddi_sav[IMAX];
// this function reads in the HELM table and stores in the above arrays
void read_helm_table() {
string data;
int i, j;
//set T and Rho (d) arrays
for (j=0; j<JMAX; j++) { t[j] = exp10(tlo + tstp*j); }
for (i=0; i<IMAX; i++) { d[i] = exp10(dlo + dstp*i); }
ifstream helm_table("helm_table.dat");
//read the Helmholtz free energy and its derivatives
for (j=0; j<JMAX; j++) {
for (i=0; i<IMAX; i++){
getline(helm_table, data);
stringstream id(data);
id >> f[i][j];
id >> fd[i][j];
id >> ft[i][j];
id >> fdd[i][j];
id >> ftt[i][j];
id >> fdt[i][j];
id >> fddt[i][j];
id >> fdtt[i][j];
id >> fddtt[i][j];
}
}
//read the pressure derivative with density
for (j=0; j<JMAX; j++) {
for (i=0; i<IMAX; i++){
getline(helm_table, data);
stringstream id(data);
id >> dpdf[i][j];
id >> dpdfd[i][j];
id >> dpdft[i][j];
id >> dpdfdt[i][j];
}
}
//read the electron chemical potential
for (j=0; j<JMAX; j++) {
for (i=0; i<IMAX; i++){
getline(helm_table, data);
stringstream id(data);
id >> ef[i][j];
id >> efd[i][j];
id >> eft[i][j];
id >> efdt[i][j];
}
}
//read the number density
for (j=0; j<JMAX; j++) {
for (i=0; i<IMAX; i++){
getline(helm_table, data);
stringstream id(data);
id >> xf[i][j];
id >> xfd[i][j];
id >> xft[i][j];
id >> xfdt[i][j];
}
}
helm_table.close(); //done reading
// construct the temperature and density deltas and their inverses
for (j=0; j<JMAX; j++) {
double dth = t[j+1] - t[j];
double dt2 = dth * dth;
double dti = 1.0/dth;
double dt2i = 1.0/dt2;
dt_sav[j] = dth;
dt2_sav[j] = dt2;
dti_sav[j] = dti;
dt2i_sav[j] = dt2i;
}
for (i=0; i<IMAX; i++) {
double dd = d[i+1] - d[i];
double dd2 = dd * dd;
double ddi = 1.0/dd;
dd_sav[i] = dd;
dd2_sav[i] = dd2;
ddi_sav[i] = ddi;
}
}
//this class stores the information for one EOS call
class EOS {
public:
//input
double T, rho, abar, zbar;
//output
double ye, etaele, xnefer; //
double ptot, pgas, prad; // pressure
double etot, egas, erad; // energy
double stot, sgas, srad; // entropy
//derivatives
double dpresdd, dpresdt, dpresda, dpresdz;
double dentrdd, dentrdt, dentrda, dentrdz;
double denerdd, denerdt, denerda, denerdz;
//these could become functions:
double chiT, chiRho, csound, grad_ad; // derived quantities
double gamma1, gamma2, gamma3, cV, cP; // derived quantities
double dse, dpe, dsp; // Maxwell relations
};
/***
this function evaluates the EOS components:
ion, radiation, electron-positron and Coulomb interaction
and returns the calculated quantities in the input
***/
void helmEOS(EOS &q) {
const double pi = 3.1415926535897932384;
const double amu = 1.66053878283e-24;
const double h = 6.6260689633e-27;
const double avo = 6.0221417930e23;
const double kerg = 1.380650424e-16;
const double kergavo = kerg*avo;
const double clight = 2.99792458e10;
const double clight2 = clight * clight;
const double ssol = 5.6704e-5;
const double asol = 4 * ssol / clight;
const double asoli3 = asol / 3;
const double sioncon = (2.0 * pi * amu * kerg)/(h*h);
const double qe = 4.8032042712e-10;
const double esqu = qe * qe;
const double third = 1./3.;
double T, rho, s, fi[36], free, abar, zbar, ytot1, ye, ptot;
double prad, srad, erad, etot, stot;
double dpraddd, dpraddt, dpradda, dpraddz;
double deraddd, deraddt, deradda, deraddz;
double dsraddd, dsraddt, dsradda, dsraddz;
double dpiondd, dpiondt, dpionda, dpiondz;
double deiondd, deiondt, deionda, deiondz;
double dsiondd, dsiondt, dsionda, dsiondz;
double dpcouldd, dpcouldt, dpcoulda, dpcouldz;
double decouldd, decouldt, decoulda, decouldz;
double dscouldd, dscouldt, dscoulda, dscouldz;
double dpgasdd, dpgasdt, dpgasda, dpgasdz;
double degasdd, degasdt, degasda, degasdz;
double dsgasdd, dsgasdt, dsgasda, dsgasdz;
double dpepdd, dpepdt, dpepda, dpepdz;
double dsepdd, dsepdt, dsepda, dsepdz;
double deepdd, deepdt, deepda, deepdz;
double dxnedd, dxnedt, dxneda, dxnedz;
double dxnidd, dxnida, detadz, detada;
double dpresdd, dpresdt, dpresda, dpresdz;
double denerdd, denerdt, denerda, denerdz;
double dentrdd, dentrdt, dentrda, dentrdz;
double xni, sion, pion, eion, x, y, z;
double pgas, egas, sgas, pele, sele, eele, pcoul, scoul, ecoul;
int jat, iat;
T = q.T; rho = q.rho; abar = q.abar; zbar = q.zbar;
ytot1 = 1/abar;
ye = zbar * ytot1;
double rhoi = 1/rho;
double Ti = 1/T;
double kT= kerg*T;
double kTinv = 1/kT;
/*** radiation ***/
prad = asoli3 * pow(T,4);
dpraddd = 0.;
dpraddt = 4 * prad * Ti;
dpradda = 0.;
dpraddz = 0.;
erad = 3 * prad * rhoi;
deraddd = -erad * rhoi;
deraddt = 3 * dpraddt * rhoi;
deradda = 0.;
deraddz = 0.;
srad = (prad*rhoi + erad)*Ti;
dsraddd = (dpraddd*rhoi - prad*rhoi*rhoi + deraddd)*Ti;
dsraddt = (dpraddt*rhoi + deraddt - srad)*Ti;
dsradda = 0.;
dsraddz = 0.;
/*** ions ***/
xni = avo * ytot1 * rho;
dxnidd = avo * ytot1;
dxnida = -xni * ytot1;
pion = xni * kT;
dpiondd = dxnidd * kT;
dpiondt = xni * kerg;
dpionda = dxnida * kT;
dpiondz = 0.;
eion = 1.5 * pion * rhoi;
deiondd = (1.5 * dpiondd - eion)*rhoi;
deiondt = 1.5 * dpiondt*rhoi;
deionda = 1.5 * dpionda*rhoi;
deiondz = 0.;
x = abar * abar * sqrt(abar) * rhoi/avo;
s = sioncon * T;
z = x * s * sqrt(s);
y = log(z);
sion = (pion * rhoi + eion) * Ti + kergavo * ytot1 * y;
dsiondd = (dpiondd*rhoi - pion*rhoi*rhoi + deiondd)*Ti
- kergavo * rhoi * ytot1;
dsiondt = (dpiondt*rhoi + deiondt)*Ti - (pion*rhoi + eion) * Ti*Ti
+ 1.5 * kergavo * Ti*ytot1;
x = avo*kerg/abar;
dsionda = (dpionda*rhoi + deionda)*Ti + kergavo*ytot1*ytot1* (2.5 - y);
dsiondz = 0.;
// electron-positron
double xnem = xni * zbar;
double din = ye*rho;
// hash the table location in T, rho:
jat = (log10(T) - tlo)*tstpi;
jat = max(0, min(jat, JMAX-1)); // bounds
iat = max(0, min(iat, IMAX-1)); // bounds
iat = (log10(din)-dlo)*dstpi;
//using the indices, access the table:
fi[0] = f[iat][jat];
fi[1] = f[iat+1][jat];
fi[2] = f[iat][jat+1];
fi[3] = f[iat+1][jat+1];
fi[4] = ft[iat][jat];
fi[5] = ft[iat+1][jat];
fi[6] = ft[iat][jat+1];
fi[7] = ft[iat+1][jat+1];
fi[8] = ftt[iat][jat];
fi[9] = ftt[iat+1][jat];
fi[10] = ftt[iat][jat+1];
fi[11] = ftt[iat+1][jat+1];
fi[12] = fd[iat][jat];
fi[13] = fd[iat+1][jat];
fi[14] = fd[iat][jat+1];
fi[15] = fd[iat+1][jat+1];
fi[16] = fdd[iat][jat];
fi[17] = fdd[iat+1][jat];
fi[18] = fdd[iat][jat+1];
fi[19] = fdd[iat+1][jat+1];
fi[20] = fdt[iat][jat];
fi[21] = fdt[iat+1][jat];
fi[22] = fdt[iat][jat+1];
fi[23] = fdt[iat+1][jat+1];
fi[24] = fddt[iat][jat];
fi[25] = fddt[iat+1][jat];
fi[26] = fddt[iat][jat+1];
fi[27] = fddt[iat+1][jat+1];
fi[28] = fdtt[iat][jat];
fi[29] = fdtt[iat+1][jat];
fi[30] = fdtt[iat][jat+1];
fi[31] = fdtt[iat+1][jat+1];
fi[32] = fddtt[iat][jat];
fi[33] = fddtt[iat+1][jat];
fi[34] = fddtt[iat][jat+1];
fi[35] = fddtt[iat+1][jat+1];
double xt = (T - t[jat]) * dti_sav[jat];
double xd = (din - d[iat]) * ddi_sav[iat];
double mxt = 1-xt;
double mxd = 1-xd;
//density and temperature basis functions
double si0t = psi0(xt);
double si1t = psi1(xt)*dt_sav[jat];
double si2t = psi2(xt)*dt2_sav[jat];
double si0mt = psi0(mxt);
double si1mt = -psi1(mxt)*dt_sav[jat];
double si2mt = psi2(mxt)*dt2_sav[jat];
double si0d = psi0(xd);
double si1d = psi1(xd)*dd_sav[iat];
double si2d = psi2(xd)*dd2_sav[iat];
double si0md = psi0(mxd);
double si1md = -psi1(mxd)*dd_sav[iat];
double si2md = psi2(mxd)*dd2_sav[iat];
// derivatives of the weight functions
double dsi0t = dpsi0(xt)*dti_sav[jat];
double dsi1t = dpsi1(xt);
double dsi2t = dpsi2(xt)*dt_sav[jat];
double dsi0mt = -dpsi0(mxt)*dti_sav[jat];
double dsi1mt = dpsi1(mxt);
double dsi2mt = -dpsi2(mxt)*dt_sav[jat];
double dsi0d = dpsi0(xd)*ddi_sav[iat];
double dsi1d = dpsi1(xd);
double dsi2d = dpsi2(xd)*dd_sav[iat];
double dsi0md = -dpsi0(mxd)*ddi_sav[iat];
double dsi1md = dpsi1(mxd);
double dsi2md = -dpsi2(mxd)*dd_sav[iat];
// second derivatives of the weight functions
double ddsi0t = ddpsi0(xt)*dt2i_sav[jat];
double ddsi1t = ddpsi1(xt)*dti_sav[jat];
double ddsi2t = ddpsi2(xt);
double ddsi0mt = ddpsi0(mxt)*dt2i_sav[jat];
double ddsi1mt = -ddpsi1(mxt)*dti_sav[jat];
double ddsi2mt = ddpsi2(mxt);
cout << " xt, mxt = " << xt << " " << mxt << endl;
cout << " dti_sav[jat] = " << dti_sav[jat] << endl;
cout << " dt2i_sav[jat] = " << dt2i_sav[jat] << endl;
cout << " ddpsi0(xt) = " << ddpsi0(xt) << endl;
cout << " ddpsi1(xt) = " << ddpsi1(xt) << endl;
cout << " ddpsi2(xt) = " << ddpsi2(xt) << endl;
cout << " ddpsi0(mxt) = " << ddpsi0(mxt) << endl;
cout << " ddpsi1(mxt) = " << ddpsi1(mxt) << endl;
cout << " ddpsi2(mxt) = " << ddpsi2(mxt) << endl;
cout << " ddsi0t = " << ddsi0t << endl;
cout << " ddsi1t = " << ddsi1t << endl;
cout << " ddsi2t = " << ddsi2t << endl;
cout << " ddsi0mt = " << ddsi0mt << endl;
cout << " ddsi1mt = " << ddsi1mt << endl;
cout << " ddsi2mt = " << ddsi2mt << endl;
ddsi0t = 0.0; ddsi1t = 0.0; ddsi2t = 1.0;
ddsi0mt = 0.0; ddsi1mt = 0.0; ddsi2mt = 0.0;
// free energy
free = h5(fi,si0t, si1t, si2t, si0mt, si1mt, si2mt,
si0d, si1d, si2d, si0md, si1md, si2md);
// derivative with respect to density
double df_d = h5(fi, si0t, si1t, si2t, si0mt, si1mt, si2mt,
dsi0d, dsi1d, dsi2d, dsi0md, dsi1md, dsi2md);
// derivative with respect to temperature
double df_t = h5(fi, dsi0t, dsi1t, dsi2t, dsi0mt, dsi1mt, dsi2mt,
si0d, si1d, si2d, si0md, si1md, si2md);
// second derivative with respect to temperature
double df_tt = h5(fi, ddsi0t, ddsi1t, ddsi2t, ddsi0mt, ddsi1mt, ddsi2mt,
si0d, si1d, si2d, si0md, si1md, si2md);
cout << "df_tt = " << df_tt << endl;
// derivative with respect to temperature and density
double df_dt = h5(fi, dsi0t, dsi1t, dsi2t, dsi0mt, dsi1mt, dsi2mt,
dsi0d, dsi1d, dsi2d, dsi0md, dsi1md, dsi2md);
// pressure derivatives
si0t = xpsi0(xt);
si1t = xpsi1(xt)*dt_sav[jat];
si0mt = xpsi0(mxt);
si1mt = -xpsi1(mxt)*dt_sav[jat];
si0d = xpsi0(xd);
si1d = xpsi1(xd)*dd_sav[iat];
si0md = xpsi0(mxd);
si1md = -xpsi1(mxd)*dd_sav[iat];
// derivatives of weight functions
dsi0t = xdpsi0(xt)*dti_sav[jat];
dsi1t = xdpsi1(xt);
dsi0mt = -xdpsi0(mxt)*dti_sav[jat];
dsi1mt = xdpsi1(mxt);
dsi0d = xdpsi0(xd)*ddi_sav[iat];
dsi1d = xdpsi1(xd);
dsi0md = -xdpsi0(mxd)*ddi_sav[iat];
dsi1md = xdpsi1(mxd);
// pressure derivative
fi[0] = dpdf[iat][jat];
fi[1] = dpdf[iat+1][jat];
fi[2] = dpdf[iat][jat+1];
fi[3] = dpdf[iat+1][jat+1];
fi[4] = dpdft[iat][jat];
fi[5] = dpdft[iat+1][jat];
fi[6] = dpdft[iat][jat+1];
fi[7] = dpdft[iat+1][jat+1];
fi[8] = dpdfd[iat][jat];
fi[9] = dpdfd[iat+1][jat];
fi[10] = dpdfd[iat][jat+1];
fi[11] = dpdfd[iat+1][jat+1];
fi[12] = dpdfdt[iat][jat];
fi[13] = dpdfdt[iat+1][jat];
fi[14] = dpdfdt[iat][jat+1];
fi[15] = dpdfdt[iat+1][jat+1];
// pressure derivative with density
dpepdd = h3(fi, si0t, si1t, si0mt, si1mt, si0d, si1d, si0md, si1md);
dpepdd = max(ye * dpepdd,1e-30);
// look in the electron chemical potential table only once
fi[0] = ef[iat][jat];
fi[1] = ef[iat+1][jat];
fi[2] = ef[iat][jat+1];
fi[3] = ef[iat+1][jat+1];
fi[4] = eft[iat][jat];
fi[5] = eft[iat+1][jat];
fi[6] = eft[iat][jat+1];
fi[7] = eft[iat+1][jat+1];
fi[8] = efd[iat][jat];
fi[9] = efd[iat+1][jat];
fi[10] = efd[iat][jat+1];
fi[11] = efd[iat+1][jat+1];
fi[12] = efdt[iat][jat];
fi[13] = efdt[iat+1][jat];
fi[14] = efdt[iat][jat+1];
fi[15] = efdt[iat+1][jat+1];
// electron chemical potential etaele
double etaele = h3(fi, si0t, si1t, si0mt, si1mt, si0d, si1d, si0md, si1md);
// derivative with respect to density
x = h3(fi, si0t, si1t, si0mt, si1mt, dsi0d, dsi1d, dsi0md, dsi1md);
double detadd = ye * x;
// derivative with respect to temperature
double detadt = h3(fi, dsi0t, dsi1t, dsi0mt, dsi1mt, si0d,
si1d, si0md, si1md);
// derivative with respect to abar and zbar
detada = -x * din * ytot1;
detadz = x * rho * ytot1;
// look in the number density table only once
fi[0] = xf[iat][jat];
fi[1] = xf[iat+1][jat];
fi[2] = xf[iat][jat+1];
fi[3] = xf[iat+1][jat+1];
fi[4] = xft[iat][jat];
fi[5] = xft[iat+1][jat];
fi[6] = xft[iat][jat+1];
fi[7] = xft[iat+1][jat+1];
fi[8] = xfd[iat][jat];
fi[9] = xfd[iat+1][jat];
fi[10] = xfd[iat][jat+1];
fi[11] = xfd[iat+1][jat+1];
fi[12] = xfdt[iat][jat];
fi[13] = xfdt[iat+1][jat];
fi[14] = xfdt[iat][jat+1];
fi[15] = xfdt[iat+1][jat+1];
// electron chemical potential etaele
double xnefer = h3(fi, si0t, si1t, si0mt, si1mt, si0d, si1d, si0md, si1md);
// derivative with respect to density
x = h3(fi, si0t, si1t, si0mt, si1mt, dsi0d, dsi1d, dsi0md, dsi1md);
x = max(x,1e-30);
dxnedd = ye * x;
// derivative with respect to temperature
dxnedt = h3(fi, dsi0t, dsi1t, dsi0mt, dsi1mt, si0d, si1d, si0md, si1md);
// derivative with respect to abar and zbar
dxneda = -x * din * ytot1;
dxnedz = x * rho * ytot1;
// the desired electron-positron thermodynamic quantities
x = din * din;
pele = x * df_d;
dpepdt = x * df_dt;
//dpepdd is set above
s = dpepdd/ye - 2 * din * df_d;
dpepda = -ytot1 * (2 * pele + s * din);
dpepdz = rho * ytot1 * (2 * din * df_d + s);
x = ye * ye;
sele = -df_t * ye;
dsepdt = -df_tt * ye;
dsepdd = -df_dt * x;
dsepda = ytot1 * (ye * df_dt * din - sele);
dsepdz = -ytot1 * (ye * df_dt * rho + df_t);
eele = ye*free + T * sele;
deepdt = T * dsepdt;
deepdd = x * df_d + T * dsepdd;
deepda = -ye * ytot1 * (free + df_d * din) + T * dsepda;
deepdz = ytot1 * (free + ye * df_d * rho) + T * dsepdz;
// coulomb
const double a1 = -0.898004;
const double b1 = 0.96786;
const double c1 = 0.220703;
const double d1 = -0.86097;
const double e1 = 2.5269;
const double a2 = 0.29561;
const double b2 = 1.9885;
const double c2 = 0.288675;
z = 4 * pi * third;
s = z * xni;
double dsdd = z * dxnidd;
double dsda = z * dxnida;
double lami = pow(s,-1./3.);
double inv_lami = 1.0/lami;
z = -lami/3.;
double lamidd = z * dsdd/s;
double lamida = z * dsda/s;
double plasg = zbar * zbar * esqu * kTinv * inv_lami;
z = -plasg * inv_lami;
double plasgdd = z * lamidd;
double plasgda = z * lamida;
double plasgdt = -plasg*kTinv * kerg;
double plasgdz = 2 * plasg/zbar;
// yakovlev & shalybkov 1989 equations 82, 85, 86, 87
if (plasg >= 1.0) {
x = pow(plasg,0.25);
y = avo * ytot1 * kerg;
ecoul = y * T * (a1*plasg + b1*x + c1/x + d1);
pcoul = third * rho * ecoul;
scoul = -y * (3*b1*x - 5*c1/x + d1 * (log(plasg) - 1) - e1);
y = avo * ytot1 * kT * (a1 + 0.25/plasg*(b1*x - c1/x));
decouldd = y * plasgdd;
decouldt = y * plasgdt + ecoul/T;
decoulda = y * plasgda - ecoul/abar;
decouldz = y * plasgdz;
y = third * rho;
dpcouldd = third * ecoul + y * decouldd;
dpcouldt = y * decouldt;
dpcoulda = y * decoulda;
dpcouldz = y * decouldz;
y = -avo * kerg / (abar*plasg) * (0.75*b1*x+1.25*c1/x+d1);
dscouldd = y * plasgdd;
dscouldt = y * plasgdt;
dscoulda = y * plasgda - scoul/abar;
dscouldz = y * plasgdz;
}
// yakovlev & shalybkov 1989 equations 102, 103, 104
else {
x = plasg*sqrt(plasg);
y = pow(plasg,b2);
z = c2 * x - third * a2 * y;
pcoul = -pion * z;
ecoul = 3 * pcoul/rho;
scoul = -(avo/abar) * kerg * (c2*x -a2*(b2-1)/b2*y);
s = 1.5*c2*x/plasg - third*a2*b2*y/plasg;
dpcouldd = -dpiondd*z - pion*s*plasgdd;
dpcouldt = -dpiondt*z - pion*s*plasgdt;
dpcoulda = -dpionda*z - pion*s*plasgda;
dpcouldz = -dpiondz*z - pion*s*plasgdz;
s = 3/rho;
decouldd = s * dpcouldd - ecoul/rho;
decouldt = s * dpcouldt;
decoulda = s * dpcoulda;
decouldz = s * dpcouldz;
s = -avo*kerg/(abar*plasg)*(1.5*c2*x - a2*(b2-1)*y);
dscouldd = s * plasgdd;
dscouldt = s * plasgdt;
dscoulda = s * plasgda - scoul/abar;
dscouldz = s * plasgdz;
}
// "bomb proof"
x = prad + pion + pele + pcoul;
y = erad + eion + eele + ecoul;
z = srad + sion + sele + scoul;
if(x <= 0 || y <= 0) {
// zero out coulomb terms
pcoul = 0; dpcouldd=0; dpcouldt=0; dpcoulda=0; dpcouldz=0;
ecoul = 0; decouldd=0; decouldt=0; decoulda=0; decouldz=0;
scoul = 0; dscouldd=0; dscouldt=0; dscoulda=0; dscouldz=0;
}
// gas subtotals
pgas = pion + pele + pcoul;
egas = eion + eele + ecoul;
sgas = sion + sele + scoul;
cout << "pgas = " << pgas << endl;
cout << "prad = " << prad << endl;
cout << "pion = " << pion << endl;
cout << "pele = " << pele << endl;
cout << "pcoul= " << pcoul << endl;
dpgasdd = dpiondd + dpepdd + dpcouldd;
dpgasdt = dpiondt + dpepdt + dpcouldt;
dpgasda = dpionda + dpepda + dpcoulda;
dpgasdz = dpiondz + dpepdz + dpcouldz;
cout << "dpgasdd = " << dpgasdd << endl;
cout << "dpraddd = " << dpraddd << endl;
cout << "dpiondd = " << dpiondd << endl;
cout << "dpeledd = " << dpepdd << endl;
cout << "dpcouldd= " << dpcouldd << endl;
degasdd = deiondd + deepdd + decouldd;
degasdt = deiondt + deepdt + decouldt;
degasda = deionda + deepda + decoulda;
degasdz = deiondz + deepdz + decouldz;
cout << "degasdd = " << degasdd << endl;
cout << "deraddd = " << deraddd << endl;
cout << "deiondd = " << deiondd << endl;
cout << "deeledd = " << deepdd << endl;
cout << "decouldd= " << decouldd << endl;
cout << "degasdt = " << degasdt << endl;
cout << "deraddt = " << deraddt << endl;
cout << "deiondt = " << deiondt << endl;
cout << "deeledt = " << deepdt << endl;
cout << "decouldt= " << decouldt << endl;
dsgasdd = dsiondd + dsepdd + dscouldd;
dsgasdt = dsiondt + dsepdt + dscouldt;
dsgasda = dsionda + dsepda + dscoulda;
dsgasdz = dsiondz + dsepdz + dscouldz;
cout << "sgas = " << sgas << endl;
cout << "srad = " << srad << endl;
cout << "sion = " << sion << endl;
cout << "sele = " << sele << endl;
cout << "scoul= " << scoul << endl;
cout << "dsgasdd = " << dsgasdd << endl;
cout << "dsraddd = " << dsraddd << endl;
cout << "dsiondd = " << dsiondd << endl;
cout << "dseledd = " << dsepdd << endl;
cout << "dscouldd= " << dscouldd << endl;
cout << "dsgasdt = " << dsgasdt << endl;
cout << "dsraddt = " << dsraddt << endl;
cout << "dsiondt = " << dsiondt << endl;
cout << "dseledt = " << dsepdt << endl;
cout << "dscouldt= " << dscouldt << endl;
//total
ptot = prad + pgas;
stot = srad + sgas;
etot = erad + egas;
dpresdd = dpraddd + dpgasdd;
dpresdt = dpraddt + dpgasdt;
dpresda = dpradda + dpgasda;
dpresdz = dpraddz + dpgasdz;
denerdd = deraddd + degasdd;
denerdt = deraddt + degasdt;
denerda = deradda + degasda;
denerdz = deraddz + degasdz;
dentrdd = dsraddd + dsgasdd;
dentrdt = dsraddt + dsgasdt;
dentrda = dsradda + dsgasda;
dentrdz = dsraddz + dsgasdz;
// derived quantities
double zz = ptot * rhoi;
double zzi = rho/ptot;
double chiT = (T/ptot)*dpresdt;
double chiRho = zzi*dpresdd;
double cV = denerdt;
x = zz * chiT/(T*cV);
double gamma3 = x + 1.0;
double gamma1 = chiT * x + chiRho;
double grad_ad = x/gamma1;
double gamma2 = 1.0/(1.0 - grad_ad);
double cP = cV * gamma1/chiRho;
z = 1.0 + (etot + clight2)*zzi;
double csound = clight * sqrt(gamma1/z);
// package in q:
q.ptot = ptot; q.etot = etot; q.stot = stot;
q.pgas = pgas; q.egas = egas; q.sgas = sgas;
q.prad = prad; q.erad = erad; q.srad = srad;
q.chiT = chiT, q.chiRho = chiRho, q.csound = csound;
q.gamma1 = gamma1, q.gamma2 = gamma2; q.gamma3 = gamma3;
q.cV = cV, q.cP = cP, q.grad_ad = grad_ad;
q.etaele = etaele, q.xnefer = xnefer, q.ye = ye;
q.dpresdd = dpresdd; q.dentrdd = dentrdd; q.denerdd = denerdd;
q.dpresdt = dpresdt; q.dentrdt = dentrdt; q.denerdt = denerdt;
q.dpresda = dpresda; q.dentrda = dentrda; q.denerda = denerda;
q.dpresdz = dpresdz; q.dentrdz = dentrdz; q.denerdz = denerdz;
// Maxwell relations
x = rho * rho;
q.dse = T * dentrdt/denerdt - 1.0;
q.dpe = (denerdd*x + T*dpresdt)/ptot - 1.0;
q.dsp = -dentrdd*x/dpresdt - 1.0;
}
int main() {
int i;
double asum, zsum;
const int nel=3;
double xmass[nel], aion[nel], zion[nel];
EOS eos1;
xmass[0] = 0.75; aion[0] = 1.0; zion[0] = 1.0;
xmass[1] = 0.23; aion[1] = 4.0; zion[1] = 2.0;
xmass[2] = 0.02; aion[2] = 12.0; zion[2] = 6.0;
eos1.T = 1.0e8;
eos1.rho = 1.0e6;
asum = 0.0;
zsum = 0.0;
for (i=0; i<nel; i++) {
asum += xmass[i]/aion[i];
zsum += xmass[i]*zion[i]/aion[i];
}
eos1.abar = 1.0/asum;
eos1.zbar = eos1.abar*zsum;
read_helm_table();
helmEOS(eos1);
cout << " T=" << eos1.T << " rho = " << eos1.rho << " ye = " << eos1.ye << endl;
cout << "helm: " << " value " << " d/drho" << " d/dT " << " d/da" << " d/dz" << std::scientific << endl;
cout << "p tot=" << " " << eos1.ptot << " " << eos1.dpresdd << " " << eos1.dpresdt << " " << eos1.dpresda << " " << eos1.dpresdz << endl;
cout << "e tot=" << " " << eos1.etot << " " << eos1.denerdd << " " << eos1.denerdt << " " << eos1.denerda << " " << eos1.denerdz << endl;
cout << "s tot=" << " " << eos1.stot << " " << eos1.dentrdd << " " << eos1.dentrdt << " " << eos1.dentrda << " " << eos1.dentrdz << endl;
cout << " chiT = " << format("{0:14.8e}", eos1.chiT) << endl;
cout << " chiRho = " << format("{0:14.8e}", eos1.chiRho) << endl;
cout << " etaele = " << format("{0:14.8e}", eos1.etaele) << endl;
cout << " xnefer = " << format("{0:14.8e}", eos1.xnefer) << endl;
cout << " cV = " << format("{0:14.8e}", eos1.cV) << endl;
cout << " cP = " << format("{0:14.8e}", eos1.cP) << endl;
cout << " gamma1 = " << format("{0:14.8e}", eos1.gamma1) << endl;
cout << " gamma2 = " << format("{0:14.8e}", eos1.gamma2) << endl;
cout << " gamma3 = " << format("{0:14.8e}", eos1.gamma3) << endl;
cout << " csound = " << format("{0:14.8e}", eos1.csound) << endl;
cout << " Maxwell Relations " << endl;
cout << " dse = " << format("{0:14.8e}", eos1.dse) << endl;
cout << " dpe = " << format("{0:14.8e}", eos1.dpe) << endl;
cout << " dsp = " << format("{0:14.8e}", eos1.dsp) << endl;
return 0;
}

359
src/eos/public/helm.h Normal file
View File

@@ -0,0 +1,359 @@
// this file contains a C++ conversion of Frank Timmes' fortran code
// helmholtz.f90, which implements the Helmholtz Free Energy EOS described
// by Timmes & Swesty (2000) doi:10.1086/313304 -- Aaron Dotter 2025
#ifndef HELM_H
#define HELM_H
#define IMAX 541
#define JMAX 201
#include <array>
#include <iostream>
#include <iomanip>
#include <sstream>
#include <string>
/**
* @brief 2D array template alias.
* @tparam T Type of the array elements.
* @tparam J Number of columns.
* @tparam I Number of rows.
*/
template <typename T, std::size_t J, std::size_t I>
using array2D = std::array<std::array<T, J>, I>;
double** heap_allocate_contiguous_2D_memory(int rows, int cols);
void heap_deallocate_contiguous_2D_memory(double **array);
namespace helmholtz
{
const double tlo = 3.0, thi = 13.0;
const double dlo = -12.0, dhi = 15.0;
const double tstp = (thi - tlo) / (JMAX - 1), tstpi = 1 / tstp;
const double dstp = (dhi - dlo) / (IMAX - 1), dstpi = 1 / dstp;
/**
* @struct HELMTable
* @brief Structure to hold the Helmholtz EOS table data.
*/
struct HELMTable
{
bool loaded = false;
const int imax = IMAX, jmax = JMAX;
double t[JMAX], d[IMAX];
double** f;
double** fd;
double** ft;
double** fdd;
double** ftt;
double** fdt;
double** fddt;
double** fdtt;
double** fddtt;
double** dpdf;
double** dpdfd;
double** dpdft;
double** dpdfdt;
double** ef;
double** efd;
double** eft;
double** efdt;
double** xf;
double** xfd;
double** xft;
double** xfdt;
double dt_sav[JMAX], dt2_sav[JMAX], dti_sav[JMAX], dt2i_sav[JMAX];
double dd_sav[IMAX], dd2_sav[IMAX], ddi_sav[IMAX];
// Constructor
HELMTable() {
f = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
fd = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
ft = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
fdd = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
ftt = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
fdt = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
fddt = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
fdtt = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
fddtt = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
dpdf = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
dpdfd = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
dpdft = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
dpdfdt = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
ef = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
efd = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
eft = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
efdt = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
xf = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
xfd = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
xft = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
xfdt = heap_allocate_contiguous_2D_memory(IMAX, JMAX);
}
// Destructor
~HELMTable() {
heap_deallocate_contiguous_2D_memory(f);
heap_deallocate_contiguous_2D_memory(fd);
heap_deallocate_contiguous_2D_memory(ft);
heap_deallocate_contiguous_2D_memory(fdd);
heap_deallocate_contiguous_2D_memory(ftt);
heap_deallocate_contiguous_2D_memory(fdt);
heap_deallocate_contiguous_2D_memory(fddt);
heap_deallocate_contiguous_2D_memory(fdtt);
heap_deallocate_contiguous_2D_memory(fddtt);
heap_deallocate_contiguous_2D_memory(dpdf);
heap_deallocate_contiguous_2D_memory(dpdfd);
heap_deallocate_contiguous_2D_memory(dpdft);
heap_deallocate_contiguous_2D_memory(dpdfdt);
heap_deallocate_contiguous_2D_memory(ef);
heap_deallocate_contiguous_2D_memory(efd);
heap_deallocate_contiguous_2D_memory(eft);
heap_deallocate_contiguous_2D_memory(efdt);
heap_deallocate_contiguous_2D_memory(xf);
heap_deallocate_contiguous_2D_memory(xfd);
heap_deallocate_contiguous_2D_memory(xft);
heap_deallocate_contiguous_2D_memory(xfdt);
}
friend std::ostream& operator<<(std::ostream& os, const helmholtz::HELMTable& table) {
if (!table.loaded) {
os << "HELMTable not loaded\n";
return os;
}
double tMin = std::numeric_limits<double>::max(), tMax = std::numeric_limits<double>::lowest();
double dMin = std::numeric_limits<double>::max(), dMax = std::numeric_limits<double>::lowest();
for (const auto& temp : table.t) {
if (temp < tMin) tMin = temp;
if (temp > tMax) tMax = temp;
}
for (const auto& dens : table.d) {
if (dens < dMin) dMin = dens;
if (dens > dMax) dMax = dens;
}
os << "HELMTable Data:\n";
os << " imax: " << table.imax << ", jmax: " << table.jmax << "\n";
os << " Temperature Range: [" << tMin << ", " << tMax << "]\n";
os << " Density Range: [" << dMin << ", " << dMax << "]\n";
return os;
}
};
/**
* @struct EOSInput
* @brief Structure to hold the input parameters for the EOS calculation.
*/
struct EOSInput
{
double T; ///< Temperature.
double rho; ///< Density.
double abar; ///< Mean atomic mass.
double zbar; ///< Mean atomic number.
friend std::ostream& operator<<(std::ostream& os, const helmholtz::EOSInput& eosInput) {
os << "EOSInput Data:\n";
os << " Temperature: " << eosInput.T << "\n";
os << " Density: " << eosInput.rho << "\n";
os << " Mean Atomic Mass: " << eosInput.abar << "\n";
os << " Mean Atomic Number: " << eosInput.zbar << "\n";
return os;
}
};
/**
* @struct EOS
* @brief Structure to hold the output parameters and derivatives of the EOS calculation.
*/
struct EOS
{
// output
double ye, etaele, xnefer; //
double ptot, pgas, prad; // pressure
double etot, egas, erad; // energy
double stot, sgas, srad; // entropy
// derivatives
double dpresdd, dpresdt, dpresda, dpresdz;
double dentrdd, dentrdt, dentrda, dentrdz;
double denerdd, denerdt, denerda, denerdz;
// these could become functions:
double chiT, chiRho, csound, grad_ad; // derived quantities
double gamma1, gamma2, gamma3, cV, cP; // derived quantities
double dse, dpe, dsp; // Maxwell relations
friend std::ostream& operator<<(std::ostream& os, const helmholtz::EOS& eos) {
os << "EOS Data:\n" << std::setw(20) << std::left;
os << " Electron Fraction: " << eos.ye << "\n";
os << " Electron Chemical Potential: " << eos.etaele << "\n";
os << " Electron Number Density: " << eos.xnefer << "\n";
os << " Total Pressure: " << eos.ptot << "\n";
os << " Gas Pressure: " << eos.pgas << "\n";
os << " Radiation Pressure: " << eos.prad << "\n";
os << " Total Energy: " << eos.etot << "\n";
os << " Gas Energy: " << eos.egas << "\n";
os << " Radiation Energy: " << eos.erad << "\n";
os << " Total Entropy: " << eos.stot << "\n";
os << " Gas Entropy: " << eos.sgas << "\n";
os << " Radiation Entropy: " << eos.srad;
return os;
}
};
// interpolating polynomial function definitions
/**
* @brief Interpolating polynomial function psi0.
* @param z Input value.
* @return Result of the polynomial.
*/
double psi0(double z);
/**
* @brief Derivative of the interpolating polynomial function psi0.
* @param z Input value.
* @return Derivative of the polynomial.
*/
double dpsi0(double z);
/**
* @brief Second derivative of the interpolating polynomial function psi0.
* @param z Input value.
* @return Second derivative of the polynomial.
*/
double ddpsi0(double z);
/**
* @brief Interpolating polynomial function psi1.
* @param z Input value.
* @return Result of the polynomial.
*/
double psi1(double z);
/**
* @brief Derivative of the interpolating polynomial function psi1.
* @param z Input value.
* @return Derivative of the polynomial.
*/
double dpsi1(double z);
/**
* @brief Second derivative of the interpolating polynomial function psi1.
* @param z Input value.
* @return Second derivative of the polynomial.
*/
double ddpsi1(double z);
/**
* @brief Interpolating polynomial function psi2.
* @param z Input value.
* @return Result of the polynomial.
*/
double psi2(double z);
/**
* @brief Derivative of the interpolating polynomial function psi2.
* @param z Input value.
* @return Derivative of the polynomial.
*/
double dpsi2(double z);
/**
* @brief Second derivative of the interpolating polynomial function psi2.
* @param z Input value.
* @return Second derivative of the polynomial.
*/
double ddpsi2(double z);
/**
* @brief Interpolating polynomial function xpsi0.
* @param z Input value.
* @return Result of the polynomial.
*/
double xpsi0(double z);
/**
* @brief Derivative of the interpolating polynomial function xpsi0.
* @param z Input value.
* @return Derivative of the polynomial.
*/
double xdpsi0(double z);
/**
* @brief Interpolating polynomial function xpsi1.
* @param z Input value.
* @return Result of the polynomial.
*/
double xpsi1(double z);
/**
* @brief Derivative of the interpolating polynomial function xpsi1.
* @param z Input value.
* @return Derivative of the polynomial.
*/
double xdpsi1(double z);
/**
* @brief Interpolating polynomial function h3.
* @param fi Array of coefficients.
* @param w0t Weight 0 for temperature.
* @param w1t Weight 1 for temperature.
* @param w0mt Weight 0 for temperature (minus).
* @param w1mt Weight 1 for temperature (minus).
* @param w0d Weight 0 for density.
* @param w1d Weight 1 for density.
* @param w0md Weight 0 for density (minus).
* @param w1md Weight 1 for density (minus).
* @return Result of the polynomial.
*/
double h3(double fi[36], double w0t, double w1t, double w0mt, double w1mt,
double w0d, double w1d, double w0md, double w1md);
/**
* @brief Interpolating polynomial function h5.
* @param fi Array of coefficients.
* @param w0t Weight 0 for temperature.
* @param w1t Weight 1 for temperature.
* @param w2t Weight 2 for temperature.
* @param w0mt Weight 0 for temperature (minus).
* @param w1mt Weight 1 for temperature (minus).
* @param w2mt Weight 2 for temperature (minus).
* @param w0d Weight 0 for density.
* @param w1d Weight 1 for density.
* @param w2d Weight 2 for density.
* @param w0md Weight 0 for density (minus).
* @param w1md Weight 1 for density (minus).
* @param w2md Weight 2 for density (minus).
* @return Result of the polynomial.
*/
double h5(double fi[36], double w0t, double w1t, double w2t, double w0mt,
double w1mt, double w2mt, double w0d, double w1d, double w2d,
double w0md, double w1md, double w2md);
/**
* @brief Read the Helmholtz EOS table from a file.
* @param filename Path to the file containing the table.
* @return HELMTable structure containing the table data.
*/
HELMTable read_helm_table(const std::string filename);
/**
* @brief Calculate the Helmholtz EOS components.
* @param q EOSInput parameters for the EOS calculation.
* @param table HELMTable structure containing the table data.
* @return EOS structure containing the calculated quantities.
*/
EOS get_helm_EOS(EOSInput &q, const HELMTable &table);
}
#endif // HELM_H