From 56f93420527912cf94da1d52bd0d5c32bd0311a9 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 31 Oct 2025 07:38:04 -0400 Subject: [PATCH] fix(AtomicWeakRate): added bool overrides for sparsity calculations these allow the atomic weak rate to be used with the current sparsity code. Note also that this commit reintroduces the full weak reaction set which does introduce an enormouse degree of stiffness which will need to be filtered out. What this means is that this commit practically cannot evolve a network due to this stiffness --- src/include/gridfire/reaction/weak/weak.h | 16 ++++- src/lib/engine/engine_graph.cpp | 18 ++--- src/lib/engine/procedures/construction.cpp | 84 +++++++++++----------- src/lib/reaction/weak/weak.cpp | 52 ++++++++++++-- 4 files changed, 113 insertions(+), 57 deletions(-) diff --git a/src/include/gridfire/reaction/weak/weak.h b/src/include/gridfire/reaction/weak/weak.h index d0865352..cec8daee 100644 --- a/src/include/gridfire/reaction/weak/weak.h +++ b/src/include/gridfire/reaction/weak/weak.h @@ -473,6 +473,20 @@ namespace gridfire::rates::weak { CppAD::vector>& st ) override; + bool for_sparse_jac( + size_t q, + const CppAD::vector &r, + CppAD::vector &s, + const CppAD::vector &x + ) override; + + bool rev_sparse_jac( + size_t q, + const CppAD::vector &rt, + CppAD::vector &st, + const CppAD::vector &x + ) override; + private: const WeakRateInterpolator& m_interpolator; const size_t m_a; @@ -536,7 +550,7 @@ namespace gridfire::rates::weak { T rateConstant = static_cast(0.0); if constexpr (std::is_same_v>) { // Case where T is an AD type - std::vector ax = {T9, log_rhoYe, mue}; + std::vector ax = {T9, log_rhoYe}; std::vector ay(2); m_atomic(ax, ay); rateConstant = static_cast(ay[0]); diff --git a/src/lib/engine/engine_graph.cpp b/src/lib/engine/engine_graph.cpp index 092a2724..b400e740 100644 --- a/src/lib/engine/engine_graph.cpp +++ b/src/lib/engine/engine_graph.cpp @@ -1407,7 +1407,7 @@ namespace gridfire { } bool GraphEngine::AtomicReverseRate::reverse( - size_t q, + const size_t q, const CppAD::vector &tx, const CppAD::vector &ty, CppAD::vector &px, @@ -1424,7 +1424,7 @@ namespace gridfire { } bool GraphEngine::AtomicReverseRate::for_sparse_jac( - size_t q, + const size_t q, const CppAD::vector> &r, CppAD::vector> &s ) { @@ -1433,7 +1433,7 @@ namespace gridfire { } bool GraphEngine::AtomicReverseRate::rev_sparse_jac( - size_t q, + const size_t q, const CppAD::vector> &rt, CppAD::vector> &st ) { @@ -1442,7 +1442,7 @@ namespace gridfire { } bool GraphEngine::AtomicReverseRate::for_sparse_jac( - size_t q, + const size_t q, const CppAD::vector &r, CppAD::vector &s, const CppAD::vector &x @@ -1450,8 +1450,8 @@ namespace gridfire { constexpr size_t n = 1; constexpr size_t m = 1; - CPPAD_ASSERT_KNOWN(r.size() == n * q, "for_sparse_jac: 'r' size is incorrect."); - CPPAD_ASSERT_KNOWN(s.size() == m * q, "for_sparse_jac: 's' size is incorrect."); + CPPAD_ASSERT_KNOWN(r.size() == n * q, "AtomicReverseRate::for_sparse_jac: 'r' size is incorrect."); + CPPAD_ASSERT_KNOWN(s.size() == m * q, "AtomicReverseRate::for_sparse_jac: 's' size is incorrect."); // S = R for (size_t j = 0; j < q; j++) { @@ -1463,7 +1463,7 @@ namespace gridfire { } bool GraphEngine::AtomicReverseRate::rev_sparse_jac( - size_t q, + const size_t q, const CppAD::vector &rt, CppAD::vector &st, const CppAD::vector &x @@ -1471,8 +1471,8 @@ namespace gridfire { constexpr size_t n = 1; constexpr size_t m = 1; - CPPAD_ASSERT_KNOWN(rt.size() == n * q, "for_sparse_jac: 'r' size is incorrect."); - CPPAD_ASSERT_KNOWN(st.size() == m * q, "for_sparse_jac: 's' size is incorrect."); + CPPAD_ASSERT_KNOWN(rt.size() == n * q, "AtomicReverseRate::for_sparse_jac: 'r' size is incorrect."); + CPPAD_ASSERT_KNOWN(st.size() == m * q, "AtomicReverseRate::for_sparse_jac: 's' size is incorrect."); // st = rt for (size_t j = 0; j < q; j++) { diff --git a/src/lib/engine/procedures/construction.cpp b/src/lib/engine/procedures/construction.cpp index dfabe8c1..f1957816 100644 --- a/src/lib/engine/procedures/construction.cpp +++ b/src/lib/engine/procedures/construction.cpp @@ -76,48 +76,48 @@ namespace gridfire { } // --- Clone all possible weak reactions into the master reaction pool --- - // for (const auto& parent_species: weakInterpolator.available_isotopes()) { - // std::expected upProduct = fourdst::atomic::az_to_species( - // parent_species.a(), - // parent_species.z() + 1 - // ); - // std::expected downProduct = fourdst::atomic::az_to_species( - // parent_species.a(), - // parent_species.z() - 1 - // ); - // if (downProduct.has_value()) { // Only add the reaction if the Species map contains the product - // master_reaction_pool.add_reaction( - // std::make_unique( - // parent_species, - // rates::weak::WeakReactionType::BETA_PLUS_DECAY, - // weakInterpolator - // ) - // ); - // master_reaction_pool.add_reaction( - // std::make_unique( - // parent_species, - // rates::weak::WeakReactionType::ELECTRON_CAPTURE, - // weakInterpolator - // ) - // ); - // } - // if (upProduct.has_value()) { // Only add the reaction if the Species map contains the product - // master_reaction_pool.add_reaction( - // std::make_unique( - // parent_species, - // rates::weak::WeakReactionType::BETA_MINUS_DECAY, - // weakInterpolator - // ) - // ); - // master_reaction_pool.add_reaction( - // std::make_unique( - // parent_species, - // rates::weak::WeakReactionType::POSITRON_CAPTURE, - // weakInterpolator - // ) - // ); - // } - // } // TODO: Remove comments, weak reactions have been disabled for testing + for (const auto& parent_species: weakInterpolator.available_isotopes()) { + std::expected upProduct = fourdst::atomic::az_to_species( + parent_species.a(), + parent_species.z() + 1 + ); + std::expected downProduct = fourdst::atomic::az_to_species( + parent_species.a(), + parent_species.z() - 1 + ); + if (downProduct.has_value()) { // Only add the reaction if the Species map contains the product + master_reaction_pool.add_reaction( + std::make_unique( + parent_species, + rates::weak::WeakReactionType::BETA_PLUS_DECAY, + weakInterpolator + ) + ); + master_reaction_pool.add_reaction( + std::make_unique( + parent_species, + rates::weak::WeakReactionType::ELECTRON_CAPTURE, + weakInterpolator + ) + ); + } + if (upProduct.has_value()) { // Only add the reaction if the Species map contains the product + master_reaction_pool.add_reaction( + std::make_unique( + parent_species, + rates::weak::WeakReactionType::BETA_MINUS_DECAY, + weakInterpolator + ) + ); + master_reaction_pool.add_reaction( + std::make_unique( + parent_species, + rates::weak::WeakReactionType::POSITRON_CAPTURE, + weakInterpolator + ) + ); + } + } // TODO: Remove comments, weak reactions have been disabled for testing // --- Step 2: Use non-owning raw pointers for the fast build algorithm --- std::vector remainingReactions; diff --git a/src/lib/reaction/weak/weak.cpp b/src/lib/reaction/weak/weak.cpp index 0d0713f2..dd4bd9b6 100644 --- a/src/lib/reaction/weak/weak.cpp +++ b/src/lib/reaction/weak/weak.cpp @@ -382,7 +382,7 @@ namespace gridfire::rates::weak { ) const { const CppAD::AD log_rhoYe = CppAD::log10(rho * Ye); - const std::vector> ax = {T9, log_rhoYe, mue}; + const std::vector> ax = {T9, log_rhoYe}; std::vector> ay(2); // 2 outputs are the reaction rate (1/s) and the neutrino loss (MeV) m_atomic(ax, ay); // Note: We needed to make m_atomic mutable to allow this call in a const method. @@ -491,7 +491,7 @@ namespace gridfire::rates::weak { return logNeutrinoLoss; } - // Note that the input vector tx is of size 3: [T9, log10(rho*Ye), mu_e] + // Note that the input vector tx is of size 2: [T9, log10(rho*Ye)] bool WeakReaction::AtomicWeakRate::forward ( const size_t p, const size_t q, @@ -549,7 +549,7 @@ namespace gridfire::rates::weak { } bool WeakReaction::AtomicWeakRate::reverse( - size_t q, + const size_t q, const CppAD::vector &tx, const CppAD::vector &ty, CppAD::vector &px, @@ -624,7 +624,7 @@ namespace gridfire::rates::weak { } bool WeakReaction::AtomicWeakRate::for_sparse_jac( - size_t q, + const size_t q, const CppAD::vector > &r, CppAD::vector > &s ) { @@ -640,7 +640,7 @@ namespace gridfire::rates::weak { } bool WeakReaction::AtomicWeakRate::rev_sparse_jac( - size_t q, + const size_t q, const CppAD::vector > &rt, CppAD::vector > &st ) { @@ -655,11 +655,53 @@ namespace gridfire::rates::weak { } + bool WeakReaction::AtomicWeakRate::for_sparse_jac( + const size_t q, + const CppAD::vector &r, + CppAD::vector &s, + const CppAD::vector &x + ) { + constexpr size_t n = 2; // Number of inputs + constexpr size_t m = 2; // Number of outputs + CPPAD_ASSERT_KNOWN(r.size() == q * n, "AtomicWeakRate::for_sparse_jac: 'r' size is incorrect!"); + CPPAD_ASSERT_KNOWN(s.size() == q * m, "AtomicWeakRate::for_sparse_jac: 's' size is incorrect!"); + // Both outputs depend on both inputs + // s[i + j*m] represents s(i,j) - output i, direction j + // r[k + j*n] represents r(k,j) - input k, direction j + for (size_t j = 0; j < q; j++) { + // s(0,j) = r(0,j) || r(1,j) --- output 0 depends on both inputs + s[0 + j*m] = r[0 + j*n] || r[1 + j*n]; + // s(1,j) = r(0,j) || r(1,j) --- output 1 depends on both inputs + s[1 + j*m] = r[0 + j*n] || r[1 + j*n]; + } + return true; + } + bool WeakReaction::AtomicWeakRate::rev_sparse_jac( + const size_t q, + const CppAD::vector &rt, + CppAD::vector &st, + const CppAD::vector &x + ) { + constexpr size_t n = 2; // Number of inputs + constexpr size_t m = 2; // Number of outputs + CPPAD_ASSERT_KNOWN(rt.size() == q * m, "AtomicWeakRate::rev_sparse_jac: 'rt' size is incorrect!"); + CPPAD_ASSERT_KNOWN(st.size() == q * n, "AtomicWeakRate::rev_sparse_jac: 'st' size is incorrect!"); + + for (size_t j = 0; j < q; j++) { + //st(0,j) = rt(0,j) || rt(1,j) --- input 0 affects both outputs + st[0 + j*n] = rt[0 + j*m] || rt[1 + j*m]; + + //st(1,j) = rt(0,j) || rt(1,j) --- input 1 affects both outputs + st[1 + j*n] = rt[0 + j*m] || rt[1 + j*m]; + } + + return true; + } }