SERiF 0.0.1a
3+1D Stellar Structure and Evolution
Loading...
Searching...
No Matches
helm.cpp
Go to the documentation of this file.
1/* ***********************************************************************
2//
3// Copyright (C) 2025 -- The 4D-STAR Collaboration
4// File Authors: Aaron Dotter, Emily Boudreaux
5// Last Modified: March 20, 2025
6//
7// 4DSSE is free software; you can use it and/or modify
8// it under the terms and restrictions the GNU General Library Public
9// License version 3 (GPLv3) as published by the Free Software Foundation.
10//
11// 4DSSE is distributed in the hope that it will be useful,
12// but WITHOUT ANY WARRANTY; without even the implied warranty of
13// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
14// See the GNU Library General Public License for more details.
15//
16// You should have received a copy of the GNU Library General Public License
17// along with this software; if not, write to the Free Software
18// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19//
20// *********************************************************************** */
21// this file contains a C++ conversion of Frank Timmes' fortran code
22// helmholtz.f90, which implements the Helmholtz Free Energy EOS described
23// by Timmes & Swesty (2000) doi:10.1086/313304
24
25#include <iostream>
26#include <fstream>
27#include <memory>
28#include <sstream>
29#include <cmath>
30#include <string>
31#include <stdexcept>
32#include <array>
33#include <numbers>
34#include <cerrno>
35
36
37
38#include "helm.h"
39#include "const.h"
40#include "probe.h"
41#include "config.h"
42#include "quill/LogMacros.h"
43
44namespace serif::constant {
45 class Constants;
46}
47
48using namespace std;
49
50// interpolating polynomila function definitions
51namespace serif::eos::helmholtz {
52 double** heap_allocate_contiguous_2D_memory(const int rows, const int cols) {
53 if (rows <= 0 || cols <= 0) {
54 throw std::invalid_argument("Rows and columns must be positive integers.");
55 }
56 const auto array = new double*[rows];
57 array[0] = new double[rows * cols];
58
59 for (int i = 1; i < rows; i++) {
60 array[i] = array[0] + i * cols;
61 }
62
63 return array;
64 }
65
66 void heap_deallocate_contiguous_2D_memory(double **array) {
67 delete[] array[0];
68 delete[] array;
69 }
70 double psi0(const double z) { return z*z*z * (z * (-6.0 * z + 15.0) - 10.0) + 1.0; }
71
72 double dpsi0(const double z) { return z*z * ( z * (-30.0 * z + 60.0) - 30.0); }
73
74 double ddpsi0(const double z) { return z * ( z * (-120.0 * z + 180.0) - 60.0); }
75
76 double psi1(const double z) { return z * ( z*z * ( z * (-3.0*z + 8.0) - 6.0) + 1.0); }
77
78 double dpsi1(const double z) { return z*z * ( z * (-15.0*z + 32.0) - 18.0) +1.0; }
79
80 double ddpsi1(const double z) { return z * (z * (-60.0*z + 96.0) - 36.0); }
81
82 double psi2(const double z) { return 0.5 * z*z * ( z* ( z * (-z + 3.0) - 3.0) + 1.0); }
83
84 double dpsi2(const double z) { return 0.5 * z * ( z * (z * (-5.0*z + 12.0) - 9.0) + 2.0); }
85
86 double ddpsi2(const double z) { return 0.5*(z*( z * (-20.*z + 36.) - 18.) + 2.); }
87
88 double xpsi0(const double z) { return z * z * (2.0 * z - 3.0) + 1.0; }
89
90 double xdpsi0(const double z) { return z * (6.0*z - 6.0); }
91
92 double xpsi1(const double z) { return z * ( z * (z - 2.0) + 1.0); }
93
94 double xdpsi1(const double z) { return z * (3.0 * z - 4.0) + 1.0; }
95
96 double h3(const std::array<double, 36> &fi, const double w0t, const double w1t, const double w0mt, const double w1mt,
97 const double w0d, const double w1d, const double w0md, const double w1md) {
98 return fi[0] * w0d * w0t + fi[1] * w0md * w0t
99 + fi[2] * w0d * w0mt + fi[3] * w0md * w0mt
100 + fi[4] * w0d * w1t + fi[5] * w0md * w1t
101 + fi[6] * w0d * w1mt + fi[7] * w0md * w1mt
102 + fi[8] * w1d * w0t + fi[9] * w1md * w0t
103 + fi[10] * w1d * w0mt + fi[11] * w1md * w0mt
104 + fi[12] * w1d * w1t + fi[13] * w1md * w1t
105 + fi[14] * w1d * w1mt + fi[15] * w1md * w1mt;
106 }
107
108 double h5(const std::array<double, 36> &fi, const double w0t, const double w1t, const double w2t, const double w0mt,
109 const double w1mt, const double w2mt, const double w0d, const double w1d, const double w2d,
110 const double w0md, const double w1md, const double w2md) {
111 return fi[0] * w0d * w0t + fi[1] * w0md * w0t
112 + fi[2] * w0d * w0mt + fi[3] * w0md * w0mt
113 + fi[4] * w0d * w1t + fi[5] * w0md * w1t
114 + fi[6] * w0d * w1mt + fi[7] * w0md * w1mt
115 + fi[8] * w0d * w2t + fi[ 9] * w0md * w2t
116 + fi[10] * w0d * w2mt + fi[11] * w0md * w2mt
117 + fi[12] * w1d * w0t + fi[13] * w1md * w0t
118 + fi[14] * w1d * w0mt + fi[15] * w1md * w0mt
119 + fi[16] * w2d * w0t + fi[17] * w2md * w0t
120 + fi[18] * w2d * w0mt + fi[19] * w2md * w0mt
121 + fi[20] * w1d * w1t + fi[21] * w1md * w1t
122 + fi[22] * w1d * w1mt + fi[23] * w1md * w1mt
123 + fi[24] * w2d * w1t + fi[25] * w2md * w1t
124 + fi[26] * w2d * w1mt + fi[27] * w2md * w1mt
125 + fi[28] * w1d * w2t + fi[29] * w1md * w2t
126 + fi[30] * w1d * w2mt + fi[31] * w1md * w2mt
127 + fi[32] * w2d * w2t + fi[33] * w2md * w2t
128 + fi[34] * w2d * w2mt + fi[35] * w2md * w2mt;
129 }
130
131 // this function reads in the HELM table and stores in the above arrays
132 std::unique_ptr<HELMTable> read_helm_table(const std::string &filename) {
133 serif::config::Config& config = serif::config::Config::getInstance();
134 auto logFile = config.get<std::string>("EOS:Helm:LogFile", "log");
135 serif::probe::LogManager& logManager = serif::probe::LogManager::getInstance();
136 quill::Logger* logger = logManager.getLogger(logFile);
137 LOG_INFO(logger, "read_helm_table : Reading HELM table from file {}", filename);
138
139 // Make a unique pointer to the HELMTable
140 auto table = std::make_unique<serif::eos::helmholtz::HELMTable>();
141 string data;
142 int i, j;
143
144 //set T and Rho (d) arrays
145 for (j=0; j<table->jmax; j++) { table->t[j] = pow(10, tlo + tstp*j); }
146
147 for (i=0; i<table->imax; i++) { table->d[i] = pow(10, dlo + dstp*i); }
148
149 ifstream helm_table(filename);
150 if (!helm_table) {
151 int errorCode = errno;
152 LOG_ERROR(logger, "read_helm_table : Error ({}) opening file {}", errorCode, filename);
153 throw std::runtime_error("Error (" + std::to_string(errorCode) + ") opening file " + filename);
154 }
155 //read the Helmholtz free energy and its derivatives
156 for (j=0; j<table->jmax; j++) {
157 for (i=0; i<table->imax; i++){
158 getline(helm_table, data);
159 stringstream id(data);
160 id >> table->f[i][j];
161 id >> table->fd[i][j];
162 id >> table->ft[i][j];
163 id >> table->fdd[i][j];
164 id >> table->ftt[i][j];
165 id >> table->fdt[i][j];
166 id >> table->fddt[i][j];
167 id >> table->fdtt[i][j];
168 id >> table->fddtt[i][j];
169 }
170 }
171
172 //read the pressure derivative with density
173 for (j=0; j<table->jmax; j++) {
174 for (i=0; i<table->imax; i++){
175 getline(helm_table, data);
176 stringstream id(data);
177 id >> table->dpdf[i][j];
178 id >> table->dpdfd[i][j];
179 id >> table->dpdft[i][j];
180 id >> table->dpdfdt[i][j];
181 }
182 }
183
184 //read the electron chemical potential
185 for (j=0; j<table->jmax; j++) {
186 for (i=0; i<table->imax; i++){
187 getline(helm_table, data);
188 stringstream id(data);
189 id >> table->ef[i][j];
190 id >> table->efd[i][j];
191 id >> table->eft[i][j];
192 id >> table->efdt[i][j];
193 }
194 }
195
196 //read the number density
197 for (j=0; j<table->jmax; j++) {
198 for (i=0; i<table->imax; i++){
199 getline(helm_table, data);
200 stringstream id(data);
201 id >> table->xf[i][j];
202 id >> table->xfd[i][j];
203 id >> table->xft[i][j];
204 id >> table->xfdt[i][j];
205 }
206 }
207
208 helm_table.close(); //done reading
209
210 // construct the temperature and density deltas and their inverses
211 for (j=0; j<table->jmax; j++) {
212 double dth = table->t[j+1] - table->t[j];
213 double dt2 = dth * dth;
214 double dti = 1.0/dth;
215 double dt2i = 1.0/dt2;
216 table->dt_sav[j] = dth;
217 table->dt2_sav[j] = dt2;
218 table->dti_sav[j] = dti;
219 table->dt2i_sav[j] = dt2i;
220 }
221
222
223 for (i=0; i<table->imax; i++) {
224 double dd = table->d[i+1] - table->d[i];
225 double dd2 = dd * dd;
226 double ddi = 1.0/dd;
227 table->dd_sav[i] = dd;
228 table->dd2_sav[i] = dd2;
229 table->ddi_sav[i] = ddi;
230 }
231
232 table->loaded = true;
233
234 return table;
235 }
236
237
238 /***
239 this function evaluates the EOS components:
240 ion, radiation, electron-positron and Coulomb interaction
241 and returns the calculated quantities in the input
242 ***/
243 serif::eos::helmholtz::HELMEOSOutput get_helm_EOS(serif::eos::helmholtz::HELMEOSInput &q, const serif::eos::helmholtz::HELMTable &table) {
244 serif::config::Config& config = serif::config::Config::getInstance();
245 auto logFile = config.get<std::string>("EOS:Helm:LogFile", "log");
246 serif::probe::LogManager& logManager = serif::probe::LogManager::getInstance();
247 quill::Logger* logger = logManager.getLogger(logFile);
248
249 serif::constant::Constants& constants = serif::constant::Constants::getInstance();
250 const double pi = std::numbers::pi;
251 const double amu = constants.get("u").value;
252 const double h = constants.get("h").value;
253 const double avo = constants.get("N_a").value;
254 const double kerg = constants.get("kB").value;
255 const double clight = constants.get("c").value;
256
257 const double kergavo = kerg*avo;
258 const double clight2 = clight * clight;
259 constexpr double ssol = 5.6704e-5;
260 const double asol = 4 * ssol / clight;
261 const double asoli3 = asol / 3;
262 const double sioncon = (2.0 * pi * amu * kerg)/(h*h);
263 constexpr double qe = 4.8032042712e-10;
264 constexpr double esqu = qe * qe;
265 constexpr double third = 1./3.;
266 std::array<double, 36> fi;
267 double T, rho, s, free, abar, zbar, ytot1, ye, ptot;
268 double prad, srad, erad, etot, stot;
269 double dpraddd, dpraddt, dpradda, dpraddz;
270 double deraddd, deraddt, deradda, deraddz;
271 double dsraddd, dsraddt, dsradda, dsraddz;
272 double dpiondd, dpiondt, dpionda, dpiondz;
273 double deiondd, deiondt, deionda, deiondz;
274 double dsiondd, dsiondt, dsionda, dsiondz;
275 double dpcouldd, dpcouldt, dpcoulda, dpcouldz;
276 double decouldd, decouldt, decoulda, decouldz;
277 double dscouldd, dscouldt, dscoulda, dscouldz;
278 double dpgasdd, dpgasdt, dpgasda, dpgasdz;
279 double degasdd, degasdt, degasda, degasdz;
280 double dsgasdd, dsgasdt, dsgasda, dsgasdz;
281 double dpepdd, dpepdt, dpepda, dpepdz;
282 double dsepdd, dsepdt, dsepda, dsepdz;
283 double deepdd, deepdt, deepda, deepdz;
284 double dxnidd, dxnida;
285 double dpresdd, dpresdt, dpresda, dpresdz;
286 double denerdd, denerdt, denerda, denerdz;
287 double dentrdd, dentrdt, dentrda, dentrdz;
288 double xni, sion, pion, eion, x, y, z;
289 double pgas, egas, sgas, pele, sele, eele, pcoul, scoul, ecoul;
290 int jat, iat;
291
292 T = q.T; rho = q.rho; abar = q.abar; zbar = q.zbar;
293
294 ytot1 = 1/abar;
295 ye = zbar * ytot1;
296
297 double rhoi = 1/rho;
298 double Ti = 1/T;
299 double kT= kerg*T;
300 double kTinv = 1/kT;
301
302 /*** radiation ***/
303 prad = asoli3 * pow(T,4);
304 dpraddd = 0.;
305 dpraddt = 4 * prad * Ti;
306 dpradda = 0.;
307 dpraddz = 0.;
308
309 erad = 3 * prad * rhoi;
310 deraddd = -erad * rhoi;
311 deraddt = 3 * dpraddt * rhoi;
312 deradda = 0.;
313 deraddz = 0.;
314
315 srad = (prad*rhoi + erad)*Ti;
316 dsraddd = (dpraddd*rhoi - prad*rhoi*rhoi + deraddd)*Ti;
317 dsraddt = (dpraddt*rhoi + deraddt - srad)*Ti;
318 dsradda = 0.;
319 dsraddz = 0.;
320
321
322 /*** ions ***/
323 xni = avo * ytot1 * rho;
324 dxnidd = avo * ytot1;
325 dxnida = -xni * ytot1;
326
327 pion = xni * kT;
328 dpiondd = dxnidd * kT;
329 dpiondt = xni * kerg;
330 dpionda = dxnida * kT;
331 dpiondz = 0.;
332
333 eion = 1.5 * pion * rhoi;
334 deiondd = (1.5 * dpiondd - eion)*rhoi;
335 deiondt = 1.5 * dpiondt*rhoi;
336 deionda = 1.5 * dpionda*rhoi;
337 deiondz = 0.;
338
339 x = abar * abar * sqrt(abar) * rhoi/avo;
340 s = sioncon * T;
341 z = x * s * sqrt(s);
342 y = log(z);
343 sion = (pion * rhoi + eion) * Ti + kergavo * ytot1 * y;
344
345 dsiondd = (dpiondd*rhoi - pion*rhoi*rhoi + deiondd)*Ti
346 - kergavo * rhoi * ytot1;
347 dsiondt = (dpiondt*rhoi + deiondt)*Ti - (pion*rhoi + eion) * Ti*Ti
348 + 1.5 * kergavo * Ti*ytot1;
349 x = avo*kerg/abar;
350 dsionda = (dpionda*rhoi + deionda)*Ti + kergavo*ytot1*ytot1* (2.5 - y);
351 dsiondz = 0.;
352
353
354 // electron-positron
355 // double xnem = xni * zbar;
356 double din = ye*rho;
357
358 // hash the table location in t, rho:
359 jat = (log10(T) - tlo)*tstpi;
360 jat = max(0, min(jat, table.jmax-1)); // bounds
361 iat = (log10(din)-dlo)*dstpi;
362 iat = max(0, min(iat, table.imax-1)); // bounds
363
364 //using the indices, access the table:
365 fi[0] = table.f[iat][jat];
366 fi[1] = table.f[iat+1][jat];
367 fi[2] = table.f[iat][jat+1];
368 fi[3] = table.f[iat+1][jat+1];
369 fi[4] = table.ft[iat][jat];
370 fi[5] = table.ft[iat+1][jat];
371 fi[6] = table.ft[iat][jat+1];
372 fi[7] = table.ft[iat+1][jat+1];
373 fi[8] = table.ftt[iat][jat];
374 fi[9] = table.ftt[iat+1][jat];
375 fi[10] = table.ftt[iat][jat+1];
376 fi[11] = table.ftt[iat+1][jat+1];
377 fi[12] = table.fd[iat][jat];
378 fi[13] = table.fd[iat+1][jat];
379 fi[14] = table.fd[iat][jat+1];
380 fi[15] = table.fd[iat+1][jat+1];
381 fi[16] = table.fdd[iat][jat];
382 fi[17] = table.fdd[iat+1][jat];
383 fi[18] = table.fdd[iat][jat+1];
384 fi[19] = table.fdd[iat+1][jat+1];
385 fi[20] = table.fdt[iat][jat];
386 fi[21] = table.fdt[iat+1][jat];
387 fi[22] = table.fdt[iat][jat+1];
388 fi[23] = table.fdt[iat+1][jat+1];
389 fi[24] = table.fddt[iat][jat];
390 fi[25] = table.fddt[iat+1][jat];
391 fi[26] = table.fddt[iat][jat+1];
392 fi[27] = table.fddt[iat+1][jat+1];
393 fi[28] = table.fdtt[iat][jat];
394 fi[29] = table.fdtt[iat+1][jat];
395 fi[30] = table.fdtt[iat][jat+1];
396 fi[31] = table.fdtt[iat+1][jat+1];
397 fi[32] = table.fddtt[iat][jat];
398 fi[33] = table.fddtt[iat+1][jat];
399 fi[34] = table.fddtt[iat][jat+1];
400 fi[35] = table.fddtt[iat+1][jat+1];
401
402 double xt = (T - table.t[jat]) * table.dti_sav[jat];
403 double xd = (din - table.d[iat]) * table.ddi_sav[iat];
404 double mxt = 1-xt;
405 double mxd = 1-xd;
406
407 //density and temperature basis functions
408 double si0t = psi0(xt);
409 double si1t = psi1(xt)*table.dt_sav[jat];
410 double si2t = psi2(xt)*table.dt2_sav[jat];
411
412 double si0mt = psi0(mxt);
413 double si1mt = -psi1(mxt)*table.dt_sav[jat];
414 double si2mt = psi2(mxt)*table.dt2_sav[jat];
415
416 double si0d = psi0(xd);
417 double si1d = psi1(xd)*table.dd_sav[iat];
418 double si2d = psi2(xd)*table.dd2_sav[iat];
419
420 double si0md = psi0(mxd);
421 double si1md = -psi1(mxd)*table.dd_sav[iat];
422 double si2md = psi2(mxd)*table.dd2_sav[iat];
423
424 // derivatives of the weight functions
425 double dsi0t = dpsi0(xt)*table.dti_sav[jat];
426 double dsi1t = dpsi1(xt);
427 double dsi2t = dpsi2(xt)*table.dt_sav[jat];
428
429 double dsi0mt = -dpsi0(mxt)*table.dti_sav[jat];
430 double dsi1mt = dpsi1(mxt);
431 double dsi2mt = -dpsi2(mxt)*table.dt_sav[jat];
432
433 double dsi0d = dpsi0(xd)*table.ddi_sav[iat];
434 double dsi1d = dpsi1(xd);
435 double dsi2d = dpsi2(xd)*table.dd_sav[iat];
436
437 double dsi0md = -dpsi0(mxd)*table.ddi_sav[iat];
438 double dsi1md = dpsi1(mxd);
439 double dsi2md = -dpsi2(mxd)*table.dd_sav[iat];
440
441 // second derivatives of the weight functions
442 double ddsi0t = ddpsi0(xt)*table.dt2i_sav[jat];
443 double ddsi1t = ddpsi1(xt)*table.dti_sav[jat];
444 double ddsi2t = ddpsi2(xt);
445
446 double ddsi0mt = ddpsi0(mxt)*table.dt2i_sav[jat];
447 double ddsi1mt = -ddpsi1(mxt)*table.dti_sav[jat];
448 double ddsi2mt = ddpsi2(mxt);
449
450 // Refactor these to use LOG_DEBUG(logger, message)
451 LOG_DEBUG(logger, " xt, mxt = {} {}", xt, mxt);
452 LOG_DEBUG(logger, " dti_sav[jat] = {}", table.dti_sav[jat]);
453 LOG_DEBUG(logger, " dt2i_sav[jat] = {}", table.dt2i_sav[jat]);
454 LOG_DEBUG(logger, " ddpsi0(xt) = {}", ddpsi0(xt));
455 LOG_DEBUG(logger, " ddpsi1(xt) = {}", ddpsi1(xt));
456 LOG_DEBUG(logger, " ddpsi2(xt) = {}", ddpsi2(xt));
457 LOG_DEBUG(logger, " ddpsi0(mxt) = {}", ddpsi0(mxt));
458 LOG_DEBUG(logger, " ddpsi1(mxt) = {}", ddpsi1(mxt));
459 LOG_DEBUG(logger, " ddpsi2(mxt) = {}", ddpsi2(mxt));
460 LOG_DEBUG(logger, " ddsi0t = {}", ddsi0t);
461 LOG_DEBUG(logger, " ddsi1t = {}", ddsi1t);
462 LOG_DEBUG(logger, " ddsi2t = {}", ddsi2t);
463 LOG_DEBUG(logger, " ddsi0mt = {}", ddsi0mt);
464 LOG_DEBUG(logger, " ddsi1mt = {}", ddsi1mt);
465 LOG_DEBUG(logger, " ddsi2mt = {}", ddsi2mt);
466
467
468 ddsi0t = 0.0; ddsi1t = 0.0; ddsi2t = 1.0;
469 ddsi0mt = 0.0; ddsi1mt = 0.0; ddsi2mt = 0.0;
470
471 // free energy
472 free = h5(fi,si0t, si1t, si2t, si0mt, si1mt, si2mt,
473 si0d, si1d, si2d, si0md, si1md, si2md);
474
475 // derivative with respect to density
476 double df_d = h5(fi, si0t, si1t, si2t, si0mt, si1mt, si2mt,
477 dsi0d, dsi1d, dsi2d, dsi0md, dsi1md, dsi2md);
478
479 // derivative with respect to temperature
480 double df_t = h5(fi, dsi0t, dsi1t, dsi2t, dsi0mt, dsi1mt, dsi2mt,
481 si0d, si1d, si2d, si0md, si1md, si2md);
482
483 // second derivative with respect to temperature
484 double df_tt = h5(fi, ddsi0t, ddsi1t, ddsi2t, ddsi0mt, ddsi1mt, ddsi2mt,
485 si0d, si1d, si2d, si0md, si1md, si2md);
486
487 LOG_DEBUG(logger, "df_tt = {}", df_tt);
488
489
490 // derivative with respect to temperature and density
491 double df_dt = h5(fi, dsi0t, dsi1t, dsi2t, dsi0mt, dsi1mt, dsi2mt,
492 dsi0d, dsi1d, dsi2d, dsi0md, dsi1md, dsi2md);
493
494
495 // pressure derivatives
496 si0t = xpsi0(xt);
497 si1t = xpsi1(xt)*table.dt_sav[jat];
498
499 si0mt = xpsi0(mxt);
500 si1mt = -xpsi1(mxt)*table.dt_sav[jat];
501
502 si0d = xpsi0(xd);
503 si1d = xpsi1(xd)*table.dd_sav[iat];
504
505 si0md = xpsi0(mxd);
506 si1md = -xpsi1(mxd)*table.dd_sav[iat];
507
508
509 // derivatives of weight functions
510 dsi0t = xdpsi0(xt)*table.dti_sav[jat];
511 dsi1t = xdpsi1(xt);
512
513 dsi0mt = -xdpsi0(mxt)*table.dti_sav[jat];
514 dsi1mt = xdpsi1(mxt);
515
516 dsi0d = xdpsi0(xd)*table.ddi_sav[iat];
517 dsi1d = xdpsi1(xd);
518
519 dsi0md = -xdpsi0(mxd)*table.ddi_sav[iat];
520 dsi1md = xdpsi1(mxd);
521
522 // pressure derivative
523 fi[0] = table.dpdf[iat][jat];
524 fi[1] = table.dpdf[iat+1][jat];
525 fi[2] = table.dpdf[iat][jat+1];
526 fi[3] = table.dpdf[iat+1][jat+1];
527 fi[4] = table.dpdft[iat][jat];
528 fi[5] = table.dpdft[iat+1][jat];
529 fi[6] = table.dpdft[iat][jat+1];
530 fi[7] = table.dpdft[iat+1][jat+1];
531 fi[8] = table.dpdfd[iat][jat];
532 fi[9] = table.dpdfd[iat+1][jat];
533 fi[10] = table.dpdfd[iat][jat+1];
534 fi[11] = table.dpdfd[iat+1][jat+1];
535 fi[12] = table.dpdfdt[iat][jat];
536 fi[13] = table.dpdfdt[iat+1][jat];
537 fi[14] = table.dpdfdt[iat][jat+1];
538 fi[15] = table.dpdfdt[iat+1][jat+1];
539
540
541 // pressure derivative with density
542 dpepdd = h3(fi, si0t, si1t, si0mt, si1mt, si0d, si1d, si0md, si1md);
543 dpepdd = max(ye * dpepdd,1e-30);
544
545
546 // look in the electron chemical potential table only once
547 fi[0] = table.ef[iat][jat];
548 fi[1] = table.ef[iat+1][jat];
549 fi[2] = table.ef[iat][jat+1];
550 fi[3] = table.ef[iat+1][jat+1];
551 fi[4] = table.eft[iat][jat];
552 fi[5] = table.eft[iat+1][jat];
553 fi[6] = table.eft[iat][jat+1];
554 fi[7] = table.eft[iat+1][jat+1];
555 fi[8] = table.efd[iat][jat];
556 fi[9] = table.efd[iat+1][jat];
557 fi[10] = table.efd[iat][jat+1];
558 fi[11] = table.efd[iat+1][jat+1];
559 fi[12] = table.efdt[iat][jat];
560 fi[13] = table.efdt[iat+1][jat];
561 fi[14] = table.efdt[iat][jat+1];
562 fi[15] = table.efdt[iat+1][jat+1];
563
564 // electron chemical potential etaele
565 double etaele = h3(fi, si0t, si1t, si0mt, si1mt, si0d, si1d, si0md, si1md);
566
567 // derivative with respect to density
568 x = h3(fi, si0t, si1t, si0mt, si1mt, dsi0d, dsi1d, dsi0md, dsi1md);
569 // double detadd = ye * x;
570
571 // derivative with respect to temperature
572 // double detadt = h3(fi, dsi0t, dsi1t, dsi0mt, dsi1mt, si0d,
573 // si1d, si0md, si1md);
574
575 // derivative with respect to abar and zbar
576 // double detada = -x * din * ytot1;
577 // double detadz = x * rho * ytot1;
578
579 // look in the number density table only once
580 fi[0] = table.xf[iat][jat];
581 fi[1] = table.xf[iat+1][jat];
582 fi[2] = table.xf[iat][jat+1];
583 fi[3] = table.xf[iat+1][jat+1];
584 fi[4] = table.xft[iat][jat];
585 fi[5] = table.xft[iat+1][jat];
586 fi[6] = table.xft[iat][jat+1];
587 fi[7] = table.xft[iat+1][jat+1];
588 fi[8] = table.xfd[iat][jat];
589 fi[9] = table.xfd[iat+1][jat];
590 fi[10] = table.xfd[iat][jat+1];
591 fi[11] = table.xfd[iat+1][jat+1];
592 fi[12] = table.xfdt[iat][jat];
593 fi[13] = table.xfdt[iat+1][jat];
594 fi[14] = table.xfdt[iat][jat+1];
595 fi[15] = table.xfdt[iat+1][jat+1];
596
597 // electron chemical potential etaele
598 double xnefer = h3(fi, si0t, si1t, si0mt, si1mt, si0d, si1d, si0md, si1md);
599
600 // derivative with respect to density
601 x = h3(fi, si0t, si1t, si0mt, si1mt, dsi0d, dsi1d, dsi0md, dsi1md);
602 x = max(x,1e-30);
603 // double dxnedd = ye * x;
604
605 // derivative with respect to temperature
606 // double dxnedt = h3(fi, dsi0t, dsi1t, dsi0mt, dsi1mt, si0d, si1d, si0md, si1md);
607
608 // derivative with respect to abar and zbar
609 // double dxneda = -x * din * ytot1;
610 // double dxnedz = x * rho * ytot1;
611
612
613 // the desired electron-positron thermodynamic quantities
614 x = din * din;
615 pele = x * df_d;
616 dpepdt = x * df_dt;
617 //dpepdd is set above
618 s = dpepdd/ye - 2 * din * df_d;
619 dpepda = -ytot1 * (2 * pele + s * din);
620 dpepdz = rho * ytot1 * (2 * din * df_d + s);
621
622 x = ye * ye;
623 sele = -df_t * ye;
624 dsepdt = -df_tt * ye;
625 dsepdd = -df_dt * x;
626 dsepda = ytot1 * (ye * df_dt * din - sele);
627 dsepdz = -ytot1 * (ye * df_dt * rho + df_t);
628
629 eele = ye*free + T * sele;
630 deepdt = T * dsepdt;
631 deepdd = x * df_d + T * dsepdd;
632 deepda = -ye * ytot1 * (free + df_d * din) + T * dsepda;
633 deepdz = ytot1 * (free + ye * df_d * rho) + T * dsepdz;
634
635 // coulomb
636
637 constexpr double a1 = -0.898004;
638 constexpr double b1 = 0.96786;
639 constexpr double c1 = 0.220703;
640 constexpr double d1 = -0.86097;
641 constexpr double e1 = 2.5269;
642 constexpr double a2 = 0.29561;
643 constexpr double b2 = 1.9885;
644 constexpr double c2 = 0.288675;
645
646 z = 4 * pi * third;
647 s = z * xni;
648 double dsdd = z * dxnidd;
649 double dsda = z * dxnida;
650
651 double lami = pow(s,-1./3.);
652 double inv_lami = 1.0/lami;
653 z = -lami/3.;
654 double lamidd = z * dsdd/s;
655 double lamida = z * dsda/s;
656
657 double plasg = zbar * zbar * esqu * kTinv * inv_lami;
658 z = -plasg * inv_lami;
659 double plasgdd = z * lamidd;
660 double plasgda = z * lamida;
661 double plasgdt = -plasg*kTinv * kerg;
662 double plasgdz = 2 * plasg/zbar;
663
664 // yakovlev & shalybkov 1989 equations 82, 85, 86, 87
665 if (plasg >= 1.0) {
666 x = pow(plasg,0.25);
667 y = avo * ytot1 * kerg;
668 ecoul = y * T * (a1*plasg + b1*x + c1/x + d1);
669 pcoul = third * rho * ecoul;
670 scoul = -y * (3*b1*x - 5*c1/x + d1 * (log(plasg) - 1) - e1);
671
672 y = avo * ytot1 * kT * (a1 + 0.25/plasg*(b1*x - c1/x));
673 decouldd = y * plasgdd;
674 decouldt = y * plasgdt + ecoul/T;
675 decoulda = y * plasgda - ecoul/abar;
676 decouldz = y * plasgdz;
677
678 y = third * rho;
679 dpcouldd = third * ecoul + y * decouldd;
680 dpcouldt = y * decouldt;
681 dpcoulda = y * decoulda;
682 dpcouldz = y * decouldz;
683
684
685 y = -avo * kerg / (abar*plasg) * (0.75*b1*x+1.25*c1/x+d1);
686 dscouldd = y * plasgdd;
687 dscouldt = y * plasgdt;
688 dscoulda = y * plasgda - scoul/abar;
689 dscouldz = y * plasgdz;
690 }
691
692 // yakovlev & shalybkov 1989 equations 102, 103, 104
693 else {
694 x = plasg*sqrt(plasg);
695 y = pow(plasg,b2);
696 z = c2 * x - third * a2 * y;
697 pcoul = -pion * z;
698 ecoul = 3 * pcoul/rho;
699 scoul = -(avo/abar) * kerg * (c2*x -a2*(b2-1)/b2*y);
700
701 s = 1.5*c2*x/plasg - third*a2*b2*y/plasg;
702 dpcouldd = -dpiondd*z - pion*s*plasgdd;
703 dpcouldt = -dpiondt*z - pion*s*plasgdt;
704 dpcoulda = -dpionda*z - pion*s*plasgda;
705 dpcouldz = -dpiondz*z - pion*s*plasgdz;
706
707 s = 3/rho;
708 decouldd = s * dpcouldd - ecoul/rho;
709 decouldt = s * dpcouldt;
710 decoulda = s * dpcoulda;
711 decouldz = s * dpcouldz;
712
713 s = -avo*kerg/(abar*plasg)*(1.5*c2*x - a2*(b2-1)*y);
714 dscouldd = s * plasgdd;
715 dscouldt = s * plasgdt;
716 dscoulda = s * plasgda - scoul/abar;
717 dscouldz = s * plasgdz;
718 }
719
720 // "bomb proof"
721 x = prad + pion + pele + pcoul;
722 y = erad + eion + eele + ecoul;
723 z = srad + sion + sele + scoul;
724
725 if(x <= 0 || y <= 0) {
726 // zero out coulomb terms
727 pcoul = 0; dpcouldd=0; dpcouldt=0; dpcoulda=0; dpcouldz=0;
728 ecoul = 0; decouldd=0; decouldt=0; decoulda=0; decouldz=0;
729 scoul = 0; dscouldd=0; dscouldt=0; dscoulda=0; dscouldz=0;
730 }
731
732
733 // gas subtotals
734 pgas = pion + pele + pcoul;
735 egas = eion + eele + ecoul;
736 sgas = sion + sele + scoul;
737
738 LOG_DEBUG(logger, "pgas = {}", pgas);
739 LOG_DEBUG(logger, "prad = {}", prad);
740 LOG_DEBUG(logger, "pion = {}", pion);
741 LOG_DEBUG(logger, "pele = {}", pele);
742 LOG_DEBUG(logger, "pcoul = {}", pcoul);
743
744 dpgasdd = dpiondd + dpepdd + dpcouldd;
745 dpgasdt = dpiondt + dpepdt + dpcouldt;
746 dpgasda = dpionda + dpepda + dpcoulda;
747 dpgasdz = dpiondz + dpepdz + dpcouldz;
748
749 LOG_DEBUG(logger, "dpgasdd = {}", dpgasdd);
750 LOG_DEBUG(logger, "dpraddd = {}", dpraddd);
751 LOG_DEBUG(logger, "dpiondd = {}", dpiondd);
752 LOG_DEBUG(logger, "dpeledd = {}", dpepdd);
753 LOG_DEBUG(logger, "dpcouldd = {}", dpcouldd);
754
755 degasdd = deiondd + deepdd + decouldd;
756 degasdt = deiondt + deepdt + decouldt;
757 degasda = deionda + deepda + decoulda;
758 degasdz = deiondz + deepdz + decouldz;
759
760 LOG_DEBUG(logger, "degasdd = {}", degasdd);
761 LOG_DEBUG(logger, "deraddd = {}", deraddd);
762 LOG_DEBUG(logger, "deiondd = {}", deiondd);
763 LOG_DEBUG(logger, "deeledd = {}", deepdd);
764 LOG_DEBUG(logger, "decouldd = {}", decouldd);
765
766 LOG_DEBUG(logger, "degasdt = {}", degasdt);
767 LOG_DEBUG(logger, "deraddt = {}", deraddt);
768 LOG_DEBUG(logger, "deiondt = {}", deiondt);
769 LOG_DEBUG(logger, "deeledt = {}", deepdt);
770 LOG_DEBUG(logger, "decouldt = {}", decouldt);
771
772 dsgasdd = dsiondd + dsepdd + dscouldd;
773 dsgasdt = dsiondt + dsepdt + dscouldt;
774 dsgasda = dsionda + dsepda + dscoulda;
775 dsgasdz = dsiondz + dsepdz + dscouldz;
776
777
778 LOG_DEBUG(logger, "sgas = {}", sgas);
779 LOG_DEBUG(logger, "srad = {}", srad);
780 LOG_DEBUG(logger, "sion = {}", sion);
781 LOG_DEBUG(logger, "sele = {}", sele);
782 LOG_DEBUG(logger, "scoul = {}", scoul);
783
784 LOG_DEBUG(logger, "dsgasdd = {}", dsgasdd);
785 LOG_DEBUG(logger, "dsraddd = {}", dsraddd);
786 LOG_DEBUG(logger, "dsiondd = {}", dsiondd);
787 LOG_DEBUG(logger, "dseledd = {}", dsepdd);
788 LOG_DEBUG(logger, "dscouldd = {}", dscouldd);
789
790 LOG_DEBUG(logger, "dsgasdt = {}", dsgasdt);
791 LOG_DEBUG(logger, "dsraddt = {}", dsraddt);
792 LOG_DEBUG(logger, "dsiondt = {}", dsiondt);
793 LOG_DEBUG(logger, "dseledt = {}", dsepdt);
794 LOG_DEBUG(logger, "dscouldt = {}", dscouldt);
795
796 //total
797 ptot = prad + pgas;
798 stot = srad + sgas;
799 etot = erad + egas;
800
801 dpresdd = dpraddd + dpgasdd;
802 dpresdt = dpraddt + dpgasdt;
803 dpresda = dpradda + dpgasda;
804 dpresdz = dpraddz + dpgasdz;
805
806 denerdd = deraddd + degasdd;
807 denerdt = deraddt + degasdt;
808 denerda = deradda + degasda;
809 denerdz = deraddz + degasdz;
810
811 dentrdd = dsraddd + dsgasdd;
812 dentrdt = dsraddt + dsgasdt;
813 dentrda = dsradda + dsgasda;
814 dentrdz = dsraddz + dsgasdz;
815
816 // derived quantities
817 double zz = ptot * rhoi;
818 double zzi = rho/ptot;
819 double chit = (T/ptot)*dpresdt;
820 double chirho = zzi*dpresdd;
821 double cv = denerdt;
822 x = zz * chit/(T*cv);
823 double gamma3 = x + 1.0;
824 double gamma1 = chit * x + chirho;
825 double grad_ad = x/gamma1;
826 double gamma2 = 1.0/(1.0 - grad_ad);
827 double cp = cv * gamma1/chirho;
828 z = 1.0 + (etot + clight2)*zzi;
829 double csound = clight * sqrt(gamma1/z);
830
831 // package in q:
832 serif::eos::helmholtz::HELMEOSOutput eos;
833 eos.ptot = ptot; eos.etot = etot; eos.stot = stot;
834 eos.pgas = pgas; eos.egas = egas; eos.sgas = sgas;
835 eos.prad = prad; eos.erad = erad; eos.srad = srad;
836 eos.chiT = chit, eos.chiRho = chirho, eos.csound = csound;
837 eos.gamma1 = gamma1, eos.gamma2 = gamma2; eos.gamma3 = gamma3;
838 eos.cV = cv, eos.cP = cp, eos.grad_ad = grad_ad;
839 eos.etaele = etaele, eos.xnefer = xnefer, eos.ye = ye;
840
841 eos.dpresdd = dpresdd; eos.dentrdd = dentrdd; eos.denerdd = denerdd;
842 eos.dpresdt = dpresdt; eos.dentrdt = dentrdt; eos.denerdt = denerdt;
843 eos.dpresda = dpresda; eos.dentrda = dentrda; eos.denerda = denerda;
844 eos.dpresdz = dpresdz; eos.dentrdz = dentrdz; eos.denerdz = denerdz;
845
846 // maxwell relations
847 x = rho * rho;
848 eos.dse = T * dentrdt/denerdt - 1.0;
849 eos.dpe = (denerdd*x + T*dpresdt)/ptot - 1.0;
850 eos.dsp = -dentrdd*x/dpresdt - 1.0;
851
852 return eos;
853 // TODO: In future this might be made more preformant by avoiding the copy.
854 }
855}
serif::constant::Constants
Class to manage a collection of constants.
Definition const.h:63
serif::constant::Constants::getInstance
static Constants & getInstance()
get instance of constants singleton
Definition const.h:99
serif::constant::Constants::get
Constant get(const std::string &key) const
Get a constant by key.
Definition const.cpp:42
serif::probe::LogManager
Class to manage logging operations.
Definition probe.h:79
serif::probe::LogManager::getInstance
static LogManager & getInstance()
Get the singleton instance of LogManager.
Definition probe.h:103
serif::eos::helmholtz::ddpsi1
double ddpsi1(const double z)
Second derivative of the interpolating polynomial function psi1.
Definition helm.cpp:80
serif::eos::helmholtz::dlo
const double dlo
Definition helm.h:54
serif::eos::helmholtz::psi0
double psi0(const double z)
Interpolating polynomial function psi0.
Definition helm.cpp:70
serif::eos::helmholtz::dpsi0
double dpsi0(const double z)
Derivative of the interpolating polynomial function psi0.
Definition helm.cpp:72
serif::eos::helmholtz::psi1
double psi1(const double z)
Interpolating polynomial function psi1.
Definition helm.cpp:76
serif::eos::helmholtz::xdpsi1
double xdpsi1(const double z)
Derivative of the interpolating polynomial function xpsi1.
Definition helm.cpp:94
serif::eos::helmholtz::tstpi
const double tstpi
Definition helm.h:55
serif::eos::helmholtz::h3
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)
Interpolating polynomial function h3.
Definition helm.cpp:96
serif::eos::helmholtz::dstp
const double dstp
Definition helm.h:56
serif::eos::helmholtz::ddpsi2
double ddpsi2(const double z)
Second derivative of the interpolating polynomial function psi2.
Definition helm.cpp:86
serif::eos::helmholtz::h5
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)
Interpolating polynomial function h5.
Definition helm.cpp:108
serif::eos::helmholtz::xdpsi0
double xdpsi0(const double z)
Derivative of the interpolating polynomial function xpsi0.
Definition helm.cpp:90
serif::eos::helmholtz::dpsi2
double dpsi2(const double z)
Derivative of the interpolating polynomial function psi2.
Definition helm.cpp:84
serif::eos::helmholtz::dpsi1
double dpsi1(const double z)
Derivative of the interpolating polynomial function psi1.
Definition helm.cpp:78
serif::eos::helmholtz::xpsi1
double xpsi1(const double z)
Interpolating polynomial function xpsi1.
Definition helm.cpp:92
serif::eos::helmholtz::ddpsi0
double ddpsi0(const double z)
Second derivative of the interpolating polynomial function psi0.
Definition helm.cpp:74
serif::eos::helmholtz::tlo
const double tlo
Definition helm.h:53
serif::eos::helmholtz::dstpi
const double dstpi
Definition helm.h:56
serif::eos::helmholtz::heap_deallocate_contiguous_2D_memory
void heap_deallocate_contiguous_2D_memory(double **array)
Definition helm.cpp:66
serif::eos::helmholtz::xpsi0
double xpsi0(const double z)
Interpolating polynomial function xpsi0.
Definition helm.cpp:88
serif::eos::helmholtz::heap_allocate_contiguous_2D_memory
double ** heap_allocate_contiguous_2D_memory(const int rows, const int cols)
Definition helm.cpp:52
serif::eos::helmholtz::get_helm_EOS
serif::eos::helmholtz::HELMEOSOutput get_helm_EOS(serif::eos::helmholtz::HELMEOSInput &q, const serif::eos::helmholtz::HELMTable &table)
Calculate the Helmholtz EOS components.
Definition helm.cpp:243
serif::eos::helmholtz::tstp
const double tstp
Definition helm.h:55
serif::eos::helmholtz::psi2
double psi2(const double z)
Interpolating polynomial function psi2.
Definition helm.cpp:82
serif::eos::helmholtz::read_helm_table
std::unique_ptr< HELMTable > read_helm_table(const std::string &filename)
Read the Helmholtz EOS table from a file.
Definition helm.cpp:132
serif::constant::Constant::value
const double value
Value of the constant.
Definition const.h:33
serif::eos::helmholtz::HELMEOSInput
Structure to hold the input parameters for the EOS calculation.
Definition helm.h:176
serif::eos::helmholtz::HELMEOSInput::abar
double abar
Mean atomic mass.
Definition helm.h:179
serif::eos::helmholtz::HELMEOSInput::zbar
double zbar
Mean atomic number.
Definition helm.h:180
serif::eos::helmholtz::HELMTable
Structure to hold the Helmholtz EOS table data.
Definition helm.h:63