diff --git a/src/include/gridfire/engine/engine_graph.h b/src/include/gridfire/engine/engine_graph.h index edbdfcbd..692f7459 100644 --- a/src/include/gridfire/engine/engine_graph.h +++ b/src/include/gridfire/engine/engine_graph.h @@ -25,6 +25,7 @@ #include "cppad/cppad.hpp" #include "cppad/utility/sparse_rc.hpp" #include "cppad/speed/sparse_jac_fun.hpp" + #include "gridfire/reaction/weak/weak_interpolator.h" #include "gridfire/reaction/weak/weak_rate_library.h" @@ -729,6 +730,8 @@ namespace gridfire { BuildDepthType depth ) override; + void lumpReactions(); + private: struct PrecomputedReaction { @@ -754,6 +757,20 @@ namespace gridfire { const double c = Constants::getInstance().get("c").value; ///< Speed of light in cm/s. const double kB = Constants::getInstance().get("kB").value; ///< Boltzmann constant in erg/K. }; + + enum class JacobianMatrixState { + UNINITIALIZED, + STALE, + READY_DENSE, + READY_SPARSE + }; + + std::unordered_map m_jacobianMatrixStateNameMap = { + {JacobianMatrixState::UNINITIALIZED, "Uninitialized"}, + {JacobianMatrixState::STALE, "Stale"}, + {JacobianMatrixState::READY_DENSE, "Ready (dense)"}, + {JacobianMatrixState::READY_SPARSE, "Ready (sparse)"}, + }; private: class AtomicReverseRate final : public CppAD::atomic_base { public: @@ -790,6 +807,18 @@ namespace gridfire { const CppAD::vector>&rt, 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 reaction::Reaction& m_reaction; @@ -811,10 +840,11 @@ namespace gridfire { std::unordered_map m_speciesToIndexMap; ///< Map from species to their index in the stoichiometry matrix. std::unordered_map m_indexToSpeciesMap; ///< Map from index to species in the stoichiometry matrix. - boost::numeric::ublas::compressed_matrix m_stoichiometryMatrix; ///< Stoichiometry matrix (species x reactions). mutable boost::numeric::ublas::compressed_matrix m_jacobianMatrix; ///< Jacobian matrix (species x species). + mutable JacobianMatrixState m_jacobianMatrixState = JacobianMatrixState::UNINITIALIZED; + mutable CppAD::ADFun m_rhsADFun; ///< CppAD function for the right-hand side of the ODE. mutable CppAD::ADFun m_epsADFun; ///< CppAD function for the energy generation rate. mutable CppAD::sparse_jac_work m_jac_work; ///< Work object for sparse Jacobian calculations. @@ -827,7 +857,6 @@ namespace gridfire { std::unique_ptr m_screeningModel = screening::selectScreeningModel(m_screeningType); bool m_usePrecomputation = true; ///< Flag to enable or disable using precomputed reactions for efficiency. Mathematically, this should not change the results. Generally end users should not need to change this. - bool m_useReverseReactions = true; ///< Flag to enable or disable reverse reactions. If false, only forward reactions are considered. BuildDepthType m_depth; @@ -978,46 +1007,8 @@ namespace gridfire { std::function(const fourdst::atomic::Species &)> speciesLookup, const std::function& reactionLookup ) const; - - // /** - // * @brief Calculates all derivatives (dY/dt) and the energy generation rate (double precision). - // * - // * @param Y Vector of molar abundances for all species in the network. - // * @param T9 Temperature in units of 10^9 K. - // * @param rho Density in g/cm^3. - // * @return StepDerivatives containing dY/dt and energy generation rate. - // * - // * This method calculates the time derivatives of all species and the - // * specific nuclear energy generation rate for the current state using - // * double precision arithmetic. - // */ - // [[nodiscard]] StepDerivatives calculateAllDerivatives( - // const std::vector& Y, - // double T9, - // double rho - // ) const; - // - // /** - // * @brief Calculates all derivatives (dY/dt) and the energy generation rate (automatic differentiation). - // * - // * @param Y Vector of molar abundances for all species in the network. - // * @param T9 Temperature in units of 10^9 K. - // * @param rho Density in g/cm^3. - // * @return StepDerivatives containing dY/dt and energy generation rate. - // * - // * This method calculates the time derivatives of all species and the - // * specific nuclear energy generation rate for the current state using - // * automatic differentiation. - // */ - // [[nodiscard]] StepDerivatives calculateAllDerivatives( - // const std::vector &Y, - // ADDouble T9, - // ADDouble rho - // ) const; }; - - template T GraphEngine::calculateReverseMolarReactionFlow( T T9, diff --git a/src/include/gridfire/exceptions/error_engine.h b/src/include/gridfire/exceptions/error_engine.h index b4fecd84..3937a760 100644 --- a/src/include/gridfire/exceptions/error_engine.h +++ b/src/include/gridfire/exceptions/error_engine.h @@ -110,4 +110,37 @@ namespace gridfire::exceptions { std::string m_message; }; + class JacobianError : public EngineError {}; + + class StaleJacobianError final : public JacobianError { + public: + explicit StaleJacobianError(std::string message) : m_message(std::move(message)) {} + [[nodiscard]] const char* what() const noexcept override { + return m_message.c_str(); + } + + private: + std::string m_message; + }; + + class UninitializedJacobianError final: public JacobianError { + public: + explicit UninitializedJacobianError(std::string message): m_message(std::move(message)) {} + [[nodiscard]] const char* what() const noexcept override { + return m_message.c_str(); + } + private: + std::string m_message; + }; + + class UnknownJacobianError final : public JacobianError { + public: + explicit UnknownJacobianError(std::string message): m_message(std::move(message)) {} + [[nodiscard]] const char* what() const noexcept override { + return m_message.c_str(); + } + private: + std::string m_message; + }; + } \ No newline at end of file diff --git a/src/include/gridfire/reaction/reaction.h b/src/include/gridfire/reaction/reaction.h index cb0401a0..c8d43484 100644 --- a/src/include/gridfire/reaction/reaction.h +++ b/src/include/gridfire/reaction/reaction.h @@ -1,5 +1,6 @@ #pragma once +#include #include #include @@ -739,7 +740,7 @@ namespace gridfire::reaction { rate.a5 * T953 + rate.a6 * logT9; sum += CppAD::exp(exponent); - // return sum; // TODO: REMOVE THIS ITS FOR TESTING ONLY + // return sum; // TODO: REMOVE OR COMMENT THIS. ITS FOR TESTING ONLY } return sum; } diff --git a/src/include/gridfire/solver/strategies/CVODE_solver_strategy.h b/src/include/gridfire/solver/strategies/CVODE_solver_strategy.h index 6917a1e7..34c891bb 100644 --- a/src/include/gridfire/solver/strategies/CVODE_solver_strategy.h +++ b/src/include/gridfire/solver/strategies/CVODE_solver_strategy.h @@ -257,7 +257,7 @@ namespace gridfire::solver { * sorted table of species with the highest error ratios; then invokes diagnostic routines to * inspect Jacobian stiffness and species balance. */ - void log_step_diagnostics(const CVODEUserData& user_data) const; + void log_step_diagnostics(const CVODEUserData& user_data, bool displayJacobianStiffness) const; private: SUNContext m_sun_ctx = nullptr; ///< SUNDIALS context (lifetime of the solver). void* m_cvode_mem = nullptr; ///< CVODE memory block. diff --git a/src/include/gridfire/utils/table_format.h b/src/include/gridfire/utils/table_format.h index 7df92e9e..c5a21848 100644 --- a/src/include/gridfire/utils/table_format.h +++ b/src/include/gridfire/utils/table_format.h @@ -82,7 +82,7 @@ namespace gridfire::utils { // --- Helper to draw horizontal border --- auto draw_border = [&]() { table_ss << "+"; - for (size_t width : col_widths) { + for (const size_t width : col_widths) { table_ss << std::string(width + 2, '-'); // +2 for padding table_ss << "+"; } diff --git a/src/lib/engine/engine_graph.cpp b/src/lib/engine/engine_graph.cpp index 14777ff2..25ff7068 100644 --- a/src/lib/engine/engine_graph.cpp +++ b/src/lib/engine/engine_graph.cpp @@ -6,6 +6,7 @@ #include "gridfire/partition/partition_ground.h" #include "gridfire/engine/procedures/construction.h" #include "gridfire/utils/hashing.h" +#include "gridfire/utils/table_format.h" #include "fourdst/composition/species.h" #include "fourdst/composition/atomicSpecies.h" @@ -164,6 +165,7 @@ namespace gridfire { } void GraphEngine::syncInternalMaps() { + LOG_INFO(m_logger, "Synchronizing internal maps for REACLIB graph network (serif::network::GraphNetwork)..."); collectNetworkSpecies(); populateReactionIDMap(); @@ -174,15 +176,18 @@ namespace gridfire { recordADTape(); // Record the AD tape for the RHS of the ODE (dY/di and dEps/di) for all independent variables i - const size_t n = m_rhsADFun.Domain(); - const size_t m = m_rhsADFun.Range(); + const size_t inputSize = m_rhsADFun.Domain(); + const size_t outputSize = m_rhsADFun.Range(); - const std::vector select_domain(n, true); - const std::vector select_range(m, true); + // Create a range x range identity pattern + CppAD::sparse_rc> patternIn(outputSize, outputSize, outputSize); + for (size_t i = 0; i < outputSize; ++i) { + patternIn.set(i, i, i); + } + + m_rhsADFun.rev_jac_sparsity(patternIn, false, false, false, m_full_jacobian_sparsity_pattern); - m_rhsADFun.subgraph_sparsity(select_domain, select_range, false, m_full_jacobian_sparsity_pattern); m_jac_work.clear(); - m_full_sparsity_set.clear(); const auto& rows = m_full_jacobian_sparsity_pattern.row(); const auto& cols = m_full_jacobian_sparsity_pattern.col(); @@ -261,6 +266,8 @@ namespace gridfire { m_jacobianMatrix.resize(numSpecies, numSpecies, false); // Sparse matrix, no initial values LOG_TRACE_L2(m_logger, "Jacobian matrix resized to {} rows and {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2()); + + m_jacobianMatrixState = JacobianMatrixState::UNINITIALIZED; } // --- Basic Accessors and Queries --- @@ -274,13 +281,12 @@ namespace gridfire { void GraphEngine::setNetworkReactions(const reaction::ReactionSet &reactions) { m_reactions = reactions; + m_jacobianMatrixState = JacobianMatrixState::STALE; syncInternalMaps(); } 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; } @@ -509,6 +515,9 @@ namespace gridfire { } void GraphEngine::setUseReverseReactions(const bool useReverse) { + if (useReverse != m_useReverseReactions) { + m_jacobianMatrixState = JacobianMatrixState::STALE; + } m_useReverseReactions = useReverse; } @@ -568,6 +577,7 @@ namespace gridfire { if (depth != m_depth) { m_depth = depth; m_reactions = build_nuclear_network(comp, m_weakRateInterpolator, m_depth, false); + m_jacobianMatrixState = JacobianMatrixState::STALE; syncInternalMaps(); // Resync internal maps after changing the depth } else { LOG_DEBUG(m_logger, "Rebuild requested with the same depth. No changes made to the network."); @@ -729,6 +739,7 @@ namespace gridfire { void GraphEngine::setScreeningModel(const screening::ScreeningType model) { m_screeningModel = screening::selectScreeningModel(model); m_screeningType = model; + m_jacobianMatrixState = JacobianMatrixState::STALE; // The screening model affects the jacobian so if its changed the jacobian must be made stale } screening::ScreeningType GraphEngine::getScreeningModel() const { @@ -777,12 +788,24 @@ namespace gridfire { const double T9, const double rho ) const { + fourdst::composition::Composition mutableComp = comp; + for (const auto& species : m_networkSpecies) { + if (!comp.hasSpecies(species)) { + mutableComp.registerSpecies(species); + mutableComp.setMassFraction(species, 0.0); + } + } + const bool didFinalize = mutableComp.finalize(false); + if (!didFinalize) { + LOG_CRITICAL(m_logger, "Could not finalize the composition used to generate the jacobian matrix!"); + throw std::runtime_error("Could not finalize the composition used to generate the jacobian matrix"); + } LOG_TRACE_L1_LIMIT_EVERY_N(1000, 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 std::vector adInput(numSpecies + 2, 0.0); // +2 for T9 and rho - const std::vector& Y_dynamic = comp.getMolarAbundanceVector(); + const std::vector& Y_dynamic = mutableComp.getMolarAbundanceVector(); for (size_t i = 0; i < numSpecies; ++i) { adInput[i] = std::max(Y_dynamic[i], 1e-99); // regularize the jacobian... } @@ -794,16 +817,22 @@ namespace gridfire { // 3. Pack jacobian vector into sparse matrix m_jacobianMatrix.clear(); + // std::vector> columns; for (size_t i = 0; i < numSpecies; ++i) { + // std::vector colData; for (size_t j = 0; j < numSpecies; ++j) { const double value = dotY[i * (numSpecies + 2) + j]; if (std::abs(value) > MIN_JACOBIAN_THRESHOLD || i == j) { // Always keep diagonal elements to avoid pathological stiffness m_jacobianMatrix(i, j) = value; - } + // colData.push_back(value); } + // columns.push_back(std::make_unique>(std::to_string(i), colData)); } + // std::cout << utils::format_table("Jacobian after dense calculation", columns) << std::endl; + // exit(0); LOG_TRACE_L1_LIMIT_EVERY_N(1000, m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2()); + m_jacobianMatrixState = JacobianMatrixState::READY_DENSE; } void GraphEngine::generateJacobianMatrix( @@ -844,6 +873,7 @@ namespace gridfire { const double rho, const SparsityPattern &sparsityPattern ) const { + //TODO: The issue now seems to be that the jacobian is returning all zeros. I need to sort out why this is SparsityPattern intersectionSparsityPattern; for (const auto& entry : sparsityPattern) { if (m_full_sparsity_set.contains(entry)) { @@ -891,6 +921,12 @@ namespace gridfire { CppAD::sparse_rcv, std::vector> jac_subset(CppAD_sparsity_pattern); + // PERF: one of *the* most pressing things that needs to be done is remove the nead for this call every + // time the jacobian is needed since coloring is expensive and we are throwing away the caching + // power of CppAD by clearing the work vector each time. We do this since we make a new subset every + // time. However, a better solution would be to make the subset stateful so it only changes if the requested + // sparsity pattern changes. This way we could reuse the work vector. + m_jac_work.clear(); m_rhsADFun.sparse_jac_rev( x, jac_subset, // Sparse Jacobian output @@ -910,16 +946,35 @@ namespace gridfire { m_jacobianMatrix(row, col) = value; // Insert into the sparse matrix } } + m_jacobianMatrixState = JacobianMatrixState::READY_SPARSE; } double GraphEngine::getJacobianMatrixEntry( const fourdst::atomic::Species& rowSpecies, const fourdst::atomic::Species& colSpecies ) const { - //PERF: There may be some way to make this more efficient - const size_t i = getSpeciesIndex(rowSpecies); - const size_t j = getSpeciesIndex(colSpecies); - return m_jacobianMatrix(i, j); + switch (m_jacobianMatrixState) { + case JacobianMatrixState::STALE: { + const std::string staleMsg = std::format("Cannot retrieve jacobian entry for row {}, column {} as jacobian matrix is stale (has not been regenerated since last network topology change)", rowSpecies.name(), colSpecies.name()); + throw exceptions::StaleJacobianError(staleMsg); + } + case JacobianMatrixState::UNINITIALIZED: { + const std::string unInitMsg = std::format("Cannot retrieve jacobian entry for row {}, column {} as jacobian matrix is uninitialized (will return all 0s)", rowSpecies.name(), colSpecies.name()); + throw exceptions::UninitializedJacobianError(unInitMsg); + } + case JacobianMatrixState::READY_DENSE: + [[fallthrough]]; + case JacobianMatrixState::READY_SPARSE: { + const size_t i = getSpeciesIndex(rowSpecies); + const size_t j = getSpeciesIndex(colSpecies); + return m_jacobianMatrix(i, j); + } + default: { + // Code should not be able to get into this state + const std::string msg = std::format("An unknown error has occurred while attempting to retrieve the jacobian element at row {}, column {}. This should be taken as a catastrophic failure and reported to GridFire developers.", rowSpecies.name(), colSpecies.name()); + throw exceptions::UnknownJacobianError(msg); + } + } } std::unordered_map GraphEngine::getNetReactionStoichiometry( @@ -1365,4 +1420,46 @@ namespace gridfire { st[0] = rt[0]; return true; } + + bool GraphEngine::AtomicReverseRate::for_sparse_jac( + size_t q, + const CppAD::vector &r, + CppAD::vector &s, + const CppAD::vector &x + ) { + 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."); + + // S = R + for (size_t j = 0; j < q; j++) { + // s(0,j) = r(0,j) + s[j*m] = r[j*n]; + } + + return true; + } + + bool GraphEngine::AtomicReverseRate::rev_sparse_jac( + size_t q, + const CppAD::vector &rt, + CppAD::vector &st, + const CppAD::vector &x + ) { + 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."); + + // st = rt + for (size_t j = 0; j < q; j++) { + // st(j, 0) = rt(j, 0) + st[j * n] = rt[j * m]; + } + + return true; + } } diff --git a/src/lib/engine/procedures/construction.cpp b/src/lib/engine/procedures/construction.cpp index 1c90fcee..194ed57f 100644 --- a/src/lib/engine/procedures/construction.cpp +++ b/src/lib/engine/procedures/construction.cpp @@ -53,48 +53,49 @@ namespace gridfire { } } - 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 - ) - ); - } - } + // --- 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 // --- Step 2: Use non-owning raw pointers for the fast build algorithm --- std::vector remainingReactions; diff --git a/src/lib/reaction/reaction.cpp b/src/lib/reaction/reaction.cpp index 8dad84f8..ef89d495 100644 --- a/src/lib/reaction/reaction.cpp +++ b/src/lib/reaction/reaction.cpp @@ -2,6 +2,7 @@ #include #include +#include #include #include #include diff --git a/src/lib/solver/strategies/CVODE_solver_strategy.cpp b/src/lib/solver/strategies/CVODE_solver_strategy.cpp index 65b6639c..7b8cd839 100644 --- a/src/lib/solver/strategies/CVODE_solver_strategy.cpp +++ b/src/lib/solver/strategies/CVODE_solver_strategy.cpp @@ -160,6 +160,7 @@ namespace gridfire::solver { check_cvode_flag(m_cvode_mem == nullptr ? -1 : 0, "CVodeCreate"); initialize_cvode_integration_resources(N, numSpecies, 0.0, equilibratedComposition, absTol, relTol, 0.0); + m_engine.generateJacobianMatrix(equilibratedComposition, T9, netIn.density); CVODEUserData user_data; user_data.solver_instance = this; @@ -192,7 +193,7 @@ namespace gridfire::solver { std::rethrow_exception(std::make_exception_ptr(*user_data.captured_exception)); } - log_step_diagnostics(user_data); + // log_step_diagnostics(user_data, false); check_cvode_flag(flag, "CVode"); long int n_steps; @@ -290,9 +291,17 @@ namespace gridfire::solver { initialize_cvode_integration_resources(N, numSpecies, current_time, currentComposition, absTol, relTol, accumulated_energy); check_cvode_flag(CVodeReInit(m_cvode_mem, current_time, m_Y), "CVodeReInit"); + + LOG_TRACE_L1(m_logger, "Regenerating jacobian matrix..."); + m_engine.generateJacobianMatrix(currentComposition, T9, netIn.density); + LOG_TRACE_L1(m_logger, "Done regenerating jacobian matrix..."); } } + + // TODO: Need a more reliable way to get the final composition out, probably some methods that bubble it or something + // aside from that this now seems to be working + LOG_TRACE_L2(m_logger, "CVODE iteration complete"); sunrealtype* y_data = N_VGetArrayPointer(m_Y); @@ -313,6 +322,7 @@ namespace gridfire::solver { speciesNames.emplace_back(species.name()); } + LOG_TRACE_L2(m_logger, "Constructing final composition= with {} species", speciesNames.size()); fourdst::composition::Composition outputComposition(speciesNames); outputComposition.setMassFraction(speciesNames, finalMassFractions); bool didFinalize = outputComposition.finalize(true); @@ -320,6 +330,7 @@ namespace gridfire::solver { LOG_ERROR(m_logger, "Failed to finalize output composition after CVODE integration. Check output mass fractions for validity."); throw std::runtime_error("Failed to finalize output composition after CVODE integration."); } + LOG_TRACE_L2(m_logger, "Final composition constructed successfully!"); LOG_TRACE_L2(m_logger, "Constructing output data..."); NetOut netOut; @@ -327,11 +338,13 @@ namespace gridfire::solver { netOut.energy = accumulated_energy; check_cvode_flag(CVodeGetNumSteps(m_cvode_mem, reinterpret_cast(&netOut.num_steps)), "CVodeGetNumSteps"); + LOG_TRACE_L2(m_logger, "generating final nuclear energy generation rate derivatives..."); auto [dEps_dT, dEps_dRho] = m_engine.calculateEpsDerivatives( outputComposition, T9, netIn.density ); + LOG_TRACE_L2(m_logger, "Found dEps/dT: {:0.3E} and dEps/dRho: {:0.3E}", dEps_dT, dEps_dRho); netOut.dEps_dT = dEps_dT; netOut.dEps_dRho = dEps_dRho; @@ -396,16 +409,19 @@ namespace gridfire::solver { const long int N = SUNDenseMatrix_Columns(J); for (size_t j = 0; j < numSpecies; ++j) { + const fourdst::atomic::Species& species_j = engine->getNetworkSpecies()[j]; for (size_t i = 0; i < numSpecies; ++i) { - const fourdst::atomic::Species& species_j = engine->getNetworkSpecies()[j]; const fourdst::atomic::Species& species_i = engine->getNetworkSpecies()[i]; // J(i,j) = d(f_i)/d(y_j) - // Column-major order format for SUNDenseMatrix: J_data[j*N + i] - J_data[j * N + i] = engine->getJacobianMatrixEntry(species_i, species_j); + // Column-major order format for SUNDenseMatrix: J_data[j*N + i] indexes J(i,j) + const double dYi_dt = engine->getJacobianMatrixEntry(species_i, species_j); + J_data[j * N + i] = dYi_dt; } } // For now assume that the energy derivatives wrt. abundances are zero + // TODO: Need a better way to build this part of the output jacobian so it properly pushes the solver + // in the right direction. Currently we effectively are doing a fixed point iteration in energy space. for (size_t i = 0; i < N; ++i) { J_data[(N - 1) * N + i] = 0.0; // df(energy_dot)/df(y_i) J_data[i * N + (N - 1)] = 0.0; // df(f_i)/df(energy_dot) @@ -423,6 +439,15 @@ namespace gridfire::solver { const size_t numSpecies = m_engine.getNetworkSpecies().size(); sunrealtype* y_data = N_VGetArrayPointer(y); + // Solver constraints should keep these values very close to 0 but floating point noise can still result in very + // small negative numbers which can result in NaN's and more immediate crashes in the composition + // finalization stage + for (size_t i = 0; i < numSpecies; ++i) { + if (y_data[i] < 0.0) { + y_data[i] = 0.0; + } + } + // PERF: The trade off of ensured index consistency is some performance here. If this becomes a bottleneck we can revisit. // The specific trade off is that we have decided to enforce that all interfaces accept composition objects rather // than raw vectors of molar abundances. This then lets any method lookup the species by name rather than relying on @@ -539,7 +564,7 @@ namespace gridfire::solver { } } - void CVODESolverStrategy::log_step_diagnostics(const CVODEUserData &user_data) const { + void CVODESolverStrategy::log_step_diagnostics(const CVODEUserData &user_data, bool displayJacobianStiffness) const { check_cvode_flag(CVodeGetEstLocalErrors(m_cvode_mem, m_YErr), "CVodeGetEstLocalErrors"); sunrealtype *y_data = N_VGetArrayPointer(m_Y); @@ -616,9 +641,11 @@ namespace gridfire::solver { std::cout << utils::format_table("Species Error Ratios", columns) << std::endl; - diagnostics::inspect_jacobian_stiffness(*user_data.engine, composition, user_data.T9, user_data.rho); - for (const auto& species : sorted_speciesNames) { - diagnostics::inspect_species_balance(*user_data.engine, species, composition, user_data.T9, user_data.rho); + if (displayJacobianStiffness) { + diagnostics::inspect_jacobian_stiffness(*user_data.engine, composition, user_data.T9, user_data.rho); + for (const auto& species : sorted_speciesNames) { + diagnostics::inspect_species_balance(*user_data.engine, species, composition, user_data.T9, user_data.rho); + } } }