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

refactor(network): updated network and network::approx8 to use composition module

Browse Source 
This is a very basic wrapper implimentation currently. This is sufficient to lock the interface down so that other code can target it. However, internally there is just a "convert" function. Eventually we should rework the code itself to use the composition module more directly.
This commit is contained in:
Emily Boudreaux
2025-06-17 09:43:43 -04:00
parent 06de84592e
commit 70f13b7222
8 changed files with 366 additions and 298 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
assets/static/atomic/meson.build
View File

@@ -1,3 +1,4 @@
species_weight_dep = declare_dependency( species_weight_dep = declare_dependency(
include_directories: include_directories('include'), include_directories: include_directories('include'),
) )
message('✅ SERiF species_weight dependency declared')

4
src/network/meson.build
View File

@@ -15,7 +15,9 @@ dependencies = [
quill_dep, quill_dep,
mfem_dep, mfem_dep,
config_dep, config_dep,
probe_dep probe_dep,
species_weight_dep,
composition_dep,
] ]
# Define the libnetwork library so it can be linked against by other parts of the build system # Define the libnetwork library so it can be linked against by other parts of the build system

447
src/network/private/approx8.cpp
View File

@@ -31,7 +31,7 @@
#include "approx8.h" #include "approx8.h"
#include "network.h" #include "network.h"
/* Nuclear reaction network in cgs units based on Frank Timmes' "aprox8". /* Nuclear reaction network in cgs units based on Frank Timmes' "approx8".
At this time it does neither screening nor neutrino losses. It includes At this time it does neither screening nor neutrino losses. It includes
the following 8 isotopes: the following 8 isotopes:
@@ -78,7 +78,7 @@ namespace serif::network::approx8{
using namespace boost::numeric::odeint; using namespace boost::numeric::odeint;
//helper functions //helper functions
// a function to multilpy two arrays and then sum the resulting elements: sum(a*b) // 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_product( const vec7 &a, const vec7 &b){
if (a.size() != b.size()) { if (a.size() != b.size()) {
throw std::runtime_error("Error: array size mismatch in sum_product"); throw std::runtime_error("Error: array size mismatch in sum_product");
@@ -96,8 +96,8 @@ namespace serif::network::approx8{
// this function returns an array of the T9 terms in that order, where T9 is the temperatures in GigaKelvin // 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 get_T9_array(const double &T) {
vec7 arr; vec7 arr;
double T9=1e-9*T; const double T9=1e-9*T;
double T913=pow(T9,1./3.); const double T913=pow(T9,1./3.);
arr[0]=1; arr[0]=1;
arr[1]=1/T9; arr[1]=1/T9;
@@ -117,125 +117,125 @@ namespace serif::network::approx8{
// p + p -> d; this, like some of the other rates, this is a composite of multiple fits // p + p -> d; this, like some of the other rates, this is a composite of multiple fits
double pp_rate(const vec7 &T9) { double pp_rate(const vec7 &T9) {
vec7 a1 = {-34.78630, 0,-3.511930, 3.100860, -0.1983140, 1.262510e-2, -1.025170}; constexpr vec7 a1 = {-34.78630, 0,-3.511930, 3.100860, -0.1983140, 1.262510e-2, -1.025170};
vec7 a2 = { -4.364990e+1,-2.460640e-3,-2.750700,-4.248770e-1,1.598700e-2,-6.908750e-4,-2.076250e-1}; 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); return rate_fit(T9,a1) + rate_fit(T9,a2);
} }
// p + d -> he3 // p + d -> he3
double dp_rate(const vec7 &T9) { double dp_rate(const vec7 &T9) {
vec7 a1 = {7.528980, 0, -3.720800, 0.8717820, 0, 0,-0.6666670}; constexpr vec7 a1 = {7.528980, 0, -3.720800, 0.8717820, 0, 0,-0.6666670};
vec7 a2 = {8.935250, 0, -3.720800, 0.1986540, 0, 0, 0.3333330}; constexpr vec7 a2 = {8.935250, 0, -3.720800, 0.1986540, 0, 0, 0.3333330};
return rate_fit(T9,a1) + rate_fit(T9,a2); return rate_fit(T9,a1) + rate_fit(T9,a2);
} }
// he3 + he3 -> he4 + 2p // he3 + he3 -> he4 + 2p
double he3he3_rate(const vec7 &T9){ double he3he3_rate(const vec7 &T9){
vec7 a = {2.477880e+01,0,-12.27700,-0.1036990,-6.499670e-02,1.681910e-02,-6.666670e-01}; 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); return rate_fit(T9,a);
} }
// he3(he3,2p)he4 // he3(he3,2p)he4
double he3he4_rate(const vec7 &T9){ double he3he4_rate(const vec7 &T9){
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 a1 = {1.560990e+01,0.000000e+00,-1.282710e+01,-3.082250e-02,-6.546850e-01,8.963310e-02,-6.666670e-01};
vec7 a2 = {1.770750e+01,0.000000e+00,-1.282710e+01,-3.812600e+00,9.422850e-02,-3.010180e-03,1.333330e+00}; 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); return rate_fit(T9,a1) + rate_fit(T9,a2);
} }
// he4 + he4 + he4 -> c12 // he4 + he4 + he4 -> c12
double triple_alpha_rate(const vec7 &T9){ double triple_alpha_rate(const vec7 &T9){
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 a1 = {-9.710520e-01,0.000000e+00,-3.706000e+01,2.934930e+01,-1.155070e+02,-1.000000e+01,-1.333330e+00};
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 a2 = {-1.178840e+01,-1.024460e+00,-2.357000e+01,2.048860e+01,-1.298820e+01,-2.000000e+01,-2.166670e+00};
vec7 a3 = {-2.435050e+01,-4.126560e+00,-1.349000e+01,2.142590e+01,-1.347690e+00,8.798160e-02,-1.316530e+01}; 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); return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
} }
// c12 + p -> n13 // c12 + p -> n13
double c12p_rate(const vec7 &T9){ double c12p_rate(const vec7 &T9){
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 a1={1.714820e+01,0.000000e+00,-1.369200e+01,-2.308810e-01,4.443620e+00,-3.158980e+00,-6.666670e-01};
vec7 a2={1.754280e+01,-3.778490e+00,-5.107350e+00,-2.241110e+00,1.488830e-01,0.000000e+00,-1.500000e+00}; 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); return rate_fit(T9,a1) + rate_fit(T9,a2);
} }
// c12 + he4 -> o16 // c12 + he4 -> o16
double c12a_rate(const vec7 &T9){ double c12a_rate(const vec7 &T9){
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 a1={6.965260e+01,-1.392540e+00,5.891280e+01,-1.482730e+02,9.083240e+00,-5.410410e-01,7.035540e+01};
vec7 a2={2.546340e+02,-1.840970e+00,1.034110e+02,-4.205670e+02,6.408740e+01,-1.246240e+01,1.373030e+02}; 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); return rate_fit(T9,a1) + rate_fit(T9,a2);
} }
// n14(p,g)o15 - o15 + p -> c12 + he4 // n14(p,g)o15 - o15 + p -> c12 + he4
double n14p_rate(const vec7 &T9){ double n14p_rate(const vec7 &T9){
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 a1={1.701000e+01,0.000000e+00,-1.519300e+01,-1.619540e-01,-7.521230e+00,-9.875650e-01,-6.666670e-01};
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 a2={2.011690e+01,0.000000e+00,-1.519300e+01,-4.639750e+00,9.734580e+00,-9.550510e+00,3.333330e-01};
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 a3={7.654440e+00,-2.998000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
vec7 a4={6.735780e+00,-4.891000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,6.820000e-02}; 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); 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 // n14(a,g)f18 assumed to go on to ne20
double n14a_rate(const vec7 &T9){ double n14a_rate(const vec7 &T9){
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 a1={2.153390e+01,0.000000e+00,-3.625040e+01,0.000000e+00,0.000000e+00,-5.000000e+00,-6.666670e-01};
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 a2={1.968380e-01,-5.160340e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
vec7 a3={1.389950e+01,-1.096560e+01,-5.622700e+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); return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
} }
// n15(p,a)c12 (CNO I) // n15(p,a)c12 (CNO I)
double n15pa_rate(const vec7 &T9){ double n15pa_rate(const vec7 &T9){
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 a1 = {2.747640e+01,0.000000e+00,-1.525300e+01,1.593180e+00,2.447900e+00,-2.197080e+00,-6.666670e-01};
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 a2 = {-4.873470e+00,-2.021170e+00,0.000000e+00,3.084970e+01,-8.504330e+00,-1.544260e+00,-1.500000e+00};
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 a3 = {2.089720e+01,-7.406000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
vec7 a4 = {-6.575220e+00,-1.163800e+00,0.000000e+00,2.271050e+01,-2.907070e+00,2.057540e-01,-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); return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3) + rate_fit(T9,a4);
} }
// n15(p,g)o16 (CNO II) // n15(p,g)o16 (CNO II)
double n15pg_rate(const vec7 &T9){ double n15pg_rate(const vec7 &T9){
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 a1 = {2.001760e+01,0.000000e+00,-1.524000e+01,3.349260e-01,4.590880e+00,-4.784680e+00,-6.666670e-01};
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 a2 = {6.590560e+00,-2.923150e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
vec7 a3 = {1.454440e+01,-1.022950e+01,0.000000e+00,0.000000e+00,4.590370e-02,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); return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
} }
double n15pg_frac(const vec7 &T9){ double n15pg_frac(const vec7 &T9){
double f1=n15pg_rate(T9); const double f1=n15pg_rate(T9);
double f2=n15pa_rate(T9); const double f2=n15pa_rate(T9);
return f1/(f1+f2); return f1/(f1+f2);
} }
// o16(p,g)f17 then f17 -> o17(p,a)n14 // o16(p,g)f17 then f17 -> o17(p,a)n14
double o16p_rate(const vec7 &T9){ double o16p_rate(const vec7 &T9){
vec7 a={1.909040e+01,0.000000e+00,-1.669600e+01,-1.162520e+00,2.677030e-01,-3.384110e-02,-6.666670e-01}; 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); return rate_fit(T9,a);
} }
// o16(a,g)ne20 // o16(a,g)ne20
double o16a_rate(const vec7 &T9){ double o16a_rate(const vec7 &T9){
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 a1={2.390300e+01,0.000000e+00,-3.972620e+01,-2.107990e-01,4.428790e-01,-7.977530e-02,-6.666670e-01};
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 a2={3.885710e+00,-1.035850e+01,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
vec7 a3={9.508480e+00,-1.276430e+01,0.000000e+00,-3.659250e+00,7.142240e-01,-1.075080e-03,-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); return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
} }
// ne20(a,g)mg24 // ne20(a,g)mg24
double ne20a_rate(const vec7 &T9){ double ne20a_rate(const vec7 &T9){
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 a1={2.450580e+01,0.000000e+00,-4.625250e+01,5.589010e+00,7.618430e+00,-3.683000e+00,-6.666670e-01};
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 a2={-3.870550e+01,-2.506050e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
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 a3={1.983070e+00,-9.220260e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
vec7 a4={-8.798270e+00,-1.278090e+01,0.000000e+00,1.692290e+01,-2.573250e+00,2.089970e-01,-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); return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3) + rate_fit(T9,a4);
} }
// c12(c12,a)ne20 // c12(c12,a)ne20
double c12c12_rate(const vec7 &T9){ double c12c12_rate(const vec7 &T9){
vec7 a={6.128630e+01,0.000000e+00,-8.416500e+01,-1.566270e+00,-7.360840e-02,-7.279700e-02,-6.666670e-01}; 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); return rate_fit(T9,a);
} }
// c12(o16,a)mg24 // c12(o16,a)mg24
double c12o16_rate(const vec7 &T9){ double c12o16_rate(const vec7 &T9){
vec7 a={4.853410e+01,3.720400e-01,-1.334130e+02,5.015720e+01,-3.159870e+00,1.782510e-02,-2.370270e+01}; 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); return rate_fit(T9,a);
} }
@@ -244,209 +244,211 @@ namespace serif::network::approx8{
// a Jacobian matrix for implicit solvers // a Jacobian matrix for implicit solvers
void Jacobian::operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) { void Jacobian::operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) const {
serif::constant::Constants& constants = serif::constant::Constants::getInstance(); serif::constant::Constants& constants = serif::constant::Constants::getInstance();
const double avo = constants.get("N_a").value; const double avo = constants.get("N_a").value;
const double clight = constants.get("c").value; const double clight = constants.get("c").value;
// EOS // EOS
vec7 T9=get_T9_array(y[Net::itemp]); const vec7 T9=get_T9_array(y[Approx8Net::iTemp]);
// evaluate rates once per call // evaluate rates once per call
double rpp=pp_rate(T9); const double rpp=pp_rate(T9);
double r33=he3he3_rate(T9); const double r33=he3he3_rate(T9);
double r34=he3he4_rate(T9); const double r34=he3he4_rate(T9);
double r3a=triple_alpha_rate(T9); const double r3a=triple_alpha_rate(T9);
double rc12p=c12p_rate(T9); const double rc12p=c12p_rate(T9);
double rc12a=c12a_rate(T9); const double rc12a=c12a_rate(T9);
double rn14p=n14p_rate(T9); const double rn14p=n14p_rate(T9);
double rn14a=n14a_rate(T9); const double rn14a=n14a_rate(T9);
double ro16p=o16p_rate(T9); const double ro16p=o16p_rate(T9);
double ro16a=o16a_rate(T9); const double ro16a=o16a_rate(T9);
double rne20a=ne20a_rate(T9); const double rne20a=ne20a_rate(T9);
double r1212=c12c12_rate(T9); const double r1212=c12c12_rate(T9);
double r1216=c12o16_rate(T9); const double r1216=c12o16_rate(T9);
double pfrac=n15pg_frac(T9); const double pFrac=n15pg_frac(T9);
double afrac=1-pfrac; const double aFrac=1-pFrac;
double yh1 = y[Net::ih1]; const double yh1 = y[Approx8Net::ih1];
double yhe3 = y[Net::ihe3]; const double yhe3 = y[Approx8Net::ihe3];
double yhe4 = y[Net::ihe4]; const double yhe4 = y[Approx8Net::ihe4];
double yc12 = y[Net::ic12]; const double yc12 = y[Approx8Net::ic12];
double yn14 = y[Net::in14]; const double yn14 = y[Approx8Net::in14];
double yo16 = y[Net::io16]; const double yo16 = y[Approx8Net::io16];
double yne20 = y[Net::ine20]; const double yne20 = y[Approx8Net::ine20];
// zero all elements to begin // zero all elements to begin
for (int i=0; i < Net::nvar; i++) { for (int i=0; i < Approx8Net::nVar; i++) {
for (int j=0; j < Net::nvar; j++) { for (int j=0; j < Approx8Net::nVar; j++) {
J(i,j)=0.0; J(i,j)=0.0;
} }
} }
// h1 jacobian elements // h1 jacobian elements
J(Net::ih1,Net::ih1) = -3*yh1*rpp - 2*yc12*rc12p -2*yn14*rn14p -2*yo16*ro16p; J(Approx8Net::ih1,Approx8Net::ih1) = -3*yh1*rpp - 2*yc12*rc12p -2*yn14*rn14p -2*yo16*ro16p;
J(Net::ih1,Net::ihe3) = 2*yhe3*r33 - yhe4*r34; J(Approx8Net::ih1,Approx8Net::ihe3) = 2*yhe3*r33 - yhe4*r34;
J(Net::ih1,Net::ihe4) = -yhe3*r34; J(Approx8Net::ih1,Approx8Net::ihe4) = -yhe3*r34;
J(Net::ih1,Net::ic12) = -2*yh1*rc12p; J(Approx8Net::ih1,Approx8Net::ic12) = -2*yh1*rc12p;
J(Net::ih1,Net::in14) = -2*yh1*rn14p; J(Approx8Net::ih1,Approx8Net::in14) = -2*yh1*rn14p;
J(Net::ih1,Net::io16) = -2*yh1*ro16p; J(Approx8Net::ih1,Approx8Net::io16) = -2*yh1*ro16p;
// he3 jacobian elements // he3 jacobian elements
J(Net::ihe3,Net::ih1) = yh1*rpp; J(Approx8Net::ihe3,Approx8Net::ih1) = yh1*rpp;
J(Net::ihe3,Net::ihe3) = -2*yhe3*r33 - yhe4*r34; J(Approx8Net::ihe3,Approx8Net::ihe3) = -2*yhe3*r33 - yhe4*r34;
J(Net::ihe3,Net::ihe4) = -yhe3*r34; J(Approx8Net::ihe3,Approx8Net::ihe4) = -yhe3*r34;
// he4 jacobian elements // he4 jacobian elements
J(Net::ihe4,Net::ih1) = yn14*afrac*rn14p + yo16*ro16p; J(Approx8Net::ihe4,Approx8Net::ih1) = yn14*aFrac*rn14p + yo16*ro16p;
J(Net::ihe4,Net::ihe3) = yhe3*r33 - yhe4*r34; J(Approx8Net::ihe4,Approx8Net::ihe3) = yhe3*r33 - yhe4*r34;
J(Net::ihe4,Net::ihe4) = yhe3*r34 - 1.5*yhe4*yhe4*r3a - yc12*rc12a - 1.5*yn14*rn14a - yo16*ro16a - yne20*rne20a; J(Approx8Net::ihe4,Approx8Net::ihe4) = yhe3*r34 - 1.5*yhe4*yhe4*r3a - yc12*rc12a - 1.5*yn14*rn14a - yo16*ro16a - yne20*rne20a;
J(Net::ihe4,Net::ic12) = -yhe4*rc12a + yc12*r1212 + yo16*r1216; J(Approx8Net::ihe4,Approx8Net::ic12) = -yhe4*rc12a + yc12*r1212 + yo16*r1216;
J(Net::ihe4,Net::in14) = yh1*afrac*rn14p - 1.5*yhe4*rn14a; J(Approx8Net::ihe4,Approx8Net::in14) = yh1*aFrac*rn14p - 1.5*yhe4*rn14a;
J(Net::ihe4,Net::io16) = yh1*ro16p - yhe4*ro16a + yc12*r1216; J(Approx8Net::ihe4,Approx8Net::io16) = yh1*ro16p - yhe4*ro16a + yc12*r1216;
J(Net::ihe4,Net::ine20) = -yhe4*rne20a; J(Approx8Net::ihe4,Approx8Net::ine20) = -yhe4*rne20a;
// c12 jacobian elements // c12 jacobian elements
J(Net::ic12,Net::ih1) = -yc12*rc12p + yn14*afrac*rn14p; J(Approx8Net::ic12,Approx8Net::ih1) = -yc12*rc12p + yn14*aFrac*rn14p;
J(Net::ic12,Net::ihe4) = 0.5*yhe4*yhe4*r3a - yhe4*rc12a; J(Approx8Net::ic12,Approx8Net::ihe4) = 0.5*yhe4*yhe4*r3a - yhe4*rc12a;
J(Net::ic12,Net::ic12) = -yh1*rc12p - yhe4*rc12a - yo16*r1216 - 2*yc12*r1212; J(Approx8Net::ic12,Approx8Net::ic12) = -yh1*rc12p - yhe4*rc12a - yo16*r1216 - 2*yc12*r1212;
J(Net::ic12,Net::in14) = yh1*yn14*afrac*rn14p; J(Approx8Net::ic12,Approx8Net::in14) = yh1*yn14*aFrac*rn14p;
J(Net::ic12,Net::io16) = -yc12*r1216; J(Approx8Net::ic12,Approx8Net::io16) = -yc12*r1216;
// n14 jacobian elements // n14 jacobian elements
J(Net::in14,Net::ih1) = yc12*rc12p - yn14*rn14p + yo16*ro16p; J(Approx8Net::in14,Approx8Net::ih1) = yc12*rc12p - yn14*rn14p + yo16*ro16p;
J(Net::in14,Net::ihe4) = -yn14*rn14a; J(Approx8Net::in14,Approx8Net::ihe4) = -yn14*rn14a;
J(Net::in14,Net::ic12) = yh1*rc12p; J(Approx8Net::in14,Approx8Net::ic12) = yh1*rc12p;
J(Net::in14,Net::in14) = -yh1*rn14p - yhe4*rn14a; J(Approx8Net::in14,Approx8Net::in14) = -yh1*rn14p - yhe4*rn14a;
J(Net::in14,Net::io16) = yo16*ro16p; J(Approx8Net::in14,Approx8Net::io16) = yo16*ro16p;
// o16 jacobian elements // o16 jacobian elements
J(Net::io16,Net::ih1) = yn14*pfrac*rn14p - yo16*ro16p; J(Approx8Net::io16,Approx8Net::ih1) = yn14*pFrac*rn14p - yo16*ro16p;
J(Net::io16,Net::ihe4) = yc12*rc12a - yo16*ro16a; J(Approx8Net::io16,Approx8Net::ihe4) = yc12*rc12a - yo16*ro16a;
J(Net::io16,Net::ic12) = yhe4*rc12a - yo16*r1216; J(Approx8Net::io16,Approx8Net::ic12) = yhe4*rc12a - yo16*r1216;
J(Net::io16,Net::in14) = yh1*pfrac*rn14p; J(Approx8Net::io16,Approx8Net::in14) = yh1*pFrac*rn14p;
J(Net::io16,Net::io16) = yh1*ro16p - yc12*r1216 -yhe4*ro16a; J(Approx8Net::io16,Approx8Net::io16) = yh1*ro16p - yc12*r1216 -yhe4*ro16a;
// ne20 jacobian elements // ne20 jacobian elements
J(Net::ine20,Net::ihe4) = yn14*rn14a + yo16*ro16a - yne20*rne20a; J(Approx8Net::ine20,Approx8Net::ihe4) = yn14*rn14a + yo16*ro16a - yne20*rne20a;
J(Net::ine20,Net::ic12) = yc12*r1212; J(Approx8Net::ine20,Approx8Net::ic12) = yc12*r1212;
J(Net::ine20,Net::in14) = yhe4*rn14a; J(Approx8Net::ine20,Approx8Net::in14) = yhe4*rn14a;
J(Net::ine20,Net::io16) = yo16*ro16a; J(Approx8Net::ine20,Approx8Net::io16) = yo16*ro16a;
J(Net::ine20,Net::ine20) = -yhe4*rne20a; J(Approx8Net::ine20,Approx8Net::ine20) = -yhe4*rne20a;
// mg24 jacobian elements // mg24 jacobian elements
J(Net::img24,Net::ihe4) = yne20*rne20a; J(Approx8Net::img24,Approx8Net::ihe4) = yne20*rne20a;
J(Net::img24,Net::ic12) = yo16*r1216; J(Approx8Net::img24,Approx8Net::ic12) = yo16*r1216;
J(Net::img24,Net::io16) = yc12*r1216; J(Approx8Net::img24,Approx8Net::io16) = yc12*r1216;
J(Net::img24,Net::ine20) = yhe4*rne20a; J(Approx8Net::img24,Approx8Net::ine20) = yhe4*rne20a;
// energy accountings // energy accounting
for (int j=0; j<Net::niso; j++) { for (int j=0; j<Approx8Net::nIso; j++) {
for (int i=0; i<Net::niso; i++) { for (int i=0; i<Approx8Net::nIso; i++) {
J(Net::iener,j) += J(i,j)*Net::mion[i]; J(Approx8Net::iEnergy,j) += J(i,j)*Approx8Net::mIon[i];
} }
J(Net::iener,j) *= -avo*clight*clight; J(Approx8Net::iEnergy,j) *= -avo*clight*clight;
} }
} }
void ODE::operator() ( const vector_type &y, vector_type &dydt, double /* t */) { void ODE::operator() ( const vector_type &y, vector_type &dydt, double /* t */) const {
serif::constant::Constants& constants = serif::constant::Constants::getInstance(); const serif::constant::Constants& constants = serif::constant::Constants::getInstance();
const double avo = constants.get("N_a").value; const double avo = constants.get("N_a").value;
const double clight = constants.get("c").value; const double clight = constants.get("c").value;
// EOS // EOS
double T=y[Net::itemp]; const double T = y[Approx8Net::iTemp];
double den=y[Net::iden]; const double den = y[Approx8Net::iDensity];
vec7 T9=get_T9_array(T); const vec7 T9=get_T9_array(T);
// rates // rates
double rpp=den*pp_rate(T9); const double rpp=den*pp_rate(T9);
double r33=den*he3he3_rate(T9); const double r33=den*he3he3_rate(T9);
double r34=den*he3he4_rate(T9); const double r34=den*he3he4_rate(T9);
double r3a=den*den*triple_alpha_rate(T9); const double r3a=den*den*triple_alpha_rate(T9);
double rc12p=den*c12p_rate(T9); const double rc12p=den*c12p_rate(T9);
double rc12a=den*c12a_rate(T9); const double rc12a=den*c12a_rate(T9);
double rn14p=den*n14p_rate(T9); const double rn14p=den*n14p_rate(T9);
double rn14a=n14a_rate(T9); const double rn14a=n14a_rate(T9);
double ro16p=den*o16p_rate(T9); const double ro16p=den*o16p_rate(T9);
double ro16a=den*o16a_rate(T9); const double ro16a=den*o16a_rate(T9);
double rne20a=den*ne20a_rate(T9); const double rne20a=den*ne20a_rate(T9);
double r1212=den*c12c12_rate(T9); const double r1212=den*c12c12_rate(T9);
double r1216=den*c12o16_rate(T9); const double r1216=den*c12o16_rate(T9);
double pfrac=n15pg_frac(T9); const double pFrac=n15pg_frac(T9);
double afrac=1-pfrac; const double aFrac=1-pFrac;
double yh1 = y[Net:: ih1]; const double yh1 = y[Approx8Net:: ih1];
double yhe3 = y[Net:: ihe3]; const double yhe3 = y[Approx8Net:: ihe3];
double yhe4 = y[Net:: ihe4]; const double yhe4 = y[Approx8Net:: ihe4];
double yc12 = y[Net:: ic12]; const double yc12 = y[Approx8Net:: ic12];
double yn14 = y[Net:: in14]; const double yn14 = y[Approx8Net:: in14];
double yo16 = y[Net:: io16]; const double yo16 = y[Approx8Net:: io16];
double yne20 = y[Net::ine20]; const double yne20 = y[Approx8Net::ine20];
dydt[Net::ih1] = -1.5*yh1*yh1*rpp; dydt[Approx8Net::ih1] = -1.5*yh1*yh1*rpp;
dydt[Net::ih1] += yhe3*yhe3*r33; dydt[Approx8Net::ih1] += yhe3*yhe3*r33;
dydt[Net::ih1] += -yhe3*yhe4*r34; dydt[Approx8Net::ih1] += -yhe3*yhe4*r34;
dydt[Net::ih1] += -2*yh1*yc12*rc12p; dydt[Approx8Net::ih1] += -2*yh1*yc12*rc12p;
dydt[Net::ih1] += -2*yh1*yn14*rn14p; dydt[Approx8Net::ih1] += -2*yh1*yn14*rn14p;
dydt[Net::ih1] += -2*yh1*yo16*ro16p; dydt[Approx8Net::ih1] += -2*yh1*yo16*ro16p;
dydt[Net::ihe3] = 0.5*yh1*yh1*rpp; dydt[Approx8Net::ihe3] = 0.5*yh1*yh1*rpp;
dydt[Net::ihe3] += -yhe3*yhe3*r33; dydt[Approx8Net::ihe3] += -yhe3*yhe3*r33;
dydt[Net::ihe3] += -yhe3*yhe4*r34; dydt[Approx8Net::ihe3] += -yhe3*yhe4*r34;
dydt[Net::ihe4] = 0.5*yhe3*yhe3*r33; dydt[Approx8Net::ihe4] = 0.5*yhe3*yhe3*r33;
dydt[Net::ihe4] += yhe3*yhe4*r34; dydt[Approx8Net::ihe4] += yhe3*yhe4*r34;
dydt[Net::ihe4] += -yhe4*yc12*rc12a; dydt[Approx8Net::ihe4] += -yhe4*yc12*rc12a;
dydt[Net::ihe4] += yh1*yn14*afrac*rn14p; dydt[Approx8Net::ihe4] += yh1*yn14*aFrac*rn14p;
dydt[Net::ihe4] += yh1*yo16*ro16p; dydt[Approx8Net::ihe4] += yh1*yo16*ro16p;
dydt[Net::ihe4] += -0.5*yhe4*yhe4*yhe4*r3a; dydt[Approx8Net::ihe4] += -0.5*yhe4*yhe4*yhe4*r3a;
dydt[Net::ihe4] += -yhe4*yo16*ro16a; dydt[Approx8Net::ihe4] += -yhe4*yo16*ro16a;
dydt[Net::ihe4] += 0.5*yc12*yc12*r1212; dydt[Approx8Net::ihe4] += 0.5*yc12*yc12*r1212;
dydt[Net::ihe4] += yc12*yo16*r1216; dydt[Approx8Net::ihe4] += yc12*yo16*r1216;
dydt[Net::ihe4] += -yhe4*yne20*rne20a; dydt[Approx8Net::ihe4] += -yhe4*yne20*rne20a;
dydt[Net::ic12] = (1./6.)*yhe4*yhe4*yhe4*r3a; dydt[Approx8Net::ic12] = (1./6.)*yhe4*yhe4*yhe4*r3a;
dydt[Net::ic12] += -yhe4*yc12*rc12a; dydt[Approx8Net::ic12] += -yhe4*yc12*rc12a;
dydt[Net::ic12] += -yh1*yc12*rc12p; dydt[Approx8Net::ic12] += -yh1*yc12*rc12p;
dydt[Net::ic12] += yh1*yn14*afrac*rn14p; dydt[Approx8Net::ic12] += yh1*yn14*aFrac*rn14p;
dydt[Net::ic12] += -yc12*yc12*r1212; dydt[Approx8Net::ic12] += -yc12*yc12*r1212;
dydt[Net::ic12] += -yc12*yo16*r1216; dydt[Approx8Net::ic12] += -yc12*yo16*r1216;
dydt[Net::in14] = yh1*yc12*rc12p; dydt[Approx8Net::in14] = yh1*yc12*rc12p;
dydt[Net::in14] += -yh1*yn14*rn14p; dydt[Approx8Net::in14] += -yh1*yn14*rn14p;
dydt[Net::in14] += yh1*yo16*ro16p; dydt[Approx8Net::in14] += yh1*yo16*ro16p;
dydt[Net::in14] += -yhe4*yn14*rn14a; dydt[Approx8Net::in14] += -yhe4*yn14*rn14a;
dydt[Net::io16] = yhe4*yc12*rc12a; dydt[Approx8Net::io16] = yhe4*yc12*rc12a;
dydt[Net::io16] += yh1*yn14*pfrac*rn14p; dydt[Approx8Net::io16] += yh1*yn14*pFrac*rn14p;
dydt[Net::io16] += -yh1*yo16*ro16p; dydt[Approx8Net::io16] += -yh1*yo16*ro16p;
dydt[Net::io16] += -yc12*yo16*r1216; dydt[Approx8Net::io16] += -yc12*yo16*r1216;
dydt[Net::io16] += -yhe4*yo16*ro16a; dydt[Approx8Net::io16] += -yhe4*yo16*ro16a;
dydt[Net::ine20] = 0.5*yc12*yc12*r1212; dydt[Approx8Net::ine20] = 0.5*yc12*yc12*r1212;
dydt[Net::ine20] += yhe4*yn14*rn14a; dydt[Approx8Net::ine20] += yhe4*yn14*rn14a;
dydt[Net::ine20] += yhe4*yo16*ro16a; dydt[Approx8Net::ine20] += yhe4*yo16*ro16a;
dydt[Net::ine20] += -yhe4*yne20*rne20a; dydt[Approx8Net::ine20] += -yhe4*yne20*rne20a;
dydt[Net::img24] = yc12*yo16*r1216; dydt[Approx8Net::img24] = yc12*yo16*r1216;
dydt[Net::img24] += yhe4*yne20*rne20a; dydt[Approx8Net::img24] += yhe4*yne20*rne20a;
dydt[Net::itemp] = 0.; dydt[Approx8Net::iTemp] = 0.;
dydt[Net::iden] = 0.; dydt[Approx8Net::iDensity] = 0.;
// energy accounting // energy accounting
double enuc = 0.; double eNuc = 0.;
for (int i=0; i<Net::niso; i++) { for (int i=0; i<Approx8Net::nIso; i++) {
enuc += Net::mion[i]*dydt[i]; eNuc += Approx8Net::mIon[i]*dydt[i];
} }
dydt[Net::iener] = -enuc*avo*clight*clight; dydt[Approx8Net::iEnergy] = -eNuc*avo*clight*clight;
} }
Approx8Network::Approx8Network() : Network(APPROX8) {}
NetOut Approx8Network::evaluate(const NetIn &netIn) { NetOut Approx8Network::evaluate(const NetIn &netIn) {
m_y = convert_netIn(netIn); m_y = convert_netIn(netIn);
m_tmax = netIn.tmax; m_tMax = netIn.tMax;
m_dt0 = netIn.dt0; m_dt0 = netIn.dt0;
const double stiff_abs_tol = m_config.get<double>("Network:Approx8:Stiff:AbsTol", 1.0e-6); const double stiff_abs_tol = m_config.get<double>("Network:Approx8:Stiff:AbsTol", 1.0e-6);
@@ -463,7 +465,7 @@ namespace serif::network::approx8{
std::make_pair(ODE(), Jacobian()), std::make_pair(ODE(), Jacobian()),
m_y, m_y,
0.0, 0.0,
m_tmax, m_tMax,
m_dt0 m_dt0
); );
@@ -474,31 +476,33 @@ namespace serif::network::approx8{
ODE(), ODE(),
m_y, m_y,
0.0, 0.0,
m_tmax, m_tMax,
m_dt0 m_dt0
); );
} }
double ysum = 0.0; double ySum = 0.0;
for (int i = 0; i < Net::niso; i++) { for (int i = 0; i < Approx8Net::nIso; i++) {
m_y[i] *= Net::aion[i]; m_y[i] *= Approx8Net::aIon[i];
ysum += m_y[i]; ySum += m_y[i];
} }
for (int i = 0; i < Net::niso; i++) { for (int i = 0; i < Approx8Net::nIso; i++) {
m_y[i] /= ysum; m_y[i] /= ySum;
} }
NetOut netOut; NetOut netOut;
std::vector<double> outComposition; std::vector<double> outComposition;
outComposition.reserve(Net::nvar); outComposition.reserve(Approx8Net::nVar);
for (int i = 0; i < Net::niso; i++) { for (int i = 0; i < Approx8Net::nIso; i++) {
outComposition.push_back(m_y[i]); outComposition.push_back(m_y[i]);
} }
netOut.energy = m_y[Net::iener]; netOut.energy = m_y[Approx8Net::iEnergy];
netOut.composition = outComposition;
netOut.num_steps = num_steps; 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 = serif::composition::Composition(symbols, outComposition);
return netOut; return netOut;
} }
@@ -507,31 +511,26 @@ namespace serif::network::approx8{
} }
vector_type Approx8Network::convert_netIn(const NetIn &netIn) { vector_type Approx8Network::convert_netIn(const NetIn &netIn) {
if (netIn.composition.size() != Net::niso) { vector_type y(Approx8Net::nVar, 0.0);
LOG_ERROR(m_logger, "Error: composition size mismatch in convert_netIn"); y[Approx8Net::ih1] = netIn.composition.getMassFraction("H-1");
throw std::runtime_error("Error: composition size mismatch in convert_netIn"); y[Approx8Net::ihe3] = netIn.composition.getMassFraction("He-3");
} y[Approx8Net::ihe4] = netIn.composition.getMassFraction("He-4");
y[Approx8Net::ic12] = netIn.composition.getMassFraction("C-12");
y[Approx8Net::in14] = netIn.composition.getMassFraction("N-14");
y[Approx8Net::io16] = netIn.composition.getMassFraction("O-16");
y[Approx8Net::ine20] = netIn.composition.getMassFraction("Ne-20");
y[Approx8Net::img24] = netIn.composition.getMassFraction("Mg-24");
y[Approx8Net::iTemp] = netIn.temperature;
y[Approx8Net::iDensity] = netIn.density;
y[Approx8Net::iEnergy] = netIn.energy;
vector_type y(Net::nvar, 0.0); double ySum = 0.0;
y[Net::ih1] = netIn.composition[0]; for (int i = 0; i < Approx8Net::nIso; i++) {
y[Net::ihe3] = netIn.composition[1]; y[i] /= Approx8Net::aIon[i];
y[Net::ihe4] = netIn.composition[2]; ySum += y[i];
y[Net::ic12] = netIn.composition[3];
y[Net::in14] = netIn.composition[4];
y[Net::io16] = netIn.composition[5];
y[Net::ine20] = netIn.composition[6];
y[Net::img24] = netIn.composition[7];
y[Net::itemp] = netIn.temperature;
y[Net::iden] = netIn.density;
y[Net::iener] = netIn.energy;
double ysum = 0.0;
for (int i = 0; i < Net::niso; i++) {
y[i] /= Net::aion[i];
ysum += y[i];
} }
for (int i = 0; i < Net::niso; i++) { for (int i = 0; i < Approx8Net::nIso; i++) {
y[i] /= ysum; y[i] /= ySum;
} }
return y; return y;

42
src/network/private/network.cpp
View File

@@ -19,18 +19,50 @@
// //
// *********************************************************************** */ // *********************************************************************** */
#include "network.h" #include "network.h"
#include "approx8.h"
#include "probe.h" #include "probe.h"
#include "quill/LogMacros.h" #include "quill/LogMacros.h"
namespace serif::network { namespace serif::network {
Network::Network() : Network::Network(const NetworkFormat format) :
m_config(serif::config::Config::getInstance()), m_config(serif::config::Config::getInstance()),
m_logManager(serif::probe::LogManager::getInstance()), m_logManager(serif::probe::LogManager::getInstance()),
m_logger(m_logManager.getLogger("log")) { m_logger(m_logManager.getLogger("log")),
m_format(format) {
if (format == NetworkFormat::UNKNOWN) {
LOG_ERROR(m_logger, "nuclearNetwork::Network::Network() called with UNKNOWN format");
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;
}
NetOut Network::evaluate(const NetIn &netIn) { NetOut Network::evaluate(const NetIn &netIn) {
// You can throw an exception here or log a warning if it should never be used. NetOut netOut;
LOG_ERROR(m_logger, "nuclearNetwork::Network::evaluate() is not implemented"); switch (m_format) {
throw std::runtime_error("nuclearNetwork::Network::evaluate() is not implemented"); case APPROX8: {
approx8::Approx8Network network;
netOut = network.evaluate(netIn);
break;
}
case UNKNOWN: {
LOG_ERROR(m_logger, "Network format {} is not implemented.", FormatStringLookup.at(m_format));
throw std::runtime_error("Network format not implemented.");
}
default: {
LOG_ERROR(m_logger, "Unknown network format.");
throw std::runtime_error("Unknown network format.");
}
}
return netOut;
} }
} }

62
src/network/public/approx8.h
View File

@@ -31,10 +31,13 @@
* @brief Header file for the Approx8 nuclear reaction network. * @brief Header file for the Approx8 nuclear reaction network.
* *
* This file contains the definitions and declarations for the Approx8 nuclear reaction network. * This file contains the definitions and declarations for the Approx8 nuclear reaction network.
* The network is based on Frank Timmes' "aprox8" and includes 8 isotopes and various nuclear reactions. * The network is based on Frank Timmes' "approx8" and includes 8 isotopes and various nuclear reactions.
* The rates are evaluated using a fitting function with coefficients from reaclib.jinaweb.org. * The rates are evaluated using a fitting function with coefficients from reaclib.jinaweb.org.
*/ */
namespace serif::network::approx8{
/** /**
* @typedef vector_type * @typedef vector_type
* @brief Alias for a vector of doubles using Boost uBLAS. * @brief Alias for a vector of doubles using Boost uBLAS.
@@ -53,32 +56,30 @@ typedef boost::numeric::ublas::matrix< double > matrix_type;
*/ */
typedef std::array<double,7> vec7; typedef std::array<double,7> vec7;
namespace serif::network::approx8{
using namespace boost::numeric::odeint; using namespace boost::numeric::odeint;
/** /**
* @struct Net * @struct Approx8Net
* @brief Contains constants and arrays related to the nuclear network. * @brief Contains constants and arrays related to the nuclear network.
*/ */
struct Net{ struct Approx8Net{
const static int ih1=0; static constexpr int ih1=0;
const static int ihe3=1; static constexpr int ihe3=1;
const static int ihe4=2; static constexpr int ihe4=2;
const static int ic12=3; static constexpr int ic12=3;
const static int in14=4; static constexpr int in14=4;
const static int io16=5; static constexpr int io16=5;
const static int ine20=6; static constexpr int ine20=6;
const static int img24=7; static constexpr int img24=7;
const static int itemp=img24+1; static constexpr int iTemp=img24+1;
const static int iden =itemp+1; static constexpr int iDensity =iTemp+1;
const static int iener=iden+1; static constexpr int iEnergy=iDensity+1;
const static int niso=img24+1; // number of isotopes static constexpr int nIso=img24+1; // number of isotopes
const static int nvar=iener+1; // number of variables static constexpr int nVar=iEnergy+1; // number of variables
static constexpr std::array<int,niso> aion = { static constexpr std::array<int,nIso> aIon = {
1, 1,
3, 3,
4, 4,
@@ -89,7 +90,7 @@ namespace serif::network::approx8{
24 24
}; };
static constexpr std::array<double,niso> mion = { static constexpr std::array<double,nIso> mIon = {
1.67262164e-24, 1.67262164e-24,
5.00641157e-24, 5.00641157e-24,
6.64465545e-24, 6.64465545e-24,
@@ -270,10 +271,9 @@ namespace serif::network::approx8{
* @brief Calculates the Jacobian matrix. * @brief Calculates the Jacobian matrix.
* @param y State vector. * @param y State vector.
* @param J Jacobian matrix. * @param J Jacobian matrix.
* @param t Time.
* @param dfdt Derivative of the state vector. * @param dfdt Derivative of the state vector.
*/ */
void operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ); void operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) const;
}; };
/** /**
@@ -285,23 +285,24 @@ namespace serif::network::approx8{
* @brief Calculates the derivatives of the state vector. * @brief Calculates the derivatives of the state vector.
* @param y State vector. * @param y State vector.
* @param dydt Derivative of the state vector. * @param dydt Derivative of the state vector.
* @param t Time.
*/ */
void operator() ( const vector_type &y, vector_type &dydt, double /* t */); void operator() ( const vector_type &y, vector_type &dydt, double /* t */) const;
}; };
/** /**
* @class Approx8Network * @class Approx8Network
* @brief Class for the Approx8 nuclear reaction network. * @brief Class for the Approx8 nuclear reaction network.
*/ */
class Approx8Network : public Network { class Approx8Network final : public Network {
public: public:
Approx8Network();
/** /**
* @brief Evaluates the nuclear network. * @brief Evaluates the nuclear network.
* @param netIn Input parameters for the network. * @param netIn Input parameters for the network.
* @return Output results from the network. * @return Output results from the network.
*/ */
virtual NetOut evaluate(const NetIn &netIn); NetOut evaluate(const NetIn &netIn) override;
/** /**
* @brief Sets whether the solver should use a stiff method. * @brief Sets whether the solver should use a stiff method.
@@ -313,11 +314,11 @@ namespace serif::network::approx8{
* @brief Checks if the solver is using a stiff method. * @brief Checks if the solver is using a stiff method.
* @return Boolean indicating if a stiff method is being used. * @return Boolean indicating if a stiff method is being used.
*/ */
bool isStiff() { return m_stiff; } bool isStiff() const { return m_stiff; }
private: private:
vector_type m_y; vector_type m_y;
double m_tmax; double m_tMax = 0;
double m_dt0; double m_dt0 = 0;
bool m_stiff = false; bool m_stiff = false;
/** /**
@@ -325,7 +326,8 @@ namespace serif::network::approx8{
* @param netIn Input parameters for the network. * @param netIn Input parameters for the network.
* @return Internal state vector. * @return Internal state vector.
*/ */
vector_type convert_netIn(const NetIn &netIn); static vector_type convert_netIn(const NetIn &netIn);
}; };
} // namespace nnApprox8 } // namespace nnApprox8

27
src/network/public/network.h
View File

@@ -25,9 +25,21 @@
#include "probe.h" #include "probe.h"
#include "config.h" #include "config.h"
#include "quill/Logger.h" #include "quill/Logger.h"
#include "composition.h"
#include <unordered_map>
namespace serif::network { namespace serif::network {
enum NetworkFormat {
APPROX8, ///< Approx8 nuclear reaction network format.
UNKNOWN,
};
static inline std::unordered_map<NetworkFormat, std::string> FormatStringLookup = {
{APPROX8, "Approx8"},
{UNKNOWN, "Unknown"}
};
/** /**
* @struct NetIn * @struct NetIn
* @brief Input structure for the network evaluation. * @brief Input structure for the network evaluation.
@@ -46,8 +58,8 @@ namespace serif::network {
* @endcode * @endcode
*/ */
struct NetIn { struct NetIn {
std::vector<double> composition; ///< Composition of the network serif::composition::Composition composition; ///< Composition of the network
double tmax; ///< Maximum time double tMax; ///< Maximum time
double dt0; ///< Initial time step double dt0; ///< Initial time step
double temperature; ///< Temperature in Kelvin double temperature; ///< Temperature in Kelvin
double density; ///< Density in g/cm^3 double density; ///< Density in g/cm^3
@@ -69,7 +81,7 @@ namespace serif::network {
* @endcode * @endcode
*/ */
struct NetOut { struct NetOut {
std::vector<double> composition; ///< Composition of the network after evaluation serif::composition::Composition composition; ///< Composition of the network after evaluation
int num_steps; ///< Number of steps taken in the evaluation int num_steps; ///< Number of steps taken in the evaluation
double energy; ///< Energy in ergs after evaluation double energy; ///< Energy in ergs after evaluation
}; };
@@ -90,9 +102,12 @@ namespace serif::network {
*/ */
class Network { class Network {
public: public:
Network(); explicit Network(const NetworkFormat format = NetworkFormat::APPROX8);
virtual ~Network() = default; virtual ~Network() = default;
NetworkFormat getFormat() const;
NetworkFormat setFormat(const NetworkFormat format);
/** /**
* @brief Evaluate the network based on the input parameters. * @brief Evaluate the network based on the input parameters.
* *
@@ -105,6 +120,10 @@ namespace serif::network {
serif::config::Config& m_config; ///< Configuration instance serif::config::Config& m_config; ///< Configuration instance
serif::probe::LogManager& m_logManager; ///< Log manager instance serif::probe::LogManager& m_logManager; ///< Log manager instance
quill::Logger* m_logger; ///< Logger instance quill::Logger* m_logger; ///< Logger instance
NetworkFormat m_format; ///< Format of the network
}; };
} // namespace nuclearNetwork } // namespace nuclearNetwork

23
tests/network/approx8Test.cpp
View File

@@ -4,6 +4,7 @@
#include "approx8.h" #include "approx8.h"
#include "config.h" #include "config.h"
#include "network.h" #include "network.h"
#include "composition.h"
#include <vector> #include <vector>
@@ -32,19 +33,31 @@ TEST_F(approx8Test, evaluate) {
serif::network::NetIn netIn; serif::network::NetIn netIn;
std::vector<double> comp = {0.708, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4}; std::vector<double> comp = {0.708, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4};
std::vector<std::string> symbols = {"H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"};
netIn.composition = comp; serif::composition::Composition composition;
composition.registerSymbol(symbols, true);
composition.setMassFraction(symbols, comp);
bool isFinalized = composition.finalize(true);
EXPECT_TRUE(isFinalized);
netIn.composition = composition;
netIn.temperature = 1e7; netIn.temperature = 1e7;
netIn.density = 1e2; netIn.density = 1e2;
netIn.energy = 0.0; netIn.energy = 0.0;
netIn.tmax = 3.15e17; netIn.tMax = 3.15e17;
netIn.dt0 = 1e12; netIn.dt0 = 1e12;
serif::network::NetOut netOut; serif::network::NetOut netOut;
EXPECT_NO_THROW(netOut = network.evaluate(netIn)); EXPECT_NO_THROW(netOut = network.evaluate(netIn));
EXPECT_DOUBLE_EQ(netOut.composition[serif::network::approx8::Net::ih1], 0.50166260916650918); double energyFraction = netOut.energy / 1.6433051127589775E+18;
EXPECT_DOUBLE_EQ(netOut.composition[serif::network::approx8::Net::ihe4],0.48172270591286032); double H1MassFraction = netOut.composition.getMassFraction("H-1")/ 0.50166262445895604;
EXPECT_DOUBLE_EQ(netOut.energy, 1.6433049870528356e+18); double He4MassFraction = netOut.composition.getMassFraction("He-4") / 0.48172273720971226;
double relError = 1e-8;
EXPECT_NEAR(H1MassFraction, 1.0, relError);
EXPECT_NEAR(He4MassFraction, 1.0, relError);
EXPECT_NEAR(energyFraction, 1.0, relError);
} }

2
tests/network/meson.build
View File

@@ -11,7 +11,7 @@ foreach test_file : test_sources
test_exe = executable( test_exe = executable(
exe_name, exe_name,
test_file, test_file,
dependencies: [gtest_dep, gtest_main, network_dep], dependencies: [gtest_dep, gtest_main, network_dep, species_weight_dep, composition_dep],
include_directories: include_directories('../../src/network/public'), include_directories: include_directories('../../src/network/public'),
link_with: libnetwork, # Link the dobj library link_with: libnetwork, # Link the dobj library
install_rpath: '@loader_path/../../src' # Ensure runtime library path resolves correctly install_rpath: '@loader_path/../../src' # Ensure runtime library path resolves correctly