feat(GridFire): SERiF now uses GridFire
This commit is contained in:
4
.gitignore
vendored
4
.gitignore
vendored
@@ -72,8 +72,12 @@ subprojects/libconstants/
|
||||
subprojects/liblogging/
|
||||
subprojects/libconfig/
|
||||
subprojects/libcomposition/
|
||||
subprojects/GridFire/
|
||||
|
||||
qhull.wrap
|
||||
quill.wrap
|
||||
yaml-cpp.wrap
|
||||
cppad.wrap
|
||||
|
||||
.vscode/
|
||||
|
||||
|
||||
2
build-config/GridFire/meson.build
Normal file
2
build-config/GridFire/meson.build
Normal file
@@ -0,0 +1,2 @@
|
||||
gridfire_p = subproject('GridFire')
|
||||
gridfire_dep = gridfire_p.get_variable('network_dep')
|
||||
@@ -1,9 +1,9 @@
|
||||
cmake = import('cmake')
|
||||
|
||||
subdir('fourdst')
|
||||
subdir('GridFire')
|
||||
|
||||
subdir('mfem')
|
||||
subdir('yaml-cpp')
|
||||
subdir('boost')
|
||||
subdir('opatIO')
|
||||
subdir('mpi')
|
||||
|
||||
@@ -1,10 +0,0 @@
|
||||
quill_cmake_options = cmake.subproject_options()
|
||||
quill_cmake_options.add_cmake_defines({
|
||||
'BUILD_SHARED_LIBS': 'ON',
|
||||
'CMAKE_SKIP_INSTALL_RULES': 'ON'
|
||||
})
|
||||
quill_sp = cmake.subproject(
|
||||
'quill',
|
||||
options: quill_cmake_options,
|
||||
)
|
||||
quill_dep = quill_sp.dependency('quill')
|
||||
@@ -1,12 +0,0 @@
|
||||
yaml_cpp_cmake_options = cmake.subproject_options()
|
||||
yaml_cpp_cmake_options.add_cmake_defines({
|
||||
'CMAKE_POLICY_VERSION_MINIMUM': '3.5',
|
||||
'BUILD_SHARED_LIBS': 'ON',
|
||||
'CMAKE_SKIP_INSTALL_RULES': 'ON',
|
||||
'YAML_CPP_BUILD_TESTS': 'OFF'
|
||||
})
|
||||
yaml_cpp_sp = cmake.subproject(
|
||||
'yaml-cpp',
|
||||
options: yaml_cpp_cmake_options,
|
||||
)
|
||||
yaml_cpp_dep = yaml_cpp_sp.dependency('yaml-cpp')
|
||||
@@ -25,7 +25,7 @@ py_mod = py_installation.extension_module(
|
||||
mfem_dep,
|
||||
polysolver_dep,
|
||||
trampoline_dep,
|
||||
network_dep,
|
||||
gridfire_dep,
|
||||
],
|
||||
cpp_args : ['-UNDEBUG'], # Example: Ensure assertions are enabled if needed
|
||||
install : true,
|
||||
|
||||
@@ -18,7 +18,7 @@
|
||||
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
|
||||
#
|
||||
# *********************************************************************** #
|
||||
project('SERiF', 'cpp', version: '0.0.1a', default_options: ['cpp_std=c++23'], meson_version: '>=1.6.0')
|
||||
project('SERiF', 'cpp', version: 'v0.0.1a', default_options: ['cpp_std=c++23'], meson_version: '>=1.5.0')
|
||||
|
||||
# Add default visibility for all C++ targets
|
||||
add_project_arguments('-fvisibility=default', language: 'cpp')
|
||||
@@ -37,6 +37,8 @@ endif
|
||||
|
||||
# Pass the DATA_DIR definition to the compiler
|
||||
add_project_arguments('-DDATA_DIR=' + data_dir, language : 'cpp')
|
||||
# We disable these because of CppAD
|
||||
add_project_arguments('-Wno-bitwise-instead-of-logical', language: 'cpp')
|
||||
|
||||
# Build external dependencies first so that all the embedded resources are available to the other targets
|
||||
subdir('build-config')
|
||||
|
||||
@@ -15,7 +15,6 @@ subdir('meshIO')
|
||||
subdir('resource')
|
||||
|
||||
# Physics Libraries
|
||||
subdir('network')
|
||||
subdir('polytrope')
|
||||
|
||||
# Python Bindings
|
||||
|
||||
@@ -1,37 +0,0 @@
|
||||
# Define the library
|
||||
network_sources = files(
|
||||
'private/network.cpp',
|
||||
'private/approx8.cpp'
|
||||
)
|
||||
|
||||
network_headers = files(
|
||||
'public/network.h',
|
||||
'public/approx8.h'
|
||||
)
|
||||
|
||||
dependencies = [
|
||||
boost_dep,
|
||||
const_dep,
|
||||
mfem_dep,
|
||||
config_dep,
|
||||
probe_dep,
|
||||
composition_dep,
|
||||
logging_dep
|
||||
]
|
||||
|
||||
# Define the libnetwork library so it can be linked against by other parts of the build system
|
||||
libnetwork = 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')
|
||||
@@ -1,542 +0,0 @@
|
||||
/* ***********************************************************************
|
||||
//
|
||||
// Copyright (C) 2025 -- The 4D-STAR Collaboration
|
||||
// File Author: Emily Boudreaux
|
||||
// Last Modified: March 21, 2025
|
||||
//
|
||||
// 4DSSE is free software; you can use it and/or modify
|
||||
// it under the terms and restrictions the GNU General Library Public
|
||||
// License version 3 (GPLv3) as published by the Free Software Foundation.
|
||||
//
|
||||
// 4DSSE is distributed in the hope that it will be useful,
|
||||
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
|
||||
// See the GNU Library General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU Library General Public License
|
||||
// along with this software; if not, write to the Free Software
|
||||
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
|
||||
//
|
||||
// *********************************************************************** */
|
||||
#include <cmath>
|
||||
#include <stdexcept>
|
||||
#include <array>
|
||||
|
||||
#include <boost/numeric/odeint.hpp>
|
||||
|
||||
#include "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' "approx8".
|
||||
At this time it does neither screening nor neutrino losses. It includes
|
||||
the following 8 isotopes:
|
||||
|
||||
h1
|
||||
he3
|
||||
he4
|
||||
c12
|
||||
n14
|
||||
o16
|
||||
ne20
|
||||
mg24
|
||||
|
||||
and the following nuclear reactions:
|
||||
|
||||
---pp chain---
|
||||
p(p,e+)d
|
||||
d(p,g)he3
|
||||
he3(he3,2p)he4
|
||||
|
||||
---CNO cycle---
|
||||
c12(p,g)n13 - n13 -> c13 + p -> n14
|
||||
n14(p,g)o15 - o15 + p -> c12 + he4
|
||||
n14(a,g)f18 - proceeds to ne20
|
||||
n15(p,a)c12 - / these two n15 reactions are \ CNO I
|
||||
n15(p,g)o16 - \ used to calculate a fraction / CNO II
|
||||
o16(p,g)f17 - f17 + e -> o17 and then o17 + p -> n14 + he4
|
||||
|
||||
---alpha captures---
|
||||
c12(a,g)o16
|
||||
triple alpha
|
||||
o16(a,g)ne20
|
||||
ne20(a,g)mg24
|
||||
c12(c12,a)ne20
|
||||
c12(o16,a)mg24
|
||||
|
||||
At present the rates are all evaluated using a fitting function.
|
||||
The coefficients to the fit are from reaclib.jinaweb.org .
|
||||
|
||||
*/
|
||||
|
||||
namespace serif::network::approx8{
|
||||
|
||||
// using namespace std;
|
||||
using namespace boost::numeric::odeint;
|
||||
|
||||
//helper functions
|
||||
// a function to multiply two arrays and then sum the resulting elements: sum(a*b)
|
||||
double sum_product( const vec7 &a, const vec7 &b){
|
||||
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;
|
||||
const double T9=1e-9*T;
|
||||
const double T913=pow(T9,1./3.);
|
||||
|
||||
arr[0]=1;
|
||||
arr[1]=1/T9;
|
||||
arr[2]=1/T913;
|
||||
arr[3]=T913;
|
||||
arr[4]=T9;
|
||||
arr[5]=pow(T9,5./3.);
|
||||
arr[6]=log(T9);
|
||||
|
||||
return arr;
|
||||
}
|
||||
|
||||
// this function uses the two preceding functions to evaluate the rate given the T9 array and the coefficients
|
||||
double rate_fit(const vec7 &T9, const vec7 &coef){
|
||||
return exp(sum_product(T9,coef));
|
||||
}
|
||||
|
||||
// p + p -> d; this, like some of the other rates, this is a composite of multiple fits
|
||||
double pp_rate(const vec7 &T9) {
|
||||
constexpr vec7 a1 = {-34.78630, 0,-3.511930, 3.100860, -0.1983140, 1.262510e-2, -1.025170};
|
||||
constexpr vec7 a2 = { -4.364990e+1,-2.460640e-3,-2.750700,-4.248770e-1,1.598700e-2,-6.908750e-4,-2.076250e-1};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2);
|
||||
}
|
||||
|
||||
// p + d -> he3
|
||||
double dp_rate(const vec7 &T9) {
|
||||
constexpr vec7 a1 = {7.528980, 0, -3.720800, 0.8717820, 0, 0,-0.6666670};
|
||||
constexpr vec7 a2 = {8.935250, 0, -3.720800, 0.1986540, 0, 0, 0.3333330};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2);
|
||||
}
|
||||
|
||||
// he3 + he3 -> he4 + 2p
|
||||
double he3he3_rate(const vec7 &T9){
|
||||
constexpr vec7 a = {2.477880e+01,0,-12.27700,-0.1036990,-6.499670e-02,1.681910e-02,-6.666670e-01};
|
||||
return rate_fit(T9,a);
|
||||
}
|
||||
|
||||
// he3(he3,2p)he4
|
||||
double he3he4_rate(const vec7 &T9){
|
||||
constexpr vec7 a1 = {1.560990e+01,0.000000e+00,-1.282710e+01,-3.082250e-02,-6.546850e-01,8.963310e-02,-6.666670e-01};
|
||||
constexpr vec7 a2 = {1.770750e+01,0.000000e+00,-1.282710e+01,-3.812600e+00,9.422850e-02,-3.010180e-03,1.333330e+00};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2);
|
||||
}
|
||||
|
||||
// he4 + he4 + he4 -> c12
|
||||
double triple_alpha_rate(const vec7 &T9){
|
||||
constexpr vec7 a1 = {-9.710520e-01,0.000000e+00,-3.706000e+01,2.934930e+01,-1.155070e+02,-1.000000e+01,-1.333330e+00};
|
||||
constexpr vec7 a2 = {-1.178840e+01,-1.024460e+00,-2.357000e+01,2.048860e+01,-1.298820e+01,-2.000000e+01,-2.166670e+00};
|
||||
constexpr vec7 a3 = {-2.435050e+01,-4.126560e+00,-1.349000e+01,2.142590e+01,-1.347690e+00,8.798160e-02,-1.316530e+01};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
|
||||
}
|
||||
|
||||
// c12 + p -> n13
|
||||
double c12p_rate(const vec7 &T9){
|
||||
constexpr vec7 a1={1.714820e+01,0.000000e+00,-1.369200e+01,-2.308810e-01,4.443620e+00,-3.158980e+00,-6.666670e-01};
|
||||
constexpr vec7 a2={1.754280e+01,-3.778490e+00,-5.107350e+00,-2.241110e+00,1.488830e-01,0.000000e+00,-1.500000e+00};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2);
|
||||
}
|
||||
|
||||
// c12 + he4 -> o16
|
||||
double c12a_rate(const vec7 &T9){
|
||||
constexpr vec7 a1={6.965260e+01,-1.392540e+00,5.891280e+01,-1.482730e+02,9.083240e+00,-5.410410e-01,7.035540e+01};
|
||||
constexpr vec7 a2={2.546340e+02,-1.840970e+00,1.034110e+02,-4.205670e+02,6.408740e+01,-1.246240e+01,1.373030e+02};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2);
|
||||
}
|
||||
|
||||
// n14(p,g)o15 - o15 + p -> c12 + he4
|
||||
double n14p_rate(const vec7 &T9){
|
||||
constexpr vec7 a1={1.701000e+01,0.000000e+00,-1.519300e+01,-1.619540e-01,-7.521230e+00,-9.875650e-01,-6.666670e-01};
|
||||
constexpr vec7 a2={2.011690e+01,0.000000e+00,-1.519300e+01,-4.639750e+00,9.734580e+00,-9.550510e+00,3.333330e-01};
|
||||
constexpr vec7 a3={7.654440e+00,-2.998000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
|
||||
constexpr vec7 a4={6.735780e+00,-4.891000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,6.820000e-02};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3) + rate_fit(T9,a4);
|
||||
}
|
||||
|
||||
// n14(a,g)f18 assumed to go on to ne20
|
||||
double n14a_rate(const vec7 &T9){
|
||||
constexpr vec7 a1={2.153390e+01,0.000000e+00,-3.625040e+01,0.000000e+00,0.000000e+00,-5.000000e+00,-6.666670e-01};
|
||||
constexpr vec7 a2={1.968380e-01,-5.160340e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
|
||||
constexpr vec7 a3={1.389950e+01,-1.096560e+01,-5.622700e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
|
||||
}
|
||||
|
||||
// n15(p,a)c12 (CNO I)
|
||||
double n15pa_rate(const vec7 &T9){
|
||||
constexpr vec7 a1 = {2.747640e+01,0.000000e+00,-1.525300e+01,1.593180e+00,2.447900e+00,-2.197080e+00,-6.666670e-01};
|
||||
constexpr vec7 a2 = {-4.873470e+00,-2.021170e+00,0.000000e+00,3.084970e+01,-8.504330e+00,-1.544260e+00,-1.500000e+00};
|
||||
constexpr vec7 a3 = {2.089720e+01,-7.406000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
|
||||
constexpr vec7 a4 = {-6.575220e+00,-1.163800e+00,0.000000e+00,2.271050e+01,-2.907070e+00,2.057540e-01,-1.500000e+00};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3) + rate_fit(T9,a4);
|
||||
}
|
||||
|
||||
// n15(p,g)o16 (CNO II)
|
||||
double n15pg_rate(const vec7 &T9){
|
||||
constexpr vec7 a1 = {2.001760e+01,0.000000e+00,-1.524000e+01,3.349260e-01,4.590880e+00,-4.784680e+00,-6.666670e-01};
|
||||
constexpr vec7 a2 = {6.590560e+00,-2.923150e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
|
||||
constexpr vec7 a3 = {1.454440e+01,-1.022950e+01,0.000000e+00,0.000000e+00,4.590370e-02,0.000000e+00,-1.500000e+00};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
|
||||
}
|
||||
|
||||
double n15pg_frac(const vec7 &T9){
|
||||
const double f1=n15pg_rate(T9);
|
||||
const double f2=n15pa_rate(T9);
|
||||
return f1/(f1+f2);
|
||||
}
|
||||
|
||||
// o16(p,g)f17 then f17 -> o17(p,a)n14
|
||||
double o16p_rate(const vec7 &T9){
|
||||
constexpr vec7 a={1.909040e+01,0.000000e+00,-1.669600e+01,-1.162520e+00,2.677030e-01,-3.384110e-02,-6.666670e-01};
|
||||
return rate_fit(T9,a);
|
||||
}
|
||||
|
||||
// o16(a,g)ne20
|
||||
double o16a_rate(const vec7 &T9){
|
||||
constexpr vec7 a1={2.390300e+01,0.000000e+00,-3.972620e+01,-2.107990e-01,4.428790e-01,-7.977530e-02,-6.666670e-01};
|
||||
constexpr vec7 a2={3.885710e+00,-1.035850e+01,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
|
||||
constexpr vec7 a3={9.508480e+00,-1.276430e+01,0.000000e+00,-3.659250e+00,7.142240e-01,-1.075080e-03,-1.500000e+00};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
|
||||
}
|
||||
|
||||
// ne20(a,g)mg24
|
||||
double ne20a_rate(const vec7 &T9){
|
||||
constexpr vec7 a1={2.450580e+01,0.000000e+00,-4.625250e+01,5.589010e+00,7.618430e+00,-3.683000e+00,-6.666670e-01};
|
||||
constexpr vec7 a2={-3.870550e+01,-2.506050e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
|
||||
constexpr vec7 a3={1.983070e+00,-9.220260e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
|
||||
constexpr vec7 a4={-8.798270e+00,-1.278090e+01,0.000000e+00,1.692290e+01,-2.573250e+00,2.089970e-01,-1.500000e+00};
|
||||
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3) + rate_fit(T9,a4);
|
||||
}
|
||||
|
||||
// c12(c12,a)ne20
|
||||
double c12c12_rate(const vec7 &T9){
|
||||
constexpr vec7 a={6.128630e+01,0.000000e+00,-8.416500e+01,-1.566270e+00,-7.360840e-02,-7.279700e-02,-6.666670e-01};
|
||||
return rate_fit(T9,a);
|
||||
}
|
||||
|
||||
// c12(o16,a)mg24
|
||||
double c12o16_rate(const vec7 &T9){
|
||||
constexpr vec7 a={4.853410e+01,3.720400e-01,-1.334130e+02,5.015720e+01,-3.159870e+00,1.782510e-02,-2.370270e+01};
|
||||
return rate_fit(T9,a);
|
||||
}
|
||||
|
||||
|
||||
// for Boost ODE solvers either a struct or a class is required
|
||||
|
||||
// a Jacobian matrix for implicit solvers
|
||||
|
||||
void Jacobian::operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) const {
|
||||
fourdst::constant::Constants& constants = fourdst::constant::Constants::getInstance();
|
||||
const double avo = constants.get("N_a").value;
|
||||
const double clight = constants.get("c").value;
|
||||
// EOS
|
||||
const vec7 T9=get_T9_array(y[Approx8Net::iTemp]);
|
||||
|
||||
// evaluate rates once per call
|
||||
const double rpp=pp_rate(T9);
|
||||
const double r33=he3he3_rate(T9);
|
||||
const double r34=he3he4_rate(T9);
|
||||
const double r3a=triple_alpha_rate(T9);
|
||||
const double rc12p=c12p_rate(T9);
|
||||
const double rc12a=c12a_rate(T9);
|
||||
const double rn14p=n14p_rate(T9);
|
||||
const double rn14a=n14a_rate(T9);
|
||||
const double ro16p=o16p_rate(T9);
|
||||
const double ro16a=o16a_rate(T9);
|
||||
const double rne20a=ne20a_rate(T9);
|
||||
const double r1212=c12c12_rate(T9);
|
||||
const double r1216=c12o16_rate(T9);
|
||||
|
||||
const double pFrac=n15pg_frac(T9);
|
||||
const double aFrac=1-pFrac;
|
||||
|
||||
const double yh1 = y[Approx8Net::ih1];
|
||||
const double yhe3 = y[Approx8Net::ihe3];
|
||||
const double yhe4 = y[Approx8Net::ihe4];
|
||||
const double yc12 = y[Approx8Net::ic12];
|
||||
const double yn14 = y[Approx8Net::in14];
|
||||
const double yo16 = y[Approx8Net::io16];
|
||||
const double yne20 = y[Approx8Net::ine20];
|
||||
|
||||
// zero all elements to begin
|
||||
for (int i=0; i < Approx8Net::nVar; i++) {
|
||||
for (int j=0; j < Approx8Net::nVar; j++) {
|
||||
J(i,j)=0.0;
|
||||
}
|
||||
}
|
||||
|
||||
// h1 jacobian elements
|
||||
J(Approx8Net::ih1,Approx8Net::ih1) = -3*yh1*rpp - 2*yc12*rc12p -2*yn14*rn14p -2*yo16*ro16p;
|
||||
J(Approx8Net::ih1,Approx8Net::ihe3) = 2*yhe3*r33 - yhe4*r34;
|
||||
J(Approx8Net::ih1,Approx8Net::ihe4) = -yhe3*r34;
|
||||
J(Approx8Net::ih1,Approx8Net::ic12) = -2*yh1*rc12p;
|
||||
J(Approx8Net::ih1,Approx8Net::in14) = -2*yh1*rn14p;
|
||||
J(Approx8Net::ih1,Approx8Net::io16) = -2*yh1*ro16p;
|
||||
|
||||
// he3 jacobian elements
|
||||
J(Approx8Net::ihe3,Approx8Net::ih1) = yh1*rpp;
|
||||
J(Approx8Net::ihe3,Approx8Net::ihe3) = -2*yhe3*r33 - yhe4*r34;
|
||||
J(Approx8Net::ihe3,Approx8Net::ihe4) = -yhe3*r34;
|
||||
|
||||
// he4 jacobian elements
|
||||
J(Approx8Net::ihe4,Approx8Net::ih1) = yn14*aFrac*rn14p + yo16*ro16p;
|
||||
J(Approx8Net::ihe4,Approx8Net::ihe3) = yhe3*r33 - yhe4*r34;
|
||||
J(Approx8Net::ihe4,Approx8Net::ihe4) = yhe3*r34 - 1.5*yhe4*yhe4*r3a - yc12*rc12a - 1.5*yn14*rn14a - yo16*ro16a - yne20*rne20a;
|
||||
J(Approx8Net::ihe4,Approx8Net::ic12) = -yhe4*rc12a + yc12*r1212 + yo16*r1216;
|
||||
J(Approx8Net::ihe4,Approx8Net::in14) = yh1*aFrac*rn14p - 1.5*yhe4*rn14a;
|
||||
J(Approx8Net::ihe4,Approx8Net::io16) = yh1*ro16p - yhe4*ro16a + yc12*r1216;
|
||||
J(Approx8Net::ihe4,Approx8Net::ine20) = -yhe4*rne20a;
|
||||
|
||||
// c12 jacobian elements
|
||||
J(Approx8Net::ic12,Approx8Net::ih1) = -yc12*rc12p + yn14*aFrac*rn14p;
|
||||
J(Approx8Net::ic12,Approx8Net::ihe4) = 0.5*yhe4*yhe4*r3a - yhe4*rc12a;
|
||||
J(Approx8Net::ic12,Approx8Net::ic12) = -yh1*rc12p - yhe4*rc12a - yo16*r1216 - 2*yc12*r1212;
|
||||
J(Approx8Net::ic12,Approx8Net::in14) = yh1*yn14*aFrac*rn14p;
|
||||
J(Approx8Net::ic12,Approx8Net::io16) = -yc12*r1216;
|
||||
|
||||
// n14 jacobian elements
|
||||
J(Approx8Net::in14,Approx8Net::ih1) = yc12*rc12p - yn14*rn14p + yo16*ro16p;
|
||||
J(Approx8Net::in14,Approx8Net::ihe4) = -yn14*rn14a;
|
||||
J(Approx8Net::in14,Approx8Net::ic12) = yh1*rc12p;
|
||||
J(Approx8Net::in14,Approx8Net::in14) = -yh1*rn14p - yhe4*rn14a;
|
||||
J(Approx8Net::in14,Approx8Net::io16) = yo16*ro16p;
|
||||
|
||||
// o16 jacobian elements
|
||||
J(Approx8Net::io16,Approx8Net::ih1) = yn14*pFrac*rn14p - yo16*ro16p;
|
||||
J(Approx8Net::io16,Approx8Net::ihe4) = yc12*rc12a - yo16*ro16a;
|
||||
J(Approx8Net::io16,Approx8Net::ic12) = yhe4*rc12a - yo16*r1216;
|
||||
J(Approx8Net::io16,Approx8Net::in14) = yh1*pFrac*rn14p;
|
||||
J(Approx8Net::io16,Approx8Net::io16) = yh1*ro16p - yc12*r1216 -yhe4*ro16a;
|
||||
|
||||
// ne20 jacobian elements
|
||||
J(Approx8Net::ine20,Approx8Net::ihe4) = yn14*rn14a + yo16*ro16a - yne20*rne20a;
|
||||
J(Approx8Net::ine20,Approx8Net::ic12) = yc12*r1212;
|
||||
J(Approx8Net::ine20,Approx8Net::in14) = yhe4*rn14a;
|
||||
J(Approx8Net::ine20,Approx8Net::io16) = yo16*ro16a;
|
||||
J(Approx8Net::ine20,Approx8Net::ine20) = -yhe4*rne20a;
|
||||
|
||||
// mg24 jacobian elements
|
||||
J(Approx8Net::img24,Approx8Net::ihe4) = yne20*rne20a;
|
||||
J(Approx8Net::img24,Approx8Net::ic12) = yo16*r1216;
|
||||
J(Approx8Net::img24,Approx8Net::io16) = yc12*r1216;
|
||||
J(Approx8Net::img24,Approx8Net::ine20) = yhe4*rne20a;
|
||||
|
||||
// energy accounting
|
||||
for (int j=0; j<Approx8Net::nIso; j++) {
|
||||
for (int i=0; i<Approx8Net::nIso; i++) {
|
||||
J(Approx8Net::iEnergy,j) += J(i,j)*Approx8Net::mIon[i];
|
||||
}
|
||||
J(Approx8Net::iEnergy,j) *= -avo*clight*clight;
|
||||
}
|
||||
}
|
||||
|
||||
void ODE::operator() ( const vector_type &y, vector_type &dydt, double /* t */) const {
|
||||
const fourdst::constant::Constants& constants = fourdst::constant::Constants::getInstance();
|
||||
const double avo = constants.get("N_a").value;
|
||||
const double clight = constants.get("c").value;
|
||||
|
||||
// EOS
|
||||
const double T = y[Approx8Net::iTemp];
|
||||
const double den = y[Approx8Net::iDensity];
|
||||
const vec7 T9=get_T9_array(T);
|
||||
|
||||
// rates
|
||||
const double rpp=den*pp_rate(T9);
|
||||
const double r33=den*he3he3_rate(T9);
|
||||
const double r34=den*he3he4_rate(T9);
|
||||
const double r3a=den*den*triple_alpha_rate(T9);
|
||||
const double rc12p=den*c12p_rate(T9);
|
||||
const double rc12a=den*c12a_rate(T9);
|
||||
const double rn14p=den*n14p_rate(T9);
|
||||
const double rn14a=n14a_rate(T9);
|
||||
const double ro16p=den*o16p_rate(T9);
|
||||
const double ro16a=den*o16a_rate(T9);
|
||||
const double rne20a=den*ne20a_rate(T9);
|
||||
const double r1212=den*c12c12_rate(T9);
|
||||
const double r1216=den*c12o16_rate(T9);
|
||||
|
||||
const double pFrac=n15pg_frac(T9);
|
||||
const double aFrac=1-pFrac;
|
||||
|
||||
const double yh1 = y[Approx8Net:: ih1];
|
||||
const double yhe3 = y[Approx8Net:: ihe3];
|
||||
const double yhe4 = y[Approx8Net:: ihe4];
|
||||
const double yc12 = y[Approx8Net:: ic12];
|
||||
const double yn14 = y[Approx8Net:: in14];
|
||||
const double yo16 = y[Approx8Net:: io16];
|
||||
const double yne20 = y[Approx8Net::ine20];
|
||||
|
||||
dydt[Approx8Net::ih1] = -1.5*yh1*yh1*rpp;
|
||||
dydt[Approx8Net::ih1] += yhe3*yhe3*r33;
|
||||
dydt[Approx8Net::ih1] += -yhe3*yhe4*r34;
|
||||
dydt[Approx8Net::ih1] += -2*yh1*yc12*rc12p;
|
||||
dydt[Approx8Net::ih1] += -2*yh1*yn14*rn14p;
|
||||
dydt[Approx8Net::ih1] += -2*yh1*yo16*ro16p;
|
||||
|
||||
dydt[Approx8Net::ihe3] = 0.5*yh1*yh1*rpp;
|
||||
dydt[Approx8Net::ihe3] += -yhe3*yhe3*r33;
|
||||
dydt[Approx8Net::ihe3] += -yhe3*yhe4*r34;
|
||||
|
||||
dydt[Approx8Net::ihe4] = 0.5*yhe3*yhe3*r33;
|
||||
dydt[Approx8Net::ihe4] += yhe3*yhe4*r34;
|
||||
dydt[Approx8Net::ihe4] += -yhe4*yc12*rc12a;
|
||||
dydt[Approx8Net::ihe4] += yh1*yn14*aFrac*rn14p;
|
||||
dydt[Approx8Net::ihe4] += yh1*yo16*ro16p;
|
||||
dydt[Approx8Net::ihe4] += -0.5*yhe4*yhe4*yhe4*r3a;
|
||||
dydt[Approx8Net::ihe4] += -yhe4*yo16*ro16a;
|
||||
dydt[Approx8Net::ihe4] += 0.5*yc12*yc12*r1212;
|
||||
dydt[Approx8Net::ihe4] += yc12*yo16*r1216;
|
||||
dydt[Approx8Net::ihe4] += -yhe4*yne20*rne20a;
|
||||
|
||||
dydt[Approx8Net::ic12] = (1./6.)*yhe4*yhe4*yhe4*r3a;
|
||||
dydt[Approx8Net::ic12] += -yhe4*yc12*rc12a;
|
||||
dydt[Approx8Net::ic12] += -yh1*yc12*rc12p;
|
||||
dydt[Approx8Net::ic12] += yh1*yn14*aFrac*rn14p;
|
||||
dydt[Approx8Net::ic12] += -yc12*yc12*r1212;
|
||||
dydt[Approx8Net::ic12] += -yc12*yo16*r1216;
|
||||
|
||||
dydt[Approx8Net::in14] = yh1*yc12*rc12p;
|
||||
dydt[Approx8Net::in14] += -yh1*yn14*rn14p;
|
||||
dydt[Approx8Net::in14] += yh1*yo16*ro16p;
|
||||
dydt[Approx8Net::in14] += -yhe4*yn14*rn14a;
|
||||
|
||||
dydt[Approx8Net::io16] = yhe4*yc12*rc12a;
|
||||
dydt[Approx8Net::io16] += yh1*yn14*pFrac*rn14p;
|
||||
dydt[Approx8Net::io16] += -yh1*yo16*ro16p;
|
||||
dydt[Approx8Net::io16] += -yc12*yo16*r1216;
|
||||
dydt[Approx8Net::io16] += -yhe4*yo16*ro16a;
|
||||
|
||||
dydt[Approx8Net::ine20] = 0.5*yc12*yc12*r1212;
|
||||
dydt[Approx8Net::ine20] += yhe4*yn14*rn14a;
|
||||
dydt[Approx8Net::ine20] += yhe4*yo16*ro16a;
|
||||
dydt[Approx8Net::ine20] += -yhe4*yne20*rne20a;
|
||||
|
||||
dydt[Approx8Net::img24] = yc12*yo16*r1216;
|
||||
dydt[Approx8Net::img24] += yhe4*yne20*rne20a;
|
||||
|
||||
dydt[Approx8Net::iTemp] = 0.;
|
||||
dydt[Approx8Net::iDensity] = 0.;
|
||||
|
||||
// energy accounting
|
||||
double eNuc = 0.;
|
||||
for (int i=0; i<Approx8Net::nIso; i++) {
|
||||
eNuc += Approx8Net::mIon[i]*dydt[i];
|
||||
}
|
||||
dydt[Approx8Net::iEnergy] = -eNuc*avo*clight*clight;
|
||||
}
|
||||
|
||||
Approx8Network::Approx8Network() : Network(APPROX8) {}
|
||||
|
||||
NetOut Approx8Network::evaluate(const NetIn &netIn) {
|
||||
m_y = convert_netIn(netIn);
|
||||
m_tMax = netIn.tMax;
|
||||
m_dt0 = netIn.dt0;
|
||||
|
||||
const double stiff_abs_tol = m_config.get<double>("Network:Approx8:Stiff:AbsTol", 1.0e-6);
|
||||
const double stiff_rel_tol = m_config.get<double>("Network:Approx8:Stiff:RelTol", 1.0e-6);
|
||||
const double nonstiff_abs_tol = m_config.get<double>("Network:Approx8:NonStiff:AbsTol", 1.0e-6);
|
||||
const double nonstiff_rel_tol = m_config.get<double>("Network:Approx8:NonStiff:RelTol", 1.0e-6);
|
||||
|
||||
int num_steps = -1;
|
||||
|
||||
if (m_stiff) {
|
||||
LOG_DEBUG(m_logger, "Using stiff solver for Approx8Network");
|
||||
num_steps = integrate_const(
|
||||
make_dense_output<rosenbrock4<double>>(stiff_abs_tol, stiff_rel_tol),
|
||||
std::make_pair(ODE(), Jacobian()),
|
||||
m_y,
|
||||
0.0,
|
||||
m_tMax,
|
||||
m_dt0
|
||||
);
|
||||
|
||||
} else {
|
||||
LOG_DEBUG(m_logger, "Using non stiff solver for Approx8Network");
|
||||
num_steps = integrate_const (
|
||||
make_dense_output<runge_kutta_dopri5<vector_type>>(nonstiff_abs_tol, nonstiff_rel_tol),
|
||||
ODE(),
|
||||
m_y,
|
||||
0.0,
|
||||
m_tMax,
|
||||
m_dt0
|
||||
);
|
||||
}
|
||||
|
||||
double ySum = 0.0;
|
||||
for (int i = 0; i < Approx8Net::nIso; i++) {
|
||||
m_y[i] *= Approx8Net::aIon[i];
|
||||
ySum += m_y[i];
|
||||
}
|
||||
for (int i = 0; i < Approx8Net::nIso; i++) {
|
||||
m_y[i] /= ySum;
|
||||
}
|
||||
|
||||
NetOut netOut;
|
||||
std::vector<double> outComposition;
|
||||
outComposition.reserve(Approx8Net::nVar);
|
||||
|
||||
for (int i = 0; i < Approx8Net::nIso; i++) {
|
||||
outComposition.push_back(m_y[i]);
|
||||
}
|
||||
netOut.energy = m_y[Approx8Net::iEnergy];
|
||||
netOut.num_steps = num_steps;
|
||||
|
||||
const std::vector<std::string> symbols = {"H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"};
|
||||
netOut.composition = fourdst::composition::Composition(symbols, outComposition);
|
||||
|
||||
return netOut;
|
||||
}
|
||||
|
||||
void Approx8Network::setStiff(bool stiff) {
|
||||
m_stiff = stiff;
|
||||
}
|
||||
|
||||
vector_type Approx8Network::convert_netIn(const NetIn &netIn) {
|
||||
vector_type y(Approx8Net::nVar, 0.0);
|
||||
y[Approx8Net::ih1] = netIn.composition.getMassFraction("H-1");
|
||||
y[Approx8Net::ihe3] = netIn.composition.getMassFraction("He-3");
|
||||
y[Approx8Net::ihe4] = netIn.composition.getMassFraction("He-4");
|
||||
y[Approx8Net::ic12] = netIn.composition.getMassFraction("C-12");
|
||||
y[Approx8Net::in14] = netIn.composition.getMassFraction("N-14");
|
||||
y[Approx8Net::io16] = netIn.composition.getMassFraction("O-16");
|
||||
y[Approx8Net::ine20] = netIn.composition.getMassFraction("Ne-20");
|
||||
y[Approx8Net::img24] = netIn.composition.getMassFraction("Mg-24");
|
||||
y[Approx8Net::iTemp] = netIn.temperature;
|
||||
y[Approx8Net::iDensity] = netIn.density;
|
||||
y[Approx8Net::iEnergy] = netIn.energy;
|
||||
|
||||
double ySum = 0.0;
|
||||
for (int i = 0; i < Approx8Net::nIso; i++) {
|
||||
y[i] /= Approx8Net::aIon[i];
|
||||
ySum += y[i];
|
||||
}
|
||||
for (int i = 0; i < Approx8Net::nIso; i++) {
|
||||
y[i] /= ySum;
|
||||
}
|
||||
|
||||
return y;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
// main program
|
||||
|
||||
@@ -1,69 +0,0 @@
|
||||
/* ***********************************************************************
|
||||
//
|
||||
// Copyright (C) 2025 -- The 4D-STAR Collaboration
|
||||
// File Authors: Aaron Dotter, Emily Boudreaux
|
||||
// Last Modified: March 21, 2025
|
||||
//
|
||||
// 4DSSE is free software; you can use it and/or modify
|
||||
// it under the terms and restrictions the GNU General Library Public
|
||||
// License version 3 (GPLv3) as published by the Free Software Foundation.
|
||||
//
|
||||
// 4DSSE is distributed in the hope that it will be useful,
|
||||
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
|
||||
// See the GNU Library General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU Library General Public License
|
||||
// along with this software; if not, write to the Free Software
|
||||
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
|
||||
//
|
||||
// *********************************************************************** */
|
||||
#include "network.h"
|
||||
#include "config.h"
|
||||
|
||||
#include "approx8.h"
|
||||
#include "logging.h"
|
||||
#include "quill/LogMacros.h"
|
||||
|
||||
namespace serif::network {
|
||||
Network::Network(const NetworkFormat format) :
|
||||
m_config(fourdst::config::Config::getInstance()),
|
||||
m_logManager(fourdst::logging::LogManager::getInstance()),
|
||||
m_logger(m_logManager.getLogger("log")),
|
||||
m_format(format) {
|
||||
if (format == NetworkFormat::UNKNOWN) {
|
||||
LOG_ERROR(m_logger, "nuclearNetwork::Network::Network() called with UNKNOWN format");
|
||||
throw std::runtime_error("nuclearNetwork::Network::Network() called with UNKNOWN format");
|
||||
}
|
||||
}
|
||||
|
||||
NetworkFormat Network::getFormat() const {
|
||||
return m_format;
|
||||
}
|
||||
|
||||
NetworkFormat Network::setFormat(const NetworkFormat format) {
|
||||
const NetworkFormat oldFormat = m_format;
|
||||
m_format = format;
|
||||
return oldFormat;
|
||||
}
|
||||
|
||||
NetOut Network::evaluate(const NetIn &netIn) {
|
||||
NetOut netOut;
|
||||
switch (m_format) {
|
||||
case APPROX8: {
|
||||
approx8::Approx8Network network;
|
||||
netOut = network.evaluate(netIn);
|
||||
break;
|
||||
}
|
||||
case UNKNOWN: {
|
||||
LOG_ERROR(m_logger, "Network format {} is not implemented.", FormatStringLookup.at(m_format));
|
||||
throw std::runtime_error("Network format not implemented.");
|
||||
}
|
||||
default: {
|
||||
LOG_ERROR(m_logger, "Unknown network format.");
|
||||
throw std::runtime_error("Unknown network format.");
|
||||
}
|
||||
}
|
||||
return netOut;
|
||||
}
|
||||
}
|
||||
@@ -1,333 +0,0 @@
|
||||
/* ***********************************************************************
|
||||
//
|
||||
// Copyright (C) 2025 -- The 4D-STAR Collaboration
|
||||
// File Author: Emily Boudreaux
|
||||
// Last Modified: March 21, 2025
|
||||
//
|
||||
// 4DSSE is free software; you can use it and/or modify
|
||||
// it under the terms and restrictions the GNU General Library Public
|
||||
// License version 3 (GPLv3) as published by the Free Software Foundation.
|
||||
//
|
||||
// 4DSSE is distributed in the hope that it will be useful,
|
||||
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
|
||||
// See the GNU Library General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU Library General Public License
|
||||
// along with this software; if not, write to the Free Software
|
||||
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
|
||||
//
|
||||
// *********************************************************************** */
|
||||
#pragma once
|
||||
|
||||
#include <array>
|
||||
|
||||
#include <boost/numeric/odeint.hpp>
|
||||
|
||||
#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' "approx8" and includes 8 isotopes and various nuclear reactions.
|
||||
* The rates are evaluated using a fitting function with coefficients from reaclib.jinaweb.org.
|
||||
*/
|
||||
|
||||
|
||||
namespace serif::network::approx8{
|
||||
|
||||
/**
|
||||
* @typedef vector_type
|
||||
* @brief Alias for a vector of doubles using Boost uBLAS.
|
||||
*/
|
||||
typedef boost::numeric::ublas::vector< double > vector_type;
|
||||
|
||||
/**
|
||||
* @typedef matrix_type
|
||||
* @brief Alias for a matrix of doubles using Boost uBLAS.
|
||||
*/
|
||||
typedef boost::numeric::ublas::matrix< double > matrix_type;
|
||||
|
||||
/**
|
||||
* @typedef vec7
|
||||
* @brief Alias for a std::array of 7 doubles.
|
||||
*/
|
||||
typedef std::array<double,7> vec7;
|
||||
|
||||
using namespace boost::numeric::odeint;
|
||||
|
||||
/**
|
||||
* @struct Approx8Net
|
||||
* @brief Contains constants and arrays related to the nuclear network.
|
||||
*/
|
||||
struct Approx8Net{
|
||||
static constexpr int ih1=0;
|
||||
static constexpr int ihe3=1;
|
||||
static constexpr int ihe4=2;
|
||||
static constexpr int ic12=3;
|
||||
static constexpr int in14=4;
|
||||
static constexpr int io16=5;
|
||||
static constexpr int ine20=6;
|
||||
static constexpr int img24=7;
|
||||
|
||||
static constexpr int iTemp=img24+1;
|
||||
static constexpr int iDensity =iTemp+1;
|
||||
static constexpr int iEnergy=iDensity+1;
|
||||
|
||||
static constexpr int nIso=img24+1; // number of isotopes
|
||||
static constexpr int nVar=iEnergy+1; // number of variables
|
||||
|
||||
static constexpr std::array<int,nIso> aIon = {
|
||||
1,
|
||||
3,
|
||||
4,
|
||||
12,
|
||||
14,
|
||||
16,
|
||||
20,
|
||||
24
|
||||
};
|
||||
|
||||
static constexpr std::array<double,nIso> 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 dfdt Derivative of the state vector.
|
||||
*/
|
||||
void operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) const;
|
||||
};
|
||||
|
||||
/**
|
||||
* @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.
|
||||
*/
|
||||
void operator() ( const vector_type &y, vector_type &dydt, double /* t */) const;
|
||||
};
|
||||
|
||||
/**
|
||||
* @class Approx8Network
|
||||
* @brief Class for the Approx8 nuclear reaction network.
|
||||
*/
|
||||
class Approx8Network final : public Network {
|
||||
public:
|
||||
Approx8Network();
|
||||
|
||||
/**
|
||||
* @brief Evaluates the nuclear network.
|
||||
* @param netIn Input parameters for the network.
|
||||
* @return Output results from the network.
|
||||
*/
|
||||
NetOut evaluate(const NetIn &netIn) override;
|
||||
|
||||
/**
|
||||
* @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() const { return m_stiff; }
|
||||
private:
|
||||
vector_type m_y;
|
||||
double m_tMax = 0;
|
||||
double m_dt0 = 0;
|
||||
bool m_stiff = false;
|
||||
|
||||
/**
|
||||
* @brief Converts the input parameters to the internal state vector.
|
||||
* @param netIn Input parameters for the network.
|
||||
* @return Internal state vector.
|
||||
*/
|
||||
static vector_type convert_netIn(const NetIn &netIn);
|
||||
};
|
||||
|
||||
|
||||
} // namespace nnApprox8
|
||||
@@ -1,129 +0,0 @@
|
||||
/* ***********************************************************************
|
||||
//
|
||||
// Copyright (C) 2025 -- The 4D-STAR Collaboration
|
||||
// File Authors: Emily Boudreaux, Aaron Dotter
|
||||
// Last Modified: March 21, 2025
|
||||
//
|
||||
// 4DSSE is free software; you can use it and/or modify
|
||||
// it under the terms and restrictions the GNU General Library Public
|
||||
// License version 3 (GPLv3) as published by the Free Software Foundation.
|
||||
//
|
||||
// 4DSSE is distributed in the hope that it will be useful,
|
||||
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
|
||||
// See the GNU Library General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU Library General Public License
|
||||
// along with this software; if not, write to the Free Software
|
||||
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
|
||||
//
|
||||
// *********************************************************************** */
|
||||
#pragma once
|
||||
|
||||
#include <vector>
|
||||
|
||||
#include "logging.h"
|
||||
#include "config.h"
|
||||
#include "quill/Logger.h"
|
||||
#include "composition.h"
|
||||
#include <unordered_map>
|
||||
|
||||
namespace serif::network {
|
||||
|
||||
enum NetworkFormat {
|
||||
APPROX8, ///< Approx8 nuclear reaction network format.
|
||||
UNKNOWN,
|
||||
};
|
||||
|
||||
static inline std::unordered_map<NetworkFormat, std::string> FormatStringLookup = {
|
||||
{APPROX8, "Approx8"},
|
||||
{UNKNOWN, "Unknown"}
|
||||
};
|
||||
|
||||
/**
|
||||
* @struct NetIn
|
||||
* @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 {
|
||||
fourdst::composition::Composition composition; ///< Composition of the network
|
||||
double tMax; ///< Maximum time
|
||||
double dt0; ///< Initial time step
|
||||
double temperature; ///< Temperature in Kelvin
|
||||
double density; ///< Density in g/cm^3
|
||||
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<double> composition = netOut.composition;
|
||||
* int steps = netOut.num_steps;
|
||||
* double energy = netOut.energy;
|
||||
* @endcode
|
||||
*/
|
||||
struct NetOut {
|
||||
fourdst::composition::Composition composition; ///< Composition of the network after evaluation
|
||||
int num_steps; ///< Number of steps taken in the evaluation
|
||||
double energy; ///< Energy in ergs after evaluation
|
||||
};
|
||||
|
||||
/**
|
||||
* @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:
|
||||
explicit Network(const NetworkFormat format = NetworkFormat::APPROX8);
|
||||
virtual ~Network() = default;
|
||||
|
||||
NetworkFormat getFormat() const;
|
||||
NetworkFormat setFormat(const NetworkFormat format);
|
||||
|
||||
/**
|
||||
* @brief Evaluate the network based on the input parameters.
|
||||
*
|
||||
* @param netIn Input parameters for the network evaluation.
|
||||
* @return NetOut Output results from the network evaluation.
|
||||
*/
|
||||
virtual NetOut evaluate(const NetIn &netIn);
|
||||
|
||||
protected:
|
||||
fourdst::config::Config& m_config; ///< Configuration instance
|
||||
fourdst::logging::LogManager& m_logManager; ///< Log manager instance
|
||||
quill::Logger* m_logger; ///< Logger instance
|
||||
|
||||
NetworkFormat m_format; ///< Format of the network
|
||||
};
|
||||
|
||||
|
||||
|
||||
} // namespace nuclearNetwork
|
||||
@@ -13,22 +13,22 @@ namespace py = pybind11;
|
||||
|
||||
|
||||
void register_network_bindings(pybind11::module &network_submodule) {
|
||||
py::enum_<serif::network::NetworkFormat>(network_submodule, "NetworkFormat")
|
||||
.value("APPROX8", serif::network::NetworkFormat::APPROX8)
|
||||
.value("UNKNOWN", serif::network::NetworkFormat::UNKNOWN)
|
||||
.def("__str__", [](const serif::network::NetworkFormat format) {
|
||||
return serif::network::FormatStringLookup[format];
|
||||
py::enum_<gridfire::NetworkFormat>(network_submodule, "NetworkFormat")
|
||||
.value("APPROX8", gridfire::NetworkFormat::APPROX8)
|
||||
.value("UNKNOWN", gridfire::NetworkFormat::UNKNOWN)
|
||||
.def("__str__", [](const gridfire::NetworkFormat format) {
|
||||
return gridfire::FormatStringLookup[format];
|
||||
});
|
||||
|
||||
py::class_<serif::network::NetIn>(network_submodule, "NetIn")
|
||||
py::class_<gridfire::NetIn>(network_submodule, "NetIn")
|
||||
.def(py::init<>())
|
||||
.def_readwrite("composition", &serif::network::NetIn::composition)
|
||||
.def_readwrite("tMax", &serif::network::NetIn::tMax)
|
||||
.def_readwrite("dt0", &serif::network::NetIn::dt0)
|
||||
.def_readwrite("temperature", &serif::network::NetIn::temperature)
|
||||
.def_readwrite("density", &serif::network::NetIn::density)
|
||||
.def_readwrite("energy", &serif::network::NetIn::energy)
|
||||
.def("__repr__", [](const serif::network::NetIn &netIn) {
|
||||
.def_readwrite("composition", &gridfire::NetIn::composition)
|
||||
.def_readwrite("tMax", &gridfire::NetIn::tMax)
|
||||
.def_readwrite("dt0", &gridfire::NetIn::dt0)
|
||||
.def_readwrite("temperature", &gridfire::NetIn::temperature)
|
||||
.def_readwrite("density", &gridfire::NetIn::density)
|
||||
.def_readwrite("energy", &gridfire::NetIn::energy)
|
||||
.def("__repr__", [](const gridfire::NetIn &netIn) {
|
||||
std::stringstream ss;
|
||||
ss << "NetIn(composition=" << netIn.composition
|
||||
<< ", tMax=" << netIn.tMax
|
||||
@@ -39,11 +39,11 @@ void register_network_bindings(pybind11::module &network_submodule) {
|
||||
return ss.str();
|
||||
});
|
||||
|
||||
py::class_<serif::network::NetOut>(network_submodule, "NetOut")
|
||||
.def_readonly("composition", &serif::network::NetOut::composition)
|
||||
.def_readonly("num_steps", &serif::network::NetOut::num_steps)
|
||||
.def_readonly("energy", &serif::network::NetOut::energy)
|
||||
.def("__repr__", [](const serif::network::NetOut &netOut) {
|
||||
py::class_<gridfire::NetOut>(network_submodule, "NetOut")
|
||||
.def_readonly("composition", &gridfire::NetOut::composition)
|
||||
.def_readonly("num_steps", &gridfire::NetOut::num_steps)
|
||||
.def_readonly("energy", &gridfire::NetOut::energy)
|
||||
.def("__repr__", [](const gridfire::NetOut &netOut) {
|
||||
std::stringstream ss;
|
||||
ss << "NetOut(composition=" << netOut.composition
|
||||
<< ", num_steps=" << netOut.num_steps
|
||||
@@ -51,14 +51,14 @@ void register_network_bindings(pybind11::module &network_submodule) {
|
||||
return ss.str();
|
||||
});
|
||||
|
||||
py::class_<serif::network::Network>(network_submodule, "Network")
|
||||
.def(py::init<serif::network::NetworkFormat>(), py::arg("format"))
|
||||
.def("evaluate", &serif::network::Network::evaluate, py::arg("netIn"))
|
||||
.def("getFormat", &serif::network::Network::getFormat)
|
||||
.def("setFormat", &serif::network::Network::setFormat, py::arg("format"))
|
||||
.def("__repr__", [](const serif::network::Network &network) {
|
||||
py::class_<gridfire::Network>(network_submodule, "Network")
|
||||
.def(py::init<gridfire::NetworkFormat>(), py::arg("format"))
|
||||
.def("evaluate", &gridfire::Network::evaluate, py::arg("netIn"))
|
||||
.def("getFormat", &gridfire::Network::getFormat)
|
||||
.def("setFormat", &gridfire::Network::setFormat, py::arg("format"))
|
||||
.def("__repr__", [](const gridfire::Network &network) {
|
||||
std::stringstream ss;
|
||||
ss << "Network(format=" << serif::network::FormatStringLookup[network.getFormat()] << ")";
|
||||
ss << "Network(format=" << gridfire::FormatStringLookup[network.getFormat()] << ")";
|
||||
return ss.str();
|
||||
});
|
||||
|
||||
@@ -68,7 +68,7 @@ void register_network_bindings(pybind11::module &network_submodule) {
|
||||
}
|
||||
|
||||
void register_approx8_bindings(pybind11::module &network_submodule) {
|
||||
using namespace serif::network::approx8;
|
||||
using namespace gridfire::approx8;
|
||||
|
||||
py::class_<vec7>(network_submodule, "vec7")
|
||||
.def("__getitem__", [](const vec7 &v, const size_t i) {
|
||||
|
||||
@@ -5,7 +5,7 @@ bindings_headers = files('bindings.h')
|
||||
dependencies = [
|
||||
python3_dep,
|
||||
pybind11_dep,
|
||||
network_dep,
|
||||
gridfire_dep,
|
||||
]
|
||||
|
||||
shared_module('py_network',
|
||||
|
||||
4
subprojects/GridFire.wrap
Normal file
4
subprojects/GridFire.wrap
Normal file
@@ -0,0 +1,4 @@
|
||||
[wrap-git]
|
||||
url = https://github.com/4D-STAR/GridFire.git
|
||||
revision = v0.0.3a
|
||||
depth = 1
|
||||
@@ -8,7 +8,6 @@ subdir('meshIO')
|
||||
subdir('probe')
|
||||
subdir('eos')
|
||||
subdir('resource')
|
||||
subdir('network')
|
||||
subdir('poly')
|
||||
|
||||
|
||||
|
||||
@@ -1,63 +0,0 @@
|
||||
#include <gtest/gtest.h>
|
||||
#include <string>
|
||||
|
||||
#include "approx8.h"
|
||||
#include "config.h"
|
||||
#include "network.h"
|
||||
#include "composition.h"
|
||||
|
||||
#include <vector>
|
||||
|
||||
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) {
|
||||
fourdst::config::Config& config = fourdst::config::Config::getInstance();
|
||||
config.loadConfig(TEST_CONFIG);
|
||||
EXPECT_NO_THROW(serif::network::approx8::Approx8Network());
|
||||
}
|
||||
|
||||
TEST_F(approx8Test, setStiff) {
|
||||
serif::network::approx8::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) {
|
||||
serif::network::approx8::Approx8Network network;
|
||||
serif::network::NetIn netIn;
|
||||
|
||||
std::vector<double> comp = {0.708, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4};
|
||||
std::vector<std::string> symbols = {"H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"};
|
||||
|
||||
fourdst::composition::Composition composition;
|
||||
composition.registerSymbol(symbols, true);
|
||||
composition.setMassFraction(symbols, comp);
|
||||
bool isFinalized = composition.finalize(true);
|
||||
EXPECT_TRUE(isFinalized);
|
||||
|
||||
netIn.composition = composition;
|
||||
netIn.temperature = 1e7;
|
||||
netIn.density = 1e2;
|
||||
netIn.energy = 0.0;
|
||||
|
||||
netIn.tMax = 3.15e17;
|
||||
netIn.dt0 = 1e12;
|
||||
|
||||
serif::network::NetOut netOut;
|
||||
EXPECT_NO_THROW(netOut = network.evaluate(netIn));
|
||||
|
||||
double energyFraction = netOut.energy / 1.6433051127589775E+18;
|
||||
double H1MassFraction = netOut.composition.getMassFraction("H-1")/ 0.50166262445895604;
|
||||
double He4MassFraction = netOut.composition.getMassFraction("He-4") / 0.48172273720971226;
|
||||
|
||||
double relError = 1e-6;
|
||||
EXPECT_NEAR(H1MassFraction, 1.0, relError);
|
||||
EXPECT_NEAR(He4MassFraction, 1.0, relError);
|
||||
EXPECT_NEAR(energyFraction, 1.0, relError);
|
||||
}
|
||||
@@ -1,25 +0,0 @@
|
||||
# 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, gtest_main, network_dep, composition_dep],
|
||||
include_directories: include_directories('../../src/network/public'),
|
||||
link_with: libnetwork, # Link the dobj library
|
||||
install_rpath: '@loader_path/../../src' # Ensure runtime library path resolves correctly
|
||||
)
|
||||
|
||||
# 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
|
||||
Reference in New Issue
Block a user