From 03c7d428a71e00744dd5682ba3f23aa5eccd784c Mon Sep 17 00:00:00 2001 From: Aaron Dotter Date: Thu, 20 Mar 2025 15:22:04 -0400 Subject: [PATCH 01/18] adding a nuclear reaction network --- src/network/private/network.cpp | 533 ++++++++++++++++++++++++++++++++ 1 file changed, 533 insertions(+) create mode 100644 src/network/private/network.cpp diff --git a/src/network/private/network.cpp b/src/network/private/network.cpp new file mode 100644 index 0000000..49a8692 --- /dev/null +++ b/src/network/private/network.cpp @@ -0,0 +1,533 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include "const.h" +//#include "probe.h" +//#include "config.h" +//#include "quill/LogMacros.h" + +/* Nuclear reaction network in cgs units based on Frank Timmes' "aprox8". + 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 . + +*/ + + +// using namespace std; +using namespace boost::numeric::odeint; +namespace phoenix = boost::phoenix; + +// these types are required by the Rosenbrock implicit solver +typedef boost::numeric::ublas::vector< double > vector_type; +typedef boost::numeric::ublas::matrix< double > matrix_type; + +// this array is used only in the nuclear reaction rate evaluations +typedef std::array vec; + +// only need a couple of constants +Constants& constants = Constants::getInstance(); +const double avo = constants.get("N_a").value; +const double clight = constants.get("c").value; + +// identify the isotopes used in the network. +const int ih1=0; +const int ihe3=1; +const int ihe4=2; +const int ic12=3; +const int in14=4; +const int io16=5; +const int ine20=6; +const int img24=7; + +// physical variables; this routine currently does not need to call EOS +// since the temperature and density are assumed constant during the burn +const int itemp=img24+1; +const int iden =itemp+1; +const int iener=iden+1; + +// number of isotopes and number of variables +const int niso=img24+1; // number of isotopes +const int nvar=iener+1; // number of variables + +// atomic stuff +std::array aion = {1,3,4,12,14,16,20,24}; +//std::array zion = {1,2,2, 6, 7, 8,10,12}; +//std::array bion = {0,7.71819,28.29603,92.16294,104.65998,127.62093,160.64788,198.25790}; +//nion = aion - zion #neutrons +//mion = nion*mn + zion*mp - bion*mev2gr #mass +std::array mion = {1.67262164e-24, 5.00641157e-24, 6.64465545e-24, 1.99209977e-23, + 2.32462686e-23, 2.65528858e-23, 3.31891077e-23, 3.98171594e-23}; + +//helper functions +// a function to multilpy two arrays and then sum the resulting elements: sum(a*b) +double sum_product( const vec &a, const vec &b){ + if (a.size() != b.size()){ + throw std::runtime_error("Error: array size mismatch in sum_product"); + } + + 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 +vec get_T9_array(const double &T) { + vec arr; + double T9=1e-9*T; + 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 vec &T9, const vec &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 vec &T9) { + vec a1 = {-34.78630, 0,-3.511930, 3.100860, -0.1983140, 1.262510e-2, -1.025170}; + vec 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 vec &T9) { + vec a1 = {7.528980, 0, -3.720800, 0.8717820, 0, 0,-0.6666670}; + vec 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 vec &T9){ + vec 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 vec &T9){ + vec a1 = {1.560990e+01,0.000000e+00,-1.282710e+01,-3.082250e-02,-6.546850e-01,8.963310e-02,-6.666670e-01}; + vec 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 vec &T9){ + vec a1 = {-9.710520e-01,0.000000e+00,-3.706000e+01,2.934930e+01,-1.155070e+02,-1.000000e+01,-1.333330e+00}; + vec a2 = {-1.178840e+01,-1.024460e+00,-2.357000e+01,2.048860e+01,-1.298820e+01,-2.000000e+01,-2.166670e+00}; + vec 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 vec &T9){ + vec a1={1.714820e+01,0.000000e+00,-1.369200e+01,-2.308810e-01,4.443620e+00,-3.158980e+00,-6.666670e-01}; + vec 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 vec &T9){ + vec a1={6.965260e+01,-1.392540e+00,5.891280e+01,-1.482730e+02,9.083240e+00,-5.410410e-01,7.035540e+01}; + vec 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 vec &T9){ + vec a1={1.701000e+01,0.000000e+00,-1.519300e+01,-1.619540e-01,-7.521230e+00,-9.875650e-01,-6.666670e-01}; + vec a2={2.011690e+01,0.000000e+00,-1.519300e+01,-4.639750e+00,9.734580e+00,-9.550510e+00,3.333330e-01}; + vec a3={7.654440e+00,-2.998000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + vec 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 vec &T9){ + vec a1={2.153390e+01,0.000000e+00,-3.625040e+01,0.000000e+00,0.000000e+00,-5.000000e+00,-6.666670e-01}; + vec a2={1.968380e-01,-5.160340e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + vec 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 vec &T9){ + vec a1 = {2.747640e+01,0.000000e+00,-1.525300e+01,1.593180e+00,2.447900e+00,-2.197080e+00,-6.666670e-01}; + vec a2 = {-4.873470e+00,-2.021170e+00,0.000000e+00,3.084970e+01,-8.504330e+00,-1.544260e+00,-1.500000e+00}; + vec a3 = {2.089720e+01,-7.406000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + vec 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 vec &T9){ + vec a1 = {2.001760e+01,0.000000e+00,-1.524000e+01,3.349260e-01,4.590880e+00,-4.784680e+00,-6.666670e-01}; + vec a2 = {6.590560e+00,-2.923150e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + vec 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 vec &T9){ + double f1=n15pg_rate(T9); + double f2=n15pa_rate(T9); + return f1/(f1+f2); +} + +// o16(p,g)f17 then f17 -> o17(p,a)n14 +double o16p_rate(const vec &T9){ + vec 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 vec &T9){ + vec a1={2.390300e+01,0.000000e+00,-3.972620e+01,-2.107990e-01,4.428790e-01,-7.977530e-02,-6.666670e-01}; + vec a2={3.885710e+00,-1.035850e+01,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + vec 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 vec &T9){ + vec a1={2.450580e+01,0.000000e+00,-4.625250e+01,5.589010e+00,7.618430e+00,-3.683000e+00,-6.666670e-01}; + vec a2={-3.870550e+01,-2.506050e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + vec a3={1.983070e+00,-9.220260e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + vec 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 vec &T9){ + vec 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 vec &T9){ + vec 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 +struct Jacobian { + + void operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) + { + // EOS + vec T9=get_T9_array(y[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); + + double pfrac=n15pg_frac(T9); + double afrac=1-pfrac; + + double yh1 = y[ ih1]; + double yhe3 = y[ ihe3]; + double yhe4 = y[ ihe4]; + double yc12 = y[ ic12]; + double yn14 = y[ in14]; + double yo16 = y[ io16]; + double yne20 = y[ine20]; + + // zero all elements to begin + for (int i=0; i < nvar; i++) { + for (int j=0; j < nvar; j++) { J(i,j)=0.0; }} + + // h1 jacobian elements + J(ih1,ih1) = -3*yh1*rpp - 2*yc12*rc12p -2*yn14*rn14p -2*yo16*ro16p; + J(ih1,ihe3) = 2*yhe3*r33 - yhe4*r34; + J(ih1,ihe4) = -yhe3*r34; + J(ih1,ic12) = -2*yh1*rc12p; + J(ih1,in14) = -2*yh1*rn14p; + J(ih1,io16) = -2*yh1*ro16p; + + // he3 jacobian elements + J(ihe3,ih1) = yh1*rpp; + J(ihe3,ihe3) = -2*yhe3*r33 - yhe4*r34; + J(ihe3,ihe4) = -yhe3*r34; + + // he4 jacobian elements + J(ihe4,ih1) = yn14*afrac*rn14p + yo16*ro16p; + J(ihe4,ihe3) = yhe3*r33 - yhe4*r34; + J(ihe4,ihe4) = yhe3*r34 - 1.5*yhe4*yhe4*r3a - yc12*rc12a - 1.5*yn14*rn14a - yo16*ro16a - yne20*rne20a; + J(ihe4,ic12) = -yhe4*rc12a + yc12*r1212 + yo16*r1216; + J(ihe4,in14) = yh1*afrac*rn14p - 1.5*yhe4*rn14a; + J(ihe4,io16) = yh1*ro16p - yhe4*ro16a + yc12*r1216; + J(ihe4,ine20) = -yhe4*rne20a; + + // c12 jacobian elements + J(ic12,ih1) = -yc12*rc12p + yn14*afrac*rn14p; + J(ic12,ihe4) = 0.5*yhe4*yhe4*r3a - yhe4*rc12a; + J(ic12,ic12) = -yh1*rc12p - yhe4*rc12a - yo16*r1216 - 2*yc12*r1212; + J(ic12,in14) = yh1*yn14*afrac*rn14p; + J(ic12,io16) = -yc12*r1216; + + // n14 jacobian elements + J(in14,ih1) = yc12*rc12p - yn14*rn14p + yo16*ro16p; + J(in14,ihe4) = -yn14*rn14a; + J(in14,ic12) = yh1*rc12p; + J(in14,in14) = -yh1*rn14p - yhe4*rn14a; + J(in14,io16) = yo16*ro16p; + + // o16 jacobian elements + J(io16,ih1) = yn14*pfrac*rn14p - yo16*ro16p; + J(io16,ihe4) = yc12*rc12a - yo16*ro16a; + J(io16,ic12) = yhe4*rc12a - yo16*r1216; + J(io16,in14) = yh1*pfrac*rn14p; + J(io16,io16) = yh1*ro16p - yc12*r1216 -yhe4*ro16a; + + // ne20 jacobian elements + J(ine20,ihe4) = yn14*rn14a + yo16*ro16a - yne20*rne20a; + J(ine20,ic12) = yc12*r1212; + J(ine20,in14) = yhe4*rn14a; + J(ine20,io16) = yo16*ro16a; + J(ine20,ine20) = -yhe4*rne20a; + + // mg24 jacobian elements + J(img24,ihe4) = yne20*rne20a; + J(img24,ic12) = yo16*r1216; + J(img24,io16) = yc12*r1216; + J(img24,ine20) = yhe4*rne20a; + + // energy accounting + for (int j=0; j>(1.0e-6,1.0e-6) , + std::make_pair( ODE(), Jacobian() ), y, 0.0, tmax, dt0); + + } else { + std::cout << " *** Explicit RK Dormand-Prince *** " << std::endl; + + num_of_steps = integrate_const( make_dense_output>(1.0e-6, 1.0e-6), + ODE(), y, 0.0, tmax, dt0); + + } + + //convert number fraction to mass fraction + ysum=0; + for (int i=0; i < niso; i++) { + y[i] *= aion[i]; + ysum+= y[i]; + } + for (int i=0; i < niso; i++) {y[i] /= ysum;} + + std::cout << " H1: " << y[ih1] << std::endl; + std::cout << " He4: " << y[ihe4] << std::endl; + std::cout << "energy: " << y[iener] << std::endl; + std::cout << " steps: " << num_of_steps << std::endl; + return 0; + +} From 5d1044a55afcddbe14d0c151cc584ecf8b004b08 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 21 Mar 2025 10:39:34 -0400 Subject: [PATCH 02/18] build(boost): incorporated boost into the build system all builds should be run with ./mk now --- .gitignore | 3 + build-config/boost/install.sh | 226 +++++++++++++++++++++++++++++++++ build-config/boost/meson.build | 1 + build-config/meson.build | 3 +- mk | 174 ++++++++++++++++++++++--- 5 files changed, 390 insertions(+), 17 deletions(-) create mode 100755 build-config/boost/install.sh create mode 100644 build-config/boost/meson.build diff --git a/.gitignore b/.gitignore index 90155ef..4e708e7 100644 --- a/.gitignore +++ b/.gitignore @@ -69,3 +69,6 @@ subprojects/googletest-1.15.2/ *.log output/ + +.boost_installed +4DSSE_logs/ diff --git a/build-config/boost/install.sh b/build-config/boost/install.sh new file mode 100755 index 0000000..ad30e59 --- /dev/null +++ b/build-config/boost/install.sh @@ -0,0 +1,226 @@ +#!/usr/bin/env bash +# +# install.sh - Interactive installation script for Boost. +# +# This script checks if Boost is installed (by looking for boost/version.hpp in common locations). +# If not found, it prompts the user to install Boost using the native package manager for: +# - FreeBSD, Ubuntu, Debian, Fedora, Arch, Mint, Manjaro, macOS, OpenSuse, and Nix. +# +# All output is colorized for clarity and logged to meson-log.txt. + +set -e + +# ANSI color codes. +RED="\033[0;31m" +GREEN="\033[0;32m" +YELLOW="\033[0;33m" +BLUE="\033[0;34m" +NC="\033[0m" # No Color + +# Log file. +LOGFILE="4DSSE-install-log.txt" +touch "$LOGFILE" + + +# Log function: prints to stdout and appends to logfile. +log() { + local message="$1" + # Print the colored message to stdout. + echo -e "$message" + # Strip ANSI escape sequences and append the cleaned message to the log file. + echo -e "$message" | sed -r 's/\x1B\[[0-9;]*[mK]//g' >> "$LOGFILE" +} + + + +check_boost_installed() { + log "${BLUE}[Info] Checking for Boost::numeric and Boost::phoenix using Meson...${NC}" + local tmpDir + if [[ -d mesonTest ]]; then + rm -rf mesonTest + fi + tmpDir=$(mkdir mesonTest) + local sourceFile="mesonTest/test.cpp" + local mesonFile="mesonTest/meson.build" + + # Write test.cpp + cat > "$sourceFile" < +#include + +int main() { + boost::numeric::ublas::vector v(3); + v[0] = 1.0; + return v[0]; +} +EOF + + # Write meson.build + cat > "$mesonFile" < /dev/null; then + log "${YELLOW}[Warning] Homebrew is not installed. It is recommended for installing Boost on macOS.${NC}" + if prompt_yes_no "[Query] Would you like to install Homebrew now? (y/n): "; then + log "${BLUE}[Info] Installing Homebrew...${NC}" + /bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)" | tee -a "$LOGFILE" + else + log "${RED}[Error] Homebrew is required for an easy Boost installation on macOS. Aborting.${NC}" + exit 1 + fi + fi + log "${BLUE}[Info] Running: brew install boost${NC}" + brew install boost | tee -a "$LOGFILE" + ;; + "Linux") + case "$DISTRO" in + "ubuntu"|"debian"|"linuxmint") + log "${BLUE}[Info] Running: sudo apt-get update && sudo apt-get install -y libboost-all-dev${NC}" + sudo apt-get update | tee -a "$LOGFILE" + sudo apt-get install -y libboost-all-dev | tee -a "$LOGFILE" + ;; + "fedora") + log "${BLUE}[Info] Running: sudo dnf install boost-devel${NC}" + sudo dnf install -y boost-devel | tee -a "$LOGFILE" + ;; + "arch"|"manjaro") + log "${BLUE}[Info] Running: sudo pacman -S boost${NC}" + sudo pacman -S --noconfirm boost | tee -a "$LOGFILE" + ;; + "opensuse") + log "${BLUE}[Info] Running: sudo zypper install libboost-devel${NC}" + sudo zypper install -y libboost-devel | tee -a "$LOGFILE" + ;; + "nixos"|"nix") + log "${BLUE}[Info] Running: nix-env -iA nixpkgs.boost${NC}" + nix-env -iA nixpkgs.boost | tee -a "$LOGFILE" + ;; + *) + log "${RED}[Error] Your Linux distribution is not recognized. Please install Boost manually.${NC}" + exit 1 + ;; + esac + ;; + "FreeBSD") + log "${BLUE}[Info] Running: sudo pkg install boost-all${NC}" + sudo pkg install -y boost-all | tee -a "$LOGFILE" + ;; + *) + log "${RED}[Error] Automatic Boost installation is not supported on OS: ${OS}. Please install Boost manually.${NC}" + exit 1 + ;; + esac + + # Verify Boost installation. + if check_boost_installed; then + log "${GREEN}[Success] Boost installation succeeded.${NC}" + else + log "${RED}[Error] Boost installation appears to have failed. Please install Boost manually.${NC}" + exit 1 + fi + + else + log "${RED}[Error] Boost is required. Please install it using the appropriate command for your system:${NC}" + case "$OS" in + "macOS") + log "${RED}[INFO] brew install boost (Install Homebrew from https://brew.sh if not present)${NC}" + ;; + "Linux") + case "$DISTRO" in + "ubuntu"|"debian"|"linuxmint") + log "${RED}[INFO] sudo apt-get install libboost-all-dev${NC}" + ;; + "fedora") + log "${RED}[INFO] sudo dnf install boost-devel${NC}" + ;; + "arch"|"manjaro") + log "${RED}[INFO] sudo pacman -S boost${NC}" + ;; + "opensuse") + log "${RED}[INFO] sudo zypper install libboost-devel${NC}" + ;; + "nixos"|"nix") + log "${RED}[INFO] nix-env -iA nixpkgs.boost${NC}" + ;; + *) + log "${RED}[INFO] Please consult your distribution's documentation for installing Boost.${NC}" + ;; + esac + ;; + "FreeBSD") + log "${RED}[INFO] sudo pkg install boost-all${NC}" + ;; + *) + log "${RED}[INFO] Please consult your operating system's documentation for installing Boost.${NC}" + ;; + esac + exit 1 + fi +fi + +check_boost_installed + +log "${GREEN}[Success] Boost installed successfully!${NC}" +log "${GREEN}[Success] 4DSSE installation can now continue.${NC}" diff --git a/build-config/boost/meson.build b/build-config/boost/meson.build new file mode 100644 index 0000000..fa83c8f --- /dev/null +++ b/build-config/boost/meson.build @@ -0,0 +1 @@ +boost_dep = dependency('boost', version: '1.87.0') diff --git a/build-config/meson.build b/build-config/meson.build index 92335cf..d800bdd 100644 --- a/build-config/meson.build +++ b/build-config/meson.build @@ -2,4 +2,5 @@ cmake = import('cmake') subdir('mfem') subdir('yaml-cpp') -subdir('quill') \ No newline at end of file +subdir('quill') +subdir('boost') \ No newline at end of file diff --git a/mk b/mk index 5cc0ba7..1c9bb3c 100755 --- a/mk +++ b/mk @@ -1,23 +1,165 @@ -#!/bin/bash +#!/usr/bin/env bash -# Check for the --noTest flag -if [[ "$1" == "--noTest" ]]; then - meson setup build -Dbuild_tests=false --buildtype=release +set -e + +# --- Defaults --- +userFlag=0 +testsFlag=0 +buildDir="build" +installDir="$HOME/.local/4DSSE" + +# --- Help Menu --- +show_help() { + cat <> "$LOGFILE" +} + +if [[ -f "$LOGFILE" ]]; then + # Get the calendar datetime of the log file creation. + # convert UNIX timestamp to human-readable date + rm "$LOGFILE" + log "${BLUE}[Info] Old log file removed.${NC}" +fi +touch "$LOGFILE" + +# --- Build 4DSSE --- +log "${BLUE}[Info] Building 4DSSE...${NC}" +if [[ $userFlag -eq 1 ]]; then + log "${BLUE}[Info] Installing 4DSSE in user mode...${NC}" +else + log "${BLUE}[Info] Installing 4DSSE for developers...${NC}" fi -# Compile the project -meson compile -C build - -# If tests are built, run them -if [[ "$1" != "--noTest" ]]; then - meson test -C build +if [[ $testsFlag -eq 1 ]]; then + log "${BLUE}[Info] Building with tests...${NC}" +else + log "${BLUE}[Info] Building without tests...${NC}" fi -# Check if --docs are to be built -if [[ "$*" == *"--docs"* ]]; then - echo "Generating documentation..." - doxygen - cd docs/latex && make +# --- First check Boost status --- +log "${BLUE}[Info] Checking Boost status...${NC}" +# if the following script exists with anything other than a 0 status the script will exit +if [[ -f ./build-config/boost/.boost_installed ]]; then + log "${BLUE}[Info] Boost already installed. Skipping...${NC}" +else + log "${BLUE}[Info] Installing Boost...${NC}" + if ! ./build-config/boost/install.sh; then + log "${RED}[Error] Boost check failed. Exiting...${NC}" + exit 1 + else + touch ./build-config/boost/.boost_installed + log "${GREEN}[Success] Boost check passed.${NC}" + fi fi + +# --- Check if the build dir has been configured with the correct flags --- +log "${BLUE}[Info] Configuring build directory...${NC}" +if [ -f "$buildDir/build.ninja" ]; then + log "${BLUE}[Info] Build directory already configured. Skipping...${NC}" +else + if [[ $userFlag -eq 1 ]]; then + meson setup "$buildDir" -Duser_mode=true --buildtype=release + else + if [[ $testsFlag -eq 1 ]]; then + meson setup "$buildDir" --buildtype=debug + else + meson setup "$buildDir" --buildtype=debug -Dbuild_tests=false + fi + fi + log "${GREEN}[Success] Build directory configured successfully.${NC}" +fi + +log "${BLUE}[Info] Building 4DSSE...${NC}" +meson compile -C "$buildDir" + +# --- Install 4DSSE --- +if [[ $userFlag -eq 1 ]]; then + log "${BLUE}[Info] Installing 4DSSE in user mode...${NC}" + meson install -C "$buildDir" --no-rebuild --prefix="$installDir" +fi + +if [[ $testsFlag -eq 1 ]]; then + log "${BLUE}[Info] Running tests...${NC}" + meson test -C "$buildDir" +fi + +if [[ $userFlag -eq 1 ]]; then + log "${GREEN}[Success] 4DSSE built and installed successfully.${NC}" +else + log "${GREEN}[Success] 4DSSE built successfully.${NC}" +fi \ No newline at end of file From 6876b87947880eff8be480e438ebc58d5befdaaa Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 21 Mar 2025 10:39:53 -0400 Subject: [PATCH 03/18] feat(network): began incoporating network --- src/meson.build | 5 ++++- src/network/meson.build | 30 ++++++++++++++++++++++++++++++ src/network/private/network.cpp | 5 ----- src/network/public/network.h | 0 4 files changed, 34 insertions(+), 6 deletions(-) create mode 100644 src/network/meson.build create mode 100644 src/network/public/network.h diff --git a/src/meson.build b/src/meson.build index d2d2399..54fbee7 100644 --- a/src/meson.build +++ b/src/meson.build @@ -1,4 +1,6 @@ # Build the main source code in the correct order +# Unless you know what you are doing, do not change the order of the subdirectories +# as there are dependencies which exist between them. # Utility Libraries subdir('misc') @@ -15,4 +17,5 @@ subdir('meshIO') # Resouce Manager Libraries subdir('resource') -# Physics Libraries \ No newline at end of file +# Physics Libraries +subdir('network') \ No newline at end of file diff --git a/src/network/meson.build b/src/network/meson.build new file mode 100644 index 0000000..72f5bd1 --- /dev/null +++ b/src/network/meson.build @@ -0,0 +1,30 @@ +# Define the library +network_sources = files( + 'private/network.cpp', +) + +network_headers = files( + 'public/network.h' +) + +dependencies = [ + boost_dep, + const_dep +] + +# Define the libnetwork library so it can be linked against by other parts of the build system +libnetwork = static_library('network', + network_sources, + include_directories: include_directories('public'), + dependencies: dependencies, + install : true) + +network_dep = declare_dependency( + include_directories: include_directories('public'), + link_with: libnetwork, + sources: network_sources, + dependencies: dependencies, +) + +# Make headers accessible +install_headers(network_headers, subdir : '4DSSE/network') \ No newline at end of file diff --git a/src/network/private/network.cpp b/src/network/private/network.cpp index 49a8692..5a9ff2b 100644 --- a/src/network/private/network.cpp +++ b/src/network/private/network.cpp @@ -1,12 +1,7 @@ #include -#include -#include #include -#include #include #include -#include -#include #include #include diff --git a/src/network/public/network.h b/src/network/public/network.h new file mode 100644 index 0000000..e69de29 From 43f57d5e738ee54b5d2b14061483f8cf80d34ab1 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 21 Mar 2025 11:14:28 -0400 Subject: [PATCH 04/18] build(mk): updated auto build script --- .gitignore | 1 + mk | 104 ++++++++++++++++++++++++++++++++++++++++++----------- 2 files changed, 84 insertions(+), 21 deletions(-) diff --git a/.gitignore b/.gitignore index 4e708e7..77b7f3d 100644 --- a/.gitignore +++ b/.gitignore @@ -72,3 +72,4 @@ output/ .boost_installed 4DSSE_logs/ +.last_build_flags diff --git a/mk b/mk index 1c9bb3c..28e2811 100755 --- a/mk +++ b/mk @@ -5,9 +5,11 @@ set -e # --- Defaults --- userFlag=0 testsFlag=0 +runTestsFlag=0 buildDir="build" installDir="$HOME/.local/4DSSE" + # --- Help Menu --- show_help() { cat < /dev/null; then + log "${RED}[Error] Meson is not installed. Exiting...${NC}" + log "${YELLOW}[Info] Please install Meson and try again.${NC}" + exit 1 +else + log "${MAGENTA}[Success] Meson is installed. Continuing...${NC}" +fi + # --- Build 4DSSE --- log "${BLUE}[Info] Building 4DSSE...${NC}" if [[ $userFlag -eq 1 ]]; then @@ -105,7 +123,7 @@ else log "${BLUE}[Info] Installing 4DSSE for developers...${NC}" fi -if [[ $testsFlag -eq 1 ]]; then +if [[ $testsFlag -eq 0 ]]; then log "${BLUE}[Info] Building with tests...${NC}" else log "${BLUE}[Info] Building without tests...${NC}" @@ -115,7 +133,7 @@ fi log "${BLUE}[Info] Checking Boost status...${NC}" # if the following script exists with anything other than a 0 status the script will exit if [[ -f ./build-config/boost/.boost_installed ]]; then - log "${BLUE}[Info] Boost already installed. Skipping...${NC}" + log "${MAGENTA}[Succsess] Boost already installed. Skipping...${NC}" else log "${BLUE}[Info] Installing Boost...${NC}" if ! ./build-config/boost/install.sh; then @@ -127,21 +145,62 @@ else fi fi -# --- Check if the build dir has been configured with the correct flags --- -log "${BLUE}[Info] Configuring build directory...${NC}" -if [ -f "$buildDir/build.ninja" ]; then - log "${BLUE}[Info] Build directory already configured. Skipping...${NC}" -else - if [[ $userFlag -eq 1 ]]; then - meson setup "$buildDir" -Duser_mode=true --buildtype=release +#check if the last build flags are the same as the current build flags +# if the flags are the same skip the configuration step +# if they are different then reconfigure the build directory + +# do not use grep -oP as it is not POSIX compliant + +doReconfigure=1 + +log "${BLUE}[Info] Checking last build flags...${NC}" +if [[ -f $buildDir/.last_build_flags ]]; then + lastUserFlag=$(grep -Eo 'userFlag=[0-9]+' $buildDir/.last_build_flags | cut -d'=' -f2) + lastTestsFlag=$(grep -Eo 'testsFlag=[0-9]+' $buildDir/.last_build_flags | cut -d'=' -f2) + if [[ $lastUserFlag -eq $userFlag && $lastTestsFlag -eq $testsFlag ]]; then + log "${MAGENTA}[Succsess] Last build flags match current build flags. Skipping configuration...${NC}" + doReconfigure=0 else - if [[ $testsFlag -eq 1 ]]; then - meson setup "$buildDir" --buildtype=debug - else - meson setup "$buildDir" --buildtype=debug -Dbuild_tests=false - fi + log "${YELLOW}[Info] Last build flags do not match current build flags. Reconfiguring...${NC}" + rm -rf $buildDir fi - log "${GREEN}[Success] Build directory configured successfully.${NC}" +else + log "${BLUE}[Info] Last build flags not found. Reconfiguring...${NC}" +fi + + +# --- Check if the build dir has been configured with the correct flags --- +if [[ $doReconfigure -eq 1 ]]; then + log "${BLUE}[Info] Configuring build directory...${NC}" + if [ -f "$buildDir/build.ninja" ]; then + log "${MAGENTA}[Succsess] Build directory already configured. Skipping...${NC}" + else + if [[ $userFlag -eq 1 ]]; then + meson setup "$buildDir" -Duser_mode=true --buildtype=release + else + if [[ $testsFlag -eq 0 ]]; then + meson setup "$buildDir" --buildtype=debug + else + meson setup "$buildDir" --buildtype=debug -Dbuild_tests=false + fi + fi + log "${GREEN}[Success] Build directory configured successfully.${NC}" + fi + log "${BLUE}[Info] Caching last build flags...${NC}" + + # --- Cache the last build flags --- + if [[ $userFlag -eq 1 ]]; then + echo "userFlag=1" > $buildDir/.last_build_flags + else + echo "userFlag=0" > $buildDir/.last_build_flags + fi + if [[ $testsFlag -eq 1 ]]; then + echo "testsFlag=1" >> $buildDir/.last_build_flags + else + echo "testsFlag=0" >> $buildDir/.last_build_flags + fi + + log "${GREEN}[Success] Last build flags cached.${NC}" fi log "${BLUE}[Info] Building 4DSSE...${NC}" @@ -150,10 +209,10 @@ meson compile -C "$buildDir" # --- Install 4DSSE --- if [[ $userFlag -eq 1 ]]; then log "${BLUE}[Info] Installing 4DSSE in user mode...${NC}" - meson install -C "$buildDir" --no-rebuild --prefix="$installDir" + meson install -C "$buildDir" fi -if [[ $testsFlag -eq 1 ]]; then +if [[ $runTestsFlag -eq 1 ]]; then log "${BLUE}[Info] Running tests...${NC}" meson test -C "$buildDir" fi @@ -162,4 +221,7 @@ if [[ $userFlag -eq 1 ]]; then log "${GREEN}[Success] 4DSSE built and installed successfully.${NC}" else log "${GREEN}[Success] 4DSSE built successfully.${NC}" -fi \ No newline at end of file +fi + + +log "${GREEN}[Success] 4DSSE mk script complete.${NC}" From 30a795e6b1b9f97305600edbe677a27a2b2b3b01 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 21 Mar 2025 11:30:24 -0400 Subject: [PATCH 05/18] build(mk): added compiler and builder checks --- mk | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/mk b/mk index 28e2811..97a1e1e 100755 --- a/mk +++ b/mk @@ -105,6 +105,16 @@ if [[ -f "$LOGFILE" ]]; then fi touch "$LOGFILE" +# --- Check if Clang or GCC is installed --- +log "${BLUE}[Info] Checking if Clang or GCC is installed...${NC}" +if ! command -v clang &> /dev/null && ! command -v gcc &> /dev/null; then + log "${RED}[Error] Clang or GCC is not installed. Exiting...${NC}" + log "${YELLOW}[Info] Please install Clang or GCC and try again.${NC}" + exit 1 +else + log "${MAGENTA}[Success] Clang or GCC is installed. Continuing...${NC}" +fi + # --- Check if MESON is installed --- log "${BLUE}[Info] Checking if Meson is installed...${NC}" if ! command -v meson &> /dev/null; then @@ -115,6 +125,16 @@ else log "${MAGENTA}[Success] Meson is installed. Continuing...${NC}" fi +# --- Check if NINJA is installed --- +log "${BLUE}[Info] Checking if Ninja is installed...${NC}" +if ! command -v ninja &> /dev/null; then + log "${RED}[Error] Ninja is not installed. Exiting...${NC}" + log "${YELLOW}[Info] Please install Ninja and try again.${NC}" + exit 1 +else + log "${MAGENTA}[Success] Ninja is installed. Continuing...${NC}" +fi + # --- Build 4DSSE --- log "${BLUE}[Info] Building 4DSSE...${NC}" if [[ $userFlag -eq 1 ]]; then From 9f524cf821b17bba30f97bf2b59489c0443d3598 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 21 Mar 2025 11:39:02 -0400 Subject: [PATCH 06/18] build(mk): added meson auto installer --- mk | 58 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 56 insertions(+), 2 deletions(-) diff --git a/mk b/mk index 97a1e1e..cf93850 100755 --- a/mk +++ b/mk @@ -119,8 +119,62 @@ fi log "${BLUE}[Info] Checking if Meson is installed...${NC}" if ! command -v meson &> /dev/null; then log "${RED}[Error] Meson is not installed. Exiting...${NC}" - log "${YELLOW}[Info] Please install Meson and try again.${NC}" - exit 1 + # Check if the user would like to try to install meson by getting input + installMeson=0 + while true; do + read -p "Would you like to try to install Meson? [y/n]: " yn + case $yn in + [Yy]* ) installMeson=1; break;; + [Nn]* ) break;; + * ) echo "Please answer yes or no.";; + esac + done + if [[ $installMeson -eq 1 ]]; then + # check if pip is installed + if ! command -v pip &> /dev/null; then + log "${RED}[Error] pip is not installed. Exiting...${NC}" + log "${YELLOW}[Info] Please install pip and try again.${NC}" + exit 1 + else + log "${MAGENTA}[Success] pip is installed. Continuing...${NC}" + fi + + # check if python3 is installed + if ! command -v python3 &> /dev/null; then + log "${RED}[Error] python3 is not installed. Exiting...${NC}" + log "${YELLOW}[Info] Please install python3 and try again.${NC}" + log "${YELLOW}[INFO] If you have python3 installed, please add it to your PATH.${NC}" + exit 1 + else + log "${MAGENTA}[Success] python3 is installed. Continuing...${NC}" + fi + + # Check if the venv ~/.4DSSE_env exists + if [[ ! -d "$HOME/.4DSSE_env" ]]; then + log "${BLUE}[Info] Creating virtual environment...${NC}" + python3 -m venv "$HOME/.4DSSE_env" + source "$HOME/.4DSSE_env/bin/activate" + log "${GREEN}[Success] Virtual environment created.${NC}" + else + log "${MAGENTA}[Succsess] Virtual environment already exists. Skipping...${NC}" + fi + + # install meson + log "${BLUE}[Info] Installing Meson...${NC}" + pip install meson + + # confirm meson is installed + if ! command -v meson &> /dev/null; then + log "${RED}[Error] Meson did not install properly. Exiting...${NC}" + log "${YELLOW}[Info] Please manually install Meson and try again.${NC}" + exit 1 + else + log "${GREEN}[Success] Meson installed.${NC}" + fi + else + log "${YELLOW}[Info] Please install Meson and try again.${NC}" + exit 1 + fi else log "${MAGENTA}[Success] Meson is installed. Continuing...${NC}" fi From 58ef0d495e528beac2ed80891a4f57b9c8e0b13f Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 21 Mar 2025 11:42:01 -0400 Subject: [PATCH 07/18] build(mk): added auto ninja installer --- mk | 59 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 57 insertions(+), 2 deletions(-) diff --git a/mk b/mk index cf93850..62c03c0 100755 --- a/mk +++ b/mk @@ -183,8 +183,63 @@ fi log "${BLUE}[Info] Checking if Ninja is installed...${NC}" if ! command -v ninja &> /dev/null; then log "${RED}[Error] Ninja is not installed. Exiting...${NC}" - log "${YELLOW}[Info] Please install Ninja and try again.${NC}" - exit 1 + + # Check if the user would like to try to install ninja by getting input + installNinja=0 + while true; do + read -p "Would you like to try to install Ninja? [y/n]: " yn + case $yn in + [Yy]* ) installNinja=1; break;; + [Nn]* ) break;; + * ) echo "Please answer yes or no.";; + esac + done + if [[ $installNinja -eq 1 ]]; then + # check if pip is installed + if ! command -v pip &> /dev/null; then + log "${RED}[Error] pip is not installed. Exiting...${NC}" + log "${YELLOW}[Info] Please install pip and try again.${NC}" + exit 1 + else + log "${MAGENTA}[Success] pip is installed. Continuing...${NC}" + fi + + # check if python3 is installed + if ! command -v python3 &> /dev/null; then + log "${RED}[Error] python3 is not installed. Exiting...${NC}" + log "${YELLOW}[Info] Please install python3 and try again.${NC}" + log "${YELLOW}[INFO] If you have python3 installed, please add it to your PATH.${NC}" + exit 1 + else + log "${MAGENTA}[Success] python3 is installed. Continuing...${NC}" + fi + + # Check if the venv ~/.4DSSE_env exists + if [[ ! -d "$HOME/.4DSSE_env" ]]; then + log "${BLUE}[Info] Creating virtual environment...${NC}" + python3 -m venv "$HOME/.4DSSE_env" + source "$HOME/.4DSSE_env/bin/activate" + log "${GREEN}[Success] Virtual environment created.${NC}" + else + log "${MAGENTA}[Succsess] Virtual environment already exists. Skipping...${NC}" + fi + + # install ninja + log "${BLUE}[Info] Installing Ninja...${NC}" + pip install ninja + + # confirm ninja is installed + if ! command -v ninja &> /dev/null; then + log "${RED}[Error] Ninja did not install properly. Exiting...${NC}" + log "${YELLOW}[Info] Please manually install Ninja and try again.${NC}" + exit 1 + else + log "${GREEN}[Success] Ninja installed.${NC}" + fi + else + log "${YELLOW}[Info] Please install Ninja and try again.${NC}" + exit 1 + fi else log "${MAGENTA}[Success] Ninja is installed. Continuing...${NC}" fi From 29c2ca3ca5f6dc6f75fd1d35293e1ceebd28151d Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 21 Mar 2025 11:50:16 -0400 Subject: [PATCH 08/18] build(mk): added auto cmake detect --- mk | 82 ++++++++++++++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 72 insertions(+), 10 deletions(-) diff --git a/mk b/mk index 62c03c0..008c900 100755 --- a/mk +++ b/mk @@ -122,7 +122,7 @@ if ! command -v meson &> /dev/null; then # Check if the user would like to try to install meson by getting input installMeson=0 while true; do - read -p "Would you like to try to install Meson? [y/n]: " yn + read -p "${YELLOW}[Query] Would you like to try to install Meson? [y/n]: ${NC}" yn case $yn in [Yy]* ) installMeson=1; break;; [Nn]* ) break;; @@ -136,7 +136,7 @@ if ! command -v meson &> /dev/null; then log "${YELLOW}[Info] Please install pip and try again.${NC}" exit 1 else - log "${MAGENTA}[Success] pip is installed. Continuing...${NC}" + log "${MAGENTA}[Succsess] pip is installed. Continuing...${NC}" fi # check if python3 is installed @@ -146,7 +146,7 @@ if ! command -v meson &> /dev/null; then log "${YELLOW}[INFO] If you have python3 installed, please add it to your PATH.${NC}" exit 1 else - log "${MAGENTA}[Success] python3 is installed. Continuing...${NC}" + log "${MAGENTA}[Succsess] python3 is installed. Continuing...${NC}" fi # Check if the venv ~/.4DSSE_env exists @@ -154,7 +154,7 @@ if ! command -v meson &> /dev/null; then log "${BLUE}[Info] Creating virtual environment...${NC}" python3 -m venv "$HOME/.4DSSE_env" source "$HOME/.4DSSE_env/bin/activate" - log "${GREEN}[Success] Virtual environment created.${NC}" + log "${GREEN}[Succsess] Virtual environment created.${NC}" else log "${MAGENTA}[Succsess] Virtual environment already exists. Skipping...${NC}" fi @@ -169,14 +169,14 @@ if ! command -v meson &> /dev/null; then log "${YELLOW}[Info] Please manually install Meson and try again.${NC}" exit 1 else - log "${GREEN}[Success] Meson installed.${NC}" + log "${GREEN}[Succsess] Meson installed.${NC}" fi else log "${YELLOW}[Info] Please install Meson and try again.${NC}" exit 1 fi else - log "${MAGENTA}[Success] Meson is installed. Continuing...${NC}" + log "${MAGENTA}[Succsess] Meson is installed. Continuing...${NC}" fi # --- Check if NINJA is installed --- @@ -187,7 +187,7 @@ if ! command -v ninja &> /dev/null; then # Check if the user would like to try to install ninja by getting input installNinja=0 while true; do - read -p "Would you like to try to install Ninja? [y/n]: " yn + read -p "${YELLOW}[Query] Would you like to try to install Ninja? [y/n]: ${NC}" yn case $yn in [Yy]* ) installNinja=1; break;; [Nn]* ) break;; @@ -201,7 +201,7 @@ if ! command -v ninja &> /dev/null; then log "${YELLOW}[Info] Please install pip and try again.${NC}" exit 1 else - log "${MAGENTA}[Success] pip is installed. Continuing...${NC}" + log "${MAGENTA}[Succsess] pip is installed. Continuing...${NC}" fi # check if python3 is installed @@ -211,7 +211,7 @@ if ! command -v ninja &> /dev/null; then log "${YELLOW}[INFO] If you have python3 installed, please add it to your PATH.${NC}" exit 1 else - log "${MAGENTA}[Success] python3 is installed. Continuing...${NC}" + log "${MAGENTA}[Succsess] python3 is installed. Continuing...${NC}" fi # Check if the venv ~/.4DSSE_env exists @@ -219,7 +219,7 @@ if ! command -v ninja &> /dev/null; then log "${BLUE}[Info] Creating virtual environment...${NC}" python3 -m venv "$HOME/.4DSSE_env" source "$HOME/.4DSSE_env/bin/activate" - log "${GREEN}[Success] Virtual environment created.${NC}" + log "${GREEN}[Succsess] Virtual environment created.${NC}" else log "${MAGENTA}[Succsess] Virtual environment already exists. Skipping...${NC}" fi @@ -244,6 +244,68 @@ else log "${MAGENTA}[Success] Ninja is installed. Continuing...${NC}" fi +# --- Check if CMake is installed --- +log "${BLUE}[Info] Checking if CMake is installed...${NC}" +if ! command -v cmake &> /dev/null; then + log "${RED}[Error] CMake is not installed. Exiting...${NC}" + + # Check if the user would like to try to install cmake by getting input + installCMake=0 + while true; do + read -p "${YELLOW}[Query] Would you like to try to install CMake? [y/n]: ${NC}" yn + case $yn in + [Yy]* ) installCMake=1; break;; + [Nn]* ) break;; + * ) echo "Please answer yes or no.";; + esac + done + if [[ $installCMake -eq 1 ]]; then + # check if pip is installed + if ! command -v pip &> /dev/null; then + log "${RED}[Error] pip is not installed. Exiting...${NC}" + log "${YELLOW}[Info] Please install pip and try again.${NC}" + exit 1 + else + log "${MAGENTA}[Succsess] pip is installed. Continuing...${NC}" + fi + + # check if python3 is installed + if ! command -v python3 &> /dev/null; then + log "${RED}[Error] python3 is not installed. Exiting...${NC}" + log "${YELLOW}[Info] Please install python3 and try again.${NC}" + log "${YELLOW}[INFO] If you have python3 installed, please add it to your PATH.${NC}" + exit 1 + else + log "${MAGENTA}[Succsess] python3 is installed. Continuing...${NC}" + fi + + # Check if the venv ~/.4DSSE_env exists + if [[ ! -d "$HOME/.4DSSE_env" ]]; then + log "${BLUE}[Info] Creating virtual environment...${NC}" + python3 -m venv "$HOME/.4DSSE_env" + source "$HOME/.4DSSE_env/bin/activate" + log "${GREEN}[Succsess] Virtual environment created.${NC}" + else + log "${MAGENTA}[Succsess] Virtual environment already exists. Skipping...${NC}" + fi + + # install cmake + log "${BLUE}[Info] Installing CMake...${NC}" + pip install cmake + + # confirm cmake is installed + if ! command -v cmake &> /dev/null; then + log "${RED}[Error] CMake did not install properly. Exiting...${NC}" + log "${YELLOW}[Info] Please manually install CMake and try again.${NC}" + exit 1 + else + log "${GREEN}[Success] CMake installed.${NC}" + fi + fi +else + log "${MAGENTA}[Success] CMake is installed. Continuing...${NC}" +fi + # --- Build 4DSSE --- log "${BLUE}[Info] Building 4DSSE...${NC}" if [[ $userFlag -eq 1 ]]; then From 7c97afc9a283c476c85d7b2a9d5db6dc37077bd1 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 21 Mar 2025 11:53:34 -0400 Subject: [PATCH 09/18] build(mk): updated query colors --- mk | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/mk b/mk index 008c900..f28f541 100755 --- a/mk +++ b/mk @@ -122,7 +122,8 @@ if ! command -v meson &> /dev/null; then # Check if the user would like to try to install meson by getting input installMeson=0 while true; do - read -p "${YELLOW}[Query] Would you like to try to install Meson? [y/n]: ${NC}" yn + echo -ne "${YELLOW}[Query] Would you like to try to install Meson? [y/n]: ${NC}" + read -p "" yn case $yn in [Yy]* ) installMeson=1; break;; [Nn]* ) break;; @@ -187,7 +188,8 @@ if ! command -v ninja &> /dev/null; then # Check if the user would like to try to install ninja by getting input installNinja=0 while true; do - read -p "${YELLOW}[Query] Would you like to try to install Ninja? [y/n]: ${NC}" yn + echo -ne "${YELLOW}[Query] Would you like to try to install Ninja? [y/n]: ${NC}" + read -p "" yn case $yn in [Yy]* ) installNinja=1; break;; [Nn]* ) break;; @@ -247,12 +249,13 @@ fi # --- Check if CMake is installed --- log "${BLUE}[Info] Checking if CMake is installed...${NC}" if ! command -v cmake &> /dev/null; then - log "${RED}[Error] CMake is not installed. Exiting...${NC}" + log "${RED}[Error] CMake is not installed....${NC}" # Check if the user would like to try to install cmake by getting input installCMake=0 while true; do - read -p "${YELLOW}[Query] Would you like to try to install CMake? [y/n]: ${NC}" yn + echo -ne "${YELLOW}[Query] Would you like to try to install CMake? [y/n]: ${NC}" + read -p "" yn case $yn in [Yy]* ) installCMake=1; break;; [Nn]* ) break;; From 2a7cd1e9e606ca3a549965994c55edbd07a7172e Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 21 Mar 2025 11:54:33 -0400 Subject: [PATCH 10/18] build(boost): lowered required version --- build-config/boost/meson.build | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/build-config/boost/meson.build b/build-config/boost/meson.build index fa83c8f..321216a 100644 --- a/build-config/boost/meson.build +++ b/build-config/boost/meson.build @@ -1 +1 @@ -boost_dep = dependency('boost', version: '1.87.0') +boost_dep = dependency('boost', version: '1.8.0') From 5d3c0d45ecaff5aa2751d6a72f3ad122af2803a8 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 21 Mar 2025 11:56:12 -0400 Subject: [PATCH 11/18] build(boost): added greater than to version for boost --- build-config/boost/meson.build | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/build-config/boost/meson.build b/build-config/boost/meson.build index 321216a..f2b581a 100644 --- a/build-config/boost/meson.build +++ b/build-config/boost/meson.build @@ -1 +1 @@ -boost_dep = dependency('boost', version: '1.8.0') +boost_dep = dependency('boost', version: '>=1.8.0') From 16613288766c1b394f1d95964cc7ecef04b50e70 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 21 Mar 2025 12:03:48 -0400 Subject: [PATCH 12/18] refactor(mk): fixed some typos --- mk | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/mk b/mk index f28f541..8a8f521 100755 --- a/mk +++ b/mk @@ -112,7 +112,7 @@ if ! command -v clang &> /dev/null && ! command -v gcc &> /dev/null; then log "${YELLOW}[Info] Please install Clang or GCC and try again.${NC}" exit 1 else - log "${MAGENTA}[Success] Clang or GCC is installed. Continuing...${NC}" + log "${MAGENTA}[Succsess] Clang or GCC is installed. Continuing...${NC}" fi # --- Check if MESON is installed --- @@ -236,14 +236,14 @@ if ! command -v ninja &> /dev/null; then log "${YELLOW}[Info] Please manually install Ninja and try again.${NC}" exit 1 else - log "${GREEN}[Success] Ninja installed.${NC}" + log "${GREEN}[Succsess] Ninja installed.${NC}" fi else log "${YELLOW}[Info] Please install Ninja and try again.${NC}" exit 1 fi else - log "${MAGENTA}[Success] Ninja is installed. Continuing...${NC}" + log "${MAGENTA}[Succsess] Ninja is installed. Continuing...${NC}" fi # --- Check if CMake is installed --- @@ -302,11 +302,11 @@ if ! command -v cmake &> /dev/null; then log "${YELLOW}[Info] Please manually install CMake and try again.${NC}" exit 1 else - log "${GREEN}[Success] CMake installed.${NC}" + log "${GREEN}[Succsess] CMake installed.${NC}" fi fi else - log "${MAGENTA}[Success] CMake is installed. Continuing...${NC}" + log "${MAGENTA}[Succsess] CMake is installed. Continuing...${NC}" fi # --- Build 4DSSE --- @@ -335,7 +335,7 @@ else exit 1 else touch ./build-config/boost/.boost_installed - log "${GREEN}[Success] Boost check passed.${NC}" + log "${GREEN}[Succsess] Boost check passed.${NC}" fi fi @@ -378,7 +378,7 @@ if [[ $doReconfigure -eq 1 ]]; then meson setup "$buildDir" --buildtype=debug -Dbuild_tests=false fi fi - log "${GREEN}[Success] Build directory configured successfully.${NC}" + log "${GREEN}[Succsess] Build directory configured successfully.${NC}" fi log "${BLUE}[Info] Caching last build flags...${NC}" @@ -394,7 +394,7 @@ if [[ $doReconfigure -eq 1 ]]; then echo "testsFlag=0" >> $buildDir/.last_build_flags fi - log "${GREEN}[Success] Last build flags cached.${NC}" + log "${GREEN}[Succsess] Last build flags cached.${NC}" fi log "${BLUE}[Info] Building 4DSSE...${NC}" @@ -412,10 +412,10 @@ if [[ $runTestsFlag -eq 1 ]]; then fi if [[ $userFlag -eq 1 ]]; then - log "${GREEN}[Success] 4DSSE built and installed successfully.${NC}" + log "${GREEN}[Succsess] 4DSSE built and installed successfully.${NC}" else - log "${GREEN}[Success] 4DSSE built successfully.${NC}" + log "${GREEN}[Succsess] 4DSSE built successfully.${NC}" fi -log "${GREEN}[Success] 4DSSE mk script complete.${NC}" +log "${GREEN}[Succsess] 4DSSE mk script complete.${NC}" From c58bd50f013d8dfa805c4a50c3ffccc7ab8aa92f Mon Sep 17 00:00:00 2001 From: Aaron Dotter Date: Fri, 21 Mar 2025 12:08:40 -0400 Subject: [PATCH 13/18] refactor(network) Added header file --- src/network/public/network.h | 526 +++++++++++++++++++++++++++++++++++ 1 file changed, 526 insertions(+) diff --git a/src/network/public/network.h b/src/network/public/network.h index e69de29..3466aa0 100644 --- a/src/network/public/network.h +++ b/src/network/public/network.h @@ -0,0 +1,526 @@ + +#include +#include +#include +#include + +#include +#include +#include + +#include "const.h" + +/* Nuclear reaction network in cgs units based on Frank Timmes' "aprox8". + 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 . + +*/ + + +// using namespace std; +using namespace boost::numeric::odeint; +namespace phoenix = boost::phoenix; + +// these types are required by the Rosenbrock implicit solver +typedef boost::numeric::ublas::vector< double > vector_type; +typedef boost::numeric::ublas::matrix< double > matrix_type; + +// this array is used only in the nuclear reaction rate evaluations +typedef std::array vec; + +// only need a couple of constants +Constants& constants = Constants::getInstance(); +const double avo = constants.get("N_a").value; +const double clight = constants.get("c").value; + +// identify the isotopes used in the network. +const int ih1=0; +const int ihe3=1; +const int ihe4=2; +const int ic12=3; +const int in14=4; +const int io16=5; +const int ine20=6; +const int img24=7; + +// physical variables; this routine currently does not need to call EOS +// since the temperature and density are assumed constant during the burn +const int itemp=img24+1; +const int iden =itemp+1; +const int iener=iden+1; + +// number of isotopes and number of variables +const int niso=img24+1; // number of isotopes +const int nvar=iener+1; // number of variables + +// atomic stuff +std::array aion = {1,3,4,12,14,16,20,24}; +//std::array zion = {1,2,2, 6, 7, 8,10,12}; +//std::array bion = {0,7.71819,28.29603,92.16294,104.65998,127.62093,160.64788,198.25790}; +//nion = aion - zion #neutrons +//mion = nion*mn + zion*mp - bion*mev2gr #mass +std::array mion = {1.67262164e-24, 5.00641157e-24, 6.64465545e-24, 1.99209977e-23, + 2.32462686e-23, 2.65528858e-23, 3.31891077e-23, 3.98171594e-23}; + +//helper functions +// a function to multilpy two arrays and then sum the resulting elements: sum(a*b) +double sum_product( const vec &a, const vec &b){ + if (a.size() != b.size()){ + throw std::runtime_error("Error: array size mismatch in sum_product"); + } + + 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 +vec get_T9_array(const double &T) { + vec arr; + double T9=1e-9*T; + 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 vec &T9, const vec &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 vec &T9) { + vec a1 = {-34.78630, 0,-3.511930, 3.100860, -0.1983140, 1.262510e-2, -1.025170}; + vec 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 vec &T9) { + vec a1 = {7.528980, 0, -3.720800, 0.8717820, 0, 0,-0.6666670}; + vec 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 vec &T9){ + vec 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 vec &T9){ + vec a1 = {1.560990e+01,0.000000e+00,-1.282710e+01,-3.082250e-02,-6.546850e-01,8.963310e-02,-6.666670e-01}; + vec 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 vec &T9){ + vec a1 = {-9.710520e-01,0.000000e+00,-3.706000e+01,2.934930e+01,-1.155070e+02,-1.000000e+01,-1.333330e+00}; + vec a2 = {-1.178840e+01,-1.024460e+00,-2.357000e+01,2.048860e+01,-1.298820e+01,-2.000000e+01,-2.166670e+00}; + vec 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 vec &T9){ + vec a1={1.714820e+01,0.000000e+00,-1.369200e+01,-2.308810e-01,4.443620e+00,-3.158980e+00,-6.666670e-01}; + vec 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 vec &T9){ + vec a1={6.965260e+01,-1.392540e+00,5.891280e+01,-1.482730e+02,9.083240e+00,-5.410410e-01,7.035540e+01}; + vec 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 vec &T9){ + vec a1={1.701000e+01,0.000000e+00,-1.519300e+01,-1.619540e-01,-7.521230e+00,-9.875650e-01,-6.666670e-01}; + vec a2={2.011690e+01,0.000000e+00,-1.519300e+01,-4.639750e+00,9.734580e+00,-9.550510e+00,3.333330e-01}; + vec a3={7.654440e+00,-2.998000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + vec 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 vec &T9){ + vec a1={2.153390e+01,0.000000e+00,-3.625040e+01,0.000000e+00,0.000000e+00,-5.000000e+00,-6.666670e-01}; + vec a2={1.968380e-01,-5.160340e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + vec 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 vec &T9){ + vec a1 = {2.747640e+01,0.000000e+00,-1.525300e+01,1.593180e+00,2.447900e+00,-2.197080e+00,-6.666670e-01}; + vec a2 = {-4.873470e+00,-2.021170e+00,0.000000e+00,3.084970e+01,-8.504330e+00,-1.544260e+00,-1.500000e+00}; + vec a3 = {2.089720e+01,-7.406000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + vec 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 vec &T9){ + vec a1 = {2.001760e+01,0.000000e+00,-1.524000e+01,3.349260e-01,4.590880e+00,-4.784680e+00,-6.666670e-01}; + vec a2 = {6.590560e+00,-2.923150e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + vec 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 vec &T9){ + double f1=n15pg_rate(T9); + double f2=n15pa_rate(T9); + return f1/(f1+f2); +} + +// o16(p,g)f17 then f17 -> o17(p,a)n14 +double o16p_rate(const vec &T9){ + vec 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 vec &T9){ + vec a1={2.390300e+01,0.000000e+00,-3.972620e+01,-2.107990e-01,4.428790e-01,-7.977530e-02,-6.666670e-01}; + vec a2={3.885710e+00,-1.035850e+01,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + vec 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 vec &T9){ + vec a1={2.450580e+01,0.000000e+00,-4.625250e+01,5.589010e+00,7.618430e+00,-3.683000e+00,-6.666670e-01}; + vec a2={-3.870550e+01,-2.506050e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + vec a3={1.983070e+00,-9.220260e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + vec 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 vec &T9){ + vec 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 vec &T9){ + vec 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 +struct Jacobian { + + void operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) + { + // EOS + vec T9=get_T9_array(y[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); + + double pfrac=n15pg_frac(T9); + double afrac=1-pfrac; + + double yh1 = y[ ih1]; + double yhe3 = y[ ihe3]; + double yhe4 = y[ ihe4]; + double yc12 = y[ ic12]; + double yn14 = y[ in14]; + double yo16 = y[ io16]; + double yne20 = y[ine20]; + + // zero all elements to begin + for (int i=0; i < nvar; i++) { + for (int j=0; j < nvar; j++) { J(i,j)=0.0; }} + + // h1 jacobian elements + J(ih1,ih1) = -3*yh1*rpp - 2*yc12*rc12p -2*yn14*rn14p -2*yo16*ro16p; + J(ih1,ihe3) = 2*yhe3*r33 - yhe4*r34; + J(ih1,ihe4) = -yhe3*r34; + J(ih1,ic12) = -2*yh1*rc12p; + J(ih1,in14) = -2*yh1*rn14p; + J(ih1,io16) = -2*yh1*ro16p; + + // he3 jacobian elements + J(ihe3,ih1) = yh1*rpp; + J(ihe3,ihe3) = -2*yhe3*r33 - yhe4*r34; + J(ihe3,ihe4) = -yhe3*r34; + + // he4 jacobian elements + J(ihe4,ih1) = yn14*afrac*rn14p + yo16*ro16p; + J(ihe4,ihe3) = yhe3*r33 - yhe4*r34; + J(ihe4,ihe4) = yhe3*r34 - 1.5*yhe4*yhe4*r3a - yc12*rc12a - 1.5*yn14*rn14a - yo16*ro16a - yne20*rne20a; + J(ihe4,ic12) = -yhe4*rc12a + yc12*r1212 + yo16*r1216; + J(ihe4,in14) = yh1*afrac*rn14p - 1.5*yhe4*rn14a; + J(ihe4,io16) = yh1*ro16p - yhe4*ro16a + yc12*r1216; + J(ihe4,ine20) = -yhe4*rne20a; + + // c12 jacobian elements + J(ic12,ih1) = -yc12*rc12p + yn14*afrac*rn14p; + J(ic12,ihe4) = 0.5*yhe4*yhe4*r3a - yhe4*rc12a; + J(ic12,ic12) = -yh1*rc12p - yhe4*rc12a - yo16*r1216 - 2*yc12*r1212; + J(ic12,in14) = yh1*yn14*afrac*rn14p; + J(ic12,io16) = -yc12*r1216; + + // n14 jacobian elements + J(in14,ih1) = yc12*rc12p - yn14*rn14p + yo16*ro16p; + J(in14,ihe4) = -yn14*rn14a; + J(in14,ic12) = yh1*rc12p; + J(in14,in14) = -yh1*rn14p - yhe4*rn14a; + J(in14,io16) = yo16*ro16p; + + // o16 jacobian elements + J(io16,ih1) = yn14*pfrac*rn14p - yo16*ro16p; + J(io16,ihe4) = yc12*rc12a - yo16*ro16a; + J(io16,ic12) = yhe4*rc12a - yo16*r1216; + J(io16,in14) = yh1*pfrac*rn14p; + J(io16,io16) = yh1*ro16p - yc12*r1216 -yhe4*ro16a; + + // ne20 jacobian elements + J(ine20,ihe4) = yn14*rn14a + yo16*ro16a - yne20*rne20a; + J(ine20,ic12) = yc12*r1212; + J(ine20,in14) = yhe4*rn14a; + J(ine20,io16) = yo16*ro16a; + J(ine20,ine20) = -yhe4*rne20a; + + // mg24 jacobian elements + J(img24,ihe4) = yne20*rne20a; + J(img24,ic12) = yo16*r1216; + J(img24,io16) = yc12*r1216; + J(img24,ine20) = yhe4*rne20a; + + // energy accounting + for (int j=0; j>(1.0e-6,1.0e-6) , + std::make_pair( ODE(), Jacobian() ), y, 0.0, tmax, dt0); + + } else { + std::cout << " *** Explicit RK Dormand-Prince *** " << std::endl; + + num_of_steps = integrate_const( make_dense_output>(1.0e-6, 1.0e-6), + ODE(), y, 0.0, tmax, dt0); + + } + + //convert number fraction to mass fraction + ysum=0; + for (int i=0; i < niso; i++) { + y[i] *= aion[i]; + ysum+= y[i]; + } + for (int i=0; i < niso; i++) {y[i] /= ysum;} + + std::cout << " H1: " << y[ih1] << std::endl; + std::cout << " He4: " << y[ihe4] << std::endl; + std::cout << "energy: " << y[iener] << std::endl; + std::cout << " steps: " << num_of_steps << std::endl; + return 0; + +} From 7c40db4b09db1c1524d71ac78ed550c2ae99ed89 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 21 Mar 2025 14:03:18 -0400 Subject: [PATCH 14/18] feat(approx8-network-integrated): added network handleing semantics and incorporated the approx8 network into them --- src/network/meson.build | 10 +- src/network/private/approx8.cpp | 525 +++++++++++++++++++++++++++ src/network/private/network.cpp | 536 +-------------------------- src/network/public/approx8.h | 316 ++++++++++++++++ src/network/public/network.h | 619 +++++--------------------------- 5 files changed, 954 insertions(+), 1052 deletions(-) create mode 100644 src/network/private/approx8.cpp create mode 100644 src/network/public/approx8.h diff --git a/src/network/meson.build b/src/network/meson.build index 72f5bd1..2add636 100644 --- a/src/network/meson.build +++ b/src/network/meson.build @@ -1,15 +1,21 @@ # Define the library network_sources = files( 'private/network.cpp', + 'private/approx8.cpp' ) network_headers = files( - 'public/network.h' + 'public/network.h', + 'public/approx8.h' ) dependencies = [ boost_dep, - const_dep + const_dep, + quill_dep, + mfem_dep, + config_dep, + probe_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 new file mode 100644 index 0000000..ec0abc8 --- /dev/null +++ b/src/network/private/approx8.cpp @@ -0,0 +1,525 @@ +#include +#include +#include + +#include +#include +#include + +#include "const.h" +#include "config.h" +#include "quill/LogMacros.h" + +#include "approx8.h" +#include "network.h" + +/* Nuclear reaction network in cgs units based on Frank Timmes' "aprox8". + 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 nnApprox8{ + + // using namespace std; + using namespace boost::numeric::odeint; + + //helper functions + // a function to multilpy 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"); + } + + 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; + double T9=1e-9*T; + 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) { + 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}; + 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}; + 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}; + 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}; + 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}; + 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}; + 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}; + 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}; + 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}; + 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}; + 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}; + 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); + 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}; + 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}; + 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}; + 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}; + 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}; + 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 ) { + Constants& constants = 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]); + + // 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); + + double pfrac=n15pg_frac(T9); + 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]; + + // 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; + } + } + + // 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; + + // 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; + + // 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; + + // 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; + + // 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; + + // 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; + + // 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; + + // 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; + + // energy accountings + for (int j=0; j("Network:Approx8:Stiff:AbsTol", 1.0e-6); + const double stiff_rel_tol = m_config.get("Network:Approx8:Stiff:RelTol", 1.0e-6); + const double nonstiff_abs_tol = m_config.get("Network:Approx8:NonStiff:AbsTol", 1.0e-6); + const double nonstiff_rel_tol = m_config.get("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>(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>(nonstiff_abs_tol, nonstiff_rel_tol), + ODE(), + m_y, + 0.0, + 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]; + } + for (int i = 0; i < Net::niso; i++) { + m_y[i] /= ysum; + } + + nuclearNetwork::NetOut netOut; + std::vector outComposition; + outComposition.reserve(Net::nvar); + + for (int i = 0; i < Net::niso; i++) { + outComposition.push_back(m_y[i]); + } + netOut.energy = m_y[Net::iener]; + netOut.composition = outComposition; + netOut.num_steps = num_steps; + + return netOut; + } + + void Approx8Network::setStiff(bool stiff) { + m_stiff = stiff; + } + + vector_type Approx8Network::convert_netIn(const nuclearNetwork::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(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]; + } + for (int i = 0; i < Net::niso; i++) { + y[i] /= ysum; + } + + return y; + } +}; + + +// main program + diff --git a/src/network/private/network.cpp b/src/network/private/network.cpp index 5a9ff2b..5ffaac7 100644 --- a/src/network/private/network.cpp +++ b/src/network/private/network.cpp @@ -1,528 +1,16 @@ -#include -#include -#include -#include +#include "network.h" +#include "probe.h" +#include "quill/LogMacros.h" -#include -#include -#include - -#include "const.h" -//#include "probe.h" -//#include "config.h" -//#include "quill/LogMacros.h" - -/* Nuclear reaction network in cgs units based on Frank Timmes' "aprox8". - 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 . - -*/ - - -// using namespace std; -using namespace boost::numeric::odeint; -namespace phoenix = boost::phoenix; - -// these types are required by the Rosenbrock implicit solver -typedef boost::numeric::ublas::vector< double > vector_type; -typedef boost::numeric::ublas::matrix< double > matrix_type; - -// this array is used only in the nuclear reaction rate evaluations -typedef std::array vec; - -// only need a couple of constants -Constants& constants = Constants::getInstance(); -const double avo = constants.get("N_a").value; -const double clight = constants.get("c").value; - -// identify the isotopes used in the network. -const int ih1=0; -const int ihe3=1; -const int ihe4=2; -const int ic12=3; -const int in14=4; -const int io16=5; -const int ine20=6; -const int img24=7; - -// physical variables; this routine currently does not need to call EOS -// since the temperature and density are assumed constant during the burn -const int itemp=img24+1; -const int iden =itemp+1; -const int iener=iden+1; - -// number of isotopes and number of variables -const int niso=img24+1; // number of isotopes -const int nvar=iener+1; // number of variables - -// atomic stuff -std::array aion = {1,3,4,12,14,16,20,24}; -//std::array zion = {1,2,2, 6, 7, 8,10,12}; -//std::array bion = {0,7.71819,28.29603,92.16294,104.65998,127.62093,160.64788,198.25790}; -//nion = aion - zion #neutrons -//mion = nion*mn + zion*mp - bion*mev2gr #mass -std::array mion = {1.67262164e-24, 5.00641157e-24, 6.64465545e-24, 1.99209977e-23, - 2.32462686e-23, 2.65528858e-23, 3.31891077e-23, 3.98171594e-23}; - -//helper functions -// a function to multilpy two arrays and then sum the resulting elements: sum(a*b) -double sum_product( const vec &a, const vec &b){ - if (a.size() != b.size()){ - throw std::runtime_error("Error: array size mismatch in sum_product"); +namespace nuclearNetwork { + Network::Network() : + m_config(Config::getInstance()), + m_logManager(Probe::LogManager::getInstance()), + m_logger(m_logManager.getLogger("log")) { } - - double sum=0; - for (size_t i=0; i < a.size(); i++) { - sum += a[i] * b[i]; + nuclearNetwork::NetOut nuclearNetwork::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"); } - 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 -vec get_T9_array(const double &T) { - vec arr; - double T9=1e-9*T; - 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 vec &T9, const vec &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 vec &T9) { - vec a1 = {-34.78630, 0,-3.511930, 3.100860, -0.1983140, 1.262510e-2, -1.025170}; - vec 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 vec &T9) { - vec a1 = {7.528980, 0, -3.720800, 0.8717820, 0, 0,-0.6666670}; - vec 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 vec &T9){ - vec 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 vec &T9){ - vec a1 = {1.560990e+01,0.000000e+00,-1.282710e+01,-3.082250e-02,-6.546850e-01,8.963310e-02,-6.666670e-01}; - vec 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 vec &T9){ - vec a1 = {-9.710520e-01,0.000000e+00,-3.706000e+01,2.934930e+01,-1.155070e+02,-1.000000e+01,-1.333330e+00}; - vec a2 = {-1.178840e+01,-1.024460e+00,-2.357000e+01,2.048860e+01,-1.298820e+01,-2.000000e+01,-2.166670e+00}; - vec 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 vec &T9){ - vec a1={1.714820e+01,0.000000e+00,-1.369200e+01,-2.308810e-01,4.443620e+00,-3.158980e+00,-6.666670e-01}; - vec 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 vec &T9){ - vec a1={6.965260e+01,-1.392540e+00,5.891280e+01,-1.482730e+02,9.083240e+00,-5.410410e-01,7.035540e+01}; - vec 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 vec &T9){ - vec a1={1.701000e+01,0.000000e+00,-1.519300e+01,-1.619540e-01,-7.521230e+00,-9.875650e-01,-6.666670e-01}; - vec a2={2.011690e+01,0.000000e+00,-1.519300e+01,-4.639750e+00,9.734580e+00,-9.550510e+00,3.333330e-01}; - vec a3={7.654440e+00,-2.998000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec 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 vec &T9){ - vec a1={2.153390e+01,0.000000e+00,-3.625040e+01,0.000000e+00,0.000000e+00,-5.000000e+00,-6.666670e-01}; - vec a2={1.968380e-01,-5.160340e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec 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 vec &T9){ - vec a1 = {2.747640e+01,0.000000e+00,-1.525300e+01,1.593180e+00,2.447900e+00,-2.197080e+00,-6.666670e-01}; - vec a2 = {-4.873470e+00,-2.021170e+00,0.000000e+00,3.084970e+01,-8.504330e+00,-1.544260e+00,-1.500000e+00}; - vec a3 = {2.089720e+01,-7.406000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec 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 vec &T9){ - vec a1 = {2.001760e+01,0.000000e+00,-1.524000e+01,3.349260e-01,4.590880e+00,-4.784680e+00,-6.666670e-01}; - vec a2 = {6.590560e+00,-2.923150e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec 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 vec &T9){ - double f1=n15pg_rate(T9); - double f2=n15pa_rate(T9); - return f1/(f1+f2); -} - -// o16(p,g)f17 then f17 -> o17(p,a)n14 -double o16p_rate(const vec &T9){ - vec 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 vec &T9){ - vec a1={2.390300e+01,0.000000e+00,-3.972620e+01,-2.107990e-01,4.428790e-01,-7.977530e-02,-6.666670e-01}; - vec a2={3.885710e+00,-1.035850e+01,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec 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 vec &T9){ - vec a1={2.450580e+01,0.000000e+00,-4.625250e+01,5.589010e+00,7.618430e+00,-3.683000e+00,-6.666670e-01}; - vec a2={-3.870550e+01,-2.506050e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec a3={1.983070e+00,-9.220260e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec 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 vec &T9){ - vec 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 vec &T9){ - vec 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 -struct Jacobian { - - void operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) - { - // EOS - vec T9=get_T9_array(y[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); - - double pfrac=n15pg_frac(T9); - double afrac=1-pfrac; - - double yh1 = y[ ih1]; - double yhe3 = y[ ihe3]; - double yhe4 = y[ ihe4]; - double yc12 = y[ ic12]; - double yn14 = y[ in14]; - double yo16 = y[ io16]; - double yne20 = y[ine20]; - - // zero all elements to begin - for (int i=0; i < nvar; i++) { - for (int j=0; j < nvar; j++) { J(i,j)=0.0; }} - - // h1 jacobian elements - J(ih1,ih1) = -3*yh1*rpp - 2*yc12*rc12p -2*yn14*rn14p -2*yo16*ro16p; - J(ih1,ihe3) = 2*yhe3*r33 - yhe4*r34; - J(ih1,ihe4) = -yhe3*r34; - J(ih1,ic12) = -2*yh1*rc12p; - J(ih1,in14) = -2*yh1*rn14p; - J(ih1,io16) = -2*yh1*ro16p; - - // he3 jacobian elements - J(ihe3,ih1) = yh1*rpp; - J(ihe3,ihe3) = -2*yhe3*r33 - yhe4*r34; - J(ihe3,ihe4) = -yhe3*r34; - - // he4 jacobian elements - J(ihe4,ih1) = yn14*afrac*rn14p + yo16*ro16p; - J(ihe4,ihe3) = yhe3*r33 - yhe4*r34; - J(ihe4,ihe4) = yhe3*r34 - 1.5*yhe4*yhe4*r3a - yc12*rc12a - 1.5*yn14*rn14a - yo16*ro16a - yne20*rne20a; - J(ihe4,ic12) = -yhe4*rc12a + yc12*r1212 + yo16*r1216; - J(ihe4,in14) = yh1*afrac*rn14p - 1.5*yhe4*rn14a; - J(ihe4,io16) = yh1*ro16p - yhe4*ro16a + yc12*r1216; - J(ihe4,ine20) = -yhe4*rne20a; - - // c12 jacobian elements - J(ic12,ih1) = -yc12*rc12p + yn14*afrac*rn14p; - J(ic12,ihe4) = 0.5*yhe4*yhe4*r3a - yhe4*rc12a; - J(ic12,ic12) = -yh1*rc12p - yhe4*rc12a - yo16*r1216 - 2*yc12*r1212; - J(ic12,in14) = yh1*yn14*afrac*rn14p; - J(ic12,io16) = -yc12*r1216; - - // n14 jacobian elements - J(in14,ih1) = yc12*rc12p - yn14*rn14p + yo16*ro16p; - J(in14,ihe4) = -yn14*rn14a; - J(in14,ic12) = yh1*rc12p; - J(in14,in14) = -yh1*rn14p - yhe4*rn14a; - J(in14,io16) = yo16*ro16p; - - // o16 jacobian elements - J(io16,ih1) = yn14*pfrac*rn14p - yo16*ro16p; - J(io16,ihe4) = yc12*rc12a - yo16*ro16a; - J(io16,ic12) = yhe4*rc12a - yo16*r1216; - J(io16,in14) = yh1*pfrac*rn14p; - J(io16,io16) = yh1*ro16p - yc12*r1216 -yhe4*ro16a; - - // ne20 jacobian elements - J(ine20,ihe4) = yn14*rn14a + yo16*ro16a - yne20*rne20a; - J(ine20,ic12) = yc12*r1212; - J(ine20,in14) = yhe4*rn14a; - J(ine20,io16) = yo16*ro16a; - J(ine20,ine20) = -yhe4*rne20a; - - // mg24 jacobian elements - J(img24,ihe4) = yne20*rne20a; - J(img24,ic12) = yo16*r1216; - J(img24,io16) = yc12*r1216; - J(img24,ine20) = yhe4*rne20a; - - // energy accounting - for (int j=0; j>(1.0e-6,1.0e-6) , - std::make_pair( ODE(), Jacobian() ), y, 0.0, tmax, dt0); - - } else { - std::cout << " *** Explicit RK Dormand-Prince *** " << std::endl; - - num_of_steps = integrate_const( make_dense_output>(1.0e-6, 1.0e-6), - ODE(), y, 0.0, tmax, dt0); - - } - - //convert number fraction to mass fraction - ysum=0; - for (int i=0; i < niso; i++) { - y[i] *= aion[i]; - ysum+= y[i]; - } - for (int i=0; i < niso; i++) {y[i] /= ysum;} - - std::cout << " H1: " << y[ih1] << std::endl; - std::cout << " He4: " << y[ihe4] << std::endl; - std::cout << "energy: " << y[iener] << std::endl; - std::cout << " steps: " << num_of_steps << std::endl; - return 0; - } diff --git a/src/network/public/approx8.h b/src/network/public/approx8.h new file mode 100644 index 0000000..22ff778 --- /dev/null +++ b/src/network/public/approx8.h @@ -0,0 +1,316 @@ +#ifndef APPROX8_H +#define APPROX8_H + +#include + +#include +#include +#include + +#include "network.h" + +/** + * @file approx8.h + * @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 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 nnApprox8{ + + using namespace boost::numeric::odeint; + + /** + * @struct Net + * @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; + + const static int itemp=img24+1; + const static int iden =itemp+1; + const static int iener=iden+1; + + const static int niso=img24+1; // number of isotopes + const static int nvar=iener+1; // number of variables + + static constexpr std::array aion = { + 1, + 3, + 4, + 12, + 14, + 16, + 20, + 24 + }; + + static constexpr std::array mion = { + 1.67262164e-24, + 5.00641157e-24, + 6.64465545e-24, + 1.99209977e-23, + 2.32462686e-23, + 2.65528858e-23, + 3.31891077e-23, + 3.98171594e-23 + }; + + }; + + /** + * @brief Multiplies two arrays and sums the resulting elements. + * @param a First array. + * @param b Second array. + * @return Sum of the product of the arrays. + * @example + * @code + * vec7 a = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0}; + * vec7 b = {0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5}; + * double result = sum_product(a, b); + * @endcode + */ + double sum_product( const vec7 &a, const vec7 &b); + + /** + * @brief Returns an array of T9 terms for the nuclear reaction rate fit. + * @param T Temperature in GigaKelvin. + * @return Array of T9 terms. + * @example + * @code + * double T = 1.5; + * vec7 T9_array = get_T9_array(T); + * @endcode + */ + vec7 get_T9_array(const double &T); + + /** + * @brief Evaluates the nuclear reaction rate given the T9 array and coefficients. + * @param T9 Array of T9 terms. + * @param coef Array of coefficients. + * @return Evaluated rate. + * @example + * @code + * vec7 T9 = get_T9_array(1.5); + * vec7 coef = {1.0, 0.1, 0.01, 0.001, 0.0001, 0.00001, 0.000001}; + * double rate = rate_fit(T9, coef); + * @endcode + */ + double rate_fit(const vec7 &T9, const vec7 &coef); + + /** + * @brief Calculates the rate for the reaction p + p -> d. + * @param T9 Array of T9 terms. + * @return Rate of the reaction. + */ + double pp_rate(const vec7 &T9); + + /** + * @brief Calculates the rate for the reaction p + d -> he3. + * @param T9 Array of T9 terms. + * @return Rate of the reaction. + */ + double dp_rate(const vec7 &T9); + + /** + * @brief Calculates the rate for the reaction he3 + he3 -> he4 + 2p. + * @param T9 Array of T9 terms. + * @return Rate of the reaction. + */ + double he3he3_rate(const vec7 &T9); + + /** + * @brief Calculates the rate for the reaction he3(he3,2p)he4. + * @param T9 Array of T9 terms. + * @return Rate of the reaction. + */ + double he3he4_rate(const vec7 &T9); + + /** + * @brief Calculates the rate for the reaction he4 + he4 + he4 -> c12. + * @param T9 Array of T9 terms. + * @return Rate of the reaction. + */ + double triple_alpha_rate(const vec7 &T9); + + /** + * @brief Calculates the rate for the reaction c12 + p -> n13. + * @param T9 Array of T9 terms. + * @return Rate of the reaction. + */ + double c12p_rate(const vec7 &T9); + + /** + * @brief Calculates the rate for the reaction c12 + he4 -> o16. + * @param T9 Array of T9 terms. + * @return Rate of the reaction. + */ + double c12a_rate(const vec7 &T9); + + /** + * @brief Calculates the rate for the reaction n14(p,g)o15 - o15 + p -> c12 + he4. + * @param T9 Array of T9 terms. + * @return Rate of the reaction. + */ + double n14p_rate(const vec7 &T9); + + /** + * @brief Calculates the rate for the reaction n14(a,g)f18 assumed to go on to ne20. + * @param T9 Array of T9 terms. + * @return Rate of the reaction. + */ + double n14a_rate(const vec7 &T9); + + /** + * @brief Calculates the rate for the reaction n15(p,a)c12 (CNO I). + * @param T9 Array of T9 terms. + * @return Rate of the reaction. + */ + double n15pa_rate(const vec7 &T9); + + /** + * @brief Calculates the rate for the reaction n15(p,g)o16 (CNO II). + * @param T9 Array of T9 terms. + * @return Rate of the reaction. + */ + double n15pg_rate(const vec7 &T9); + + /** + * @brief Calculates the fraction for the reaction n15(p,g)o16. + * @param T9 Array of T9 terms. + * @return Fraction of the reaction. + */ + double n15pg_frac(const vec7 &T9); + + /** + * @brief Calculates the rate for the reaction o16(p,g)f17 then f17 -> o17(p,a)n14. + * @param T9 Array of T9 terms. + * @return Rate of the reaction. + */ + double o16p_rate(const vec7 &T9); + + /** + * @brief Calculates the rate for the reaction o16(a,g)ne20. + * @param T9 Array of T9 terms. + * @return Rate of the reaction. + */ + double o16a_rate(const vec7 &T9); + + /** + * @brief Calculates the rate for the reaction ne20(a,g)mg24. + * @param T9 Array of T9 terms. + * @return Rate of the reaction. + */ + double ne20a_rate(const vec7 &T9); + + /** + * @brief Calculates the rate for the reaction c12(c12,a)ne20. + * @param T9 Array of T9 terms. + * @return Rate of the reaction. + */ + double c12c12_rate(const vec7 &T9); + + /** + * @brief Calculates the rate for the reaction c12(o16,a)mg24. + * @param T9 Array of T9 terms. + * @return Rate of the reaction. + */ + double c12o16_rate(const vec7 &T9); + + /** + * @struct Jacobian + * @brief Functor to calculate the Jacobian matrix for implicit solvers. + */ + struct Jacobian { + /** + * @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 ); + }; + + /** + * @struct ODE + * @brief Functor to calculate the derivatives for the ODE solver. + */ + struct ODE { + /** + * @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 */); + }; + + /** + * @class Approx8Network + * @brief Class for the Approx8 nuclear reaction network. + */ + class Approx8Network : public nuclearNetwork::Network { + public: + /** + * @brief Evaluates the nuclear network. + * @param netIn Input parameters for the network. + * @return Output results from the network. + */ + virtual nuclearNetwork::NetOut evaluate(const nuclearNetwork::NetIn &netIn); + + /** + * @brief Sets whether the solver should use a stiff method. + * @param stiff Boolean indicating if a stiff method should be used. + */ + void setStiff(bool stiff); + + /** + * @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; } + private: + vector_type m_y; + double m_tmax; + double m_dt0; + bool m_stiff; + + /** + * @brief Converts the input parameters to the internal state vector. + * @param netIn Input parameters for the network. + * @return Internal state vector. + */ + vector_type convert_netIn(const nuclearNetwork::NetIn &netIn); + }; + +} // namespace nnApprox8 + +#endif \ No newline at end of file diff --git a/src/network/public/network.h b/src/network/public/network.h index 3466aa0..9cb6a3c 100644 --- a/src/network/public/network.h +++ b/src/network/public/network.h @@ -1,526 +1,93 @@ - -#include -#include -#include -#include - -#include -#include -#include - -#include "const.h" - -/* Nuclear reaction network in cgs units based on Frank Timmes' "aprox8". - 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 . - -*/ - - -// using namespace std; -using namespace boost::numeric::odeint; -namespace phoenix = boost::phoenix; - -// these types are required by the Rosenbrock implicit solver -typedef boost::numeric::ublas::vector< double > vector_type; -typedef boost::numeric::ublas::matrix< double > matrix_type; - -// this array is used only in the nuclear reaction rate evaluations -typedef std::array vec; - -// only need a couple of constants -Constants& constants = Constants::getInstance(); -const double avo = constants.get("N_a").value; -const double clight = constants.get("c").value; - -// identify the isotopes used in the network. -const int ih1=0; -const int ihe3=1; -const int ihe4=2; -const int ic12=3; -const int in14=4; -const int io16=5; -const int ine20=6; -const int img24=7; - -// physical variables; this routine currently does not need to call EOS -// since the temperature and density are assumed constant during the burn -const int itemp=img24+1; -const int iden =itemp+1; -const int iener=iden+1; - -// number of isotopes and number of variables -const int niso=img24+1; // number of isotopes -const int nvar=iener+1; // number of variables - -// atomic stuff -std::array aion = {1,3,4,12,14,16,20,24}; -//std::array zion = {1,2,2, 6, 7, 8,10,12}; -//std::array bion = {0,7.71819,28.29603,92.16294,104.65998,127.62093,160.64788,198.25790}; -//nion = aion - zion #neutrons -//mion = nion*mn + zion*mp - bion*mev2gr #mass -std::array mion = {1.67262164e-24, 5.00641157e-24, 6.64465545e-24, 1.99209977e-23, - 2.32462686e-23, 2.65528858e-23, 3.31891077e-23, 3.98171594e-23}; - -//helper functions -// a function to multilpy two arrays and then sum the resulting elements: sum(a*b) -double sum_product( const vec &a, const vec &b){ - if (a.size() != b.size()){ - throw std::runtime_error("Error: array size mismatch in sum_product"); - } - - 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 -vec get_T9_array(const double &T) { - vec arr; - double T9=1e-9*T; - 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 vec &T9, const vec &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 vec &T9) { - vec a1 = {-34.78630, 0,-3.511930, 3.100860, -0.1983140, 1.262510e-2, -1.025170}; - vec 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 vec &T9) { - vec a1 = {7.528980, 0, -3.720800, 0.8717820, 0, 0,-0.6666670}; - vec 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 vec &T9){ - vec 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 vec &T9){ - vec a1 = {1.560990e+01,0.000000e+00,-1.282710e+01,-3.082250e-02,-6.546850e-01,8.963310e-02,-6.666670e-01}; - vec 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 vec &T9){ - vec a1 = {-9.710520e-01,0.000000e+00,-3.706000e+01,2.934930e+01,-1.155070e+02,-1.000000e+01,-1.333330e+00}; - vec a2 = {-1.178840e+01,-1.024460e+00,-2.357000e+01,2.048860e+01,-1.298820e+01,-2.000000e+01,-2.166670e+00}; - vec 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 vec &T9){ - vec a1={1.714820e+01,0.000000e+00,-1.369200e+01,-2.308810e-01,4.443620e+00,-3.158980e+00,-6.666670e-01}; - vec 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 vec &T9){ - vec a1={6.965260e+01,-1.392540e+00,5.891280e+01,-1.482730e+02,9.083240e+00,-5.410410e-01,7.035540e+01}; - vec 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 vec &T9){ - vec a1={1.701000e+01,0.000000e+00,-1.519300e+01,-1.619540e-01,-7.521230e+00,-9.875650e-01,-6.666670e-01}; - vec a2={2.011690e+01,0.000000e+00,-1.519300e+01,-4.639750e+00,9.734580e+00,-9.550510e+00,3.333330e-01}; - vec a3={7.654440e+00,-2.998000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec 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 vec &T9){ - vec a1={2.153390e+01,0.000000e+00,-3.625040e+01,0.000000e+00,0.000000e+00,-5.000000e+00,-6.666670e-01}; - vec a2={1.968380e-01,-5.160340e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec 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 vec &T9){ - vec a1 = {2.747640e+01,0.000000e+00,-1.525300e+01,1.593180e+00,2.447900e+00,-2.197080e+00,-6.666670e-01}; - vec a2 = {-4.873470e+00,-2.021170e+00,0.000000e+00,3.084970e+01,-8.504330e+00,-1.544260e+00,-1.500000e+00}; - vec a3 = {2.089720e+01,-7.406000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec 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 vec &T9){ - vec a1 = {2.001760e+01,0.000000e+00,-1.524000e+01,3.349260e-01,4.590880e+00,-4.784680e+00,-6.666670e-01}; - vec a2 = {6.590560e+00,-2.923150e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec 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 vec &T9){ - double f1=n15pg_rate(T9); - double f2=n15pa_rate(T9); - return f1/(f1+f2); -} - -// o16(p,g)f17 then f17 -> o17(p,a)n14 -double o16p_rate(const vec &T9){ - vec 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 vec &T9){ - vec a1={2.390300e+01,0.000000e+00,-3.972620e+01,-2.107990e-01,4.428790e-01,-7.977530e-02,-6.666670e-01}; - vec a2={3.885710e+00,-1.035850e+01,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec 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 vec &T9){ - vec a1={2.450580e+01,0.000000e+00,-4.625250e+01,5.589010e+00,7.618430e+00,-3.683000e+00,-6.666670e-01}; - vec a2={-3.870550e+01,-2.506050e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec a3={1.983070e+00,-9.220260e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec 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 vec &T9){ - vec 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 vec &T9){ - vec 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 -struct Jacobian { - - void operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) - { - // EOS - vec T9=get_T9_array(y[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); - - double pfrac=n15pg_frac(T9); - double afrac=1-pfrac; - - double yh1 = y[ ih1]; - double yhe3 = y[ ihe3]; - double yhe4 = y[ ihe4]; - double yc12 = y[ ic12]; - double yn14 = y[ in14]; - double yo16 = y[ io16]; - double yne20 = y[ine20]; - - // zero all elements to begin - for (int i=0; i < nvar; i++) { - for (int j=0; j < nvar; j++) { J(i,j)=0.0; }} - - // h1 jacobian elements - J(ih1,ih1) = -3*yh1*rpp - 2*yc12*rc12p -2*yn14*rn14p -2*yo16*ro16p; - J(ih1,ihe3) = 2*yhe3*r33 - yhe4*r34; - J(ih1,ihe4) = -yhe3*r34; - J(ih1,ic12) = -2*yh1*rc12p; - J(ih1,in14) = -2*yh1*rn14p; - J(ih1,io16) = -2*yh1*ro16p; - - // he3 jacobian elements - J(ihe3,ih1) = yh1*rpp; - J(ihe3,ihe3) = -2*yhe3*r33 - yhe4*r34; - J(ihe3,ihe4) = -yhe3*r34; - - // he4 jacobian elements - J(ihe4,ih1) = yn14*afrac*rn14p + yo16*ro16p; - J(ihe4,ihe3) = yhe3*r33 - yhe4*r34; - J(ihe4,ihe4) = yhe3*r34 - 1.5*yhe4*yhe4*r3a - yc12*rc12a - 1.5*yn14*rn14a - yo16*ro16a - yne20*rne20a; - J(ihe4,ic12) = -yhe4*rc12a + yc12*r1212 + yo16*r1216; - J(ihe4,in14) = yh1*afrac*rn14p - 1.5*yhe4*rn14a; - J(ihe4,io16) = yh1*ro16p - yhe4*ro16a + yc12*r1216; - J(ihe4,ine20) = -yhe4*rne20a; - - // c12 jacobian elements - J(ic12,ih1) = -yc12*rc12p + yn14*afrac*rn14p; - J(ic12,ihe4) = 0.5*yhe4*yhe4*r3a - yhe4*rc12a; - J(ic12,ic12) = -yh1*rc12p - yhe4*rc12a - yo16*r1216 - 2*yc12*r1212; - J(ic12,in14) = yh1*yn14*afrac*rn14p; - J(ic12,io16) = -yc12*r1216; - - // n14 jacobian elements - J(in14,ih1) = yc12*rc12p - yn14*rn14p + yo16*ro16p; - J(in14,ihe4) = -yn14*rn14a; - J(in14,ic12) = yh1*rc12p; - J(in14,in14) = -yh1*rn14p - yhe4*rn14a; - J(in14,io16) = yo16*ro16p; - - // o16 jacobian elements - J(io16,ih1) = yn14*pfrac*rn14p - yo16*ro16p; - J(io16,ihe4) = yc12*rc12a - yo16*ro16a; - J(io16,ic12) = yhe4*rc12a - yo16*r1216; - J(io16,in14) = yh1*pfrac*rn14p; - J(io16,io16) = yh1*ro16p - yc12*r1216 -yhe4*ro16a; - - // ne20 jacobian elements - J(ine20,ihe4) = yn14*rn14a + yo16*ro16a - yne20*rne20a; - J(ine20,ic12) = yc12*r1212; - J(ine20,in14) = yhe4*rn14a; - J(ine20,io16) = yo16*ro16a; - J(ine20,ine20) = -yhe4*rne20a; - - // mg24 jacobian elements - J(img24,ihe4) = yne20*rne20a; - J(img24,ic12) = yo16*r1216; - J(img24,io16) = yc12*r1216; - J(img24,ine20) = yhe4*rne20a; - - // energy accounting - for (int j=0; j>(1.0e-6,1.0e-6) , - std::make_pair( ODE(), Jacobian() ), y, 0.0, tmax, dt0); - - } else { - std::cout << " *** Explicit RK Dormand-Prince *** " << std::endl; - - num_of_steps = integrate_const( make_dense_output>(1.0e-6, 1.0e-6), - ODE(), y, 0.0, tmax, dt0); - - } - - //convert number fraction to mass fraction - ysum=0; - for (int i=0; i < niso; i++) { - y[i] *= aion[i]; - ysum+= y[i]; - } - for (int i=0; i < niso; i++) {y[i] /= ysum;} - - std::cout << " H1: " << y[ih1] << std::endl; - std::cout << " He4: " << y[ihe4] << std::endl; - std::cout << "energy: " << y[iener] << std::endl; - std::cout << " steps: " << num_of_steps << std::endl; - return 0; - -} +#ifndef NETWORK_H +#define NETWORK_H + +#include + +#include "probe.h" +#include "config.h" +#include "quill/Logger.h" + +namespace nuclearNetwork { + + /** + * @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; + * netIn.composition = {1.0, 0.0, 0.0}; + * netIn.tmax = 1.0e6; + * netIn.dt0 = 1.0e-3; + * netIn.temperature = 1.0e8; + * netIn.density = 1.0e5; + * netIn.energy = 1.0e12; + * @endcode + */ + struct NetIn { + std::vector 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 + double energy; ///< Energy in ergs + }; + + /** + * @struct NetOut + * @brief Output structure for the network evaluation. + * + * This structure holds the output results from the network evaluation. + * + * Example usage: + * @code + * nuclearNetwork::NetOut netOut = network.evaluate(netIn); + * std::vector composition = netOut.composition; + * int steps = netOut.num_steps; + * double energy = netOut.energy; + * @endcode + */ + struct NetOut { + std::vector composition; ///< Composition of the network after evaluation + int num_steps; ///< Number of steps taken in the evaluation + double energy; ///< Energy in ergs after evaluation + }; + + /** + * @class Network + * @brief Class for network evaluation. + * + * This class provides methods to evaluate the network based on the input parameters. + * + * Example usage: + * @code + * nuclearNetwork::Network network; + * nuclearNetwork::NetIn netIn; + * // Set netIn parameters... + * nuclearNetwork::NetOut netOut = network.evaluate(netIn); + * @endcode + */ + class Network { + public: + Network(); + virtual ~Network() = default; + + /** + * @brief Evaluate the network based on the input parameters. + * + * @param netIn Input parameters for the network evaluation. + * @return NetOut Output results from the network evaluation. + */ + virtual NetOut evaluate(const NetIn &netIn); + + protected: + Config& m_config; ///< Configuration instance + Probe::LogManager& m_logManager; ///< Log manager instance + quill::Logger* m_logger; ///< Logger instance + }; + +} // namespace nuclearNetwork + +#endif // NETWORK_H \ No newline at end of file From 8fa950835d9e09ab13f64b02dd0a4b4a1a77511e Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 21 Mar 2025 14:03:46 -0400 Subject: [PATCH 15/18] test(network): added approx8 network tests --- tests/meson.build | 1 + tests/network/approx8Test.cpp | 50 +++++++++++++++++++++++++++++++++++ tests/network/meson.build | 25 ++++++++++++++++++ 3 files changed, 76 insertions(+) create mode 100644 tests/network/approx8Test.cpp create mode 100644 tests/network/meson.build diff --git a/tests/meson.build b/tests/meson.build index 0d89e07..0c045fa 100644 --- a/tests/meson.build +++ b/tests/meson.build @@ -12,6 +12,7 @@ subdir('config') subdir('probe') subdir('eos') subdir('resource') +subdir('network') # Subdirectories for sandbox tests subdir('dobj_sandbox') diff --git a/tests/network/approx8Test.cpp b/tests/network/approx8Test.cpp new file mode 100644 index 0000000..a944c1f --- /dev/null +++ b/tests/network/approx8Test.cpp @@ -0,0 +1,50 @@ +#include +#include + +#include "approx8.h" +#include "config.h" +#include "network.h" + +#include + +std::string TEST_CONFIG = std::string(getenv("MESON_SOURCE_ROOT")) + "/tests/testsConfig.yaml"; +class approx8Test : public ::testing::Test {}; + +/** + * @brief Test the constructor of the Config class. + */ +TEST_F(approx8Test, constructor) { + Config& config = Config::getInstance(); + config.loadConfig(TEST_CONFIG); + EXPECT_NO_THROW(nnApprox8::Approx8Network()); +} + +TEST_F(approx8Test, setStiff) { + nnApprox8::Approx8Network network; + EXPECT_NO_THROW(network.setStiff(true)); + EXPECT_TRUE(network.isStiff()); + EXPECT_NO_THROW(network.setStiff(false)); + EXPECT_FALSE(network.isStiff()); +} + +TEST_F(approx8Test, evaluate) { + nnApprox8::Approx8Network network; + nuclearNetwork::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}; + + netIn.composition = comp; + netIn.temperature = 1e7; + netIn.density = 1e2; + netIn.energy = 0.0; + + netIn.tmax = 3.15e17; + netIn.dt0 = 1e12; + + nuclearNetwork::NetOut netOut; + EXPECT_NO_THROW(netOut = network.evaluate(netIn)); + + EXPECT_DOUBLE_EQ(netOut.composition[nnApprox8::Net::ih1], 0.50166260916650918); + EXPECT_DOUBLE_EQ(netOut.composition[nnApprox8::Net::ihe4],0.48172270591286032); + EXPECT_DOUBLE_EQ(netOut.energy, -1.6433049870528356e+18); +} diff --git a/tests/network/meson.build b/tests/network/meson.build new file mode 100644 index 0000000..5e316fd --- /dev/null +++ b/tests/network/meson.build @@ -0,0 +1,25 @@ +# Test files for network +test_sources = [ + 'approx8Test.cpp', +] + +foreach test_file : test_sources + exe_name = test_file.split('.')[0] + message('Building test: ' + exe_name) + + # Create an executable target for each test + test_exe = executable( + exe_name, + test_file, + dependencies: [gtest_dep, network_dep, gtest_main], + 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 + ) + + # Add the executable as a test + test( + exe_name, + test_exe, + env: ['MESON_SOURCE_ROOT=' + meson.project_source_root(), 'MESON_BUILD_ROOT=' + meson.project_build_root()]) +endforeach From 3109c198cfdb846242a6da3b54ea32d17c58f6af Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 21 Mar 2025 14:35:18 -0400 Subject: [PATCH 16/18] refactor(network): removed unused phoenix dependency and defauled stiff=false to use explicit solver --- src/network/private/approx8.cpp | 2 -- src/network/public/approx8.h | 4 +--- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/src/network/private/approx8.cpp b/src/network/private/approx8.cpp index ec0abc8..c80e3ac 100644 --- a/src/network/private/approx8.cpp +++ b/src/network/private/approx8.cpp @@ -3,8 +3,6 @@ #include #include -#include -#include #include "const.h" #include "config.h" diff --git a/src/network/public/approx8.h b/src/network/public/approx8.h index 22ff778..75d88e8 100644 --- a/src/network/public/approx8.h +++ b/src/network/public/approx8.h @@ -4,8 +4,6 @@ #include #include -#include -#include #include "network.h" @@ -301,7 +299,7 @@ namespace nnApprox8{ vector_type m_y; double m_tmax; double m_dt0; - bool m_stiff; + bool m_stiff = false; /** * @brief Converts the input parameters to the internal state vector. From 59d4c290ba5384bdfcd470bc705f8395e7073aee Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 21 Mar 2025 14:44:35 -0400 Subject: [PATCH 17/18] fix(network): added missing negative sign in energy calculation --- src/network/private/approx8.cpp | 2 +- tests/network/approx8Test.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/network/private/approx8.cpp b/src/network/private/approx8.cpp index c80e3ac..8d08ff5 100644 --- a/src/network/private/approx8.cpp +++ b/src/network/private/approx8.cpp @@ -421,7 +421,7 @@ namespace nnApprox8{ for (int i=0; i Date: Fri, 21 Mar 2025 14:46:22 -0400 Subject: [PATCH 18/18] fix(network): fixed missing negative in jacobian energy accounting --- src/network/private/approx8.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/network/private/approx8.cpp b/src/network/private/approx8.cpp index 8d08ff5..c9e597a 100644 --- a/src/network/private/approx8.cpp +++ b/src/network/private/approx8.cpp @@ -325,7 +325,7 @@ namespace nnApprox8{ for (int i=0; i