From cd191cff231fb9501e8abc254ac2bf0be2aa958a Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 26 Jun 2025 15:13:46 -0400 Subject: [PATCH] feat(GridFire): major design changes Switching to an Engine + solver design. Also brought xxHash and Eigen in. Working on QSE and Culling. --- .../static/reaclib/include/gridfire/reaclib.h | 254 -------- .../reaclib/include/gridfire/reactions.h | 5 +- assets/static/reaclib/meson.build | 2 - build-config/eigen/meson.build | 5 + build-config/meson.build | 2 + build-config/xxHash/LICENSE | 33 + build-config/xxHash/include/constexpr-xxh3.h | 355 +++++++++++ build-config/xxHash/include/xxhash64.h | 202 +++++++ build-config/xxHash/meson.build | 4 + meson.build | 25 + meson_options.txt | 21 +- .../include/gridfire/engine/engine_abstract.h | 66 ++ .../{approx8.h => engine/engine_approx8.h} | 0 .../include/gridfire/engine/engine_culled.h | 55 ++ .../include/gridfire/engine/engine_graph.h | 277 +++++++++ src/network/include/gridfire/netgraph.h | 571 ------------------ src/network/include/gridfire/network.h | 111 +--- .../include/gridfire/reaction/reaction.h | 270 +++++++++ src/network/include/gridfire/solver/solver.h | 256 ++++++++ .../engine_approx8.cpp} | 34 +- .../{netgraph.cpp => engine/engine_graph.cpp} | 402 ++++++------ src/network/lib/network.cpp | 322 +++------- src/network/lib/reaction/reaction.cpp | 427 +++++++++++++ src/network/lib/solver/solver.cpp | 384 ++++++++++++ src/network/meson.build | 19 +- subprojects/eigen.wrap | 13 + subprojects/libcomposition.wrap | 2 +- tests/graphnet_sandbox/approx8 | 0 tests/graphnet_sandbox/main.cpp | 48 ++ tests/graphnet_sandbox/meson.build | 5 + tests/meson.build | 1 + tests/network/approx8Test.cpp | 7 +- 32 files changed, 2737 insertions(+), 1441 deletions(-) delete mode 100644 assets/static/reaclib/include/gridfire/reaclib.h create mode 100644 build-config/eigen/meson.build create mode 100644 build-config/xxHash/LICENSE create mode 100644 build-config/xxHash/include/constexpr-xxh3.h create mode 100644 build-config/xxHash/include/xxhash64.h create mode 100644 build-config/xxHash/meson.build create mode 100644 src/network/include/gridfire/engine/engine_abstract.h rename src/network/include/gridfire/{approx8.h => engine/engine_approx8.h} (100%) create mode 100644 src/network/include/gridfire/engine/engine_culled.h create mode 100644 src/network/include/gridfire/engine/engine_graph.h delete mode 100644 src/network/include/gridfire/netgraph.h create mode 100644 src/network/include/gridfire/reaction/reaction.h create mode 100644 src/network/include/gridfire/solver/solver.h rename src/network/lib/{approx8.cpp => engine/engine_approx8.cpp} (95%) rename src/network/lib/{netgraph.cpp => engine/engine_graph.cpp} (53%) create mode 100644 src/network/lib/reaction/reaction.cpp create mode 100644 src/network/lib/solver/solver.cpp create mode 100644 subprojects/eigen.wrap create mode 100644 tests/graphnet_sandbox/approx8 create mode 100644 tests/graphnet_sandbox/main.cpp create mode 100644 tests/graphnet_sandbox/meson.build diff --git a/assets/static/reaclib/include/gridfire/reaclib.h b/assets/static/reaclib/include/gridfire/reaclib.h deleted file mode 100644 index be2ffcc7..00000000 --- a/assets/static/reaclib/include/gridfire/reaclib.h +++ /dev/null @@ -1,254 +0,0 @@ -#pragma once -#include -#include -#include -#include -#include -#include -#include "fourdst/composition/atomicSpecies.h" - -#include "cppad/cppad.hpp" - -namespace gridfire::reaclib { - /** - * @struct RateFitSet - * @brief Holds the seven fitting parameters for a single REACLIB rate set. - * @details The thermonuclear reaction rate for a single set is calculated as: - * rate = exp(a0 + a1/T9 + a2/T9^(-1/3) + a3*T9^(1/3) + a4*T9 + a5*T9^(5/3) + a6*ln(T9)) - * where T9 is the temperature in billions of Kelvin. The total rate for a - * reaction is the sum of the rates from all its sets. - */ - struct RateFitSet { - double a0, a1, a2, a3, a4, a5, a6; - }; - - /** - * @struct REACLIBReaction - * @brief Represents a single nuclear reaction from the JINA REACLIB database. - * @details This struct is designed to be constructed at compile time (constexpr) from - * the data parsed by the Python generation script. It stores all necessary - * information to identify a reaction and calculate its rate. - */ - class REACLIBReaction { - private: - std::string m_id; ///< Unique identifier for the reaction, generated by the Python script. - std::string_view m_peName; ///< Name of the reaction in (projectile, ejectile) notation (e.g. p(p, g)d) - int m_chapter; ///< Chapter number from the REACLIB database, defining the reaction structure. - std::vector m_reactantNames; ///< Names of the reactant species involved in the reaction. - std::vector m_productNames; ///< Names of the product species produced by the reaction. - double m_qValue; ///< Q-value of the reaction in MeV, representing the energy released or absorbed. - std::string m_sourceLabel; ///< Source label for the rate data, indicating the origin of the rate coefficients (e.g., "wc12w", "st08"). - RateFitSet m_rateSets; ///< Fitting coefficients for the reaction rate, containing the seven parameters used in the rate calculation. - bool m_reverse = false; ///< indicates if the reaction is reversed - - friend bool operator==(const REACLIBReaction& lhs, const REACLIBReaction& rhs); - friend bool operator!=(const REACLIBReaction& lhs, const REACLIBReaction& rhs); - public: - /** - * @brief Constructs a REACLIBReaction object at compile time. - * @param id A unique string identifier generated by the Python script. - * @param chapter The REACLIB chapter number, defining the reaction structure. - * @param reactants A vector of strings with the names of the reactant species. - * @param products A vector of strings with the names of the product species. - * @param qValue The Q-value of the reaction in MeV. - * @param label The source label for the rate data (e.g., "wc12w", "st08"). - * @param sets A vector of RateFitSet, containing the fitting coefficients for the rate. - * @param reverse A boolean indicating if the reaction is reversed (default is false). - */ - REACLIBReaction( - std::string id, - std::string_view peName, - int chapter, - std::vector reactants, - std::vector products, - double qValue, - std::string label, - RateFitSet sets, - bool reverse = false) - : m_id(std::move(id)), - m_peName(peName), - m_chapter(chapter), - m_reactantNames(std::move(reactants)), - m_productNames(std::move(products)), - m_qValue(qValue), - m_sourceLabel(std::move(label)), - m_rateSets(sets), - m_reverse(reverse) {} - - template - [[nodiscard]] GeneralScalarType calculate_rate(const GeneralScalarType T9) const { - const GeneralScalarType T913 = CppAD::pow(T9, 1.0/3.0); - const GeneralScalarType rateExponent = m_rateSets.a0 + - m_rateSets.a1 / T9 + - m_rateSets.a2 / T913 + - m_rateSets.a3 * T913 + - m_rateSets.a4 * T9 + - m_rateSets.a5 * CppAD::pow(T9, 5.0/3.0) + - m_rateSets.a6 * CppAD::log(T9); - return CppAD::exp(rateExponent); - } - - [[nodiscard]] const std::string& id() const { return m_id; } - - [[nodiscard]] std::string_view peName() const { return m_peName; } - - [[nodiscard]] int chapter() const { return m_chapter; } - - [[nodiscard]] const std::vector& reactants() const { return m_reactantNames; } - - [[nodiscard]] const std::vector& products() const { return m_productNames; } - - [[nodiscard]] double qValue() const { return m_qValue; } - - [[nodiscard]] std::string sourceLabel() const { return m_sourceLabel; } - - [[nodiscard]] bool is_reverse() const { return m_reverse; } - - [[nodiscard]] RateFitSet rateFits() const { return m_rateSets; } - - friend std::ostream& operator<<(std::ostream& os, const REACLIBReaction& reaction) { - os << "REACLIBReaction(" << reaction.m_peName << ", " - << "Chapter: " << reaction.m_chapter << ")"; - return os; - } - - friend bool operator==(const REACLIBReaction& lhs, const REACLIBReaction& rhs); - - friend bool operator!=(const REACLIBReaction& lhs, const REACLIBReaction& rhs); - }; - - class REACLIBReactionSet { - private: - std::vector m_reactions; - public: - REACLIBReactionSet() = default; - - explicit REACLIBReactionSet(std::vector r) : m_reactions(std::move(r)) {} - - void add_reaction(const REACLIBReaction& reaction) { - if (not contains(reaction.id())) { - m_reactions.push_back(reaction); - } - } - - void remove_reaction(const REACLIBReaction& reaction) { - if (contains(reaction.id())) { - m_reactions.erase(std::remove_if(m_reactions.begin(), m_reactions.end(), - [&reaction](const REACLIBReaction& r) { return r.id() == reaction.id(); }), m_reactions.end()); - } - } - - [[nodiscard]] bool contains(const std::string& id) const { - for (const auto& reaction : m_reactions) { - if (reaction.id() == id) { - return true; - } - } - return false; - } - - [[nodiscard]] bool contains(const REACLIBReaction& reaction) const { - for (const auto& r : m_reactions) { - if (r == reaction) { - return true; - } - } - return false; - } - - [[nodiscard]] size_t size() const { - return m_reactions.size(); - } - - friend std::ostream& operator<<(std::ostream& os, const REACLIBReactionSet& set) { - os << "REACLIBReactionSet("; - int reactionNo = 0; - for (const auto& reaction : set.m_reactions) { - reactionNo++; - os << reaction.id(); - if (reactionNo < set.m_reactions.size()) { - os << ", "; - } - } - os << ")"; - return os; - } - - void sort(double T9=1.0) { - // Sort based on the evaluation of each reaction's rate function - std::sort(m_reactions.begin(), m_reactions.end(), - [T9](const REACLIBReaction& a, const REACLIBReaction& b) { - return a.calculate_rate(T9) > b.calculate_rate(T9); - }); - } - - auto begin() { - return m_reactions.begin(); - } - - auto end() { - return m_reactions.end(); - } - - bool containsSpecies(const fourdst::atomic::Species &species) const { - for (const auto& reaction : m_reactions) { - if (std::ranges::find(reaction.reactants(), species) != reaction.reactants().end() || - std::ranges::find(reaction.products(), species) != reaction.products().end()) { - return true; - } - } - return false; - } - - [[nodiscard]] const REACLIBReaction& operator[](size_t index) const { - if (index >= m_reactions.size()) { - throw std::out_of_range("Index out of range in REACLIBReactionSet."); - } - return m_reactions[index]; - } - - [[nodiscard]] std::vector get_reactions() const { - return m_reactions; - } - - [[nodiscard]] auto begin() const { - return m_reactions.cbegin(); - } - - [[nodiscard]] auto end() const { - return m_reactions.cend(); - } - - }; - static std::unordered_map s_all_reaclib_reactions; - static bool s_initialized = false; - - inline bool operator==(const REACLIBReaction& lhs, const REACLIBReaction& rhs) { - return lhs.m_id == rhs.m_id; - } - inline bool operator!=(const REACLIBReaction& lhs, const REACLIBReaction& rhs) { - return !(lhs == rhs); - } - - inline bool operator==(const REACLIBReactionSet& lhs, const REACLIBReactionSet& rhs) { - if (lhs.size() != rhs.size()) return false; - for (const auto& reaction : lhs) { - if (!rhs.contains(reaction)) return false; - } - return true; - } - - inline bool operator!=(const REACLIBReactionSet& lhs, const REACLIBReactionSet& rhs) { - return !(lhs == rhs); - } - -} - -namespace std { - template<> - struct hash { - size_t operator()(const gridfire::reaclib::REACLIBReaction& r) const noexcept { - return std::hash()(r.id()); - } - }; -} // namespace std diff --git a/assets/static/reaclib/include/gridfire/reactions.h b/assets/static/reaclib/include/gridfire/reactions.h index 7876fbf1..8df8f7cf 100644 --- a/assets/static/reaclib/include/gridfire/reactions.h +++ b/assets/static/reaclib/include/gridfire/reactions.h @@ -9,12 +9,15 @@ #pragma once #include "fourdst/composition/atomicSpecies.h" #include "fourdst/composition/species.h" -#include "reaclib.h" +#include "gridfire/reaction/reaction.h" namespace gridfire::reaclib { + static std::unordered_map s_all_reaclib_reactions; + static bool s_initialized = false; inline void initializeAllReaclibReactions() { + using namespace reaction; if (s_initialized) return; // already initialized s_initialized = true; s_all_reaclib_reactions.clear(); diff --git a/assets/static/reaclib/meson.build b/assets/static/reaclib/meson.build index 80d086ad..11bcf24c 100644 --- a/assets/static/reaclib/meson.build +++ b/assets/static/reaclib/meson.build @@ -1,5 +1,4 @@ required_headers = [ - 'gridfire/reaclib.h', 'gridfire/reactions.h', ] @@ -14,7 +13,6 @@ reaclib_reactions_dep = declare_dependency( message('✅ GridFire reaclib_reactions dependency declared') to_install_headers = [ - 'include/gridfire/reaclib.h', 'include/gridfire/reactions.h', ] install_headers(to_install_headers, subdir: 'gridfire/gridfire') diff --git a/build-config/eigen/meson.build b/build-config/eigen/meson.build new file mode 100644 index 00000000..ce160f25 --- /dev/null +++ b/build-config/eigen/meson.build @@ -0,0 +1,5 @@ +eigen_dep = dependency( + 'eigen3', + version : '>=3.3', + fallback : ['eigen', 'eigen_dep'] +) \ No newline at end of file diff --git a/build-config/meson.build b/build-config/meson.build index 7721f292..6087921d 100644 --- a/build-config/meson.build +++ b/build-config/meson.build @@ -3,5 +3,7 @@ cmake = import('cmake') subdir('fourdst') subdir('boost') subdir('cppad') +subdir('xxHash') +subdir('eigen') diff --git a/build-config/xxHash/LICENSE b/build-config/xxHash/LICENSE new file mode 100644 index 00000000..94092e1f --- /dev/null +++ b/build-config/xxHash/LICENSE @@ -0,0 +1,33 @@ +xxHash - Extremely Fast Hash algorithm +Header File +Copyright (C) 2012-2021 Yann Collet + +BSD 2-Clause License (https://www.opensource.org/licenses/bsd-license.php) + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following disclaimer + in the documentation and/or other materials provided with the + distribution. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +You can contact the author at: + - xxHash homepage: https://www.xxhash.com + - xxHash source repository: https://github.com/Cyan4973/xxHash + diff --git a/build-config/xxHash/include/constexpr-xxh3.h b/build-config/xxHash/include/constexpr-xxh3.h new file mode 100644 index 00000000..c5840bc8 --- /dev/null +++ b/build-config/xxHash/include/constexpr-xxh3.h @@ -0,0 +1,355 @@ +/* +BSD 2-Clause License + +constexpr-xxh3 - C++20 constexpr implementation of the XXH3 64-bit variant of xxHash +Copyright (c) 2021-2023, chys +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + +/* +This file uses code from Yann Collet's xxHash implementation. +Original xxHash copyright notice: + +xxHash - Extremely Fast Hash algorithm +Header File +Copyright (C) 2012-2020 Yann Collet +*/ + +#pragma once + +#include +#include +#include // for std::data, std::size +#include +#include + +namespace constexpr_xxh3 { + +template +concept ByteType = (std::is_integral_v && sizeof(T) == 1) +#if defined __cpp_lib_byte && __cpp_lib_byte >= 201603 + || std::is_same_v +#endif + ; + +template +concept BytePtrType = requires (T ptr) { + requires std::is_pointer_v; + requires ByteType>; +}; + +template +concept BytesType = requires (const T& bytes) { + { std::data(bytes) }; + requires BytePtrType; + // -> std::convertible_to is not supported widely enough + { static_cast(std::size(bytes)) }; +}; + +inline constexpr uint32_t swap32(uint32_t x) noexcept { + return ((x << 24) & 0xff000000) | ((x << 8) & 0x00ff0000) | + ((x >> 8) & 0x0000ff00) | ((x >> 24) & 0x000000ff); +} + +template +inline constexpr uint32_t readLE32(const T* ptr) noexcept { + return uint8_t(ptr[0]) | uint32_t(uint8_t(ptr[1])) << 8 | + uint32_t(uint8_t(ptr[2])) << 16 | uint32_t(uint8_t(ptr[3])) << 24; +} + +inline constexpr uint64_t swap64(uint64_t x) noexcept { + return ((x << 56) & 0xff00000000000000ULL) | + ((x << 40) & 0x00ff000000000000ULL) | + ((x << 24) & 0x0000ff0000000000ULL) | + ((x << 8) & 0x000000ff00000000ULL) | + ((x >> 8) & 0x00000000ff000000ULL) | + ((x >> 24) & 0x0000000000ff0000ULL) | + ((x >> 40) & 0x000000000000ff00ULL) | + ((x >> 56) & 0x00000000000000ffULL); +} + +template +inline constexpr uint64_t readLE64(const T* ptr) noexcept { + return readLE32(ptr) | uint64_t(readLE32(ptr + 4)) << 32; +} + +inline constexpr void writeLE64(uint8_t* dst, uint64_t v) noexcept { + for (int i = 0; i < 8; ++i) dst[i] = uint8_t(v >> (i * 8)); +} + +inline constexpr uint32_t PRIME32_1 = 0x9E3779B1U; +inline constexpr uint32_t PRIME32_2 = 0x85EBCA77U; +inline constexpr uint32_t PRIME32_3 = 0xC2B2AE3DU; + +inline constexpr uint64_t PRIME64_1 = 0x9E3779B185EBCA87ULL; +inline constexpr uint64_t PRIME64_2 = 0xC2B2AE3D27D4EB4FULL; +inline constexpr uint64_t PRIME64_3 = 0x165667B19E3779F9ULL; +inline constexpr uint64_t PRIME64_4 = 0x85EBCA77C2B2AE63ULL; +inline constexpr uint64_t PRIME64_5 = 0x27D4EB2F165667C5ULL; + +inline constexpr size_t SECRET_DEFAULT_SIZE = 192; +inline constexpr size_t SECRET_SIZE_MIN = 136; + +inline constexpr uint8_t kSecret[SECRET_DEFAULT_SIZE]{ + 0xb8, 0xfe, 0x6c, 0x39, 0x23, 0xa4, 0x4b, 0xbe, 0x7c, 0x01, 0x81, 0x2c, + 0xf7, 0x21, 0xad, 0x1c, 0xde, 0xd4, 0x6d, 0xe9, 0x83, 0x90, 0x97, 0xdb, + 0x72, 0x40, 0xa4, 0xa4, 0xb7, 0xb3, 0x67, 0x1f, 0xcb, 0x79, 0xe6, 0x4e, + 0xcc, 0xc0, 0xe5, 0x78, 0x82, 0x5a, 0xd0, 0x7d, 0xcc, 0xff, 0x72, 0x21, + 0xb8, 0x08, 0x46, 0x74, 0xf7, 0x43, 0x24, 0x8e, 0xe0, 0x35, 0x90, 0xe6, + 0x81, 0x3a, 0x26, 0x4c, 0x3c, 0x28, 0x52, 0xbb, 0x91, 0xc3, 0x00, 0xcb, + 0x88, 0xd0, 0x65, 0x8b, 0x1b, 0x53, 0x2e, 0xa3, 0x71, 0x64, 0x48, 0x97, + 0xa2, 0x0d, 0xf9, 0x4e, 0x38, 0x19, 0xef, 0x46, 0xa9, 0xde, 0xac, 0xd8, + 0xa8, 0xfa, 0x76, 0x3f, 0xe3, 0x9c, 0x34, 0x3f, 0xf9, 0xdc, 0xbb, 0xc7, + 0xc7, 0x0b, 0x4f, 0x1d, 0x8a, 0x51, 0xe0, 0x4b, 0xcd, 0xb4, 0x59, 0x31, + 0xc8, 0x9f, 0x7e, 0xc9, 0xd9, 0x78, 0x73, 0x64, 0xea, 0xc5, 0xac, 0x83, + 0x34, 0xd3, 0xeb, 0xc3, 0xc5, 0x81, 0xa0, 0xff, 0xfa, 0x13, 0x63, 0xeb, + 0x17, 0x0d, 0xdd, 0x51, 0xb7, 0xf0, 0xda, 0x49, 0xd3, 0x16, 0x55, 0x26, + 0x29, 0xd4, 0x68, 0x9e, 0x2b, 0x16, 0xbe, 0x58, 0x7d, 0x47, 0xa1, 0xfc, + 0x8f, 0xf8, 0xb8, 0xd1, 0x7a, 0xd0, 0x31, 0xce, 0x45, 0xcb, 0x3a, 0x8f, + 0x95, 0x16, 0x04, 0x28, 0xaf, 0xd7, 0xfb, 0xca, 0xbb, 0x4b, 0x40, 0x7e, +}; + +inline constexpr std::pair mult64to128( + uint64_t lhs, uint64_t rhs) noexcept { + uint64_t lo_lo = uint64_t(uint32_t(lhs)) * uint32_t(rhs); + uint64_t hi_lo = (lhs >> 32) * uint32_t(rhs); + uint64_t lo_hi = uint32_t(lhs) * (rhs >> 32); + uint64_t hi_hi = (lhs >> 32) * (rhs >> 32); + uint64_t cross = (lo_lo >> 32) + uint32_t(hi_lo) + lo_hi; + uint64_t upper = (hi_lo >> 32) + (cross >> 32) + hi_hi; + uint64_t lower = (cross << 32) | uint32_t(lo_lo); + return {lower, upper}; +} + +inline constexpr uint64_t mul128_fold64(uint64_t lhs, uint64_t rhs) noexcept { +#if defined __GNUC__ && __WORDSIZE >= 64 + // It appears both GCC and Clang support evaluating __int128 as constexpr + auto product = static_cast(lhs) * rhs; + return uint64_t(product >> 64) ^ uint64_t(product); +#else + auto product = mult64to128(lhs, rhs); + return product.first ^ product.second; +#endif +} + +inline constexpr uint64_t XXH64_avalanche(uint64_t h) noexcept { + h = (h ^ (h >> 33)) * PRIME64_2; + h = (h ^ (h >> 29)) * PRIME64_3; + return h ^ (h >> 32); +} + +inline constexpr uint64_t XXH3_avalanche(uint64_t h) noexcept { + h = (h ^ (h >> 37)) * 0x165667919E3779F9ULL; + return h ^ (h >> 32); +} + +inline constexpr uint64_t rrmxmx(uint64_t h, uint64_t len) noexcept { + h ^= ((h << 49) | (h >> 15)) ^ ((h << 24) | (h >> 40)); + h *= 0x9FB21C651E98DF25ULL; + h ^= (h >> 35) + len; + h *= 0x9FB21C651E98DF25ULL; + return h ^ (h >> 28); +} + +template +constexpr uint64_t mix16B(const T* input, const S* secret, + uint64_t seed) noexcept { + return mul128_fold64(readLE64(input) ^ (readLE64(secret) + seed), + readLE64(input + 8) ^ (readLE64(secret + 8) - seed)); +} + +inline constexpr size_t STRIPE_LEN = 64; +inline constexpr size_t SECRET_CONSUME_RATE = 8; +inline constexpr size_t ACC_NB = STRIPE_LEN / sizeof(uint64_t); + +template +constexpr void accumulate_512(uint64_t* acc, const T* input, + const S* secret) noexcept { + for (size_t i = 0; i < ACC_NB; i++) { + uint64_t data_val = readLE64(input + 8 * i); + uint64_t data_key = data_val ^ readLE64(secret + i * 8); + acc[i ^ 1] += data_val; + acc[i] += uint32_t(data_key) * (data_key >> 32); + } +} + +template +constexpr uint64_t hashLong_64b_internal(const T* input, size_t len, + const S* secret, + size_t secretSize) noexcept { + uint64_t acc[ACC_NB]{PRIME32_3, PRIME64_1, PRIME64_2, PRIME64_3, + PRIME64_4, PRIME32_2, PRIME64_5, PRIME32_1}; + size_t nbStripesPerBlock = (secretSize - STRIPE_LEN) / SECRET_CONSUME_RATE; + size_t block_len = STRIPE_LEN * nbStripesPerBlock; + size_t nb_blocks = (len - 1) / block_len; + + for (size_t n = 0; n < nb_blocks; n++) { + for (size_t i = 0; i < nbStripesPerBlock; i++) + accumulate_512(acc, input + n * block_len + i * STRIPE_LEN, + secret + i * SECRET_CONSUME_RATE); + for (size_t i = 0; i < ACC_NB; i++) + acc[i] = (acc[i] ^ (acc[i] >> 47) ^ + readLE64(secret + secretSize - STRIPE_LEN + 8 * i)) * + PRIME32_1; + } + + size_t nbStripes = ((len - 1) - (block_len * nb_blocks)) / STRIPE_LEN; + for (size_t i = 0; i < nbStripes; i++) + accumulate_512(acc, input + nb_blocks * block_len + i * STRIPE_LEN, + secret + i * SECRET_CONSUME_RATE); + accumulate_512(acc, input + len - STRIPE_LEN, + secret + secretSize - STRIPE_LEN - 7); + uint64_t result = len * PRIME64_1; + for (size_t i = 0; i < 4; i++) + result += + mul128_fold64(acc[2 * i] ^ readLE64(secret + 11 + 16 * i), + acc[2 * i + 1] ^ readLE64(secret + 11 + 16 * i + 8)); + return XXH3_avalanche(result); +} + +template +constexpr uint64_t XXH3_64bits_internal(const T* input, size_t len, + uint64_t seed, const S* secret, + size_t secretLen, + HashLong f_hashLong) noexcept { + if (len == 0) { + return XXH64_avalanche(seed ^ + (readLE64(secret + 56) ^ readLE64(secret + 64))); + } else if (len < 4) { + uint64_t keyed = ((uint32_t(uint8_t(input[0])) << 16) | + (uint32_t(uint8_t(input[len >> 1])) << 24) | + uint8_t(input[len - 1]) | (uint32_t(len) << 8)) ^ + ((readLE32(secret) ^ readLE32(secret + 4)) + seed); + return XXH64_avalanche(keyed); + } else if (len <= 8) { + uint64_t keyed = + (readLE32(input + len - 4) + (uint64_t(readLE32(input)) << 32)) ^ + ((readLE64(secret + 8) ^ readLE64(secret + 16)) - + (seed ^ (uint64_t(swap32(uint32_t(seed))) << 32))); + return rrmxmx(keyed, len); + } else if (len <= 16) { + uint64_t input_lo = + readLE64(input) ^ + ((readLE64(secret + 24) ^ readLE64(secret + 32)) + seed); + uint64_t input_hi = + readLE64(input + len - 8) ^ + ((readLE64(secret + 40) ^ readLE64(secret + 48)) - seed); + uint64_t acc = + len + swap64(input_lo) + input_hi + mul128_fold64(input_lo, input_hi); + return XXH3_avalanche(acc); + } else if (len <= 128) { + uint64_t acc = len * PRIME64_1; + size_t secret_off = 0; + for (size_t i = 0, j = len; j > i; i += 16, j -= 16) { + acc += mix16B(input + i, secret + secret_off, seed); + acc += mix16B(input + j - 16, secret + secret_off + 16, seed); + secret_off += 32; + } + return XXH3_avalanche(acc); + } else if (len <= 240) { + uint64_t acc = len * PRIME64_1; + for (size_t i = 0; i < 128; i += 16) + acc += mix16B(input + i, secret + i, seed); + acc = XXH3_avalanche(acc); + for (size_t i = 128; i < len / 16 * 16; i += 16) + acc += mix16B(input + i, secret + (i - 128) + 3, seed); + acc += mix16B(input + len - 16, secret + SECRET_SIZE_MIN - 17, seed); + return XXH3_avalanche(acc); + } else { + return f_hashLong(input, len, seed, secret, secretLen); + } +} + +template +constexpr size_t bytes_size(const Bytes& bytes) noexcept { + return std::size(bytes); +} + +template +constexpr size_t bytes_size(T (&)[N]) noexcept { + return (N ? N - 1 : 0); +} + +/// Basic interfaces + +template +consteval uint64_t XXH3_64bits_const(const T* input, size_t len) noexcept { + return XXH3_64bits_internal( + input, len, 0, kSecret, sizeof(kSecret), + [](const T* input, size_t len, uint64_t, const void*, + size_t) constexpr noexcept { + return hashLong_64b_internal(input, len, kSecret, sizeof(kSecret)); + }); +} + +template +consteval uint64_t XXH3_64bits_withSecret_const(const T* input, size_t len, + const S* secret, + size_t secretSize) noexcept { + return XXH3_64bits_internal( + input, len, 0, secret, secretSize, + [](const T* input, size_t len, uint64_t, const S* secret, + size_t secretLen) constexpr noexcept { + return hashLong_64b_internal(input, len, secret, secretLen); + }); +} + +template +consteval uint64_t XXH3_64bits_withSeed_const(const T* input, size_t len, + uint64_t seed) noexcept { + if (seed == 0) return XXH3_64bits_const(input, len); + return XXH3_64bits_internal( + input, len, seed, kSecret, sizeof(kSecret), + [](const T* input, size_t len, uint64_t seed, const void*, + size_t) constexpr noexcept { + uint8_t secret[SECRET_DEFAULT_SIZE]; + for (size_t i = 0; i < SECRET_DEFAULT_SIZE; i += 16) { + writeLE64(secret + i, readLE64(kSecret + i) + seed); + writeLE64(secret + i + 8, readLE64(kSecret + i + 8) - seed); + } + return hashLong_64b_internal(input, len, secret, sizeof(secret)); + }); +} + +/// Convenient interfaces + +template +consteval uint64_t XXH3_64bits_const(const Bytes& input) noexcept { + return XXH3_64bits_const(std::data(input), bytes_size(input)); +} + +template +consteval uint64_t XXH3_64bits_withSecret_const(const Bytes& input, + const Secret& secret) noexcept { + return XXH3_64bits_withSecret_const(std::data(input), bytes_size(input), + std::data(secret), bytes_size(secret)); +} + +template +consteval uint64_t XXH3_64bits_withSeed_const(const Bytes& input, + uint64_t seed) noexcept { + return XXH3_64bits_withSeed_const(std::data(input), bytes_size(input), seed); +} + +} // namespace constexpr_xxh3 diff --git a/build-config/xxHash/include/xxhash64.h b/build-config/xxHash/include/xxhash64.h new file mode 100644 index 00000000..4d0bbc5d --- /dev/null +++ b/build-config/xxHash/include/xxhash64.h @@ -0,0 +1,202 @@ +// ////////////////////////////////////////////////////////// +// xxhash64.h +// Copyright (c) 2016 Stephan Brumme. All rights reserved. +// see http://create.stephan-brumme.com/disclaimer.html +// + +#pragma once +#include // for uint32_t and uint64_t + +/// XXHash (64 bit), based on Yann Collet's descriptions, see http://cyan4973.github.io/xxHash/ +/** How to use: + uint64_t myseed = 0; + XXHash64 myhash(myseed); + myhash.add(pointerToSomeBytes, numberOfBytes); + myhash.add(pointerToSomeMoreBytes, numberOfMoreBytes); // call add() as often as you like to ... + // and compute hash: + uint64_t result = myhash.hash(); + + // or all of the above in one single line: + uint64_t result2 = XXHash64::hash(mypointer, numBytes, myseed); + + Note: my code is NOT endian-aware ! +**/ +class XXHash64 +{ +public: + /// create new XXHash (64 bit) + /** @param seed your seed value, even zero is a valid seed **/ + explicit XXHash64(uint64_t seed) + { + state[0] = seed + Prime1 + Prime2; + state[1] = seed + Prime2; + state[2] = seed; + state[3] = seed - Prime1; + bufferSize = 0; + totalLength = 0; + } + + /// add a chunk of bytes + /** @param input pointer to a continuous block of data + @param length number of bytes + @return false if parameters are invalid / zero **/ + bool add(const void* input, uint64_t length) + { + // no data ? + if (!input || length == 0) + return false; + + totalLength += length; + // byte-wise access + const unsigned char* data = (const unsigned char*)input; + + // unprocessed old data plus new data still fit in temporary buffer ? + if (bufferSize + length < MaxBufferSize) + { + // just add new data + while (length-- > 0) + buffer[bufferSize++] = *data++; + return true; + } + + // point beyond last byte + const unsigned char* stop = data + length; + const unsigned char* stopBlock = stop - MaxBufferSize; + + // some data left from previous update ? + if (bufferSize > 0) + { + // make sure temporary buffer is full (16 bytes) + while (bufferSize < MaxBufferSize) + buffer[bufferSize++] = *data++; + + // process these 32 bytes (4x8) + process(buffer, state[0], state[1], state[2], state[3]); + } + + // copying state to local variables helps optimizer A LOT + uint64_t s0 = state[0], s1 = state[1], s2 = state[2], s3 = state[3]; + // 32 bytes at once + while (data <= stopBlock) + { + // local variables s0..s3 instead of state[0]..state[3] are much faster + process(data, s0, s1, s2, s3); + data += 32; + } + // copy back + state[0] = s0; state[1] = s1; state[2] = s2; state[3] = s3; + + // copy remainder to temporary buffer + bufferSize = stop - data; + for (uint64_t i = 0; i < bufferSize; i++) + buffer[i] = data[i]; + + // done + return true; + } + + /// get current hash + /** @return 64 bit XXHash **/ + uint64_t hash() const + { + // fold 256 bit state into one single 64 bit value + uint64_t result; + if (totalLength >= MaxBufferSize) + { + result = rotateLeft(state[0], 1) + + rotateLeft(state[1], 7) + + rotateLeft(state[2], 12) + + rotateLeft(state[3], 18); + result = (result ^ processSingle(0, state[0])) * Prime1 + Prime4; + result = (result ^ processSingle(0, state[1])) * Prime1 + Prime4; + result = (result ^ processSingle(0, state[2])) * Prime1 + Prime4; + result = (result ^ processSingle(0, state[3])) * Prime1 + Prime4; + } + else + { + // internal state wasn't set in add(), therefore original seed is still stored in state2 + result = state[2] + Prime5; + } + + result += totalLength; + + // process remaining bytes in temporary buffer + const unsigned char* data = buffer; + // point beyond last byte + const unsigned char* stop = data + bufferSize; + + // at least 8 bytes left ? => eat 8 bytes per step + for (; data + 8 <= stop; data += 8) + result = rotateLeft(result ^ processSingle(0, *(uint64_t*)data), 27) * Prime1 + Prime4; + + // 4 bytes left ? => eat those + if (data + 4 <= stop) + { + result = rotateLeft(result ^ (*(uint32_t*)data) * Prime1, 23) * Prime2 + Prime3; + data += 4; + } + + // take care of remaining 0..3 bytes, eat 1 byte per step + while (data != stop) + result = rotateLeft(result ^ (*data++) * Prime5, 11) * Prime1; + + // mix bits + result ^= result >> 33; + result *= Prime2; + result ^= result >> 29; + result *= Prime3; + result ^= result >> 32; + return result; + } + + + /// combine constructor, add() and hash() in one static function (C style) + /** @param input pointer to a continuous block of data + @param length number of bytes + @param seed your seed value, e.g. zero is a valid seed + @return 64 bit XXHash **/ + static uint64_t hash(const void* input, uint64_t length, uint64_t seed) + { + XXHash64 hasher(seed); + hasher.add(input, length); + return hasher.hash(); + } + +private: + /// magic constants :-) + static const uint64_t Prime1 = 11400714785074694791ULL; + static const uint64_t Prime2 = 14029467366897019727ULL; + static const uint64_t Prime3 = 1609587929392839161ULL; + static const uint64_t Prime4 = 9650029242287828579ULL; + static const uint64_t Prime5 = 2870177450012600261ULL; + + /// temporarily store up to 31 bytes between multiple add() calls + static const uint64_t MaxBufferSize = 31+1; + + uint64_t state[4]; + unsigned char buffer[MaxBufferSize]; + uint64_t bufferSize; + uint64_t totalLength; + + /// rotate bits, should compile to a single CPU instruction (ROL) + static inline uint64_t rotateLeft(uint64_t x, unsigned char bits) + { + return (x << bits) | (x >> (64 - bits)); + } + + /// process a single 64 bit value + static inline uint64_t processSingle(uint64_t previous, uint64_t input) + { + return rotateLeft(previous + input * Prime2, 31) * Prime1; + } + + /// process a block of 4x4 bytes, this is the main part of the XXHash32 algorithm + static inline void process(const void* data, uint64_t& state0, uint64_t& state1, uint64_t& state2, uint64_t& state3) + { + const uint64_t* block = (const uint64_t*) data; + state0 = processSingle(state0, block[0]); + state1 = processSingle(state1, block[1]); + state2 = processSingle(state2, block[2]); + state3 = processSingle(state3, block[3]); + } +}; diff --git a/build-config/xxHash/meson.build b/build-config/xxHash/meson.build new file mode 100644 index 00000000..a8551b51 --- /dev/null +++ b/build-config/xxHash/meson.build @@ -0,0 +1,4 @@ + +xxhash_dep = declare_dependency( + include_directories: include_directories('include') +) \ No newline at end of file diff --git a/meson.build b/meson.build index 30f2a160..621a63d8 100644 --- a/meson.build +++ b/meson.build @@ -24,6 +24,31 @@ project('GridFire', 'cpp', version: '0.0.1a', default_options: ['cpp_std=c++23'] add_project_arguments('-fvisibility=default', language: 'cpp') # We disable these because of CppAD add_project_arguments('-Wno-bitwise-instead-of-logical', language: 'cpp') +llevel = get_option('log_level') + +logbase='QUILL_COMPILE_ACTIVE_LOG_LEVEL_' + +if (llevel == 'traceL3') + log_argument = logbase + 'TRACE_L3' +elif (llevel == 'traceL2') + log_argument = logbase + 'TRACE_L2' +elif (llevel == 'traceL1') + log_argument = logbase + 'TRACE_L1' +elif (llevel == 'debug') + log_argument = logbase + 'DEBUG' +elif (llevel == 'info') + log_argument = logbase + 'INFO' +elif (llevel == 'warning') + log_argument = logbase + 'WARNING' +elif (llevel == 'error') + log_argument = logbase + 'ERROR' +elif (llevel == 'critical') + log_argument = logbase + 'CRITICAL' +endif + +log_argument = '-DQUILL_COMPILE_ACTIVE_LOG_LEVEL=' + log_argument +add_project_arguments(log_argument, language: 'cpp') + cpp = meson.get_compiler('cpp') subdir('build-config') diff --git a/meson_options.txt b/meson_options.txt index f5bc8c17..efbdc46b 100644 --- a/meson_options.txt +++ b/meson_options.txt @@ -1,20 +1 @@ -# Build options for the project - -# Usage: -# -Doption_name=value - -# -Dbuild_tests=true (default) build tests -# -Duser_mode=false (default) enable user mode (set mode = 0) If user mode is enabled then the optimization level is set to 3 and the build type is set to release -# -Dbuild_python=true (default) build Python bindings -option('build_tests', type: 'boolean', value: true, description: 'Build tests') -option('user_mode', type: 'boolean', value: false, description: 'Enable user mode (set mode = 0)') -option('build_python', type: 'boolean', value: true, description: 'Build Python bindings') -option( - 'config_error_handling', - type: 'combo', - choices: [ 'none', 'warn', 'harsh' ], - value: 'none', - description: 'What to do if a config file fails to load: silent (none), warning (warn), or error (harsh)' -) -option('build_post_run_utils', type: 'boolean', value: true, description: 'Build Helper Utilities') -option('build_debug_utils', type: 'boolean', value: true, description: 'Build Debug Utilities') +option('log_level', type: 'combo', choices: ['traceL3', 'traceL2', 'traceL1', 'debug', 'info', 'warning', 'error', 'critial'], value: 'info', description: 'Set the log level for the GridFire library') \ No newline at end of file diff --git a/src/network/include/gridfire/engine/engine_abstract.h b/src/network/include/gridfire/engine/engine_abstract.h new file mode 100644 index 00000000..a62e7c88 --- /dev/null +++ b/src/network/include/gridfire/engine/engine_abstract.h @@ -0,0 +1,66 @@ +#pragma once + +#include "gridfire/network.h" // For NetIn, NetOut +#include "../reaction/reaction.h" + +#include "fourdst/composition/composition.h" +#include "fourdst/config/config.h" +#include "fourdst/logging/logging.h" + +#include +#include + +namespace gridfire { + + template + concept IsArithmeticOrAD = std::is_same_v || std::is_same_v>; + + template + struct StepDerivatives { + std::vector dydt; ///< Derivatives of abundances. + T nuclearEnergyGenerationRate = T(0.0); ///< Specific energy generation rate. + }; + + class Engine { + public: + virtual ~Engine() = default; + virtual const std::vector& getNetworkSpecies() const = 0; + virtual StepDerivatives calculateRHSAndEnergy( + const std::vector& Y, + double T9, + double rho + ) const = 0; + }; + + class DynamicEngine : public Engine { + public: + virtual void generateJacobianMatrix( + const std::vector& Y, + double T9, double rho + ) = 0; + virtual double getJacobianMatrixEntry( + int i, + int j + ) const = 0; + + virtual void generateStoichiometryMatrix() = 0; + virtual int getStoichiometryMatrixEntry( + int speciesIndex, + int reactionIndex + ) const = 0; + + virtual double calculateMolarReactionFlow( + const reaction::Reaction& reaction, + const std::vector& Y, + double T9, + double rho + ) const = 0; + + virtual const reaction::REACLIBLogicalReactionSet& getNetworkReactions() const = 0; + virtual std::unordered_map getSpeciesTimescales( + const std::vector& Y, + double T9, + double rho + ) const = 0; + }; +} \ No newline at end of file diff --git a/src/network/include/gridfire/approx8.h b/src/network/include/gridfire/engine/engine_approx8.h similarity index 100% rename from src/network/include/gridfire/approx8.h rename to src/network/include/gridfire/engine/engine_approx8.h diff --git a/src/network/include/gridfire/engine/engine_culled.h b/src/network/include/gridfire/engine/engine_culled.h new file mode 100644 index 00000000..c7d02826 --- /dev/null +++ b/src/network/include/gridfire/engine/engine_culled.h @@ -0,0 +1,55 @@ +#pragma once +#include "gridfire/engine/engine_abstract.h" + +#include "fourdst/composition/atomicSpecies.h" + +namespace gridfire { + class CulledEngine final : public DynamicEngine { + public: + explicit CulledEngine(DynamicEngine& baseEngine); + + const std::vector& getNetworkSpecies() const override; + StepDerivatives calculateRHSAndEnergy( + const std::vector &Y, + double T9, + double rho + ) const override; + + void generateJacobianMatrix( + const std::vector &Y, + double T9, + double rho + ) override; + double getJacobianMatrixEntry( + int i, + int j + ) const override; + + void generateStoichiometryMatrix() override; + int getStoichiometryMatrixEntry( + int speciesIndex, + int reactionIndex + ) const override; + + double calculateMolarReactionFlow( + const reaction::Reaction &reaction, + const std::vector &Y, + double T9, + double rho + ) const override; + + const reaction::REACLIBLogicalReactionSet& getNetworkReactions() const override; + std::unordered_map getSpeciesTimescales( + const std::vector &Y, + double T9, + double rho + ) const override; + private: + DynamicEngine& m_baseEngine; + std::unordered_map m_fullSpeciesToIndexMap; + std::vector m_culledSpecies; + private: + std::unordered_map populatedSpeciesToIndexMap; + std::vector determineCullableSpecies; + }; +} diff --git a/src/network/include/gridfire/engine/engine_graph.h b/src/network/include/gridfire/engine/engine_graph.h new file mode 100644 index 00000000..ac444186 --- /dev/null +++ b/src/network/include/gridfire/engine/engine_graph.h @@ -0,0 +1,277 @@ +#pragma once + +#include "fourdst/composition/atomicSpecies.h" +#include "fourdst/composition/composition.h" +#include "fourdst/logging/logging.h" +#include "fourdst/config/config.h" + +#include "gridfire/network.h" +#include "gridfire/reaction/reaction.h" +#include "gridfire/engine/engine_abstract.h" + +#include +#include +#include + +#include + +#include "cppad/cppad.hpp" + +// PERF: The function getNetReactionStoichiometry returns a map of species to their stoichiometric coefficients for a given reaction. +// this makes extra copies of the species, which is not ideal and could be optimized further. +// Even more relevant is the member m_reactionIDMap which makes copies of a REACLIBReaction for each reaction ID. +// REACLIBReactions are quite large data structures, so this could be a performance bottleneck. + +namespace gridfire { + typedef CppAD::AD ADDouble; ///< Alias for CppAD AD type for double precision. + + using fourdst::config::Config; + using fourdst::logging::LogManager; + using fourdst::constant::Constants; + + static constexpr double MIN_DENSITY_THRESHOLD = 1e-18; + static constexpr double MIN_ABUNDANCE_THRESHOLD = 1e-18; + static constexpr double MIN_JACOBIAN_THRESHOLD = 1e-24; + + class GraphEngine final : public DynamicEngine{ + public: + explicit GraphEngine(const fourdst::composition::Composition &composition); + explicit GraphEngine(reaction::REACLIBLogicalReactionSet reactions); + + StepDerivatives calculateRHSAndEnergy( + const std::vector& Y, + const double T9, + const double rho + ) const override; + + void generateJacobianMatrix( + const std::vector& Y, + const double T9, + const double rho + ) override; + + void generateStoichiometryMatrix() override; + + double calculateMolarReactionFlow( + const reaction::Reaction& reaction, + const std::vector&Y, + const double T9, + const double rho + ) const override; + + [[nodiscard]] const std::vector& getNetworkSpecies() const override; + [[nodiscard]] const reaction::REACLIBLogicalReactionSet& getNetworkReactions() const override; + [[nodiscard]] double getJacobianMatrixEntry( + const int i, + const int j + ) const override; + [[nodiscard]] std::unordered_map getNetReactionStoichiometry( + const reaction::Reaction& reaction + ) const; + [[nodiscard]] int getStoichiometryMatrixEntry( + const int speciesIndex, + const int reactionIndex + ) const override; + [[nodiscard]] std::unordered_map getSpeciesTimescales( + const std::vector& Y, + double T9, + double rho + ) const override; + + [[nodiscard]] bool involvesSpecies( + const fourdst::atomic::Species& species + ) const; + + void exportToDot( + const std::string& filename + ) const; + void exportToCSV( + const std::string& filename + ) const; + + + private: + reaction::REACLIBLogicalReactionSet m_reactions; ///< Set of REACLIB reactions in the network. + std::unordered_map m_reactionIDMap; ///< Map from reaction ID to REACLIBReaction. //PERF: This makes copies of REACLIBReaction and could be a performance bottleneck. + + std::vector m_networkSpecies; ///< Vector of unique species in the network. + std::unordered_map m_networkSpeciesMap; ///< Map from species name to Species object. + std::unordered_map m_speciesToIndexMap; ///< Map from species to their index in the stoichiometry matrix. + + boost::numeric::ublas::compressed_matrix m_stoichiometryMatrix; ///< Stoichiometry matrix (species x reactions). + boost::numeric::ublas::compressed_matrix m_jacobianMatrix; ///< Jacobian matrix (species x species). + + CppAD::ADFun m_rhsADFun; ///< CppAD function for the right-hand side of the ODE. + + Config& m_config = Config::getInstance(); + Constants& m_constants = Constants::getInstance(); ///< Access to physical constants. + quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); + + private: + void syncInternalMaps(); + void collectNetworkSpecies(); + void populateReactionIDMap(); + void populateSpeciesToIndexMap(); + void reserveJacobianMatrix(); + void recordADTape(); + + [[nodiscard]] bool validateConservation() const; + void validateComposition( + const fourdst::composition::Composition &composition, + double culling, + double T9 + ); + + template + T calculateMolarReactionFlow( + const reaction::Reaction &reaction, + const std::vector &Y, + const T T9, + const T rho + ) const; + + template + StepDerivatives calculateAllDerivatives( + const std::vector &Y_in, + T T9, + T rho + ) const; + + StepDerivatives calculateAllDerivatives( + const std::vector& Y_in, + const double T9, + const double rho + ) const; + + StepDerivatives calculateAllDerivatives( + const std::vector& Y_in, + const ADDouble T9, + const ADDouble rho + ) const; + }; + + + template + StepDerivatives GraphEngine::calculateAllDerivatives( + const std::vector &Y_in, T T9, T rho) const { + + // --- Setup output derivatives structure --- + StepDerivatives result; + result.dydt.resize(m_networkSpecies.size(), static_cast(0.0)); + + // --- AD Pre-setup (flags to control conditionals in an AD safe / branch aware manner) --- + // ----- Constants for AD safe calculations --- + const T zero = static_cast(0.0); + const T one = static_cast(1.0); + + // ----- Initialize variables for molar concentration product and thresholds --- + // Note: the logic here is that we use CppAD::CondExprLt to test thresholds and if they are less we set the flag + // to zero so that the final returned reaction flow is 0. This is as opposed to standard if statements + // which create branches that break the AD tape. + const T rho_threshold = static_cast(MIN_DENSITY_THRESHOLD); + + // --- Check if the density is below the threshold where we ignore reactions --- + T threshold_flag = CppAD::CondExpLt(rho, rho_threshold, zero, one); // If rho < threshold, set flag to 0 + + std::vector Y = Y_in; + for (size_t i = 0; i < m_networkSpecies.size(); ++i) { + // We use CppAD::CondExpLt to handle AD taping and prevent branching + // Note that while this is syntactically more complex this is equivalent to + // if (Y[i] < 0) {Y[i] = 0;} + // The issue is that this would introduce a branch which would require the auto diff tape to be re-recorded + // each timestep, which is very inefficient. + Y[i] = CppAD::CondExpLt(Y[i], zero, zero, Y[i]); // Ensure no negative abundances + } + + const T u = static_cast(m_constants.get("u").value); // Atomic mass unit in grams + const T N_A = static_cast(m_constants.get("N_a").value); // Avogadro's number in mol^-1 + const T c = static_cast(m_constants.get("c").value); // Speed of light in cm/s + + // --- SINGLE LOOP OVER ALL REACTIONS --- + for (size_t reactionIndex = 0; reactionIndex < m_reactions.size(); ++reactionIndex) { + const auto& reaction = m_reactions[reactionIndex]; + + // 1. Calculate reaction rate + const T molarReactionFlow = calculateMolarReactionFlow(reaction, Y, T9, rho); + + // 2. Use the rate to update all relevant species derivatives (dY/dt) + for (size_t speciesIndex = 0; speciesIndex < m_networkSpecies.size(); ++speciesIndex) { + const T nu_ij = static_cast(m_stoichiometryMatrix(speciesIndex, reactionIndex)); + result.dydt[speciesIndex] += threshold_flag * nu_ij * molarReactionFlow / rho; + } + } + + T massProductionRate = static_cast(0.0); // [mol][s^-1] + for (const auto& [species, index] : m_speciesToIndexMap) { + massProductionRate += result.dydt[index] * species.mass() * u; + } + + result.nuclearEnergyGenerationRate = -massProductionRate * N_A * c * c; // [cm^2][s^-3] = [erg][s^-1][g^-1] + + return result; + } + + + template + T GraphEngine::calculateMolarReactionFlow( + const reaction::Reaction &reaction, + const std::vector &Y, + const T T9, + const T rho + ) const { + + // --- Pre-setup (flags to control conditionals in an AD safe / branch aware manner) --- + // ----- Constants for AD safe calculations --- + const T zero = static_cast(0.0); + const T one = static_cast(1.0); + + // ----- Initialize variables for molar concentration product and thresholds --- + // Note: the logic here is that we use CppAD::CondExprLt to test thresholds and if they are less we set the flag + // to zero so that the final returned reaction flow is 0. This is as opposed to standard if statements + // which create branches that break the AD tape. + const T Y_threshold = static_cast(MIN_ABUNDANCE_THRESHOLD); + T threshold_flag = one; + + // --- Calculate the molar reaction rate (in units of [s^-1][cm^3(N-1)][mol^(1-N)] for N reactants) --- + const T k_reaction = reaction.calculate_rate(T9); + + // --- Cound the number of each reactant species to account for species multiplicity --- + std::unordered_map reactant_counts; + reactant_counts.reserve(reaction.reactants().size()); + for (const auto& reactant : reaction.reactants()) { + reactant_counts[std::string(reactant.name())]++; + } + + // --- Accumulator for the molar concentration --- + auto molar_concentration_product = static_cast(1.0); + + // --- Loop through each unique reactant species and calculate the molar concentration for that species then multiply that into the accumulator --- + for (const auto& [species_name, count] : reactant_counts) { + // --- Resolve species to molar abundance --- + // PERF: Could probably optimize out this lookup + const auto species_it = m_speciesToIndexMap.find(m_networkSpeciesMap.at(species_name)); + const size_t species_index = species_it->second; + const T Yi = Y[species_index]; + + // --- Check if the species abundance is below the threshold where we ignore reactions --- + threshold_flag *= CppAD::CondExpLt(Yi, Y_threshold, zero, one); + + // --- Convert from molar abundance to molar concentration --- + T molar_concentration = Yi * rho; + + // --- If count is > 1 , we need to raise the molar concentration to the power of count since there are really count bodies in that reaction --- + molar_concentration_product *= CppAD::pow(molar_concentration, static_cast(count)); // ni^count + + // --- Apply factorial correction for identical reactions --- + if (count > 1) { + molar_concentration_product /= static_cast(std::tgamma(static_cast(count + 1))); // Gamma function for factorial + } + } + // --- Final reaction flow calculation [mol][s^-1][cm^-3] --- + // Note: If the threshold flag ever gets set to zero this will return zero. + // This will result basically in multiple branches being written to the AD tape, which will make + // the tape more expensive to record, but it will also mean that we only need to record it once for + // the entire network. + return molar_concentration_product * k_reaction * threshold_flag; + } +}; \ No newline at end of file diff --git a/src/network/include/gridfire/netgraph.h b/src/network/include/gridfire/netgraph.h deleted file mode 100644 index 5d5aa438..00000000 --- a/src/network/include/gridfire/netgraph.h +++ /dev/null @@ -1,571 +0,0 @@ -#pragma once - -#include "fourdst/composition/atomicSpecies.h" -#include "fourdst/composition/composition.h" - -#include "gridfire/network.h" -#include "gridfire/reaclib.h" - -#include -#include -#include - -#include - -#include "cppad/cppad.hpp" - -#include "quill/LogMacros.h" - -// PERF: The function getNetReactionStoichiometry returns a map of species to their stoichiometric coefficients for a given reaction. -// this makes extra copies of the species, which is not ideal and could be optimized further. -// Even more relevant is the member m_reactionIDMap which makes copies of a REACLIBReaction for each reaction ID. -// REACLIBReactions are quite large data structures, so this could be a performance bottleneck. - -/** - * @file netgraph.h - * @brief Defines the GraphNetwork class for representing and evolving nuclear reaction networks as a heavily connected graph. - * - * This file provides the GraphNetwork class, which implements a graph-based nuclear reaction network - * using REACLIB reactions and supports both stiff and non-stiff ODE integration. The network is constructed - * from a composition and can be queried for species, reactions, and stoichiometry. - * - * This is a general nuclear reaction network; however, it does not currently manage reverse reactions, weak reactions, - * or other reactions which become much more relevant for more extreme astrophysical sources. - * - * @see gridfire::GraphNetwork - */ - -namespace gridfire { - - /** - * @brief Concept to check if a type is either double or CppAD::AD. - * - * Used to enable templated functions for both arithmetic and automatic differentiation types. - */ - template - concept IsArithmeticOrAD = std::is_same_v || std::is_same_v>; - - /// Minimum density threshold (g/cm^3) below which reactions are ignored. - static constexpr double MIN_DENSITY_THRESHOLD = 1e-18; - /// Minimum abundance threshold below which reactions are ignored. - static constexpr double MIN_ABUNDANCE_THRESHOLD = 1e-18; - /// Minimum Jacobian entry threshold for sparsity. - static constexpr double MIN_JACOBIAN_THRESHOLD = 1e-24; - - /** - * @brief Graph-based nuclear reaction network using REACLIB reactions. - * - * GraphNetwork constructs a reaction network from a given composition, builds the associated - * stoichiometry and Jacobian matrices, and provides methods for evaluating the network's evolution - * using ODE solvers. It supports both stiff and non-stiff integration and can be queried for - * network species, reactions, and stoichiometry. - * - * @note Currently does not handle reverse reactions, weak reactions, or other complex reaction types. These will be added in future versions. - * - * Example usage: - * @code - * serif::composition::Composition comp = ...; - * serif::network::GraphNetwork net(comp); - * serif::network::NetIn input; - * input.composition = comp; - * input.temperature = 1.5e9; - * input.density = 1e6; - * input.tMax = 1.0; - * input.dt0 = 1e-3; - * serif::network::NetOut result = net.evaluate(input); - * @endcode - */ - class GraphNetwork final : public Network { - public: - /** - * @brief Construct a GraphNetwork from a composition. - * @param composition The composition specifying the initial isotopic abundances. - */ - explicit GraphNetwork(const fourdst::composition::Composition &composition); - - /** - * @brief Construct a GraphNetwork from a composition with culling and temperature. - * @param composition The composition specifying the initial isotopic abundances. - * @param cullingThreshold specific reaction rate threshold to cull reactions below. - * @param T9 Temperature in units of 10^9 K where culling is evaluated at. - * - * @see serif::network::build_reaclib_nuclear_network - */ - explicit GraphNetwork(const fourdst::composition::Composition &composition, - const double cullingThreshold, const double T9); - - explicit GraphNetwork(const reaclib::REACLIBReactionSet& reactions); - - /** - * @brief Evolve the network for the given input conditions. - * - * This is the primary entry point for users of GridFire. This function integrates the network ODEs using either a stiff or non-stiff solver, - * depending on the detected stiffness of the system. - * - * @param netIn The input structure containing composition, temperature, density, and integration parameters. - * @return NetOut The output structure containing the final composition, number of steps, and energy. - * - * Example: - * @code - * serif::network::NetIn input; - * // ... set up input ... - * serif::network::NetOut result = net.evaluate(input); - * @endcode - */ - NetOut evaluate(const NetIn &netIn) override; - - /** - * @brief Get the vector of unique species in the network. - * @return Reference to the vector of species. - * - * Example: - * @code - * const auto& species = net.getNetworkSpecies(); - * for (const auto& sp : species) { - * std::cout << sp.name() << std::endl; - * } - * @endcode - */ - [[nodiscard]] const std::vector& getNetworkSpecies() const; - - /** - * @brief Get the set of REACLIB reactions in the network. - * @return Reference to the set of reactions. - */ - [[nodiscard]] const reaclib::REACLIBReactionSet& getNetworkReactions() const; - - /** - * @brief Get the net stoichiometric coefficients for all species in a reaction. - * - * Returns a map from species to their net stoichiometric coefficient (products minus reactants). - * - * @param reaction The REACLIB reaction to analyze. - * @return Map from species to stoichiometric coefficient. - * - * @throws std::runtime_error If a species in the reaction is not found in the network. - * - * Example: - * @code - * for (const auto& reaction : net.getNetworkReactions()) { - * auto stoichiometry = net.getNetReactionStoichiometry(reaction); - * // ... - * } - * @endcode - */ - [[nodiscard]] std::unordered_map getNetReactionStoichiometry( - const reaclib::REACLIBReaction &reaction) const; - - /** - * @brief Check if a species is present in the network. - * @param species The species to check. - * @return True if the species is present, false otherwise. - * - * Example: - * @code - * if (net.involvesSpecies(mySpecies)) { ... } - * @endcode - */ - [[nodiscard]] bool involvesSpecies(const fourdst::atomic::Species& species) const; - - /** - * @brief Detect cycles in the reaction network (not implemented). - * @return Vector of cycles, each represented as a vector of reaction IDs. - * - * @note This is not yet implemented; however, it will allow for easy detection of things like the CNO cycle. - * @deprecated Not implemented. - */ - [[deprecated("not implimented")]] [[nodiscard]] std::vector> detectCycles() const = delete; - - /** - * @brief Perform a topological sort of the reactions (not implemented). - * @return Vector of reaction IDs in topologically sorted order. - * @deprecated Not implemented. - */ - [[deprecated("not implimented")]] [[nodiscard]] std::vector topologicalSortReactions() const = delete; - - /** - * @brief Export the network to a DOT file for visualization. - * @param filename The name of the output DOT file. - */ - void exportToDot(const std::string& filename) const; - - private: - reaclib::REACLIBReactionSet m_reactions; ///< Set of REACLIB reactions in the network. - std::unordered_map m_reactionIDMap; ///< Map from reaction ID to REACLIBReaction. //PERF: This makes copies of REACLIBReaction and could be a performance bottleneck. - - std::vector m_networkSpecies; ///< Vector of unique species in the network. - std::unordered_map m_networkSpeciesMap; ///< Map from species name to Species object. - std::unordered_map m_speciesToIndexMap; ///< Map from species to their index in the stoichiometry matrix. - - boost::numeric::ublas::compressed_matrix m_stoichiometryMatrix; ///< Stoichiometry matrix (species x reactions). - boost::numeric::ublas::compressed_matrix m_jacobianMatrix; ///< Jacobian matrix (species x species). - - CppAD::ADFun m_rhsADFun; ///< CppAD function for the right-hand side of the ODE. - - /** - * @brief Functor for ODE right-hand side evaluation. - * - * Used by ODE solvers to compute dY/dt and energy generation rate. This is the only functor used in the non-NSE case. - */ - struct ODETerm { - GraphNetwork& m_network; ///< Reference to the network - const double m_T9; ///< Temperature in 10^9 K - const double m_rho; ///< Density in g/cm^3 - const size_t m_numSpecies; ///< Number of species - - /** - * @brief Construct an ODETerm functor. - * @param network Reference to the GraphNetwork. - * @param T9 Temperature in 10^9 K. - * @param rho Density in g/cm^3. - */ - ODETerm(GraphNetwork& network, const double T9, double rho) : - m_network(network), m_T9(T9), m_rho(rho), m_numSpecies(m_network.m_networkSpecies.size()) {} - - /** - * @brief Compute dY/dt and energy rate for the ODE solver. - * @param Y Input vector of abundances. - * @param dYdt Output vector for derivatives (resized). - * - * @note this functor does not need auto differentiation to it called the version of calculateAllDerivatives. - */ - void operator()(const boost::numeric::ublas::vector&Y, boost::numeric::ublas::vector& dYdt, double /* t */) const { - const std::vector y(Y.begin(), m_numSpecies + Y.begin()); - - StepDerivatives derivatives = m_network.calculateAllDerivatives(y, m_T9, m_rho); - - dYdt.resize(m_numSpecies + 1); - std::ranges::copy(derivatives.dydt, dYdt.begin()); - dYdt(m_numSpecies) = derivatives.specificEnergyRate; // Last element is the specific energy rate - } - }; - - /** - * @brief Functor for Jacobian evaluation for stiff ODE solvers. (used in the NSE case) - */ - struct JacobianTerm { - GraphNetwork& m_network; ///< Reference to the network - const double m_T9; ///< Temperature in 10^9 K - const double m_rho; ///< Density in g/cm^3 - const size_t m_numSpecies; ///< Number of species - - /** - * @brief Construct a JacobianTerm functor. - * @param network Reference to the GraphNetwork. - * @param T9 Temperature in 10^9 K. - * @param rho Density in g/cm^3. - */ - JacobianTerm(GraphNetwork& network, const double T9, const double rho) : - m_network(network), m_T9(T9), m_rho(rho), m_numSpecies(m_network.m_networkSpecies.size()) {} - - /** - * @brief Compute the Jacobian matrix for the ODE solver. - * @param Y Input vector of abundances. - * @param J Output matrix for the Jacobian (resized). - */ - void operator()(const boost::numeric::ublas::vector& Y, boost::numeric::ublas::matrix& J, double /* t */, boost::numeric::ublas::vector& /* dfdt */) const { - const std::vector y_species(Y.begin(), Y.begin() + m_numSpecies); - - m_network.generateJacobianMatrix(y_species, m_T9, m_rho); - - J.resize(m_numSpecies + 1, m_numSpecies + 1); - J.clear(); // Zero out the matrix - - // PERF: Could probably be optimized - for (auto it1 = m_network.m_jacobianMatrix.begin1(); it1 != m_network.m_jacobianMatrix.end1(); ++it1) { - for (auto it2 = it1.begin(); it2 != it1.end(); ++it2) { - J(it2.index1(), it2.index2()) = *it2; - } - } - } - }; - - /** - * @brief Struct holding derivatives for a single ODE step. - * @tparam T Scalar type (double or CppAD::AD). - */ - template - struct StepDerivatives { - std::vector dydt; ///< Derivatives of abundances. - T specificEnergyRate = T(0.0); ///< Specific energy generation rate. - }; - - private: - /** - * @brief Synchronize all internal maps and matrices after network changes. - * - * Called after the reaction set is updated to rebuild all derived data structures. - * - * @note This method must be called any time the network topology changes. - */ - void syncInternalMaps(); - - /** - * @brief Collect all unique species from the reaction set. - * - * Populates m_networkSpecies and m_networkSpeciesMap. - * @throws std::runtime_error If a species is not found in the global atomic species database. - * - * @note Called by syncInternalMaps() to ensure the species list is up-to-date. - */ - void collectNetworkSpecies(); - - /** - * @brief Populate the reaction ID map from the reaction set. - * - * @note Called by syncInternalMaps() to ensure the reaction ID map is up-to-date. - */ - void populateReactionIDMap(); - - /** - * @brief Populate the species-to-index map for matrix construction. - * - * @note Called by syncInternalMaps() to ensure the species-to-index map is up-to-date. - */ - void populateSpeciesToIndexMap(); - - /** - * @brief Reserve and resize the Jacobian matrix. - * - * @note Called by syncInternalMaps() to ensure the Jacobian matrix is ready for use. - */ - void reserveJacobianMatrix(); - - /** - * @brief Record the CppAD tape for automatic differentiation. - * @throws std::runtime_error If there are no species in the network. - * - * @note Called by syncInternalMaps() to ensure the AD tape is recorded for the current network state. - */ - void recordADTape(); - - // --- Validation methods --- - - /** - * @brief Validate mass and charge conservation for all reactions. - * @return True if all reactions conserve mass and charge, false otherwise. - * - * @note Logs errors for any violations. - */ - [[nodiscard]] bool validateConservation() const; - - /** - * @brief Validate and update the network composition if needed. - * - * If the composition or culling/temperature has changed, rebuilds the reaction set and synchronizes maps. - * This is primarily used to enable caching in dynamic network situations where the composition, temperature, and density - * may change but the relevant reaction set remains equivalent. - * - * @param composition The composition to validate. - * @param culling Culling threshold. - * @param T9 Temperature in 10^9 K. - */ - void validateComposition(const fourdst::composition::Composition &composition, double culling, double T9); - - // --- Simple Derivative Calculations --- - - /** - * @brief Calculate all derivatives (dY/dt and energy rate) for the current state. - * @tparam T Scalar type (double or CppAD::AD). - * @param Y Vector of abundances. - * @param T9 Temperature in 10^9 K. - * @param rho Density in g/cm^3. - * @return StepDerivatives containing dY/dt and energy rate. - */ - template - StepDerivatives calculateAllDerivatives(const std::vector& Y, T T9, T rho) const; - - // --- Generate the system matrices --- - - /** - * @brief Generate the stoichiometry matrix for the network. - * - * Populates m_stoichiometryMatrix based on the current reactions and species. - * @throws std::runtime_error If a species is not found in the species-to-index map. - */ - void generateStoichiometryMatrix(); - - /** - * @brief Generate the Jacobian matrix for the network. - * - * Populates m_jacobianMatrix using automatic differentiation via the precomputed CppAD tape. - * - * @param Y Vector of abundances. - * @param T9 Temperature in 10^9 K. - * @param rho Density in g/cm^3. - */ - void generateJacobianMatrix(const std::vector& Y, double T9, const double rho); - - /** - * @brief Calculate the right-hand side (dY/dt) for the ODE system. - * @tparam GeneralScalarType Scalar type (double or CppAD::AD). - * @param Y Vector of abundances. - * @param T9 Temperature in 10^9 K. - * @param rho Density in g/cm^3. - * @return Vector of dY/dt values. - */ - template - std::vector calculateRHS(const std::vector &Y, const GeneralScalarType T9, - GeneralScalarType rho) const; - - /** - * @brief Calculate the reaction rate for a given reaction in dimensions of particles per cm^3 per second. - * @tparam GeneralScalarType Scalar type (double or CppAD::AD). - * @param reaction The REACLIB reaction. - * @param Y Vector of abundances. - * @param T9 Temperature in 10^9 K. - * @param rho Density in g/cm^3. - * @return Reaction rate. - * @throws std::runtime_error If a reactant species is not found in the species-to-index map. - * - * @note reaclib uses molar units that vary depending on the number of reactants, It's pretty straightforward - * to convert these into particle based units. Specifically, we just need to divide the final result by - * Avogadro's number raised to the number of reactants - 1; - */ - template - GeneralScalarType calculateReactionRate(const reaclib::REACLIBReaction &reaction, - const std::vector &Y, const GeneralScalarType T9, - const GeneralScalarType rho) const; - - /** - * @brief Detect if the network is stiff and select the appropriate ODE solver. - * - * Heuristically determines stiffness based on the ratio of timescales. The stiffness heuristic is as follows: - * - * 1. For each species, compute the timescale as |Y_i / (dY_i/dt)|, where Y_i is the abundance and dY_i/dt is its time derivative. - * 2. Find the minimum and maximum timescales across all species. - * 3. Compute the stiffness ratio as (maximum timescale) / (minimum timescale). - * 4. If the stiffness ratio exceeds a threshold (default: 1e6), the system is considered stiff and a stiff ODE solver is used. - * Otherwise, a non-stiff ODE solver is used. - * - * This heuristic is based on the observation that stiff systems have widely varying timescales, which can cause explicit - * solvers to become unstable or inefficient. The threshold can be tuned based on the characteristics of the network. - * - * @param netIn The input structure. - * @param T9 Temperature in 10^9 K. - * @param numSpecies Number of species. - * @param Y Vector of abundances. - */ - void detectStiff(const NetIn &netIn, double T9, unsigned long numSpecies, const boost::numeric::ublas::vector& Y); - - }; - - - template - GraphNetwork::StepDerivatives GraphNetwork::calculateAllDerivatives( - const std::vector &Y, T T9, T rho) const { - StepDerivatives result; - result.dydt.resize(m_networkSpecies.size(), static_cast(0.0)); - - if (rho < MIN_DENSITY_THRESHOLD) { - return result; // Return zero for everything if density is too low - } - - const T u = static_cast(m_constants.get("u").value); // Atomic mass unit in grams - const T MeV_to_erg = static_cast(m_constants.get("MeV_to_erg").value); - - T volumetricEnergyRate = static_cast(0.0); // Accumulator for erg / (cm^3 * s) - - // --- SINGLE LOOP OVER ALL REACTIONS --- - for (size_t reactionIndex = 0; reactionIndex < m_reactions.size(); ++reactionIndex) { - const auto& reaction = m_reactions[reactionIndex]; - - // 1. Calculate reaction rate - const T reactionRate = calculateReactionRate(reaction, Y, T9, rho); - - // 2. Use the rate to update all relevant species derivatives (dY/dt) - for (size_t speciesIndex = 0; speciesIndex < m_networkSpecies.size(); ++speciesIndex) { - const T nu_ij = static_cast(m_stoichiometryMatrix(speciesIndex, reactionIndex)); - if (nu_ij != 0) { - const T speciesAtomicMassAMU = static_cast(m_networkSpecies[speciesIndex].mass()); - const T speciesAtomicMassGrams = speciesAtomicMassAMU * u; - result.dydt[speciesIndex] += (nu_ij * reactionRate * speciesAtomicMassGrams) / rho; - } - } - - // 3. Use the same rate to update the energy generation rate - const T q_value_ergs = static_cast(reaction.qValue()) * MeV_to_erg; - volumetricEnergyRate += reactionRate * q_value_ergs; - } - - result.specificEnergyRate = volumetricEnergyRate / rho; - - return result; - } - - - template - GeneralScalarType GraphNetwork::calculateReactionRate(const reaclib::REACLIBReaction &reaction, const std::vector &Y, - const GeneralScalarType T9, const GeneralScalarType rho) const { - const auto &constants = fourdst::constant::Constants::getInstance(); - - const auto u = constants.get("u"); // Atomic mass unit in g/mol - const auto uValue = static_cast(u.value); // Convert to double for calculations - const GeneralScalarType k_reaction = reaction.calculate_rate(T9); - - auto reactant_product_or_particle_densities = static_cast(1.0); - - std::unordered_map reactant_counts; - reactant_counts.reserve(reaction.reactants().size()); - for (const auto& reactant : reaction.reactants()) { - reactant_counts[std::string(reactant.name())]++; - } - - const auto minAbundanceThreshold = static_cast(MIN_ABUNDANCE_THRESHOLD); - const auto minDensityThreshold = static_cast(MIN_DENSITY_THRESHOLD); - - if (rho < minDensityThreshold) { - // LOG_INFO(m_logger, "Density is below the minimum threshold ({} g/cm^3), returning zero reaction rate for reaction '{}'.", - // CppAD::Value(rho), reaction.id()); // CppAD::Value causes a compilation error here, not clear why... - return static_cast(0.0); // If density is below a threshold, return zero rate. - } - - for (const auto& [species_name, count] : reactant_counts) { - auto species_it = m_speciesToIndexMap.find(m_networkSpeciesMap.at(species_name)); - if (species_it == m_speciesToIndexMap.end()) { - LOG_ERROR(m_logger, "Reactant species '{}' not found in species to index map for reaction '{}'.", - species_name, reaction.id()); - throw std::runtime_error("Reactant species not found in species to index map: " + species_name); - } - const size_t species_index = species_it->second; - const GeneralScalarType Yi = Y[species_index]; - - if (Yi < minAbundanceThreshold) { - return static_cast(0.0); // If any reactant is below a threshold, return zero rate. - } - - auto atomicMassAMU = static_cast(m_networkSpecies[species_index].mass()); - - // Convert to number density - GeneralScalarType ni; - const GeneralScalarType denominator = atomicMassAMU * uValue; - assert (denominator > 0.0); - ni = (Yi * rho) / (denominator); - - reactant_product_or_particle_densities *= ni; - - // Apply factorial correction for identical reactions - if (count > 1) { - reactant_product_or_particle_densities /= static_cast(std::tgamma(static_cast(count + 1))); // Gamma function for factorial - } - } - const GeneralScalarType Na = static_cast(constants.get("N_a").value); // Avogadro's number in mol^-1 - const int numReactants = static_cast(reaction.reactants().size()); - auto molarCorrectionFactor = static_cast(1.0); // No correction needed for single reactant reactions - if (numReactants > 1) { - molarCorrectionFactor = CppAD::pow(Na, static_cast(reaction.reactants().size() - 1)); - } - return (reactant_product_or_particle_densities * k_reaction) / molarCorrectionFactor; // reaction rate in per volume per time (particles/cm^3/s) - } - - template - std::vector GraphNetwork::calculateRHS( - const std::vector& Y, - const GeneralScalarType T9, - const GeneralScalarType rho) const { - - auto derivatives = calculateAllDerivatives(Y, T9, rho); - return derivatives.dydt; - } - -}; \ No newline at end of file diff --git a/src/network/include/gridfire/network.h b/src/network/include/gridfire/network.h index 057a9066..a9d56283 100644 --- a/src/network/include/gridfire/network.h +++ b/src/network/include/gridfire/network.h @@ -25,12 +25,15 @@ #include "fourdst/logging/logging.h" #include "fourdst/config/config.h" #include "fourdst/composition/species.h" -#include "quill/Logger.h" #include "fourdst/composition/composition.h" -#include "gridfire/reaclib.h" +#include "fourdst/constants/const.h" + +#include "gridfire/reaction/reaction.h" + +#include "quill/Logger.h" + #include -#include "fourdst/constants/const.h" namespace gridfire { @@ -127,7 +130,7 @@ namespace gridfire { * @param netIn Input parameters for the network evaluation. * @return NetOut Output results from the network evaluation. */ - virtual NetOut evaluate(const NetIn &netIn); + virtual NetOut evaluate(const NetIn &netIn) = 0; virtual bool isStiff() const { return m_stiff; } virtual void setStiff(const bool stiff) { m_stiff = stiff; } @@ -143,105 +146,9 @@ namespace gridfire { bool m_stiff = false; ///< Flag indicating if the network is stiff }; - class LogicalReaction { - public: - explicit LogicalReaction(const std::vector &reactions); - explicit LogicalReaction(const reaclib::REACLIBReaction &reaction); - void add_reaction(const reaclib::REACLIBReaction& reaction); - template - [[nodiscard]] T calculate_rate(const T T9) const { - T sum = static_cast(0.0); - const T T913 = CppAD::pow(T9, 1.0/3.0); - const T T953 = CppAD::pow(T9, 5.0/3.0); - const T logT9 = CppAD::log(T9); - for (const auto& rate : m_rates) { - const T exponent = rate.a0 + - rate.a1 / T9 + - rate.a2 / T913 + - rate.a3 * T913 + - rate.a4 * T9 + - rate.a5 * T953 + - rate.a6 * logT9; - sum += CppAD::exp(exponent); - } - return sum; - } - - [[nodiscard]] const std::string& id() const {return std::string(m_peID); } - - [[nodiscard]] std::string_view peName() const { return m_peID; } - - [[nodiscard]] int chapter() const { return m_chapter; } - - [[nodiscard]] const std::vector& reactants() const { return m_reactants; } - - [[nodiscard]] const std::vector& products() const { return m_products; } - - [[nodiscard]] double qValue() const { return m_qValue; } - - [[nodiscard]] bool is_reverse() const { return m_reverse; } - - - auto begin(); - auto begin() const; - auto end(); - auto end() const; - - private: - const quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); - std::string_view m_peID; - std::vector m_reactants; ///< Reactants of the reaction - std::vector m_products; ///< Products of the reaction - double m_qValue = 0.0; ///< Q-value of the reaction - bool m_reverse = false; ///< True if the reaction is reverse - int m_chapter; - - std::vector m_rates; - - }; - - class LogicalReactionSet { - public: - LogicalReactionSet() = default; - explicit LogicalReactionSet(const std::vector& reactions); - explicit LogicalReactionSet(const std::vector& reactions); - explicit LogicalReactionSet(const reaclib::REACLIBReactionSet& reactionSet); - - void add_reaction(const LogicalReaction& reaction); - void add_reaction(const reaclib::REACLIBReaction& reaction); - - void remove_reaction(const LogicalReaction& reaction); - - [[nodiscard]] bool contains(const std::string_view& id) const; - [[nodiscard]] bool contains(const LogicalReaction& reactions) const; - [[nodiscard]] bool contains(const reaclib::REACLIBReaction& reaction) const; - - [[nodiscard]] size_t size() const; - - void sort(double T9=1.0); - - bool contains_species(const fourdst::atomic::Species &species) const; - bool contains_reactant(const fourdst::atomic::Species &species) const; - bool contains_product(const fourdst::atomic::Species &species) const; - - [[nodiscard]] const LogicalReaction& operator[](size_t index) const; - [[nodiscard]] const LogicalReaction& operator[](const std::string_view& id) const; - - auto begin(); - auto begin() const; - auto end(); - auto end() const; - - - private: - const quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); - std::vector m_reactions; - std::unordered_map m_reactionNameMap; - }; - - LogicalReactionSet build_reaclib_nuclear_network(const fourdst::composition::Composition &composition); - LogicalReactionSet build_reaclib_nuclear_network(const fourdst::composition::Composition &composition, double culling, double T9 = 1.0); + reaction::REACLIBLogicalReactionSet build_reaclib_nuclear_network(const fourdst::composition::Composition &composition, bool reverse); + reaction::REACLIBLogicalReactionSet build_reaclib_nuclear_network_from_file(const std::string& filename, bool reverse); } // namespace nuclearNetwork diff --git a/src/network/include/gridfire/reaction/reaction.h b/src/network/include/gridfire/reaction/reaction.h new file mode 100644 index 00000000..250b32f9 --- /dev/null +++ b/src/network/include/gridfire/reaction/reaction.h @@ -0,0 +1,270 @@ +#pragma once + +#include + +#include "fourdst/composition/atomicSpecies.h" +#include "fourdst/logging/logging.h" +#include "quill/Logger.h" +#include +#include +#include +#include +#include + + +#include "cppad/cppad.hpp" + +namespace gridfire::reaction { + /** + * @struct REACLIBRateCoefficientSet + * @brief Holds the seven fitting parameters for a single REACLIB rate set. + * @details The thermonuclear reaction rate for a single set is calculated as: + * rate = exp(a0 + a1/T9 + a2/T9^(-1/3) + a3*T9^(1/3) + a4*T9 + a5*T9^(5/3) + a6*ln(T9)) + * where T9 is the temperature in billions of Kelvin. The total rate for a + * reaction is the sum of the rates from all its sets. + */ + + struct REACLIBRateCoefficientSet { + const double a0; + const double a1; + const double a2; + const double a3; + const double a4; + const double a5; + const double a6; + + friend std::ostream& operator<<(std::ostream& os, const REACLIBRateCoefficientSet& r) { + os << "[" << r.a0 << ", " << r.a1 << ", " << r.a2 << ", " + << r.a3 << ", " << r.a4 << ", " << r.a5 << ", " << r.a6 << "]"; + return os; + } + }; + + class Reaction { + public: + Reaction( + const std::string_view id, + const double qValue, + const std::vector& reactants, + const std::vector& products, + const bool reverse = false + ); + virtual ~Reaction() = default; + + virtual std::unique_ptr clone() const = 0; + + virtual double calculate_rate(double T9) const = 0; + virtual CppAD::AD calculate_rate(const CppAD::AD T9) const = 0; + + virtual std::string_view peName() const { return ""; } + + [[nodiscard]] bool contains(const fourdst::atomic::Species& species) const; + [[nodiscard]] bool contains_reactant(const fourdst::atomic::Species& species) const; + [[nodiscard]] bool contains_product(const fourdst::atomic::Species& species) const; + + std::unordered_set all_species() const; + std::unordered_set reactant_species() const; + std::unordered_set product_species() const; + + size_t num_species() const; + + int stoichiometry(const fourdst::atomic::Species& species) const; + std::unordered_map stoichiometry() const; + + std::string_view id() const { return m_id; } + double qValue() const { return m_qValue; } + const std::vector& reactants() const { return m_reactants; } + const std::vector& products() const { return m_products; } + bool is_reverse() const { return m_reverse; } + + double excess_energy() const; + + bool operator==(const Reaction& other) const { return m_id == other.m_id; } + bool operator!=(const Reaction& other) const { return !(*this == other); } + [[nodiscard]] uint64_t hash(uint64_t seed) const; + + protected: + quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); + std::string m_id; + double m_qValue = 0.0; ///< Q-value of the reaction + std::vector m_reactants; ///< Reactants of the reaction + std::vector m_products; ///< Products of the reaction + bool m_reverse = false; + }; + + class ReactionSet { + public: + explicit ReactionSet(std::vector> reactions); + ReactionSet(const ReactionSet& other); + ReactionSet& operator=(const ReactionSet& other); + virtual ~ReactionSet() = default; + + virtual void add_reaction(std::unique_ptr reaction); + virtual void remove_reaction(const std::unique_ptr& reaction); + + bool contains(const std::string_view& id) const; + bool contains(const Reaction& reaction) const; + + size_t size() const { return m_reactions.size(); } + + void sort(double T9=1.0); + + bool contains_species(const fourdst::atomic::Species& species) const; + bool contains_reactant(const fourdst::atomic::Species& species) const; + bool contains_product(const fourdst::atomic::Species& species) const; + + [[nodiscard]] const Reaction& operator[](size_t index) const; + [[nodiscard]] const Reaction& operator[](const std::string_view& id) const; + + bool operator==(const ReactionSet& other) const; + bool operator!=(const ReactionSet& other) const; + + [[nodiscard]] uint64_t hash(uint64_t seed = 0) const; + + auto begin() { return m_reactions.begin(); } + auto begin() const { return m_reactions.cbegin(); } + auto end() { return m_reactions.end(); } + auto end() const { return m_reactions.cend(); } + protected: + quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); + std::vector> m_reactions; + std::string m_id; + std::unordered_map m_reactionNameMap; ///< Maps reaction IDs to Reaction objects for quick lookup + + }; + + /** + * @struct REACLIBReaction + * @brief Represents a single nuclear reaction from the JINA REACLIB database. + * @details This struct is designed to be constructed at compile time (constexpr) from + * the data parsed by the Python generation script. It stores all necessary + * information to identify a reaction and calculate its rate. + */ + class REACLIBReaction final : public Reaction { + public: + /** + * @brief Constructs a REACLIBReaction object at compile time. + * @param id A unique string identifier generated by the Python script. + * @param peName + * @param chapter The REACLIB chapter number, defining the reaction structure. + * @param reactants A vector of strings with the names of the reactant species. + * @param products A vector of strings with the names of the product species. + * @param qValue The Q-value of the reaction in MeV. + * @param label The source label for the rate data (e.g., "wc12w", "st08"). + * @param sets A vector of RateFitSet, containing the fitting coefficients for the rate. + * @param reverse A boolean indicating if the reaction is reversed (default is false). + */ + REACLIBReaction( + const std::string_view id, + const std::string_view peName, + const int chapter, + const std::vector &reactants, + const std::vector &products, + const double qValue, + const std::string_view label, + const REACLIBRateCoefficientSet &sets, + const bool reverse = false); + + [[nodiscard]] std::unique_ptr clone() const override; + + [[nodiscard]] double calculate_rate(const double T9) const override; + [[nodiscard]] CppAD::AD calculate_rate(const CppAD::AD T9) const override; + + template + [[nodiscard]] GeneralScalarType calculate_rate(const GeneralScalarType T9) const { + const GeneralScalarType T913 = CppAD::pow(T9, 1.0/3.0); + const GeneralScalarType rateExponent = m_rateCoefficients.a0 + + m_rateCoefficients.a1 / T9 + + m_rateCoefficients.a2 / T913 + + m_rateCoefficients.a3 * T913 + + m_rateCoefficients.a4 * T9 + + m_rateCoefficients.a5 * CppAD::pow(T9, 5.0/3.0) + + m_rateCoefficients.a6 * CppAD::log(T9); + return CppAD::exp(rateExponent); + } + + [[nodiscard]] std::string_view peName() const override { return m_peName; } + + [[nodiscard]] int chapter() const { return m_chapter; } + + [[nodiscard]] std::string_view sourceLabel() const { return m_sourceLabel; } + + [[nodiscard]] const REACLIBRateCoefficientSet& rateCoefficients() const { return m_rateCoefficients; } + + friend std::ostream& operator<<(std::ostream& os, const REACLIBReaction& reaction); + private: + std::string m_peName; ///< Name of the reaction in (projectile, ejectile) notation (e.g. p(p, g)d) + int m_chapter; ///< Chapter number from the REACLIB database, defining the reaction structure. + std::string m_sourceLabel; ///< Source label for the rate data, indicating the origin of the rate coefficients (e.g., "wc12w", "st08"). + REACLIBRateCoefficientSet m_rateCoefficients; + }; + + class REACLIBReactionSet final : public ReactionSet { + public: + explicit REACLIBReactionSet(std::vector); + std::unordered_set peNames() const; + friend std::ostream& operator<<(std::ostream& os, const REACLIBReactionSet& set); + }; + + class REACLIBLogicalReaction final : public Reaction { + public: + explicit REACLIBLogicalReaction(const std::vector &reactions); + explicit REACLIBLogicalReaction(const REACLIBReaction &reaction); + void add_reaction(const REACLIBReaction& reaction); + + [[nodiscard]] std::unique_ptr clone() const override; + + [[nodiscard]] std::string_view peName() const override { return m_id; }; + [[nodiscard]] size_t size() const { return m_rates.size(); } + + [[nodiscard]] std::vector sources() const { return m_sources; } + + [[nodiscard]] double calculate_rate(const double T9) const override; + + [[nodiscard]] CppAD::AD calculate_rate(const CppAD::AD T9) const override; + + template + [[nodiscard]] T calculate_rate(const T T9) const { + T sum = static_cast(0.0); + const T T913 = CppAD::pow(T9, 1.0/3.0); + const T T953 = CppAD::pow(T9, 5.0/3.0); + const T logT9 = CppAD::log(T9); + // ReSharper disable once CppUseStructuredBinding + for (const auto& rate : m_rates) { + const T exponent = rate.a0 + + rate.a1 / T9 + + rate.a2 / T913 + + rate.a3 * T913 + + rate.a4 * T9 + + rate.a5 * T953 + + rate.a6 * logT9; + sum += CppAD::exp(exponent); + } + return sum; + } + + [[nodiscard]] int chapter() const { return m_chapter; } + + auto begin() { return m_rates.begin(); } + auto begin() const { return m_rates.cbegin(); } + auto end() { return m_rates.end(); } + auto end() const { return m_rates.cend(); } + + private: + int m_chapter; + std::vector m_sources; + std::vector m_rates; + + }; + + class REACLIBLogicalReactionSet final : public ReactionSet { + public: + REACLIBLogicalReactionSet() = delete; + explicit REACLIBLogicalReactionSet(const REACLIBReactionSet& reactionSet); + + [[nodiscard]] std::unordered_set peNames() const; + private: + std::unordered_set m_peNames; + }; + +} \ No newline at end of file diff --git a/src/network/include/gridfire/solver/solver.h b/src/network/include/gridfire/solver/solver.h new file mode 100644 index 00000000..50181e20 --- /dev/null +++ b/src/network/include/gridfire/solver/solver.h @@ -0,0 +1,256 @@ +#pragma once + +#include "gridfire/engine/engine_graph.h" +#include "gridfire/engine/engine_abstract.h" +#include "gridfire/network.h" + +#include "fourdst/logging/logging.h" +#include "fourdst/config/config.h" + +#include "quill/Logger.h" + +#include "Eigen/Dense" + +#include + +namespace gridfire::solver { + + struct dynamicQSESpeciesIndices { + std::vector dynamicSpeciesIndices; // Slow species that are not in QSE + std::vector QSESpeciesIndices; // Fast species that are in QSE + }; + + template + class NetworkSolverStrategy { + public: + explicit NetworkSolverStrategy(EngineT& engine) : m_engine(engine) {}; + virtual ~NetworkSolverStrategy() = default; + virtual NetOut evaluate(const NetIn& netIn) = 0; + protected: + EngineT& m_engine; + }; + + using DynamicNetworkSolverStrategy = NetworkSolverStrategy; + using StaticNetworkSolverStrategy = NetworkSolverStrategy; + + + class QSENetworkSolver final : public DynamicNetworkSolverStrategy { + public: + using DynamicNetworkSolverStrategy::DynamicNetworkSolverStrategy; + NetOut evaluate(const NetIn& netIn) override; + private: // methods + dynamicQSESpeciesIndices packSpeciesTypeIndexVectors( + const std::vector& Y, + const double T9, + const double rho + ) const; + Eigen::VectorXd calculateSteadyStateAbundances( + const std::vector& Y, + const double T9, + const double rho, + const dynamicQSESpeciesIndices& indices + ) const; + NetOut initializeNetworkWithShortIgnition( + const NetIn& netIn + ) const; + private: // Nested functors for ODE integration + struct RHSFunctor { + DynamicEngine& m_engine; + const std::vector& m_dynamicSpeciesIndices; + const std::vector& m_QSESpeciesIndices; + const Eigen::VectorXd& m_Y_QSE; + const double m_T9; + const double m_rho; + + RHSFunctor( + DynamicEngine& engine, + const std::vector& dynamicSpeciesIndices, + const std::vector& QSESpeciesIndices, + const Eigen::VectorXd& Y_QSE, + const double T9, + const double rho + ) : + m_engine(engine), + m_dynamicSpeciesIndices(dynamicSpeciesIndices), + m_QSESpeciesIndices(QSESpeciesIndices), + m_Y_QSE(Y_QSE), + m_T9(T9), + m_rho(rho) {} + + void operator()( + const boost::numeric::ublas::vector& YDynamic, + boost::numeric::ublas::vector& dYdtDynamic, + double t + ) const; + + }; + + struct JacobianFunctor { + DynamicEngine& m_engine; + const std::vector& m_dynamicSpeciesIndices; + const std::vector& m_QSESpeciesIndices; + const double m_T9; + const double m_rho; + + JacobianFunctor( + DynamicEngine& engine, + const std::vector& dynamicSpeciesIndices, + const std::vector& QSESpeciesIndices, + const double T9, + const double rho + ) : + m_engine(engine), + m_dynamicSpeciesIndices(dynamicSpeciesIndices), + m_QSESpeciesIndices(QSESpeciesIndices), + m_T9(T9), + m_rho(rho) {} + + void operator()( + const boost::numeric::ublas::vector& YDynamic, + boost::numeric::ublas::matrix& JDynamic, + double t, + boost::numeric::ublas::vector& dfdt + ) const; + }; + + template + struct EigenFunctor { + using InputType = Eigen::Matrix; + using OutputType = Eigen::Matrix; + using JacobianType = Eigen::Matrix; + + DynamicEngine& m_engine; + const std::vector& m_YDynamic; + const std::vector& m_dynamicSpeciesIndices; + const std::vector& m_QSESpeciesIndices; + const double m_T9; + const double m_rho; + + EigenFunctor( + DynamicEngine& engine, + const std::vector& YDynamic, + const std::vector& dynamicSpeciesIndices, + const std::vector& QSESpeciesIndices, + const double T9, + const double rho + ) : + m_engine(engine), + m_YDynamic(YDynamic), + m_dynamicSpeciesIndices(dynamicSpeciesIndices), + m_QSESpeciesIndices(QSESpeciesIndices), + m_T9(T9), + m_rho(rho) {} + + int operator()(const InputType& v_QSE, OutputType& f_QSE) const; + int df(const InputType& v_QSE, JacobianType& J_QSE) const; + }; + private: + quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); + fourdst::config::Config& m_config = fourdst::config::Config::getInstance(); + }; + + class DirectNetworkSolver final : public DynamicNetworkSolverStrategy { + public: + using DynamicNetworkSolverStrategy::DynamicNetworkSolverStrategy; + NetOut evaluate(const NetIn& netIn) override; + private: + struct RHSFunctor { + DynamicEngine& m_engine; + const double m_T9; + const double m_rho; + const size_t m_numSpecies; + + RHSFunctor( + DynamicEngine& engine, + const double T9, + const double rho + ) : + m_engine(engine), + m_T9(T9), + m_rho(rho), + m_numSpecies(engine.getNetworkSpecies().size()) {} + + void operator()( + const boost::numeric::ublas::vector& Y, + boost::numeric::ublas::vector& dYdt, + double t + ) const; + }; + struct JacobianFunctor { + DynamicEngine& m_engine; + const double m_T9; + const double m_rho; + const size_t m_numSpecies; + + JacobianFunctor( + DynamicEngine& engine, + const double T9, + const double rho + ) : + m_engine(engine), + m_T9(T9), + m_rho(rho), + m_numSpecies(engine.getNetworkSpecies().size()) {} + + void operator()( + const boost::numeric::ublas::vector& Y, + boost::numeric::ublas::matrix& J, + double t, + boost::numeric::ublas::vector& dfdt + ) const; + + }; + + private: + quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); + fourdst::config::Config& m_config = fourdst::config::Config::getInstance(); + }; + + template + int QSENetworkSolver::EigenFunctor::operator()(const InputType &v_QSE, OutputType &f_QSE) const { + std::vector YFull(m_engine.getNetworkSpecies().size(), 0.0); + Eigen::VectorXd Y_QSE(v_QSE.array().exp()); + for (size_t i = 0; i < m_dynamicSpeciesIndices.size(); ++i) { + YFull[m_dynamicSpeciesIndices[i]] = m_YDynamic[i]; + } + for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { + YFull[m_QSESpeciesIndices[i]] = Y_QSE(i); + } + const auto [full_dYdt, specificEnergyGenerationRate] = m_engine.calculateRHSAndEnergy(YFull, m_T9, m_rho); + f_QSE.resize(m_QSESpeciesIndices.size()); + for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { + f_QSE(i) = full_dYdt[m_QSESpeciesIndices[i]]; + } + + return 0; // Success + } + + template + int QSENetworkSolver::EigenFunctor::df(const InputType& v_QSE, JacobianType& J_QSE) const { + std::vector YFull(m_engine.getNetworkSpecies().size(), 0.0); + Eigen::VectorXd Y_QSE(v_QSE.array().exp()); + for (size_t i = 0; i < m_dynamicSpeciesIndices.size(); ++i) { + YFull[m_dynamicSpeciesIndices[i]] = m_YDynamic[i]; + } + for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { + YFull[m_QSESpeciesIndices[i]] = Y_QSE(i); + } + + m_engine.generateJacobianMatrix(YFull, m_T9, m_rho); + + Eigen::MatrixXd J_orig(m_QSESpeciesIndices.size(), m_QSESpeciesIndices.size()); + for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { + for (size_t j = 0; j < m_QSESpeciesIndices.size(); ++j) { + J_orig(i, j) = m_engine.getJacobianMatrixEntry(m_QSESpeciesIndices[i], m_QSESpeciesIndices[j]); + } + } + + J_QSE = J_orig; + for (long j = 0; j < J_QSE.cols(); ++j) { + J_QSE.col(j) *= Y_QSE(j); // Chain rule for log space + } + return 0; // Success + } + + +} \ No newline at end of file diff --git a/src/network/lib/approx8.cpp b/src/network/lib/engine/engine_approx8.cpp similarity index 95% rename from src/network/lib/approx8.cpp rename to src/network/lib/engine/engine_approx8.cpp index 2e5c4b23..a4b58169 100644 --- a/src/network/lib/approx8.cpp +++ b/src/network/lib/engine/engine_approx8.cpp @@ -28,7 +28,7 @@ #include "fourdst/config/config.h" #include "quill/LogMacros.h" -#include "gridfire/approx8.h" +#include "gridfire/engine/engine_approx8.h" #include "gridfire/network.h" /* Nuclear reaction network in cgs units based on Frank Timmes' "approx8". @@ -131,14 +131,14 @@ namespace gridfire::approx8{ return rate_fit(T9,a); } - // he3(he3,2p)he4 + // he3(he3,2p)he4 ** (missing both coefficients but have a reaction) 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 + // he4 + he4 + he4 -> c12 ** (missing middle coefficient but have other two) 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}; @@ -153,7 +153,7 @@ namespace gridfire::approx8{ return rate_fit(T9,a1) + rate_fit(T9,a2); } - // c12 + he4 -> o16 + // c12 + he4 -> o16 ** (missing first coefficient but have the second) 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}; @@ -508,27 +508,19 @@ namespace gridfire::approx8{ 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::ih1] = netIn.composition.getNumberFraction("H-1"); + std::cout << "Approx8::convert_netIn -> H-1 fraction: " << y[Approx8Net::ih1] << std::endl; + y[Approx8Net::ihe3] = netIn.composition.getNumberFraction("He-3"); + y[Approx8Net::ihe4] = netIn.composition.getNumberFraction("He-4"); + y[Approx8Net::ic12] = netIn.composition.getNumberFraction("C-12"); + y[Approx8Net::in14] = netIn.composition.getNumberFraction("N-14"); + y[Approx8Net::io16] = netIn.composition.getNumberFraction("O-16"); + y[Approx8Net::ine20] = netIn.composition.getNumberFraction("Ne-20"); + y[Approx8Net::img24] = netIn.composition.getNumberFraction("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; } }; diff --git a/src/network/lib/netgraph.cpp b/src/network/lib/engine/engine_graph.cpp similarity index 53% rename from src/network/lib/netgraph.cpp rename to src/network/lib/engine/engine_graph.cpp index 2deeb780..5d0bfafd 100644 --- a/src/network/lib/netgraph.cpp +++ b/src/network/lib/engine/engine_graph.cpp @@ -1,9 +1,9 @@ -#include "gridfire/netgraph.h" -#include "fourdst/composition/atomicSpecies.h" -#include "fourdst/constants/const.h" +#include "gridfire/engine/engine_graph.h" +#include "gridfire/reaction/reaction.h" #include "gridfire/network.h" -#include "gridfire/reaclib.h" + #include "fourdst/composition/species.h" +#include "fourdst/composition/atomicSpecies.h" #include "quill/LogMacros.h" @@ -14,6 +14,7 @@ #include #include #include +#include #include #include @@ -22,31 +23,28 @@ namespace gridfire { - GraphNetwork::GraphNetwork( + GraphEngine::GraphEngine( const fourdst::composition::Composition &composition ): - Network(REACLIB), - m_reactions(build_reaclib_nuclear_network(composition)) { + m_reactions(build_reaclib_nuclear_network(composition, false)) { syncInternalMaps(); } - GraphNetwork::GraphNetwork( - const fourdst::composition::Composition &composition, - const double cullingThreshold, - const double T9 - ): - Network(REACLIB), - m_reactions(build_reaclib_nuclear_network(composition, cullingThreshold, T9)) { - syncInternalMaps(); - } - - GraphNetwork::GraphNetwork(const reaclib::REACLIBReactionSet &reactions) : - Network(REACLIB), - m_reactions(reactions) { + GraphEngine::GraphEngine(reaction::REACLIBLogicalReactionSet reactions) : + m_reactions(std::move(reactions)) { syncInternalMaps(); } - void GraphNetwork::syncInternalMaps() { + StepDerivatives GraphEngine::calculateRHSAndEnergy( + const std::vector &Y, + const double T9, + const double rho + ) const { + return calculateAllDerivatives(Y, T9, rho); + } + + + void GraphEngine::syncInternalMaps() { collectNetworkSpecies(); populateReactionIDMap(); populateSpeciesToIndexMap(); @@ -56,17 +54,17 @@ namespace gridfire { } // --- Network Graph Construction Methods --- - void GraphNetwork::collectNetworkSpecies() { + void GraphEngine::collectNetworkSpecies() { m_networkSpecies.clear(); m_networkSpeciesMap.clear(); std::set uniqueSpeciesNames; for (const auto& reaction: m_reactions) { - for (const auto& reactant: reaction.reactants()) { + for (const auto& reactant: reaction->reactants()) { uniqueSpeciesNames.insert(reactant.name()); } - for (const auto& product: reaction.products()) { + for (const auto& product: reaction->products()) { uniqueSpeciesNames.insert(product.name()); } } @@ -84,84 +82,55 @@ namespace gridfire { } - void GraphNetwork::populateReactionIDMap() { - LOG_INFO(m_logger, "Populating reaction ID map for REACLIB graph network (serif::network::GraphNetwork)..."); + void GraphEngine::populateReactionIDMap() { + LOG_TRACE_L1(m_logger, "Populating reaction ID map for REACLIB graph network (serif::network::GraphNetwork)..."); m_reactionIDMap.clear(); for (const auto& reaction: m_reactions) { - m_reactionIDMap.insert({reaction.id(), reaction}); + m_reactionIDMap.emplace(reaction->id(), reaction.get()); } - LOG_INFO(m_logger, "Populated {} reactions in the reaction ID map.", m_reactionIDMap.size()); + LOG_TRACE_L1(m_logger, "Populated {} reactions in the reaction ID map.", m_reactionIDMap.size()); } - void GraphNetwork::populateSpeciesToIndexMap() { + void GraphEngine::populateSpeciesToIndexMap() { m_speciesToIndexMap.clear(); for (size_t i = 0; i < m_networkSpecies.size(); ++i) { m_speciesToIndexMap.insert({m_networkSpecies[i], i}); } } - void GraphNetwork::reserveJacobianMatrix() { + void GraphEngine::reserveJacobianMatrix() { // The implementation of this function (and others) constrains this nuclear network to a constant temperature and density during // each evaluation. size_t numSpecies = m_networkSpecies.size(); m_jacobianMatrix.clear(); m_jacobianMatrix.resize(numSpecies, numSpecies, false); // Sparse matrix, no initial values - LOG_INFO(m_logger, "Jacobian matrix resized to {} rows and {} columns.", + LOG_TRACE_L2(m_logger, "Jacobian matrix resized to {} rows and {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2()); } // --- Basic Accessors and Queries --- - const std::vector& GraphNetwork::getNetworkSpecies() const { + const std::vector& GraphEngine::getNetworkSpecies() const { // Returns a constant reference to the vector of unique species in the network. LOG_DEBUG(m_logger, "Providing access to network species vector. Size: {}.", m_networkSpecies.size()); return m_networkSpecies; } - const reaclib::REACLIBReactionSet& GraphNetwork::getNetworkReactions() const { + const reaction::REACLIBLogicalReactionSet& GraphEngine::getNetworkReactions() const { // Returns a constant reference to the set of reactions in the network. LOG_DEBUG(m_logger, "Providing access to network reactions set. Size: {}.", m_reactions.size()); return m_reactions; } - bool GraphNetwork::involvesSpecies(const fourdst::atomic::Species& species) const { + bool GraphEngine::involvesSpecies(const fourdst::atomic::Species& species) const { // Checks if a given species is present in the network's species map for efficient lookup. const bool found = m_networkSpeciesMap.contains(species.name()); LOG_DEBUG(m_logger, "Checking if species '{}' is involved in the network: {}.", species.name(), found ? "Yes" : "No"); return found; } - std::unordered_map GraphNetwork::getNetReactionStoichiometry(const reaclib::REACLIBReaction& reaction) const { - // Calculates the net stoichiometric coefficients for species in a given reaction. - std::unordered_map stoichiometry; - - // Iterate through reactants, decrementing their counts - for (const auto& reactant : reaction.reactants()) { - auto it = m_networkSpeciesMap.find(reactant.name()); - if (it != m_networkSpeciesMap.end()) { - stoichiometry[it->second]--; // Copy Species by value (PERF: Future performance improvements by using pointers or references (std::reference_wrapper) or something like that) - } else { - LOG_WARNING(m_logger, "Reactant species '{}' in reaction '{}' not found in network species map during stoichiometry calculation.", - reactant.name(), reaction.id()); - } - } - - // Iterate through products, incrementing their counts - for (const auto& product : reaction.products()) { - auto it = m_networkSpeciesMap.find(product.name()); - if (it != m_networkSpeciesMap.end()) { - stoichiometry[it->second]++; // Copy Species by value (PERF: Future performance improvements by using pointers or references (std::reference_wrapper) or something like that) - } else { - LOG_WARNING(m_logger, "Product species '{}' in reaction '{}' not found in network species map during stoichiometry calculation.", - product.name(), reaction.id()); - } - } - LOG_DEBUG(m_logger, "Calculated net stoichiometry for reaction '{}'. Total unique species in stoichiometry: {}.", reaction.id(), stoichiometry.size()); - return stoichiometry; - } - // --- Validation Methods --- - bool GraphNetwork::validateConservation() const { - LOG_INFO(m_logger, "Validating mass (A) and charge (Z) conservation across all reactions in the network."); + bool GraphEngine::validateConservation() const { + LOG_TRACE_L1(m_logger, "Validating mass (A) and charge (Z) conservation across all reactions in the network."); for (const auto& reaction : m_reactions) { uint64_t totalReactantA = 0; @@ -170,7 +139,7 @@ namespace gridfire { uint64_t totalProductZ = 0; // Calculate total A and Z for reactants - for (const auto& reactant : reaction.reactants()) { + for (const auto& reactant : reaction->reactants()) { auto it = m_networkSpeciesMap.find(reactant.name()); if (it != m_networkSpeciesMap.end()) { totalReactantA += it->second.a(); @@ -179,13 +148,13 @@ namespace gridfire { // This scenario indicates a severe data integrity issue: // a reactant is part of a reaction but not in the network's species map. LOG_ERROR(m_logger, "CRITICAL ERROR: Reactant species '{}' in reaction '{}' not found in network species map during conservation validation.", - reactant.name(), reaction.id()); + reactant.name(), reaction->id()); return false; } } // Calculate total A and Z for products - for (const auto& product : reaction.products()) { + for (const auto& product : reaction->products()) { auto it = m_networkSpeciesMap.find(product.name()); if (it != m_networkSpeciesMap.end()) { totalProductA += it->second.a(); @@ -193,7 +162,7 @@ namespace gridfire { } else { // Similar critical error for product species LOG_ERROR(m_logger, "CRITICAL ERROR: Product species '{}' in reaction '{}' not found in network species map during conservation validation.", - product.name(), reaction.id()); + product.name(), reaction->id()); return false; } } @@ -201,25 +170,24 @@ namespace gridfire { // Compare totals for conservation if (totalReactantA != totalProductA) { LOG_ERROR(m_logger, "Mass number (A) not conserved for reaction '{}': Reactants A={} vs Products A={}.", - reaction.id(), totalReactantA, totalProductA); + reaction->id(), totalReactantA, totalProductA); return false; } if (totalReactantZ != totalProductZ) { LOG_ERROR(m_logger, "Atomic number (Z) not conserved for reaction '{}': Reactants Z={} vs Products Z={}.", - reaction.id(), totalReactantZ, totalProductZ); + reaction->id(), totalReactantZ, totalProductZ); return false; } } - LOG_INFO(m_logger, "Mass (A) and charge (Z) conservation validated successfully for all reactions."); + LOG_TRACE_L1(m_logger, "Mass (A) and charge (Z) conservation validated successfully for all reactions."); return true; // All reactions passed the conservation check } - void GraphNetwork::validateComposition(const fourdst::composition::Composition &composition, double culling, double T9) { - + void GraphEngine::validateComposition(const fourdst::composition::Composition &composition, double culling, double T9) { // Check if the requested network has already been cached. // PERF: Rebuilding this should be pretty fast but it may be a good point of optimization in the future. - const reaclib::REACLIBReactionSet validationReactionSet = build_reaclib_nuclear_network(composition, culling, T9); + const reaction::REACLIBLogicalReactionSet validationReactionSet = build_reaclib_nuclear_network(composition, false); // TODO: need some more robust method here to // A. Build the basic network from the composition's species with non zero mass fractions. // B. rebuild a new composition from the reaction set's reactants + products (with the mass fractions from the things that are only products set to 0) @@ -230,22 +198,22 @@ namespace gridfire { // This allows for dynamic network modification while retaining caching for networks which are very similar. if (validationReactionSet != m_reactions) { - LOG_INFO(m_logger, "Reaction set not cached. Rebuilding the reaction set for T9={} and culling={}.", T9, culling); + LOG_DEBUG(m_logger, "Reaction set not cached. Rebuilding the reaction set for T9={} and culling={}.", T9, culling); m_reactions = validationReactionSet; syncInternalMaps(); // Re-sync internal maps after updating reactions. Note this will also retrace the AD tape. } } // --- Generate Stoichiometry Matrix --- - void GraphNetwork::generateStoichiometryMatrix() { - LOG_INFO(m_logger, "Generating stoichiometry matrix..."); + void GraphEngine::generateStoichiometryMatrix() { + LOG_TRACE_L1(m_logger, "Generating stoichiometry matrix..."); // Task 1: Set dimensions and initialize the matrix size_t numSpecies = m_networkSpecies.size(); size_t numReactions = m_reactions.size(); m_stoichiometryMatrix.resize(numSpecies, numReactions, false); - LOG_INFO(m_logger, "Stoichiometry matrix initialized with dimensions: {} rows (species) x {} columns (reactions).", + LOG_TRACE_L1(m_logger, "Stoichiometry matrix initialized with dimensions: {} rows (species) x {} columns (reactions).", numSpecies, numReactions); // Task 2: Populate the stoichiometry matrix @@ -253,13 +221,10 @@ namespace gridfire { size_t reactionColumnIndex = 0; for (const auto& reaction : m_reactions) { // Get the net stoichiometry for the current reaction - std::unordered_map netStoichiometry = getNetReactionStoichiometry(reaction); + std::unordered_map netStoichiometry = reaction->stoichiometry(); // Iterate through the species and their coefficients in the stoichiometry map - for (const auto& pair : netStoichiometry) { - const fourdst::atomic::Species& species = pair.first; // The Species object - const int coefficient = pair.second; // The stoichiometric coefficient - + for (const auto& [species, coefficient] : netStoichiometry) { // Find the row index for this species auto it = m_speciesToIndexMap.find(species); if (it != m_speciesToIndexMap.end()) { @@ -269,21 +234,49 @@ namespace gridfire { } else { // This scenario should ideally not happen if m_networkSpeciesMap and m_speciesToIndexMap are correctly synced LOG_ERROR(m_logger, "CRITICAL ERROR: Species '{}' from reaction '{}' stoichiometry not found in species to index map.", - species.name(), reaction.id()); + species.name(), reaction->id()); throw std::runtime_error("Species not found in species to index map: " + std::string(species.name())); } } reactionColumnIndex++; // Move to the next column for the next reaction } - LOG_INFO(m_logger, "Stoichiometry matrix population complete. Number of non-zero elements: {}.", + LOG_TRACE_L1(m_logger, "Stoichiometry matrix population complete. Number of non-zero elements: {}.", m_stoichiometryMatrix.nnz()); // Assuming nnz() exists for compressed_matrix } - void GraphNetwork::generateJacobianMatrix(const std::vector &Y, const double T9, - const double rho) { + StepDerivatives GraphEngine::calculateAllDerivatives( + const std::vector &Y_in, + const double T9, + const double rho + ) const { + return calculateAllDerivatives(Y_in, T9, rho); + } - LOG_INFO(m_logger, "Generating jacobian matrix for T9={}, rho={}..", T9, rho); + StepDerivatives GraphEngine::calculateAllDerivatives( + const std::vector &Y_in, + const ADDouble T9, + const ADDouble rho + ) const { + return calculateAllDerivatives(Y_in, T9, rho); + } + + double GraphEngine::calculateMolarReactionFlow( + const reaction::Reaction &reaction, + const std::vector &Y, + const double T9, + const double rho + ) const { + return calculateMolarReactionFlow(reaction, Y, T9, rho); + } + + void GraphEngine::generateJacobianMatrix( + const std::vector &Y, + const double T9, + const double rho + ) { + + LOG_TRACE_L1(m_logger, "Generating jacobian matrix for T9={}, rho={}..", T9, rho); const size_t numSpecies = m_networkSpecies.size(); // 1. Pack the input variables into a vector for CppAD @@ -307,51 +300,28 @@ namespace gridfire { } } } - LOG_INFO(m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2()); + LOG_DEBUG(m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2()); } - void GraphNetwork::detectStiff(const NetIn &netIn, const double T9, const unsigned long numSpecies, const boost::numeric::ublas::vector& Y) { - // --- Heuristic for automatic stiffness detection --- - const std::vector initial_y_stl(Y.begin(), Y.begin() + numSpecies); - const auto [dydt, specificEnergyRate] = calculateAllDerivatives(initial_y_stl, T9, netIn.density); - const std::vector& initial_dotY = dydt; - - double min_destruction_timescale = std::numeric_limits::max(); - - for (size_t i = 0; i < numSpecies; ++i) { - if (Y(i) > MIN_ABUNDANCE_THRESHOLD && initial_dotY[i] < 0.0) { - const double timescale = std::abs(Y(i) / initial_dotY[i]); - if (timescale < min_destruction_timescale) { - min_destruction_timescale = timescale; - } - } - } - - // If no species are being destroyed, the system is not stiff. - if (min_destruction_timescale == std::numeric_limits::max()) { - LOG_INFO(m_logger, "No species are undergoing net destruction. Network is considered non-stiff."); - m_stiff = false; - return; - } - - constexpr double saftey_factor = 10; - const bool is_stiff = (netIn.dt0 > saftey_factor * min_destruction_timescale); - - LOG_INFO(m_logger, "Fastest destruction timescale: {}. Initial dt0: {}. Stiffness detected: {}.", - min_destruction_timescale, netIn.dt0, is_stiff ? "Yes" : "No"); - - if (is_stiff) { - m_stiff = true; - LOG_INFO(m_logger, "Network is detected as stiff."); - } else { - m_stiff = false; - LOG_INFO(m_logger, "Network is detected as non-stiff."); - } - + double GraphEngine::getJacobianMatrixEntry(const int i, const int j) const { + return m_jacobianMatrix(i, j); } - void GraphNetwork::exportToDot(const std::string &filename) const { - LOG_INFO(m_logger, "Exporting network graph to DOT file: {}", filename); + std::unordered_map GraphEngine::getNetReactionStoichiometry( + const reaction::Reaction &reaction + ) const { + return reaction.stoichiometry(); + } + + int GraphEngine::getStoichiometryMatrixEntry( + const int speciesIndex, + const int reactionIndex + ) const { + return m_stoichiometryMatrix(speciesIndex, reactionIndex); + } + + void GraphEngine::exportToDot(const std::string &filename) const { + LOG_TRACE_L1(m_logger, "Exporting network graph to DOT file: {}", filename); std::ofstream dotFile(filename); if (!dotFile.is_open()) { @@ -375,118 +345,104 @@ namespace gridfire { dotFile << " // --- Reaction Edges ---\n"; for (const auto& reaction : m_reactions) { // Create a unique ID for the reaction node - std::string reactionNodeId = "reaction_" + std::string(reaction.id()); + std::string reactionNodeId = "reaction_" + std::string(reaction->id()); // Define the reaction node (small, black dot) dotFile << " \"" << reactionNodeId << "\" [shape=point, fillcolor=black, width=0.1, height=0.1, label=\"\"];\n"; // Draw edges from reactants to the reaction node - for (const auto& reactant : reaction.reactants()) { + for (const auto& reactant : reaction->reactants()) { dotFile << " \"" << reactant.name() << "\" -> \"" << reactionNodeId << "\";\n"; } // Draw edges from the reaction node to products - for (const auto& product : reaction.products()) { - dotFile << " \"" << reactionNodeId << "\" -> \"" << product.name() << "\" [label=\"" << reaction.qValue() << " MeV\"];\n"; + for (const auto& product : reaction->products()) { + dotFile << " \"" << reactionNodeId << "\" -> \"" << product.name() << "\" [label=\"" << reaction->qValue() << " MeV\"];\n"; } dotFile << "\n"; } dotFile << "}\n"; dotFile.close(); - LOG_INFO(m_logger, "Successfully exported network to {}", filename); + LOG_TRACE_L1(m_logger, "Successfully exported network to {}", filename); } - NetOut GraphNetwork::evaluate(const NetIn &netIn) { - namespace ublas = boost::numeric::ublas; - namespace odeint = boost::numeric::odeint; + void GraphEngine::exportToCSV(const std::string &filename) const { + LOG_TRACE_L1(m_logger, "Exporting network graph to CSV file: {}", filename); - const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9) - // validateComposition(netIn.composition, netIn.culling, T9); - - const unsigned long numSpecies = m_networkSpecies.size(); - constexpr double abs_tol = 1.0e-8; - constexpr double rel_tol = 1.0e-8; - - size_t stepCount = 0; - - // TODO: Pull these out into configuration options - - ODETerm rhs_functor(*this, T9, netIn.density); - - - ublas::vector Y(numSpecies + 1); - for (size_t i = 0; i < numSpecies; ++i) { - const auto& species = m_networkSpecies[i]; - // Get the mass fraction for this specific species from the input object - try { - Y(i) = netIn.composition.getMassFraction(std::string(species.name())); - } catch (const std::runtime_error &e) { - LOG_INFO(m_logger, "Species {} not in base composition, adding...", species.name()); - Y(i) = 0.0; // If the species is not in the composition, set its mass fraction to + std::ofstream csvFile(filename, std::ios::out | std::ios::trunc); + if (!csvFile.is_open()) { + LOG_ERROR(m_logger, "Failed to open file for writing: {}", filename); + throw std::runtime_error("Failed to open file for writing: " + filename); + } + csvFile << "Reaction;Reactants;Products;Q-value;sources;rates\n"; + for (const auto& reaction : m_reactions) { + // Dynamic cast to REACLIBReaction to access specific properties + csvFile << reaction->id() << ";"; + // Reactants + int count = 0; + for (const auto& reactant : reaction->reactants()) { + csvFile << reactant.name(); + if (++count < reaction->reactants().size()) { + csvFile << ","; + } } + csvFile << ";"; + count = 0; + for (const auto& product : reaction->products()) { + csvFile << product.name(); + if (++count < reaction->products().size()) { + csvFile << ","; + } + } + csvFile << ";" << reaction->qValue() << ";"; + // Reaction coefficients + auto* reaclibReaction = dynamic_cast(reaction.get()); + if (!reaclibReaction) { + LOG_ERROR(m_logger, "Failed to cast Reaction to REACLIBLogicalReaction in GraphNetwork::exportToCSV()."); + throw std::runtime_error("Failed to cast Reaction to REACLIBLogicalReaction in GraphNetwork::exportToCSV(). This should not happen, please check your reaction setup."); + } + auto sources = reaclibReaction->sources(); + count = 0; + for (const auto& source : sources) { + csvFile << source; + if (++count < sources.size()) { + csvFile << ","; + } + } + csvFile << ";"; + // Reaction coefficients + count = 0; + for (const auto& rates : *reaclibReaction) { + csvFile << rates; + if (++count < reaclibReaction->size()) { + csvFile << ","; + } + } + csvFile << "\n"; } - Y(numSpecies) = 0; // initial specific energy rate, will be updated later - - detectStiff(netIn, T9, numSpecies, Y); - m_stiff = false; - - if (m_stiff) { - JacobianTerm jacobian_functor(*this, T9, netIn.density); - LOG_INFO(m_logger, "Making use of stiff ODE solver for network evaluation."); - auto stepper = odeint::make_controlled>(abs_tol, rel_tol); - stepCount = odeint::integrate_adaptive( - stepper, - std::make_pair(rhs_functor, jacobian_functor), - Y, - 0.0, // Start time - netIn.tMax, - netIn.dt0 - ); - - } else { - LOG_INFO(m_logger, "Making use of ODE solver (non-stiff) for network evaluation."); - using state_type = ublas::vector; - auto stepper = odeint::make_controlled>(abs_tol, rel_tol); - stepCount = odeint::integrate_adaptive( - stepper, - rhs_functor, - Y, - 0.0, // Start time - netIn.tMax, - netIn.dt0 - ); - - - } - - double sumY = 0.0; - for (int i = 0; i < numSpecies; ++i) { sumY += Y(i); } - for (int i = 0; i < numSpecies; ++i) { Y(i) /= sumY; } - - // --- Marshall output variables --- - // PERF: Im sure this step could be tuned to avoid so many copies, that is a job for another day - std::vector speciesNames; - speciesNames.reserve(numSpecies); - for (const auto& species : m_networkSpecies) { - speciesNames.push_back(std::string(species.name())); - } - - std::vector finalAbundances(Y.begin(), Y.begin() + numSpecies); - fourdst::composition::Composition outputComposition(speciesNames, finalAbundances); - outputComposition.finalize(true); - - NetOut netOut; - netOut.composition = outputComposition; - netOut.num_steps = stepCount; - netOut.energy = Y(numSpecies); // The last element in Y is the specific energy rate - - return netOut; - + csvFile.close(); + LOG_TRACE_L1(m_logger, "Successfully exported network graph to {}", filename); } - void GraphNetwork::recordADTape() { - LOG_INFO(m_logger, "Recording AD tape for the RHS calculation..."); + std::unordered_map GraphEngine::getSpeciesTimescales(const std::vector &Y, const double T9, + const double rho) const { + auto [dydt, _] = calculateAllDerivatives(Y, T9, rho); + std::unordered_map speciesTimescales; + speciesTimescales.reserve(m_networkSpecies.size()); + for (size_t i = 0; i < m_networkSpecies.size(); ++i) { + double timescale = std::numeric_limits::infinity(); + const auto species = m_networkSpecies[i]; + if (std::abs(dydt[i]) > 0.0) { + timescale = std::abs(Y[i] / dydt[i]); + } + speciesTimescales.emplace(species, timescale); + } + return speciesTimescales; + } + + void GraphEngine::recordADTape() { + LOG_TRACE_L1(m_logger, "Recording AD tape for the RHS calculation..."); // Task 1: Set dimensions and initialize the matrix const size_t numSpecies = m_networkSpecies.size(); @@ -521,11 +477,11 @@ namespace gridfire { // 5. Call the actual templated function // We let T9 and rho be constant, so we pass them as fixed values. - auto derivatives = calculateAllDerivatives>(adY, adT9, adRho); + auto [dydt, nuclearEnergyGenerationRate] = calculateAllDerivatives>(adY, adT9, adRho); - m_rhsADFun.Dependent(adInput, derivatives.dydt); + m_rhsADFun.Dependent(adInput, dydt); - LOG_INFO(m_logger, "AD tape recorded successfully for the RHS calculation. Number of independent variables: {}.", + LOG_TRACE_L1(m_logger, "AD tape recorded successfully for the RHS calculation. Number of independent variables: {}.", adInput.size()); } } diff --git a/src/network/lib/network.cpp b/src/network/lib/network.cpp index fc3bafa6..6793cf90 100644 --- a/src/network/lib/network.cpp +++ b/src/network/lib/network.cpp @@ -19,13 +19,13 @@ // // *********************************************************************** */ #include "gridfire/network.h" +#include "gridfire/reactions.h" +#include "../include/gridfire/reaction/reaction.h" #include +#include -#include "gridfire/approx8.h" #include "quill/LogMacros.h" -#include "gridfire/reaclib.h" -#include "gridfire/reactions.h" namespace gridfire { Network::Network(const NetworkFormat format) : @@ -50,231 +50,20 @@ namespace gridfire { 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; - } - - LogicalReaction::LogicalReaction(const std::vector &reactions) { - if (reactions.empty()) { - LOG_ERROR(m_logger, "Empty reaction set provided to LogicalReaction constructor."); - throw std::runtime_error("Empty reaction set provided to LogicalReaction constructor."); - } - - const auto& first_reaction = reactions.front(); - m_peID = first_reaction.peName(); - m_reactants = first_reaction.reactants(); - m_products = first_reaction.products(); - m_qValue = first_reaction.qValue(); - m_reverse = first_reaction.is_reverse(); - m_chapter = first_reaction.chapter(); - m_rates.reserve(reactions.size()); - for (const auto& reaction : reactions) { - m_rates.push_back(reaction.rateFits()); - if (std::abs(reaction.qValue() - m_qValue) > 1e-6) { - LOG_ERROR(m_logger, "Inconsistent Q-values in reactions {}. All reactions must have the same Q-value.", m_peID); - throw std::runtime_error("Inconsistent Q-values in reactions. All reactions must have the same Q-value."); - } - } - } - - LogicalReaction::LogicalReaction(const reaclib::REACLIBReaction &first_reaction) { - m_peID = first_reaction.peName(); - m_reactants = first_reaction.reactants(); - m_products = first_reaction.products(); - m_qValue = first_reaction.qValue(); - m_reverse = first_reaction.is_reverse(); - m_chapter = first_reaction.chapter(); - m_rates.reserve(1); - m_rates.push_back(first_reaction.rateFits()); - - } - - void LogicalReaction::add_reaction(const reaclib::REACLIBReaction &reaction) { - if (std::abs(reaction.qValue() - m_qValue > 1e-6)) { - LOG_ERROR(m_logger, "Inconsistent Q-values in reactions {}. All reactions must have the same Q-value.", m_peID); - throw std::runtime_error("Inconsistent Q-values in reactions. All reactions must have the same Q-value."); - } - m_rates.push_back(reaction.rateFits()); - } - - auto LogicalReaction::begin() { - return m_rates.begin(); - } - auto LogicalReaction::begin() const { - return m_rates.cbegin(); - } - auto LogicalReaction::end() { - return m_rates.end(); - } - auto LogicalReaction::end() const { - return m_rates.cend(); - } - - LogicalReactionSet::LogicalReactionSet(const std::vector &reactions) { - if (reactions.empty()) { - LOG_ERROR(m_logger, "Empty reaction set provided to LogicalReactionSet constructor."); - throw std::runtime_error("Empty reaction set provided to LogicalReactionSet constructor."); - } - - for (const auto& reaction : reactions) { - m_reactions.push_back(reaction); - } - - m_reactionNameMap.reserve(m_reactions.size()); - for (const auto& reaction : m_reactions) { - if (m_reactionNameMap.contains(reaction.id())) { - LOG_ERROR(m_logger, "Duplicate reaction ID '{}' found in LogicalReactionSet.", reaction.id()); - throw std::runtime_error("Duplicate reaction ID found in LogicalReactionSet: " + std::string(reaction.id())); - } - m_reactionNameMap[reaction.id()] = reaction; - } - } - - LogicalReactionSet::LogicalReactionSet(const std::vector &reactions) { - std::vector uniquePeNames; - for (const auto& reaction: reactions) { - if (uniquePeNames.empty()) { - uniquePeNames.emplace_back(reaction.peName()); - } else if (std::ranges::find(uniquePeNames, reaction.peName()) == uniquePeNames.end()) { - uniquePeNames.emplace_back(reaction.peName()); - } - } - - for (const auto& peName : uniquePeNames) { - std::vector reactionsForPe; - for (const auto& reaction : reactions) { - if (reaction.peName() == peName) { - reactionsForPe.push_back(reaction); - } - } - m_reactions.emplace_back(reactionsForPe); - } - - m_reactionNameMap.reserve(m_reactions.size()); - for (const auto& reaction : m_reactions) { - if (m_reactionNameMap.contains(reaction.id())) { - LOG_ERROR(m_logger, "Duplicate reaction ID '{}' found in LogicalReactionSet.", reaction.id()); - throw std::runtime_error("Duplicate reaction ID found in LogicalReactionSet: " + std::string(reaction.id())); - } - m_reactionNameMap[reaction.id()] = reaction; - } - } - - LogicalReactionSet::LogicalReactionSet(const reaclib::REACLIBReactionSet &reactionSet) { - LogicalReactionSet(reactionSet.get_reactions()); - } - - void LogicalReactionSet::add_reaction(const LogicalReaction& reaction) { - if (m_reactionNameMap.contains(reaction.id())) { - LOG_WARNING(m_logger, "Reaction {} already exists in the set. Not adding again.", reaction.id()); - std::cerr << "Warning: Reaction " << reaction.id() << " already exists in the set. Not adding again." << std::endl; - return; - } - m_reactions.push_back(reaction); - m_reactionNameMap[reaction.id()] = reaction; - } - - void LogicalReactionSet::add_reaction(const reaclib::REACLIBReaction& reaction) { - if (m_reactionNameMap.contains(reaction.id())) { - m_reactionNameMap[reaction.id()].add_reaction(reaction); - } else { - const LogicalReaction logicalReaction(reaction); - m_reactions.push_back(logicalReaction); - m_reactionNameMap[reaction.id()] = logicalReaction; - } - } - - bool LogicalReactionSet::contains(const std::string_view &id) const { - for (const auto& reaction : m_reactions) { - if (reaction.id() == id) { - return true; - } - } - return false; - } - - bool LogicalReactionSet::contains(const LogicalReaction &reactions) const { - for (const auto& reaction : m_reactions) { - if (reaction.id() == reactions.id()) { - return true; - } - } - return false; - } - - bool LogicalReactionSet::contains_species(const fourdst::atomic::Species &species) const { - return contains_reactant(species) || contains_product(species); - } - - bool LogicalReactionSet::contains_reactant(const fourdst::atomic::Species &species) const { - for (const auto& reaction : m_reactions) { - if (std::ranges::find(reaction.reactants(), species) != reaction.reactants().end()) { - return true; - } - } - return false; - } - - bool LogicalReactionSet::contains_product(const fourdst::atomic::Species &species) const { - for (const auto& reaction : m_reactions) { - if (std::ranges::find(reaction.products(), species) != reaction.products().end()) { - return true; - } - } - return false; - } - - const LogicalReaction & LogicalReactionSet::operator[](size_t index) const { - return m_reactions.at(index); - } - - const LogicalReaction & LogicalReactionSet::operator[](const std::string_view &id) const { - return m_reactionNameMap.at(id); - } - - auto LogicalReactionSet::begin() { - return m_reactions.begin(); - } - - auto LogicalReactionSet::begin() const { - return m_reactions.cbegin(); - } - - auto LogicalReactionSet::end() { - return m_reactions.end(); - } - - auto LogicalReactionSet::end() const { - return m_reactions.cend(); - } - - LogicalReactionSet build_reaclib_nuclear_network(const fourdst::composition::Composition &composition) { - using namespace reaclib; - REACLIBReactionSet reactions; + reaction::REACLIBLogicalReactionSet build_reaclib_nuclear_network(const fourdst::composition::Composition &composition, bool reverse) { + using namespace reaction; + std::vector reaclibReactions; auto logger = fourdst::logging::LogManager::getInstance().getLogger("log"); - if (!s_initialized) { - LOG_INFO(logger, "REACLIB reactions not initialized. Calling initializeAllReaclibReactions()..."); - initializeAllReaclibReactions(); + if (!reaclib::s_initialized) { + LOG_DEBUG(logger, "REACLIB reactions not initialized. Calling initializeAllReaclibReactions()..."); + reaclib::initializeAllReaclibReactions(); } - for (const auto &reaction: s_all_reaclib_reactions | std::views::values) { + for (const auto &reaction: reaclib::s_all_reaclib_reactions | std::views::values) { + if (reaction.is_reverse() != reverse) { + continue; // Skip reactions that do not match the requested direction + } bool gotReaction = true; const auto& reactants = reaction.reactants(); for (const auto& reactant : reactants) { @@ -284,24 +73,85 @@ namespace gridfire { } } if (gotReaction) { - LOG_INFO(logger, "Adding reaction {} to REACLIB reaction set.", reaction.peName()); - reactions.add_reaction(reaction); + LOG_TRACE_L3(logger, "Adding reaction {} to REACLIB reaction set.", reaction.peName()); + reaclibReactions.push_back(reaction); } } - reactions.sort(); - return LogicalReactionSet(reactions); + const REACLIBReactionSet reactionSet(reaclibReactions); + return REACLIBLogicalReactionSet(reactionSet); } - LogicalReactionSet build_reaclib_nuclear_network(const fourdst::composition::Composition &composition, const double culling, const double T9) { - using namespace reaclib; - LogicalReactionSet allReactions = build_reaclib_nuclear_network(composition); - LogicalReactionSet reactions; - for (const auto& reaction : allReactions) { - if (reaction.calculate_rate(T9) >= culling) { - reactions.add_reaction(reaction); + // Trim whitespace from both ends of a string + std::string trim_whitespace(const std::string& str) { + auto startIt = str.begin(); + auto endIt = str.end(); + + while (startIt != endIt && std::isspace(static_cast(*startIt))) { + ++startIt; + } + if (startIt == endIt) { + return ""; + } + auto ritr = std::find_if(str.rbegin(), std::string::const_reverse_iterator(startIt), + [](unsigned char ch){ return !std::isspace(ch); }); + return std::string(startIt, ritr.base()); + } + + reaction::REACLIBLogicalReactionSet build_reaclib_nuclear_network_from_file(const std::string& filename, bool reverse) { + const auto logger = fourdst::logging::LogManager::getInstance().getLogger("log"); + std::vector reactionPENames; + std::ifstream infile(filename, std::ios::in); + if (!infile.is_open()) { + LOG_ERROR(logger, "Failed to open network definition file: {}", filename); + throw std::runtime_error("Failed to open network definition file: " + filename); + } + std::string line; + while (std::getline(infile, line)) { + std::string trimmedLine = trim_whitespace(line); + if (trimmedLine.empty() || trimmedLine.starts_with('#')) { + continue; // Skip empty lines and comments + } + reactionPENames.push_back(trimmedLine); + } + infile.close(); + std::vector reaclibReactions; + + if (reactionPENames.empty()) { + LOG_ERROR(logger, "No reactions found in the network definition file: {}", filename); + throw std::runtime_error("No reactions found in the network definition file: " + filename); + } + + if (!reaclib::s_initialized) { + LOG_DEBUG(logger, "REACLIB reactions not initialized. Calling initializeAllReaclibReactions()..."); + reaclib::initializeAllReaclibReactions(); + } else { + LOG_DEBUG(logger, "REACLIB reactions already initialized."); + } + + for (const auto& peName : reactionPENames) { + bool found = false; + for (const auto& reaction : reaclib::s_all_reaclib_reactions | std::views::values) { + if (reaction.peName() == peName && reaction.is_reverse() == reverse) { + reaclibReactions.push_back(reaction); + found = true; + LOG_TRACE_L3(logger, "Found reaction {} (version {}) in REACLIB database.", peName, reaction.sourceLabel()); + } + } + if (!found) { + LOG_ERROR(logger, "Reaction {} not found in REACLIB database. Skipping...", peName); + throw std::runtime_error("Reaction not found in REACLIB database."); } } - return reactions; + + + // Log any reactions that were not found in the REACLIB database + for (const auto& peName : reactionPENames) { + if (std::ranges::find(reaclibReactions, peName, &reaction::REACLIBReaction::peName) == reaclibReactions.end()) { + LOG_WARNING(logger, "Reaction {} not found in REACLIB database.", peName); + } + } + const reaction::REACLIBReactionSet reactionSet(reaclibReactions); + return reaction::REACLIBLogicalReactionSet(reactionSet); } } diff --git a/src/network/lib/reaction/reaction.cpp b/src/network/lib/reaction/reaction.cpp new file mode 100644 index 00000000..ce68e9bf --- /dev/null +++ b/src/network/lib/reaction/reaction.cpp @@ -0,0 +1,427 @@ +#include "gridfire/reaction/reaction.h" + +#include +#include +#include +#include +#include +#include + +#include "quill/LogMacros.h" + +#include "fourdst/composition/atomicSpecies.h" + +#include "xxhash64.h" + +namespace gridfire::reaction { + using namespace fourdst::atomic; + + Reaction::Reaction( + const std::string_view id, + const double qValue, + const std::vector& reactants, + const std::vector& products, + const bool reverse) : + m_id(id), + m_qValue(qValue), + m_reactants(std::move(reactants)), + m_products(std::move(products)), + m_reverse(reverse) {} + + bool Reaction::contains(const Species &species) const { + return contains_reactant(species) || contains_product(species); + } + + + bool Reaction::contains_reactant(const Species& species) const { + for (const auto& reactant : m_reactants) { + if (reactant == species) { + return true; + } + } + return false; + } + + bool Reaction::contains_product(const Species& species) const { + for (const auto& product : m_products) { + if (product == species) { + return true; + } + } + return false; + } + + std::unordered_set Reaction::all_species() const { + auto rs = reactant_species(); + auto ps = product_species(); + rs.insert(ps.begin(), ps.end()); + return rs; + } + + std::unordered_set Reaction::reactant_species() const { + std::unordered_set reactantsSet; + for (const auto& reactant : m_reactants) { + reactantsSet.insert(reactant); + } + return reactantsSet; + } + + std::unordered_set Reaction::product_species() const { + std::unordered_set productsSet; + for (const auto& product : m_products) { + productsSet.insert(product); + } + return productsSet; + } + + int Reaction::stoichiometry(const Species& species) const { + int s = 0; + for (const auto& reactant : m_reactants) { + if (reactant == species) { + s--; + } + } + for (const auto& product : m_products) { + if (product == species) { + s++; + } + } + return s; + } + + size_t Reaction::num_species() const { + return all_species().size(); + } + + std::unordered_map Reaction::stoichiometry() const { + std::unordered_map stoichiometryMap; + for (const auto& reactant : m_reactants) { + stoichiometryMap[reactant]--; + } + for (const auto& product : m_products) { + stoichiometryMap[product]++; + } + return stoichiometryMap; + } + + double Reaction::excess_energy() const { + double reactantMass = 0.0; + double productMass = 0.0; + constexpr double AMU2MeV = 931.494893; // Conversion factor from atomic mass unit to MeV + for (const auto& reactant : m_reactants) { + reactantMass += reactant.mass(); + } + for (const auto& product : m_products) { + productMass += product.mass(); + } + return (reactantMass - productMass) * AMU2MeV; + } + + uint64_t Reaction::hash(uint64_t seed) const { + return XXHash64::hash(m_id.data(), m_id.size(), seed); + } + + ReactionSet::ReactionSet( + std::vector> reactions) : + m_reactions(std::move(reactions)) { + if (m_reactions.empty()) { + return; // Case where the reactions will be added later. + } + m_reactionNameMap.reserve(reactions.size()); + for (const auto& reaction : m_reactions) { + m_id += reaction->id(); + m_reactionNameMap.emplace(reaction->id(), reaction.get()); + } + } + + ReactionSet::ReactionSet(const ReactionSet &other) { + m_reactions.reserve(other.m_reactions.size()); + for (const auto& reaction_ptr: other.m_reactions) { + m_reactions.push_back(reaction_ptr->clone()); + } + + m_reactionNameMap.reserve(other.m_reactionNameMap.size()); + for (const auto& reaction_ptr : m_reactions) { + m_reactionNameMap.emplace(reaction_ptr->id(), reaction_ptr.get()); + } + } + + ReactionSet & ReactionSet::operator=(const ReactionSet &other) { + if (this != &other) { + ReactionSet temp(other); + std::swap(m_reactions, temp.m_reactions); + std::swap(m_reactionNameMap, temp.m_reactionNameMap); + } + return *this; + } + + void ReactionSet::add_reaction(std::unique_ptr reaction) { + m_reactions.emplace_back(std::move(reaction)); + m_id += m_reactions.back()->id(); + m_reactionNameMap.emplace(m_reactions.back()->id(), m_reactions.back().get()); + } + + void ReactionSet::remove_reaction(const std::unique_ptr& reaction) { + if (!m_reactionNameMap.contains(std::string(reaction->id()))) { + // LOG_INFO(m_logger, "Attempted to remove reaction {} that does not exist in ReactionSet. Skipping...", reaction->id()); + return; + } + + m_reactionNameMap.erase(std::string(reaction->id())); + + std::erase_if(m_reactions, [&reaction](const std::unique_ptr& r) { + return *r == *reaction; + }); + } + + bool ReactionSet::contains(const std::string_view& id) const { + for (const auto& reaction : m_reactions) { + if (reaction->id() == id) { + return true; + } + } + return false; + } + + bool ReactionSet::contains(const Reaction& reaction) const { + for (const auto& r : m_reactions) { + if (*r == reaction) { + return true; + } + } + return false; + } + + void ReactionSet::sort(double T9) { + std::ranges::sort(m_reactions, + [&T9](const std::unique_ptr& r1, const std::unique_ptr& r2) { + return r1->calculate_rate(T9) < r2->calculate_rate(T9); + }); + } + + bool ReactionSet::contains_species(const Species& species) const { + for (const auto& reaction : m_reactions) { + if (reaction->contains(species)) { + return true; + } + } + return false; + } + + bool ReactionSet::contains_reactant(const Species& species) const { + for (const auto& r : m_reactions) { + if (r->contains_reactant(species)) { + return true; + } + } + return false; + } + + bool ReactionSet::contains_product(const Species& species) const { + for (const auto& r : m_reactions) { + if (r->contains_product(species)) { + return true; + } + } + return false; + } + + const Reaction& ReactionSet::operator[](const size_t index) const { + if (index >= m_reactions.size()) { + throw std::out_of_range("Index" + std::to_string(index) + " out of range for ReactionSet of size " + std::to_string(m_reactions.size()) + "."); + } + return *m_reactions[index]; + } + + const Reaction& ReactionSet::operator[](const std::string_view& id) const { + if (auto it = m_reactionNameMap.find(std::string(id)); it != m_reactionNameMap.end()) { + return *it->second; + } + throw std::out_of_range("Species " + std::string(id) + " does not exist in ReactionSet."); + } + + bool ReactionSet::operator==(const ReactionSet& other) const { + if (size() != other.size()) { + return false; + } + return hash() == other.hash(); + } + + bool ReactionSet::operator!=(const ReactionSet& other) const { + return !(*this == other); + } + + uint64_t ReactionSet::hash(uint64_t seed) const { + if (m_reactions.empty()) { + return XXHash64::hash(nullptr, 0, seed); + } + std::vector individualReactionHashes; + individualReactionHashes.reserve(m_reactions.size()); + for (const auto& reaction : m_reactions) { + individualReactionHashes.push_back(reaction->hash(seed)); + } + + std::ranges::sort(individualReactionHashes); + + const void* data = static_cast(individualReactionHashes.data()); + size_t sizeInBytes = individualReactionHashes.size() * sizeof(uint64_t); + return XXHash64::hash(data, sizeInBytes, seed); + } + + REACLIBReaction::REACLIBReaction( + const std::string_view id, + const std::string_view peName, + const int chapter, + const std::vector &reactants, + const std::vector &products, + const double qValue, + const std::string_view label, + const REACLIBRateCoefficientSet &sets, + const bool reverse) : + Reaction(id, qValue, reactants, products, reverse), + m_peName(peName), + m_chapter(chapter), + m_sourceLabel(label), + m_rateCoefficients(sets) {} + + std::unique_ptr REACLIBReaction::clone() const { + return std::make_unique(*this); + } + + + double REACLIBReaction::calculate_rate(const double T9) const { + return calculate_rate(T9); + } + + CppAD::AD REACLIBReaction::calculate_rate(const CppAD::AD T9) const { + return calculate_rate>(T9); + } + + REACLIBReactionSet::REACLIBReactionSet(std::vector reactions) : + ReactionSet(std::vector>()) { + // Convert REACLIBReaction to unique_ptr and store in m_reactions + m_reactions.reserve(reactions.size()); + m_reactionNameMap.reserve(reactions.size()); + for (auto& reaction : reactions) { + m_reactions.emplace_back(std::make_unique(std::move(reaction))); + m_reactionNameMap.emplace(std::string(reaction.id()), m_reactions.back().get()); + } + } + + std::unordered_set REACLIBReactionSet::peNames() const { + std::unordered_set peNames; + for (const auto& reactionPtr: m_reactions) { + if (const auto* reaction = dynamic_cast(reactionPtr.get())) { + peNames.insert(std::string(reaction->peName())); + } else { + // LOG_ERROR(m_logger, "Failed to cast Reaction to REACLIBReaction in REACLIBReactionSet::peNames()."); + throw std::runtime_error("Failed to cast Reaction to REACLIBReaction in REACLIBReactionSet::peNames(). This should not happen, please check your reaction setup."); + } + } + return peNames; + } + + REACLIBLogicalReaction::REACLIBLogicalReaction(const std::vector& reactants) : + Reaction(reactants.front().peName(), + reactants.front().qValue(), + reactants.front().reactants(), + reactants.front().products(), + reactants.front().is_reverse()), + m_chapter(reactants.front().chapter()) { + + m_sources.reserve(reactants.size()); + m_rates.reserve(reactants.size()); + for (const auto& reaction : reactants) { + if (std::abs(reaction.qValue() - m_qValue) > 1e-6) { + LOG_ERROR(m_logger, "REACLIBLogicalReaction constructed with reactions having different Q-values. Expected {} got {}.", m_qValue, reaction.qValue()); + throw std::runtime_error("REACLIBLogicalReaction constructed with reactions having different Q-values. Expected " + std::to_string(m_qValue) + " got " + std::to_string(reaction.qValue()) + "."); + } + m_sources.push_back(std::string(reaction.sourceLabel())); + m_rates.push_back(reaction.rateCoefficients()); + } + } + + REACLIBLogicalReaction::REACLIBLogicalReaction(const REACLIBReaction& reaction) : + Reaction(reaction.peName(), + reaction.qValue(), + reaction.reactants(), + reaction.products(), + reaction.is_reverse()), + m_chapter(reaction.chapter()) { + m_sources.reserve(1); + m_rates.reserve(1); + m_sources.push_back(std::string(reaction.sourceLabel())); + m_rates.push_back(reaction.rateCoefficients()); + } + + void REACLIBLogicalReaction::add_reaction(const REACLIBReaction& reaction) { + if (reaction.peName() != m_id) { + LOG_ERROR(m_logger, "Cannot add reaction with different peName to REACLIBLogicalReaction. Expected {} got {}.", m_id, reaction.peName()); + throw std::runtime_error("Cannot add reaction with different peName to REACLIBLogicalReaction. Expected " + std::string(m_id) + " got " + std::string(reaction.peName()) + "."); + } + for (const auto& source : m_sources) { + if (source == reaction.sourceLabel()) { + LOG_ERROR(m_logger, "Cannot add reaction with duplicate source label {} to REACLIBLogicalReaction.", reaction.sourceLabel()); + throw std::runtime_error("Cannot add reaction with duplicate source label " + std::string(reaction.sourceLabel()) + " to REACLIBLogicalReaction."); + } + } + if (std::abs(reaction.qValue() - m_qValue) > 1e-6) { + LOG_ERROR(m_logger, "REACLIBLogicalReaction constructed with reactions having different Q-values. Expected {} got {}.", m_qValue, reaction.qValue()); + throw std::runtime_error("REACLIBLogicalReaction constructed with reactions having different Q-values. Expected " + std::to_string(m_qValue) + " got " + std::to_string(reaction.qValue()) + "."); + } + m_sources.push_back(std::string(reaction.sourceLabel())); + m_rates.push_back(reaction.rateCoefficients()); + } + + std::unique_ptr REACLIBLogicalReaction::clone() const { + return std::make_unique(*this); + } + + double REACLIBLogicalReaction::calculate_rate(const double T9) const { + return calculate_rate(T9); + } + + CppAD::AD REACLIBLogicalReaction::calculate_rate(const CppAD::AD T9) const { + return calculate_rate>(T9); + } + + REACLIBLogicalReactionSet::REACLIBLogicalReactionSet(const REACLIBReactionSet &reactionSet) : + ReactionSet(std::vector>()) { + + std::unordered_map> grouped_reactions; + + for (const auto& reaction_ptr : reactionSet) { + if (const auto* reaclib_reaction = dynamic_cast(reaction_ptr.get())) { + grouped_reactions[reaclib_reaction->peName()].push_back(*reaclib_reaction); + } + } + m_reactions.reserve(grouped_reactions.size()); + m_reactionNameMap.reserve(grouped_reactions.size()); + for (const auto& [peName, reactions_for_peName] : grouped_reactions) { + m_peNames.insert(std::string(peName)); + auto logical_reaction = std::make_unique(reactions_for_peName); + m_reactionNameMap.emplace(logical_reaction->id(), logical_reaction.get()); + m_reactions.push_back(std::move(logical_reaction)); + } + } + + std::unordered_set REACLIBLogicalReactionSet::peNames() const { + return m_peNames; + } +} + +namespace std { + template<> + struct hash { + size_t operator()(const gridfire::reaction::Reaction& r) const noexcept { + return r.hash(0); + } + }; + + template<> + struct hash { + size_t operator()(const gridfire::reaction::ReactionSet& s) const noexcept { + return s.hash(0); + } + }; +} // namespace std diff --git a/src/network/lib/solver/solver.cpp b/src/network/lib/solver/solver.cpp new file mode 100644 index 00000000..620787de --- /dev/null +++ b/src/network/lib/solver/solver.cpp @@ -0,0 +1,384 @@ +#include "gridfire/solver/solver.h" +#include "gridfire/engine/engine_graph.h" +#include "gridfire/network.h" + +#include "fourdst/composition/atomicSpecies.h" +#include "fourdst/composition/composition.h" +#include "fourdst/config/config.h" + +#include "Eigen/Dense" +#include "unsupported/Eigen/NonLinearOptimization" + +#include + +#include +#include +#include +#include + +#include "quill/LogMacros.h" + +namespace gridfire::solver { + + NetOut QSENetworkSolver::evaluate(const NetIn &netIn) { + using state_type = boost::numeric::ublas::vector; + using namespace boost::numeric::odeint; + + NetOut postIgnition = initializeNetworkWithShortIgnition(netIn); + + constexpr double abundance_floor = 1.0e-30; + std::vectorY_sanitized_initial; + Y_sanitized_initial.reserve(m_engine.getNetworkSpecies().size()); + + LOG_DEBUG(m_logger, "Sanitizing initial abundances with a floor of {:0.3E}...", abundance_floor); + for (const auto& species : m_engine.getNetworkSpecies()) { + double molar_abundance = 0.0; + if (postIgnition.composition.contains(species)) { + molar_abundance = postIgnition.composition.getMolarAbundance(std::string(species.name())); + } + + if (molar_abundance < abundance_floor) { + molar_abundance = abundance_floor; + } + Y_sanitized_initial.push_back(molar_abundance); + } + const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9) + const double rho = netIn.density; // Density in g/cm^3 + + const auto indices = packSpeciesTypeIndexVectors(Y_sanitized_initial, T9, rho); + Eigen::VectorXd Y_QSE; + try { + Y_QSE = calculateSteadyStateAbundances(Y_sanitized_initial, T9, rho, indices); + } catch (const std::runtime_error& e) { + LOG_ERROR(m_logger, "Failed to calculate steady state abundances. Aborting QSE evaluation."); + m_logger->flush_log(); + throw std::runtime_error("Failed to calculate steady state abundances: " + std::string(e.what())); + } + + state_type YDynamic_ublas(indices.dynamicSpeciesIndices.size() + 1); + for (size_t i = 0; i < indices.dynamicSpeciesIndices.size(); ++i) { + YDynamic_ublas(i) = Y_sanitized_initial[indices.dynamicSpeciesIndices[i]]; + } + YDynamic_ublas(indices.dynamicSpeciesIndices.size()) = 0.0; // Placeholder for specific energy rate + + const RHSFunctor rhs_functor(m_engine, indices.dynamicSpeciesIndices, indices.QSESpeciesIndices, Y_QSE, T9, rho); + const auto stepper = make_controlled>(1.0e-8, 1.0e-8); + + size_t stepCount = integrate_adaptive( + stepper, + rhs_functor, + YDynamic_ublas, + 0.0, // Start time + netIn.tMax, + netIn.dt0 + ); + + std::vector YFinal = Y_sanitized_initial; + for (size_t i = 0; i speciesNames(m_engine.getNetworkSpecies().size()); + std::vector finalMassFractions(m_engine.getNetworkSpecies().size()); + double massFractionSum = 0.0; + + for (size_t i = 0; i < speciesNames.size(); ++i) { + const auto& species = m_engine.getNetworkSpecies()[i]; + speciesNames[i] = species.name(); + finalMassFractions[i] = YFinal[i] * species.mass(); // Convert from molar abundance to mass fraction + massFractionSum += finalMassFractions[i]; + } + for (auto& mf : finalMassFractions) { + mf /= massFractionSum; // Normalize to get mass fractions + } + + fourdst::composition::Composition outputComposition(speciesNames, finalMassFractions); + NetOut netOut; + netOut.composition = outputComposition; + netOut.energy = finalSpecificEnergyRate; // Specific energy rate + netOut.num_steps = stepCount; + return netOut; + } + + dynamicQSESpeciesIndices QSENetworkSolver::packSpeciesTypeIndexVectors( + const std::vector& Y, + const double T9, + const double rho + ) const { + constexpr double timescaleCutoff = 1.0e-5; + constexpr double abundanceCutoff = 1.0e-15; + + LOG_INFO(m_logger, "Partitioning species using T9={:0.2f} and ρ={:0.2e}", T9, rho); + LOG_INFO(m_logger, "Timescale Cutoff: {:.1e} s, Abundance Cutoff: {:.1e}", timescaleCutoff, abundanceCutoff); + + std::vectordynamicSpeciesIndices; // Slow species that are not in QSE + std::vectorQSESpeciesIndices; // Fast species that are in QSE + + std::unordered_map speciesTimescale = m_engine.getSpeciesTimescales(Y, T9, rho); + + for (size_t i = 0; i < m_engine.getNetworkSpecies().size(); ++i) { + const auto& species = m_engine.getNetworkSpecies()[i]; + const double timescale = speciesTimescale[species]; + const double abundance = Y[i]; + + if (std::isinf(timescale) || abundance < abundanceCutoff || timescale <= timescaleCutoff) { + QSESpeciesIndices.push_back(i); + } else { + dynamicSpeciesIndices.push_back(i); + } + } + LOG_INFO(m_logger, "Partitioning complete. Dynamical species: {}, QSE species: {}.", dynamicSpeciesIndices.size(), QSESpeciesIndices.size()); + LOG_INFO(m_logger, "Dynamic species: {}", [dynamicSpeciesIndices](const DynamicEngine& engine_wrapper) -> std::string { + std::string result; + int count = 0; + for (const auto& i : dynamicSpeciesIndices) { + result += std::string(engine_wrapper.getNetworkSpecies()[i].name()); + if (count < dynamicSpeciesIndices.size() - 2) { + result += ", "; + } else if (count == dynamicSpeciesIndices.size() - 2) { + result += " and "; + } + count++; + } + return result; + }(m_engine)); + LOG_INFO(m_logger, "QSE species: {}", [QSESpeciesIndices](const DynamicEngine& engine_wrapper) -> std::string { + std::string result; + int count = 0; + for (const auto& i : QSESpeciesIndices) { + result += std::string(engine_wrapper.getNetworkSpecies()[i].name()); + if (count < QSESpeciesIndices.size() - 2) { + result += ", "; + } else if (count == QSESpeciesIndices.size() - 2) { + result += " and "; + } + count++; + } + return result; + }(m_engine)); + return {dynamicSpeciesIndices, QSESpeciesIndices}; + } + + Eigen::VectorXd QSENetworkSolver::calculateSteadyStateAbundances( + const std::vector &Y, + const double T9, + const double rho, + const dynamicQSESpeciesIndices &indices + ) const { + std::vector Y_dyn; + Eigen::VectorXd Y_qse_initial(indices.QSESpeciesIndices.size()); + for (const auto& i : indices.dynamicSpeciesIndices) { Y_dyn.push_back(Y[i]); } + for (size_t i = 0; i < indices.QSESpeciesIndices.size(); ++i) { + Y_qse_initial(i) = Y[indices.QSESpeciesIndices[i]]; + if (Y_qse_initial(i) < 1.0e-99) { Y_qse_initial(i) = 1.0e-99; } + } + + Eigen::VectorXd v_qse = Y_qse_initial.array().log(); + + EigenFunctor qse_problem(m_engine, Y_dyn, indices.dynamicSpeciesIndices, indices.QSESpeciesIndices, T9, rho); + LOG_INFO(m_logger, "--- QSE Pre-Solve Diagnostics ---"); + + Eigen::VectorXd f_initial(indices.QSESpeciesIndices.size()); + qse_problem(v_qse, f_initial); + LOG_INFO(m_logger, "Initial Guess ||f||: {:0.4e}", f_initial.norm()); + + Eigen::MatrixXd J_initial(indices.QSESpeciesIndices.size(), indices.QSESpeciesIndices.size()); + qse_problem.df(v_qse, J_initial); + const Eigen::JacobiSVD svd(J_initial); + double cond = svd.singularValues().maxCoeff() / svd.singularValues().minCoeff(); + LOG_INFO(m_logger, "Initial Jacobian Condition Number: {:0.4e}", cond); + LOG_INFO(m_logger, "Starting QSE solve..."); + + Eigen::HybridNonLinearSolver> solver(qse_problem); + solver.parameters.xtol = 1.0e-8; // Set tolerance + + // 5. Run the solver. It will modify v_qse in place. + const int eigenStatus = solver.solve(v_qse); + + // 6. Check for convergence and return the result + if(eigenStatus != Eigen::Success) { + LOG_WARNING(m_logger, "--- QSE SOLVER FAILED ---"); + LOG_WARNING(m_logger, "Eigen status code: {}", eigenStatus); + LOG_WARNING(m_logger, "Iterations performed: {}", solver.iter); + + // Log the final state that caused the failure + Eigen::VectorXd Y_qse_final_fail = v_qse.array().exp(); + for(long i=0; i Y_qse[{}]: {:0.4e}", i, v_qse(i), i, Y_qse_final_fail(i)); + } + + // Log the residual at the final state + Eigen::VectorXd f_final(indices.QSESpeciesIndices.size()); + qse_problem(v_qse, f_final); + LOG_WARNING(m_logger, "Final ||f||: {:0.4e}", f_final.norm()); + + throw std::runtime_error("Eigen QSE solver did not converge."); + } + LOG_INFO(m_logger, "Eigen QSE solver converged in {} iterations.", solver.iter); + + return v_qse.array().exp(); + + } + + NetOut QSENetworkSolver::initializeNetworkWithShortIgnition(const NetIn &netIn) const { + const auto ignitionTemperature = m_config.get( + "gridfire:solver:QSE:ignition:temperature", + 2e8 + ); // 0.2 GK + const auto ignitionDensity = m_config.get( + "gridfire:solver:QSE:ignition:density", + 1e6 + ); // 1e6 g/cm^3 + const auto ignitionTime = m_config.get( + "gridfire:solver:QSE:ignition:tMax", + 1e-7 + ); // 0.1 μs + const auto ignitionStepSize = m_config.get( + "gridfire:solver:QSE:ignition:dt0", + 1e-15 + ); // 1e-15 seconds + + LOG_INFO( + m_logger, + "Igniting network with T={:<5.3E}, ρ={:<5.3E}, tMax={:<5.3E}, dt0={:<5.3E}...", + ignitionTemperature, + ignitionDensity, + ignitionTime, + ignitionStepSize + ); + + NetIn preIgnition = netIn; + preIgnition.temperature = ignitionTemperature; + preIgnition.density = ignitionDensity; + preIgnition.tMax = ignitionTime; + preIgnition.dt0 = ignitionStepSize; + + DirectNetworkSolver ignitionSolver(m_engine); + NetOut postIgnition = ignitionSolver.evaluate(preIgnition); + LOG_INFO(m_logger, "Network ignition completed in {} steps.", postIgnition.num_steps); + return postIgnition; + } + + void QSENetworkSolver::RHSFunctor::operator()( + const boost::numeric::ublas::vector &YDynamic, + boost::numeric::ublas::vector &dYdtDynamic, + double t + ) const { + // --- Populate the slow / dynamic species vector --- + std::vector YFull(m_engine.getNetworkSpecies().size()); + for (size_t i = 0; i < m_dynamicSpeciesIndices.size(); ++i) { + YFull[m_dynamicSpeciesIndices[i]] = YDynamic(i); + } + + // --- Populate the QSE species vector --- + for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { + YFull[m_QSESpeciesIndices[i]] = m_Y_QSE(i); + } + + auto [full_dYdt, specificEnergyRate] = m_engine.calculateRHSAndEnergy(YFull, m_T9, m_rho); + + dYdtDynamic.resize(m_dynamicSpeciesIndices.size() + 1); + for (size_t i = 0; i < m_dynamicSpeciesIndices.size(); ++i) { + dYdtDynamic(i) = full_dYdt[m_dynamicSpeciesIndices[i]]; + } + dYdtDynamic[m_dynamicSpeciesIndices.size()] = specificEnergyRate; + } + + NetOut DirectNetworkSolver::evaluate(const NetIn &netIn) { + namespace ublas = boost::numeric::ublas; + namespace odeint = boost::numeric::odeint; + using fourdst::composition::Composition; + + const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9) + const unsigned long numSpecies = m_engine.getNetworkSpecies().size(); + + const auto absTol = m_config.get("gridfire:solver:DirectNetworkSolver:absTol", 1.0e-8); + const auto relTol = m_config.get("gridfire:solver:DirectNetworkSolver:relTol", 1.0e-8); + + size_t stepCount = 0; + + RHSFunctor rhsFunctor(m_engine, T9, netIn.density); + JacobianFunctor jacobianFunctor(m_engine, T9, netIn.density); + + ublas::vector Y(numSpecies + 1); + + for (size_t i = 0; i < numSpecies; ++i) { + const auto& species = m_engine.getNetworkSpecies()[i]; + try { + Y(i) = netIn.composition.getMolarAbundance(std::string(species.name())); + } catch (const std::runtime_error) { + LOG_DEBUG(m_logger, "Species '{}' not found in composition. Setting abundance to 0.0.", species.name()); + Y(i) = 0.0; + } + } + Y(numSpecies) = 0.0; + + const auto stepper = odeint::make_controlled>(absTol, relTol); + stepCount = odeint::integrate_adaptive( + stepper, + std::make_pair(rhsFunctor, jacobianFunctor), + Y, + 0.0, + netIn.tMax, + netIn.dt0 + ); + + std::vector finalMassFractions(numSpecies); + for (size_t i = 0; i < numSpecies; ++i) { + const double molarMass = m_engine.getNetworkSpecies()[i].mass(); + finalMassFractions[i] = Y(i) * molarMass; // Convert from molar abundance to mass fraction + if (finalMassFractions[i] < MIN_ABUNDANCE_THRESHOLD) { + finalMassFractions[i] = 0.0; + } + } + + std::vector speciesNames; + speciesNames.reserve(numSpecies); + for (const auto& species : m_engine.getNetworkSpecies()) { + speciesNames.push_back(std::string(species.name())); + } + + Composition outputComposition(speciesNames, finalMassFractions); + outputComposition.finalize(true); + + NetOut netOut; + netOut.composition = std::move(outputComposition); + netOut.energy = Y(numSpecies); // Specific energy rate + netOut.num_steps = stepCount; + + return netOut; + } + + void DirectNetworkSolver::RHSFunctor::operator()( + const boost::numeric::ublas::vector &Y, + boost::numeric::ublas::vector &dYdt, + double t + ) const { + const std::vector y(Y.begin(), m_numSpecies + Y.begin()); + auto [dydt, eps] = m_engine.calculateRHSAndEnergy(y, m_T9, m_rho); + dYdt.resize(m_numSpecies + 1); + std::ranges::copy(dydt, dydt.begin()); + dYdt(m_numSpecies) = eps; + } + + void DirectNetworkSolver::JacobianFunctor::operator()( + const boost::numeric::ublas::vector &Y, + boost::numeric::ublas::matrix &J, + double t, + boost::numeric::ublas::vector &dfdt + ) const { + J.resize(m_numSpecies + 1, m_numSpecies + 1); + J.clear(); + for (int i = 0; i < m_numSpecies + 1; ++i) { + for (int j = 0; j < m_numSpecies + 1; ++j) { + J(i, j) = m_engine.getJacobianMatrixEntry(i, j); + } + } + } +} \ No newline at end of file diff --git a/src/network/meson.build b/src/network/meson.build index fb5be9a8..5b68487f 100644 --- a/src/network/meson.build +++ b/src/network/meson.build @@ -1,8 +1,10 @@ # Define the library network_sources = files( 'lib/network.cpp', - 'lib/approx8.cpp', - 'lib/netgraph.cpp' + 'lib/engine/engine_approx8.cpp', + 'lib/engine/engine_graph.cpp', + 'lib/reaction/reaction.cpp', + 'lib/solver/solver.cpp', ) @@ -13,7 +15,9 @@ dependencies = [ composition_dep, reaclib_reactions_dep, cppad_dep, - log_dep + log_dep, + xxhash_dep, + eigen_dep, ] # Define the libnetwork library so it can be linked against by other parts of the build system @@ -33,7 +37,10 @@ network_dep = declare_dependency( # Make headers accessible network_headers = files( 'include/gridfire/network.h', - 'include/gridfire/approx8.h', - 'include/gridfire/netgraph.h', + 'include/gridfire/engine/engine_abstract.h', + 'include/gridfire/engine/engine_approx8.h', + 'include/gridfire/engine/engine_graph.h', + 'include/gridfire/reaction/reaction.h', + 'include/gridfire/solver/solver.h', ) -install_headers(network_headers, subdir : 'gridfire/gridfire') +install_headers(network_headers, subdir : 'gridfire') diff --git a/subprojects/eigen.wrap b/subprojects/eigen.wrap new file mode 100644 index 00000000..becc4767 --- /dev/null +++ b/subprojects/eigen.wrap @@ -0,0 +1,13 @@ +[wrap-file] +directory = eigen-3.4.0 +source_url = https://gitlab.com/libeigen/eigen/-/archive/3.4.0/eigen-3.4.0.tar.bz2 +source_filename = eigen-3.4.0.tar.bz2 +source_hash = b4c198460eba6f28d34894e3a5710998818515104d6e74e5cc331ce31e46e626 +patch_filename = eigen_3.4.0-2_patch.zip +patch_url = https://wrapdb.mesonbuild.com/v2/eigen_3.4.0-2/get_patch +patch_hash = cb764fd9fec02d94aaa2ec673d473793c0d05da4f4154c142f76ef923ea68178 +source_fallback_url = https://github.com/mesonbuild/wrapdb/releases/download/eigen_3.4.0-2/eigen-3.4.0.tar.bz2 +wrapdb_version = 3.4.0-2 + +[provide] +eigen3 = eigen_dep diff --git a/subprojects/libcomposition.wrap b/subprojects/libcomposition.wrap index 2fb8e98f..766eea18 100644 --- a/subprojects/libcomposition.wrap +++ b/subprojects/libcomposition.wrap @@ -1,4 +1,4 @@ [wrap-git] url = https://github.com/4D-STAR/libcomposition.git -revision = v1.0.7 +revision = v1.0.9 depth = 1 \ No newline at end of file diff --git a/tests/graphnet_sandbox/approx8 b/tests/graphnet_sandbox/approx8 new file mode 100644 index 00000000..e69de29b diff --git a/tests/graphnet_sandbox/main.cpp b/tests/graphnet_sandbox/main.cpp new file mode 100644 index 00000000..5f06ded9 --- /dev/null +++ b/tests/graphnet_sandbox/main.cpp @@ -0,0 +1,48 @@ +#define QUILL_COMPILE_ACTIVE_LOG_LEVEL QUILL_COMPILE_ACTIVE_LOG_LEVEL_DEBUG +#include +#include "../../src/network/include/gridfire/engine/engine_graph.h" +#include "gridfire/network.h" +#include "../../src/network/include/gridfire/engine/engine_approx8.h" +#include "../../src/network/include/gridfire/solver/solver.h" + +#include "fourdst/composition/composition.h" + +int main() { + using namespace gridfire; + // reaction::REACLIBLogicalReactionSet reactions = build_reaclib_nuclear_network_from_file("approx8", false); + // GraphNetwork network(reactions); + const std::vector comp = {0.708, 0.0, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4}; + const std::vector symbols = {"H-1", "H-2", "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); + composition.finalize(true); + + + NetIn netIn; + netIn.composition = composition; + netIn.temperature = 1e7; + netIn.density = 1e2; + netIn.energy = 0.0; + + netIn.tMax = 3.15e17; + netIn.dt0 = 1e12; + + NetOut netOut; + approx8::Approx8Network approx8Network; + netOut = approx8Network.evaluate(netIn); + std::cout << "Approx8 Network Output: " << netOut << std::endl; + + // // Network Ignition + netIn.dt0 = 1e-15; + // netIn.temperature = 2.0e8; + // netIn.density = 1.0e6; + // netIn.tMax = 1e-1; + + GraphEngine ReaclibEngine(composition); + // netOut = network.evaluate(netIn); + solver::QSENetworkSolver solver(ReaclibEngine); + netOut = solver.evaluate(netIn); + std::cout << "QSE Graph Network Output: " << netOut << std::endl; +} \ No newline at end of file diff --git a/tests/graphnet_sandbox/meson.build b/tests/graphnet_sandbox/meson.build new file mode 100644 index 00000000..1c823c61 --- /dev/null +++ b/tests/graphnet_sandbox/meson.build @@ -0,0 +1,5 @@ +executable( + 'graphnet_sandbox', + 'main.cpp', + dependencies: [network_dep, composition_dep], +) diff --git a/tests/meson.build b/tests/meson.build index 11444fea..720bd28e 100644 --- a/tests/meson.build +++ b/tests/meson.build @@ -5,3 +5,4 @@ gtest_nomain_dep = dependency('gtest', main: false, required : true) # Subdirectories for unit and integration tests subdir('network') +subdir('graphnet_sandbox') diff --git a/tests/network/approx8Test.cpp b/tests/network/approx8Test.cpp index 7388ef97..1848e096 100644 --- a/tests/network/approx8Test.cpp +++ b/tests/network/approx8Test.cpp @@ -3,10 +3,9 @@ #include "fourdst/composition/composition.h" #include "fourdst/config/config.h" -#include "gridfire/approx8.h" -#include "gridfire/netgraph.h" +#include "../../src/network/include/gridfire/engine/engine_approx8.h" +#include "../../src/network/include/gridfire/engine/engine_graph.h" #include "gridfire/network.h" -#include "gridfire/reaclib.h" #include @@ -86,7 +85,7 @@ TEST_F(approx8Test, reaclib) { netIn.tMax = 3.15e17; netIn.dt0 = 1e12; - GraphNetwork network(composition); + GraphEngine network(composition); NetOut netOut; netOut = network.evaluate(netIn); std::cout << netOut << std::endl;