diff --git a/meson.build b/meson.build index e58bca8e..2dae91f6 100644 --- a/meson.build +++ b/meson.build @@ -79,9 +79,19 @@ subdir('build-config') subdir('src') -subdir('build-python') +if get_option('build-python') + message('Configuring Python bindings...') + subdir('src-pybind') +else + message('Skipping Python bindings...') +endif -subdir('tests') +if get_option('build-tests') + message('Setting up tests for GridFire...') + subdir('tests') +else + message('Skipping tests for GridFire...') +endif if get_option('pkg-config') diff --git a/meson_options.txt b/meson_options.txt index 82a8d3db..c42ba6fa 100644 --- a/meson_options.txt +++ b/meson_options.txt @@ -1,2 +1,4 @@ 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') -option('pkg-config', type: 'boolean', value: true, description: 'generate pkg-config file for GridFire (gridfire.pc)') \ No newline at end of file +option('pkg-config', type: 'boolean', value: true, description: 'generate pkg-config file for GridFire (gridfire.pc)') +option('build-python', type: 'boolean', value: true, description: 'build the python bindings so you can use GridFire from python') +option('build-tests', type: 'boolean', value: true, description: 'build the test suite') diff --git a/src/include/gridfire/engine/diagnostics/dynamic_engine_diagnostics.h b/src/include/gridfire/engine/diagnostics/dynamic_engine_diagnostics.h new file mode 100644 index 00000000..a65d9bd5 --- /dev/null +++ b/src/include/gridfire/engine/diagnostics/dynamic_engine_diagnostics.h @@ -0,0 +1,34 @@ +#pragma once + +#include "gridfire/engine/engine_abstract.h" + +#include +#include + +namespace gridfire::diagnostics { + void report_limiting_species( + const DynamicEngine& engine, + const std::vector& Y_full, + const std::vector& E_full, + const std::vector& dydt_full, + double relTol, + double absTol, + size_t top_n = 10 + ); + + void inspect_species_balance( + const DynamicEngine& engine, + const std::string& species_name, + const std::vector& Y_full, + double T9, + double rho + ); + + void inspect_jacobian_stiffness( + const DynamicEngine& engine, + const std::vector& Y_full, + double T9, + double rho + ); + +} \ No newline at end of file diff --git a/src/include/gridfire/engine/engine_abstract.h b/src/include/gridfire/engine/engine_abstract.h index e03efa9b..c9658426 100644 --- a/src/include/gridfire/engine/engine_abstract.h +++ b/src/include/gridfire/engine/engine_abstract.h @@ -64,6 +64,16 @@ namespace gridfire { using SparsityPattern = std::vector>; + struct EnergyDerivatives { + double dEps_dT = 0.0; ///< Partial derivative of energy generation rate with respect to temperature. + double dEps_dRho = 0.0;///< Partial derivative of energy generation rate with respect to density. + + friend std::ostream& operator<<(std::ostream& os, const EnergyDerivatives& ed) { + os << ""; + return os; + } + }; + /** * @brief Abstract base class for a reaction network engine. * @@ -210,6 +220,23 @@ namespace gridfire { double rho ) const = 0; + /** + * @brief Calculate the derivatives of the energy generation rate with respect to T and rho. + * + * @param Y Vector of current abundances. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @return EnergyDerivatives containing dEps/dT and dEps/dRho. + * + * This method computes the partial derivatives of the specific nuclear energy + * generation rate with respect to temperature and density for the current state. + */ + [[nodiscard]] virtual EnergyDerivatives calculateEpsDerivatives( + const std::vector& Y, + const double T9, + const double rho + ) const = 0; + /** * @brief Get the set of logical reactions in the network. * diff --git a/src/include/gridfire/engine/engine_graph.h b/src/include/gridfire/engine/engine_graph.h index b02736d4..0c918858 100644 --- a/src/include/gridfire/engine/engine_graph.h +++ b/src/include/gridfire/engine/engine_graph.h @@ -68,6 +68,7 @@ namespace gridfire { static constexpr double MIN_JACOBIAN_THRESHOLD = 1e-24; + /** * @class GraphEngine * @brief A reaction network engine that uses a graph-based representation. @@ -143,6 +144,12 @@ namespace gridfire { const double rho ) const override; + [[nodiscard]] EnergyDerivatives calculateEpsDerivatives( + const std::vector& Y, + const double T9, + const double rho + ) const override; + /** * @brief Generates the Jacobian matrix for the current state. * @@ -581,6 +588,7 @@ namespace gridfire { mutable boost::numeric::ublas::compressed_matrix m_jacobianMatrix; ///< Jacobian matrix (species x species). 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. CppAD::sparse_rc> m_full_jacobian_sparsity_pattern; ///< Full sparsity pattern for the Jacobian matrix. @@ -652,6 +660,8 @@ namespace gridfire { */ void recordADTape() const; + void recordEpsADTape() const; + void collectAtomicReverseRateAtomicBases(); void precomputeNetwork(); diff --git a/src/include/gridfire/engine/views/engine_adaptive.h b/src/include/gridfire/engine/views/engine_adaptive.h index a65dd59a..fd7bf4bd 100644 --- a/src/include/gridfire/engine/views/engine_adaptive.h +++ b/src/include/gridfire/engine/views/engine_adaptive.h @@ -107,6 +107,12 @@ namespace gridfire { const double rho ) const override; + [[nodiscard]] EnergyDerivatives calculateEpsDerivatives( + const std::vector &Y_culled, + const double T9, + const double rho + ) const override; + /** * @brief Generates the Jacobian matrix for the active species. * diff --git a/src/include/gridfire/engine/views/engine_defined.h b/src/include/gridfire/engine/views/engine_defined.h index c47ddab2..92de79c4 100644 --- a/src/include/gridfire/engine/views/engine_defined.h +++ b/src/include/gridfire/engine/views/engine_defined.h @@ -42,6 +42,13 @@ namespace gridfire{ const double T9, const double rho ) const override; + + [[nodiscard]] EnergyDerivatives calculateEpsDerivatives( + const std::vector &Y, + const double T9, + const double rho + ) const override; + /** * @brief Generates the Jacobian matrix for the active species. * diff --git a/src/include/gridfire/engine/views/engine_multiscale.h b/src/include/gridfire/engine/views/engine_multiscale.h index 4ca7b0dc..cf4ff91c 100644 --- a/src/include/gridfire/engine/views/engine_multiscale.h +++ b/src/include/gridfire/engine/views/engine_multiscale.h @@ -236,6 +236,12 @@ namespace gridfire { double rho ) const override; + [[nodiscard]] EnergyDerivatives calculateEpsDerivatives( + const std::vector &Y, + const double T9, + const double rho + ) const override; + /** * @brief Generates the Jacobian matrix for the current state. * @@ -754,6 +760,10 @@ namespace gridfire { * @return True if the sets of species indices are not identical. */ bool operator!=(const QSEGroup& other) const; + + std::string toString() const; + + std::string toString(DynamicEngine &engine) const; }; /** @@ -1056,9 +1066,11 @@ namespace gridfire { * flux exceeds a configurable threshold, the group is considered valid and is added * to the returned vector. */ - std::vector validateGroupsWithFluxAnalysis( + std::pair, std::vector> validateGroupsWithFluxAnalysis( const std::vector &candidate_groups, - const std::vector& Y, + const std::vector &Y, double T9, double rho ) const; diff --git a/src/include/gridfire/network.h b/src/include/gridfire/network.h index fb614aa8..036db3e2 100644 --- a/src/include/gridfire/network.h +++ b/src/include/gridfire/network.h @@ -57,16 +57,18 @@ namespace gridfire { double energy; ///< Energy in ergs double culling = 0.0; ///< Culling threshold for reactions (default is 0.0, meaning no culling) - std::vector MolarAbundance() const; + [[nodiscard]] std::vector MolarAbundance() const; }; struct NetOut { fourdst::composition::Composition composition; ///< Composition of the network after evaluation int num_steps; ///< Number of steps taken in the evaluation double energy; ///< Energy in ergs after evaluation + double dEps_dT; ///< Partial derivative of energy generation rate with respect to temperature + double dEps_dRho; ///< Partial derivative of energy generation rate with respect to density friend std::ostream& operator<<(std::ostream& os, const NetOut& netOut) { - os << "NetOut(composition=" << netOut.composition << ", num_steps=" << netOut.num_steps << ", energy=" << netOut.energy << ")"; + os << "NetOut(composition=" << netOut.composition << ", num_steps=" << netOut.num_steps << ", ε=" << netOut.energy << ", dε/dT=" << netOut.dEps_dT << ", dε/dρ=" << netOut.dEps_dRho << ")"; return os; } }; diff --git a/src/include/gridfire/partition/partition_abstract.h b/src/include/gridfire/partition/partition_abstract.h index 8230d86f..b64a8d29 100644 --- a/src/include/gridfire/partition/partition_abstract.h +++ b/src/include/gridfire/partition/partition_abstract.h @@ -4,7 +4,6 @@ #include namespace gridfire::partition { - /** * @class PartitionFunction * @brief Abstract interface for evaluating nuclear partition functions. diff --git a/src/include/gridfire/solver/strategies/CVODE_solver_strategy.h b/src/include/gridfire/solver/strategies/CVODE_solver_strategy.h index 4b2f7e92..6a05c48b 100644 --- a/src/include/gridfire/solver/strategies/CVODE_solver_strategy.h +++ b/src/include/gridfire/solver/strategies/CVODE_solver_strategy.h @@ -98,16 +98,18 @@ namespace gridfire::solver { DynamicEngine* engine; double T9; double rho; + double energy; const std::vector* networkSpecies; std::unique_ptr captured_exception = nullptr; }; private: Config& m_config = Config::getInstance(); + quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); static int cvode_rhs_wrapper(sunrealtype t, N_Vector y, N_Vector ydot, void *user_data); static int cvode_jac_wrapper(sunrealtype t, N_Vector y, N_Vector ydot, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); - int calculate_rhs(sunrealtype t, N_Vector y, N_Vector ydot, const CVODEUserData* data) const; + void calculate_rhs(sunrealtype t, N_Vector y, N_Vector ydot, const CVODEUserData* data) const; void initialize_cvode_integration_resources( uint64_t N, @@ -120,10 +122,13 @@ namespace gridfire::solver { ); void cleanup_cvode_resources(bool memFree); + + void log_step_diagnostics(const CVODEUserData& user_data) const; private: SUNContext m_sun_ctx = nullptr; void* m_cvode_mem = nullptr; N_Vector m_Y = nullptr; + N_Vector m_YErr = nullptr; SUNMatrix m_J = nullptr; SUNLinearSolver m_LS = nullptr; diff --git a/src/include/gridfire/utils/logging.h b/src/include/gridfire/utils/logging.h index c1b8d0b4..e08cf85a 100644 --- a/src/include/gridfire/utils/logging.h +++ b/src/include/gridfire/utils/logging.h @@ -4,6 +4,12 @@ #include #include +#include +#include +#include +#include +#include + namespace gridfire::utils { /** @@ -62,4 +68,6 @@ namespace gridfire::utils { const double T9, const double rho ); + + } \ No newline at end of file diff --git a/src/include/gridfire/utils/table_format.h b/src/include/gridfire/utils/table_format.h new file mode 100644 index 00000000..0e8d0b00 --- /dev/null +++ b/src/include/gridfire/utils/table_format.h @@ -0,0 +1,121 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + + +namespace gridfire::utils { + class ColumnBase { + public: + virtual ~ColumnBase() = default; + // Gets the string representation of the data at a given row + virtual std::string getCellData(size_t rowIndex) const = 0; + // Gets the header text for the column + virtual std::string getHeader() const = 0; + // Gets the number of data rows in the column + virtual size_t getRowCount() const = 0; + }; + + template + class Column final : public ColumnBase { + public: + Column(const std::string& header, const std::vector& data) + : m_header(header), m_data(data) {} + + std::string getCellData(size_t rowIndex) const override { + std::stringstream ss; + if (rowIndex < m_data.size()) { + ss << m_data[rowIndex]; + } + return ss.str(); + } + + std::string getHeader() const override { + return m_header; + } + + size_t getRowCount() const override { + return m_data.size(); + } + private: + std::string m_header; + std::vector m_data; + + }; + + inline std::string format_table(const std::string& tableName, const std::vector>& columns) { + // --- 1. Handle Empty Table --- + if (columns.empty()) { + return tableName + "\n(Table has no columns)\n"; + } + + // --- 2. Determine dimensions and calculate column widths --- + size_t num_cols = columns.size(); + size_t num_rows = 0; + for(const auto& col : columns) { + num_rows = std::max(num_rows, col->getRowCount()); + } + + std::vector col_widths(num_cols); + for (size_t j = 0; j < num_cols; ++j) { + col_widths[j] = columns[j]->getHeader().length(); + for (size_t i = 0; i < num_rows; ++i) { + col_widths[j] = std::max(col_widths[j], columns[j]->getCellData(i).length()); + } + } + + // --- 3. Build the table string using stringstream --- + std::stringstream table_ss; + + // --- Table Title --- + size_t total_width = std::accumulate(col_widths.begin(), col_widths.end(), 0) + (num_cols * 3) + 1; + size_t title_padding = (total_width > tableName.length()) ? (total_width - tableName.length()) / 2 : 0; + table_ss << std::string(title_padding, ' ') << tableName << "\n"; + + // --- Helper to draw horizontal border --- + auto draw_border = [&]() { + table_ss << "+"; + for (size_t width : col_widths) { + table_ss << std::string(width + 2, '-'); // +2 for padding + table_ss << "+"; + } + table_ss << "\n"; + }; + + // --- Draw Top Border --- + draw_border(); + + // --- Draw Header Row --- + table_ss << "|"; + for (size_t j = 0; j < num_cols; ++j) { + table_ss << " " << std::left << std::setw(col_widths[j]) << columns[j]->getHeader() << " |"; + } + table_ss << "\n"; + + // --- Draw Separator --- + draw_border(); + + // --- Draw Data Rows --- + for (size_t i = 0; i < num_rows; ++i) { + table_ss << "|"; + for (size_t j = 0; j < num_cols; ++j) { + table_ss << " " << std::left << std::setw(col_widths[j]) << columns[j]->getCellData(i) << " |"; + } + table_ss << "\n"; + } + + // --- Draw Bottom Border --- + draw_border(); + + return table_ss.str(); + } + + + +} \ No newline at end of file diff --git a/src/lib/engine/diagnostics/dynamic_engine_diagnostics.cpp b/src/lib/engine/diagnostics/dynamic_engine_diagnostics.cpp new file mode 100644 index 00000000..bbe434c7 --- /dev/null +++ b/src/lib/engine/diagnostics/dynamic_engine_diagnostics.cpp @@ -0,0 +1,168 @@ +#include "gridfire/engine/diagnostics/dynamic_engine_diagnostics.h" +#include "gridfire/engine/engine_abstract.h" +#include "gridfire/utils/table_format.h" + +#include +#include +#include +#include +#include + +namespace gridfire::diagnostics { + void report_limiting_species( + const DynamicEngine& engine, + const std::vector& Y_full, + const std::vector& E_full, + const std::vector& dydt_full, + const double relTol, + const double absTol, + const size_t top_n + ) { + struct SpeciesError { + std::string name; + double ratio; + double abundance; + double dydt; + }; + + const auto& species_list = engine.getNetworkSpecies(); + std::vector errors; + + for (size_t i = 0; i < species_list.size(); ++i) { + const double weight = relTol * std::abs(Y_full[i]) + absTol; + if (weight > 1e-99) { // Avoid division by zero for zero-abundance species + const double ratio = std::abs(E_full[i]) / weight; + errors.push_back({ + std::string(species_list[i].name()), + ratio, + Y_full[i], + dydt_full[i] + }); + } + } + + // Sort by error ratio in descending order + std::ranges::sort( + errors, + [](const auto& a, const auto& b) { + return a.ratio > b.ratio; + } + ); + + std::vector sorted_speciesNames; + std::vector sorted_err_ratios; + std::vector sorted_abundances; + std::vector sorted_dydt; + + for (size_t i = 0; i < std::min(top_n, errors.size()); ++i) { + sorted_speciesNames.push_back(errors[i].name); + sorted_err_ratios.push_back(errors[i].ratio); + sorted_abundances.push_back(errors[i].abundance); + sorted_dydt.push_back(errors[i].dydt); + } + + std::vector> columns; + columns.push_back(std::make_unique>("Species", sorted_speciesNames)); + columns.push_back(std::make_unique>("Error Ratio", sorted_err_ratios)); + columns.push_back(std::make_unique>("Abundance", sorted_abundances)); + columns.push_back(std::make_unique>("dY/dt", sorted_dydt)); + + std::cout << utils::format_table("Timestep Limiting Species", columns) << std::endl; + } + + void inspect_species_balance( + const DynamicEngine& engine, + const std::string& species_name, + const std::vector& Y_full, + const double T9, + const double rho + ) { + const auto& species_obj = fourdst::atomic::species.at(species_name); + + std::vector creation_ids, destruction_ids; + std::vector creation_stoichiometry, destruction_stoichiometry; + std::vector creation_flows, destruction_flows; + double total_creation_flow = 0.0; + double total_destruction_flow = 0.0; + + for (const auto& reaction : engine.getNetworkReactions()) { + const int stoichiometry = reaction->stoichiometry(species_obj); + if (stoichiometry == 0) continue; + + const double flow = engine.calculateMolarReactionFlow(*reaction, Y_full, T9, rho); + + if (stoichiometry > 0) { + creation_ids.push_back(std::string(reaction->id())); + creation_stoichiometry.push_back(stoichiometry); + creation_flows.push_back(flow); + total_creation_flow += stoichiometry * flow; + } else { + destruction_ids.push_back(std::string(reaction->id())); + destruction_stoichiometry.push_back(stoichiometry); + destruction_flows.push_back(flow); + total_destruction_flow += std::abs(stoichiometry) * flow; + } + } + + { + std::vector> columns; + columns.push_back(std::make_unique>("Reaction ID", creation_ids)); + columns.push_back(std::make_unique>("Stoichiometry", creation_stoichiometry)); + columns.push_back(std::make_unique>("Molar Flow", creation_flows)); + std::cout << utils::format_table("Creation Reactions for " + species_name, columns) << std::endl; + } + + { + std::vector> columns; + columns.push_back(std::make_unique>("Reaction ID", destruction_ids)); + columns.push_back(std::make_unique>("Stoichiometry", destruction_stoichiometry)); + columns.push_back(std::make_unique>("Molar Flow", destruction_flows)); + std::cout << utils::format_table("Destruction Reactions for " + species_name, columns) << std::endl; + } + + std::cout << "--- Balance Summary for " << species_name << " ---" << std::endl; + std::cout << " Total Creation Rate: " << std::scientific << total_creation_flow << " [mol/g/s]" << std::endl; + std::cout << " Total Destruction Rate: " << std::scientific << total_destruction_flow << " [mol/g/s]" << std::endl; + std::cout << " Net dY/dt: " << std::scientific << (total_creation_flow - total_destruction_flow) << std::endl; + std::cout << "-----------------------------------" << std::endl; + } + + void inspect_jacobian_stiffness( + const DynamicEngine& engine, + const std::vector& Y_full, + const double T9, + const double rho + ) { + engine.generateJacobianMatrix(Y_full, T9, rho); + const auto& species_list = engine.getNetworkSpecies(); + + double max_diag = 0.0; + double max_off_diag = 0.0; + int max_diag_idx = -1; + int max_off_diag_i = -1, max_off_diag_j = -1; + + for (size_t i = 0; i < species_list.size(); ++i) { + for (size_t j = 0; j < species_list.size(); ++j) { + const double val = std::abs(engine.getJacobianMatrixEntry(i, j)); + if (i == j) { + if (val > max_diag) { max_diag = val; max_diag_idx = i; } + } else { + if (val > max_off_diag) { max_off_diag = val; max_off_diag_i = i; max_off_diag_j = j; } + } + } + } + + std::cout << "\n--- Jacobian Stiffness Report ---" << std::endl; + if (max_diag_idx != -1) { + std::cout << " Largest Diagonal Element (d(dYi/dt)/dYi): " << std::scientific << max_diag + << " for species " << species_list[max_diag_idx].name() << std::endl; + } + if (max_off_diag_i != -1) { + std::cout << " Largest Off-Diagonal Element (d(dYi/dt)/dYj): " << std::scientific << max_off_diag + << " for d(" << species_list[max_off_diag_i].name() + << ")/d(" << species_list[max_off_diag_j].name() << ")" << std::endl; + } + std::cout << "---------------------------------" << std::endl; + } + +} \ No newline at end of file diff --git a/src/lib/engine/engine_graph.cpp b/src/lib/engine/engine_graph.cpp index c7c2c42f..352586db 100644 --- a/src/lib/engine/engine_graph.cpp +++ b/src/lib/engine/engine_graph.cpp @@ -78,6 +78,49 @@ namespace gridfire { } } + EnergyDerivatives GraphEngine::calculateEpsDerivatives( + const std::vector &Y, + const double T9, + const double rho + ) const { + const size_t numSpecies = m_networkSpecies.size(); + const size_t numADInputs = numSpecies + 2; // +2 for T9 and rho + + if (Y.size() != numSpecies) { + LOG_ERROR(m_logger, "Input abundance vector size ({}) does not match number of species in the network ({}).", + Y.size(), numSpecies); + throw std::invalid_argument("Input abundance vector size does not match number of species in the network."); + } + + std::vector x(numADInputs); + for (size_t i = 0; i < numSpecies; ++i) { + x[i] = Y[i]; + } + + x[numSpecies] = T9; + x[numSpecies + 1] = rho; + + // Use reverse mode to get the gradient. W selects which dependent variable we care about, the Eps AD tape only has eps as a dependent variable so we just select set the 0th element to 1. + std::vector w(1); + w[0] = 1.0; // We want the derivative of the energy generation rate + + // Sweep the tape forward to record the function value at x + m_epsADFun.Forward(0, x); + + // Extract the gradient at the previously evaluated point x using reverse mode + const std::vector eps_derivatives = m_epsADFun.Reverse(1, w); + + const double dEps_dT9 = eps_derivatives[numSpecies]; + const double dEps_dRho = eps_derivatives[numSpecies + 1]; + + // Chain rule to scale from deps/dT9 to deps/dT + // dT9/dT = 1e-9 + const double dEps_dT = dEps_dT9 * 1e-9; + + + return {dEps_dT, dEps_dRho}; + } + void GraphEngine::syncInternalMaps() { LOG_INFO(m_logger, "Synchronizing internal maps for REACLIB graph network (serif::network::GraphNetwork)..."); collectNetworkSpecies(); @@ -86,7 +129,8 @@ namespace gridfire { collectAtomicReverseRateAtomicBases(); generateStoichiometryMatrix(); reserveJacobianMatrix(); - recordADTape(); + recordADTape(); // Record the AD tape for the RHS function + recordEpsADTape(); // Record the AD tape for the energy generation rate function const size_t n = m_rhsADFun.Domain(); const size_t m = m_rhsADFun.Range(); @@ -161,13 +205,13 @@ namespace gridfire { // --- Basic Accessors and Queries --- const std::vector& GraphEngine::getNetworkSpecies() const { // Returns a constant reference to the vector of unique species in the network. - LOG_TRACE_L3(m_logger, "Providing access to network species vector. Size: {}.", m_networkSpecies.size()); + // LOG_TRACE_L3(m_logger, "Providing access to network species vector. Size: {}.", m_networkSpecies.size()); return m_networkSpecies; } const reaction::ReactionSet& GraphEngine::getNetworkReactions() const { // Returns a constant reference to the set of reactions in the network. - LOG_TRACE_L3(m_logger, "Providing access to network reactions set. Size: {}.", m_reactions.size()); + // LOG_TRACE_L3(m_logger, "Providing access to network reactions set. Size: {}.", m_reactions.size()); return m_reactions; } @@ -931,6 +975,43 @@ namespace gridfire { adInput.size()); } + void GraphEngine::recordEpsADTape() const { + LOG_TRACE_L1(m_logger, "Recording AD tape for the EPS calculation..."); + + const size_t numSpecies = m_networkSpecies.size(); + if (numSpecies == 0) { + LOG_ERROR(m_logger, "Cannot record EPS tape: No species in the network."); + throw std::runtime_error("Cannot record EPS AD tape: No species in the network."); + } + + const size_t numADInputs = numSpecies + 2; // Y + T9 + rho + + std::vector> adInput(numADInputs, 0.0); + + for (size_t i = 0; i < numSpecies; ++i) { + adInput[i] = static_cast>(1.0 / static_cast(numSpecies)); // Uniform distribution + } + adInput[numSpecies] = 1.0; // Dummy T9 + adInput[numSpecies + 1] = 1.0; // Dummy rho + + CppAD::Independent(adInput); + + const std::vector> adY(adInput.begin(), adInput.begin() + numSpecies); + const CppAD::AD adT9 = adInput[numSpecies]; + const CppAD::AD adRho = adInput[numSpecies + 1]; + + StepDerivatives> derivatives = calculateAllDerivatives>(adY, adT9, adRho); + + std::vector> adOutput(1); + adOutput[0] = derivatives.nuclearEnergyGenerationRate; + m_epsADFun.Dependent(adInput, adOutput); + // m_epsADFun.optimize(); + + LOG_TRACE_L1(m_logger, "AD tape recorded successfully for the EPS calculation. Number of independent variables: {}.", adInput.size()); + + + } + void GraphEngine::collectAtomicReverseRateAtomicBases() { m_atomicReverseRates.clear(); m_atomicReverseRates.reserve(m_reactions.size()); diff --git a/src/lib/engine/procedures/priming.cpp b/src/lib/engine/procedures/priming.cpp index 7f3bfd79..8e8600fc 100644 --- a/src/lib/engine/procedures/priming.cpp +++ b/src/lib/engine/procedures/priming.cpp @@ -35,10 +35,42 @@ namespace gridfire { return dominateReaction; } - + /** + * @brief Primes absent species in the network to their equilibrium abundances using a robust, two-stage approach. + * + * @details This function implements a robust network priming algorithm that avoids the pitfalls of + * sequential, one-by-one priming. The previous, brittle method could allow an early priming + * reaction to consume all of a shared reactant, starving later reactions. This new, two-stage + * method ensures that all priming reactions are considered collectively, competing for the + * same limited pool of initial reactants in a physically consistent manner. + * + * The algorithm proceeds in three main stages: + * 1. **Calculation Stage:** It first loops through all species that need priming. For each one, + * it calculates its theoretical equilibrium mass fraction and identifies the dominant + * creation channel. Crucially, it *does not* modify any abundances at this stage. Instead, + * it stores these calculations as a list of "mass transfer requests". + * + * 2. **Collective Scaling Stage:** It then processes the full list of requests to determine the + * total "debit" required from each reactant. By comparing these total debits to the + * initially available mass of each reactant, it calculates a single, global `scalingFactor`. + * If any reactant is overdrawn, this factor will be less than 1.0, ensuring that no + * reactant's abundance can go negative. + * + * 3. **Application Stage:** Finally, it loops through the requests again, applying the mass + * transfers. Each calculated equilibrium mass fraction and corresponding reactant debit is + * multiplied by the global `scalingFactor` before being applied to the final composition. + * This ensures that if resources are limited, all primed species are scaled down proportionally. + * + * @param netIn Input network data containing initial composition, temperature, and density. + * @param engine DynamicEngine used to build and evaluate the reaction network. + * @return PrimingReport encapsulating the results of the priming operation, including the new + * robustly primed composition. + */ PrimingReport primeNetwork(const NetIn& netIn, DynamicEngine& engine) { auto logger = LogManager::getInstance().getLogger("log"); + // --- Initial Setup --- + // Identify all species with zero initial mass fraction that need to be primed. std::vector speciesToPrime; for (const auto &entry: netIn.composition | std::views::values) { if (entry.mass_fraction() == 0.0) { @@ -47,6 +79,7 @@ namespace gridfire { } LOG_DEBUG(logger, "Priming {} species in the network.", speciesToPrime.size()); + // If no species need priming, return immediately. PrimingReport report; if (speciesToPrime.empty()) { report.primedComposition = netIn.composition; @@ -58,43 +91,46 @@ namespace gridfire { const double T9 = netIn.temperature / 1e9; const double rho = netIn.density; const auto initialReactionSet = engine.getNetworkReactions(); - report.status = PrimingReportStatus::FULL_SUCCESS; report.success = true; - // --- 1: pack composition into internal map --- + // Create a mutable map of the mass fractions that we will modify. std::unordered_map currentMassFractions; for (const auto& entry : netIn.composition | std::views::values) { currentMassFractions[entry.isotope()] = entry.mass_fraction(); } + // Ensure all species to be primed exist in the map, initialized to zero. for (const auto& entry : speciesToPrime) { - currentMassFractions[entry] = 0.0; // Initialize priming species with 0 mass fraction + currentMassFractions[entry] = 0.0; } - std::unordered_map totalMassFractionChanges; - + // Rebuild the engine with the full network to ensure all possible creation channels are available. engine.rebuild(netIn.composition, NetworkBuildDepth::Full); - for (const auto& primingSpecies : speciesToPrime) { - LOG_TRACE_L3(logger, "Priming species: {}", primingSpecies.name()); + // --- STAGE 1: Calculation and Bookkeeping (No Modifications) --- + // In this stage, we calculate all required mass transfers but do not apply them yet. - // Create a temporary composition from the current internal state for the primer + // A struct to hold the result of each individual priming calculation. + struct MassTransferRequest { + Species species_to_prime; + double equilibrium_mass_fraction; + std::vector reactants; + }; + std::vector requests; + + for (const auto& primingSpecies : speciesToPrime) { + // Create a temporary composition reflecting the current state for rate calculations. Composition tempComp; for(const auto& [sp, mf] : currentMassFractions) { tempComp.registerSymbol(std::string(sp.name())); - if (mf < 0.0 && std::abs(mf) < 1e-16) { - tempComp.setMassFraction(sp, 0.0); // Avoid negative mass fractions - } else { - tempComp.setMassFraction(sp, mf); - } + tempComp.setMassFraction(sp, std::max(0.0, mf)); } - tempComp.finalize(true); // Finalize with normalization + tempComp.finalize(true); NetIn tempNetIn = netIn; tempNetIn.composition = tempComp; NetworkPrimingEngineView primer(primingSpecies, engine); - if (primer.getNetworkReactions().size() == 0) { LOG_ERROR(logger, "No priming reactions found for species {}.", primingSpecies.name()); report.success = false; @@ -106,60 +142,87 @@ namespace gridfire { const double destructionRateConstant = calculateDestructionRateConstant(primer, primingSpecies, Y, T9, rho); if (destructionRateConstant > 1e-99) { - double equilibriumMassFraction = 0.0; const double creationRate = calculateCreationRate(primer, primingSpecies, Y, T9, rho); - equilibriumMassFraction = (creationRate / destructionRateConstant) * primingSpecies.mass(); - if (std::isnan(equilibriumMassFraction)) { - LOG_WARNING(logger, "Equilibrium mass fraction for {} is NaN. Setting to 0.0. This is likely not an issue. It probably originates from all reactions leading to creation and destruction being frozen out. In that case 0.0 should be a good approximation. Hint: This happens often when the network temperature is very the low. ", primingSpecies.name()); - equilibriumMassFraction = 0.0; - } - LOG_TRACE_L3(logger, "Found equilibrium for {}: X_eq = {:.4e}", primingSpecies.name(), equilibriumMassFraction); + double equilibriumMassFraction = (creationRate / destructionRateConstant) * primingSpecies.mass(); + if (std::isnan(equilibriumMassFraction)) equilibriumMassFraction = 0.0; if (const reaction::Reaction* dominantChannel = findDominantCreationChannel(primer, primingSpecies, Y, T9, rho)) { - LOG_TRACE_L3(logger, "Dominant creation channel for {}: {}", primingSpecies.name(), dominantChannel->id()); - - double totalReactantMass = 0.0; - for (const auto& reactant : dominantChannel->reactants()) { - totalReactantMass += reactant.mass(); - } - - double scalingFactor = 1.0; - for (const auto& reactant : dominantChannel->reactants()) { - const double massToSubtract = equilibriumMassFraction * (reactant.mass() / totalReactantMass); - double availableMass = 0.0; - if (currentMassFractions.contains(reactant)) { - availableMass = currentMassFractions.at(reactant); - } - if (massToSubtract > availableMass && availableMass > 0) { - scalingFactor = std::min(scalingFactor, availableMass / massToSubtract); - } - } - - if (scalingFactor < 1.0) { - LOG_WARNING(logger, "Priming for {} was limited by reactant availability. Scaling transfer by {:.4e}", primingSpecies.name(), scalingFactor); - equilibriumMassFraction *= scalingFactor; - } - - // Update the internal mass fraction map and accumulate total changes - totalMassFractionChanges[primingSpecies] += equilibriumMassFraction; - currentMassFractions[primingSpecies] += equilibriumMassFraction; - - for (const auto& reactant : dominantChannel->reactants()) { - const double massToSubtract = equilibriumMassFraction * (reactant.mass() / totalReactantMass); - totalMassFractionChanges[reactant] -= massToSubtract; - currentMassFractions[reactant] -= massToSubtract; - } + // Store the request instead of applying it immediately. + requests.push_back({primingSpecies, equilibriumMassFraction, dominantChannel->reactants()}); } else { - LOG_ERROR(logger, "Failed to find dominant creation channel for {}.", primingSpecies.name()); - report.status = PrimingReportStatus::FAILED_TO_FIND_CREATION_CHANNEL; - totalMassFractionChanges[primingSpecies] += 1e-40; - currentMassFractions[primingSpecies] += 1e-40; + LOG_ERROR(logger, "Failed to find dominant creation channel for {}.", primingSpecies.name()); + report.status = PrimingReportStatus::FAILED_TO_FIND_CREATION_CHANNEL; } } else { LOG_WARNING(logger, "No destruction channel found for {}. Using fallback abundance.", primingSpecies.name()); - totalMassFractionChanges[primingSpecies] += 1e-40; - currentMassFractions[primingSpecies] += 1e-40; - report.status = PrimingReportStatus::BASE_NETWORK_TOO_SHALLOW; + // For species with no destruction, we can't calculate an equilibrium. + // We add a request with a tiny fallback abundance to ensure it exists in the network. + requests.push_back({primingSpecies, 1e-40, {}}); + } + } + + // --- STAGE 2: Collective Scaling Based on Reactant Availability --- + // Now, we determine the total demand for each reactant and find a global scaling factor. + + std::unordered_map total_mass_debits; + for (const auto& req : requests) { + if (req.reactants.empty()) continue; // Skip fallbacks which don't consume reactants. + + double totalReactantMass = 0.0; + for (const auto& reactant : req.reactants) { + totalReactantMass += reactant.mass(); + } + if (totalReactantMass == 0.0) continue; + + for (const auto& reactant : req.reactants) { + const double massToSubtract = req.equilibrium_mass_fraction * (reactant.mass() / totalReactantMass); + total_mass_debits[reactant] += massToSubtract; + } + } + + double globalScalingFactor = 1.0; + for (const auto& [reactant, total_debit] : total_mass_debits) { + double availableMass; + if (currentMassFractions.contains(reactant)) { + availableMass = currentMassFractions.at(reactant); + } else { + availableMass = 0.0; + } + if (total_debit > availableMass && availableMass > 0) { + globalScalingFactor = std::min(globalScalingFactor, availableMass / total_debit); + } + } + + if (globalScalingFactor < 1.0) { + LOG_WARNING(logger, "Priming was limited by reactant availability. All transfers will be scaled by {:.4e}", globalScalingFactor); + } + + // --- STAGE 3: Application of Scaled Mass Transfers --- + // Finally, apply all the transfers, scaled by our global factor. + + std::unordered_map totalMassFractionChanges; + for (const auto&[species_to_prime, equilibrium_mass_fraction, reactants] : requests) { + const double scaled_equilibrium_mf = equilibrium_mass_fraction * globalScalingFactor; + + // Add the scaled mass to the primed species. + currentMassFractions.at(species_to_prime) += scaled_equilibrium_mf; + totalMassFractionChanges[species_to_prime] += scaled_equilibrium_mf; + + // Subtract the scaled mass from the reactants. + if (!reactants.empty()) { + double totalReactantMass = 0.0; + for (const auto& reactant : reactants) { + totalReactantMass += reactant.mass(); + } + if (totalReactantMass == 0.0) continue; + + for (const auto& reactant : reactants) { + const double massToSubtract = scaled_equilibrium_mf * (reactant.mass() / totalReactantMass); + if (massToSubtract != 0) { + currentMassFractions.at(reactant) -= massToSubtract; + totalMassFractionChanges[reactant] -= massToSubtract; + } + } } } @@ -168,14 +231,9 @@ namespace gridfire { std::vector final_mass_fractions; for(const auto& [species, mass_fraction] : currentMassFractions) { final_symbols.emplace_back(species.name()); - if (mass_fraction < 0.0 && std::abs(mass_fraction) < 1e-16) { - final_mass_fractions.push_back(0.0); // Avoid negative mass fractions - } else { - final_mass_fractions.push_back(mass_fraction); - } + final_mass_fractions.push_back(std::max(0.0, mass_fraction)); // Ensure no negative mass fractions. } - // Create the final composition object from the pre-normalized mass fractions Composition primedComposition(final_symbols, final_mass_fractions, true); report.primedComposition = primedComposition; @@ -183,10 +241,166 @@ namespace gridfire { report.massFractionChanges.emplace_back(species, change); } + // Restore the engine to its original, smaller network state. engine.setNetworkReactions(initialReactionSet); return report; } + + + // PrimingReport primeNetwork(const NetIn& netIn, DynamicEngine& engine) { + // auto logger = LogManager::getInstance().getLogger("log"); + // + // std::vector speciesToPrime; + // for (const auto &entry: netIn.composition | std::views::values) { + // std::cout << "Checking species: " << entry.isotope().name() << " with mass fraction: " << entry.mass_fraction() << std::endl; + // if (entry.mass_fraction() == 0.0) { + // speciesToPrime.push_back(entry.isotope()); + // } + // } + // LOG_DEBUG(logger, "Priming {} species in the network.", speciesToPrime.size()); + // + // PrimingReport report; + // if (speciesToPrime.empty()) { + // report.primedComposition = netIn.composition; + // report.success = true; + // report.status = PrimingReportStatus::NO_SPECIES_TO_PRIME; + // return report; + // } + // + // const double T9 = netIn.temperature / 1e9; + // const double rho = netIn.density; + // const auto initialReactionSet = engine.getNetworkReactions(); + // + // report.status = PrimingReportStatus::FULL_SUCCESS; + // report.success = true; + // + // // --- 1: pack composition into internal map --- + // std::unordered_map currentMassFractions; + // for (const auto& entry : netIn.composition | std::views::values) { + // currentMassFractions[entry.isotope()] = entry.mass_fraction(); + // } + // for (const auto& entry : speciesToPrime) { + // currentMassFractions[entry] = 0.0; // Initialize priming species with 0 mass fraction + // } + // + // std::unordered_map totalMassFractionChanges; + // + // engine.rebuild(netIn.composition, NetworkBuildDepth::Full); + // + // for (const auto& primingSpecies : speciesToPrime) { + // LOG_TRACE_L3(logger, "Priming species: {}", primingSpecies.name()); + // + // // Create a temporary composition from the current internal state for the primer + // Composition tempComp; + // for(const auto& [sp, mf] : currentMassFractions) { + // tempComp.registerSymbol(std::string(sp.name())); + // if (mf < 0.0 && std::abs(mf) < 1e-16) { + // tempComp.setMassFraction(sp, 0.0); // Avoid negative mass fractions + // } else { + // tempComp.setMassFraction(sp, mf); + // } + // } + // tempComp.finalize(true); // Finalize with normalization + // + // NetIn tempNetIn = netIn; + // tempNetIn.composition = tempComp; + // + // NetworkPrimingEngineView primer(primingSpecies, engine); + // + // if (primer.getNetworkReactions().size() == 0) { + // LOG_ERROR(logger, "No priming reactions found for species {}.", primingSpecies.name()); + // report.success = false; + // report.status = PrimingReportStatus::FAILED_TO_FIND_PRIMING_REACTIONS; + // continue; + // } + // + // const auto Y = primer.mapNetInToMolarAbundanceVector(tempNetIn); + // const double destructionRateConstant = calculateDestructionRateConstant(primer, primingSpecies, Y, T9, rho); + // + // if (destructionRateConstant > 1e-99) { + // double equilibriumMassFraction = 0.0; + // const double creationRate = calculateCreationRate(primer, primingSpecies, Y, T9, rho); + // equilibriumMassFraction = (creationRate / destructionRateConstant) * primingSpecies.mass(); + // if (std::isnan(equilibriumMassFraction)) { + // LOG_WARNING(logger, "Equilibrium mass fraction for {} is NaN. Setting to 0.0. This is likely not an issue. It probably originates from all reactions leading to creation and destruction being frozen out. In that case 0.0 should be a good approximation. Hint: This happens often when the network temperature is very the low. ", primingSpecies.name()); + // equilibriumMassFraction = 0.0; + // } + // LOG_TRACE_L3(logger, "Found equilibrium for {}: X_eq = {:.4e}", primingSpecies.name(), equilibriumMassFraction); + // + // if (const reaction::Reaction* dominantChannel = findDominantCreationChannel(primer, primingSpecies, Y, T9, rho)) { + // LOG_TRACE_L3(logger, "Dominant creation channel for {}: {}", primingSpecies.name(), dominantChannel->id()); + // + // double totalReactantMass = 0.0; + // for (const auto& reactant : dominantChannel->reactants()) { + // totalReactantMass += reactant.mass(); + // } + // + // double scalingFactor = 1.0; + // for (const auto& reactant : dominantChannel->reactants()) { + // const double massToSubtract = equilibriumMassFraction * (reactant.mass() / totalReactantMass); + // double availableMass = 0.0; + // if (currentMassFractions.contains(reactant)) { + // availableMass = currentMassFractions.at(reactant); + // } + // if (massToSubtract > availableMass && availableMass > 0) { + // scalingFactor = std::min(scalingFactor, availableMass / massToSubtract); + // } + // } + // + // if (scalingFactor < 1.0) { + // LOG_WARNING(logger, "Priming for {} was limited by reactant availability. Scaling transfer by {:.4e}", primingSpecies.name(), scalingFactor); + // equilibriumMassFraction *= scalingFactor; + // } + // + // // Update the internal mass fraction map and accumulate total changes + // totalMassFractionChanges[primingSpecies] += equilibriumMassFraction; + // currentMassFractions.at(primingSpecies) += equilibriumMassFraction; + // + // for (const auto& reactant : dominantChannel->reactants()) { + // const double massToSubtract = equilibriumMassFraction * (reactant.mass() / totalReactantMass); + // std::cout << "[Priming: " << primingSpecies.name() << ", Channel: " << dominantChannel->id() << "] Subtracting " << massToSubtract << " from reactant " << reactant.name() << std::endl; + // totalMassFractionChanges[reactant] -= massToSubtract; + // currentMassFractions[reactant] -= massToSubtract; + // } + // } else { + // LOG_ERROR(logger, "Failed to find dominant creation channel for {}.", primingSpecies.name()); + // report.status = PrimingReportStatus::FAILED_TO_FIND_CREATION_CHANNEL; + // totalMassFractionChanges[primingSpecies] += 1e-40; + // currentMassFractions.at(primingSpecies) += 1e-40; + // } + // } else { + // LOG_WARNING(logger, "No destruction channel found for {}. Using fallback abundance.", primingSpecies.name()); + // totalMassFractionChanges.at(primingSpecies) += 1e-40; + // currentMassFractions.at(primingSpecies) += 1e-40; + // report.status = PrimingReportStatus::BASE_NETWORK_TOO_SHALLOW; + // } + // } + // + // // --- Final Composition Construction --- + // std::vector final_symbols; + // std::vector final_mass_fractions; + // for(const auto& [species, mass_fraction] : currentMassFractions) { + // final_symbols.emplace_back(species.name()); + // if (mass_fraction < 0.0 && std::abs(mass_fraction) < 1e-16) { + // final_mass_fractions.push_back(0.0); // Avoid negative mass fractions + // } else { + // final_mass_fractions.push_back(mass_fraction); + // } + // } + // + // // Create the final composition object from the pre-normalized mass fractions + // Composition primedComposition(final_symbols, final_mass_fractions, true); + // + // report.primedComposition = primedComposition; + // for (const auto& [species, change] : totalMassFractionChanges) { + // report.massFractionChanges.emplace_back(species, change); + // } + // + // engine.setNetworkReactions(initialReactionSet); + // return report; + // } + double calculateDestructionRateConstant( const DynamicEngine& engine, const fourdst::atomic::Species& species, diff --git a/src/lib/engine/views/engine_adaptive.cpp b/src/lib/engine/views/engine_adaptive.cpp index 34d830ca..7b906088 100644 --- a/src/lib/engine/views/engine_adaptive.cpp +++ b/src/lib/engine/views/engine_adaptive.cpp @@ -167,6 +167,17 @@ namespace gridfire { return culledResults; } + EnergyDerivatives AdaptiveEngineView::calculateEpsDerivatives( + const std::vector &Y_culled, + const double T9, + const double rho + ) const { + validateState(); + const auto Y_full = mapCulledToFull(Y_culled); + + return m_baseEngine.calculateEpsDerivatives(Y_full, T9, rho); + } + void AdaptiveEngineView::generateJacobianMatrix( const std::vector &Y_dynamic, const double T9, diff --git a/src/lib/engine/views/engine_defined.cpp b/src/lib/engine/views/engine_defined.cpp index c05d36a9..5c30d628 100644 --- a/src/lib/engine/views/engine_defined.cpp +++ b/src/lib/engine/views/engine_defined.cpp @@ -48,6 +48,17 @@ namespace gridfire { return definedResults; } + EnergyDerivatives DefinedEngineView::calculateEpsDerivatives( + const std::vector &Y_dynamic, + const double T9, + const double rho + ) const { + validateNetworkState(); + + const auto Y_full = mapViewToFull(Y_dynamic); + return m_baseEngine.calculateEpsDerivatives(Y_full, T9, rho); + } + void DefinedEngineView::generateJacobianMatrix( const std::vector &Y_dynamic, const double T9, diff --git a/src/lib/engine/views/engine_multiscale.cpp b/src/lib/engine/views/engine_multiscale.cpp index 13309ffe..0c5e5113 100644 --- a/src/lib/engine/views/engine_multiscale.cpp +++ b/src/lib/engine/views/engine_multiscale.cpp @@ -1,5 +1,6 @@ #include "gridfire/engine/views/engine_multiscale.h" #include "gridfire/exceptions/error_engine.h" +#include "gridfire/engine/procedures/priming.h" #include #include @@ -11,6 +12,7 @@ #include #include +#include "gridfire/utils/logging.h" #include "quill/LogMacros.h" #include "quill/Logger.h" @@ -200,6 +202,14 @@ namespace gridfire { return deriv; } + EnergyDerivatives MultiscalePartitioningEngineView::calculateEpsDerivatives( + const std::vector &Y, + const double T9, + const double rho + ) const { + return m_baseEngine.calculateEpsDerivatives(Y, T9, rho); + } + void MultiscalePartitioningEngineView::generateJacobianMatrix( const std::vector &Y_full, const double T9, @@ -228,12 +238,12 @@ namespace gridfire { const int i_full, const int j_full ) const { - // // Check if the species we are differentiating with respect to is algebraic or dynamic. If it is algebraic we can reduce the work significantly... - // if (std::ranges::contains(m_algebraic_species_indices, j_full)) { - // const auto& species = m_baseEngine.getNetworkSpecies()[j_full]; - // // If j is algebraic, we can return 0.0 since the Jacobian entry for algebraic species is always zero. - // return 0.0; - // } + // Check if the species we are differentiating with respect to is algebraic or dynamic. If it is algebraic we can reduce the work significantly... + if (std::ranges::contains(m_algebraic_species_indices, j_full)) { + // const auto& species = m_baseEngine.getNetworkSpecies()[j_full]; + // If j is algebraic, we can return 0.0 since the Jacobian entry for algebraic species is always zero. + return 0.0; + } // Otherwise we need to query the full jacobian return m_baseEngine.getJacobianMatrixEntry(i_full, j_full); } @@ -481,9 +491,42 @@ namespace gridfire { LOG_TRACE_L1(m_logger, "Identifying potential seed species for candidate pools..."); const std::vector candidate_groups = constructCandidateGroups(connected_pools, Y, T9, rho); LOG_TRACE_L1(m_logger, "Found {} candidate QSE groups for further analysis", candidate_groups.size()); + LOG_TRACE_L2( + m_logger, + "{}", + [&]() -> std::string { + std::stringstream ss; + int j = 0; + for (const auto& group : candidate_groups) { + ss << "CandidateQSEGroup(Algebraic: {"; + int i = 0; + for (const auto& index : group.algebraic_indices) { + ss << m_baseEngine.getNetworkSpecies()[index].name(); + if (i < group.algebraic_indices.size() - 1) { + ss << ", "; + } + } + ss << "}, Seed: {"; + i = 0; + for (const auto& index : group.seed_indices) { + ss << m_baseEngine.getNetworkSpecies()[index].name(); + if (i < group.seed_indices.size() - 1) { + ss << ", "; + } + i++; + } + ss << "})"; + if (j < candidate_groups.size() - 1) { + ss << ", "; + } + j++; + } + return ss.str(); + }() + ); LOG_TRACE_L1(m_logger, "Validating candidate groups with flux analysis..."); - const std::vector validated_groups = validateGroupsWithFluxAnalysis(candidate_groups, Y, T9, rho); + const auto [validated_groups, invalidate_groups] = validateGroupsWithFluxAnalysis(candidate_groups, Y, T9, rho); LOG_TRACE_L1( m_logger, "Validated {} group(s) QSE groups. {}", @@ -507,6 +550,13 @@ namespace gridfire { }() ); + // Push the invalidated groups' species into the dynamic set + for (const auto& group : invalidate_groups) { + for (const auto& index : group.algebraic_indices) { + m_dynamic_species.push_back(m_baseEngine.getNetworkSpecies()[index]); + } + } + m_qse_groups = validated_groups; LOG_TRACE_L1(m_logger, "Identified {} QSE groups.", m_qse_groups.size()); @@ -542,6 +592,10 @@ namespace gridfire { m_qse_groups.size() == 1 ? "" : "s" ); + // throw std::runtime_error( + // "Partitioning complete. Throwing an error to end the program during debugging. This error should not be caught by the caller. " + // ); + } void MultiscalePartitioningEngineView::partitionNetwork( @@ -846,6 +900,8 @@ namespace gridfire { return m_baseEngine.getSpeciesIndex(species); } + + std::vector> MultiscalePartitioningEngineView::partitionByTimescale( const std::vector& Y_full, const double T9, @@ -853,19 +909,31 @@ namespace gridfire { ) const { LOG_TRACE_L1(m_logger, "Partitioning by timescale..."); const auto result= m_baseEngine.getSpeciesDestructionTimescales(Y_full, T9, rho); + const auto netTimescale = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho); + std::cout << gridfire::utils::formatNuclearTimescaleLogString(m_baseEngine, Y_full, T9, rho) << std::endl; if (!result) { - LOG_ERROR(m_logger, "Failed to get species timescales due to stale engine state"); + LOG_ERROR(m_logger, "Failed to get species destruction timescales due to stale engine state"); m_logger->flush_log(); - throw exceptions::StaleEngineError("Failed to get species timescales due to stale engine state"); + throw exceptions::StaleEngineError("Failed to get species destruction timescales due to stale engine state"); + } + if (!netTimescale) { + LOG_ERROR(m_logger, "Failed to get net species timescales due to stale engine state"); + m_logger->flush_log(); + throw exceptions::StaleEngineError("Failed to get net species timescales due to stale engine state"); } const std::unordered_map& all_timescales = result.value(); + const std::unordered_map& net_timescales = netTimescale.value(); const auto& all_species = m_baseEngine.getNetworkSpecies(); std::vector> sorted_timescales; for (size_t i = 0; i < all_species.size(); ++i) { double timescale = all_timescales.at(all_species[i]); + double net_timescale = net_timescales.at(all_species[i]); if (std::isfinite(timescale) && timescale > 0) { + LOG_TRACE_L3(m_logger, "Species {} has finite destruction timescale: destruction: {} s, net: {} s", all_species[i].name(), timescale, net_timescale); sorted_timescales.emplace_back(timescale, i); + } else { + LOG_TRACE_L3(m_logger, "Species {} has infinite or negative destruction timescale: destruction: {} s, net: {} s", all_species[i].name(), timescale, net_timescale); } } @@ -971,98 +1039,120 @@ namespace gridfire { } - std::vector + std::pair, std::vector> MultiscalePartitioningEngineView::validateGroupsWithFluxAnalysis( const std::vector &candidate_groups, const std::vector &Y, const double T9, const double rho ) const { - constexpr double FLUX_RATIO_THRESHOLD = 100; - std::vector validated_groups = candidate_groups; - for (auto& group : validated_groups) { - double internal_flux = 0.0; - double external_flux = 0.0; - - const std::unordered_set group_members( - group.species_indices.begin(), - group.species_indices.end() + constexpr double FLUX_RATIO_THRESHOLD = 5; + std::vector validated_groups; + std::vector invalidated_groups; + validated_groups.reserve(candidate_groups.size()); + for (auto& group : candidate_groups) { + const std::unordered_set algebraic_group_members( + group.algebraic_indices.begin(), + group.algebraic_indices.end() ); + const std::unordered_set seed_group_members( + group.seed_indices.begin(), + group.seed_indices.end() + ); + + double coupling_flux = 0.0; + double leakage_flux = 0.0; + for (const auto& reaction: m_baseEngine.getNetworkReactions()) { const double flow = std::abs(m_baseEngine.calculateMolarReactionFlow(*reaction, Y, T9, rho)); if (flow == 0.0) { continue; // Skip reactions with zero flow } - bool has_internal_reactant = false; - bool has_external_reactant = false; + bool has_internal_algebraic_reactant = false; for (const auto& reactant : reaction->reactants()) { - if (group_members.contains(m_baseEngine.getSpeciesIndex(reactant))) { - has_internal_reactant = true; - } else { - has_external_reactant = true; + if (algebraic_group_members.contains(m_baseEngine.getSpeciesIndex(reactant))) { + has_internal_algebraic_reactant = true; } } - bool has_internal_product = false; - bool has_external_product = false; + bool has_internal_algebraic_product = false; for (const auto& product : reaction->products()) { - if (group_members.contains(m_baseEngine.getSpeciesIndex(product))) { - has_internal_product = true; - } else { - has_external_product = true; + if (algebraic_group_members.contains(m_baseEngine.getSpeciesIndex(product))) { + has_internal_algebraic_product = true; } } - // Classify the reaction based on its participants - if ((has_internal_reactant && has_internal_product) && !(has_external_reactant || has_external_product)) { - LOG_TRACE_L3( - m_logger, - "Reaction {} is internal to the group containing {} and contributes to internal flux by {}", - reaction->id(), - [&]() -> std::string { - std::stringstream ss; - int count = 0; - for (const auto& idx : group.algebraic_indices) { - ss << m_baseEngine.getNetworkSpecies()[idx].name(); - if (count < group.species_indices.size() - 1) { - ss << ", "; - } - count++; - } - return ss.str(); - }(), - flow - ); - internal_flux += flow; // Internal flux within the group - } else if ((has_internal_reactant || has_internal_product) && (has_external_reactant || has_external_product)) { - LOG_TRACE_L3( - m_logger, - "Reaction {} is external to the group containing {} and contributes to external flux by {}", - reaction->id(), - [&]() -> std::string { - std::stringstream ss; - int count = 0; - for (const auto& idx : group.algebraic_indices) { - ss << m_baseEngine.getNetworkSpecies()[idx].name(); - if (count < group.species_indices.size() - 1) { - ss << ", "; - } - count++; - } - return ss.str(); - }(), - flow - ); - external_flux += flow; // External flux to/from the group + if (!has_internal_algebraic_product && !has_internal_algebraic_reactant) { + LOG_TRACE_L3(m_logger, "{}: Skipping reaction {} as it has no internal algebraic species in reactants or products.", group.toString(m_baseEngine), reaction->id()); + continue; } - // otherwise the reaction is fully decoupled from the QSE group and can be ignored. + + double algebraic_participants = 0; + double seed_participants = 0; + double external_participants = 0; + + std::unordered_set participants; + for(const auto& p : reaction->reactants()) participants.insert(p); + for(const auto& p : reaction->products()) participants.insert(p); + + for (const auto& species : participants) { + const size_t species_idx = m_baseEngine.getSpeciesIndex(species); + if (algebraic_group_members.contains(species_idx)) { + LOG_TRACE_L3(m_logger, "{}: Species {} is an algebraic participant in reaction {}.", group.toString(m_baseEngine), species.name(), reaction->id()); + algebraic_participants++; + } else if (seed_group_members.contains(species_idx)) { + LOG_TRACE_L3(m_logger, "{}: Species {} is a seed participant in reaction {}.", group.toString(m_baseEngine), species.name(), reaction->id()); + seed_participants++; + } else { + LOG_TRACE_L3(m_logger, "{}: Species {} is an external participant in reaction {}.", group.toString(m_baseEngine), species.name(), reaction->id()); + external_participants++; + } + } + + const double total_participants = algebraic_participants + seed_participants + external_participants; + + if (total_participants == 0) { + LOG_CRITICAL(m_logger, "Some catastrophic error has occurred. Reaction {} has no participants.", reaction->id()); + throw std::runtime_error("Some catastrophic error has occurred. Reaction " + std::string(reaction->id()) + " has no participants."); + } + + const double leakage_fraction = external_participants / total_participants; + const double coupling_fraction = (algebraic_participants + seed_participants) / total_participants; + + leakage_flux += flow * leakage_fraction; + coupling_flux += flow * coupling_fraction; } - if (internal_flux > FLUX_RATIO_THRESHOLD * external_flux) { + + // if (leakage_flux < 1e-99) { + // LOG_TRACE_L1( + // m_logger, + // "Group containing {} is in equilibrium due to vanishing leakage: leakage flux = {}, coupling flux = {}, ratio = {}", + // [&]() -> std::string { + // std::stringstream ss; + // int count = 0; + // for (const auto& idx : group.algebraic_indices) { + // ss << m_baseEngine.getNetworkSpecies()[idx].name(); + // if (count < group.species_indices.size() - 1) { + // ss << ", "; + // } + // count++; + // } + // return ss.str(); + // }(), + // leakage_flux, + // coupling_flux, + // coupling_flux / leakage_flux + // ); + // validated_groups.emplace_back(group); + // validated_groups.back().is_in_equilibrium = true; + // } else if ((coupling_flux / leakage_flux ) > FLUX_RATIO_THRESHOLD) { + if ((coupling_flux / leakage_flux ) > FLUX_RATIO_THRESHOLD) { LOG_TRACE_L1( m_logger, - "Group containing {} is in equilibrium: internal flux = {}, external flux = {}, ratio = {}", + "Group containing {} is in equilibrium due to high coupling flux threshold: leakage flux = {}, coupling flux = {}, ratio = {} (Threshold: {})", [&]() -> std::string { std::stringstream ss; int count = 0; @@ -1075,15 +1165,17 @@ namespace gridfire { } return ss.str(); }(), - internal_flux, - external_flux, - internal_flux / external_flux + leakage_flux, + coupling_flux, + coupling_flux / leakage_flux, + FLUX_RATIO_THRESHOLD ); - group.is_in_equilibrium = true; // This group is in equilibrium if internal flux is significantly larger than external flux. + validated_groups.emplace_back(group); + validated_groups.back().is_in_equilibrium = true; } else { LOG_TRACE_L1( m_logger, - "Group containing {} is NOT in equilibrium: internal flux = {}, external flux = {}, ratio = {}", + "Group containing {} is NOT in equilibrium: leakage flux = {}, coupling flux = {}, ratio = {} (Threshold: {})", [&]() -> std::string { std::stringstream ss; int count = 0; @@ -1096,14 +1188,16 @@ namespace gridfire { } return ss.str(); }(), - internal_flux, - external_flux, - internal_flux / external_flux + leakage_flux, + coupling_flux, + coupling_flux / leakage_flux, + FLUX_RATIO_THRESHOLD ); - group.is_in_equilibrium = false; + invalidated_groups.emplace_back(group); + invalidated_groups.back().is_in_equilibrium = false; } } - return validated_groups; + return {validated_groups, invalidated_groups}; } std::vector MultiscalePartitioningEngineView::solveQSEAbundances( @@ -1543,10 +1637,11 @@ namespace gridfire { hash_combine(seed, bin(m_T9, m_cacheConfig.T9_tol)); hash_combine(seed, bin(m_rho, m_cacheConfig.rho_tol)); + double negThresh = 1e-10; // Threshold for considering a value as negative. for (double Yi : m_Y) { - if (Yi < 0.0 && std::abs(Yi) < 1e-20) { + if (Yi < 0.0 && std::abs(Yi) < negThresh) { Yi = 0.0; // Avoid negative abundances - } else if (Yi < 0.0 && std::abs(Yi) >= 1e-20) { + } else if (Yi < 0.0 && std::abs(Yi) >= negThresh) { throw std::invalid_argument("Yi should be positive for valid hashing (expected Yi > 0, received " + std::to_string(Yi) + ")"); } hash_combine(seed, bin(Yi, m_cacheConfig.Yi_tol)); @@ -1580,6 +1675,55 @@ namespace gridfire { return !(*this == other); } + std::string MultiscalePartitioningEngineView::QSEGroup::toString() const { + std::stringstream ss; + ss << "QSEGroup(Algebraic: ["; + size_t count = 0; + for (const auto& idx : algebraic_indices) { + ss << idx; + if (count < algebraic_indices.size() - 1) { + ss << ", "; + } + count++; + } + ss << "], Seed: ["; + count = 0; + for (const auto& idx : seed_indices) { + ss << idx; + if (count < seed_indices.size() - 1) { + ss << ", "; + } + count++; + } + ss << "], Mean Timescale: " << mean_timescale << ", Is In Equilibrium: " << (is_in_equilibrium ? "True" : "False") << ")"; + return ss.str(); + + } + + std::string MultiscalePartitioningEngineView::QSEGroup::toString(DynamicEngine &engine) const { + std::stringstream ss; + ss << "QSEGroup(Algebraic: ["; + size_t count = 0; + for (const auto& idx : algebraic_indices) { + ss << engine.getNetworkSpecies()[idx].name(); + if (count < algebraic_indices.size() - 1) { + ss << ", "; + } + count++; + } + ss << "], Seed: ["; + count = 0; + for (const auto& idx : seed_indices) { + ss << engine.getNetworkSpecies()[idx].name(); + if (count < seed_indices.size() - 1) { + ss << ", "; + } + count++; + } + ss << "], Mean Timescale: " << mean_timescale << ", Is In Equilibrium: " << (is_in_equilibrium ? "True" : "False") << ")"; + return ss.str(); + } + void MultiscalePartitioningEngineView::CacheStats::hit(const operators op) { if (op == operators::All) { throw std::invalid_argument("Cannot use 'ALL' as an operator for a hit"); diff --git a/src/lib/solver/solver.cpp b/src/lib/solver/solver.cpp index 33385af2..72d64c5f 100644 --- a/src/lib/solver/solver.cpp +++ b/src/lib/solver/solver.cpp @@ -137,7 +137,7 @@ namespace gridfire::solver { std::vector speciesNames; speciesNames.reserve(numSpecies); for (const auto& species : m_engine.getNetworkSpecies()) { - speciesNames.push_back(std::string(species.name())); + speciesNames.emplace_back(species.name()); } Composition outputComposition(speciesNames); @@ -145,10 +145,19 @@ namespace gridfire::solver { outputComposition.finalize(true); NetOut netOut; - netOut.composition = std::move(outputComposition); + netOut.composition = outputComposition; netOut.energy = accumulated_energy; // Specific energy rate netOut.num_steps = manager.m_num_steps; + auto [dEps_dT, dEps_dRho] = m_engine.calculateEpsDerivatives( + std::vector(Y.begin(), Y.begin() + numSpecies), // TODO: This narrowing should probably be solved. Its possible unforeseen bugs will arise from this + T9, + netIn.density + ); + + netOut.dEps_dT = dEps_dT; + netOut.dEps_dRho = dEps_dRho; + return netOut; } diff --git a/src/lib/solver/strategies/CVODE_solver_strategy.cpp b/src/lib/solver/strategies/CVODE_solver_strategy.cpp index 22fb48b4..df990fbb 100644 --- a/src/lib/solver/strategies/CVODE_solver_strategy.cpp +++ b/src/lib/solver/strategies/CVODE_solver_strategy.cpp @@ -1,6 +1,10 @@ #include "gridfire/solver/strategies/CVODE_solver_strategy.h" #include "gridfire/network.h" +#include "gridfire/utils/table_format.h" +#include "gridfire/engine/diagnostics/dynamic_engine_diagnostics.h" + +#include "quill/LogMacros.h" #include "fourdst/composition/composition.h" @@ -10,6 +14,9 @@ #include #include #include +#include + +#include "fourdst/composition/exceptions/exceptions_composition.h" namespace { @@ -75,11 +82,12 @@ namespace gridfire::solver { // TODO: In order to support MPI this function must be changed const int flag = SUNContext_Create(SUN_COMM_NULL, &m_sun_ctx); if (flag < 0) { - throw std::runtime_error("Failed to create SUNDIALS context (Errno: " + std::to_string(flag) + ")"); + throw std::runtime_error("Failed to create SUNDIALS context (SUNDIALS Errno: " + std::to_string(flag) + ")"); } } CVODESolverStrategy::~CVODESolverStrategy() { + std::cout << "Cleaning up CVODE resources..." << std::endl; cleanup_cvode_resources(true); if (m_sun_ctx) { @@ -111,7 +119,7 @@ namespace gridfire::solver { [[maybe_unused]] double last_callback_time = 0; m_num_steps = 0; double accumulated_energy = 0.0; - + int total_update_stages_triggered = 0; while (current_time < netIn.tMax) { try { user_data.T9 = T9; @@ -138,22 +146,47 @@ namespace gridfire::solver { double last_step_size; CVodeGetNumSteps(m_cvode_mem, &n_steps); CVodeGetLastStep(m_cvode_mem, &last_step_size); - std::cout << std::scientific << std::setprecision(3) << "Step: " << std::setw(6) << n_steps << " | Time: " << current_time << " | Last Step Size: " << last_step_size << std::endl; + long int nliters, nlcfails; + CVodeGetNumNonlinSolvIters(m_cvode_mem, &nliters); + CVodeGetNumNonlinSolvConvFails(m_cvode_mem, &nlcfails); + + sunrealtype* y_data = N_VGetArrayPointer(m_Y); + const double current_energy = y_data[numSpecies]; // Specific energy rate + std::cout << std::scientific << std::setprecision(3) + << "Step: " << std::setw(6) << n_steps + << " | Time: " << current_time << " [s]" + << " | Last Step Size: " << last_step_size + << " | Accumulated Energy: " << current_energy << " [erg/g]" + << " | NonlinIters: " << std::setw(2) << nliters + << " | ConvFails: " << std::setw(2) << nlcfails + << std::endl; + + // if (n_steps % 50 == 0) { + // std::cout << "Logging step diagnostics at step " << n_steps << "..." << std::endl; + // log_step_diagnostics(user_data); + // } + // if (n_steps == 300) { + // log_step_diagnostics(user_data); + // exit(0); + // } + // log_step_diagnostics(user_data); } catch (const exceptions::StaleEngineTrigger& e) { exceptions::StaleEngineTrigger::state staleState = e.getState(); accumulated_energy += e.energy(); // Add the specific energy rate to the accumulated energy - // total_update_stages_triggered++; + LOG_INFO( + m_logger, + "Engine Update Triggered due to StaleEngineTrigger exception at time {} ({} update{} triggered). Current total specific energy {} [erg/g]", + current_time, + total_update_stages_triggered, + total_update_stages_triggered == 1 ? "" : "s", + accumulated_energy); + total_update_stages_triggered++; fourdst::composition::Composition temp_comp; std::vector mass_fractions; size_t num_species_at_stop = e.numSpecies(); - if (num_species_at_stop != m_engine.getNetworkSpecies().size()) { - throw std::runtime_error( - "StaleEngineError state has a different number of species than the engine. This should not happen." - ); - } mass_fractions.reserve(num_species_at_stop); for (size_t i = 0; i < num_species_at_stop; ++i) { @@ -170,13 +203,22 @@ namespace gridfire::solver { netInTemp.composition = temp_comp; fourdst::composition::Composition currentComposition = m_engine.update(netInTemp); + LOG_INFO( + m_logger, + "Due to a triggered stale engine the composition was updated from size {} to {} species.", + num_species_at_stop, + m_engine.getNetworkSpecies().size() + ); numSpecies = m_engine.getNetworkSpecies().size(); N = numSpecies + 1; - initialize_cvode_integration_resources(N, numSpecies, current_time, temp_comp, absTol, relTol, accumulated_energy); + 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"); + } catch (fourdst::composition::exceptions::InvalidCompositionError& e) { + log_step_diagnostics(user_data); + std::rethrow_exception(std::make_exception_ptr(e)); } } @@ -195,7 +237,7 @@ namespace gridfire::solver { std::vector speciesNames; speciesNames.reserve(numSpecies); for (const auto& species : m_engine.getNetworkSpecies()) { - speciesNames.push_back(std::string(species.name())); + speciesNames.emplace_back(species.name()); } fourdst::composition::Composition outputComposition(speciesNames); @@ -203,9 +245,21 @@ namespace gridfire::solver { outputComposition.finalize(true); NetOut netOut; - netOut.composition = std::move(outputComposition); + netOut.composition = outputComposition; netOut.energy = accumulated_energy; check_cvode_flag(CVodeGetNumSteps(m_cvode_mem, reinterpret_cast(&netOut.num_steps)), "CVodeGetNumSteps"); + + outputComposition.setCompositionMode(false); // set to number fraction mode + std::vector Y = outputComposition.getNumberFractionVector(); // TODO need to ensure that the canonical vector representation is used throughout the code to make sure tracking does not get messed up + auto [dEps_dT, dEps_dRho] = m_engine.calculateEpsDerivatives( + std::vector(Y.begin(), Y.begin() + numSpecies), // TODO: This narrowing should probably be solved. Its possible unforeseen bugs will arise from this + T9, + netIn.density + ); + + netOut.dEps_dT = dEps_dT; + netOut.dEps_dRho = dEps_dRho; + return netOut; } @@ -232,10 +286,11 @@ namespace gridfire::solver { void *user_data ) { auto* data = static_cast(user_data); - auto* instance = data->solver_instance; + const auto* instance = data->solver_instance; try { - return instance->calculate_rhs(t, y, ydot, data); + instance->calculate_rhs(t, y, ydot, data); + return 0; } catch (const exceptions::StaleEngineTrigger& e) { data->captured_exception = std::make_unique(e); return 1; // 1 Indicates a recoverable error, CVODE will retry the step @@ -279,7 +334,7 @@ namespace gridfire::solver { return 0; } - int CVODESolverStrategy::calculate_rhs( + void CVODESolverStrategy::calculate_rhs( const sunrealtype t, const N_Vector y, const N_Vector ydot, @@ -304,7 +359,6 @@ namespace gridfire::solver { ydot_data[i] = dydt[i]; } ydot_data[numSpecies] = nuclearEnergyGenerationRate; // Set the last element to the specific energy rate - return 0; } void CVODESolverStrategy::initialize_cvode_integration_resources( @@ -319,6 +373,7 @@ namespace gridfire::solver { cleanup_cvode_resources(false); // Cleanup any existing resources before initializing new ones m_Y = init_sun_vector(N, m_sun_ctx); + m_YErr = N_VClone(m_Y); sunrealtype *y_data = N_VGetArrayPointer(m_Y); for (size_t i = 0; i < numSpecies; i++) { @@ -335,6 +390,8 @@ namespace gridfire::solver { check_cvode_flag(CVodeInit(m_cvode_mem, cvode_rhs_wrapper, current_time, m_Y), "CVodeInit"); check_cvode_flag(CVodeSStolerances(m_cvode_mem, relTol, absTol), "CVodeSStolerances"); + check_cvode_flag(CVodeSetMaxStep(m_cvode_mem, 1.0e20), "CVodeSetMaxStep"); + m_J = SUNDenseMatrix(static_cast(N), static_cast(N), m_sun_ctx); check_cvode_flag(m_J == nullptr ? -1 : 0, "SUNDenseMatrix"); m_LS = SUNLinSol_Dense(m_Y, m_J, m_sun_ctx); @@ -348,10 +405,12 @@ namespace gridfire::solver { if (m_LS) SUNLinSolFree(m_LS); if (m_J) SUNMatDestroy(m_J); if (m_Y) N_VDestroy(m_Y); + if (m_YErr) N_VDestroy(m_YErr); m_LS = nullptr; m_J = nullptr; m_Y = nullptr; + m_YErr = nullptr; if (memFree) { if (m_cvode_mem) CVodeFree(&m_cvode_mem); @@ -359,4 +418,79 @@ namespace gridfire::solver { } } + void CVODESolverStrategy::log_step_diagnostics(const CVODEUserData &user_data) const { + check_cvode_flag(CVodeGetEstLocalErrors(m_cvode_mem, m_YErr), "CVodeGetEstLocalErrors"); + + sunrealtype *y_data = N_VGetArrayPointer(m_Y); + sunrealtype *y_err_data = N_VGetArrayPointer(m_YErr); + + + std::vector err_ratios; + std::vector speciesNames; + + const auto absTol = m_config.get("gridfire:solver:CVODESolverStrategy:absTol", 1.0e-8); + const auto relTol = m_config.get("gridfire:solver:CVODESolverStrategy:relTol", 1.0e-8); + + const size_t num_components = N_VGetLength(m_Y); + err_ratios.resize(num_components - 1); + + std::vector Y_full(y_data, y_data + num_components - 1); + + std::ranges::replace_if( + Y_full, + [](const double val) { + return val < 0.0 && val > -1.0e-16; + }, + 0.0 + ); + + for (size_t i = 0; i < num_components - 1; i++) { + const double weight = relTol * std::abs(y_data[i]) + absTol; + if (weight == 0.0) continue; // Skip components with zero weight + + const double err_ratio = std::abs(y_err_data[i]) / weight; + + err_ratios[i] = err_ratio; + speciesNames.push_back(std::string(user_data.networkSpecies->at(i).name())); + } + + if (err_ratios.empty()) { + return; + } + + std::vector indices(speciesNames.size()); + for (size_t i = 0; i < indices.size(); ++i) { + indices[i] = i; + } + + std::ranges::sort( + indices, + [&err_ratios](const size_t i1, const size_t i2) { + return err_ratios[i1] > err_ratios[i2]; + } + ); + + std::vector sorted_speciesNames; + std::vector sorted_err_ratios; + + sorted_speciesNames.reserve(indices.size()); + sorted_err_ratios.reserve(indices.size()); + + for (const auto idx: indices) { + sorted_speciesNames.push_back(speciesNames[idx]); + sorted_err_ratios.push_back(err_ratios[idx]); + } + + + + std::vector> columns; + columns.push_back(std::make_unique>("Species", sorted_speciesNames)); + columns.push_back(std::make_unique>("Error Ratio", sorted_err_ratios)); + + std::cout << utils::format_table("Species Error Ratios", columns) << std::endl; + diagnostics::inspect_jacobian_stiffness(*user_data.engine, Y_full, user_data.T9, user_data.rho); + diagnostics::inspect_species_balance(*user_data.engine, "N-14", Y_full, user_data.T9, user_data.rho); + diagnostics::inspect_species_balance(*user_data.engine, "n-1", Y_full, user_data.T9, user_data.rho); + + } } diff --git a/src/meson.build b/src/meson.build index 76320efd..adff336f 100644 --- a/src/meson.build +++ b/src/meson.build @@ -9,6 +9,7 @@ gridfire_sources = files( 'lib/engine/views/engine_priming.cpp', 'lib/engine/procedures/priming.cpp', 'lib/engine/procedures/construction.cpp', + 'lib/engine/diagnostics/dynamic_engine_diagnostics.cpp', 'lib/reaction/reaction.cpp', 'lib/reaction/reaclib.cpp', 'lib/io/network_file.cpp', @@ -53,4 +54,9 @@ gridfire_dep = declare_dependency( install_subdir('include/gridfire', install_dir: get_option('includedir')) -subdir('python') +if get_option('build-python') + message('Configuring Python bindings...') + subdir('src-pybind') +else + message('Skipping Python bindings...') +endif \ No newline at end of file diff --git a/subprojects/fourdst.wrap b/subprojects/fourdst.wrap index bf00980d..fc133ce0 100644 --- a/subprojects/fourdst.wrap +++ b/subprojects/fourdst.wrap @@ -1,4 +1,4 @@ [wrap-git] url = https://github.com/4D-STAR/fourdst -revision = v0.7.0 +revision = v0.8.0 depth = 1 diff --git a/tests/graphnet_sandbox/main.cpp b/tests/graphnet_sandbox/main.cpp index 5bba4689..d641fae2 100644 --- a/tests/graphnet_sandbox/main.cpp +++ b/tests/graphnet_sandbox/main.cpp @@ -99,7 +99,7 @@ int main(int argc, char* argv[]){ g_previousHandler = std::set_terminate(quill_terminate_handler); quill::Logger* logger = fourdst::logging::LogManager::getInstance().getLogger("log"); - logger->set_log_level(quill::LogLevel::TraceL1); + logger->set_log_level(quill::LogLevel::TraceL3); LOG_INFO(logger, "Starting Adaptive Engine View Example..."); using namespace gridfire; @@ -122,12 +122,15 @@ int main(int argc, char* argv[]){ netIn.temperature = 1.5e7; netIn.density = 1.6e2; netIn.energy = 0; - netIn.tMax = 5e17; + netIn.tMax = 1e2; // netIn.tMax = 1e-14; netIn.dt0 = 1e-12; GraphEngine ReaclibEngine(composition, partitionFunction, NetworkBuildDepth::SecondOrder); + ReaclibEngine.setUseReverseReactions(false); + ReaclibEngine.setPrecomputation(false); + MultiscalePartitioningEngineView partitioningView(ReaclibEngine); AdaptiveEngineView adaptiveView(partitioningView); @@ -144,5 +147,6 @@ int main(int argc, char* argv[]){ double finalHydrogen = netOut.composition.getMassFraction("H-1"); double fractionalConsumedHydrogen = (initialHydrogen - finalHydrogen) / initialHydrogen * 100.0; std::cout << "Fractional consumed hydrogen: " << fractionalConsumedHydrogen << "%" << std::endl; + std::cout << netOut << std::endl; } \ No newline at end of file