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; + } }