From 4cb1a4841aa1b2c038bfa45082a636530708feb7 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Sat, 21 Jun 2025 13:50:04 -0400 Subject: [PATCH] feat(GridFire): SERiF now uses GridFire --- .gitignore | 4 + build-config/GridFire/meson.build | 2 + build-config/meson.build | 2 +- build-config/quill/meson.build | 10 - build-config/yaml-cpp/meson.build | 12 - build-python/meson.build | 2 +- meson.build | 4 +- src/meson.build | 1 - src/network/meson.build | 37 -- src/network/private/approx8.cpp | 542 ------------------------------ src/network/private/network.cpp | 69 ---- src/network/public/approx8.h | 333 ------------------ src/network/public/network.h | 129 ------- src/python/network/bindings.cpp | 52 +-- src/python/network/meson.build | 2 +- subprojects/GridFire.wrap | 4 + tests/meson.build | 1 - tests/network/approx8Test.cpp | 63 ---- tests/network/meson.build | 25 -- 19 files changed, 42 insertions(+), 1252 deletions(-) create mode 100644 build-config/GridFire/meson.build delete mode 100644 build-config/quill/meson.build delete mode 100644 build-config/yaml-cpp/meson.build delete mode 100644 src/network/meson.build delete mode 100644 src/network/private/approx8.cpp delete mode 100644 src/network/private/network.cpp delete mode 100644 src/network/public/approx8.h delete mode 100644 src/network/public/network.h create mode 100644 subprojects/GridFire.wrap delete mode 100644 tests/network/approx8Test.cpp delete mode 100644 tests/network/meson.build diff --git a/.gitignore b/.gitignore index 8e615c6..f00408a 100644 --- a/.gitignore +++ b/.gitignore @@ -72,8 +72,12 @@ subprojects/libconstants/ subprojects/liblogging/ subprojects/libconfig/ subprojects/libcomposition/ +subprojects/GridFire/ qhull.wrap +quill.wrap +yaml-cpp.wrap +cppad.wrap .vscode/ diff --git a/build-config/GridFire/meson.build b/build-config/GridFire/meson.build new file mode 100644 index 0000000..63e4f80 --- /dev/null +++ b/build-config/GridFire/meson.build @@ -0,0 +1,2 @@ +gridfire_p = subproject('GridFire') +gridfire_dep = gridfire_p.get_variable('network_dep') diff --git a/build-config/meson.build b/build-config/meson.build index 9a48c5a..9a7fc1b 100644 --- a/build-config/meson.build +++ b/build-config/meson.build @@ -1,9 +1,9 @@ cmake = import('cmake') subdir('fourdst') +subdir('GridFire') subdir('mfem') -subdir('yaml-cpp') subdir('boost') subdir('opatIO') subdir('mpi') diff --git a/build-config/quill/meson.build b/build-config/quill/meson.build deleted file mode 100644 index a069255..0000000 --- a/build-config/quill/meson.build +++ /dev/null @@ -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') \ No newline at end of file diff --git a/build-config/yaml-cpp/meson.build b/build-config/yaml-cpp/meson.build deleted file mode 100644 index 72f76c0..0000000 --- a/build-config/yaml-cpp/meson.build +++ /dev/null @@ -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') \ No newline at end of file diff --git a/build-python/meson.build b/build-python/meson.build index 6d03a5a..29db376 100644 --- a/build-python/meson.build +++ b/build-python/meson.build @@ -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, diff --git a/meson.build b/meson.build index e0b9301..e25c551 100644 --- a/meson.build +++ b/meson.build @@ -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') diff --git a/src/meson.build b/src/meson.build index 8436ab7..e9cedf9 100644 --- a/src/meson.build +++ b/src/meson.build @@ -15,7 +15,6 @@ subdir('meshIO') subdir('resource') # Physics Libraries -subdir('network') subdir('polytrope') # Python Bindings diff --git a/src/network/meson.build b/src/network/meson.build deleted file mode 100644 index 042a2c3..0000000 --- a/src/network/meson.build +++ /dev/null @@ -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') \ No newline at end of file diff --git a/src/network/private/approx8.cpp b/src/network/private/approx8.cpp deleted file mode 100644 index b2da919..0000000 --- a/src/network/private/approx8.cpp +++ /dev/null @@ -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 -#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' "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("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 < 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 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 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 - diff --git a/src/network/private/network.cpp b/src/network/private/network.cpp deleted file mode 100644 index f377c9b..0000000 --- a/src/network/private/network.cpp +++ /dev/null @@ -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; - } -} diff --git a/src/network/public/approx8.h b/src/network/public/approx8.h deleted file mode 100644 index 634225e..0000000 --- a/src/network/public/approx8.h +++ /dev/null @@ -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 - -#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' "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 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 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 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 diff --git a/src/network/public/network.h b/src/network/public/network.h deleted file mode 100644 index 349747a..0000000 --- a/src/network/public/network.h +++ /dev/null @@ -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 - -#include "logging.h" -#include "config.h" -#include "quill/Logger.h" -#include "composition.h" -#include - -namespace serif::network { - - enum NetworkFormat { - APPROX8, ///< Approx8 nuclear reaction network format. - UNKNOWN, - }; - - static inline std::unordered_map FormatStringLookup = { - {APPROX8, "Approx8"}, - {UNKNOWN, "Unknown"} - }; - - /** - * @struct NetIn - * @brief Input structure for the network evaluation. - * - * This structure holds the input parameters required for the network evaluation. - * - * Example usage: - * @code - * nuclearNetwork::NetIn netIn; - * 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 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 diff --git a/src/python/network/bindings.cpp b/src/python/network/bindings.cpp index 8d1c34d..7c97302 100644 --- a/src/python/network/bindings.cpp +++ b/src/python/network/bindings.cpp @@ -13,22 +13,22 @@ namespace py = pybind11; void register_network_bindings(pybind11::module &network_submodule) { - py::enum_(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_(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_(network_submodule, "NetIn") + py::class_(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_(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_(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_(network_submodule, "Network") - .def(py::init(), 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_(network_submodule, "Network") + .def(py::init(), 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_(network_submodule, "vec7") .def("__getitem__", [](const vec7 &v, const size_t i) { diff --git a/src/python/network/meson.build b/src/python/network/meson.build index 47ea8a1..cf594be 100644 --- a/src/python/network/meson.build +++ b/src/python/network/meson.build @@ -5,7 +5,7 @@ bindings_headers = files('bindings.h') dependencies = [ python3_dep, pybind11_dep, - network_dep, + gridfire_dep, ] shared_module('py_network', diff --git a/subprojects/GridFire.wrap b/subprojects/GridFire.wrap new file mode 100644 index 0000000..0315252 --- /dev/null +++ b/subprojects/GridFire.wrap @@ -0,0 +1,4 @@ +[wrap-git] +url = https://github.com/4D-STAR/GridFire.git +revision = v0.0.3a +depth = 1 \ No newline at end of file diff --git a/tests/meson.build b/tests/meson.build index d3ea172..ec25dcb 100644 --- a/tests/meson.build +++ b/tests/meson.build @@ -8,7 +8,6 @@ subdir('meshIO') subdir('probe') subdir('eos') subdir('resource') -subdir('network') subdir('poly') diff --git a/tests/network/approx8Test.cpp b/tests/network/approx8Test.cpp deleted file mode 100644 index d837275..0000000 --- a/tests/network/approx8Test.cpp +++ /dev/null @@ -1,63 +0,0 @@ -#include -#include - -#include "approx8.h" -#include "config.h" -#include "network.h" -#include "composition.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) { - 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 comp = {0.708, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4}; - std::vector symbols = {"H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"}; - - 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); -} diff --git a/tests/network/meson.build b/tests/network/meson.build deleted file mode 100644 index 72a4a55..0000000 --- a/tests/network/meson.build +++ /dev/null @@ -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