tboudreaux/SERiF
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
3117360d494db8db08203460f7d4990630980955
SERiF/src/eos/private/helm.cpp
Emily Boudreaux 3117360d49
Some checks failed
Build and Test / build-and-test-ubuntu (ubuntu-24.04) (push) Has been cancelled
Details
fix(build): started bringing SERiF back up to speed with 4D-STAR C++
2025-12-15 13:12:47 -05:00

850 lines
29 KiB
C++
Raw Blame History

/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Authors: Aaron Dotter, Emily Boudreaux
// Last Modified: March 20, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
// 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
#include <iostream>
#include <fstream>
#include <memory>
#include <sstream>
#include <cmath>
#include <string>
#include <stdexcept>
#include <array>
#include <numbers>
#include <cerrno>
#include "helm.h"
#include "fourdst/constants/const.h"
#include "fourdst/config/config.h"
#include "quill/LogMacros.h"
#include "quill/Logger.h"
#include "fourdst/logging/logging.h"
using namespace std;
// interpolating polynomila function definitions
namespace serif::eos::helmholtz {
double** heap_allocate_contiguous_2D_memory(const int rows, const int cols) {
if (rows <= 0 || cols <= 0) {
throw std::invalid_argument("Rows and columns must be positive integers.");
}
const auto 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;
}
double psi0(const double z) { return z*z*z * (z * (-6.0 * z + 15.0) - 10.0) + 1.0; }
double dpsi0(const double z) { return z*z * ( z * (-30.0 * z + 60.0) - 30.0); }
double ddpsi0(const double z) { return z * ( z * (-120.0 * z + 180.0) - 60.0); }
double psi1(const double z) { return z * ( z*z * ( z * (-3.0*z + 8.0) - 6.0) + 1.0); }
double dpsi1(const double z) { return z*z * ( z * (-15.0*z + 32.0) - 18.0) +1.0; }
double ddpsi1(const double z) { return z * (z * (-60.0*z + 96.0) - 36.0); }
double psi2(const double z) { return 0.5 * z*z * ( z* ( z * (-z + 3.0) - 3.0) + 1.0); }
double dpsi2(const double z) { return 0.5 * z * ( z * (z * (-5.0*z + 12.0) - 9.0) + 2.0); }
double ddpsi2(const double z) { return 0.5*(z*( z * (-20.*z + 36.) - 18.) + 2.); }
double xpsi0(const double z) { return z * z * (2.0 * z - 3.0) + 1.0; }
double xdpsi0(const double z) { return z * (6.0*z - 6.0); }
double xpsi1(const double z) { return z * ( z * (z - 2.0) + 1.0); }
double xdpsi1(const double z) { return z * (3.0 * z - 4.0) + 1.0; }
double h3(const std::array<double, 36> &fi, const double w0t, const double w1t, const double w0mt, const double w1mt,
const double w0d, const double w1d, const double w0md, const 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(const std::array<double, 36> &fi, const double w0t, const double w1t, const double w2t, const double w0mt,
const double w1mt, const double w2mt, const double w0d, const double w1d, const double w2d,
const double w0md, const double w1md, const 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
std::unique_ptr<HELMTable> read_helm_table(const std::string &filename) {
fourdst::logging::LogManager& logManager = fourdst::logging::LogManager::getInstance();
quill::Logger* logger = logManager.getLogger("log");
LOG_INFO(logger, "read_helm_table : Reading HELM table from file {}", filename);
// Make a unique pointer to the HELMTable
auto table = std::make_unique<serif::eos::helmholtz::HELMTable>();
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
***/
serif::eos::helmholtz::HELMEOSOutput get_helm_EOS(serif::eos::helmholtz::HELMEOSInput &q, const serif::eos::helmholtz::HELMTable &table) {
fourdst::logging::LogManager& logManager = fourdst::logging::LogManager::getInstance();
quill::Logger* logger = logManager.getLogger("log");
fourdst::constant::Constants& constants = fourdst::constant::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;
constexpr 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);
constexpr double qe = 4.8032042712e-10;
constexpr double esqu = qe * qe;
constexpr 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 = (log10(din)-dlo)*dstpi;
iat = max(0, min(iat, table.imax-1)); // bounds
//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
constexpr double a1 = -0.898004;
constexpr double b1 = 0.96786;
constexpr double c1 = 0.220703;
constexpr double d1 = -0.86097;
constexpr double e1 = 2.5269;
constexpr double a2 = 0.29561;
constexpr double b2 = 1.9885;
constexpr 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:
serif::eos::helmholtz::HELMEOSOutput 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.
}
}
Reference in New Issue View Git Blame Copy Permalink