diff --git a/assets/static/atomic/meson.build b/assets/static/atomic/meson.build index 9e024aab..b5b79b71 100644 --- a/assets/static/atomic/meson.build +++ b/assets/static/atomic/meson.build @@ -1,3 +1,4 @@ species_weight_dep = declare_dependency( include_directories: include_directories('include'), -) \ No newline at end of file +) +message('✅ SERiF species_weight dependency declared') \ No newline at end of file diff --git a/src/network/meson.build b/src/network/meson.build index 41abf3fe..97d410ee 100644 --- a/src/network/meson.build +++ b/src/network/meson.build @@ -15,7 +15,9 @@ dependencies = [ quill_dep, mfem_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 diff --git a/src/network/private/approx8.cpp b/src/network/private/approx8.cpp index c485765e..f46bb7eb 100644 --- a/src/network/private/approx8.cpp +++ b/src/network/private/approx8.cpp @@ -31,7 +31,7 @@ #include "approx8.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 the following 8 isotopes: @@ -78,7 +78,7 @@ namespace serif::network::approx8{ using namespace boost::numeric::odeint; //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){ if (a.size() != b.size()) { 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 vec7 get_T9_array(const double &T) { vec7 arr; - double T9=1e-9*T; - double T913=pow(T9,1./3.); + const double T9=1e-9*T; + const double T913=pow(T9,1./3.); arr[0]=1; 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 double pp_rate(const vec7 &T9) { - 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 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) { - 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 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){ - 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); } // he3(he3,2p)he4 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}; - 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 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 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}; - 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 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){ - 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 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 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}; - 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 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){ - 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}; - 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 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){ - 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}; - 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 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){ - 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}; - 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 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){ - 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}; - 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 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){ - double f1=n15pg_rate(T9); - double f2=n15pa_rate(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){ - 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); } // o16(a,g)ne20 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}; - 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 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){ - 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}; - 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 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){ - 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); } // c12(o16,a)mg24 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); } @@ -244,209 +244,211 @@ namespace serif::network::approx8{ // 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(); const double avo = constants.get("N_a").value; const double clight = constants.get("c").value; // EOS - vec7 T9=get_T9_array(y[Net::itemp]); + const vec7 T9=get_T9_array(y[Approx8Net::iTemp]); // evaluate rates once per call - double rpp=pp_rate(T9); - double r33=he3he3_rate(T9); - double r34=he3he4_rate(T9); - double r3a=triple_alpha_rate(T9); - double rc12p=c12p_rate(T9); - double rc12a=c12a_rate(T9); - double rn14p=n14p_rate(T9); - double rn14a=n14a_rate(T9); - double ro16p=o16p_rate(T9); - double ro16a=o16a_rate(T9); - double rne20a=ne20a_rate(T9); - double r1212=c12c12_rate(T9); - double r1216=c12o16_rate(T9); + 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); - double pfrac=n15pg_frac(T9); - double afrac=1-pfrac; + const double pFrac=n15pg_frac(T9); + const double aFrac=1-pFrac; - double yh1 = y[Net::ih1]; - double yhe3 = y[Net::ihe3]; - double yhe4 = y[Net::ihe4]; - double yc12 = y[Net::ic12]; - double yn14 = y[Net::in14]; - double yo16 = y[Net::io16]; - double yne20 = y[Net::ine20]; + 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 < Net::nvar; i++) { - for (int j=0; j < Net::nvar; j++) { - J(i,j)=0.0; + 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(Net::ih1,Net::ih1) = -3*yh1*rpp - 2*yc12*rc12p -2*yn14*rn14p -2*yo16*ro16p; - J(Net::ih1,Net::ihe3) = 2*yhe3*r33 - yhe4*r34; - J(Net::ih1,Net::ihe4) = -yhe3*r34; - J(Net::ih1,Net::ic12) = -2*yh1*rc12p; - J(Net::ih1,Net::in14) = -2*yh1*rn14p; - J(Net::ih1,Net::io16) = -2*yh1*ro16p; + 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(Net::ihe3,Net::ih1) = yh1*rpp; - J(Net::ihe3,Net::ihe3) = -2*yhe3*r33 - yhe4*r34; - J(Net::ihe3,Net::ihe4) = -yhe3*r34; - + 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(Net::ihe4,Net::ih1) = yn14*afrac*rn14p + yo16*ro16p; - J(Net::ihe4,Net::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(Net::ihe4,Net::ic12) = -yhe4*rc12a + yc12*r1212 + yo16*r1216; - J(Net::ihe4,Net::in14) = yh1*afrac*rn14p - 1.5*yhe4*rn14a; - J(Net::ihe4,Net::io16) = yh1*ro16p - yhe4*ro16a + yc12*r1216; - J(Net::ihe4,Net::ine20) = -yhe4*rne20a; + 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(Net::ic12,Net::ih1) = -yc12*rc12p + yn14*afrac*rn14p; - J(Net::ic12,Net::ihe4) = 0.5*yhe4*yhe4*r3a - yhe4*rc12a; - J(Net::ic12,Net::ic12) = -yh1*rc12p - yhe4*rc12a - yo16*r1216 - 2*yc12*r1212; - J(Net::ic12,Net::in14) = yh1*yn14*afrac*rn14p; - J(Net::ic12,Net::io16) = -yc12*r1216; + 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(Net::in14,Net::ih1) = yc12*rc12p - yn14*rn14p + yo16*ro16p; - J(Net::in14,Net::ihe4) = -yn14*rn14a; - J(Net::in14,Net::ic12) = yh1*rc12p; - J(Net::in14,Net::in14) = -yh1*rn14p - yhe4*rn14a; - J(Net::in14,Net::io16) = yo16*ro16p; - + 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(Net::io16,Net::ih1) = yn14*pfrac*rn14p - yo16*ro16p; - J(Net::io16,Net::ihe4) = yc12*rc12a - yo16*ro16a; - J(Net::io16,Net::ic12) = yhe4*rc12a - yo16*r1216; - J(Net::io16,Net::in14) = yh1*pfrac*rn14p; - J(Net::io16,Net::io16) = yh1*ro16p - yc12*r1216 -yhe4*ro16a; + 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(Net::ine20,Net::ihe4) = yn14*rn14a + yo16*ro16a - yne20*rne20a; - J(Net::ine20,Net::ic12) = yc12*r1212; - J(Net::ine20,Net::in14) = yhe4*rn14a; - J(Net::ine20,Net::io16) = yo16*ro16a; - J(Net::ine20,Net::ine20) = -yhe4*rne20a; + 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(Net::img24,Net::ihe4) = yne20*rne20a; - J(Net::img24,Net::ic12) = yo16*r1216; - J(Net::img24,Net::io16) = yc12*r1216; - J(Net::img24,Net::ine20) = yhe4*rne20a; + 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 accountings - for (int j=0; j("Network:Approx8:Stiff:AbsTol", 1.0e-6); @@ -463,7 +465,7 @@ namespace serif::network::approx8{ std::make_pair(ODE(), Jacobian()), m_y, 0.0, - m_tmax, + m_tMax, m_dt0 ); @@ -474,31 +476,33 @@ namespace serif::network::approx8{ ODE(), m_y, 0.0, - m_tmax, + m_tMax, m_dt0 ); } - double ysum = 0.0; - for (int i = 0; i < Net::niso; i++) { - m_y[i] *= Net::aion[i]; - ysum += m_y[i]; + 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 < Net::niso; i++) { - m_y[i] /= ysum; + for (int i = 0; i < Approx8Net::nIso; i++) { + m_y[i] /= ySum; } NetOut netOut; std::vector 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]); } - netOut.energy = m_y[Net::iener]; - netOut.composition = outComposition; + netOut.energy = m_y[Approx8Net::iEnergy]; netOut.num_steps = num_steps; + const std::vector 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; } @@ -507,31 +511,26 @@ namespace serif::network::approx8{ } vector_type Approx8Network::convert_netIn(const NetIn &netIn) { - if (netIn.composition.size() != Net::niso) { - LOG_ERROR(m_logger, "Error: composition size mismatch in convert_netIn"); - throw std::runtime_error("Error: composition size mismatch in convert_netIn"); - } + vector_type y(Approx8Net::nVar, 0.0); + y[Approx8Net::ih1] = netIn.composition.getMassFraction("H-1"); + 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); - y[Net::ih1] = netIn.composition[0]; - y[Net::ihe3] = netIn.composition[1]; - y[Net::ihe4] = netIn.composition[2]; - 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]; + double ySum = 0.0; + for (int i = 0; i < Approx8Net::nIso; i++) { + y[i] /= Approx8Net::aIon[i]; + ySum += y[i]; } - for (int i = 0; i < Net::niso; i++) { - y[i] /= ysum; + for (int i = 0; i < Approx8Net::nIso; i++) { + y[i] /= ySum; } return y; diff --git a/src/network/private/network.cpp b/src/network/private/network.cpp index d3b3b25e..b0e3e2d7 100644 --- a/src/network/private/network.cpp +++ b/src/network/private/network.cpp @@ -19,18 +19,50 @@ // // *********************************************************************** */ #include "network.h" + +#include "approx8.h" #include "probe.h" #include "quill/LogMacros.h" namespace serif::network { - Network::Network() : + Network::Network(const NetworkFormat format) : m_config(serif::config::Config::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) { - // You can throw an exception here or log a warning if it should never be used. - LOG_ERROR(m_logger, "nuclearNetwork::Network::evaluate() is not implemented"); - throw std::runtime_error("nuclearNetwork::Network::evaluate() is not implemented"); + NetOut netOut; + switch (m_format) { + 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; } } diff --git a/src/network/public/approx8.h b/src/network/public/approx8.h index 711d2a1e..634225e0 100644 --- a/src/network/public/approx8.h +++ b/src/network/public/approx8.h @@ -31,54 +31,55 @@ * @brief Header file 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. */ -/** - * @typedef vector_type - * @brief Alias for a vector of doubles using Boost uBLAS. - */ -typedef boost::numeric::ublas::vector< double > vector_type; - -/** - * @typedef matrix_type - * @brief Alias for a matrix of doubles using Boost uBLAS. - */ -typedef boost::numeric::ublas::matrix< double > matrix_type; - -/** - * @typedef vec7 - * @brief Alias for a std::array of 7 doubles. - */ -typedef std::array vec7; namespace serif::network::approx8{ + /** + * @typedef vector_type + * @brief Alias for a vector of doubles using Boost uBLAS. + */ + typedef boost::numeric::ublas::vector< double > vector_type; + + /** + * @typedef matrix_type + * @brief Alias for a matrix of doubles using Boost uBLAS. + */ + typedef boost::numeric::ublas::matrix< double > matrix_type; + + /** + * @typedef vec7 + * @brief Alias for a std::array of 7 doubles. + */ + typedef std::array vec7; + using namespace boost::numeric::odeint; /** - * @struct Net + * @struct Approx8Net * @brief Contains constants and arrays related to the nuclear network. */ - struct Net{ - const static int ih1=0; - const static int ihe3=1; - const static int ihe4=2; - const static int ic12=3; - const static int in14=4; - const static int io16=5; - const static int ine20=6; - const static int img24=7; + struct Approx8Net{ + static constexpr int ih1=0; + static constexpr int ihe3=1; + static constexpr int ihe4=2; + static constexpr int ic12=3; + static constexpr int in14=4; + static constexpr int io16=5; + static constexpr int ine20=6; + static constexpr int img24=7; - const static int itemp=img24+1; - const static int iden =itemp+1; - const static int iener=iden+1; + static constexpr int iTemp=img24+1; + static constexpr int iDensity =iTemp+1; + static constexpr int iEnergy=iDensity+1; - const static int niso=img24+1; // number of isotopes - const static int nvar=iener+1; // number of variables + static constexpr int nIso=img24+1; // number of isotopes + static constexpr int nVar=iEnergy+1; // number of variables - static constexpr std::array aion = { + static constexpr std::array aIon = { 1, 3, 4, @@ -89,7 +90,7 @@ namespace serif::network::approx8{ 24 }; - static constexpr std::array mion = { + static constexpr std::array mIon = { 1.67262164e-24, 5.00641157e-24, 6.64465545e-24, @@ -270,10 +271,9 @@ namespace serif::network::approx8{ * @brief Calculates the Jacobian matrix. * @param y State vector. * @param J Jacobian matrix. - * @param t Time. * @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. * @param y 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 * @brief Class for the Approx8 nuclear reaction network. */ - class Approx8Network : public Network { + class Approx8Network final : public Network { public: + Approx8Network(); + /** * @brief Evaluates the nuclear network. * @param netIn Input parameters for 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. @@ -313,11 +314,11 @@ namespace serif::network::approx8{ * @brief Checks if the solver is using a stiff method. * @return Boolean indicating if a stiff method is being used. */ - bool isStiff() { return m_stiff; } + bool isStiff() const { return m_stiff; } private: vector_type m_y; - double m_tmax; - double m_dt0; + double m_tMax = 0; + double m_dt0 = 0; bool m_stiff = false; /** @@ -325,7 +326,8 @@ namespace serif::network::approx8{ * @param netIn Input parameters for the network. * @return Internal state vector. */ - vector_type convert_netIn(const NetIn &netIn); + static vector_type convert_netIn(const NetIn &netIn); }; + } // namespace nnApprox8 diff --git a/src/network/public/network.h b/src/network/public/network.h index 0ae46779..0d7692be 100644 --- a/src/network/public/network.h +++ b/src/network/public/network.h @@ -25,15 +25,27 @@ #include "probe.h" #include "config.h" #include "quill/Logger.h" +#include "composition.h" +#include namespace serif::network { + enum NetworkFormat { + APPROX8, ///< Approx8 nuclear reaction network format. + UNKNOWN, + }; + + static inline std::unordered_map FormatStringLookup = { + {APPROX8, "Approx8"}, + {UNKNOWN, "Unknown"} + }; + /** * @struct NetIn * @brief Input structure for the network evaluation. - * + * * This structure holds the input parameters required for the network evaluation. - * + * * Example usage: * @code * nuclearNetwork::NetIn netIn; @@ -46,8 +58,8 @@ namespace serif::network { * @endcode */ struct NetIn { - std::vector composition; ///< Composition of the network - double tmax; ///< Maximum time + serif::composition::Composition composition; ///< Composition of the network + double tMax; ///< Maximum time double dt0; ///< Initial time step double temperature; ///< Temperature in Kelvin double density; ///< Density in g/cm^3 @@ -69,7 +81,7 @@ namespace serif::network { * @endcode */ struct NetOut { - std::vector 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 double energy; ///< Energy in ergs after evaluation }; @@ -90,9 +102,12 @@ namespace serif::network { */ class Network { public: - Network(); + explicit Network(const NetworkFormat format = NetworkFormat::APPROX8); virtual ~Network() = default; + NetworkFormat getFormat() const; + NetworkFormat setFormat(const NetworkFormat format); + /** * @brief Evaluate the network based on the input parameters. * @@ -105,6 +120,10 @@ namespace serif::network { serif::config::Config& m_config; ///< Configuration instance serif::probe::LogManager& m_logManager; ///< Log manager instance quill::Logger* m_logger; ///< Logger instance + + NetworkFormat m_format; ///< Format of the network }; + + } // namespace nuclearNetwork diff --git a/tests/network/approx8Test.cpp b/tests/network/approx8Test.cpp index 2bfe7fa2..899eaf68 100644 --- a/tests/network/approx8Test.cpp +++ b/tests/network/approx8Test.cpp @@ -4,6 +4,7 @@ #include "approx8.h" #include "config.h" #include "network.h" +#include "composition.h" #include @@ -32,19 +33,31 @@ TEST_F(approx8Test, evaluate) { serif::network::NetIn netIn; std::vector comp = {0.708, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4}; + std::vector 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.density = 1e2; netIn.energy = 0.0; - netIn.tmax = 3.15e17; + netIn.tMax = 3.15e17; netIn.dt0 = 1e12; serif::network::NetOut netOut; EXPECT_NO_THROW(netOut = network.evaluate(netIn)); - EXPECT_DOUBLE_EQ(netOut.composition[serif::network::approx8::Net::ih1], 0.50166260916650918); - EXPECT_DOUBLE_EQ(netOut.composition[serif::network::approx8::Net::ihe4],0.48172270591286032); - EXPECT_DOUBLE_EQ(netOut.energy, 1.6433049870528356e+18); + double energyFraction = netOut.energy / 1.6433051127589775E+18; + double H1MassFraction = netOut.composition.getMassFraction("H-1")/ 0.50166262445895604; + 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); } diff --git a/tests/network/meson.build b/tests/network/meson.build index 7b840631..1e6ccf3e 100644 --- a/tests/network/meson.build +++ b/tests/network/meson.build @@ -11,7 +11,7 @@ foreach test_file : test_sources test_exe = executable( exe_name, 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'), link_with: libnetwork, # Link the dobj library install_rpath: '@loader_path/../../src' # Ensure runtime library path resolves correctly