feat(python): added robust python bindings covering the entire codebase
This commit is contained in:
529
src/lib/engine/engine_approx8.cpp
Normal file
529
src/lib/engine/engine_approx8.cpp
Normal file
@@ -0,0 +1,529 @@
|
||||
/* ***********************************************************************
|
||||
//
|
||||
// Copyright (C) 2025 -- The 4D-STAR Collaboration
|
||||
// File Author: Emily Boudreaux
|
||||
// Last Modified: March 21, 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
|
||||
//
|
||||
// *********************************************************************** */
|
||||
#include <cmath>
|
||||
#include <stdexcept>
|
||||
#include <array>
|
||||
|
||||
#include <boost/numeric/odeint.hpp>
|
||||
|
||||
#include "fourdst/constants/const.h"
|
||||
#include "fourdst/config/config.h"
|
||||
#include "quill/LogMacros.h"
|
||||
|
||||
#include "gridfire/engine/engine_approx8.h"
|
||||
#include "gridfire/network.h"
|
||||
|
||||
/* Nuclear reaction network in cgs units based on Frank Timmes' "approx8".
|
||||
At this time it does neither screening nor neutrino losses. It includes
|
||||
the following 8 isotopes:
|
||||
|
||||
h1
|
||||
he3
|
||||
he4
|
||||
c12
|
||||
n14
|
||||
o16
|
||||
ne20
|
||||
mg24
|
||||
|
||||
and the following nuclear reactions:
|
||||
|
||||
---pp chain---
|
||||
p(p,e+)d
|
||||
d(p,g)he3
|
||||
he3(he3,2p)he4
|
||||
|
||||
---CNO cycle---
|
||||
c12(p,g)n13 - n13 -> c13 + p -> n14
|
||||
n14(p,g)o15 - o15 + p -> c12 + he4
|
||||
n14(a,g)f18 - proceeds to ne20
|
||||
n15(p,a)c12 - / these two n15 reactions are \ CNO I
|
||||
n15(p,g)o16 - \ used to calculate a fraction / CNO II
|
||||
o16(p,g)f17 - f17 + e -> o17 and then o17 + p -> n14 + he4
|
||||
|
||||
---alpha captures---
|
||||
c12(a,g)o16
|
||||
triple alpha
|
||||
o16(a,g)ne20
|
||||
ne20(a,g)mg24
|
||||
c12(c12,a)ne20
|
||||
c12(o16,a)mg24
|
||||
|
||||
At present the rates are all evaluated using a fitting function.
|
||||
The coefficients to the fit are from reaclib.jinaweb.org .
|
||||
|
||||
*/
|
||||
|
||||
namespace gridfire::approx8{
|
||||
|
||||
// using namespace std;
|
||||
using namespace boost::numeric::odeint;
|
||||
|
||||
//helper functions
|
||||
// a function to multiply two arrays and then sum the resulting elements: sum(a*b)
|
||||
double sum_product( const vec7 &a, const vec7 &b){
|
||||
double sum=0;
|
||||
for (size_t i=0; i < a.size(); i++) {
|
||||
sum += a[i] * b[i];
|
||||
}
|
||||
return sum;
|
||||
}
|
||||
|
||||
// the fit to the nuclear reaction rates is of the form:
|
||||
// exp( a0 + a1/T9 + a2/T9^(1/3) + a3*T9^(1/3) + a4*T9 + a5*T9^(5/3) + log(T9) )
|
||||
// this function returns an array of the T9 terms in that order, where T9 is the temperatures in GigaKelvin
|
||||
vec7 get_T9_array(const double &T) {
|
||||
vec7 arr;
|
||||
const double T9=1e-9*T;
|
||||
const double T913=pow(T9,1./3.);
|
||||
|
||||
arr[0]=1;
|
||||
arr[1]=1/T9;
|
||||
arr[2]=1/T913;
|
||||
arr[3]=T913;
|
||||
arr[4]=T9;
|
||||
arr[5]=pow(T9,5./3.);
|
||||
arr[6]=log(T9);
|
||||
|
||||
return arr;
|
||||
}
|
||||
|
||||
// this function uses the two preceding functions to evaluate the rate given the T9 array and the coefficients
|
||||
double rate_fit(const vec7 &T9, const vec7 &coef){
|
||||
return exp(sum_product(T9,coef));
|
||||
}
|
||||
|
||||
// p + p -> d; this, like some of the other rates, this is a composite of multiple fits
|
||||
double pp_rate(const vec7 &T9) {
|
||||
constexpr vec7 a1 = {-34.78630, 0,-3.511930, 3.100860, -0.1983140, 1.262510e-2, -1.025170};
|
||||
constexpr vec7 a2 = { -4.364990e+1,-2.460640e-3,-2.750700,-4.248770e-1,1.598700e-2,-6.908750e-4,-2.076250e-1};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2);
|
||||
}
|
||||
|
||||
// p + d -> he3
|
||||
double dp_rate(const vec7 &T9) {
|
||||
constexpr vec7 a1 = {7.528980, 0, -3.720800, 0.8717820, 0, 0,-0.6666670};
|
||||
constexpr vec7 a2 = {8.935250, 0, -3.720800, 0.1986540, 0, 0, 0.3333330};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2);
|
||||
}
|
||||
|
||||
// he3 + he3 -> he4 + 2p
|
||||
double he3he3_rate(const vec7 &T9){
|
||||
constexpr vec7 a = {2.477880e+01,0,-12.27700,-0.1036990,-6.499670e-02,1.681910e-02,-6.666670e-01};
|
||||
return rate_fit(T9,a);
|
||||
}
|
||||
|
||||
// he3(he3,2p)he4 ** (missing both coefficients but have a reaction)
|
||||
double he3he4_rate(const vec7 &T9){
|
||||
constexpr vec7 a1 = {1.560990e+01,0.000000e+00,-1.282710e+01,-3.082250e-02,-6.546850e-01,8.963310e-02,-6.666670e-01};
|
||||
constexpr vec7 a2 = {1.770750e+01,0.000000e+00,-1.282710e+01,-3.812600e+00,9.422850e-02,-3.010180e-03,1.333330e+00};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2);
|
||||
}
|
||||
|
||||
// he4 + he4 + he4 -> c12 ** (missing middle coefficient but have other two)
|
||||
double triple_alpha_rate(const vec7 &T9){
|
||||
constexpr vec7 a1 = {-9.710520e-01,0.000000e+00,-3.706000e+01,2.934930e+01,-1.155070e+02,-1.000000e+01,-1.333330e+00};
|
||||
constexpr vec7 a2 = {-1.178840e+01,-1.024460e+00,-2.357000e+01,2.048860e+01,-1.298820e+01,-2.000000e+01,-2.166670e+00};
|
||||
constexpr vec7 a3 = {-2.435050e+01,-4.126560e+00,-1.349000e+01,2.142590e+01,-1.347690e+00,8.798160e-02,-1.316530e+01};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
|
||||
}
|
||||
|
||||
// c12 + p -> n13
|
||||
double c12p_rate(const vec7 &T9){
|
||||
constexpr vec7 a1={1.714820e+01,0.000000e+00,-1.369200e+01,-2.308810e-01,4.443620e+00,-3.158980e+00,-6.666670e-01};
|
||||
constexpr vec7 a2={1.754280e+01,-3.778490e+00,-5.107350e+00,-2.241110e+00,1.488830e-01,0.000000e+00,-1.500000e+00};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2);
|
||||
}
|
||||
|
||||
// c12 + he4 -> o16 ** (missing first coefficient but have the second)
|
||||
double c12a_rate(const vec7 &T9){
|
||||
constexpr vec7 a1={6.965260e+01,-1.392540e+00,5.891280e+01,-1.482730e+02,9.083240e+00,-5.410410e-01,7.035540e+01};
|
||||
constexpr vec7 a2={2.546340e+02,-1.840970e+00,1.034110e+02,-4.205670e+02,6.408740e+01,-1.246240e+01,1.373030e+02};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2);
|
||||
}
|
||||
|
||||
// n14(p,g)o15 - o15 + p -> c12 + he4
|
||||
double n14p_rate(const vec7 &T9){
|
||||
constexpr vec7 a1={1.701000e+01,0.000000e+00,-1.519300e+01,-1.619540e-01,-7.521230e+00,-9.875650e-01,-6.666670e-01};
|
||||
constexpr vec7 a2={2.011690e+01,0.000000e+00,-1.519300e+01,-4.639750e+00,9.734580e+00,-9.550510e+00,3.333330e-01};
|
||||
constexpr vec7 a3={7.654440e+00,-2.998000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
|
||||
constexpr vec7 a4={6.735780e+00,-4.891000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,6.820000e-02};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3) + rate_fit(T9,a4);
|
||||
}
|
||||
|
||||
// n14(a,g)f18 assumed to go on to ne20
|
||||
double n14a_rate(const vec7 &T9){
|
||||
constexpr vec7 a1={2.153390e+01,0.000000e+00,-3.625040e+01,0.000000e+00,0.000000e+00,-5.000000e+00,-6.666670e-01};
|
||||
constexpr vec7 a2={1.968380e-01,-5.160340e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
|
||||
constexpr vec7 a3={1.389950e+01,-1.096560e+01,-5.622700e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
|
||||
}
|
||||
|
||||
// n15(p,a)c12 (CNO I)
|
||||
double n15pa_rate(const vec7 &T9){
|
||||
constexpr vec7 a1 = {2.747640e+01,0.000000e+00,-1.525300e+01,1.593180e+00,2.447900e+00,-2.197080e+00,-6.666670e-01};
|
||||
constexpr vec7 a2 = {-4.873470e+00,-2.021170e+00,0.000000e+00,3.084970e+01,-8.504330e+00,-1.544260e+00,-1.500000e+00};
|
||||
constexpr vec7 a3 = {2.089720e+01,-7.406000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
|
||||
constexpr vec7 a4 = {-6.575220e+00,-1.163800e+00,0.000000e+00,2.271050e+01,-2.907070e+00,2.057540e-01,-1.500000e+00};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3) + rate_fit(T9,a4);
|
||||
}
|
||||
|
||||
// n15(p,g)o16 (CNO II)
|
||||
double n15pg_rate(const vec7 &T9){
|
||||
constexpr vec7 a1 = {2.001760e+01,0.000000e+00,-1.524000e+01,3.349260e-01,4.590880e+00,-4.784680e+00,-6.666670e-01};
|
||||
constexpr vec7 a2 = {6.590560e+00,-2.923150e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
|
||||
constexpr vec7 a3 = {1.454440e+01,-1.022950e+01,0.000000e+00,0.000000e+00,4.590370e-02,0.000000e+00,-1.500000e+00};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
|
||||
}
|
||||
|
||||
double n15pg_frac(const vec7 &T9){
|
||||
const double f1=n15pg_rate(T9);
|
||||
const double f2=n15pa_rate(T9);
|
||||
return f1/(f1+f2);
|
||||
}
|
||||
|
||||
// o16(p,g)f17 then f17 -> o17(p,a)n14
|
||||
double o16p_rate(const vec7 &T9){
|
||||
constexpr vec7 a={1.909040e+01,0.000000e+00,-1.669600e+01,-1.162520e+00,2.677030e-01,-3.384110e-02,-6.666670e-01};
|
||||
return rate_fit(T9,a);
|
||||
}
|
||||
|
||||
// o16(a,g)ne20
|
||||
double o16a_rate(const vec7 &T9){
|
||||
constexpr vec7 a1={2.390300e+01,0.000000e+00,-3.972620e+01,-2.107990e-01,4.428790e-01,-7.977530e-02,-6.666670e-01};
|
||||
constexpr vec7 a2={3.885710e+00,-1.035850e+01,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
|
||||
constexpr vec7 a3={9.508480e+00,-1.276430e+01,0.000000e+00,-3.659250e+00,7.142240e-01,-1.075080e-03,-1.500000e+00};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
|
||||
}
|
||||
|
||||
// ne20(a,g)mg24
|
||||
double ne20a_rate(const vec7 &T9){
|
||||
constexpr vec7 a1={2.450580e+01,0.000000e+00,-4.625250e+01,5.589010e+00,7.618430e+00,-3.683000e+00,-6.666670e-01};
|
||||
constexpr vec7 a2={-3.870550e+01,-2.506050e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
|
||||
constexpr vec7 a3={1.983070e+00,-9.220260e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
|
||||
constexpr vec7 a4={-8.798270e+00,-1.278090e+01,0.000000e+00,1.692290e+01,-2.573250e+00,2.089970e-01,-1.500000e+00};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3) + rate_fit(T9,a4);
|
||||
}
|
||||
|
||||
// c12(c12,a)ne20
|
||||
double c12c12_rate(const vec7 &T9){
|
||||
constexpr vec7 a={6.128630e+01,0.000000e+00,-8.416500e+01,-1.566270e+00,-7.360840e-02,-7.279700e-02,-6.666670e-01};
|
||||
return rate_fit(T9,a);
|
||||
}
|
||||
|
||||
// c12(o16,a)mg24
|
||||
double c12o16_rate(const vec7 &T9){
|
||||
constexpr vec7 a={4.853410e+01,3.720400e-01,-1.334130e+02,5.015720e+01,-3.159870e+00,1.782510e-02,-2.370270e+01};
|
||||
return rate_fit(T9,a);
|
||||
}
|
||||
|
||||
|
||||
// for Boost ODE solvers either a struct or a class is required
|
||||
|
||||
// a Jacobian matrix for implicit solvers
|
||||
|
||||
void Jacobian::operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) const {
|
||||
fourdst::constant::Constants& constants = fourdst::constant::Constants::getInstance();
|
||||
const double avo = constants.get("N_a").value;
|
||||
const double clight = constants.get("c").value;
|
||||
// EOS
|
||||
const vec7 T9=get_T9_array(y[Approx8Net::iTemp]);
|
||||
|
||||
// evaluate rates once per call
|
||||
const double rpp=pp_rate(T9);
|
||||
const double r33=he3he3_rate(T9);
|
||||
const double r34=he3he4_rate(T9);
|
||||
const double r3a=triple_alpha_rate(T9);
|
||||
const double rc12p=c12p_rate(T9);
|
||||
const double rc12a=c12a_rate(T9);
|
||||
const double rn14p=n14p_rate(T9);
|
||||
const double rn14a=n14a_rate(T9);
|
||||
const double ro16p=o16p_rate(T9);
|
||||
const double ro16a=o16a_rate(T9);
|
||||
const double rne20a=ne20a_rate(T9);
|
||||
const double r1212=c12c12_rate(T9);
|
||||
const double r1216=c12o16_rate(T9);
|
||||
|
||||
const double pFrac=n15pg_frac(T9);
|
||||
const double aFrac=1-pFrac;
|
||||
|
||||
const double yh1 = y[Approx8Net::ih1];
|
||||
const double yhe3 = y[Approx8Net::ihe3];
|
||||
const double yhe4 = y[Approx8Net::ihe4];
|
||||
const double yc12 = y[Approx8Net::ic12];
|
||||
const double yn14 = y[Approx8Net::in14];
|
||||
const double yo16 = y[Approx8Net::io16];
|
||||
const double yne20 = y[Approx8Net::ine20];
|
||||
|
||||
// zero all elements to begin
|
||||
for (int i=0; i < Approx8Net::nVar; i++) {
|
||||
for (int j=0; j < Approx8Net::nVar; j++) {
|
||||
J(i,j)=0.0;
|
||||
}
|
||||
}
|
||||
|
||||
// h1 jacobian elements
|
||||
J(Approx8Net::ih1,Approx8Net::ih1) = -3*yh1*rpp - 2*yc12*rc12p -2*yn14*rn14p -2*yo16*ro16p;
|
||||
J(Approx8Net::ih1,Approx8Net::ihe3) = 2*yhe3*r33 - yhe4*r34;
|
||||
J(Approx8Net::ih1,Approx8Net::ihe4) = -yhe3*r34;
|
||||
J(Approx8Net::ih1,Approx8Net::ic12) = -2*yh1*rc12p;
|
||||
J(Approx8Net::ih1,Approx8Net::in14) = -2*yh1*rn14p;
|
||||
J(Approx8Net::ih1,Approx8Net::io16) = -2*yh1*ro16p;
|
||||
|
||||
// he3 jacobian elements
|
||||
J(Approx8Net::ihe3,Approx8Net::ih1) = yh1*rpp;
|
||||
J(Approx8Net::ihe3,Approx8Net::ihe3) = -2*yhe3*r33 - yhe4*r34;
|
||||
J(Approx8Net::ihe3,Approx8Net::ihe4) = -yhe3*r34;
|
||||
|
||||
// he4 jacobian elements
|
||||
J(Approx8Net::ihe4,Approx8Net::ih1) = yn14*aFrac*rn14p + yo16*ro16p;
|
||||
J(Approx8Net::ihe4,Approx8Net::ihe3) = yhe3*r33 - yhe4*r34;
|
||||
J(Approx8Net::ihe4,Approx8Net::ihe4) = yhe3*r34 - 1.5*yhe4*yhe4*r3a - yc12*rc12a - 1.5*yn14*rn14a - yo16*ro16a - yne20*rne20a;
|
||||
J(Approx8Net::ihe4,Approx8Net::ic12) = -yhe4*rc12a + yc12*r1212 + yo16*r1216;
|
||||
J(Approx8Net::ihe4,Approx8Net::in14) = yh1*aFrac*rn14p - 1.5*yhe4*rn14a;
|
||||
J(Approx8Net::ihe4,Approx8Net::io16) = yh1*ro16p - yhe4*ro16a + yc12*r1216;
|
||||
J(Approx8Net::ihe4,Approx8Net::ine20) = -yhe4*rne20a;
|
||||
|
||||
// c12 jacobian elements
|
||||
J(Approx8Net::ic12,Approx8Net::ih1) = -yc12*rc12p + yn14*aFrac*rn14p;
|
||||
J(Approx8Net::ic12,Approx8Net::ihe4) = 0.5*yhe4*yhe4*r3a - yhe4*rc12a;
|
||||
J(Approx8Net::ic12,Approx8Net::ic12) = -yh1*rc12p - yhe4*rc12a - yo16*r1216 - 2*yc12*r1212;
|
||||
J(Approx8Net::ic12,Approx8Net::in14) = yh1*yn14*aFrac*rn14p;
|
||||
J(Approx8Net::ic12,Approx8Net::io16) = -yc12*r1216;
|
||||
|
||||
// n14 jacobian elements
|
||||
J(Approx8Net::in14,Approx8Net::ih1) = yc12*rc12p - yn14*rn14p + yo16*ro16p;
|
||||
J(Approx8Net::in14,Approx8Net::ihe4) = -yn14*rn14a;
|
||||
J(Approx8Net::in14,Approx8Net::ic12) = yh1*rc12p;
|
||||
J(Approx8Net::in14,Approx8Net::in14) = -yh1*rn14p - yhe4*rn14a;
|
||||
J(Approx8Net::in14,Approx8Net::io16) = yo16*ro16p;
|
||||
|
||||
// o16 jacobian elements
|
||||
J(Approx8Net::io16,Approx8Net::ih1) = yn14*pFrac*rn14p - yo16*ro16p;
|
||||
J(Approx8Net::io16,Approx8Net::ihe4) = yc12*rc12a - yo16*ro16a;
|
||||
J(Approx8Net::io16,Approx8Net::ic12) = yhe4*rc12a - yo16*r1216;
|
||||
J(Approx8Net::io16,Approx8Net::in14) = yh1*pFrac*rn14p;
|
||||
J(Approx8Net::io16,Approx8Net::io16) = yh1*ro16p - yc12*r1216 -yhe4*ro16a;
|
||||
|
||||
// ne20 jacobian elements
|
||||
J(Approx8Net::ine20,Approx8Net::ihe4) = yn14*rn14a + yo16*ro16a - yne20*rne20a;
|
||||
J(Approx8Net::ine20,Approx8Net::ic12) = yc12*r1212;
|
||||
J(Approx8Net::ine20,Approx8Net::in14) = yhe4*rn14a;
|
||||
J(Approx8Net::ine20,Approx8Net::io16) = yo16*ro16a;
|
||||
J(Approx8Net::ine20,Approx8Net::ine20) = -yhe4*rne20a;
|
||||
|
||||
// mg24 jacobian elements
|
||||
J(Approx8Net::img24,Approx8Net::ihe4) = yne20*rne20a;
|
||||
J(Approx8Net::img24,Approx8Net::ic12) = yo16*r1216;
|
||||
J(Approx8Net::img24,Approx8Net::io16) = yc12*r1216;
|
||||
J(Approx8Net::img24,Approx8Net::ine20) = yhe4*rne20a;
|
||||
|
||||
// energy accounting
|
||||
for (int j=0; j<Approx8Net::nIso; j++) {
|
||||
for (int i=0; i<Approx8Net::nIso; i++) {
|
||||
J(Approx8Net::iEnergy,j) += J(i,j)*Approx8Net::mIon[i];
|
||||
}
|
||||
J(Approx8Net::iEnergy,j) *= -avo*clight*clight;
|
||||
}
|
||||
}
|
||||
|
||||
void ODE::operator() ( const vector_type &y, vector_type &dydt, double /* t */) const {
|
||||
const fourdst::constant::Constants& constants = fourdst::constant::Constants::getInstance();
|
||||
const double avo = constants.get("N_a").value;
|
||||
const double clight = constants.get("c").value;
|
||||
|
||||
// EOS
|
||||
const double T = y[Approx8Net::iTemp];
|
||||
const double den = y[Approx8Net::iDensity];
|
||||
const vec7 T9=get_T9_array(T);
|
||||
|
||||
// rates
|
||||
const double rpp=den*pp_rate(T9);
|
||||
const double r33=den*he3he3_rate(T9);
|
||||
const double r34=den*he3he4_rate(T9);
|
||||
const double r3a=den*den*triple_alpha_rate(T9);
|
||||
const double rc12p=den*c12p_rate(T9);
|
||||
const double rc12a=den*c12a_rate(T9);
|
||||
const double rn14p=den*n14p_rate(T9);
|
||||
const double rn14a=n14a_rate(T9);
|
||||
const double ro16p=den*o16p_rate(T9);
|
||||
const double ro16a=den*o16a_rate(T9);
|
||||
const double rne20a=den*ne20a_rate(T9);
|
||||
const double r1212=den*c12c12_rate(T9);
|
||||
const double r1216=den*c12o16_rate(T9);
|
||||
|
||||
const double pFrac=n15pg_frac(T9);
|
||||
const double aFrac=1-pFrac;
|
||||
|
||||
const double yh1 = y[Approx8Net:: ih1];
|
||||
const double yhe3 = y[Approx8Net:: ihe3];
|
||||
const double yhe4 = y[Approx8Net:: ihe4];
|
||||
const double yc12 = y[Approx8Net:: ic12];
|
||||
const double yn14 = y[Approx8Net:: in14];
|
||||
const double yo16 = y[Approx8Net:: io16];
|
||||
const double yne20 = y[Approx8Net::ine20];
|
||||
|
||||
dydt[Approx8Net::ih1] = -1.5*yh1*yh1*rpp;
|
||||
dydt[Approx8Net::ih1] += yhe3*yhe3*r33;
|
||||
dydt[Approx8Net::ih1] += -yhe3*yhe4*r34;
|
||||
dydt[Approx8Net::ih1] += -2*yh1*yc12*rc12p;
|
||||
dydt[Approx8Net::ih1] += -2*yh1*yn14*rn14p;
|
||||
dydt[Approx8Net::ih1] += -2*yh1*yo16*ro16p;
|
||||
|
||||
dydt[Approx8Net::ihe3] = 0.5*yh1*yh1*rpp;
|
||||
dydt[Approx8Net::ihe3] += -yhe3*yhe3*r33;
|
||||
dydt[Approx8Net::ihe3] += -yhe3*yhe4*r34;
|
||||
|
||||
dydt[Approx8Net::ihe4] = 0.5*yhe3*yhe3*r33;
|
||||
dydt[Approx8Net::ihe4] += yhe3*yhe4*r34;
|
||||
dydt[Approx8Net::ihe4] += -yhe4*yc12*rc12a;
|
||||
dydt[Approx8Net::ihe4] += yh1*yn14*aFrac*rn14p;
|
||||
dydt[Approx8Net::ihe4] += yh1*yo16*ro16p;
|
||||
dydt[Approx8Net::ihe4] += -0.5*yhe4*yhe4*yhe4*r3a;
|
||||
dydt[Approx8Net::ihe4] += -yhe4*yo16*ro16a;
|
||||
dydt[Approx8Net::ihe4] += 0.5*yc12*yc12*r1212;
|
||||
dydt[Approx8Net::ihe4] += yc12*yo16*r1216;
|
||||
dydt[Approx8Net::ihe4] += -yhe4*yne20*rne20a;
|
||||
|
||||
dydt[Approx8Net::ic12] = (1./6.)*yhe4*yhe4*yhe4*r3a;
|
||||
dydt[Approx8Net::ic12] += -yhe4*yc12*rc12a;
|
||||
dydt[Approx8Net::ic12] += -yh1*yc12*rc12p;
|
||||
dydt[Approx8Net::ic12] += yh1*yn14*aFrac*rn14p;
|
||||
dydt[Approx8Net::ic12] += -yc12*yc12*r1212;
|
||||
dydt[Approx8Net::ic12] += -yc12*yo16*r1216;
|
||||
|
||||
dydt[Approx8Net::in14] = yh1*yc12*rc12p;
|
||||
dydt[Approx8Net::in14] += -yh1*yn14*rn14p;
|
||||
dydt[Approx8Net::in14] += yh1*yo16*ro16p;
|
||||
dydt[Approx8Net::in14] += -yhe4*yn14*rn14a;
|
||||
|
||||
dydt[Approx8Net::io16] = yhe4*yc12*rc12a;
|
||||
dydt[Approx8Net::io16] += yh1*yn14*pFrac*rn14p;
|
||||
dydt[Approx8Net::io16] += -yh1*yo16*ro16p;
|
||||
dydt[Approx8Net::io16] += -yc12*yo16*r1216;
|
||||
dydt[Approx8Net::io16] += -yhe4*yo16*ro16a;
|
||||
|
||||
dydt[Approx8Net::ine20] = 0.5*yc12*yc12*r1212;
|
||||
dydt[Approx8Net::ine20] += yhe4*yn14*rn14a;
|
||||
dydt[Approx8Net::ine20] += yhe4*yo16*ro16a;
|
||||
dydt[Approx8Net::ine20] += -yhe4*yne20*rne20a;
|
||||
|
||||
dydt[Approx8Net::img24] = yc12*yo16*r1216;
|
||||
dydt[Approx8Net::img24] += yhe4*yne20*rne20a;
|
||||
|
||||
dydt[Approx8Net::iTemp] = 0.;
|
||||
dydt[Approx8Net::iDensity] = 0.;
|
||||
|
||||
// energy accounting
|
||||
double eNuc = 0.;
|
||||
for (int i=0; i<Approx8Net::nIso; i++) {
|
||||
eNuc += Approx8Net::mIon[i]*dydt[i];
|
||||
}
|
||||
dydt[Approx8Net::iEnergy] = -eNuc*avo*clight*clight;
|
||||
}
|
||||
|
||||
Approx8Network::Approx8Network() : Network(APPROX8) {}
|
||||
|
||||
NetOut Approx8Network::evaluate(const NetIn &netIn) {
|
||||
m_y = convert_netIn(netIn);
|
||||
m_tMax = netIn.tMax;
|
||||
m_dt0 = netIn.dt0;
|
||||
|
||||
const double stiff_abs_tol = m_config.get<double>("Network:Approx8:Stiff:AbsTol", 1.0e-6);
|
||||
const double stiff_rel_tol = m_config.get<double>("Network:Approx8:Stiff:RelTol", 1.0e-6);
|
||||
const double nonstiff_abs_tol = m_config.get<double>("Network:Approx8:NonStiff:AbsTol", 1.0e-6);
|
||||
const double nonstiff_rel_tol = m_config.get<double>("Network:Approx8:NonStiff:RelTol", 1.0e-6);
|
||||
|
||||
int num_steps = -1;
|
||||
|
||||
if (m_stiff) {
|
||||
LOG_DEBUG(m_logger, "Using stiff solver for Approx8Network");
|
||||
num_steps = integrate_const(
|
||||
make_dense_output<rosenbrock4<double>>(stiff_abs_tol, stiff_rel_tol),
|
||||
std::make_pair(ODE(), Jacobian()),
|
||||
m_y,
|
||||
0.0,
|
||||
m_tMax,
|
||||
m_dt0
|
||||
);
|
||||
|
||||
} else {
|
||||
LOG_DEBUG(m_logger, "Using non stiff solver for Approx8Network");
|
||||
num_steps = integrate_const (
|
||||
make_dense_output<runge_kutta_dopri5<vector_type>>(nonstiff_abs_tol, nonstiff_rel_tol),
|
||||
ODE(),
|
||||
m_y,
|
||||
0.0,
|
||||
m_tMax,
|
||||
m_dt0
|
||||
);
|
||||
}
|
||||
|
||||
double ySum = 0.0;
|
||||
for (int i = 0; i < Approx8Net::nIso; i++) {
|
||||
m_y[i] *= Approx8Net::aIon[i];
|
||||
ySum += m_y[i];
|
||||
}
|
||||
for (int i = 0; i < Approx8Net::nIso; i++) {
|
||||
m_y[i] /= ySum;
|
||||
}
|
||||
|
||||
NetOut netOut;
|
||||
std::vector<double> outComposition;
|
||||
outComposition.reserve(Approx8Net::nVar);
|
||||
|
||||
for (int i = 0; i < Approx8Net::nIso; i++) {
|
||||
outComposition.push_back(m_y[i]);
|
||||
}
|
||||
netOut.energy = m_y[Approx8Net::iEnergy];
|
||||
netOut.num_steps = num_steps;
|
||||
|
||||
const std::vector<std::string> symbols = {"H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"};
|
||||
netOut.composition = fourdst::composition::Composition(symbols, outComposition);
|
||||
|
||||
return netOut;
|
||||
}
|
||||
|
||||
void Approx8Network::setStiff(bool stiff) {
|
||||
m_stiff = stiff;
|
||||
}
|
||||
|
||||
vector_type Approx8Network::convert_netIn(const NetIn &netIn) {
|
||||
vector_type y(Approx8Net::nVar, 0.0);
|
||||
y[Approx8Net::ih1] = netIn.composition.getNumberFraction("H-1");
|
||||
y[Approx8Net::ihe3] = netIn.composition.getNumberFraction("He-3");
|
||||
y[Approx8Net::ihe4] = netIn.composition.getNumberFraction("He-4");
|
||||
y[Approx8Net::ic12] = netIn.composition.getNumberFraction("C-12");
|
||||
y[Approx8Net::in14] = netIn.composition.getNumberFraction("N-14");
|
||||
y[Approx8Net::io16] = netIn.composition.getNumberFraction("O-16");
|
||||
y[Approx8Net::ine20] = netIn.composition.getNumberFraction("Ne-20");
|
||||
y[Approx8Net::img24] = netIn.composition.getNumberFraction("Mg-24");
|
||||
y[Approx8Net::iTemp] = netIn.temperature;
|
||||
y[Approx8Net::iDensity] = netIn.density;
|
||||
y[Approx8Net::iEnergy] = netIn.energy;
|
||||
|
||||
return y;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
// main program
|
||||
|
||||
1147
src/lib/engine/engine_graph.cpp
Normal file
1147
src/lib/engine/engine_graph.cpp
Normal file
File diff suppressed because it is too large
Load Diff
113
src/lib/engine/procedures/construction.cpp
Normal file
113
src/lib/engine/procedures/construction.cpp
Normal file
@@ -0,0 +1,113 @@
|
||||
#include "gridfire/engine/procedures/construction.h"
|
||||
|
||||
#include <ranges>
|
||||
#include <stdexcept>
|
||||
|
||||
#include "gridfire/reaction/reaction.h"
|
||||
#include "gridfire/reaction/reaclib.h"
|
||||
|
||||
#include "fourdst/composition/composition.h"
|
||||
|
||||
#include "fourdst/logging/logging.h"
|
||||
|
||||
#include "quill/Logger.h"
|
||||
#include "quill/LogMacros.h"
|
||||
|
||||
namespace gridfire {
|
||||
using reaction::LogicalReactionSet;
|
||||
using reaction::ReactionSet;
|
||||
using reaction::Reaction;
|
||||
using fourdst::composition::Composition;
|
||||
using fourdst::atomic::Species;
|
||||
|
||||
|
||||
LogicalReactionSet build_reaclib_nuclear_network(
|
||||
const Composition &composition,
|
||||
BuildDepthType maxLayers,
|
||||
bool reverse
|
||||
) {
|
||||
int depth;
|
||||
if (std::holds_alternative<NetworkBuildDepth>(maxLayers)) {
|
||||
depth = static_cast<int>(std::get<NetworkBuildDepth>(maxLayers));
|
||||
} else {
|
||||
depth = std::get<int>(maxLayers);
|
||||
}
|
||||
auto logger = fourdst::logging::LogManager::getInstance().getLogger("log");
|
||||
if (depth == 0) {
|
||||
LOG_ERROR(logger, "Network build depth is set to 0. No reactions will be collected.");
|
||||
throw std::logic_error("Network build depth is set to 0. No reactions will be collected.");
|
||||
}
|
||||
|
||||
const auto allReactions = reaclib::get_all_reactions();
|
||||
std::vector<Reaction> remainingReactions;
|
||||
for (const auto& reaction : allReactions) {
|
||||
if (reaction.is_reverse() == reverse) {
|
||||
remainingReactions.push_back(reaction);
|
||||
}
|
||||
}
|
||||
|
||||
if (depth == static_cast<int>(NetworkBuildDepth::Full)) {
|
||||
LOG_INFO(logger, "Building full nuclear network with a total of {} reactions.", allReactions.size());
|
||||
const ReactionSet reactionSet(remainingReactions);
|
||||
return reaction::packReactionSetToLogicalReactionSet(reactionSet);
|
||||
}
|
||||
|
||||
std::unordered_set<Species> availableSpecies;
|
||||
for (const auto &entry: composition | std::views::values) {
|
||||
if (entry.mass_fraction() > 0.0) {
|
||||
availableSpecies.insert(entry.isotope());
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
std::vector<Reaction> collectedReactions;
|
||||
|
||||
LOG_INFO(logger, "Starting network construction with {} available species.", availableSpecies.size());
|
||||
for (int layer = 0; layer < depth && !remainingReactions.empty(); ++layer) {
|
||||
LOG_TRACE_L1(logger, "Collecting reactions for layer {} with {} remaining reactions. Currently there are {} available species", layer, remainingReactions.size(), availableSpecies.size());
|
||||
std::vector<Reaction> reactionsForNextPass;
|
||||
std::unordered_set<Species> newProductsThisLayer;
|
||||
bool newReactionsAdded = false;
|
||||
|
||||
reactionsForNextPass.reserve(remainingReactions.size());
|
||||
|
||||
for (const auto &reaction : remainingReactions) {
|
||||
bool allReactantsAvailable = true;
|
||||
for (const auto& reactant : reaction.reactants()) {
|
||||
if (!availableSpecies.contains(reactant)) {
|
||||
allReactantsAvailable = false;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if (allReactantsAvailable) {
|
||||
collectedReactions.push_back(reaction);
|
||||
newReactionsAdded = true;
|
||||
|
||||
for (const auto& product : reaction.products()) {
|
||||
newProductsThisLayer.insert(product);
|
||||
}
|
||||
} else {
|
||||
reactionsForNextPass.push_back(reaction);
|
||||
}
|
||||
}
|
||||
|
||||
if (!newReactionsAdded) {
|
||||
LOG_INFO(logger, "No new reactions added in layer {}. Stopping network construction with {} reactions collected.", layer, collectedReactions.size());
|
||||
break;
|
||||
}
|
||||
|
||||
LOG_TRACE_L1(logger, "Layer {}: Collected {} reactions. New products this layer: {}", layer, collectedReactions.size(), newProductsThisLayer.size());
|
||||
availableSpecies.insert(newProductsThisLayer.begin(), newProductsThisLayer.end());
|
||||
|
||||
remainingReactions = std::move(reactionsForNextPass);
|
||||
}
|
||||
|
||||
LOG_INFO(logger, "Network construction completed with {} reactions collected.", collectedReactions.size());
|
||||
const ReactionSet reactionSet(collectedReactions);
|
||||
return reaction::packReactionSetToLogicalReactionSet(reactionSet);
|
||||
}
|
||||
|
||||
|
||||
|
||||
}
|
||||
229
src/lib/engine/procedures/priming.cpp
Normal file
229
src/lib/engine/procedures/priming.cpp
Normal file
@@ -0,0 +1,229 @@
|
||||
#include "gridfire/engine/procedures/priming.h"
|
||||
#include "gridfire/engine/views/engine_priming.h"
|
||||
#include "gridfire/engine/procedures/construction.h"
|
||||
#include "gridfire/solver/solver.h"
|
||||
|
||||
#include "gridfire/engine/engine_abstract.h"
|
||||
#include "gridfire/network.h"
|
||||
|
||||
#include "fourdst/logging/logging.h"
|
||||
#include "quill/Logger.h"
|
||||
#include "quill/LogMacros.h"
|
||||
|
||||
namespace gridfire {
|
||||
using fourdst::composition::Composition;
|
||||
using fourdst::atomic::Species;
|
||||
|
||||
const reaction::Reaction* findDominantCreationChannel (
|
||||
const DynamicEngine& engine,
|
||||
const Species& species,
|
||||
const std::vector<double>& Y,
|
||||
const double T9,
|
||||
const double rho
|
||||
) {
|
||||
const reaction::Reaction* dominateReaction = nullptr;
|
||||
double maxFlow = -1.0;
|
||||
for (const auto& reaction : engine.getNetworkReactions()) {
|
||||
if (reaction.contains(species) && reaction.stoichiometry(species) > 0) {
|
||||
const double flow = engine.calculateMolarReactionFlow(reaction, Y, T9, rho);
|
||||
if (flow > maxFlow) {
|
||||
maxFlow = flow;
|
||||
dominateReaction = &reaction;
|
||||
}
|
||||
}
|
||||
}
|
||||
return dominateReaction;
|
||||
}
|
||||
|
||||
|
||||
PrimingReport primeNetwork(const NetIn& netIn, DynamicEngine& engine) {
|
||||
auto logger = LogManager::getInstance().getLogger("log");
|
||||
|
||||
std::vector<Species> speciesToPrime;
|
||||
for (const auto &entry: netIn.composition | std::views::values) {
|
||||
if (entry.mass_fraction() == 0.0) {
|
||||
speciesToPrime.push_back(entry.isotope());
|
||||
}
|
||||
}
|
||||
LOG_DEBUG(logger, "Priming {} species in the network.", speciesToPrime.size());
|
||||
|
||||
PrimingReport report;
|
||||
if (speciesToPrime.empty()) {
|
||||
report.primedComposition = netIn.composition;
|
||||
report.success = true;
|
||||
report.status = PrimingReportStatus::NO_SPECIES_TO_PRIME;
|
||||
return report;
|
||||
}
|
||||
|
||||
const double T9 = netIn.temperature / 1e9;
|
||||
const double rho = netIn.density;
|
||||
const auto initialReactionSet = engine.getNetworkReactions();
|
||||
|
||||
report.status = PrimingReportStatus::FULL_SUCCESS;
|
||||
report.success = true;
|
||||
|
||||
// --- 1: pack composition into internal map ---
|
||||
std::unordered_map<Species, double> currentMassFractions;
|
||||
for (const auto& entry : netIn.composition | std::views::values) {
|
||||
currentMassFractions[entry.isotope()] = entry.mass_fraction();
|
||||
}
|
||||
for (const auto& entry : speciesToPrime) {
|
||||
currentMassFractions[entry] = 0.0; // Initialize priming species with 0 mass fraction
|
||||
}
|
||||
|
||||
std::unordered_map<Species, double> totalMassFractionChanges;
|
||||
|
||||
engine.rebuild(netIn.composition, NetworkBuildDepth::Full);
|
||||
|
||||
for (const auto& primingSpecies : speciesToPrime) {
|
||||
LOG_TRACE_L3(logger, "Priming species: {}", primingSpecies.name());
|
||||
|
||||
// Create a temporary composition from the current internal state for the primer
|
||||
Composition tempComp;
|
||||
for(const auto& [sp, mf] : currentMassFractions) {
|
||||
tempComp.registerSymbol(std::string(sp.name()));
|
||||
if (mf < 0.0 && std::abs(mf) < 1e-16) {
|
||||
tempComp.setMassFraction(sp, 0.0); // Avoid negative mass fractions
|
||||
} else {
|
||||
tempComp.setMassFraction(sp, mf);
|
||||
}
|
||||
}
|
||||
tempComp.finalize(true); // Finalize with normalization
|
||||
|
||||
NetIn tempNetIn = netIn;
|
||||
tempNetIn.composition = tempComp;
|
||||
|
||||
NetworkPrimingEngineView primer(primingSpecies, engine);
|
||||
|
||||
if (primer.getNetworkReactions().size() == 0) {
|
||||
LOG_ERROR(logger, "No priming reactions found for species {}.", primingSpecies.name());
|
||||
report.success = false;
|
||||
report.status = PrimingReportStatus::FAILED_TO_FIND_PRIMING_REACTIONS;
|
||||
continue;
|
||||
}
|
||||
|
||||
const auto Y = primer.mapNetInToMolarAbundanceVector(tempNetIn);
|
||||
const double destructionRateConstant = calculateDestructionRateConstant(primer, primingSpecies, Y, T9, rho);
|
||||
|
||||
double equilibriumMassFraction = 0.0;
|
||||
|
||||
if (destructionRateConstant > 1e-99) {
|
||||
const double creationRate = calculateCreationRate(primer, primingSpecies, Y, T9, rho);
|
||||
equilibriumMassFraction = (creationRate / destructionRateConstant) * primingSpecies.mass();
|
||||
if (std::isnan(equilibriumMassFraction)) {
|
||||
LOG_WARNING(logger, "Equilibrium mass fraction for {} is NaN. Setting to 0.0. This is likely not an issue. It probably originates from all reactions leading to creation and destruction being frozen out. In that case 0.0 should be a good approximation. Hint: This happens often when the network temperature is very the low. ", primingSpecies.name());
|
||||
equilibriumMassFraction = 0.0;
|
||||
}
|
||||
LOG_TRACE_L3(logger, "Found equilibrium for {}: X_eq = {:.4e}", primingSpecies.name(), equilibriumMassFraction);
|
||||
|
||||
const reaction::Reaction* dominantChannel = findDominantCreationChannel(primer, primingSpecies, Y, T9, rho);
|
||||
|
||||
if (dominantChannel) {
|
||||
LOG_TRACE_L3(logger, "Dominant creation channel for {}: {}", primingSpecies.name(), dominantChannel->peName());
|
||||
|
||||
double totalReactantMass = 0.0;
|
||||
for (const auto& reactant : dominantChannel->reactants()) {
|
||||
totalReactantMass += reactant.mass();
|
||||
}
|
||||
|
||||
double scalingFactor = 1.0;
|
||||
for (const auto& reactant : dominantChannel->reactants()) {
|
||||
const double massToSubtract = equilibriumMassFraction * (reactant.mass() / totalReactantMass);
|
||||
double availableMass = 0.0;
|
||||
if (currentMassFractions.contains(reactant)) {
|
||||
availableMass = currentMassFractions.at(reactant);
|
||||
}
|
||||
if (massToSubtract > availableMass && availableMass > 0) {
|
||||
scalingFactor = std::min(scalingFactor, availableMass / massToSubtract);
|
||||
}
|
||||
}
|
||||
|
||||
if (scalingFactor < 1.0) {
|
||||
LOG_WARNING(logger, "Priming for {} was limited by reactant availability. Scaling transfer by {:.4e}", primingSpecies.name(), scalingFactor);
|
||||
equilibriumMassFraction *= scalingFactor;
|
||||
}
|
||||
|
||||
// Update the internal mass fraction map and accumulate total changes
|
||||
totalMassFractionChanges[primingSpecies] += equilibriumMassFraction;
|
||||
currentMassFractions[primingSpecies] += equilibriumMassFraction;
|
||||
|
||||
for (const auto& reactant : dominantChannel->reactants()) {
|
||||
const double massToSubtract = equilibriumMassFraction * (reactant.mass() / totalReactantMass);
|
||||
totalMassFractionChanges[reactant] -= massToSubtract;
|
||||
currentMassFractions[reactant] -= massToSubtract;
|
||||
}
|
||||
} else {
|
||||
LOG_ERROR(logger, "Failed to find dominant creation channel for {}.", primingSpecies.name());
|
||||
report.status = PrimingReportStatus::FAILED_TO_FIND_CREATION_CHANNEL;
|
||||
totalMassFractionChanges[primingSpecies] += 1e-40;
|
||||
currentMassFractions[primingSpecies] += 1e-40;
|
||||
}
|
||||
} else {
|
||||
LOG_WARNING(logger, "No destruction channel found for {}. Using fallback abundance.", primingSpecies.name());
|
||||
totalMassFractionChanges[primingSpecies] += 1e-40;
|
||||
currentMassFractions[primingSpecies] += 1e-40;
|
||||
report.status = PrimingReportStatus::BASE_NETWORK_TOO_SHALLOW;
|
||||
}
|
||||
}
|
||||
|
||||
// --- Final Composition Construction ---
|
||||
std::vector<std::string> final_symbols;
|
||||
std::vector<double> final_mass_fractions;
|
||||
for(const auto& [species, mass_fraction] : currentMassFractions) {
|
||||
final_symbols.push_back(std::string(species.name()));
|
||||
if (mass_fraction < 0.0 && std::abs(mass_fraction) < 1e-16) {
|
||||
final_mass_fractions.push_back(0.0); // Avoid negative mass fractions
|
||||
} else {
|
||||
final_mass_fractions.push_back(mass_fraction);
|
||||
}
|
||||
}
|
||||
|
||||
// Create the final composition object from the pre-normalized mass fractions
|
||||
Composition primedComposition(final_symbols, final_mass_fractions, true);
|
||||
|
||||
report.primedComposition = primedComposition;
|
||||
for (const auto& [species, change] : totalMassFractionChanges) {
|
||||
report.massFractionChanges.emplace_back(species, change);
|
||||
}
|
||||
|
||||
engine.setNetworkReactions(initialReactionSet);
|
||||
return report;
|
||||
}
|
||||
|
||||
double calculateDestructionRateConstant(
|
||||
const DynamicEngine& engine,
|
||||
const fourdst::atomic::Species& species,
|
||||
const std::vector<double>& Y,
|
||||
const double T9,
|
||||
const double rho
|
||||
) {
|
||||
const int speciesIndex = engine.getSpeciesIndex(species);
|
||||
std::vector<double> Y_scaled(Y.begin(), Y.end());
|
||||
Y_scaled[speciesIndex] = 1.0; // Set the abundance of the species to 1.0 for rate constant calculation
|
||||
double destructionRateConstant = 0.0;
|
||||
for (const auto& reaction: engine.getNetworkReactions()) {
|
||||
if (reaction.contains_reactant(species)) {
|
||||
const int stoichiometry = reaction.stoichiometry(species);
|
||||
destructionRateConstant += std::abs(stoichiometry) * engine.calculateMolarReactionFlow(reaction, Y_scaled, T9, rho);
|
||||
}
|
||||
}
|
||||
return destructionRateConstant;
|
||||
}
|
||||
|
||||
double calculateCreationRate(
|
||||
const DynamicEngine& engine,
|
||||
const fourdst::atomic::Species& species,
|
||||
const std::vector<double>& Y,
|
||||
const double T9,
|
||||
const double rho
|
||||
) {
|
||||
double creationRate = 0.0;
|
||||
for (const auto& reaction: engine.getNetworkReactions()) {
|
||||
const int stoichiometry = reaction.stoichiometry(species);
|
||||
if (stoichiometry > 0) {
|
||||
creationRate += stoichiometry * engine.calculateMolarReactionFlow(reaction, Y, T9, rho);
|
||||
}
|
||||
}
|
||||
return creationRate;
|
||||
}
|
||||
}
|
||||
673
src/lib/engine/views/engine_adaptive.cpp
Normal file
673
src/lib/engine/views/engine_adaptive.cpp
Normal file
@@ -0,0 +1,673 @@
|
||||
#include "gridfire/engine/views/engine_adaptive.h"
|
||||
|
||||
#include <ranges>
|
||||
#include <queue>
|
||||
#include <algorithm>
|
||||
|
||||
|
||||
#include "gridfire/network.h"
|
||||
#include "gridfire/exceptions/error_engine.h"
|
||||
|
||||
#include "quill/LogMacros.h"
|
||||
#include "quill/Logger.h"
|
||||
|
||||
namespace gridfire {
|
||||
using fourdst::atomic::Species;
|
||||
AdaptiveEngineView::AdaptiveEngineView(
|
||||
DynamicEngine &baseEngine
|
||||
) :
|
||||
m_baseEngine(baseEngine),
|
||||
m_activeSpecies(baseEngine.getNetworkSpecies()),
|
||||
m_activeReactions(baseEngine.getNetworkReactions()),
|
||||
m_speciesIndexMap(constructSpeciesIndexMap()),
|
||||
m_reactionIndexMap(constructReactionIndexMap())
|
||||
{
|
||||
}
|
||||
|
||||
std::vector<size_t> AdaptiveEngineView::constructSpeciesIndexMap() const {
|
||||
LOG_TRACE_L1(m_logger, "Constructing species index map for adaptive engine view...");
|
||||
std::unordered_map<Species, size_t> fullSpeciesReverseMap;
|
||||
const auto& fullSpeciesList = m_baseEngine.getNetworkSpecies();
|
||||
|
||||
fullSpeciesReverseMap.reserve(fullSpeciesList.size());
|
||||
|
||||
for (size_t i = 0; i < fullSpeciesList.size(); ++i) {
|
||||
fullSpeciesReverseMap[fullSpeciesList[i]] = i;
|
||||
}
|
||||
|
||||
std::vector<size_t> speciesIndexMap;
|
||||
speciesIndexMap.reserve(m_activeSpecies.size());
|
||||
|
||||
for (const auto& active_species : m_activeSpecies) {
|
||||
auto it = fullSpeciesReverseMap.find(active_species);
|
||||
if (it != fullSpeciesReverseMap.end()) {
|
||||
speciesIndexMap.push_back(it->second);
|
||||
} else {
|
||||
LOG_ERROR(m_logger, "Species '{}' not found in full species map.", active_species.name());
|
||||
m_logger -> flush_log();
|
||||
throw std::runtime_error("Species not found in full species map: " + std::string(active_species.name()));
|
||||
}
|
||||
}
|
||||
LOG_TRACE_L1(m_logger, "Species index map constructed with {} entries.", speciesIndexMap.size());
|
||||
return speciesIndexMap;
|
||||
|
||||
}
|
||||
|
||||
std::vector<size_t> AdaptiveEngineView::constructReactionIndexMap() const {
|
||||
LOG_TRACE_L1(m_logger, "Constructing reaction index map for adaptive engine view...");
|
||||
|
||||
// --- Step 1: Create a reverse map using the reaction's unique ID as the key. ---
|
||||
std::unordered_map<std::string_view, size_t> fullReactionReverseMap;
|
||||
const auto& fullReactionSet = m_baseEngine.getNetworkReactions();
|
||||
fullReactionReverseMap.reserve(fullReactionSet.size());
|
||||
|
||||
for (size_t i_full = 0; i_full < fullReactionSet.size(); ++i_full) {
|
||||
fullReactionReverseMap[fullReactionSet[i_full].id()] = i_full;
|
||||
}
|
||||
|
||||
// --- Step 2: Build the final index map using the active reaction set. ---
|
||||
std::vector<size_t> reactionIndexMap;
|
||||
reactionIndexMap.reserve(m_activeReactions.size());
|
||||
|
||||
for (const auto& active_reaction_ptr : m_activeReactions) {
|
||||
auto it = fullReactionReverseMap.find(active_reaction_ptr.id());
|
||||
|
||||
if (it != fullReactionReverseMap.end()) {
|
||||
reactionIndexMap.push_back(it->second);
|
||||
} else {
|
||||
LOG_ERROR(m_logger, "Active reaction '{}' not found in base engine during reaction index map construction.", active_reaction_ptr.id());
|
||||
m_logger->flush_log();
|
||||
throw std::runtime_error("Mismatch between active reactions and base engine.");
|
||||
}
|
||||
}
|
||||
|
||||
LOG_TRACE_L1(m_logger, "Reaction index map constructed with {} entries.", reactionIndexMap.size());
|
||||
return reactionIndexMap;
|
||||
}
|
||||
|
||||
fourdst::composition::Composition AdaptiveEngineView::update(const NetIn &netIn) {
|
||||
fourdst::composition::Composition baseUpdatedComposition = m_baseEngine.update(netIn);
|
||||
NetIn updatedNetIn = netIn;
|
||||
|
||||
// for (const auto &entry: netIn.composition | std::views::values) {
|
||||
// if (baseUpdatedComposition.contains(entry.isotope())) {
|
||||
// updatedNetIn.composition.setMassFraction(entry.isotope(), baseUpdatedComposition.getMassFraction(entry.isotope()));
|
||||
// }
|
||||
// }
|
||||
updatedNetIn.composition = baseUpdatedComposition;
|
||||
|
||||
updatedNetIn.composition.finalize(false);
|
||||
|
||||
LOG_TRACE_L1(m_logger, "Updating AdaptiveEngineView with new network input...");
|
||||
|
||||
std::vector<double> Y_Full;
|
||||
std::vector<ReactionFlow> allFlows = calculateAllReactionFlows(updatedNetIn, Y_Full);
|
||||
|
||||
double maxFlow = 0.0;
|
||||
|
||||
for (const auto&[reactionPtr, flowRate]: allFlows) {
|
||||
if (flowRate > maxFlow) {
|
||||
maxFlow = flowRate;
|
||||
}
|
||||
}
|
||||
LOG_DEBUG(m_logger, "Maximum reaction flow rate in adaptive engine view: {:0.3E} [mol/s]", maxFlow);
|
||||
|
||||
const std::unordered_set<Species> reachableSpecies = findReachableSpecies(updatedNetIn);
|
||||
LOG_DEBUG(m_logger, "Found {} reachable species in adaptive engine view.", reachableSpecies.size());
|
||||
|
||||
const std::vector<const reaction::LogicalReaction*> finalReactions = cullReactionsByFlow(allFlows, reachableSpecies, Y_Full, maxFlow);
|
||||
|
||||
finalizeActiveSet(finalReactions);
|
||||
|
||||
auto [rescuedReactions, rescuedSpecies] = rescueEdgeSpeciesDestructionChannel(Y_Full, netIn.temperature/1e9, netIn.density, m_activeSpecies, m_activeReactions);
|
||||
|
||||
for (const auto& reactionPtr : rescuedReactions) {
|
||||
m_activeReactions.add_reaction(*reactionPtr);
|
||||
}
|
||||
|
||||
for (const auto& species : rescuedSpecies) {
|
||||
if (!std::ranges::contains(m_activeSpecies, species)) {
|
||||
m_activeSpecies.push_back(species);
|
||||
}
|
||||
}
|
||||
|
||||
m_speciesIndexMap = constructSpeciesIndexMap();
|
||||
m_reactionIndexMap = constructReactionIndexMap();
|
||||
|
||||
m_isStale = false;
|
||||
|
||||
LOG_INFO(m_logger, "AdaptiveEngineView updated successfully with {} active species and {} active reactions.", m_activeSpecies.size(), m_activeReactions.size());
|
||||
|
||||
return updatedNetIn.composition;
|
||||
}
|
||||
|
||||
bool AdaptiveEngineView::isStale(const NetIn &netIn) {
|
||||
return m_isStale || m_baseEngine.isStale(netIn);
|
||||
}
|
||||
|
||||
const std::vector<Species> & AdaptiveEngineView::getNetworkSpecies() const {
|
||||
return m_activeSpecies;
|
||||
}
|
||||
|
||||
std::expected<StepDerivatives<double>, expectations::StaleEngineError> AdaptiveEngineView::calculateRHSAndEnergy(
|
||||
const std::vector<double> &Y_culled,
|
||||
const double T9,
|
||||
const double rho
|
||||
) const {
|
||||
validateState();
|
||||
|
||||
const auto Y_full = mapCulledToFull(Y_culled);
|
||||
|
||||
auto result = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho);
|
||||
|
||||
if (!result) {
|
||||
return std::unexpected{result.error()};
|
||||
}
|
||||
|
||||
const auto [dydt, nuclearEnergyGenerationRate] = result.value();
|
||||
StepDerivatives<double> culledResults;
|
||||
culledResults.nuclearEnergyGenerationRate = nuclearEnergyGenerationRate;
|
||||
culledResults.dydt = mapFullToCulled(dydt);
|
||||
|
||||
return culledResults;
|
||||
}
|
||||
|
||||
void AdaptiveEngineView::generateJacobianMatrix(
|
||||
const std::vector<double> &Y_dynamic,
|
||||
const double T9,
|
||||
const double rho
|
||||
) const {
|
||||
validateState();
|
||||
const auto Y_full = mapCulledToFull(Y_dynamic);
|
||||
|
||||
m_baseEngine.generateJacobianMatrix(Y_full, T9, rho);
|
||||
}
|
||||
|
||||
double AdaptiveEngineView::getJacobianMatrixEntry(
|
||||
const int i_culled,
|
||||
const int j_culled
|
||||
) const {
|
||||
validateState();
|
||||
const size_t i_full = mapCulledToFullSpeciesIndex(i_culled);
|
||||
const size_t j_full = mapCulledToFullSpeciesIndex(j_culled);
|
||||
|
||||
return m_baseEngine.getJacobianMatrixEntry(i_full, j_full);
|
||||
}
|
||||
|
||||
void AdaptiveEngineView::generateStoichiometryMatrix() {
|
||||
validateState();
|
||||
m_baseEngine.generateStoichiometryMatrix();
|
||||
}
|
||||
|
||||
int AdaptiveEngineView::getStoichiometryMatrixEntry(
|
||||
const int speciesIndex_culled,
|
||||
const int reactionIndex_culled
|
||||
) const {
|
||||
validateState();
|
||||
const size_t speciesIndex_full = mapCulledToFullSpeciesIndex(speciesIndex_culled);
|
||||
const size_t reactionIndex_full = mapCulledToFullReactionIndex(reactionIndex_culled);
|
||||
return m_baseEngine.getStoichiometryMatrixEntry(speciesIndex_full, reactionIndex_full);
|
||||
}
|
||||
|
||||
double AdaptiveEngineView::calculateMolarReactionFlow(
|
||||
const reaction::Reaction &reaction,
|
||||
const std::vector<double> &Y_culled,
|
||||
const double T9,
|
||||
const double rho
|
||||
) const {
|
||||
validateState();
|
||||
if (!m_activeReactions.contains(reaction)) {
|
||||
LOG_ERROR(m_logger, "Reaction '{}' is not part of the active reactions in the adaptive engine view.", reaction.id());
|
||||
m_logger -> flush_log();
|
||||
throw std::runtime_error("Reaction not found in active reactions: " + std::string(reaction.id()));
|
||||
}
|
||||
const auto Y = mapCulledToFull(Y_culled);
|
||||
|
||||
return m_baseEngine.calculateMolarReactionFlow(reaction, Y, T9, rho);
|
||||
}
|
||||
|
||||
const reaction::LogicalReactionSet & AdaptiveEngineView::getNetworkReactions() const {
|
||||
return m_activeReactions;
|
||||
}
|
||||
|
||||
void AdaptiveEngineView::setNetworkReactions(const reaction::LogicalReactionSet &reactions) {
|
||||
LOG_CRITICAL(m_logger, "AdaptiveEngineView does not support setting network reactions directly. Use update() with NetIn instead. Perhaps you meant to call this on the base engine?");
|
||||
throw exceptions::UnableToSetNetworkReactionsError("AdaptiveEngineView does not support setting network reactions directly. Use update() with NetIn instead. Perhaps you meant to call this on the base engine?");
|
||||
}
|
||||
|
||||
std::expected<std::unordered_map<Species, double>, expectations::StaleEngineError> AdaptiveEngineView::getSpeciesTimescales(
|
||||
const std::vector<double> &Y_culled,
|
||||
const double T9,
|
||||
const double rho
|
||||
) const {
|
||||
validateState();
|
||||
const auto Y_full = mapCulledToFull(Y_culled);
|
||||
const auto result = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho);
|
||||
|
||||
if (!result) {
|
||||
return std::unexpected{result.error()};
|
||||
}
|
||||
|
||||
const std::unordered_map<Species, double> fullTimescales = result.value();
|
||||
|
||||
|
||||
std::unordered_map<Species, double> culledTimescales;
|
||||
culledTimescales.reserve(m_activeSpecies.size());
|
||||
for (const auto& active_species : m_activeSpecies) {
|
||||
if (fullTimescales.contains(active_species)) {
|
||||
culledTimescales[active_species] = fullTimescales.at(active_species);
|
||||
}
|
||||
}
|
||||
return culledTimescales;
|
||||
|
||||
}
|
||||
|
||||
std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError>
|
||||
AdaptiveEngineView::getSpeciesDestructionTimescales(
|
||||
const std::vector<double> &Y,
|
||||
double T9,
|
||||
double rho
|
||||
) const {
|
||||
validateState();
|
||||
|
||||
const auto Y_full = mapCulledToFull(Y);
|
||||
const auto result = m_baseEngine.getSpeciesDestructionTimescales(Y_full, T9, rho);
|
||||
|
||||
if (!result) {
|
||||
return std::unexpected{result.error()};
|
||||
}
|
||||
|
||||
const std::unordered_map<Species, double> destructionTimescales = result.value();
|
||||
|
||||
std::unordered_map<Species, double> culledTimescales;
|
||||
culledTimescales.reserve(m_activeSpecies.size());
|
||||
for (const auto& active_species : m_activeSpecies) {
|
||||
if (destructionTimescales.contains(active_species)) {
|
||||
culledTimescales[active_species] = destructionTimescales.at(active_species);
|
||||
}
|
||||
}
|
||||
return culledTimescales;
|
||||
}
|
||||
|
||||
void AdaptiveEngineView::setScreeningModel(const screening::ScreeningType model) {
|
||||
m_baseEngine.setScreeningModel(model);
|
||||
}
|
||||
|
||||
screening::ScreeningType AdaptiveEngineView::getScreeningModel() const {
|
||||
return m_baseEngine.getScreeningModel();
|
||||
}
|
||||
|
||||
std::vector<double> AdaptiveEngineView::mapNetInToMolarAbundanceVector(const NetIn &netIn) const {
|
||||
std::vector<double> Y(m_activeSpecies.size(), 0.0); // Initialize with zeros
|
||||
for (const auto& [symbol, entry] : netIn.composition) {
|
||||
Y[getSpeciesIndex(entry.isotope())] = netIn.composition.getMolarAbundance(symbol); // Map species to their molar abundance
|
||||
}
|
||||
return Y; // Return the vector of molar abundances
|
||||
}
|
||||
|
||||
PrimingReport AdaptiveEngineView::primeEngine(const NetIn &netIn) {
|
||||
return m_baseEngine.primeEngine(netIn);
|
||||
}
|
||||
|
||||
int AdaptiveEngineView::getSpeciesIndex(const fourdst::atomic::Species &species) const {
|
||||
const auto it = std::ranges::find(m_activeSpecies, species);
|
||||
if (it != m_activeSpecies.end()) {
|
||||
return static_cast<int>(std::distance(m_activeSpecies.begin(), it));
|
||||
} else {
|
||||
LOG_ERROR(m_logger, "Species '{}' not found in active species list.", species.name());
|
||||
m_logger->flush_log();
|
||||
throw std::runtime_error("Species not found in active species list: " + std::string(species.name()));
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<double> AdaptiveEngineView::mapCulledToFull(const std::vector<double>& culled) const {
|
||||
std::vector<double> full(m_baseEngine.getNetworkSpecies().size(), 0.0);
|
||||
for (size_t i_culled = 0; i_culled < culled.size(); ++i_culled) {
|
||||
const size_t i_full = m_speciesIndexMap[i_culled];
|
||||
full[i_full] += culled[i_culled];
|
||||
}
|
||||
return full;
|
||||
}
|
||||
|
||||
std::vector<double> AdaptiveEngineView::mapFullToCulled(const std::vector<double>& full) const {
|
||||
std::vector<double> culled(m_activeSpecies.size(), 0.0);
|
||||
for (size_t i_culled = 0; i_culled < m_activeSpecies.size(); ++i_culled) {
|
||||
const size_t i_full = m_speciesIndexMap[i_culled];
|
||||
culled[i_culled] = full[i_full];
|
||||
}
|
||||
return culled;
|
||||
}
|
||||
|
||||
size_t AdaptiveEngineView::mapCulledToFullSpeciesIndex(size_t culledSpeciesIndex) const {
|
||||
if (culledSpeciesIndex < 0 || culledSpeciesIndex >= static_cast<int>(m_speciesIndexMap.size())) {
|
||||
LOG_ERROR(m_logger, "Culled index {} is out of bounds for species index map of size {}.", culledSpeciesIndex, m_speciesIndexMap.size());
|
||||
m_logger->flush_log();
|
||||
throw std::out_of_range("Culled index " + std::to_string(culledSpeciesIndex) + " is out of bounds for species index map of size " + std::to_string(m_speciesIndexMap.size()) + ".");
|
||||
}
|
||||
return m_speciesIndexMap[culledSpeciesIndex];
|
||||
}
|
||||
|
||||
size_t AdaptiveEngineView::mapCulledToFullReactionIndex(size_t culledReactionIndex) const {
|
||||
if (culledReactionIndex < 0 || culledReactionIndex >= static_cast<int>(m_reactionIndexMap.size())) {
|
||||
LOG_ERROR(m_logger, "Culled index {} is out of bounds for reaction index map of size {}.", culledReactionIndex, m_reactionIndexMap.size());
|
||||
m_logger->flush_log();
|
||||
throw std::out_of_range("Culled index " + std::to_string(culledReactionIndex) + " is out of bounds for reaction index map of size " + std::to_string(m_reactionIndexMap.size()) + ".");
|
||||
}
|
||||
return m_reactionIndexMap[culledReactionIndex];
|
||||
}
|
||||
|
||||
void AdaptiveEngineView::validateState() const {
|
||||
if (m_isStale) {
|
||||
LOG_ERROR(m_logger, "AdaptiveEngineView is stale. Please call update() before calculating RHS and energy.");
|
||||
m_logger->flush_log();
|
||||
throw std::runtime_error("AdaptiveEngineView is stale. Please call update() before calculating RHS and energy.");
|
||||
}
|
||||
}
|
||||
|
||||
// TODO: Change this to use a return value instead of an output parameter.
|
||||
std::vector<AdaptiveEngineView::ReactionFlow> AdaptiveEngineView::calculateAllReactionFlows(
|
||||
const NetIn &netIn,
|
||||
std::vector<double> &out_Y_Full
|
||||
) const {
|
||||
const auto& fullSpeciesList = m_baseEngine.getNetworkSpecies();
|
||||
out_Y_Full.clear();
|
||||
out_Y_Full.reserve(fullSpeciesList.size());
|
||||
|
||||
for (const auto& species: fullSpeciesList) {
|
||||
if (netIn.composition.contains(species)) {
|
||||
out_Y_Full.push_back(netIn.composition.getMolarAbundance(std::string(species.name())));
|
||||
} else {
|
||||
LOG_TRACE_L2(m_logger, "Species '{}' not found in composition. Setting abundance to 0.0.", species.name());
|
||||
out_Y_Full.push_back(0.0);
|
||||
}
|
||||
}
|
||||
|
||||
const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9)
|
||||
const double rho = netIn.density; // Density in g/cm^3
|
||||
|
||||
std::vector<ReactionFlow> reactionFlows;
|
||||
const auto& fullReactionSet = m_baseEngine.getNetworkReactions();
|
||||
reactionFlows.reserve(fullReactionSet.size());
|
||||
for (const auto& reaction : fullReactionSet) {
|
||||
const double flow = m_baseEngine.calculateMolarReactionFlow(reaction, out_Y_Full, T9, rho);
|
||||
reactionFlows.push_back({&reaction, flow});
|
||||
LOG_TRACE_L1(m_logger, "Reaction '{}' has flow rate: {:0.3E} [mol/s/g]", reaction.id(), flow);
|
||||
}
|
||||
return reactionFlows;
|
||||
}
|
||||
|
||||
std::unordered_set<Species> AdaptiveEngineView::findReachableSpecies(
|
||||
const NetIn &netIn
|
||||
) const {
|
||||
std::unordered_set<Species> reachable;
|
||||
std::queue<Species> to_vist;
|
||||
|
||||
constexpr double ABUNDANCE_FLOOR = 1e-12; // Abundance floor for a species to be considered part of the initial fuel
|
||||
for (const auto& species: m_baseEngine.getNetworkSpecies()) {
|
||||
if (netIn.composition.contains(species) && netIn.composition.getMassFraction(std::string(species.name())) > ABUNDANCE_FLOOR) {
|
||||
if (!reachable.contains(species)) {
|
||||
to_vist.push(species);
|
||||
reachable.insert(species);
|
||||
LOG_TRACE_L2(m_logger, "Network Connectivity Analysis: Species '{}' is part of the initial fuel.", species.name());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
bool new_species_found_in_pass = true;
|
||||
while (new_species_found_in_pass) {
|
||||
new_species_found_in_pass = false;
|
||||
for (const auto& reaction: m_baseEngine.getNetworkReactions()) {
|
||||
bool all_reactants_reachable = true;
|
||||
for (const auto& reactant: reaction.reactants()) {
|
||||
if (!reachable.contains(reactant)) {
|
||||
all_reactants_reachable = false;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (all_reactants_reachable) {
|
||||
for (const auto& product: reaction.products()) {
|
||||
if (!reachable.contains(product)) {
|
||||
reachable.insert(product);
|
||||
new_species_found_in_pass = true;
|
||||
LOG_TRACE_L2(m_logger, "Network Connectivity Analysis: Species '{}' is reachable via reaction '{}'.", product.name(), reaction.id());
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return reachable;
|
||||
}
|
||||
|
||||
std::vector<const reaction::LogicalReaction *> AdaptiveEngineView::cullReactionsByFlow(
|
||||
const std::vector<ReactionFlow> &allFlows,
|
||||
const std::unordered_set<fourdst::atomic::Species> &reachableSpecies,
|
||||
const std::vector<double> &Y_full,
|
||||
const double maxFlow
|
||||
) const {
|
||||
LOG_TRACE_L1(m_logger, "Culling reactions based on flow rates...");
|
||||
const auto relative_culling_threshold = m_config.get<double>("gridfire:AdaptiveEngineView:RelativeCullingThreshold", 1e-75);
|
||||
double absoluteCullingThreshold = relative_culling_threshold * maxFlow;
|
||||
LOG_DEBUG(m_logger, "Relative culling threshold: {:0.3E} ({})", relative_culling_threshold, absoluteCullingThreshold);
|
||||
std::vector<const reaction::LogicalReaction*> culledReactions;
|
||||
for (const auto& [reactionPtr, flowRate]: allFlows) {
|
||||
bool keepReaction = false;
|
||||
if (flowRate > absoluteCullingThreshold) {
|
||||
LOG_TRACE_L2(m_logger, "Maintaining reaction '{}' with relative (abs) flow rate: {:0.3E} ({:0.3E} [mol/s])", reactionPtr->id(), flowRate/maxFlow, flowRate);
|
||||
keepReaction = true;
|
||||
} else {
|
||||
bool zero_flow_due_to_reachable_reactants = false;
|
||||
if (flowRate < 1e-99 && flowRate > 0.0) {
|
||||
for (const auto& reactant: reactionPtr->reactants()) {
|
||||
const auto it = std::ranges::find(m_baseEngine.getNetworkSpecies(), reactant);
|
||||
const size_t index = std::distance(m_baseEngine.getNetworkSpecies().begin(), it);
|
||||
if (Y_full[index] < 1e-99 && reachableSpecies.contains(reactant)) {
|
||||
LOG_TRACE_L1(m_logger, "Maintaining reaction '{}' with low flow ({:0.3E} [mol/s/g]) due to reachable reactant '{}'.", reactionPtr->id(), flowRate, reactant.name());
|
||||
zero_flow_due_to_reachable_reactants = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
if (zero_flow_due_to_reachable_reactants) {
|
||||
keepReaction = true;
|
||||
}
|
||||
}
|
||||
if (keepReaction) {
|
||||
culledReactions.push_back(reactionPtr);
|
||||
} else {
|
||||
LOG_TRACE_L1(m_logger, "Culling reaction '{}' due to low flow rate or lack of connectivity.", reactionPtr->id());
|
||||
}
|
||||
}
|
||||
LOG_DEBUG(m_logger, "Selected {} (total: {}, culled: {}) reactions based on flow rates.", culledReactions.size(), allFlows.size(), allFlows.size() - culledReactions.size());
|
||||
return culledReactions;
|
||||
}
|
||||
|
||||
AdaptiveEngineView::RescueSet AdaptiveEngineView::rescueEdgeSpeciesDestructionChannel(
|
||||
const std::vector<double> &Y_full,
|
||||
const double T9,
|
||||
const double rho,
|
||||
const std::vector<Species> &activeSpecies,
|
||||
const reaction::LogicalReactionSet &activeReactions
|
||||
) const {
|
||||
const auto result = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho);
|
||||
if (!result) {
|
||||
LOG_ERROR(m_logger, "Failed to get species timescales due to stale engine state.");
|
||||
throw exceptions::StaleEngineError("Failed to get species timescales");
|
||||
}
|
||||
std::unordered_map<Species, double> timescales = result.value();
|
||||
std::set<Species> onlyProducedSpecies;
|
||||
for (const auto& reaction : activeReactions) {
|
||||
const std::vector<Species> products = reaction.products();
|
||||
onlyProducedSpecies.insert(products.begin(), products.end());
|
||||
}
|
||||
|
||||
// Remove species that are consumed by any one of the active reactions.
|
||||
std::erase_if(
|
||||
onlyProducedSpecies,
|
||||
[&](const Species &species) {
|
||||
for (const auto& reaction : activeReactions) {
|
||||
if (reaction.contains_reactant(species)) {
|
||||
return true; // If any active reaction consumes the species then erase it from the set.
|
||||
}
|
||||
}
|
||||
return false;
|
||||
}
|
||||
);
|
||||
|
||||
// Remove species that have a non-zero timescale (these are expected to be +inf as they should be the equilibrium species if using with a MultiscalePartitioningEngineView)
|
||||
std::erase_if(
|
||||
onlyProducedSpecies,
|
||||
[&](const Species &species) {
|
||||
return std::isinf(timescales.at(species));
|
||||
}
|
||||
);
|
||||
|
||||
LOG_TRACE_L1(
|
||||
m_logger,
|
||||
"Found {} species {}that are produced but not consumed in any reaction and do not have an inf timescale (those are expected to be equilibrium species and do not contribute to the stiffness of the network).",
|
||||
onlyProducedSpecies.size(),
|
||||
[&]() -> std::string {
|
||||
std::ostringstream ss;
|
||||
if (onlyProducedSpecies.empty()) {
|
||||
return "";
|
||||
}
|
||||
int count = 0;
|
||||
ss << "(";
|
||||
for (const auto& species : onlyProducedSpecies) {
|
||||
ss << species.name();
|
||||
if (count < onlyProducedSpecies.size() - 1) {
|
||||
ss << ", ";
|
||||
}
|
||||
count++;
|
||||
}
|
||||
ss << ") ";
|
||||
return ss.str();
|
||||
}()
|
||||
);
|
||||
|
||||
std::unordered_map<Species, const reaction::LogicalReaction*> reactionsToRescue;
|
||||
for (const auto& species : onlyProducedSpecies) {
|
||||
double maxSpeciesConsumptionRate = 0.0;
|
||||
for (const auto& reaction : m_baseEngine.getNetworkReactions()) {
|
||||
const bool speciesToCheckIsConsumed = reaction.contains_reactant(species);
|
||||
if (!speciesToCheckIsConsumed) {
|
||||
continue; // If the species is not consumed by this reaction, skip it.
|
||||
}
|
||||
bool allOtherReactantsAreAvailable = true;
|
||||
for (const auto& reactant : reaction.reactants()) {
|
||||
const bool reactantIsAvailable = std::ranges::contains(activeSpecies, reactant);
|
||||
if (!reactantIsAvailable && reactant != species) {
|
||||
allOtherReactantsAreAvailable = false;
|
||||
}
|
||||
}
|
||||
if (allOtherReactantsAreAvailable && speciesToCheckIsConsumed) {
|
||||
double rate = reaction.calculate_rate(T9);
|
||||
if (rate > maxSpeciesConsumptionRate) {
|
||||
maxSpeciesConsumptionRate = rate;
|
||||
reactionsToRescue[species] = &reaction;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
LOG_TRACE_L1(
|
||||
m_logger,
|
||||
"Rescuing {} {}reactions",
|
||||
reactionsToRescue.size(),
|
||||
[&]() -> std::string {
|
||||
std::ostringstream ss;
|
||||
if (reactionsToRescue.empty()) {
|
||||
return "";
|
||||
}
|
||||
int count = 0;
|
||||
ss << "(";
|
||||
for (const auto &reaction : reactionsToRescue | std::views::values) {
|
||||
ss << reaction->id();
|
||||
if (count < reactionsToRescue.size() - 1) {
|
||||
ss << ", ";
|
||||
}
|
||||
count++;
|
||||
}
|
||||
ss << ") ";
|
||||
return ss.str();
|
||||
}()
|
||||
);
|
||||
|
||||
LOG_TRACE_L1(
|
||||
m_logger,
|
||||
"Timescale adjustments due to reaction rescue: {}",
|
||||
[&]() -> std::string {
|
||||
std::stringstream ss;
|
||||
if (reactionsToRescue.empty()) {
|
||||
return "No reactions rescued...";
|
||||
}
|
||||
int count = 0;
|
||||
for (const auto& [species, reaction] : reactionsToRescue) {
|
||||
ss << "(Species: " << species.name() << " started with a timescale of " << timescales.at(species);
|
||||
ss << ", rescued by reaction: " << reaction->id();
|
||||
ss << " whose product timescales are --- [";
|
||||
int icount = 0;
|
||||
for (const auto& product : reaction->products()) {
|
||||
ss << product.name() << ": " << timescales.at(product);
|
||||
if (icount < reaction->products().size() - 1) {
|
||||
ss << ", ";
|
||||
}
|
||||
icount++;
|
||||
}
|
||||
ss << "])";
|
||||
|
||||
if (count < reactionsToRescue.size() - 1) {
|
||||
ss << ", ";
|
||||
}
|
||||
count++;
|
||||
}
|
||||
|
||||
return ss.str();
|
||||
}()
|
||||
);
|
||||
|
||||
RescueSet rescueSet;
|
||||
std::unordered_set<const reaction::LogicalReaction*> newReactions;
|
||||
std::unordered_set<Species> newSpecies;
|
||||
|
||||
for (const auto &reactionPtr: reactionsToRescue | std::views::values) {
|
||||
newReactions.insert(reactionPtr);
|
||||
for (const auto& product : reactionPtr->products()) {
|
||||
newSpecies.insert(product);
|
||||
}
|
||||
}
|
||||
return {std::move(newReactions), std::move(newSpecies)};
|
||||
}
|
||||
|
||||
void AdaptiveEngineView::finalizeActiveSet(
|
||||
const std::vector<const reaction::LogicalReaction *> &finalReactions
|
||||
) {
|
||||
std::unordered_set<Species>finalSpeciesSet;
|
||||
m_activeReactions.clear();
|
||||
for (const auto* reactionPtr: finalReactions) {
|
||||
m_activeReactions.add_reaction(*reactionPtr);
|
||||
for (const auto& reactant : reactionPtr->reactants()) {
|
||||
if (!finalSpeciesSet.contains(reactant)) {
|
||||
LOG_TRACE_L1(m_logger, "Adding reactant '{}' to active species set through reaction {}.", reactant.name(), reactionPtr->id());
|
||||
} else {
|
||||
LOG_TRACE_L1(m_logger, "Reactant '{}' already in active species set through another reaction.", reactant.name());
|
||||
}
|
||||
finalSpeciesSet.insert(reactant);
|
||||
}
|
||||
for (const auto& product : reactionPtr->products()) {
|
||||
if (!finalSpeciesSet.contains(product)) {
|
||||
LOG_TRACE_L1(m_logger, "Adding product '{}' to active species set through reaction {}.", product.name(), reactionPtr->id());
|
||||
} else {
|
||||
LOG_TRACE_L1(m_logger, "Product '{}' already in active species set through another reaction.", product.name());
|
||||
}
|
||||
finalSpeciesSet.insert(product);
|
||||
}
|
||||
}
|
||||
m_activeSpecies.clear();
|
||||
m_activeSpecies = std::vector<Species>(finalSpeciesSet.begin(), finalSpeciesSet.end());
|
||||
std::ranges::sort(
|
||||
m_activeSpecies,
|
||||
[](const Species &a, const Species &b) { return a.mass() < b.mass(); }
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
387
src/lib/engine/views/engine_defined.cpp
Normal file
387
src/lib/engine/views/engine_defined.cpp
Normal file
@@ -0,0 +1,387 @@
|
||||
#include "gridfire/engine/views/engine_defined.h"
|
||||
|
||||
#include <ranges>
|
||||
|
||||
#include "quill/LogMacros.h"
|
||||
|
||||
#include <string>
|
||||
#include <vector>
|
||||
#include <unordered_set>
|
||||
#include <stdexcept>
|
||||
#include <unordered_map>
|
||||
#include <utility>
|
||||
|
||||
namespace gridfire {
|
||||
using fourdst::atomic::Species;
|
||||
|
||||
DefinedEngineView::DefinedEngineView(const std::vector<std::string>& peNames, DynamicEngine& baseEngine) :
|
||||
m_baseEngine(baseEngine) {
|
||||
collect(peNames);
|
||||
}
|
||||
|
||||
const DynamicEngine & DefinedEngineView::getBaseEngine() const {
|
||||
return m_baseEngine;
|
||||
}
|
||||
|
||||
const std::vector<Species> & DefinedEngineView::getNetworkSpecies() const {
|
||||
return m_activeSpecies;
|
||||
}
|
||||
|
||||
std::expected<StepDerivatives<double>, expectations::StaleEngineError> DefinedEngineView::calculateRHSAndEnergy(
|
||||
const std::vector<double> &Y_defined,
|
||||
const double T9,
|
||||
const double rho
|
||||
) const {
|
||||
validateNetworkState();
|
||||
|
||||
const auto Y_full = mapViewToFull(Y_defined);
|
||||
const auto result = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho);
|
||||
|
||||
if (!result) {
|
||||
return std::unexpected{result.error()};
|
||||
}
|
||||
|
||||
const auto [dydt, nuclearEnergyGenerationRate] = result.value();
|
||||
StepDerivatives<double> definedResults;
|
||||
definedResults.nuclearEnergyGenerationRate = nuclearEnergyGenerationRate;
|
||||
definedResults.dydt = mapFullToView(dydt);
|
||||
return definedResults;
|
||||
}
|
||||
|
||||
void DefinedEngineView::generateJacobianMatrix(
|
||||
const std::vector<double> &Y_dynamic,
|
||||
const double T9,
|
||||
const double rho
|
||||
) const {
|
||||
validateNetworkState();
|
||||
|
||||
const auto Y_full = mapViewToFull(Y_dynamic);
|
||||
m_baseEngine.generateJacobianMatrix(Y_full, T9, rho);
|
||||
}
|
||||
|
||||
double DefinedEngineView::getJacobianMatrixEntry(
|
||||
const int i_defined,
|
||||
const int j_defined
|
||||
) const {
|
||||
validateNetworkState();
|
||||
|
||||
const size_t i_full = mapViewToFullSpeciesIndex(i_defined);
|
||||
const size_t j_full = mapViewToFullSpeciesIndex(j_defined);
|
||||
|
||||
return m_baseEngine.getJacobianMatrixEntry(i_full, j_full);
|
||||
}
|
||||
|
||||
void DefinedEngineView::generateStoichiometryMatrix() {
|
||||
validateNetworkState();
|
||||
|
||||
m_baseEngine.generateStoichiometryMatrix();
|
||||
}
|
||||
|
||||
int DefinedEngineView::getStoichiometryMatrixEntry(
|
||||
const int speciesIndex_defined,
|
||||
const int reactionIndex_defined
|
||||
) const {
|
||||
validateNetworkState();
|
||||
|
||||
const size_t i_full = mapViewToFullSpeciesIndex(speciesIndex_defined);
|
||||
const size_t j_full = mapViewToFullReactionIndex(reactionIndex_defined);
|
||||
return m_baseEngine.getStoichiometryMatrixEntry(i_full, j_full);
|
||||
}
|
||||
|
||||
double DefinedEngineView::calculateMolarReactionFlow(
|
||||
const reaction::Reaction &reaction,
|
||||
const std::vector<double> &Y_defined,
|
||||
const double T9,
|
||||
const double rho
|
||||
) const {
|
||||
validateNetworkState();
|
||||
|
||||
if (!m_activeReactions.contains(reaction)) {
|
||||
LOG_ERROR(m_logger, "Reaction '{}' is not part of the active reactions in the DefinedEngineView.", reaction.id());
|
||||
m_logger -> flush_log();
|
||||
throw std::runtime_error("Reaction not found in active reactions: " + std::string(reaction.id()));
|
||||
}
|
||||
const auto Y_full = mapViewToFull(Y_defined);
|
||||
return m_baseEngine.calculateMolarReactionFlow(reaction, Y_full, T9, rho);
|
||||
}
|
||||
|
||||
const reaction::LogicalReactionSet & DefinedEngineView::getNetworkReactions() const {
|
||||
validateNetworkState();
|
||||
|
||||
return m_activeReactions;
|
||||
}
|
||||
|
||||
void DefinedEngineView::setNetworkReactions(const reaction::LogicalReactionSet &reactions) {
|
||||
std::vector<std::string> peNames;
|
||||
for (const auto& reaction : reactions) {
|
||||
peNames.push_back(std::string(reaction.id()));
|
||||
}
|
||||
collect(peNames);
|
||||
}
|
||||
|
||||
std::expected<std::unordered_map<Species, double>, expectations::StaleEngineError> DefinedEngineView::getSpeciesTimescales(
|
||||
const std::vector<double> &Y_defined,
|
||||
const double T9,
|
||||
const double rho
|
||||
) const {
|
||||
validateNetworkState();
|
||||
|
||||
const auto Y_full = mapViewToFull(Y_defined);
|
||||
const auto result = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho);
|
||||
if (!result) {
|
||||
return std::unexpected{result.error()};
|
||||
}
|
||||
const auto& fullTimescales = result.value();
|
||||
|
||||
std::unordered_map<Species, double> definedTimescales;
|
||||
for (const auto& active_species : m_activeSpecies) {
|
||||
if (fullTimescales.contains(active_species)) {
|
||||
definedTimescales[active_species] = fullTimescales.at(active_species);
|
||||
}
|
||||
}
|
||||
return definedTimescales;
|
||||
}
|
||||
|
||||
std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError>
|
||||
DefinedEngineView::getSpeciesDestructionTimescales(
|
||||
const std::vector<double> &Y_defined,
|
||||
const double T9,
|
||||
const double rho
|
||||
) const {
|
||||
validateNetworkState();
|
||||
|
||||
const auto Y_full = mapViewToFull(Y_defined);
|
||||
const auto result = m_baseEngine.getSpeciesDestructionTimescales(Y_full, T9, rho);
|
||||
|
||||
if (!result) {
|
||||
return std::unexpected{result.error()};
|
||||
}
|
||||
|
||||
const auto& destructionTimescales = result.value();
|
||||
|
||||
std::unordered_map<Species, double> definedTimescales;
|
||||
for (const auto& active_species : m_activeSpecies) {
|
||||
if (destructionTimescales.contains(active_species)) {
|
||||
definedTimescales[active_species] = destructionTimescales.at(active_species);
|
||||
}
|
||||
}
|
||||
return definedTimescales;
|
||||
}
|
||||
|
||||
fourdst::composition::Composition DefinedEngineView::update(const NetIn &netIn) {
|
||||
return m_baseEngine.update(netIn);
|
||||
}
|
||||
|
||||
bool DefinedEngineView::isStale(const NetIn &netIn) {
|
||||
return m_baseEngine.isStale(netIn);
|
||||
}
|
||||
|
||||
|
||||
void DefinedEngineView::setScreeningModel(const screening::ScreeningType model) {
|
||||
m_baseEngine.setScreeningModel(model);
|
||||
}
|
||||
|
||||
screening::ScreeningType DefinedEngineView::getScreeningModel() const {
|
||||
return m_baseEngine.getScreeningModel();
|
||||
}
|
||||
|
||||
int DefinedEngineView::getSpeciesIndex(const Species &species) const {
|
||||
validateNetworkState();
|
||||
|
||||
const auto it = std::ranges::find(m_activeSpecies, species);
|
||||
if (it != m_activeSpecies.end()) {
|
||||
return static_cast<int>(std::distance(m_activeSpecies.begin(), it));
|
||||
} else {
|
||||
LOG_ERROR(m_logger, "Species '{}' not found in active species list.", species.name());
|
||||
m_logger->flush_log();
|
||||
throw std::runtime_error("Species not found in active species list: " + std::string(species.name()));
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<double> DefinedEngineView::mapNetInToMolarAbundanceVector(const NetIn &netIn) const {
|
||||
std::vector<double> Y(m_activeSpecies.size(), 0.0); // Initialize with zeros
|
||||
for (const auto& [symbol, entry] : netIn.composition) {
|
||||
auto it = std::ranges::find(m_activeSpecies, entry.isotope());
|
||||
if (it != m_activeSpecies.end()) {
|
||||
Y[getSpeciesIndex(entry.isotope())] = netIn.composition.getMolarAbundance(symbol); // Map species to their molar abundance
|
||||
}
|
||||
}
|
||||
return Y; // Return the vector of molar abundances
|
||||
}
|
||||
|
||||
PrimingReport DefinedEngineView::primeEngine(const NetIn &netIn) {
|
||||
return m_baseEngine.primeEngine(netIn);
|
||||
}
|
||||
|
||||
std::vector<size_t> DefinedEngineView::constructSpeciesIndexMap() const {
|
||||
LOG_TRACE_L3(m_logger, "Constructing species index map for DefinedEngineView...");
|
||||
std::unordered_map<Species, size_t> fullSpeciesReverseMap;
|
||||
const auto& fullSpeciesList = m_baseEngine.getNetworkSpecies();
|
||||
|
||||
fullSpeciesReverseMap.reserve(fullSpeciesList.size());
|
||||
|
||||
for (size_t i = 0; i < fullSpeciesList.size(); ++i) {
|
||||
fullSpeciesReverseMap[fullSpeciesList[i]] = i;
|
||||
}
|
||||
|
||||
std::vector<size_t> speciesIndexMap;
|
||||
speciesIndexMap.reserve(m_activeSpecies.size());
|
||||
|
||||
for (const auto& active_species : m_activeSpecies) {
|
||||
auto it = fullSpeciesReverseMap.find(active_species);
|
||||
if (it != fullSpeciesReverseMap.end()) {
|
||||
speciesIndexMap.push_back(it->second);
|
||||
} else {
|
||||
LOG_ERROR(m_logger, "Species '{}' not found in full species map.", active_species.name());
|
||||
m_logger -> flush_log();
|
||||
throw std::runtime_error("Species not found in full species map: " + std::string(active_species.name()));
|
||||
}
|
||||
}
|
||||
LOG_TRACE_L3(m_logger, "Species index map constructed with {} entries.", speciesIndexMap.size());
|
||||
return speciesIndexMap;
|
||||
|
||||
}
|
||||
|
||||
std::vector<size_t> DefinedEngineView::constructReactionIndexMap() const {
|
||||
LOG_TRACE_L3(m_logger, "Constructing reaction index map for DefinedEngineView...");
|
||||
|
||||
// --- Step 1: Create a reverse map using the reaction's unique ID as the key. ---
|
||||
std::unordered_map<std::string_view, size_t> fullReactionReverseMap;
|
||||
const auto& fullReactionSet = m_baseEngine.getNetworkReactions();
|
||||
fullReactionReverseMap.reserve(fullReactionSet.size());
|
||||
|
||||
for (size_t i_full = 0; i_full < fullReactionSet.size(); ++i_full) {
|
||||
fullReactionReverseMap[fullReactionSet[i_full].id()] = i_full;
|
||||
}
|
||||
|
||||
// --- Step 2: Build the final index map using the active reaction set. ---
|
||||
std::vector<size_t> reactionIndexMap;
|
||||
reactionIndexMap.reserve(m_activeReactions.size());
|
||||
|
||||
for (const auto& active_reaction_ptr : m_activeReactions) {
|
||||
auto it = fullReactionReverseMap.find(active_reaction_ptr.id());
|
||||
|
||||
if (it != fullReactionReverseMap.end()) {
|
||||
reactionIndexMap.push_back(it->second);
|
||||
} else {
|
||||
LOG_ERROR(m_logger, "Active reaction '{}' not found in base engine during reaction index map construction.", active_reaction_ptr.id());
|
||||
m_logger->flush_log();
|
||||
throw std::runtime_error("Mismatch between active reactions and base engine.");
|
||||
}
|
||||
}
|
||||
|
||||
LOG_TRACE_L3(m_logger, "Reaction index map constructed with {} entries.", reactionIndexMap.size());
|
||||
return reactionIndexMap;
|
||||
}
|
||||
|
||||
std::vector<double> DefinedEngineView::mapViewToFull(const std::vector<double>& culled) const {
|
||||
std::vector<double> full(m_baseEngine.getNetworkSpecies().size(), 0.0);
|
||||
for (size_t i_culled = 0; i_culled < culled.size(); ++i_culled) {
|
||||
const size_t i_full = m_speciesIndexMap[i_culled];
|
||||
full[i_full] += culled[i_culled];
|
||||
}
|
||||
return full;
|
||||
}
|
||||
|
||||
std::vector<double> DefinedEngineView::mapFullToView(const std::vector<double>& full) const {
|
||||
std::vector<double> culled(m_activeSpecies.size(), 0.0);
|
||||
for (size_t i_culled = 0; i_culled < m_activeSpecies.size(); ++i_culled) {
|
||||
const size_t i_full = m_speciesIndexMap[i_culled];
|
||||
culled[i_culled] = full[i_full];
|
||||
}
|
||||
return culled;
|
||||
}
|
||||
|
||||
size_t DefinedEngineView::mapViewToFullSpeciesIndex(size_t culledSpeciesIndex) const {
|
||||
if (culledSpeciesIndex < 0 || culledSpeciesIndex >= static_cast<int>(m_speciesIndexMap.size())) {
|
||||
LOG_ERROR(m_logger, "Defined index {} is out of bounds for species index map of size {}.", culledSpeciesIndex, m_speciesIndexMap.size());
|
||||
m_logger->flush_log();
|
||||
throw std::out_of_range("Defined index " + std::to_string(culledSpeciesIndex) + " is out of bounds for species index map of size " + std::to_string(m_speciesIndexMap.size()) + ".");
|
||||
}
|
||||
return m_speciesIndexMap[culledSpeciesIndex];
|
||||
}
|
||||
|
||||
size_t DefinedEngineView::mapViewToFullReactionIndex(size_t culledReactionIndex) const {
|
||||
if (culledReactionIndex < 0 || culledReactionIndex >= static_cast<int>(m_reactionIndexMap.size())) {
|
||||
LOG_ERROR(m_logger, "Defined index {} is out of bounds for reaction index map of size {}.", culledReactionIndex, m_reactionIndexMap.size());
|
||||
m_logger->flush_log();
|
||||
throw std::out_of_range("Defined index " + std::to_string(culledReactionIndex) + " is out of bounds for reaction index map of size " + std::to_string(m_reactionIndexMap.size()) + ".");
|
||||
}
|
||||
return m_reactionIndexMap[culledReactionIndex];
|
||||
}
|
||||
|
||||
void DefinedEngineView::validateNetworkState() const {
|
||||
if (m_isStale) {
|
||||
LOG_ERROR(m_logger, "DefinedEngineView is stale. Please call update() with a valid NetIn object.");
|
||||
m_logger->flush_log();
|
||||
throw std::runtime_error("DefinedEngineView is stale. Please call update() with a valid NetIn object.");
|
||||
}
|
||||
}
|
||||
|
||||
void DefinedEngineView::collect(const std::vector<std::string> &peNames) {
|
||||
std::unordered_set<Species> seenSpecies;
|
||||
|
||||
const auto& fullNetworkReactionSet = m_baseEngine.getNetworkReactions();
|
||||
for (const auto& peName : peNames) {
|
||||
if (!fullNetworkReactionSet.contains(peName)) {
|
||||
LOG_ERROR(m_logger, "Reaction with name '{}' not found in the base engine's network reactions. Aborting...", peName);
|
||||
m_logger->flush_log();
|
||||
throw std::runtime_error("Reaction with name '" + std::string(peName) + "' not found in the base engine's network reactions.");
|
||||
}
|
||||
auto reaction = fullNetworkReactionSet[peName];
|
||||
for (const auto& reactant : reaction.reactants()) {
|
||||
if (!seenSpecies.contains(reactant)) {
|
||||
seenSpecies.insert(reactant);
|
||||
m_activeSpecies.push_back(reactant);
|
||||
}
|
||||
}
|
||||
for (const auto& product : reaction.products()) {
|
||||
if (!seenSpecies.contains(product)) {
|
||||
seenSpecies.insert(product);
|
||||
m_activeSpecies.push_back(product);
|
||||
}
|
||||
}
|
||||
m_activeReactions.add_reaction(reaction);
|
||||
}
|
||||
LOG_TRACE_L3(m_logger, "DefinedEngineView built with {} active species and {} active reactions.", m_activeSpecies.size(), m_activeReactions.size());
|
||||
LOG_TRACE_L3(m_logger, "Active species: {}", [this]() -> std::string {
|
||||
std::string result;
|
||||
for (const auto& species : m_activeSpecies) {
|
||||
result += std::string(species.name()) + ", ";
|
||||
}
|
||||
if (!result.empty()) {
|
||||
result.pop_back(); // Remove last space
|
||||
result.pop_back(); // Remove last comma
|
||||
}
|
||||
return result;
|
||||
}());
|
||||
LOG_TRACE_L3(m_logger, "Active reactions: {}", [this]() -> std::string {
|
||||
std::string result;
|
||||
for (const auto& reaction : m_activeReactions) {
|
||||
result += std::string(reaction.id()) + ", ";
|
||||
}
|
||||
if (!result.empty()) {
|
||||
result.pop_back(); // Remove last space
|
||||
result.pop_back(); // Remove last comma
|
||||
}
|
||||
return result;
|
||||
}());
|
||||
m_speciesIndexMap = constructSpeciesIndexMap();
|
||||
m_reactionIndexMap = constructReactionIndexMap();
|
||||
m_isStale = false;
|
||||
}
|
||||
|
||||
|
||||
////////////////////////////////////////////
|
||||
/// FileDefinedEngineView Implementation ///
|
||||
/////////////////////////////////////////////
|
||||
|
||||
FileDefinedEngineView::FileDefinedEngineView(
|
||||
DynamicEngine &baseEngine,
|
||||
const std::string &fileName,
|
||||
const io::NetworkFileParser &parser
|
||||
):
|
||||
DefinedEngineView(parser.parse(fileName), baseEngine),
|
||||
m_fileName(fileName),
|
||||
m_parser(parser) {}
|
||||
}
|
||||
1626
src/lib/engine/views/engine_multiscale.cpp
Normal file
1626
src/lib/engine/views/engine_multiscale.cpp
Normal file
File diff suppressed because it is too large
Load Diff
73
src/lib/engine/views/engine_priming.cpp
Normal file
73
src/lib/engine/views/engine_priming.cpp
Normal file
@@ -0,0 +1,73 @@
|
||||
#include "gridfire/engine/views/engine_priming.h"
|
||||
#include "gridfire/solver/solver.h"
|
||||
|
||||
#include "fourdst/composition/species.h"
|
||||
#include "fourdst/logging/logging.h"
|
||||
|
||||
#include "quill/LogMacros.h"
|
||||
#include "quill/Logger.h"
|
||||
|
||||
#include <vector>
|
||||
#include <string>
|
||||
#include <unordered_set>
|
||||
#include <stdexcept>
|
||||
#include <unordered_map>
|
||||
#include <utility>
|
||||
#include <ranges>
|
||||
#include <cmath>
|
||||
|
||||
|
||||
namespace gridfire {
|
||||
using fourdst::atomic::species;
|
||||
|
||||
NetworkPrimingEngineView::NetworkPrimingEngineView(
|
||||
const std::string &primingSymbol,
|
||||
DynamicEngine &baseEngine
|
||||
) :
|
||||
DefinedEngineView(
|
||||
constructPrimingReactionSet(
|
||||
species.at(primingSymbol),
|
||||
baseEngine
|
||||
),
|
||||
baseEngine
|
||||
),
|
||||
m_primingSpecies(species.at(primingSymbol)) {}
|
||||
|
||||
NetworkPrimingEngineView::NetworkPrimingEngineView(
|
||||
const fourdst::atomic::Species &primingSpecies,
|
||||
DynamicEngine &baseEngine
|
||||
) :
|
||||
DefinedEngineView(
|
||||
constructPrimingReactionSet(
|
||||
primingSpecies,
|
||||
baseEngine
|
||||
),
|
||||
baseEngine
|
||||
),
|
||||
m_primingSpecies(primingSpecies) {
|
||||
}
|
||||
|
||||
|
||||
std::vector<std::string> NetworkPrimingEngineView::constructPrimingReactionSet(
|
||||
const fourdst::atomic::Species &primingSpecies,
|
||||
const DynamicEngine &baseEngine
|
||||
) const {
|
||||
std::unordered_set<std::string> primeReactions;
|
||||
for (const auto &reaction : baseEngine.getNetworkReactions()) {
|
||||
if (reaction.contains(primingSpecies)) {
|
||||
primeReactions.insert(std::string(reaction.peName()));
|
||||
}
|
||||
}
|
||||
if (primeReactions.empty()) {
|
||||
LOG_ERROR(m_logger, "No priming reactions found for species '{}'.", primingSpecies.name());
|
||||
m_logger->flush_log();
|
||||
throw std::runtime_error("No priming reactions found for species '" + std::string(primingSpecies.name()) + "'.");
|
||||
}
|
||||
std::vector<std::string> primingReactionSet(primeReactions.begin(), primeReactions.end());
|
||||
// LOG_INFO(m_logger, "Constructed priming reaction set with {} reactions for species '{}'.", primingReactionSet.size(), primingSpecies.name());
|
||||
return primingReactionSet;
|
||||
}
|
||||
|
||||
|
||||
|
||||
}
|
||||
77
src/lib/io/network_file.cpp
Normal file
77
src/lib/io/network_file.cpp
Normal file
@@ -0,0 +1,77 @@
|
||||
#include "gridfire/io/network_file.h"
|
||||
|
||||
#include <string>
|
||||
#include <vector>
|
||||
#include <algorithm>
|
||||
#include <fstream>
|
||||
#include <stdexcept>
|
||||
|
||||
#include "quill/LogMacros.h"
|
||||
|
||||
namespace gridfire::io {
|
||||
namespace {
|
||||
inline void ltrim(std::string &s) {
|
||||
s.erase(
|
||||
s.begin(),
|
||||
std::ranges::find_if(s,
|
||||
[](const unsigned char ch) {
|
||||
return !std::isspace(ch);
|
||||
})
|
||||
);
|
||||
}
|
||||
|
||||
inline void rtrim(std::string &s) {
|
||||
s.erase(
|
||||
std::find_if(
|
||||
s.rbegin(),
|
||||
s.rend(),
|
||||
[](const unsigned char ch) {
|
||||
return !std::isspace(ch);
|
||||
}).base(),
|
||||
s.end()
|
||||
);
|
||||
}
|
||||
|
||||
inline void trim(std::string &s) {
|
||||
ltrim(s);
|
||||
rtrim(s);
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
SimpleReactionListFileParser::SimpleReactionListFileParser() {}
|
||||
|
||||
ParsedNetworkData SimpleReactionListFileParser::parse(const std::string& filename) const {
|
||||
LOG_TRACE_L1(m_logger, "Parsing simple reaction list file: {}", filename);
|
||||
|
||||
std::ifstream file(filename);
|
||||
if (!file.is_open()) {
|
||||
LOG_ERROR(m_logger, "Failed to open file: {}", filename);
|
||||
m_logger -> flush_log();
|
||||
throw std::runtime_error("Could not open file: " + filename);
|
||||
}
|
||||
|
||||
ParsedNetworkData parsed;
|
||||
std::string line;
|
||||
int line_number = 0;
|
||||
while (std::getline(file, line)) {
|
||||
line_number++;
|
||||
LOG_TRACE_L3(m_logger, "Parsing reaction list file {}, line {}: {}", filename, line_number, line);
|
||||
|
||||
const size_t comment_pos = line.find('#');
|
||||
if (comment_pos != std::string::npos) {
|
||||
line = line.substr(0, comment_pos);
|
||||
}
|
||||
|
||||
trim(line);
|
||||
|
||||
if (line.empty()) {
|
||||
continue; // Skip empty lines
|
||||
}
|
||||
parsed.push_back(line);
|
||||
}
|
||||
LOG_TRACE_L1(m_logger, "Parsed {} reactions from file: {}", parsed.size(), filename);
|
||||
return parsed;
|
||||
}
|
||||
|
||||
}
|
||||
80
src/lib/network.cpp
Normal file
80
src/lib/network.cpp
Normal file
@@ -0,0 +1,80 @@
|
||||
/* ***********************************************************************
|
||||
//
|
||||
// Copyright (C) 2025 -- The 4D-STAR Collaboration
|
||||
// File Authors: Aaron Dotter, Emily Boudreaux
|
||||
// Last Modified: March 21, 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
|
||||
//
|
||||
// *********************************************************************** */
|
||||
#include "gridfire/network.h"
|
||||
#include "gridfire/reaction/reaclib.h"
|
||||
#include "gridfire/reaction/reaction.h"
|
||||
|
||||
#include <ranges>
|
||||
|
||||
#include "quill/LogMacros.h"
|
||||
|
||||
namespace gridfire {
|
||||
std::vector<double> NetIn::MolarAbundance() const {
|
||||
std::vector <double> y;
|
||||
y.reserve(composition.getRegisteredSymbols().size());
|
||||
const auto [fst, snd] = composition.getComposition();
|
||||
for (const auto &name: fst | std::views::keys) {
|
||||
y.push_back(composition.getMolarAbundance(name));
|
||||
}
|
||||
return y;
|
||||
}
|
||||
|
||||
|
||||
Network::Network(const NetworkFormat format) :
|
||||
m_config(fourdst::config::Config::getInstance()),
|
||||
m_logManager(fourdst::logging::LogManager::getInstance()),
|
||||
m_logger(m_logManager.getLogger("log")),
|
||||
m_format(format),
|
||||
m_constants(fourdst::constant::Constants::getInstance()){
|
||||
if (format == NetworkFormat::UNKNOWN) {
|
||||
LOG_ERROR(m_logger, "nuclearNetwork::Network::Network() called with UNKNOWN format");
|
||||
m_logger->flush_log();
|
||||
throw std::runtime_error("nuclearNetwork::Network::Network() called with UNKNOWN format");
|
||||
}
|
||||
}
|
||||
|
||||
NetworkFormat Network::getFormat() const {
|
||||
return m_format;
|
||||
}
|
||||
|
||||
NetworkFormat Network::setFormat(const NetworkFormat format) {
|
||||
const NetworkFormat oldFormat = m_format;
|
||||
m_format = format;
|
||||
return oldFormat;
|
||||
}
|
||||
|
||||
// Trim whitespace from both ends of a string
|
||||
std::string trim_whitespace(const std::string& str) {
|
||||
auto startIt = str.begin();
|
||||
const auto endIt = str.end();
|
||||
|
||||
while (startIt != endIt && std::isspace(static_cast<unsigned char>(*startIt))) {
|
||||
++startIt;
|
||||
}
|
||||
if (startIt == endIt) {
|
||||
return "";
|
||||
}
|
||||
const auto ritr = std::find_if(str.rbegin(), std::string::const_reverse_iterator(startIt),
|
||||
[](const unsigned char ch){ return !std::isspace(ch); });
|
||||
return std::string(startIt, ritr.base());
|
||||
}
|
||||
|
||||
}
|
||||
108
src/lib/partition/composite/partition_composite.cpp
Normal file
108
src/lib/partition/composite/partition_composite.cpp
Normal file
@@ -0,0 +1,108 @@
|
||||
#include "gridfire/partition/composite/partition_composite.h"
|
||||
|
||||
#include <vector>
|
||||
#include <set>
|
||||
|
||||
#include "gridfire/partition/partition_ground.h"
|
||||
#include "gridfire/partition/partition_rauscher_thielemann.h"
|
||||
#include "quill/LogMacros.h"
|
||||
|
||||
namespace gridfire::partition {
|
||||
CompositePartitionFunction::CompositePartitionFunction(
|
||||
const std::vector<BasePartitionType>& partitionFunctions
|
||||
) {
|
||||
for (const auto& type : partitionFunctions) {
|
||||
LOG_TRACE_L2(m_logger, "Adding partition function of type: {}", basePartitionTypeToString[type]);
|
||||
m_partitionFunctions.push_back(selectPartitionFunction(type));
|
||||
}
|
||||
}
|
||||
|
||||
CompositePartitionFunction::CompositePartitionFunction(const CompositePartitionFunction &other) {
|
||||
m_partitionFunctions.reserve(other.m_partitionFunctions.size());
|
||||
for (const auto& pf : other.m_partitionFunctions) {
|
||||
m_partitionFunctions.push_back(pf->clone());
|
||||
}
|
||||
}
|
||||
|
||||
double CompositePartitionFunction::evaluate(int z, int a, double T9) const {
|
||||
LOG_TRACE_L3(m_logger, "Evaluating partition function for Z={} A={} T9={}", z, a, T9);
|
||||
for (const auto& partitionFunction : m_partitionFunctions) {
|
||||
if (partitionFunction->supports(z, a)) {
|
||||
LOG_TRACE_L3(m_logger, "Partition function of type {} supports Z={} A={}", partitionFunction->type(), z, a);
|
||||
return partitionFunction->evaluate(z, a, T9);
|
||||
} else {
|
||||
LOG_TRACE_L3(m_logger, "Partition function of type {} does not support Z={} A={}", partitionFunction->type(), z, a);
|
||||
}
|
||||
}
|
||||
LOG_ERROR(
|
||||
m_logger,
|
||||
"No partition function supports Z={} A={} T9={}. Tried: {}",
|
||||
z,
|
||||
a,
|
||||
T9,
|
||||
type()
|
||||
);
|
||||
throw std::runtime_error("No partition function supports the given Z, A, and T9 values.");
|
||||
}
|
||||
|
||||
double CompositePartitionFunction::evaluateDerivative(int z, int a, double T9) const {
|
||||
for (const auto& partitionFunction : m_partitionFunctions) {
|
||||
if (partitionFunction->supports(z, a)) {
|
||||
LOG_TRACE_L3(m_logger, "Evaluating derivative of partition function for Z={} A={} T9={}", z, a, T9);
|
||||
return partitionFunction->evaluateDerivative(z, a, T9);
|
||||
}
|
||||
}
|
||||
LOG_ERROR(
|
||||
m_logger,
|
||||
"No partition function supports Z={} A={} T9={}. Tried: {}",
|
||||
z,
|
||||
a,
|
||||
T9,
|
||||
type()
|
||||
);
|
||||
throw std::runtime_error("No partition function supports the given Z, A, and T9 values.");
|
||||
}
|
||||
|
||||
bool CompositePartitionFunction::supports(int z, int a) const {
|
||||
for (const auto& partitionFunction : m_partitionFunctions) {
|
||||
if (partitionFunction->supports(z, a)) {
|
||||
LOG_TRACE_L2(m_logger, "Partition function supports Z={} A={}", z, a);
|
||||
return true;
|
||||
}
|
||||
}
|
||||
return false;
|
||||
}
|
||||
|
||||
std::string CompositePartitionFunction::type() const {
|
||||
std::stringstream ss;
|
||||
ss << "CompositePartitionFunction(";
|
||||
int count = 0;
|
||||
for (const auto& partitionFunction : m_partitionFunctions) {
|
||||
ss << partitionFunction->type();
|
||||
if (count < m_partitionFunctions.size() - 1) {
|
||||
ss << ", ";
|
||||
}
|
||||
count++;
|
||||
}
|
||||
ss << ")";
|
||||
std::string types = ss.str();
|
||||
return types;
|
||||
}
|
||||
|
||||
std::unique_ptr<PartitionFunction> CompositePartitionFunction::selectPartitionFunction(
|
||||
const BasePartitionType type
|
||||
) const {
|
||||
switch (type) {
|
||||
case RauscherThielemann: {
|
||||
return std::make_unique<RauscherThielemannPartitionFunction>();
|
||||
}
|
||||
case GroundState: {
|
||||
return std::make_unique<GroundStatePartitionFunction>();
|
||||
}
|
||||
default: {
|
||||
LOG_ERROR(m_logger, "Unknown partition function type");
|
||||
throw std::runtime_error("Unknown partition function type");
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
50
src/lib/partition/partition_ground.cpp
Normal file
50
src/lib/partition/partition_ground.cpp
Normal file
@@ -0,0 +1,50 @@
|
||||
#include "gridfire/partition/partition_ground.h"
|
||||
|
||||
#include <ranges>
|
||||
|
||||
#include "fourdst/logging/logging.h"
|
||||
#include "fourdst/composition/atomicSpecies.h"
|
||||
#include "fourdst/composition/species.h"
|
||||
|
||||
#include "quill/LogMacros.h"
|
||||
|
||||
namespace gridfire::partition {
|
||||
GroundStatePartitionFunction::GroundStatePartitionFunction() {
|
||||
for (const auto &isotope: fourdst::atomic::species | std::views::values) {
|
||||
m_ground_state_spin[make_key(isotope.z(), isotope.a())] = isotope.spin();
|
||||
}
|
||||
}
|
||||
|
||||
double GroundStatePartitionFunction::evaluate(
|
||||
const int z,
|
||||
const int a,
|
||||
const double T9
|
||||
) const {
|
||||
LOG_TRACE_L2(m_logger, "Evaluating ground state partition function for Z={} A={} T9={}", z, a, T9);
|
||||
const int key = make_key(z, a);
|
||||
const double spin = m_ground_state_spin.at(key);
|
||||
return (2.0 * spin) + 1.0;
|
||||
}
|
||||
|
||||
double GroundStatePartitionFunction::evaluateDerivative(
|
||||
const int z,
|
||||
const int a,
|
||||
const double T9
|
||||
) const {
|
||||
LOG_TRACE_L2(m_logger, "Evaluating derivative of ground state partition function for Z={} A={} T9={}", z, a, T9);
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
bool GroundStatePartitionFunction::supports(
|
||||
const int z,
|
||||
const int a
|
||||
) const {
|
||||
return m_ground_state_spin.contains(make_key(z, a));
|
||||
}
|
||||
constexpr int GroundStatePartitionFunction::make_key(
|
||||
const int z,
|
||||
const int a
|
||||
) {
|
||||
return z * 1000 + a; // Simple key generation for Z and A
|
||||
}
|
||||
}
|
||||
151
src/lib/partition/partition_rauscher_thielemann.cpp
Normal file
151
src/lib/partition/partition_rauscher_thielemann.cpp
Normal file
@@ -0,0 +1,151 @@
|
||||
#include "gridfire/partition/partition_rauscher_thielemann.h"
|
||||
#include "gridfire/partition/rauscher_thielemann_partition_data.h"
|
||||
#include "gridfire/partition/rauscher_thielemann_partition_data_record.h"
|
||||
|
||||
#include "fourdst/logging/logging.h"
|
||||
#include "quill/LogMacros.h"
|
||||
|
||||
#include <stdexcept>
|
||||
#include <algorithm>
|
||||
#include <vector>
|
||||
#include <array>
|
||||
#include <iostream>
|
||||
|
||||
namespace gridfire::partition {
|
||||
static constexpr std::array<double, 24> RT_TEMPERATURE_GRID_T9 = {
|
||||
0.01, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5,
|
||||
2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0
|
||||
};
|
||||
|
||||
RauscherThielemannPartitionFunction::RauscherThielemannPartitionFunction() {
|
||||
constexpr size_t numRecords = rauscher_thielemann_partition_data_len / sizeof(record::RauscherThielemannPartitionDataRecord);
|
||||
m_partitionData.reserve(numRecords);
|
||||
const auto* records = reinterpret_cast<const record::RauscherThielemannPartitionDataRecord*>(rauscher_thielemann_partition_data);
|
||||
for (size_t i = 0; i < numRecords; ++i) {
|
||||
const auto&[z, a, ground_state_spin, normalized_g_values] = records[i];
|
||||
IsotopeData data;
|
||||
data.ground_state_spin = ground_state_spin;
|
||||
std::ranges::copy(normalized_g_values, data.normalized_g_values.begin());
|
||||
const int key = make_key(z, a);
|
||||
LOG_TRACE_L3_LIMIT_EVERY_N(
|
||||
100,
|
||||
m_logger,
|
||||
"(EVERY 100) Adding Rauscher-Thielemann partition data for Z={} A={} (key={})",
|
||||
z,
|
||||
a,
|
||||
key
|
||||
);
|
||||
|
||||
m_partitionData[key] = std::move(data);
|
||||
}
|
||||
}
|
||||
|
||||
double RauscherThielemannPartitionFunction::evaluate(
|
||||
const int z,
|
||||
const int a,
|
||||
const double T9
|
||||
) const {
|
||||
LOG_TRACE_L2(m_logger, "Evaluating Rauscher-Thielemann partition function for Z={} A={} T9={}", z, a, T9);
|
||||
|
||||
const auto [bound, data, upperIndex, lowerIndex] = find(z, a, T9);
|
||||
|
||||
switch (bound) {
|
||||
case FRONT: {
|
||||
LOG_TRACE_L2(m_logger, "Using FRONT bound for Z={} A={} T9={}", z, a, T9);
|
||||
return data.normalized_g_values.front() * (2.0 * data.ground_state_spin + 1.0);
|
||||
}
|
||||
case BACK: {
|
||||
LOG_TRACE_L2(m_logger, "Using BACK bound for Z={} A={} T9={}", z, a, T9);
|
||||
return data.normalized_g_values.back() * (2.0 * data.ground_state_spin + 1.0);
|
||||
}
|
||||
case MIDDLE: {
|
||||
LOG_TRACE_L2(m_logger, "Using MIDDLE bound for Z={} A={} T9={}", z, a, T9);
|
||||
}
|
||||
}
|
||||
|
||||
const auto [T9_high, G_norm_high, T9_low, G_norm_low] = get_interpolation_points(
|
||||
upperIndex,
|
||||
lowerIndex,
|
||||
data.normalized_g_values
|
||||
);
|
||||
|
||||
const double frac = (T9 - T9_low) / (T9_high - T9_low);
|
||||
const double interpolated_g_norm = G_norm_low + frac * (G_norm_high - G_norm_low);
|
||||
|
||||
return interpolated_g_norm * (2.0 * data.ground_state_spin + 1.0);
|
||||
}
|
||||
|
||||
double RauscherThielemannPartitionFunction::evaluateDerivative(
|
||||
const int z,
|
||||
const int a,
|
||||
const double T9
|
||||
) const {
|
||||
LOG_TRACE_L2(m_logger, "Evaluating derivative of Rauscher-Thielemann partition function for Z={} A={} T9={}", z, a, T9);
|
||||
const auto [bound, data, upperIndex, lowerIndex] = find(z, a, T9);
|
||||
if (bound == FRONT || bound == BACK) {
|
||||
LOG_TRACE_L2(m_logger, "Derivative is zero for Z={} A={} T9={} (bound: {})", z, a, T9, bound == FRONT ? "FRONT" : "BACK");
|
||||
return 0.0; // Derivative is zero at the boundaries
|
||||
}
|
||||
const auto [T9_high, G_norm_high, T9_low, G_norm_low] = get_interpolation_points(
|
||||
upperIndex,
|
||||
lowerIndex,
|
||||
data.normalized_g_values
|
||||
);
|
||||
const double slope_g_norm = (G_norm_high - G_norm_low) / (T9_high - T9_low);
|
||||
return slope_g_norm * (2.0 * data.ground_state_spin + 1.0);
|
||||
}
|
||||
|
||||
bool RauscherThielemannPartitionFunction::supports(
|
||||
const int z,
|
||||
const int a
|
||||
) const {
|
||||
return m_partitionData.contains(make_key(z, a));
|
||||
}
|
||||
|
||||
RauscherThielemannPartitionFunction::InterpolationPoints RauscherThielemannPartitionFunction::get_interpolation_points(
|
||||
const size_t upper_index,
|
||||
const size_t lower_index,
|
||||
const std::array<double, 24>& normalized_g_values
|
||||
) {
|
||||
const double T_high = RT_TEMPERATURE_GRID_T9[upper_index];
|
||||
const double G_norm_high = normalized_g_values[upper_index];
|
||||
const double T_low = RT_TEMPERATURE_GRID_T9[lower_index];
|
||||
const double G_norm_low = normalized_g_values[lower_index];
|
||||
return {T_high, G_norm_high, T_low, G_norm_low};
|
||||
}
|
||||
|
||||
RauscherThielemannPartitionFunction::IdentifiedIsotope RauscherThielemannPartitionFunction::find(
|
||||
const int z,
|
||||
const int a,
|
||||
const double T9
|
||||
) const {
|
||||
const auto key = make_key(z, a);
|
||||
const auto it = m_partitionData.find(key);
|
||||
|
||||
if (it == m_partitionData.end()) {
|
||||
LOG_ERROR(m_logger, "Rauscher-Thielemann partition function data for Z={} A={} not found.", z, a);
|
||||
throw std::out_of_range("Partition function data not found for Z=" + std::to_string(z) + " A=" + std::to_string(a));
|
||||
}
|
||||
|
||||
const IsotopeData& data = it->second;
|
||||
const auto upper_it = std::ranges::lower_bound(RT_TEMPERATURE_GRID_T9, T9);
|
||||
Bounds bound;
|
||||
if (upper_it == RT_TEMPERATURE_GRID_T9.begin()) {
|
||||
bound = FRONT; // T9 is below the first grid point
|
||||
} else if (upper_it == RT_TEMPERATURE_GRID_T9.end()) {
|
||||
bound = BACK; // T9 is above the last grid point
|
||||
} else {
|
||||
bound = MIDDLE; // T9 is within the grid
|
||||
}
|
||||
const size_t upper_index = std::distance(RT_TEMPERATURE_GRID_T9.begin(), upper_it);
|
||||
const size_t lower_index = upper_index - 1;
|
||||
return {bound, data, upper_index, lower_index};
|
||||
}
|
||||
|
||||
constexpr int RauscherThielemannPartitionFunction::make_key(
|
||||
const int z,
|
||||
const int a
|
||||
) {
|
||||
return z * 1000 + a; // Simple key generation for Z and A
|
||||
}
|
||||
}
|
||||
148
src/lib/reaction/reaclib.cpp
Normal file
148
src/lib/reaction/reaclib.cpp
Normal file
@@ -0,0 +1,148 @@
|
||||
#include "fourdst/composition/atomicSpecies.h"
|
||||
#include "fourdst/composition/species.h"
|
||||
|
||||
#include "gridfire/reaction/reaclib.h"
|
||||
#include "gridfire/reaction/reactions_data.h"
|
||||
#include "gridfire/network.h"
|
||||
|
||||
#include <stdexcept>
|
||||
#include <sstream>
|
||||
#include <vector>
|
||||
#include <string>
|
||||
|
||||
std::string trim_whitespace(const std::string& str) {
|
||||
auto startIt = str.begin();
|
||||
const auto endIt = str.end();
|
||||
|
||||
while (startIt != endIt && std::isspace(static_cast<unsigned char>(*startIt))) {
|
||||
++startIt;
|
||||
}
|
||||
if (startIt == endIt) {
|
||||
return "";
|
||||
}
|
||||
const auto ritr = std::find_if(str.rbegin(), std::string::const_reverse_iterator(startIt),
|
||||
[](const unsigned char ch){ return !std::isspace(ch); });
|
||||
return std::string(startIt, ritr.base());
|
||||
}
|
||||
|
||||
namespace gridfire::reaclib {
|
||||
static reaction::LogicalReactionSet* s_all_reaclib_reactions_ptr = nullptr;
|
||||
|
||||
#pragma pack(push, 1)
|
||||
struct ReactionRecord {
|
||||
int32_t chapter;
|
||||
double qValue;
|
||||
double coeffs[7];
|
||||
bool reverse;
|
||||
char label[8];
|
||||
char rpName[64];
|
||||
char reactants_str[128];
|
||||
char products_str[128];
|
||||
};
|
||||
#pragma pack(pop)
|
||||
|
||||
std::ostream& operator<<(std::ostream& os, const ReactionRecord& r) {
|
||||
os << "Chapter: " << r.chapter
|
||||
<< ", Q-value: " << r.qValue
|
||||
<< ", Coefficients: [" << r.coeffs[0] << ", " << r.coeffs[1] << ", "
|
||||
<< r.coeffs[2] << ", " << r.coeffs[3] << ", " << r.coeffs[4] << ", "
|
||||
<< r.coeffs[5] << ", " << r.coeffs[6] << "]"
|
||||
<< ", Reverse: " << (r.reverse ? "true" : "false")
|
||||
<< ", Label: '" << std::string(r.label, strnlen(r.label, sizeof(r.label))) << "'"
|
||||
<< ", RP Name: '" << std::string(r.rpName, strnlen(r.rpName, sizeof(r.rpName))) << "'"
|
||||
<< ", Reactants: '" << std::string(r.reactants_str, strnlen(r.reactants_str, sizeof(r.reactants_str))) << "'"
|
||||
<< ", Products: '" << std::string(r.products_str, strnlen(r.products_str, sizeof(r.products_str))) << "'";
|
||||
return os;
|
||||
}
|
||||
|
||||
static std::vector<fourdst::atomic::Species> parseSpeciesString(const std::string_view str) {
|
||||
std::vector<fourdst::atomic::Species> result;
|
||||
std::stringstream ss{std::string(str)};
|
||||
std::string name;
|
||||
|
||||
while (ss >> name) {
|
||||
// Trim whitespace that might be left over from the fixed-width char arrays
|
||||
const auto trimmed_name = trim_whitespace(name);
|
||||
if (trimmed_name.empty()) continue;
|
||||
|
||||
auto it = fourdst::atomic::species.find(trimmed_name);
|
||||
if (it != fourdst::atomic::species.end()) {
|
||||
result.push_back(it->second);
|
||||
} else {
|
||||
// If a species is not found, it's a critical data error.
|
||||
throw std::runtime_error("Unknown species in reaction data: " + std::string(trimmed_name));
|
||||
}
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
static void initializeAllReaclibReactions() {
|
||||
if (s_initialized) {
|
||||
return;
|
||||
}
|
||||
|
||||
// Cast the raw byte data to our structured record format.
|
||||
const auto* records = reinterpret_cast<const ReactionRecord*>(raw_reactions_data);
|
||||
constexpr size_t num_reactions = raw_reactions_data_len / sizeof(ReactionRecord);
|
||||
|
||||
std::vector<reaction::Reaction> reaction_list;
|
||||
reaction_list.reserve(num_reactions);
|
||||
|
||||
for (size_t i = 0; i < num_reactions; ++i) {
|
||||
const auto&[chapter, qValue, coeffs, reverse, label, rpName, reactants_str, products_str] = records[i];
|
||||
|
||||
// The char arrays from the binary are not guaranteed to be null-terminated
|
||||
// if the string fills the entire buffer. We create null-terminated string_views.
|
||||
const std::string_view label_sv(label, strnlen(label, sizeof(label)));
|
||||
const std::string_view rpName_sv(rpName, strnlen(rpName, sizeof(rpName)));
|
||||
const std::string_view reactants_sv(reactants_str, strnlen(reactants_str, sizeof(reactants_str)));
|
||||
const std::string_view products_sv(products_str, strnlen(products_str, sizeof(products_str)));
|
||||
|
||||
auto reactants = parseSpeciesString(reactants_sv);
|
||||
auto products = parseSpeciesString(products_sv);
|
||||
|
||||
const reaction::RateCoefficientSet rate_coeffs = {
|
||||
coeffs[0], coeffs[1], coeffs[2],
|
||||
coeffs[3], coeffs[4], coeffs[5],
|
||||
coeffs[6]
|
||||
};
|
||||
|
||||
// Construct the Reaction object. We use rpName for both the unique ID and the human-readable name.
|
||||
reaction_list.emplace_back(
|
||||
rpName_sv,
|
||||
rpName_sv,
|
||||
chapter,
|
||||
reactants,
|
||||
products,
|
||||
qValue,
|
||||
label_sv,
|
||||
rate_coeffs,
|
||||
reverse
|
||||
);
|
||||
}
|
||||
|
||||
// The ReactionSet takes the vector of all individual reactions.
|
||||
const reaction::ReactionSet reaction_set(std::move(reaction_list));
|
||||
|
||||
// The LogicalReactionSet groups reactions by their peName, which is what we want.
|
||||
s_all_reaclib_reactions_ptr = new reaction::LogicalReactionSet(
|
||||
reaction::packReactionSetToLogicalReactionSet(reaction_set)
|
||||
);
|
||||
|
||||
s_initialized = true;
|
||||
}
|
||||
|
||||
|
||||
// --- Public Interface Implementation ---
|
||||
|
||||
const reaction::LogicalReactionSet& get_all_reactions() {
|
||||
// This ensures that the initialization happens only on the first call.
|
||||
if (!s_initialized) {
|
||||
initializeAllReaclibReactions();
|
||||
}
|
||||
if (s_all_reaclib_reactions_ptr == nullptr) {
|
||||
throw std::runtime_error("Reaclib reactions have not been initialized.");
|
||||
}
|
||||
return *s_all_reaclib_reactions_ptr;
|
||||
}
|
||||
} // namespace gridfire::reaclib
|
||||
312
src/lib/reaction/reaction.cpp
Normal file
312
src/lib/reaction/reaction.cpp
Normal file
@@ -0,0 +1,312 @@
|
||||
#include "gridfire/reaction/reaction.h"
|
||||
|
||||
#include<string_view>
|
||||
#include<string>
|
||||
#include<vector>
|
||||
#include<unordered_set>
|
||||
#include<algorithm>
|
||||
#include <ranges>
|
||||
|
||||
#include "quill/LogMacros.h"
|
||||
|
||||
#include "fourdst/composition/atomicSpecies.h"
|
||||
|
||||
#include "xxhash64.h"
|
||||
|
||||
namespace gridfire::reaction {
|
||||
using namespace fourdst::atomic;
|
||||
|
||||
Reaction::Reaction(
|
||||
const std::string_view id,
|
||||
const std::string_view peName,
|
||||
const int chapter,
|
||||
const std::vector<Species>& reactants,
|
||||
const std::vector<Species>& products,
|
||||
const double qValue,
|
||||
const std::string_view label,
|
||||
const RateCoefficientSet& sets,
|
||||
const bool reverse) :
|
||||
m_id(id),
|
||||
m_peName(peName),
|
||||
m_chapter(chapter),
|
||||
m_qValue(qValue),
|
||||
m_reactants(reactants),
|
||||
m_products(products),
|
||||
m_sourceLabel(label),
|
||||
m_rateCoefficients(sets),
|
||||
m_reverse(reverse) {}
|
||||
|
||||
double Reaction::calculate_rate(const double T9) const {
|
||||
return calculate_rate<double>(T9);
|
||||
}
|
||||
|
||||
CppAD::AD<double> Reaction::calculate_rate(const CppAD::AD<double> T9) const {
|
||||
return calculate_rate<CppAD::AD<double>>(T9);
|
||||
}
|
||||
|
||||
double Reaction::calculate_forward_rate_log_derivative(const double T9) const {
|
||||
constexpr double r_p13 = 1.0 / 3.0;
|
||||
constexpr double r_p53 = 5.0 / 3.0;
|
||||
constexpr double r_p23 = 2.0 / 3.0;
|
||||
constexpr double r_p43 = 4.0 / 3.0;
|
||||
|
||||
const double T9_m1 = 1.0 / T9;
|
||||
const double T9_m23 = std::pow(T9, -r_p23);
|
||||
const double T9_m43 = std::pow(T9, -r_p43);
|
||||
|
||||
const double d_log_k_fwd_dT9 =
|
||||
-m_rateCoefficients.a1 * T9_m1 * T9_m1
|
||||
- r_p13 * m_rateCoefficients.a2 * T9_m43
|
||||
+ r_p13 * m_rateCoefficients.a3 * T9_m23
|
||||
+ m_rateCoefficients.a4
|
||||
+ r_p53 * m_rateCoefficients.a5 * std::pow(T9, r_p23)
|
||||
+ m_rateCoefficients.a6 * T9_m1;
|
||||
|
||||
return d_log_k_fwd_dT9; // Return the derivative of the log rate with respect to T9
|
||||
}
|
||||
|
||||
bool Reaction::contains(const Species &species) const {
|
||||
return contains_reactant(species) || contains_product(species);
|
||||
}
|
||||
|
||||
|
||||
bool Reaction::contains_reactant(const Species& species) const {
|
||||
for (const auto& reactant : m_reactants) {
|
||||
if (reactant == species) {
|
||||
return true;
|
||||
}
|
||||
}
|
||||
return false;
|
||||
}
|
||||
|
||||
bool Reaction::contains_product(const Species& species) const {
|
||||
for (const auto& product : m_products) {
|
||||
if (product == species) {
|
||||
return true;
|
||||
}
|
||||
}
|
||||
return false;
|
||||
}
|
||||
|
||||
std::unordered_set<Species> Reaction::all_species() const {
|
||||
auto rs = reactant_species();
|
||||
auto ps = product_species();
|
||||
rs.insert(ps.begin(), ps.end());
|
||||
return rs;
|
||||
}
|
||||
|
||||
std::unordered_set<Species> Reaction::reactant_species() const {
|
||||
std::unordered_set<Species> reactantsSet;
|
||||
for (const auto& reactant : m_reactants) {
|
||||
reactantsSet.insert(reactant);
|
||||
}
|
||||
return reactantsSet;
|
||||
}
|
||||
|
||||
std::unordered_set<Species> Reaction::product_species() const {
|
||||
std::unordered_set<Species> productsSet;
|
||||
for (const auto& product : m_products) {
|
||||
productsSet.insert(product);
|
||||
}
|
||||
return productsSet;
|
||||
}
|
||||
|
||||
int Reaction::stoichiometry(const Species& species) const {
|
||||
int s = 0;
|
||||
for (const auto& reactant : m_reactants) {
|
||||
if (reactant == species) {
|
||||
s--;
|
||||
}
|
||||
}
|
||||
for (const auto& product : m_products) {
|
||||
if (product == species) {
|
||||
s++;
|
||||
}
|
||||
}
|
||||
return s;
|
||||
}
|
||||
|
||||
size_t Reaction::num_species() const {
|
||||
return all_species().size();
|
||||
}
|
||||
|
||||
std::unordered_map<Species, int> Reaction::stoichiometry() const {
|
||||
std::unordered_map<Species, int> stoichiometryMap;
|
||||
for (const auto& reactant : m_reactants) {
|
||||
stoichiometryMap[reactant]--;
|
||||
}
|
||||
for (const auto& product : m_products) {
|
||||
stoichiometryMap[product]++;
|
||||
}
|
||||
return stoichiometryMap;
|
||||
}
|
||||
|
||||
double Reaction::excess_energy() const {
|
||||
double reactantMass = 0.0;
|
||||
double productMass = 0.0;
|
||||
constexpr double AMU2MeV = 931.494893; // Conversion factor from atomic mass unit to MeV
|
||||
for (const auto& reactant : m_reactants) {
|
||||
reactantMass += reactant.mass();
|
||||
}
|
||||
for (const auto& product : m_products) {
|
||||
productMass += product.mass();
|
||||
}
|
||||
return (reactantMass - productMass) * AMU2MeV;
|
||||
}
|
||||
|
||||
uint64_t Reaction::hash(const uint64_t seed) const {
|
||||
return XXHash64::hash(m_id.data(), m_id.size(), seed);
|
||||
}
|
||||
|
||||
|
||||
|
||||
LogicalReaction::LogicalReaction(const std::vector<Reaction>& reactants) :
|
||||
Reaction(reactants.front().peName(),
|
||||
reactants.front().peName(),
|
||||
reactants.front().chapter(),
|
||||
reactants.front().reactants(),
|
||||
reactants.front().products(),
|
||||
reactants.front().qValue(),
|
||||
reactants.front().sourceLabel(),
|
||||
reactants.front().rateCoefficients(),
|
||||
reactants.front().is_reverse()) {
|
||||
|
||||
m_sources.reserve(reactants.size());
|
||||
m_rates.reserve(reactants.size());
|
||||
for (const auto& reaction : reactants) {
|
||||
if (std::abs(std::abs(reaction.qValue()) - std::abs(m_qValue)) > 1e-6) {
|
||||
LOG_ERROR(
|
||||
m_logger,
|
||||
"LogicalReaction constructed with reactions having different Q-values. Expected {} got {}.",
|
||||
m_qValue,
|
||||
reaction.qValue()
|
||||
);
|
||||
m_logger -> flush_log();
|
||||
throw std::runtime_error("LogicalReaction constructed with reactions having different Q-values. Expected " + std::to_string(m_qValue) + " got " + std::to_string(reaction.qValue()) + " (difference : " + std::to_string(std::abs(reaction.qValue() - m_qValue)) + ").");
|
||||
}
|
||||
m_sources.push_back(std::string(reaction.sourceLabel()));
|
||||
m_rates.push_back(reaction.rateCoefficients());
|
||||
}
|
||||
}
|
||||
|
||||
void LogicalReaction::add_reaction(const Reaction& reaction) {
|
||||
if (reaction.peName() != m_id) {
|
||||
LOG_ERROR(m_logger, "Cannot add reaction with different peName to LogicalReaction. Expected {} got {}.", m_id, reaction.peName());
|
||||
m_logger -> flush_log();
|
||||
throw std::runtime_error("Cannot add reaction with different peName to LogicalReaction. Expected " + std::string(m_id) + " got " + std::string(reaction.peName()) + ".");
|
||||
}
|
||||
for (const auto& source : m_sources) {
|
||||
if (source == reaction.sourceLabel()) {
|
||||
LOG_ERROR(m_logger, "Cannot add reaction with duplicate source label {} to LogicalReaction.", reaction.sourceLabel());
|
||||
m_logger -> flush_log();
|
||||
throw std::runtime_error("Cannot add reaction with duplicate source label " + std::string(reaction.sourceLabel()) + " to LogicalReaction.");
|
||||
}
|
||||
}
|
||||
if (std::abs(reaction.qValue() - m_qValue) > 1e-6) {
|
||||
LOG_ERROR(m_logger, "LogicalReaction constructed with reactions having different Q-values. Expected {} got {}.", m_qValue, reaction.qValue());
|
||||
m_logger -> flush_log();
|
||||
throw std::runtime_error("LogicalReaction constructed with reactions having different Q-values. Expected " + std::to_string(m_qValue) + " got " + std::to_string(reaction.qValue()) + ".");
|
||||
}
|
||||
m_sources.push_back(std::string(reaction.sourceLabel()));
|
||||
m_rates.push_back(reaction.rateCoefficients());
|
||||
}
|
||||
|
||||
double LogicalReaction::calculate_rate(const double T9) const {
|
||||
return calculate_rate<double>(T9);
|
||||
}
|
||||
|
||||
double LogicalReaction::calculate_forward_rate_log_derivative(const double T9) const {
|
||||
constexpr double r_p13 = 1.0 / 3.0;
|
||||
constexpr double r_p53 = 5.0 / 3.0;
|
||||
constexpr double r_p23 = 2.0 / 3.0;
|
||||
constexpr double r_p43 = 4.0 / 3.0;
|
||||
|
||||
double totalRate = 0.0;
|
||||
double totalRateDerivative = 0.0;
|
||||
|
||||
|
||||
const double T9_m1 = 1.0 / T9;
|
||||
const double T913 = std::pow(T9, r_p13);
|
||||
const double T953 = std::pow(T9, r_p53);
|
||||
const double logT9 = std::log(T9);
|
||||
|
||||
const double T9_m1_sq = T9_m1 * T9_m1;
|
||||
const double T9_m23 = std::pow(T9, -r_p23);
|
||||
const double T9_m43 = std::pow(T9, -r_p43);
|
||||
const double T9_p23 = std::pow(T9, r_p23);
|
||||
|
||||
|
||||
for (const auto& coeffs : m_rates) {
|
||||
const double exponent = coeffs.a0 +
|
||||
coeffs.a1 * T9_m1 +
|
||||
coeffs.a2 / T913 +
|
||||
coeffs.a3 * T913 +
|
||||
coeffs.a4 * T9 +
|
||||
coeffs.a5 * T953 +
|
||||
coeffs.a6 * logT9;
|
||||
const double individualRate = std::exp(exponent);
|
||||
|
||||
const double d_exponent_T9 =
|
||||
-coeffs.a1 * T9_m1_sq
|
||||
- r_p13 * coeffs.a2 * T9_m43
|
||||
+ r_p13 * coeffs.a3 * T9_m23
|
||||
+ coeffs.a4
|
||||
+ r_p53 * coeffs.a5 * T9_p23
|
||||
+ coeffs.a6 * T9_m1;
|
||||
|
||||
const double individualRateDerivative = individualRate * d_exponent_T9;
|
||||
totalRate += individualRate;
|
||||
totalRateDerivative += individualRateDerivative;
|
||||
}
|
||||
|
||||
if (totalRate == 0.0) {
|
||||
return 0.0; // Avoid division by zero
|
||||
}
|
||||
|
||||
return totalRateDerivative / totalRate;
|
||||
}
|
||||
|
||||
CppAD::AD<double> LogicalReaction::calculate_rate(const CppAD::AD<double> T9) const {
|
||||
return calculate_rate<CppAD::AD<double>>(T9);
|
||||
}
|
||||
|
||||
LogicalReactionSet packReactionSetToLogicalReactionSet(const ReactionSet& reactionSet) {
|
||||
std::unordered_map<std::string_view, std::vector<Reaction>> groupedReactions;
|
||||
|
||||
for (const auto& reaction: reactionSet) {
|
||||
groupedReactions[reaction.peName()].push_back(reaction);
|
||||
}
|
||||
|
||||
std::vector<LogicalReaction> reactions;
|
||||
reactions.reserve(groupedReactions.size());
|
||||
|
||||
for (const auto &reactionsGroup: groupedReactions | std::views::values) {
|
||||
LogicalReaction logicalReaction(reactionsGroup);
|
||||
reactions.push_back(logicalReaction);
|
||||
}
|
||||
return LogicalReactionSet(std::move(reactions));
|
||||
}
|
||||
}
|
||||
|
||||
namespace std {
|
||||
template<>
|
||||
struct hash<gridfire::reaction::Reaction> {
|
||||
size_t operator()(const gridfire::reaction::Reaction& r) const noexcept {
|
||||
return r.hash(0);
|
||||
}
|
||||
};
|
||||
|
||||
template<>
|
||||
struct hash<gridfire::reaction::ReactionSet> {
|
||||
size_t operator()(const gridfire::reaction::ReactionSet& s) const noexcept {
|
||||
return s.hash(0);
|
||||
}
|
||||
};
|
||||
|
||||
template<>
|
||||
struct hash<gridfire::reaction::LogicalReactionSet> {
|
||||
size_t operator()(const gridfire::reaction::LogicalReactionSet& s) const noexcept {
|
||||
return s.hash(0);
|
||||
}
|
||||
};
|
||||
} // namespace std
|
||||
31
src/lib/screening/screening_bare.cpp
Normal file
31
src/lib/screening/screening_bare.cpp
Normal file
@@ -0,0 +1,31 @@
|
||||
#include "gridfire/screening/screening_bare.h"
|
||||
|
||||
#include "fourdst/composition/atomicSpecies.h"
|
||||
|
||||
#include "cppad/cppad.hpp"
|
||||
|
||||
#include <vector>
|
||||
|
||||
|
||||
namespace gridfire::screening {
|
||||
using ADDouble = CppAD::AD<double>;
|
||||
std::vector<ADDouble> BareScreeningModel::calculateScreeningFactors(
|
||||
const reaction::LogicalReactionSet &reactions,
|
||||
const std::vector<fourdst::atomic::Species>& species,
|
||||
const std::vector<ADDouble> &Y,
|
||||
const ADDouble T9,
|
||||
const ADDouble rho
|
||||
) const {
|
||||
return calculateFactors_impl<ADDouble>(reactions, species, Y, T9, rho);
|
||||
}
|
||||
|
||||
std::vector<double> BareScreeningModel::calculateScreeningFactors(
|
||||
const reaction::LogicalReactionSet &reactions,
|
||||
const std::vector<fourdst::atomic::Species>& species,
|
||||
const std::vector<double> &Y,
|
||||
const double T9,
|
||||
const double rho
|
||||
) const {
|
||||
return calculateFactors_impl<double>(reactions, species, Y, T9, rho);
|
||||
}
|
||||
}
|
||||
19
src/lib/screening/screening_types.cpp
Normal file
19
src/lib/screening/screening_types.cpp
Normal file
@@ -0,0 +1,19 @@
|
||||
#include "gridfire/screening/screening_abstract.h"
|
||||
#include "gridfire/screening/screening_types.h"
|
||||
#include "gridfire/screening/screening_weak.h"
|
||||
#include "gridfire/screening/screening_bare.h"
|
||||
|
||||
#include <memory>
|
||||
|
||||
namespace gridfire::screening {
|
||||
std::unique_ptr<ScreeningModel> selectScreeningModel(const ScreeningType type) {
|
||||
switch (type) {
|
||||
case ScreeningType::WEAK:
|
||||
return std::make_unique<WeakScreeningModel>();
|
||||
case ScreeningType::BARE:
|
||||
return std::make_unique<BareScreeningModel>();
|
||||
default:
|
||||
return std::make_unique<BareScreeningModel>();
|
||||
}
|
||||
}
|
||||
}
|
||||
31
src/lib/screening/screening_weak.cpp
Normal file
31
src/lib/screening/screening_weak.cpp
Normal file
@@ -0,0 +1,31 @@
|
||||
#include "gridfire/screening/screening_weak.h"
|
||||
|
||||
#include "fourdst/composition/atomicSpecies.h"
|
||||
|
||||
#include "cppad/cppad.hpp"
|
||||
|
||||
#include <vector>
|
||||
|
||||
|
||||
namespace gridfire::screening {
|
||||
using ADDouble = CppAD::AD<double>;
|
||||
std::vector<ADDouble> WeakScreeningModel::calculateScreeningFactors(
|
||||
const reaction::LogicalReactionSet &reactions,
|
||||
const std::vector<fourdst::atomic::Species>& species,
|
||||
const std::vector<ADDouble> &Y,
|
||||
const ADDouble T9,
|
||||
const ADDouble rho
|
||||
) const {
|
||||
return calculateFactors_impl<ADDouble>(reactions, species, Y, T9, rho);
|
||||
}
|
||||
|
||||
std::vector<double> WeakScreeningModel::calculateScreeningFactors(
|
||||
const reaction::LogicalReactionSet &reactions,
|
||||
const std::vector<fourdst::atomic::Species>& species,
|
||||
const std::vector<double> &Y,
|
||||
const double T9,
|
||||
const double rho
|
||||
) const {
|
||||
return calculateFactors_impl<double>(reactions, species, Y, T9, rho);
|
||||
}
|
||||
}
|
||||
231
src/lib/solver/solver.cpp
Normal file
231
src/lib/solver/solver.cpp
Normal file
@@ -0,0 +1,231 @@
|
||||
#include "gridfire/solver/solver.h"
|
||||
#include "gridfire/engine/engine_graph.h"
|
||||
#include "gridfire/network.h"
|
||||
#include "gridfire/exceptions/error_engine.h"
|
||||
|
||||
#include "fourdst/composition/atomicSpecies.h"
|
||||
#include "fourdst/composition/composition.h"
|
||||
#include "fourdst/config/config.h"
|
||||
|
||||
#include "unsupported/Eigen/NonLinearOptimization"
|
||||
|
||||
#include <boost/numeric/odeint.hpp>
|
||||
|
||||
#include <vector>
|
||||
#include <string>
|
||||
#include <stdexcept>
|
||||
#include <iomanip>
|
||||
|
||||
#include "quill/LogMacros.h"
|
||||
|
||||
namespace gridfire::solver {
|
||||
NetOut DirectNetworkSolver::evaluate(const NetIn &netIn) {
|
||||
namespace ublas = boost::numeric::ublas;
|
||||
namespace odeint = boost::numeric::odeint;
|
||||
using fourdst::composition::Composition;
|
||||
|
||||
|
||||
const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9)
|
||||
|
||||
const auto absTol = m_config.get<double>("gridfire:solver:DirectNetworkSolver:absTol", 1.0e-8);
|
||||
const auto relTol = m_config.get<double>("gridfire:solver:DirectNetworkSolver:relTol", 1.0e-8);
|
||||
|
||||
Composition equilibratedComposition = m_engine.update(netIn);
|
||||
size_t numSpecies = m_engine.getNetworkSpecies().size();
|
||||
ublas::vector<double> Y(numSpecies + 1);
|
||||
|
||||
RHSManager manager(m_engine, T9, netIn.density);
|
||||
JacobianFunctor jacobianFunctor(m_engine, T9, netIn.density);
|
||||
|
||||
auto populateY = [&](const Composition& comp) {
|
||||
const size_t numSpeciesInternal = m_engine.getNetworkSpecies().size();
|
||||
Y.resize(numSpeciesInternal + 1);
|
||||
for (size_t i = 0; i < numSpeciesInternal; i++) {
|
||||
const auto& species = m_engine.getNetworkSpecies()[i];
|
||||
if (!comp.contains(species)) {
|
||||
double lim = std::numeric_limits<double>::min();
|
||||
LOG_DEBUG(m_logger, "Species '{}' not found in composition. Setting abundance to {:0.3E}.", species.name(), lim);
|
||||
Y(i) = lim; // Species not in the composition, set to zero
|
||||
} else {
|
||||
Y(i) = comp.getMolarAbundance(species);
|
||||
}
|
||||
}
|
||||
|
||||
// TODO: a good starting point to make the temperature, density, and energy self consistent would be to turn this into an accumulator
|
||||
Y(numSpeciesInternal) = 0.0; // Specific energy rate, initialized to zero
|
||||
};
|
||||
|
||||
// This is a quick debug that can be turned on. For solar code input parameters (T~1.5e7K, ρ~1.5e3 g/cm^3) this should be near 8e-17
|
||||
// std::cout << "D/H: " << equilibratedComposition.getMolarAbundance("H-2") / equilibratedComposition.getMolarAbundance("H-1") << std::endl;
|
||||
|
||||
populateY(equilibratedComposition);
|
||||
const auto stepper = odeint::make_controlled<odeint::rosenbrock4<double>>(absTol, relTol);
|
||||
|
||||
double current_time = 0.0;
|
||||
double current_initial_timestep = netIn.dt0;
|
||||
double accumulated_energy = 0.0;
|
||||
// size_t total_update_stages_triggered = 0;
|
||||
|
||||
while (current_time < netIn.tMax) {
|
||||
try {
|
||||
odeint::integrate_adaptive(
|
||||
stepper,
|
||||
std::make_pair(manager, jacobianFunctor),
|
||||
Y,
|
||||
current_time,
|
||||
netIn.tMax,
|
||||
current_initial_timestep,
|
||||
[&](const auto& state, double t) {
|
||||
current_time = t;
|
||||
manager.observe(state, t);
|
||||
}
|
||||
);
|
||||
current_time = netIn.tMax;
|
||||
} catch (const exceptions::StaleEngineTrigger &e) {
|
||||
LOG_INFO(m_logger, "Catching StaleEngineTrigger at t = {:0.3E} with T9 = {:0.3E}, rho = {:0.3E}. Triggering update stage (last stage took {} steps).", current_time, T9, netIn.density, e.totalSteps());
|
||||
exceptions::StaleEngineTrigger::state staleState = e.getState();
|
||||
accumulated_energy += e.energy(); // Add the specific energy rate to the accumulated energy
|
||||
// total_update_stages_triggered++;
|
||||
|
||||
Composition temp_comp;
|
||||
std::vector<double> mass_fractions;
|
||||
size_t num_species_at_stop = e.numSpecies();
|
||||
|
||||
if (num_species_at_stop != m_engine.getNetworkSpecies().size()) {
|
||||
throw std::runtime_error(
|
||||
"StaleEngineError state has a different number of species than the engine. This should not happen."
|
||||
);
|
||||
}
|
||||
mass_fractions.reserve(num_species_at_stop);
|
||||
|
||||
for (size_t i = 0; i < num_species_at_stop; ++i) {
|
||||
const auto& species = m_engine.getNetworkSpecies()[i];
|
||||
temp_comp.registerSpecies(species);
|
||||
mass_fractions.push_back(e.getMolarAbundance(i) * species.mass()); // Convert from molar abundance to mass fraction
|
||||
}
|
||||
temp_comp.setMassFraction(m_engine.getNetworkSpecies(), mass_fractions);
|
||||
temp_comp.finalize(true);
|
||||
|
||||
NetIn netInTemp = netIn;
|
||||
netInTemp.temperature = e.temperature();
|
||||
netInTemp.density = e.density();
|
||||
netInTemp.composition = std::move(temp_comp);
|
||||
|
||||
Composition currentComposition = m_engine.update(netInTemp);
|
||||
populateY(currentComposition);
|
||||
Y(Y.size() - 1) = e.energy(); // Set the specific energy rate from the stale state
|
||||
numSpecies = m_engine.getNetworkSpecies().size();
|
||||
|
||||
// current_initial_timestep = 0.001 * manager.m_last_step_time; // set the new timestep to the last successful timestep before repartitioning
|
||||
}
|
||||
}
|
||||
|
||||
accumulated_energy += Y(Y.size() - 1); // Add the specific energy rate to the accumulated energy
|
||||
|
||||
std::vector<double> finalMassFractions(numSpecies);
|
||||
for (size_t i = 0; i < numSpecies; ++i) {
|
||||
const double molarMass = m_engine.getNetworkSpecies()[i].mass();
|
||||
finalMassFractions[i] = Y(i) * molarMass; // Convert from molar abundance to mass fraction
|
||||
if (finalMassFractions[i] < MIN_ABUNDANCE_THRESHOLD) {
|
||||
finalMassFractions[i] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<std::string> speciesNames;
|
||||
speciesNames.reserve(numSpecies);
|
||||
for (const auto& species : m_engine.getNetworkSpecies()) {
|
||||
speciesNames.push_back(std::string(species.name()));
|
||||
}
|
||||
|
||||
Composition outputComposition(speciesNames);
|
||||
outputComposition.setMassFraction(speciesNames, finalMassFractions);
|
||||
outputComposition.finalize(true);
|
||||
|
||||
NetOut netOut;
|
||||
netOut.composition = std::move(outputComposition);
|
||||
netOut.energy = accumulated_energy; // Specific energy rate
|
||||
netOut.num_steps = manager.m_num_steps;
|
||||
|
||||
return netOut;
|
||||
}
|
||||
|
||||
void DirectNetworkSolver::RHSManager::operator()(
|
||||
const boost::numeric::ublas::vector<double> &Y,
|
||||
boost::numeric::ublas::vector<double> &dYdt,
|
||||
const double t
|
||||
) const {
|
||||
const size_t numSpecies = m_engine.getNetworkSpecies().size();
|
||||
if (t != m_cached_time || !m_cached_result.has_value() || m_cached_result.value().dydt.size() != numSpecies + 1) {
|
||||
compute_and_cache(Y, t);
|
||||
}
|
||||
const auto&[dydt, nuclearEnergyGenerationRate] = m_cached_result.value();
|
||||
dYdt.resize(numSpecies + 1);
|
||||
std::ranges::copy(dydt, dYdt.begin());
|
||||
dYdt(numSpecies) = nuclearEnergyGenerationRate; // Set the last element to the specific energy rate
|
||||
}
|
||||
|
||||
void DirectNetworkSolver::RHSManager::observe(
|
||||
const boost::numeric::ublas::vector<double> &state,
|
||||
const double t
|
||||
) const {
|
||||
double dt = t - m_last_observed_time;
|
||||
compute_and_cache(state, t);
|
||||
LOG_INFO(
|
||||
m_logger,
|
||||
"(Step {}) Observed state at t = {:0.3E} (dt = {:0.3E})",
|
||||
m_num_steps,
|
||||
t,
|
||||
dt
|
||||
);
|
||||
std::ostringstream oss;
|
||||
oss << std::scientific << std::setprecision(3);
|
||||
oss << "(Step: " << std::setw(10) << m_num_steps << ") t = " << t << " (dt = " << dt << ", eps_nuc: " << state(state.size() - 1) << " [erg])\n";
|
||||
std::cout << oss.str();
|
||||
m_last_observed_time = t;
|
||||
m_last_step_time = dt;
|
||||
|
||||
}
|
||||
|
||||
void DirectNetworkSolver::RHSManager::compute_and_cache(
|
||||
const boost::numeric::ublas::vector<double> &state,
|
||||
double t
|
||||
) const {
|
||||
std::vector<double> y_vec(state.begin(), state.end() - 1);
|
||||
std::ranges::replace_if(
|
||||
y_vec,
|
||||
[](const double yi){
|
||||
return yi < 0.0;
|
||||
},
|
||||
0.0 // Avoid negative abundances
|
||||
);
|
||||
|
||||
const auto result = m_engine.calculateRHSAndEnergy(y_vec, m_T9, m_rho);
|
||||
if (!result) {
|
||||
LOG_INFO(m_logger,
|
||||
"Triggering update stage due to stale engine detected at t = {:0.3E} with T9 = {:0.3E}, rho = {:0.3E}, y_vec (size: {})",
|
||||
t, m_T9, m_rho, y_vec.size());
|
||||
throw exceptions::StaleEngineTrigger({m_T9, m_rho, y_vec, t, m_num_steps, state(state.size() - 1)});
|
||||
}
|
||||
m_cached_result = result.value();
|
||||
m_cached_time = t;
|
||||
|
||||
m_num_steps++;
|
||||
}
|
||||
|
||||
void DirectNetworkSolver::JacobianFunctor::operator()(
|
||||
const boost::numeric::ublas::vector<double> &Y,
|
||||
boost::numeric::ublas::matrix<double> &J,
|
||||
double t,
|
||||
boost::numeric::ublas::vector<double> &dfdt
|
||||
) const {
|
||||
size_t numSpecies = m_engine.getNetworkSpecies().size();
|
||||
J.resize(numSpecies+1, numSpecies+1);
|
||||
J.clear();
|
||||
for (int i = 0; i < numSpecies; ++i) {
|
||||
for (int j = 0; j < numSpecies; ++j) {
|
||||
J(i, j) = m_engine.getJacobianMatrixEntry(i, j);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
66
src/lib/utils/logging.cpp
Normal file
66
src/lib/utils/logging.cpp
Normal file
@@ -0,0 +1,66 @@
|
||||
#include "gridfire/utils/logging.h"
|
||||
#include "gridfire/engine/engine_abstract.h"
|
||||
|
||||
#include <sstream>
|
||||
#include <iomanip>
|
||||
#include <algorithm>
|
||||
#include <ranges>
|
||||
#include <string_view>
|
||||
#include <string>
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
|
||||
std::string gridfire::utils::formatNuclearTimescaleLogString(
|
||||
const DynamicEngine& engine,
|
||||
std::vector<double> const& Y,
|
||||
const double T9,
|
||||
const double rho
|
||||
) {
|
||||
auto const& result = engine.getSpeciesTimescales(Y, T9, rho);
|
||||
if (!result) {
|
||||
std::ostringstream ss;
|
||||
ss << "Failed to get species timescales: " << result.error();
|
||||
return ss.str();
|
||||
}
|
||||
const std::unordered_map<fourdst::atomic::Species, double>& timescales = result.value();
|
||||
|
||||
// Figure out how wide the "Species" column needs to be:
|
||||
std::size_t maxNameLen = std::string_view("Species").size();
|
||||
for (const auto &key: timescales | std::views::keys) {
|
||||
std::string_view name = key.name();
|
||||
maxNameLen = std::max(maxNameLen, name.size());
|
||||
}
|
||||
|
||||
// Pick a fixed width for the timescale column:
|
||||
constexpr int timescaleWidth = 12;
|
||||
|
||||
std::ostringstream ss;
|
||||
ss << "== Timescales (s) ==\n";
|
||||
|
||||
// Header row
|
||||
ss << std::left << std::setw(static_cast<int>(maxNameLen) + 2) << "Species"
|
||||
<< std::right << std::setw(timescaleWidth) << "Timescale (s)" << "\n";
|
||||
|
||||
// Underline
|
||||
ss << std::string(static_cast<int>(maxNameLen) + 2 + timescaleWidth, '=') << "\n";
|
||||
|
||||
ss << std::scientific;
|
||||
|
||||
// Data rows
|
||||
for (auto const& [species, timescale] : timescales) {
|
||||
const std::string_view name = species.name();
|
||||
ss << std::left << std::setw(static_cast<int>(maxNameLen) + 2) << name;
|
||||
|
||||
if (std::isinf(timescale)) {
|
||||
ss << std::right << std::setw(timescaleWidth) << "inf" << "\n";
|
||||
} else {
|
||||
ss << std::right << std::setw(timescaleWidth)
|
||||
<< std::scientific << std::setprecision(3) << timescale << "\n";
|
||||
}
|
||||
}
|
||||
|
||||
// Footer underline
|
||||
ss << std::string(static_cast<int>(maxNameLen) + 2 + timescaleWidth, '=') << "\n";
|
||||
|
||||
return ss.str();
|
||||
}
|
||||
Reference in New Issue
Block a user