Merge branch 'main' into feature/pointwisePolytrope
This commit is contained in:
434964
assets/eos/helm_table.dat
Normal file
434964
assets/eos/helm_table.dat
Normal file
File diff suppressed because it is too large
Load Diff
@@ -1,13 +1,5 @@
|
||||
patchFile = files('disable_mfem_selfcheck.patch')
|
||||
|
||||
patch_check = run_command('grep', '-q', 'MFEM_CHECK_TARGET_NAME', 'subprojects/mfem/CMakeLists.txt', check: false)
|
||||
if patch_check.returncode() == 0
|
||||
message('Patching MFEM CMakeLists.txt to remove self checks')
|
||||
run_command('patch', '-p4', '-d', '../../subprojects/mfem', '-i', patchFile[0].full_path(), check: true)
|
||||
else
|
||||
message('MFEM CMakeLists.txt already patched')
|
||||
endif
|
||||
|
||||
mfem_cmake_options = cmake.subproject_options()
|
||||
mfem_cmake_options.add_cmake_defines({
|
||||
'MFEM_ENABLE_EXAMPLES': 'OFF',
|
||||
|
||||
@@ -12,8 +12,13 @@ libconst = library('const',
|
||||
const_sources,
|
||||
include_directories: include_directories('public'),
|
||||
cpp_args: ['-fvisibility=default'],
|
||||
dependencies: [const_dep],
|
||||
dependencies: [const_data_dep],
|
||||
install : true)
|
||||
|
||||
const_dep = declare_dependency(
|
||||
include_directories: include_directories('public'),
|
||||
link_with: libconst,
|
||||
)
|
||||
|
||||
# Make headers accessible
|
||||
install_headers(const_headers, subdir : '4DSSE/const')
|
||||
@@ -5,6 +5,7 @@
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
#include <set>
|
||||
#include <map>
|
||||
|
||||
/**
|
||||
|
||||
23
src/eos/meson.build
Normal file
23
src/eos/meson.build
Normal 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')
|
||||
825
src/eos/private/helm.cpp
Normal file
825
src/eos/private/helm.cpp
Normal file
@@ -0,0 +1,825 @@
|
||||
// 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 <cerrno>
|
||||
|
||||
|
||||
#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] = pow(10, tlo + tstp*j); }
|
||||
|
||||
for (i=0; i<table.imax; i++) { table.d[i] = pow(10, dlo + dstp*i); }
|
||||
|
||||
ifstream helm_table(filename);
|
||||
if (!helm_table) {
|
||||
int errorCode = errno;
|
||||
LOG_ERROR(logger, "read_helm_table : Error ({}) opening file {}", errorCode, filename);
|
||||
throw std::runtime_error("Error (" + std::to_string(errorCode) + ") 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.
|
||||
}
|
||||
}
|
||||
366
src/eos/public/helm.h
Normal file
366
src/eos/public/helm.h
Normal file
@@ -0,0 +1,366 @@
|
||||
// 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>
|
||||
#include <format>
|
||||
|
||||
/**
|
||||
* @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: " << std::format("{0:24.16e}",eos.ye) << "\n";
|
||||
os << " Electron Chemical Potential: " << std::format("{0:24.16e}",eos.etaele) << "\n";
|
||||
os << " Electron Number Density: " << std::format("{0:24.16e}",eos.xnefer) << "\n";
|
||||
os << " Total Pressure: " << std::format("{0:24.16e}",eos.ptot) << "\n";
|
||||
os << " dPres/dRho: " << std::format("{0:24.16e}",eos.dpresdd) << "\n";
|
||||
os << " dPres/dT: " << std::format("{0:24.16e}",eos.dpresdt) << "\n";
|
||||
os << " Gas Pressure: " << std::format("{0:24.16e}",eos.pgas) << "\n";
|
||||
os << " Radiation Pressure: " << std::format("{0:24.16e}",eos.prad) << "\n";
|
||||
os << " Total Energy: " << std::format("{0:24.16e}",eos.etot) << "\n";
|
||||
os << " dEner/dRho: " << std::format("{0:24.16e}",eos.denerdd) << "\n";
|
||||
os << " dEner/dT: " << std::format("{0:24.16e}",eos.denerdt) << "\n";
|
||||
os << " Gas Energy: " << std::format("{0:24.16e}",eos.egas) << "\n";
|
||||
os << " Radiation Energy: " << std::format("{0:24.16e}",eos.erad) << "\n";
|
||||
os << " Total Entropy: " << std::format("{0:24.16e}",eos.stot) << "\n";
|
||||
os << " dEntr/dRho: " << std::format("{0:24.16e}",eos.dentrdd) << "\n";
|
||||
os << " dEntr/dT: " << std::format("{0:24.16e}",eos.dentrdt) << "\n";
|
||||
os << " Gas Entropy: " << std::format("{0:24.16e}",eos.sgas) << "\n";
|
||||
os << " Radiation Entropy: " << std::format("{0:24.16e}",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
|
||||
@@ -7,5 +7,7 @@ subdir('dobj')
|
||||
subdir('const')
|
||||
subdir('opatIO')
|
||||
subdir('meshIO')
|
||||
subdir('config')
|
||||
subdir('probe')
|
||||
subdir('eos')
|
||||
subdir('poly')
|
||||
|
||||
@@ -1,9 +1,6 @@
|
||||
data_file = files('const.dat')
|
||||
command_file = files('format.sh')
|
||||
output_file = meson.current_build_dir() + '/embedded_constants.h'
|
||||
message('Data file absolute path: ' + data_file[0].full_path())
|
||||
message('Meson source directory: ' + meson.current_source_dir())
|
||||
message('Meson build directory: ' + meson.current_build_dir())
|
||||
|
||||
embedded_constants_h = custom_target('embed_constants',
|
||||
input: data_file,
|
||||
@@ -12,9 +9,9 @@ embedded_constants_h = custom_target('embed_constants',
|
||||
)
|
||||
|
||||
# Ensure the generated header is included
|
||||
const_header = include_directories('.')
|
||||
const_data_header = include_directories('.')
|
||||
|
||||
const_dep = declare_dependency(
|
||||
include_directories: const_header,
|
||||
const_data_dep = declare_dependency(
|
||||
include_directories: const_data_header,
|
||||
sources: embedded_constants_h
|
||||
)
|
||||
16
subprojects/gtest.wrap
Normal file
16
subprojects/gtest.wrap
Normal file
@@ -0,0 +1,16 @@
|
||||
[wrap-file]
|
||||
directory = googletest-1.15.2
|
||||
source_url = https://github.com/google/googletest/archive/refs/tags/v1.15.2.tar.gz
|
||||
source_filename = gtest-1.15.2.tar.gz
|
||||
source_hash = 7b42b4d6ed48810c5362c265a17faebe90dc2373c885e5216439d37927f02926
|
||||
patch_filename = gtest_1.15.2-2_patch.zip
|
||||
patch_url = https://wrapdb.mesonbuild.com/v2/gtest_1.15.2-2/get_patch
|
||||
patch_hash = 641a16b33c96cd32a593537bc30eb7d853c5cc361fa1ee96884f0e2fca21e2d3
|
||||
source_fallback_url = https://github.com/mesonbuild/wrapdb/releases/download/gtest_1.15.2-2/gtest-1.15.2.tar.gz
|
||||
wrapdb_version = 1.15.2-2
|
||||
|
||||
[provide]
|
||||
gtest = gtest_dep
|
||||
gtest_main = gtest_main_dep
|
||||
gmock = gmock_dep
|
||||
gmock_main = gmock_main_dep
|
||||
@@ -1,5 +1,6 @@
|
||||
[wrap-git]
|
||||
url = https://github.com/mfem/mfem.git
|
||||
revision = master
|
||||
diff_files = mfem/disable_mfem_selfcheck.patch
|
||||
|
||||
[cmake]
|
||||
BIN
subprojects/packagecache/gtest-1.15.2.tar.gz
Normal file
BIN
subprojects/packagecache/gtest-1.15.2.tar.gz
Normal file
Binary file not shown.
BIN
subprojects/packagecache/gtest_1.15.2-2_patch.zip
Normal file
BIN
subprojects/packagecache/gtest_1.15.2-2_patch.zip
Normal file
Binary file not shown.
@@ -1,4 +1,4 @@
|
||||
--- subprojects/mfem/CMakeLists.txt 2025-02-12 15:54:52.454728232 -0500
|
||||
--- mfem/CMakeLists.txt 2025-02-12 15:54:52.454728232 -0500
|
||||
+++ CMakeLists.txt.bak 2025-02-12 16:08:06.654542689 -0500
|
||||
@@ -765,7 +765,7 @@
|
||||
if (MFEM_ENABLE_EXAMPLES)
|
||||
@@ -11,7 +11,7 @@ foreach test_file : test_sources
|
||||
test_exe = executable(
|
||||
exe_name,
|
||||
test_file,
|
||||
dependencies: [gtest_dep, config_dep],
|
||||
dependencies: [gtest_dep, config_dep, gtest_main],
|
||||
include_directories: include_directories('../../src/config/public'),
|
||||
link_with: libconst, # Link the dobj library
|
||||
install_rpath: '@loader_path/../../src' # Ensure runtime library path resolves correctly
|
||||
|
||||
@@ -11,7 +11,7 @@ foreach test_file : test_sources
|
||||
test_exe = executable(
|
||||
exe_name,
|
||||
test_file,
|
||||
dependencies: [gtest_dep, const_dep],
|
||||
dependencies: [gtest_dep, const_dep, gtest_main],
|
||||
include_directories: include_directories('../../src/const/public'),
|
||||
link_with: libconst, # Link the dobj library
|
||||
install_rpath: '@loader_path/../../src' # Ensure runtime library path resolves correctly
|
||||
|
||||
@@ -13,7 +13,7 @@ foreach test_file : test_sources
|
||||
test_exe = executable(
|
||||
exe_name,
|
||||
test_file,
|
||||
dependencies: gtest_dep,
|
||||
dependencies: [gtest_dep, gtest_main],
|
||||
include_directories: include_directories('../../src/dobj/public'),
|
||||
link_with: libdobj, # Link the dobj library
|
||||
install_rpath: '@loader_path/../../src' # Ensure runtime library path resolves correctly
|
||||
|
||||
85
tests/eos/eosTest.cpp
Normal file
85
tests/eos/eosTest.cpp
Normal file
@@ -0,0 +1,85 @@
|
||||
#include <gtest/gtest.h>
|
||||
#include <iostream>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
#include <set>
|
||||
#include <sstream>
|
||||
|
||||
#include "helm.h"
|
||||
|
||||
/**
|
||||
* @file constTest.cpp
|
||||
* @brief Unit tests for the const class.
|
||||
*/
|
||||
|
||||
/**
|
||||
* @brief Test suite for the const class.
|
||||
*/
|
||||
class eosTest : public ::testing::Test {};
|
||||
|
||||
std::string HELM_FILENAME = std::string(getenv("MESON_SOURCE_ROOT")) + "/assets/eos/helm_table.dat";
|
||||
/**
|
||||
* @test Verify default constructor initializes correctly.
|
||||
*/
|
||||
TEST_F(eosTest, constructor) {
|
||||
using namespace helmholtz;
|
||||
EXPECT_NO_THROW(HELMTable table = read_helm_table(HELM_FILENAME));
|
||||
}
|
||||
|
||||
TEST_F(eosTest, read_helm_table) {
|
||||
using namespace helmholtz;
|
||||
HELMTable table = read_helm_table(HELM_FILENAME);
|
||||
// Capture the << operator output as a string
|
||||
std::stringstream ss;
|
||||
ss << table;
|
||||
EXPECT_EQ(ss.str(), "HELMTable Data:\n imax: 541, jmax: 201\n Temperature Range: [1000, 1e+13]\n Density Range: [1e-12, 1e+15]\n");
|
||||
}
|
||||
|
||||
TEST_F(eosTest, get_helm_EOS) {
|
||||
using namespace helmholtz;
|
||||
|
||||
const int nel=3;
|
||||
double xmass[nel], aion[nel], zion[nel];
|
||||
EOSInput 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;
|
||||
|
||||
double asum = 0.0;
|
||||
double zsum = 0.0;
|
||||
for (int 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;
|
||||
|
||||
HELMTable table = read_helm_table(HELM_FILENAME);
|
||||
EOS eos = get_helm_EOS(eos1, table);
|
||||
// std::cout << eos << std::endl;
|
||||
|
||||
//Check composition info
|
||||
EXPECT_DOUBLE_EQ( eos.ye, 8.75e-01);
|
||||
|
||||
//Check E, P, S and derivatives of each wrt Rho and T
|
||||
EXPECT_DOUBLE_EQ( eos.etaele, 2.3043348231021554e+01);
|
||||
EXPECT_DOUBLE_EQ( eos.etot, 1.1586558190936826e+17);
|
||||
EXPECT_DOUBLE_EQ(eos.denerdd, 6.1893000468285858e+10);
|
||||
EXPECT_DOUBLE_EQ(eos.denerdt, 1.2129708972542575e+08);
|
||||
EXPECT_DOUBLE_EQ( eos.ptot, 6.9610135220017030e+22);
|
||||
EXPECT_DOUBLE_EQ(eos.dpresdd, 1.0296440482849070e+17);
|
||||
EXPECT_DOUBLE_EQ(eos.dpresdt, 7.7171347517311625e+13);
|
||||
EXPECT_DOUBLE_EQ( eos.stot, 6.0647461970567346e+08);
|
||||
EXPECT_DOUBLE_EQ(eos.dentrdd,-7.7171347517311645e+01);
|
||||
EXPECT_DOUBLE_EQ(eos.dentrdt, 1.2129708972542577e+00);
|
||||
|
||||
const double abs_err = 1.0e-12;
|
||||
// Maxwell relations, should always be zero
|
||||
EXPECT_NEAR( eos.dse, 0, abs_err);
|
||||
EXPECT_NEAR( eos.dpe, 0, abs_err);
|
||||
EXPECT_NEAR( eos.dsp, 0, abs_err);
|
||||
}
|
||||
23
tests/eos/meson.build
Normal file
23
tests/eos/meson.build
Normal file
@@ -0,0 +1,23 @@
|
||||
# Test files for const
|
||||
test_sources = [
|
||||
'eosTest.cpp',
|
||||
]
|
||||
|
||||
foreach test_file : test_sources
|
||||
exe_name = test_file.split('.')[0]
|
||||
message('Building test: ' + exe_name)
|
||||
|
||||
# Create an executable target for each test
|
||||
test_exe = executable(
|
||||
exe_name,
|
||||
test_file,
|
||||
dependencies: [gtest_dep, eos_dep, gtest_main],
|
||||
install_rpath: '@loader_path/../../src' # Ensure runtime library path resolves correctly
|
||||
)
|
||||
|
||||
# Add the executable as a test
|
||||
test(
|
||||
exe_name,
|
||||
test_exe,
|
||||
env: ['MESON_SOURCE_ROOT=' + meson.project_source_root(), 'MESON_BUILD_ROOT=' + meson.project_build_root()])
|
||||
endforeach
|
||||
@@ -11,7 +11,7 @@ foreach test_file : test_sources
|
||||
test_exe = executable(
|
||||
exe_name,
|
||||
test_file,
|
||||
dependencies: [gtest_dep, mfem_dep],
|
||||
dependencies: [gtest_dep, mfem_dep, gtest_main],
|
||||
include_directories: include_directories('../../src/meshIO/public'),
|
||||
link_with: libmeshIO, # Link the dobj library
|
||||
install_rpath: '@loader_path/../../src' # Ensure runtime library path resolves correctly
|
||||
|
||||
@@ -1,5 +1,6 @@
|
||||
# Google Test dependency
|
||||
gtest_dep = dependency('gtest', main: true, required : true)
|
||||
gtest_main = dependency('gtest_main', required: true)
|
||||
gtest_nomain_dep = dependency('gtest', main: false, required : true)
|
||||
|
||||
# Subdirectories for unit and integration tests
|
||||
@@ -8,8 +9,9 @@ subdir('const')
|
||||
subdir('opatIO')
|
||||
subdir('meshIO')
|
||||
subdir('probe')
|
||||
subdir('poly')
|
||||
subdir('config')
|
||||
subdir('eos')
|
||||
subdir('poly')
|
||||
|
||||
# Subdirectories for sandbox tests
|
||||
subdir('dobj_sandbox')
|
||||
|
||||
@@ -13,7 +13,7 @@ foreach test_file : test_sources
|
||||
test_exe = executable(
|
||||
exe_name,
|
||||
test_file,
|
||||
dependencies: [gtest_dep, picosha2_dep],
|
||||
dependencies: [gtest_dep, picosha2_dep, gtest_main],
|
||||
include_directories: include_directories('../../src/opatIO/public'),
|
||||
link_with: libopatIO, # Link the dobj library
|
||||
install_rpath: '@loader_path/../../src' # Ensure runtime library path resolves correctly
|
||||
|
||||
@@ -11,7 +11,7 @@ foreach test_file : test_sources
|
||||
test_exe = executable(
|
||||
exe_name,
|
||||
test_file,
|
||||
dependencies: [gtest_dep, probe_dep, mfem_dep, quill_dep],
|
||||
dependencies: [gtest_dep, probe_dep, mfem_dep, quill_dep, gtest_main],
|
||||
install_rpath: '@loader_path/../../src' # Ensure runtime library path resolves correctly
|
||||
)
|
||||
|
||||
|
||||
Reference in New Issue
Block a user