From 8a0b5b2c3611efbd0f7f625575e9233aa821ffe7 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Tue, 7 Oct 2025 15:16:03 -0400 Subject: [PATCH] feat(weak): major weak rate progress Major weak rate progress which includes: A refactor of many of the public interfaces for GridFire Engines to use composition objects as opposed to raw abundance vectors. This helps prevent index mismatch errors. Further, the weak reaction class has been expanded with the majority of an implimentation, including an atomic_base derived class to allow for proper auto diff tracking of the interpolated table results. Some additional changes are that the version of fourdst and libcomposition have been bumped to versions with smarter caching of intermediate vectors and a few bug fixes. --- .gitignore | 1 + .../diagnostics/dynamic_engine_diagnostics.h | 4 +- src/include/gridfire/engine/engine_abstract.h | 48 +- src/include/gridfire/engine/engine_graph.h | 223 +++--- .../gridfire/engine/procedures/construction.h | 9 +- .../gridfire/engine/procedures/priming.h | 8 +- .../gridfire/engine/views/engine_adaptive.h | 87 +-- .../gridfire/engine/views/engine_defined.h | 58 +- .../gridfire/engine/views/engine_multiscale.h | 153 ++-- .../gridfire/exceptions/error_engine.h | 2 +- .../gridfire/expectations/expected_engine.h | 5 +- src/include/gridfire/io/network_file.h | 4 +- src/include/gridfire/network.h | 2 +- .../partition/composite/partition_composite.h | 2 +- src/include/gridfire/reaction/reaction.h | 86 ++- src/include/gridfire/reaction/weak/weak.h | 275 ++++++- .../reaction/weak/weak_interpolator.h | 44 ++ .../reaction/weak/weak_rate_library.h | 19 +- .../gridfire/reaction/weak/weak_types.h | 133 ++++ .../gridfire/screening/screening_abstract.h | 20 +- .../gridfire/screening/screening_bare.h | 28 +- .../screening/screening_intermediate.h | 2 +- .../gridfire/screening/screening_weak.h | 24 +- src/include/gridfire/solver/solver.h | 251 ------- .../solver/strategies/CVODE_solver_strategy.h | 10 +- .../triggers/engine_partitioning_trigger.h | 6 +- .../trigger/procedures/trigger_pprint.h | 2 +- .../gridfire/trigger/trigger_abstract.h | 11 +- src/include/gridfire/utils/logging.h | 11 +- src/include/gridfire/utils/table_format.h | 25 +- .../dynamic_engine_diagnostics.cpp | 39 +- src/lib/engine/engine_graph.cpp | 201 ++++-- src/lib/engine/procedures/construction.cpp | 106 ++- src/lib/engine/procedures/priming.cpp | 32 +- src/lib/engine/views/engine_adaptive.cpp | 119 ++- src/lib/engine/views/engine_defined.cpp | 57 +- src/lib/engine/views/engine_multiscale.cpp | 681 +++++++----------- src/lib/io/network_file.cpp | 2 +- src/lib/network.cpp | 3 +- .../partition_rauscher_thielemann.cpp | 1 - src/lib/reaction/reaction.cpp | 56 +- src/lib/reaction/weak/weak.cpp | 450 +++++++++++- src/lib/reaction/weak/weak_interpolator.cpp | 308 ++++++++ src/lib/screening/screening_bare.cpp | 2 +- src/lib/screening/screening_weak.cpp | 2 +- src/lib/solver/solver.cpp | 356 --------- .../strategies/CVODE_solver_strategy.cpp | 49 +- .../triggers/engine_partitioning_trigger.cpp | 21 +- src/lib/utils/logging.cpp | 4 +- src/meson.build | 3 +- src/python/reaction/bindings.cpp | 4 +- subprojects/fourdst.wrap | 2 +- tests/graphnet_sandbox/main.cpp | 18 - 53 files changed, 2310 insertions(+), 1759 deletions(-) create mode 100644 src/include/gridfire/reaction/weak/weak_interpolator.h create mode 100644 src/include/gridfire/reaction/weak/weak_types.h create mode 100644 src/lib/reaction/weak/weak_interpolator.cpp delete mode 100644 src/lib/solver/solver.cpp diff --git a/.gitignore b/.gitignore index feedac37..c6f3fa87 100644 --- a/.gitignore +++ b/.gitignore @@ -80,6 +80,7 @@ subprojects/libplugin/ subprojects/minizip-ng-4.0.10/ subprojects/cvode-*/ *.fbundle +*.wraplock *.csv *.dot diff --git a/src/include/gridfire/engine/diagnostics/dynamic_engine_diagnostics.h b/src/include/gridfire/engine/diagnostics/dynamic_engine_diagnostics.h index a65d9bd5..e08c990a 100644 --- a/src/include/gridfire/engine/diagnostics/dynamic_engine_diagnostics.h +++ b/src/include/gridfire/engine/diagnostics/dynamic_engine_diagnostics.h @@ -19,14 +19,14 @@ namespace gridfire::diagnostics { void inspect_species_balance( const DynamicEngine& engine, const std::string& species_name, - const std::vector& Y_full, + const fourdst::composition::Composition &comp, double T9, double rho ); void inspect_jacobian_stiffness( const DynamicEngine& engine, - const std::vector& Y_full, + const fourdst::composition::Composition &comp, double T9, double rho ); diff --git a/src/include/gridfire/engine/engine_abstract.h b/src/include/gridfire/engine/engine_abstract.h index c9658426..d74d2368 100644 --- a/src/include/gridfire/engine/engine_abstract.h +++ b/src/include/gridfire/engine/engine_abstract.h @@ -58,8 +58,10 @@ namespace gridfire { */ template struct StepDerivatives { - std::vector dydt; ///< Derivatives of abundances (dY/dt for each species). + std::map dydt; ///< Derivatives of abundances (dY/dt for each species). T nuclearEnergyGenerationRate = T(0.0); ///< Specific energy generation rate (e.g., erg/g/s). + + StepDerivatives() : dydt(), nuclearEnergyGenerationRate(T(0.0)) {} }; using SparsityPattern = std::vector>; @@ -107,7 +109,7 @@ namespace gridfire { /** * @brief Calculate the right-hand side (dY/dt) and energy generation. * - * @param Y Vector of current abundances for all species. + * @param comp Composition object containing current abundances. * @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. @@ -117,7 +119,7 @@ namespace gridfire { * rate for the current state. */ [[nodiscard]] virtual std::expected, expectations::StaleEngineError> calculateRHSAndEnergy( - const std::vector& Y, + const fourdst::composition::Composition& comp, double T9, double rho ) const = 0; @@ -142,7 +144,7 @@ namespace gridfire { /** * @brief Generate the Jacobian matrix for the current state. * - * @param Y_dynamic Vector of current abundances. + * @param comp Composition object containing current abundances. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @@ -150,13 +152,13 @@ namespace gridfire { * for the current state. The matrix can then be accessed via getJacobianMatrixEntry(). */ virtual void generateJacobianMatrix( - const std::vector& Y_dynamic, + const fourdst::composition::Composition& comp, double T9, double rho ) const = 0; virtual void generateJacobianMatrix( - const std::vector& Y_dynamic, + const fourdst::composition::Composition& comp, double T9, double rho, const SparsityPattern& sparsityPattern @@ -167,15 +169,15 @@ namespace gridfire { /** * @brief Get an entry from the previously generated Jacobian matrix. * - * @param i Row index (species index). - * @param j Column index (species index). + * @param rowSpecies The species corresponding to the row index (i) + * @param colSpecies The species corresponding to the column index (j) * @return Value of the Jacobian matrix at (i, j). * * The Jacobian must have been generated by generateJacobianMatrix() before calling this. */ [[nodiscard]] virtual double getJacobianMatrixEntry( - int i, - int j + const fourdst::atomic::Species& rowSpecies, + const fourdst::atomic::Species& colSpecies ) const = 0; @@ -190,22 +192,22 @@ namespace gridfire { /** * @brief Get an entry from the stoichiometry matrix. * - * @param speciesIndex Index of the species. - * @param reactionIndex Index of the reaction. + * @param species species to look up stoichiometry for. + * @param reaction reaction to find * @return Stoichiometric coefficient for the species in the reaction. * * The stoichiometry matrix must have been generated by generateStoichiometryMatrix(). */ [[nodiscard]] virtual int getStoichiometryMatrixEntry( - int speciesIndex, - int reactionIndex + const fourdst::atomic::Species& species, + const reaction::Reaction& reaction ) const = 0; /** * @brief Calculate the molar reaction flow for a given reaction. * * @param reaction The reaction for which to calculate the flow. - * @param Y Vector of current abundances. + * @param comp Composition object containing current abundances. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return Molar flow rate for the reaction (e.g., mol/g/s). @@ -215,7 +217,7 @@ namespace gridfire { */ [[nodiscard]] virtual double calculateMolarReactionFlow( const reaction::Reaction& reaction, - const std::vector& Y, + const fourdst::composition::Composition& comp, double T9, double rho ) const = 0; @@ -223,7 +225,7 @@ namespace gridfire { /** * @brief Calculate the derivatives of the energy generation rate with respect to T and rho. * - * @param Y Vector of current abundances. + * @param comp Composition object containing 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. @@ -232,9 +234,9 @@ namespace gridfire { * 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 fourdst::composition::Composition& comp, + double T9, + double rho ) const = 0; /** @@ -249,7 +251,7 @@ namespace gridfire { /** * @brief Compute timescales for all species in the network. * - * @param Y Vector of current abundances. + * @param comp Composition object containing current abundances. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return Map from Species to their characteristic timescales (s). @@ -258,13 +260,13 @@ namespace gridfire { * which can be used for timestep control, diagnostics, and reaction network culling. */ [[nodiscard]] virtual std::expected, expectations::StaleEngineError> getSpeciesTimescales( - const std::vector& Y, + const fourdst::composition::Composition& comp, double T9, double rho ) const = 0; [[nodiscard]] virtual std::expected, expectations::StaleEngineError> getSpeciesDestructionTimescales( - const std::vector& Y, + const fourdst::composition::Composition& comp, double T9, double rho ) const = 0; diff --git a/src/include/gridfire/engine/engine_graph.h b/src/include/gridfire/engine/engine_graph.h index 0c918858..ef972995 100644 --- a/src/include/gridfire/engine/engine_graph.h +++ b/src/include/gridfire/engine/engine_graph.h @@ -24,6 +24,8 @@ #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" // PERF: The function getNetReactionStoichiometry returns a map of species to their stoichiometric coefficients for a given reaction. // this makes extra copies of the species, which is not ideal and could be optimized further. @@ -106,13 +108,13 @@ namespace gridfire { */ explicit GraphEngine( const fourdst::composition::Composition &composition, - const BuildDepthType = NetworkBuildDepth::Full + BuildDepthType = NetworkBuildDepth::Full ); explicit GraphEngine( const fourdst::composition::Composition &composition, const partition::PartitionFunction& partitionFunction, - const BuildDepthType buildDepth = NetworkBuildDepth::Full + BuildDepthType buildDepth = NetworkBuildDepth::Full ); /** @@ -128,7 +130,7 @@ namespace gridfire { /** * @brief Calculates the right-hand side (dY/dt) and energy generation rate. * - * @param Y Vector of current abundances for all species. + * @param comp Composition object containing current abundances. * @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. @@ -139,21 +141,21 @@ namespace gridfire { * @see StepDerivatives */ [[nodiscard]] std::expected, expectations::StaleEngineError> calculateRHSAndEnergy( - const std::vector& Y, - const double T9, - const double rho + const fourdst::composition::Composition& comp, + double T9, + double rho ) const override; [[nodiscard]] EnergyDerivatives calculateEpsDerivatives( - const std::vector& Y, - const double T9, - const double rho + const fourdst::composition::Composition& comp, + double T9, + double rho ) const override; /** * @brief Generates the Jacobian matrix for the current state. * - * @param Y_dynamic Vector of current abundances. + * @param comp Composition object containing current abundances. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @@ -164,13 +166,13 @@ namespace gridfire { * @see getJacobianMatrixEntry() */ void generateJacobianMatrix( - const std::vector& Y_dynamic, - const double T9, - const double rho + const fourdst::composition::Composition& comp, + double T9, + double rho ) const override; void generateJacobianMatrix( - const std::vector &Y_dynamic, + const fourdst::composition::Composition& comp, double T9, double rho, const SparsityPattern &sparsityPattern @@ -188,7 +190,7 @@ namespace gridfire { * @brief Calculates the molar reaction flow for a given reaction. * * @param reaction The reaction for which to calculate the flow. - * @param Y Vector of current abundances. + * @param comp Composition object containing current abundances. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return Molar flow rate for the reaction (e.g., mol/g/s). @@ -198,9 +200,9 @@ namespace gridfire { */ [[nodiscard]] double calculateMolarReactionFlow( const reaction::Reaction& reaction, - const std::vector&Y, - const double T9, - const double rho + const fourdst::composition::Composition& comp, + double T9, + double rho ) const override; /** @@ -220,8 +222,8 @@ namespace gridfire { /** * @brief Gets an entry from the previously generated Jacobian matrix. * - * @param i Row index (species index). - * @param j Column index (species index). + * @param rowSpecies Species corresponding to the row index (i). + * @param colSpecies Species corresponding to the column index (j). * @return Value of the Jacobian matrix at (i, j). * * The Jacobian must have been generated by `generateJacobianMatrix()` before calling this. @@ -229,8 +231,8 @@ namespace gridfire { * @see generateJacobianMatrix() */ [[nodiscard]] double getJacobianMatrixEntry( - const int i, - const int j + const fourdst::atomic::Species& rowSpecies, + const fourdst::atomic::Species& colSpecies ) const override; /** @@ -246,8 +248,8 @@ namespace gridfire { /** * @brief Gets an entry from the stoichiometry matrix. * - * @param speciesIndex Index of the species. - * @param reactionIndex Index of the reaction. + * @param species Species to look up stoichiometry for. + * @param reaction Reaction to find. * @return Stoichiometric coefficient for the species in the reaction. * * The stoichiometry matrix must have been generated by `generateStoichiometryMatrix()`. @@ -255,14 +257,14 @@ namespace gridfire { * @see generateStoichiometryMatrix() */ [[nodiscard]] int getStoichiometryMatrixEntry( - const int speciesIndex, - const int reactionIndex + const fourdst::atomic::Species& species, + const reaction::Reaction& reaction ) const override; /** * @brief Computes timescales for all species in the network. * - * @param Y Vector of current abundances. + * @param comp Composition object containing current abundances. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return Map from Species to their characteristic timescales (s). @@ -271,13 +273,13 @@ namespace gridfire { * which can be used for timestep control or diagnostics. */ [[nodiscard]] std::expected, expectations::StaleEngineError> getSpeciesTimescales( - const std::vector& Y, + const fourdst::composition::Composition& comp, double T9, double rho ) const override; [[nodiscard]] std::expected, expectations::StaleEngineError> getSpeciesDestructionTimescales( - const std::vector& Y, + const fourdst::composition::Composition& comp, double T9, double rho ) const override; @@ -396,7 +398,7 @@ namespace gridfire { * @param reaction The reaction for which to calculate the reverse rate. * @param T9 Temperature in units of 10^9 K. * @param rho - * @param Y + * @param comp Composition object containing current abundances. * @return Reverse rate for the reaction (e.g., mol/g/s). * * This method computes the reverse rate based on the forward rate and @@ -404,7 +406,9 @@ namespace gridfire { */ [[nodiscard]] double calculateReverseRate( const reaction::Reaction &reaction, - double T9, double rho, const std::vector &Y + double T9, + double rho, + const fourdst::composition::Composition& comp ) const; /** @@ -428,8 +432,10 @@ namespace gridfire { [[nodiscard]] double calculateReverseRateTwoBodyDerivative( const reaction::Reaction &reaction, - const double T9, - double rho, const std::vector &Y, const double reverseRate + double T9, + double rho, + const fourdst::composition::Composition& comp, + double reverseRate ) const; /** @@ -577,12 +583,16 @@ namespace gridfire { constants m_constants; + rates::weak::WeakRateInterpolator m_weakRateInterpolator; ///< Interpolator for weak reaction rates. + reaction::ReactionSet m_reactions; ///< Set of REACLIB reactions in the network. std::unordered_map m_reactionIDMap; ///< Map from reaction ID to REACLIBReaction. //PERF: This makes copies of REACLIBReaction and could be a performance bottleneck. std::vector m_networkSpecies; ///< Vector of unique species in the network. std::unordered_map m_networkSpeciesMap; ///< Map from species name to Species object. std::unordered_map m_speciesToIndexMap; ///< Map from species to their index in the stoichiometry matrix. + 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). @@ -680,10 +690,11 @@ namespace gridfire { [[nodiscard]] StepDerivatives calculateAllDerivativesUsingPrecomputation( - const std::vector &Y_in, + const fourdst::composition::Composition &comp, const std::vector& bare_rates, const std::vector &bare_reverse_rates, - double T9, double rho + double T9, + double rho ) const; /** @@ -691,9 +702,11 @@ namespace gridfire { * * @tparam T The numeric type to use for the calculation. * @param reaction The reaction for which to calculate the flow. - * @param Y Vector of current abundances. + * @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. + * @param Ye + * @param mue * @return Molar flow rate for the reaction (e.g., mol/g/s). * * This method computes the net rate at which the given reaction proceeds @@ -702,17 +715,17 @@ namespace gridfire { template T calculateMolarReactionFlow( const reaction::Reaction &reaction, - const std::vector &Y, + const std::vector& Y, const T T9, - const T rho + const T rho, T Ye, T mue ) const; template T calculateReverseMolarReactionFlow( - T T9, - T rho, + const T T9, + const T rho, std::vector screeningFactors, - std::vector Y, + const std::vector& Y, size_t reactionIndex, const reaction::Reaction &reaction ) const; @@ -721,9 +734,11 @@ namespace gridfire { * @brief Calculates all derivatives (dY/dt) and the energy generation rate. * * @tparam T The numeric type to use for the calculation. - * @param Y_in Vector of current abundances for all species. + * @param Y_in 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. + * @param Ye + * @param mue * @return StepDerivatives containing dY/dt and energy generation rate. * * This method calculates the time derivatives of all species and the @@ -731,46 +746,46 @@ namespace gridfire { */ template [[nodiscard]] StepDerivatives calculateAllDerivatives( - const std::vector &Y_in, + const std::vector& Y_in, T T9, - T rho + T rho, T Ye, T mue ) const; - /** - * @brief Calculates all derivatives (dY/dt) and the energy generation rate (double precision). - * - * @param Y_in Vector of current abundances for all species. - * @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_in, - const double T9, - const double rho - ) const; - - /** - * @brief Calculates all derivatives (dY/dt) and the energy generation rate (automatic differentiation). - * - * @param Y_in Vector of current abundances for all species. - * @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_in, - const ADDouble &T9, - const ADDouble &rho - ) 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; }; @@ -780,7 +795,7 @@ namespace gridfire { T T9, T rho, std::vector screeningFactors, - std::vector Y, + const std::vector& Y, size_t reactionIndex, const reaction::Reaction &reaction ) const { @@ -804,9 +819,15 @@ namespace gridfire { } else { return reverseMolarFlow; // If no atomic function is available, return zero } - } else { + } else { // The case where T is of type double // A,B If not calling with an AD type, calculate the reverse rate directly - reverseRateConstant = calculateReverseRate(reaction, T9, 0, {}); + std::vector symbols; + symbols.reserve(m_networkSpecies.size()); + for (const auto& species : m_networkSpecies) { + symbols.emplace_back(species.name()); + } + fourdst::composition::Composition comp(symbols, Y); + reverseRateConstant = calculateReverseRate(reaction, T9, rho, comp); } // C. Get product multiplicities @@ -824,7 +845,7 @@ namespace gridfire { // E. Calculate the abundance term T productAbundanceTerm = static_cast(1.0); for (const auto& [species, count] : productCounts) { - const unsigned long speciesIndex = m_speciesToIndexMap.at(species); + const size_t speciesIndex = m_speciesToIndexMap.at(species); productAbundanceTerm *= CppAD::pow(Y[speciesIndex], count); } @@ -844,7 +865,12 @@ namespace gridfire { template StepDerivatives GraphEngine::calculateAllDerivatives( - const std::vector &Y_in, T T9, T rho) const { + const std::vector& Y_in, + const T T9, + const T rho, + const T Ye, + const T mue + ) const { std::vector screeningFactors = m_screeningModel->calculateScreeningFactors( m_reactions, m_networkSpecies, @@ -854,8 +880,10 @@ namespace gridfire { ); // --- Setup output derivatives structure --- - StepDerivatives result; - result.dydt.resize(m_networkSpecies.size(), static_cast(0.0)); + StepDerivatives result{}; + for (const auto& species : m_networkSpecies) { + result.dydt[species] = static_cast(0.0); // default the change in abundance to zero + } // --- AD Pre-setup (flags to control conditionals in an AD safe / branch aware manner) --- // ----- Constants for AD safe calculations --- @@ -885,13 +913,16 @@ namespace gridfire { const T N_A = static_cast(m_constants.Na); // Avogadro's number in mol^-1 const T c = static_cast(m_constants.c); // Speed of light in cm/s + // TODO: It may be prudent to introduce assertions here which validate units but will be removed in release builds (to ensure that unit inconsistencies do not creep in during future development) + // libconstants already has units built in so this should be straightforward. + // --- SINGLE LOOP OVER ALL REACTIONS --- for (size_t reactionIndex = 0; reactionIndex < m_reactions.size(); ++reactionIndex) { const auto& reaction = m_reactions[reactionIndex]; // 1. Calculate forward reaction rate const T forwardMolarReactionFlow = screeningFactors[reactionIndex] * - calculateMolarReactionFlow(reaction, Y, T9, rho); + calculateMolarReactionFlow(reaction, Y, T9, rho, Ye, mue); // 2. Calculate reverse reaction rate T reverseMolarFlow = static_cast(0.0); @@ -910,15 +941,15 @@ namespace gridfire { const T molarReactionFlow = forwardMolarReactionFlow - reverseMolarFlow; // Net molar reaction flow // 3. Use the rate to update all relevant species derivatives (dY/dt) - for (size_t speciesIndex = 0; speciesIndex < m_networkSpecies.size(); ++speciesIndex) { - const T nu_ij = static_cast(m_stoichiometryMatrix(speciesIndex, reactionIndex)); - result.dydt[speciesIndex] += threshold_flag * nu_ij * molarReactionFlow; + for (const auto& species: m_networkSpecies) { + const T nu_ij = static_cast(reaction.stoichiometry(species)); + result.dydt[species] += threshold_flag * nu_ij * molarReactionFlow; } } T massProductionRate = static_cast(0.0); // [mol][s^-1] - for (const auto& [species, index] : m_speciesToIndexMap) { - massProductionRate += result.dydt[index] * species.mass() * u; + for (const auto &species: m_speciesToIndexMap | std::views::keys) { + massProductionRate += result.dydt.at(species) * species.mass() * u; } result.nuclearEnergyGenerationRate = -massProductionRate * N_A * c * c; // [cm^2][s^-3] = [erg][s^-1][g^-1] @@ -930,9 +961,11 @@ namespace gridfire { template T GraphEngine::calculateMolarReactionFlow( const reaction::Reaction &reaction, - const std::vector &Y, + const std::vector& Y, const T T9, - const T rho + const T rho, + const T Ye, + const T mue ) const { // --- Pre-setup (flags to control conditionals in an AD safe / branch aware manner) --- @@ -940,7 +973,7 @@ namespace gridfire { const T zero = static_cast(0.0); // --- Calculate the molar reaction rate (in units of [s^-1][cm^3(N-1)][mol^(1-N)] for N reactants) --- - const T k_reaction = reaction.calculate_rate(T9, rho, Y); + const T k_reaction = reaction.calculate_rate(T9, rho, Ye, mue, Y, m_indexToSpeciesMap); // --- Cound the number of each reactant species to account for species multiplicity --- std::unordered_map reactant_counts; diff --git a/src/include/gridfire/engine/procedures/construction.h b/src/include/gridfire/engine/procedures/construction.h index 6dd50a2d..2a942912 100644 --- a/src/include/gridfire/engine/procedures/construction.h +++ b/src/include/gridfire/engine/procedures/construction.h @@ -7,6 +7,8 @@ #include +#include "gridfire/reaction/weak/weak_interpolator.h" + namespace gridfire { /** @@ -21,6 +23,7 @@ namespace gridfire { * * @param composition Mapping of isotopic species to their mass fractions; species with positive * mass fraction seed the network. + * @param weakInterpolator * @param maxLayers Variant specifying either a predefined NetworkBuildDepth or a custom integer depth; * negative depth (Full) collects all reactions, zero is invalid. * @param reverse If true, collects reverse reactions (decays or back-reactions); if false, uses forward reactions. @@ -32,7 +35,7 @@ namespace gridfire { */ reaction::ReactionSet build_reaclib_nuclear_network( const fourdst::composition::Composition &composition, - BuildDepthType maxLayers = NetworkBuildDepth::Full, - bool reverse = false + const rates::weak::WeakRateInterpolator &weakInterpolator, + BuildDepthType maxLayers = NetworkBuildDepth::Full, bool reverse = false ); -} \ No newline at end of file +} diff --git a/src/include/gridfire/engine/procedures/priming.h b/src/include/gridfire/engine/procedures/priming.h index 42e790dd..19a686da 100644 --- a/src/include/gridfire/engine/procedures/priming.h +++ b/src/include/gridfire/engine/procedures/priming.h @@ -37,7 +37,7 @@ namespace gridfire { * * @param engine Engine providing the current set of network reactions and flow calculations. * @param species The atomic species whose destruction rate is computed. - * @param Y Vector of molar abundances for all species in the engine. + * @param composition Current composition providing abundances for all species. * @param T9 Temperature in units of 10^9 K. * @param rho Density of the medium. * @pre Y.size() matches engine.getNetworkReactions().size() mapping species order. @@ -47,7 +47,7 @@ namespace gridfire { double calculateDestructionRateConstant( const DynamicEngine& engine, const fourdst::atomic::Species& species, - const std::vector& Y, + const fourdst::composition::Composition& composition, double T9, double rho ); @@ -60,7 +60,7 @@ namespace gridfire { * * @param engine Engine providing the current set of network reactions and flow calculations. * @param species The atomic species whose creation rate is computed. - * @param Y Vector of molar abundances for all species in the engine. + * @param composition Composition object containing current abundances. * @param T9 Temperature in units of 10^9 K. * @param rho Density of the medium. * @pre Y.size() matches engine.getNetworkReactions().size() mapping species order. @@ -70,7 +70,7 @@ namespace gridfire { double calculateCreationRate( const DynamicEngine& engine, const fourdst::atomic::Species& species, - const std::vector& Y, + const fourdst::composition::Composition& composition, double T9, double rho ); diff --git a/src/include/gridfire/engine/views/engine_adaptive.h b/src/include/gridfire/engine/views/engine_adaptive.h index fd7bf4bd..a65e9808 100644 --- a/src/include/gridfire/engine/views/engine_adaptive.h +++ b/src/include/gridfire/engine/views/engine_adaptive.h @@ -88,7 +88,7 @@ namespace gridfire { /** * @brief Calculates the right-hand side (dY/dt) and energy generation for the active species. * - * @param Y_culled A vector of abundances for the active species. + * @param comp The current composition of the system. * @param T9 The temperature in units of 10^9 K. * @param rho The density in g/cm^3. * @return A StepDerivatives struct containing the derivatives of the active species and the @@ -102,21 +102,29 @@ namespace gridfire { * @see AdaptiveEngineView::update() */ [[nodiscard]] std::expected, expectations::StaleEngineError> calculateRHSAndEnergy( - const std::vector &Y_culled, - const double T9, - const double rho + const fourdst::composition::Composition& comp, + double T9, + double rho ) const override; + + /** + * + * @param comp The current composition of the system. + * @param T9 The temperature in units of 10^9 K. + * @param rho The density in g/cm^3. + * @return A struct containing the derivatives of the energy generation rate with respect to temperature and density. + */ [[nodiscard]] EnergyDerivatives calculateEpsDerivatives( - const std::vector &Y_culled, - const double T9, - const double rho + const fourdst::composition::Composition& comp, + double T9, + double rho ) const override; /** * @brief Generates the Jacobian matrix for the active species. * - * @param Y_dynamic A vector of abundances for the active species. + * @param comp The current composition of the system. * @param T9 The temperature in units of 10^9 K. * @param rho The density in g/cm^3. * @@ -127,16 +135,16 @@ namespace gridfire { * @see AdaptiveEngineView::update() */ void generateJacobianMatrix( - const std::vector &Y_dynamic, - const double T9, - const double rho + const fourdst::composition::Composition &comp, + double T9, + double rho ) const override; /** * @brief Gets an entry from the Jacobian matrix for the active species. * - * @param i_culled The row index (species index) in the culled matrix. - * @param j_culled The column index (species index) in the culled matrix. + * @param rowSpecies The species corresponding to the row index in the culled species list. + * @param colSpecies The species corresponding to the column index in the culled species list * @return The value of the Jacobian matrix at (i_culled, j_culled). * * This method maps the culled indices to the full network indices and calls the base engine @@ -147,8 +155,8 @@ namespace gridfire { * @see AdaptiveEngineView::update() */ [[nodiscard]] double getJacobianMatrixEntry( - const int i_culled, - const int j_culled + const fourdst::atomic::Species& rowSpecies, + const fourdst::atomic::Species& colSpecies ) const override; /** @@ -165,8 +173,8 @@ namespace gridfire { /** * @brief Gets an entry from the stoichiometry matrix for the active species and reactions. * - * @param speciesIndex_culled The index of the species in the culled species list. - * @param reactionIndex_culled The index of the reaction in the culled reaction list. + * @param species The species for which to get the stoichiometric coefficient. + * @param reaction The reaction for which to get the stoichiometric coefficient. * @return The stoichiometric coefficient for the given species and reaction. * * This method maps the culled indices to the full network indices and calls the base engine @@ -177,15 +185,15 @@ namespace gridfire { * @see AdaptiveEngineView::update() */ [[nodiscard]] int getStoichiometryMatrixEntry( - const int speciesIndex_culled, - const int reactionIndex_culled + const fourdst::atomic::Species& species, + const reaction::Reaction& reaction ) const override; /** * @brief Calculates the molar reaction flow for a given reaction in the active network. * * @param reaction The reaction for which to calculate the flow. - * @param Y_culled Vector of current abundances for the active species. + * @param comp Composition object containing current abundances. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return Molar flow rate for the reaction (e.g., mol/g/s). @@ -198,7 +206,7 @@ namespace gridfire { */ [[nodiscard]] double calculateMolarReactionFlow( const reaction::Reaction &reaction, - const std::vector &Y_culled, + const fourdst::composition::Composition &comp, double T9, double rho ) const override; @@ -215,7 +223,7 @@ namespace gridfire { /** * @brief Computes timescales for all active species in the network. * - * @param Y_culled Vector of current abundances for the active species. + * @param comp Composition object containing current abundances. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return Map from Species to their characteristic timescales (s). @@ -226,13 +234,13 @@ namespace gridfire { * @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called). */ [[nodiscard]] std::expected, expectations::StaleEngineError> getSpeciesTimescales( - const std::vector &Y_culled, + const fourdst::composition::Composition &comp, double T9, double rho ) const override; [[nodiscard]] std::expected, expectations::StaleEngineError> getSpeciesDestructionTimescales( - const std::vector &Y, + const fourdst::composition::Composition &comp, double T9, double rho ) const override; @@ -391,21 +399,18 @@ namespace gridfire { * full network. * * @param netIn The current network input, containing temperature, density, and composition. - * @param out_Y_Full A vector that will be populated with the molar abundances of all species in the full network. - * @return A vector of ReactionFlow structs, each containing a pointer to a reaction and its calculated flow rate. + * @return A pair with the first element a vector of ReactionFlow structs, each containing a pointer to a reaction and its calculated flow rate and the second being a composition object where species which were not present in netIn but are present in the definition of the base engine are registered but have 0 mass fraction. * * @par Algorithm: - * 1. Clears and reserves space in `out_Y_Full`. - * 2. Iterates through all species in the base engine's network. - * 3. For each species, it retrieves the molar abundance from `netIn.composition`. If the species is not found, its abundance is set to 0.0. - * 4. Converts the temperature from Kelvin to T9. - * 5. Iterates through all reactions in the base engine's network. - * 6. For each reaction, it calls the base engine's `calculateMolarReactionFlow` to get the flow rate. - * 7. Stores the reaction pointer and its flow rate in a `ReactionFlow` struct and adds it to the returned vector. + * 1. Iterates through all species in the base engine's network. + * 2. For each species, it retrieves the molar abundance from `netIn.composition`. If the species is not found, its abundance is set to 0.0. + * 3. Converts the temperature from Kelvin to T9. + * 4. Iterates through all reactions in the base engine's network. + * 5. For each reaction, it calls the base engine's `calculateMolarReactionFlow` to get the flow rate. + * 6. Stores the reaction pointer and its flow rate in a `ReactionFlow` struct and adds it to the returned vector. */ - std::vector calculateAllReactionFlows( - const NetIn& netIn, - std::vector& out_Y_Full + std::pair, fourdst::composition::Composition> calculateAllReactionFlows( + const NetIn& netIn ) const; /** * @brief Finds all species that are reachable from the initial fuel through the reaction network. @@ -436,7 +441,7 @@ namespace gridfire { * * @param allFlows A vector of all reactions and their flow rates. * @param reachableSpecies A set of all species reachable from the initial fuel. - * @param Y_full A vector of molar abundances for all species in the full network. + * @param comp The current composition of the system. * @param maxFlow The maximum reaction flow rate in the network. * @return A vector of pointers to the reactions that have been kept after culling. * @@ -450,15 +455,15 @@ namespace gridfire { [[nodiscard]] std::vector cullReactionsByFlow( const std::vector& allFlows, const std::unordered_set& reachableSpecies, - const std::vector& Y_full, + const fourdst::composition::Composition& comp, double maxFlow ) const; typedef std::pair, std::unordered_set> RescueSet; [[nodiscard]] RescueSet rescueEdgeSpeciesDestructionChannel( - const std::vector& Y_full, - const double T9, - const double rho, + const fourdst::composition::Composition& comp, + double T9, + double rho, const std::vector& activeSpecies, const reaction::ReactionSet& activeReactions ) const; diff --git a/src/include/gridfire/engine/views/engine_defined.h b/src/include/gridfire/engine/views/engine_defined.h index 92de79c4..af788cb7 100644 --- a/src/include/gridfire/engine/views/engine_defined.h +++ b/src/include/gridfire/engine/views/engine_defined.h @@ -29,7 +29,7 @@ namespace gridfire{ /** * @brief Calculates the right-hand side (dY/dt) and energy generation for the active species. * - * @param Y_defined A vector of abundances for the active species. + * @param comp A Composition object containing the current composition of the system * @param T9 The temperature in units of 10^9 K. * @param rho The density in g/cm^3. * @return A StepDerivatives struct containing the derivatives of the active species and the @@ -38,44 +38,44 @@ namespace gridfire{ * @throws std::runtime_error If the view is stale (i.e., `update()` has not been called after `setNetworkFile()`). */ [[nodiscard]] std::expected, expectations::StaleEngineError> calculateRHSAndEnergy( - const std::vector& Y_defined, - const double T9, - const double rho + const fourdst::composition::Composition& comp, + double T9, + double rho ) const override; [[nodiscard]] EnergyDerivatives calculateEpsDerivatives( - const std::vector &Y, - const double T9, - const double rho + const fourdst::composition::Composition& comp, + double T9, + double rho ) const override; /** * @brief Generates the Jacobian matrix for the active species. * - * @param Y_dynamic A vector of abundances for the active species. + * @param comp A Composition object containing the current composition of the system * @param T9 The temperature in units of 10^9 K. * @param rho The density in g/cm^3. * * @throws std::runtime_error If the view is stale. */ void generateJacobianMatrix( - const std::vector& Y_dynamic, + const fourdst::composition::Composition& comp, const double T9, const double rho ) const override; /** * @brief Gets an entry from the Jacobian matrix for the active species. * - * @param i_defined The row index (species index) in the defined matrix. - * @param j_defined The column index (species index) in the defined matrix. - * @return The value of the Jacobian matrix at (i_defined, j_defined). + * @param rowSpecies The species corresponding to the row index. + * @param colSpecies The species corresponding to the column index. + * @return The value of the Jacobian matrix at (row species index, col species index). * * @throws std::runtime_error If the view is stale. * @throws std::out_of_range If an index is out of bounds. */ [[nodiscard]] double getJacobianMatrixEntry( - const int i_defined, - const int j_defined + const fourdst::atomic::Species& rowSpecies, + const fourdst::atomic::Species& colSpecies ) const override; /** * @brief Generates the stoichiometry matrix for the active reactions and species. @@ -86,22 +86,22 @@ namespace gridfire{ /** * @brief Gets an entry from the stoichiometry matrix for the active species and reactions. * - * @param speciesIndex_defined The index of the species in the defined species list. - * @param reactionIndex_defined The index of the reaction in the defined reaction list. + * @param species The species for which to get the stoichiometric coefficient. + * @param reaction The reaction for which to get the stoichiometric coefficient. * @return The stoichiometric coefficient for the given species and reaction. * * @throws std::runtime_error If the view is stale. * @throws std::out_of_range If an index is out of bounds. */ [[nodiscard]] int getStoichiometryMatrixEntry( - const int speciesIndex_defined, - const int reactionIndex_defined + const fourdst::atomic::Species& species, + const reaction::Reaction& reaction ) const override; /** * @brief Calculates the molar reaction flow for a given reaction in the active network. * * @param reaction The reaction for which to calculate the flow. - * @param Y_defined Vector of current abundances for the active species. + * @param comp A Composition object containing the current composition of the system * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return Molar flow rate for the reaction (e.g., mol/g/s). @@ -110,9 +110,9 @@ namespace gridfire{ */ [[nodiscard]] double calculateMolarReactionFlow( const reaction::Reaction& reaction, - const std::vector& Y_defined, - const double T9, - const double rho + const fourdst::composition::Composition& comp, + double T9, + double rho ) const override; /** * @brief Gets the set of active logical reactions in the network. @@ -127,7 +127,7 @@ namespace gridfire{ /** * @brief Computes timescales for all active species in the network. * - * @param Y_defined Vector of current abundances for the active species. + * @param comp A Composition object containing the current composition of the system * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return Map from Species to their characteristic timescales (s). @@ -135,15 +135,15 @@ namespace gridfire{ * @throws std::runtime_error If the view is stale. */ [[nodiscard]] std::expected, expectations::StaleEngineError> getSpeciesTimescales( - const std::vector& Y_defined, - const double T9, - const double rho + const fourdst::composition::Composition& comp, + double T9, + double rho ) const override; [[nodiscard]] std::expected, expectations::StaleEngineError> getSpeciesDestructionTimescales( - const std::vector& Y_defined, - const double T9, - const double rho + const fourdst::composition::Composition& comp, + double T9, + double rho ) const override; /** diff --git a/src/include/gridfire/engine/views/engine_multiscale.h b/src/include/gridfire/engine/views/engine_multiscale.h index 04179302..7c949660 100644 --- a/src/include/gridfire/engine/views/engine_multiscale.h +++ b/src/include/gridfire/engine/views/engine_multiscale.h @@ -206,7 +206,7 @@ namespace gridfire { /** * @brief Calculates the right-hand side (dY/dt) and energy generation. * - * @param Y_full Vector of current molar abundances for all species in the base engine. + * @param comp The current composition. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return A `std::expected` containing `StepDerivatives` on success, or a @@ -231,21 +231,21 @@ namespace gridfire { * (T9, rho, Y_full). This indicates `update()` was not called recently enough. */ [[nodiscard]] std::expected, expectations::StaleEngineError> calculateRHSAndEnergy( - const std::vector &Y_full, + const fourdst::composition::Composition& comp, double T9, double rho ) const override; [[nodiscard]] EnergyDerivatives calculateEpsDerivatives( - const std::vector &Y, - const double T9, - const double rho + const fourdst::composition::Composition& comp, + double T9, + double rho ) const override; /** * @brief Generates the Jacobian matrix for the current state. * - * @param Y_full Vector of current molar abundances. + * @param comp The current composition. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @@ -266,7 +266,7 @@ namespace gridfire { * without a valid partition. */ void generateJacobianMatrix( - const std::vector &Y_full, + const fourdst::composition::Composition& comp, double T9, double rho ) const override; @@ -274,8 +274,8 @@ namespace gridfire { /** * @brief Gets an entry from the previously generated Jacobian matrix. * - * @param i_full Row index (species index) in the full network. - * @param j_full Column index (species index) in the full network. + * @param rowSpecies Species corresponding to the row index (i_full). + * @param colSpecies Species corresponding to the column index (j_full). * @return Value of the Jacobian matrix at (i_full, j_full). * * @par Purpose @@ -289,8 +289,8 @@ namespace gridfire { * @pre `generateJacobianMatrix()` must have been called for the current state. */ [[nodiscard]] double getJacobianMatrixEntry( - int i_full, - int j_full + const fourdst::atomic::Species& rowSpecies, + const fourdst::atomic::Species& colSpecies ) const override; /** @@ -308,8 +308,8 @@ namespace gridfire { /** * @brief Gets an entry from the stoichiometry matrix. * - * @param speciesIndex Index of the species in the full network. - * @param reactionIndex Index of the reaction in the full network. + * @param species Species to look up stoichiometry for. + * @param reaction Reaction to find. * @return Stoichiometric coefficient for the species in the reaction. * * @par Purpose @@ -321,15 +321,15 @@ namespace gridfire { * @pre `generateStoichiometryMatrix()` must have been called. */ [[nodiscard]] int getStoichiometryMatrixEntry( - int speciesIndex, - int reactionIndex + const fourdst::atomic::Species& species, + const reaction::Reaction& reaction ) const override; /** * @brief Calculates the molar reaction flow for a given reaction. * * @param reaction The reaction for which to calculate the flow. - * @param Y_full Vector of current molar abundances for the full network. + * @param comp The current composition. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return Molar flow rate for the reaction (e.g., mol/g/s). @@ -349,7 +349,7 @@ namespace gridfire { */ [[nodiscard]] double calculateMolarReactionFlow( const reaction::Reaction &reaction, - const std::vector &Y_full, + const fourdst::composition::Composition& comp, double T9, double rho ) const override; @@ -385,7 +385,7 @@ namespace gridfire { /** * @brief Computes timescales for all species in the network. * - * @param Y Vector of current molar abundances for the full network. + * @param comp The current composition. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return A `std::expected` containing a map from `Species` to their characteristic @@ -403,7 +403,7 @@ namespace gridfire { * @throws StaleEngineError If the QSE cache misses. */ [[nodiscard]] std::expected, expectations::StaleEngineError> getSpeciesTimescales( - const std::vector &Y, + const fourdst::composition::Composition& comp, double T9, double rho ) const override; @@ -411,7 +411,7 @@ namespace gridfire { /** * @brief Computes destruction timescales for all species in the network. * - * @param Y Vector of current molar abundances for the full network. + * @param comp The current composition. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return A `std::expected` containing a map from `Species` to their characteristic @@ -429,7 +429,7 @@ namespace gridfire { * @throws StaleEngineError If the QSE cache misses. */ [[nodiscard]] std::expected, expectations::StaleEngineError> getSpeciesDestructionTimescales( - const std::vector &Y, + const fourdst::composition::Composition& comp, double T9, double rho ) const override; @@ -516,7 +516,7 @@ namespace gridfire { * * @param timescale_pools A vector of vectors of species indices, where each inner vector * represents a timescale pool. - * @param Y Vector of current molar abundances for the full network. + * @param comp Vector of current molar abundances for the full network. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return A vector of vectors of species indices, where each inner vector represents a @@ -531,9 +531,9 @@ namespace gridfire { * It then finds the connected components within that graph using a Breadth-First Search (BFS). * The resulting components from all pools are collected and returned. */ - std::vector> analyzeTimescalePoolConnectivity( - const std::vector> ×cale_pools, - const std::vector &Y, + std::vector> analyzeTimescalePoolConnectivity( + const std::vector> ×cale_pools, + const fourdst::composition::Composition &comp, double T9, double rho ) const; @@ -541,7 +541,7 @@ namespace gridfire { /** * @brief Partitions the network into dynamic and algebraic (QSE) groups based on timescales. * - * @param Y Vector of current molar abundances for the full network. + * @param comp Vector of current molar abundances for the full network. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @@ -567,7 +567,7 @@ namespace gridfire { * partitioning. */ void partitionNetwork( - const std::vector& Y, + const fourdst::composition::Composition &comp, double T9, double rho ); @@ -605,9 +605,9 @@ namespace gridfire { */ void exportToDot( const std::string& filename, - const std::vector& Y, - const double T9, - const double rho + const fourdst::composition::Composition &Y, + double T9, + double rho ) const; /** @@ -679,7 +679,7 @@ namespace gridfire { /** * @brief Equilibrates the network by partitioning and solving for QSE abundances. * - * @param Y Vector of current molar abundances for the full network. + * @param comp Vector of current molar abundances for the full network. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return A new composition object with the equilibrated abundances. @@ -698,7 +698,7 @@ namespace gridfire { * @post The engine's internal partition is updated. A new composition object is returned. */ fourdst::composition::Composition equilibrateNetwork( - const std::vector &Y, + const fourdst::composition::Composition &comp, double T9, double rho ); @@ -730,11 +730,10 @@ namespace gridfire { * in quasi-steady-state equilibrium with each other. */ struct QSEGroup { - std::set species_indices; ///< Indices of all species in this group. - bool is_in_equilibrium = false; ///< Flag set by flux analysis. - std::set algebraic_indices; ///< Indices of algebraic species in this group. - std::set seed_indices; ///< Indices of dynamic species in this group. - double mean_timescale; ///< Mean timescale of the group. + bool is_in_equilibrium = false; ///< Flag set by flux analysis. + std::set algebraic_species; ///< Algebraic species in this group. + std::set seed_species; ///< Dynamic species in this group. + double mean_timescale; ///< Mean timescale of the group. /** * @brief Less-than operator for QSEGroup, used for sorting. @@ -761,9 +760,7 @@ namespace gridfire { */ bool operator!=(const QSEGroup& other) const; - std::string toString() const; - - std::string toString(const DynamicEngine &engine) const; + [[nodiscard]] [[maybe_unused]] std::string toString() const; }; /** @@ -802,11 +799,11 @@ namespace gridfire { /** * @brief Indices of the species to solve for in the QSE group. */ - const std::vector& m_qse_solve_indices; + const std::set& m_qse_solve_species; /** * @brief Initial abundances of all species in the full network. */ - const std::vector& m_Y_full_initial; + const fourdst::composition::Composition& m_initial_comp; /** * @brief Temperature in units of 10^9 K. */ @@ -820,41 +817,49 @@ namespace gridfire { */ const Eigen::VectorXd& m_Y_scale; + /** + * @brief Mapping from species to their indices in the QSE solve vector. + */ + const std::unordered_map m_qse_solve_species_index_map; + /** * @brief Constructs an EigenFunctor. * * @param view The MultiscalePartitioningEngineView instance. - * @param qse_solve_indices Indices of the species to solve for in the QSE group. - * @param Y_full_initial Initial abundances of all species. + * @param qse_solve_species Species to solve for in the QSE group. + * @param initial_comp Initial abundances of all species in the full network. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @param Y_scale Scaling factors for the species abundances. + * @param qse_solve_species_index_map Mapping from species to their indices in the QSE solve vector. */ EigenFunctor( MultiscalePartitioningEngineView& view, - const std::vector& qse_solve_indices, - const std::vector& Y_full_initial, + const std::set& qse_solve_species, + const fourdst::composition::Composition& initial_comp, const double T9, const double rho, - const Eigen::VectorXd& Y_scale + const Eigen::VectorXd& Y_scale, + const std::unordered_map& qse_solve_species_index_map ) : m_view(&view), - m_qse_solve_indices(qse_solve_indices), - m_Y_full_initial(Y_full_initial), + m_qse_solve_species(qse_solve_species), + m_initial_comp(initial_comp), m_T9(T9), m_rho(rho), - m_Y_scale(Y_scale) {} + m_Y_scale(Y_scale), + m_qse_solve_species_index_map(qse_solve_species_index_map){} /** * @brief Gets the number of output values from the functor (size of the residual vector). * @return The number of algebraic species being solved. */ - [[nodiscard]] size_t values() const { return m_qse_solve_indices.size(); } + [[nodiscard]] size_t values() const { return m_qse_solve_species.size(); } /** * @brief Gets the number of input values to the functor (size of the variable vector). * @return The number of algebraic species being solved. */ - [[nodiscard]] size_t inputs() const { return m_qse_solve_indices.size(); } + [[nodiscard]] size_t inputs() const { return m_qse_solve_species.size(); } /** * @brief Evaluates the functor's residual vector `f_qse = dY_alg/dt`. @@ -981,23 +986,15 @@ namespace gridfire { * @brief The simplified set of species presented to the solver (the "slow" species). */ std::vector m_dynamic_species; - /** - * @brief Indices mapping the dynamic species back to the base engine's full species list. - */ - std::vector m_dynamic_species_indices; /** * @brief Species that are treated as algebraic (in QSE) in the QSE groups. */ std::vector m_algebraic_species; /** - * @breif Stateful storage of the current algebraic species abundances. This is updated every time the update method is called. + * @brief Map from species to their calculated abundances in the QSE state. */ - std::vector m_Y_algebraic; - /** - * @brief Indices of algebraic species in the full network. - */ - std::vector m_algebraic_species_indices; + std::unordered_map m_algebraic_abundances; /** * @brief Indices of all species considered active in the current partition (dynamic + algebraic). @@ -1029,7 +1026,7 @@ namespace gridfire { /** * @brief Partitions the network by timescale. * - * @param Y_full Vector of current molar abundances for all species. + * @param comp Vector of current molar abundances for all species. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return A vector of vectors of species indices, where each inner vector represents a @@ -1044,8 +1041,8 @@ namespace gridfire { * a gap between consecutive timescales that is larger than a predefined threshold * (e.g., a factor of 100). */ - std::vector> partitionByTimescale( - const std::vector &Y_full, + std::vector> partitionByTimescale( + const fourdst::composition::Composition &comp, double T9, double rho ) const; @@ -1054,7 +1051,7 @@ namespace gridfire { * @brief Validates candidate QSE groups using flux analysis. * * @param candidate_groups A vector of candidate QSE groups. - * @param Y Vector of current molar abundances for the full network. + * @param comp Vector of current molar abundances for the full network. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return A vector of validated QSE groups that meet the flux criteria. @@ -1074,7 +1071,7 @@ namespace gridfire { :: QSEGroup>> validateGroupsWithFluxAnalysis( const std::vector &candidate_groups, - const std::vector &Y, + const fourdst::composition::Composition &comp, double T9, double rho ) const; @@ -1082,7 +1079,7 @@ namespace gridfire { /** * @brief Solves for the QSE abundances of the algebraic species in a given state. * - * @param Y_full Vector of current molar abundances for all species in the base engine. + * @param comp Vector of current molar abundances for all species in the base engine. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return A vector of molar abundances for the algebraic species. @@ -1099,8 +1096,8 @@ namespace gridfire { * @pre The input state (Y_full, T9, rho) must be a valid physical state. * @post The algebraic species in the QSE cache are updated with the new equilibrium abundances. */ - std::vector solveQSEAbundances( - const std::vector &Y_full, + fourdst::composition::Composition solveQSEAbundances( + const fourdst::composition::Composition &comp, double T9, double rho ); @@ -1110,7 +1107,7 @@ namespace gridfire { * * @param pools A vector of vectors of species indices, where each inner vector represents a * timescale pool. - * @param Y Vector of current molar abundances for the full network. + * @param comp Vector of current molar abundances for the full network. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return The index of the pool with the largest (slowest) mean destruction timescale. @@ -1123,8 +1120,8 @@ namespace gridfire { * pool and returns the index of the pool with the maximum mean timescale. */ size_t identifyMeanSlowestPool( - const std::vector>& pools, - const std::vector &Y, + const std::vector>& pools, + const fourdst::composition::Composition &comp, double T9, double rho ) const; @@ -1144,8 +1141,8 @@ namespace gridfire { * and one as a product), it adds edges between all reactants and products from * that reaction that are also in the pool. */ - std::unordered_map> buildConnectivityGraph( - const std::vector& species_pool + std::unordered_map> buildConnectivityGraph( + const std::vector& species_pool ) const; /** @@ -1153,7 +1150,7 @@ namespace gridfire { * * @param candidate_pools A vector of vectors of species indices, where each inner vector * represents a connected pool of species with similar fast timescales. - * @param Y Vector of current molar abundances. + * @param comp Vector of current molar abundances. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @return A vector of `QSEGroup` structs, ready for flux validation. @@ -1168,8 +1165,8 @@ namespace gridfire { * @post A list of candidate `QSEGroup` objects is returned. */ std::vector constructCandidateGroups( - const std::vector>& candidate_pools, - const std::vector& Y, + const std::vector>& candidate_pools, + const fourdst::composition::Composition &comp, double T9, double rho ) const; diff --git a/src/include/gridfire/exceptions/error_engine.h b/src/include/gridfire/exceptions/error_engine.h index 46beddb4..b7516b0f 100644 --- a/src/include/gridfire/exceptions/error_engine.h +++ b/src/include/gridfire/exceptions/error_engine.h @@ -2,7 +2,7 @@ #include #include -#include +#include namespace gridfire::exceptions { class EngineError : public std::exception {}; diff --git a/src/include/gridfire/expectations/expected_engine.h b/src/include/gridfire/expectations/expected_engine.h index 3d885ced..6c11daf6 100644 --- a/src/include/gridfire/expectations/expected_engine.h +++ b/src/include/gridfire/expectations/expected_engine.h @@ -2,6 +2,7 @@ #include #include +#include namespace gridfire::expectations { enum class EngineErrorTypes { @@ -19,8 +20,8 @@ namespace gridfire::expectations { std::string m_message; const EngineErrorTypes type = EngineErrorTypes::FAILURE; - explicit EngineError(const std::string &message, const EngineErrorTypes type) - : m_message(message), type(type) {} + explicit EngineError(std::string message, const EngineErrorTypes type) + : m_message(std::move(message)), type(type) {} virtual ~EngineError() = default; diff --git a/src/include/gridfire/io/network_file.h b/src/include/gridfire/io/network_file.h index c9dd276d..6f25e312 100644 --- a/src/include/gridfire/io/network_file.h +++ b/src/include/gridfire/io/network_file.h @@ -99,7 +99,7 @@ namespace gridfire::io { * ParsedNetworkData data = parser.parse("reactions.txt"); * @endcode */ - ParsedNetworkData parse(const std::string& filename) const override; + [[nodiscard]] ParsedNetworkData parse(const std::string& filename) const override; private: using Config = fourdst::config::Config; using LogManager = fourdst::logging::LogManager; @@ -139,7 +139,7 @@ namespace gridfire::io { * @throws std::runtime_error If the file cannot be opened or if it * contains formatting errors. */ - ParsedNetworkData parse(const std::string& filename) const override; + [[nodiscard]] ParsedNetworkData parse(const std::string& filename) const override; private: using Config = fourdst::config::Config; using LogManager = fourdst::logging::LogManager; diff --git a/src/include/gridfire/network.h b/src/include/gridfire/network.h index 036db3e2..62ca0662 100644 --- a/src/include/gridfire/network.h +++ b/src/include/gridfire/network.h @@ -89,7 +89,7 @@ namespace gridfire { */ virtual NetOut evaluate(const NetIn &netIn) = 0; - virtual bool isStiff() const { return m_stiff; } + [[nodiscard]] virtual bool isStiff() const { return m_stiff; } virtual void setStiff(const bool stiff) { m_stiff = stiff; } protected: diff --git a/src/include/gridfire/partition/composite/partition_composite.h b/src/include/gridfire/partition/composite/partition_composite.h index a6c33a24..585b775a 100644 --- a/src/include/gridfire/partition/composite/partition_composite.h +++ b/src/include/gridfire/partition/composite/partition_composite.h @@ -99,6 +99,6 @@ namespace gridfire::partition { * @return Unique pointer to a new PartitionFunction instance of the given type. * @throws std::runtime_error If the given type is not recognized. + */ - std::unique_ptr selectPartitionFunction(const BasePartitionType type) const; + [[nodiscard]] std::unique_ptr selectPartitionFunction(const BasePartitionType type) const; }; } \ No newline at end of file diff --git a/src/include/gridfire/reaction/reaction.h b/src/include/gridfire/reaction/reaction.h index 9d0d74b3..ab9a71a8 100644 --- a/src/include/gridfire/reaction/reaction.h +++ b/src/include/gridfire/reaction/reaction.h @@ -11,6 +11,7 @@ #include "cppad/cppad.hpp" +#include "fourdst/composition/composition.h" /** * @file reaction.h @@ -78,8 +79,18 @@ namespace gridfire::reaction { public: virtual ~Reaction() = default; - [[nodiscard]] virtual double calculate_rate(double T9, double rho, const std::vector& Y) const = 0; - [[nodiscard]] virtual CppAD::AD calculate_rate(CppAD::AD T9, CppAD::AD rho, const std::vector>& Y) const = 0; + [[nodiscard]] virtual double calculate_rate( + double T9, + double rho, + double Ye, + double mue, const std::vector &Y, const std::unordered_map& index_to_species_map + ) const = 0; + [[nodiscard]] virtual CppAD::AD calculate_rate( + CppAD::AD T9, + CppAD::AD rho, + CppAD::AD Ye, + CppAD::AD mue, const std::vector>& Y, const std::unordered_map& index_to_species_map + ) const = 0; [[nodiscard]] virtual std::string_view id() const = 0; [[nodiscard]] virtual const std::vector& reactants() const = 0; @@ -103,7 +114,26 @@ namespace gridfire::reaction { [[nodiscard]] virtual uint64_t hash(uint64_t seed) const = 0; [[nodiscard]] virtual double qValue() const = 0; - [[nodiscard]] virtual double calculate_forward_rate_log_derivative(double T9, double rho, const std::vector& Y) const = 0; + + [[nodiscard]] virtual double calculate_energy_generation_rate( + double T9, + double rho, + double Ye, + double mue, const std::vector& Y, const std::unordered_map& index_to_species_map + ) const { + return calculate_rate(T9, rho, 0, 0, Y, index_to_species_map) * qValue(); + } + + [[nodiscard]] virtual CppAD::AD calculate_energy_generation_rate( + const CppAD::AD& T9, + const CppAD::AD& rho, + const CppAD::AD &Ye, + const CppAD::AD &mue, const std::vector>& Y, const std::unordered_map& index_to_species_map + ) const { + return calculate_rate(T9, rho, {}, {}, Y, index_to_species_map) * qValue(); + } + + [[nodiscard]] virtual double calculate_forward_rate_log_derivative(double T9, double rho, double Ye, double mue, const fourdst::composition::Composition& comp) const = 0; [[nodiscard]] virtual ReactionType type() const = 0; @@ -145,21 +175,37 @@ namespace gridfire::reaction { * @brief Calculates the reaction rate for a given temperature. * @param T9 The temperature in units of 10^9 K. * @param rho Density [Not used in this implementation]. - * @param Y Molar abundances of species [Not used in this implementation]. + * @param Ye + * @param mue + * @param Y + * @param index_to_species_map * @return The calculated reaction rate. */ - [[nodiscard]] double calculate_rate(double T9, double rho, const std::vector& Y) const override; + [[nodiscard]] double calculate_rate( + double T9, + double rho, + double Ye, + double mue, const std::vector &Y, const std::unordered_map& index_to_species_map + ) const override; /** * @brief Calculates the reaction rate for a given temperature using CppAD types. * @param T9 The temperature in units of 10^9 K, as a CppAD::AD. * @param rho Density, as a CppAD::AD [Not used in this implementation]. + * @param Ye + * @param mue * @param Y Molar abundances of species, as a vector of CppAD::AD [Not used in this implementation]. + * @param index_to_species_map * @return The calculated reaction rate, as a CppAD::AD. */ - [[nodiscard]] CppAD::AD calculate_rate(CppAD::AD T9, CppAD::AD rho, const std::vector>& Y) const override; + [[nodiscard]] CppAD::AD calculate_rate( + CppAD::AD T9, + CppAD::AD rho, + CppAD::AD Ye, + CppAD::AD mue, const std::vector>& Y, const std::unordered_map& index_to_species_map + ) const override; - [[nodiscard]] double calculate_forward_rate_log_derivative(double T9, double rho, const std::vector& Y) const override; + [[nodiscard]] double calculate_forward_rate_log_derivative(double T9, double rho, double Ye, double mue, const fourdst::composition::Composition& comp) const override; /** * @brief Gets the reaction name in (projectile, ejectile) notation. @@ -389,12 +435,24 @@ namespace gridfire::reaction { * @brief Calculates the total reaction rate by summing all source rates. * @param T9 The temperature in units of 10^9 K. * @param rho + * @param Ye + * @param mue * @param Y + * @param index_to_species_map * @return The total calculated reaction rate. */ - [[nodiscard]] double calculate_rate(double T9, double rho, const std::vector& Y) const override; + [[nodiscard]] double calculate_rate( + double T9, + double rho, + double Ye, + double mue, const std::vector &Y, const std::unordered_map& index_to_species_map + ) const override; - [[nodiscard]] double calculate_forward_rate_log_derivative(double T9, double rho, const std::vector& Y) const override; + [[nodiscard]] double calculate_forward_rate_log_derivative( + double T9, + double rho, + double Ye, double mue, const fourdst::composition::Composition& comp + ) const override; [[nodiscard]] ReactionType type() const override { return ReactionType::LOGICAL_REACLIB; } @@ -404,10 +462,18 @@ namespace gridfire::reaction { * @brief Calculates the total reaction rate using CppAD types. * @param T9 The temperature in units of 10^9 K, as a CppAD::AD. * @param rho + * @param Ye + * @param mue * @param Y + * @param index_to_species_map * @return The total calculated reaction rate, as a CppAD::AD. */ - [[nodiscard]] CppAD::AD calculate_rate(CppAD::AD T9, CppAD::AD rho, const std::vector>& Y) const override; + [[nodiscard]] CppAD::AD calculate_rate( + CppAD::AD T9, + CppAD::AD rho, + CppAD::AD Ye, + CppAD::AD mue, const std::vector>& Y, const std::unordered_map& index_to_species_map + ) const override; /** @name Iterators * Provides iterators to loop over the rate coefficient sets. diff --git a/src/include/gridfire/reaction/weak/weak.h b/src/include/gridfire/reaction/weak/weak.h index 55ceee0a..fbffffd4 100644 --- a/src/include/gridfire/reaction/weak/weak.h +++ b/src/include/gridfire/reaction/weak/weak.h @@ -1,57 +1,254 @@ #pragma once -#include "fourdst/composition/atomicSpecies.h" +#define GRIDFIRE_WEAK_REACTION_LIB_SENTINEL -60.0 +#include "gridfire/reaction/reaction.h" +#include "gridfire/reaction/weak/weak_types.h" +#include "gridfire/reaction/weak/weak_interpolator.h" + +#include "gridfire/engine/engine_abstract.h" + +#include "fourdst/composition/atomicSpecies.h" +#include "fourdst/constants/const.h" + +#include "cppad/cppad.hpp" + +#include #include #include +#include +#include +#include +#include + + namespace gridfire::rates::weak { - enum class WeakReactionType { - BETA_PLUS_DECAY, - BETA_MINUS_DECAY, - ELECTRON_CAPTURE, - POSITRON_CAPTURE, - }; - - inline std::unordered_map WeakReactionTypeNames = { - {WeakReactionType::BETA_PLUS_DECAY, "β+ Decay"}, - {WeakReactionType::BETA_MINUS_DECAY, "β- Decay"}, - {WeakReactionType::ELECTRON_CAPTURE, "e- Capture"}, - {WeakReactionType::POSITRON_CAPTURE, "e+ Capture"}, - }; - - struct WeakReaction { - WeakReactionType type; - float T9; - float log_rhoYe; - float mu_e; - float log_rate; - float log_neutrino_loss; - - friend std::ostream& operator<<(std::ostream& os, const WeakReaction& reaction) { - os << "WeakReaction(type=" << WeakReactionTypeNames[reaction.type] - << ", T9=" << reaction.T9 - << ", log_rhoYe=" << reaction.log_rhoYe - << ", mu_e=" << reaction.mu_e - << ", log_rate=" << reaction.log_rate - << ", log_neutrino_loss=" << reaction.log_neutrino_loss - << ")"; - return os; - } - }; - class WeakReactionMap { public: WeakReactionMap(); ~WeakReactionMap() = default; - std::vector get_all_reactions() const; + [[nodiscard]] std::vector get_all_reactions() const; - std::expected, bool> get_species_reactions(const fourdst::atomic::Species &species) const; - std::expected, bool> get_species_reactions(const std::string& species_name) const; + [[nodiscard]] std::expected, WeakMapError> get_species_reactions( + const fourdst::atomic::Species &species) const; + + [[nodiscard]] std::expected, WeakMapError> get_species_reactions( + const std::string &species_name) const; private: - std::unordered_map> m_weak_network; + std::unordered_map> m_weak_network; }; + + class WeakReaction final : public reaction::Reaction { + public: + explicit WeakReaction( + const fourdst::atomic::Species &species, + WeakReactionType type, + const WeakRateInterpolator& interpolator + ); + [[nodiscard]] double calculate_rate( + double T9, + double rho, + double Ye, + double mue, + const std::vector &Y, + const std::unordered_map& index_to_species_map + ) const override; + [[nodiscard]] CppAD::AD calculate_rate( + CppAD::AD T9, + CppAD::AD rho, + CppAD::AD Ye, + CppAD::AD mue, + const std::vector> &Y, + const std::unordered_map& index_to_species_map + ) const override; + [[nodiscard]] std::string_view id() const override { return m_id; } + [[nodiscard]] const std::vector &reactants() const override { return {m_reactant}; } + [[nodiscard]] const std::vector &products() const override { return {m_product}; } + [[nodiscard]] bool contains(const fourdst::atomic::Species &species) const override; + [[nodiscard]] bool contains_reactant(const fourdst::atomic::Species &species) const override; + [[nodiscard]] bool contains_product(const fourdst::atomic::Species &species) const override; + [[nodiscard]] std::unordered_set all_species() const override; + [[nodiscard]] std::unordered_set reactant_species() const override; + [[nodiscard]] std::unordered_set product_species() const override; + [[nodiscard]] size_t num_species() const override { return 2; } // Always 2 for weak reactions + [[nodiscard]] std::unordered_map stoichiometry() const override; + [[nodiscard]] int stoichiometry(const fourdst::atomic::Species &species) const override; + [[nodiscard]] uint64_t hash(uint64_t seed) const override; + [[nodiscard]] double qValue() const override; + [[nodiscard]] double calculate_energy_generation_rate( + double T9, + double rho, + double Ye, + double mue, + const std::vector& Y, + const std::unordered_map& index_to_species_map + ) const override; + [[nodiscard]] CppAD::AD calculate_energy_generation_rate( + const CppAD::AD &T9, + const CppAD::AD &rho, + const CppAD::AD &Ye, + const CppAD::AD &mue, + const std::vector> &Y, + const std::unordered_map &index_to_species_map + ) const override; + [[nodiscard]] double calculate_forward_rate_log_derivative( + double T9, + double rho, + double Ye, + double mue, + const fourdst::composition::Composition& composition + ) const override; + [[nodiscard]] reaction::ReactionType type() const override { return reaction::ReactionType::WEAK; } + [[nodiscard]] std::unique_ptr clone() const override; + [[nodiscard]] bool is_reverse() const override { return false; }; + + private: + template + T calculate_rate( + T T9, + T rho, + T Ye, + T mue, + const std::vector &Y, + const std::unordered_map& index_to_species_map + ) const; + + double get_log_rate_from_payload(const WeakRatePayload& payload) const; + + private: + class AtomicWeakRate final : public CppAD::atomic_base { + public: + AtomicWeakRate( + const WeakRateInterpolator& interpolator, + const size_t a, + const size_t z, + const WeakReactionType type + ) : + CppAD::atomic_base(std::format("atomic-{}-{}-weak-rate", a, z)), + m_interpolator(interpolator), + m_a(a), + m_z(z) , + m_type(type) {} + + bool forward( + size_t p, + size_t q, + const CppAD::vector& vx, + CppAD::vector& vy, + const CppAD::vector& tx, + CppAD::vector& ty + ) override; + bool reverse( + size_t q, + const CppAD::vector& tx, + const CppAD::vector& ty, + CppAD::vector& px, + const CppAD::vector& py + ) override; + bool for_sparse_jac( + size_t q, + const CppAD::vector>&r, + CppAD::vector>& s + ) override; + bool rev_sparse_jac( + size_t q, + const CppAD::vector>&rt, + CppAD::vector>& st + ) override; + + private: + const WeakRateInterpolator& m_interpolator; + const size_t m_a; + const size_t m_z; + const WeakReactionType m_type; + }; + + struct constants { + const fourdst::constant::Constants& c = fourdst::constant::Constants::getInstance(); + fourdst::constant::Constant neutronMassG = c.get("mN"); + fourdst::constant::Constant protonMassG = c.get("mP"); + fourdst::constant::Constant electronMassG = c.get("mE"); + fourdst::constant::Constant speedOfLight = c.get("c"); + fourdst::constant::Constant eVgRelation = c.get("eV_kg"); // note that despite the symbol this is in g NOT kg + fourdst::constant::Constant MeV2Erg = c.get("MeV_to_erg"); + fourdst::constant::Constant amu = c.get("u"); + + double MeVgRelation = eVgRelation.value * 1.0e6; + double MeVPerGraph = 1.0/MeVgRelation; + + double neutronMassMeV = neutronMassG.value * MeVgRelation; + double protonMassMeV = protonMassG.value * MeVgRelation; + double electronMassMeV = electronMassG.value * MeVgRelation; + + double u_to_MeV = (amu.value * speedOfLight.value * speedOfLight.value)/MeV2Erg.value; + }; + + private: + const constants m_constants; + fourdst::atomic::Species m_reactant; + fourdst::atomic::Species m_product; + + size_t m_reactant_a; + size_t m_reactant_z; + size_t m_product_a; + size_t m_product_z; + + std::string m_id; + WeakReactionType m_type; + + const WeakRateInterpolator& m_interpolator; + + AtomicWeakRate m_atomic; + + }; + + template + T WeakReaction::calculate_rate( + T T9, + T rho, + T Ye, + T mue, + const std::vector &Y, + const std::unordered_map& index_to_species_map + ) const { + T log_rhoYe = CppAD::log10(rho * Ye); + + T rateConstant = static_cast(0.0); + if constexpr (std::is_same_v>) { // Case where T is an AD type + std::vector ax = {T9, log_rhoYe, mue}; + std::vector ay(1); + m_atomic(ax, ay); + rateConstant = static_cast(ay[0]); + } else { // The case where T is of type double + const std::expected result = m_interpolator.get_rates( + static_cast(m_reactant_a), + static_cast(m_reactant_z), + T9, + log_rhoYe, + mue + ); + + if (!result.has_value()) { + const InterpolationErrorType type = result.error().type; + const std::string msg = std::format( + "Failed to interpolate weak rate for (A={}, Z={}) at T9={}, log10(rho*Ye)={}, mu_e={} with error: {}", + m_reactant.name(), m_reactant_a, m_reactant_z, T9, log_rhoYe, mue, InterpolationErrorTypeMap.at(type) + ); + throw std::runtime_error(msg); + } + const WeakRatePayload payload = result.value(); + const double logRate = get_log_rate_from_payload(payload); + if (logRate <= GRIDFIRE_WEAK_REACTION_LIB_SENTINEL) { + rateConstant = static_cast(0.0); + } else { + rateConstant = CppAD::pow(10.0, logRate); + } + + } + return rateConstant; + } } diff --git a/src/include/gridfire/reaction/weak/weak_interpolator.h b/src/include/gridfire/reaction/weak/weak_interpolator.h new file mode 100644 index 00000000..92e9c797 --- /dev/null +++ b/src/include/gridfire/reaction/weak/weak_interpolator.h @@ -0,0 +1,44 @@ +#pragma once + +#include "gridfire/reaction/weak/weak_types.h" +#include "fourdst/composition/atomicSpecies.h" + +#include +#include +#include +#include + + + +namespace gridfire::rates::weak { + class WeakRateInterpolator { + public: + using RowDataTable = std::array; // Total number of entries in the weak rate table NOTE: THIS MUST EQUAL THE VALUE IN weak_rate_library.h + + explicit WeakRateInterpolator(const RowDataTable& raw_data); + + [[nodiscard]] std::vector available_isotopes() const; + + [[nodiscard]] std::expected get_rates( + uint16_t A, + uint8_t Z, + double t9, + double log_rhoYe, + double mu_e + ) const; + + [[nodiscard]] std::expected get_rate_derivatives( + uint16_t A, + uint8_t Z, + double t9, + double log_rhoYe, + double mu_e + ) const; + private: + static uint32_t pack_isotope_id(uint16_t A, uint8_t Z); + + std::unordered_map m_rate_table; + }; + + +} \ No newline at end of file diff --git a/src/include/gridfire/reaction/weak/weak_rate_library.h b/src/include/gridfire/reaction/weak/weak_rate_library.h index 6dc4b662..16f9010c 100644 --- a/src/include/gridfire/reaction/weak/weak_rate_library.h +++ b/src/include/gridfire/reaction/weak/weak_rate_library.h @@ -8,25 +8,10 @@ #include -#include +#include "gridfire/reaction/weak/weak.h" + namespace gridfire::rates::weak { - - // Represents a single row from the unified weak rate table - struct RateDataRow { - uint16_t A; - uint8_t Z; - float t9; - float log_rhoye; - float mu_e; - float log_beta_plus; - float log_electron_capture; - float log_neutrino_loss_ec; - float log_beta_minus; - float log_positron_capture; - float log_antineutrino_loss_bd; - }; - // The complete, pre-processed weak rate table data static constexpr std::array UNIFIED_WEAK_DATA = { RateDataRow(6, 3, 0.01, 5.0, 0.485, -60.0, -60.0, -60.0, -60.0, -60.0, -60.0), diff --git a/src/include/gridfire/reaction/weak/weak_types.h b/src/include/gridfire/reaction/weak/weak_types.h new file mode 100644 index 00000000..7f0ccd9e --- /dev/null +++ b/src/include/gridfire/reaction/weak/weak_types.h @@ -0,0 +1,133 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +namespace gridfire::rates::weak { + struct RateDataRow { + uint16_t A; + uint8_t Z; + float t9; + float log_rhoye; + float mu_e; + float log_beta_plus; + float log_electron_capture; + float log_neutrino_loss_ec; + float log_beta_minus; + float log_positron_capture; + float log_antineutrino_loss_bd; + }; + + enum class WeakReactionType { + BETA_PLUS_DECAY, + BETA_MINUS_DECAY, + ELECTRON_CAPTURE, + POSITRON_CAPTURE, + }; + + enum class NeutrinoTypes { + ELECTRON_NEUTRINO, + ELECTRON_ANTINEUTRINO, + MUON_NEUTRINO, + MUON_ANTINEUTRINO, + TAU_NEUTRINO, + TAU_ANTINEUTRINO + }; + + enum class WeakMapError { + SPECIES_NOT_FOUND, + UNKNOWN_ERROR + }; + + struct WeakRatePayload { + double log_beta_plus; + double log_electron_capture; + double log_neutrino_loss_ec; + double log_beta_minus; + double log_positron_capture; + double log_antineutrino_loss_bd; + }; + + struct WeakRateDerivatives { + // Each array holds [d/dT9, d/dlogRhoYe, d/dMuE] + std::array d_log_beta_plus; + std::array d_log_electron_capture; + std::array d_log_neutrino_loss_ec; + std::array d_log_beta_minus; + std::array d_log_positron_capture; + std::array d_log_antineutrino_loss_bd; + }; + + enum class InterpolationErrorType { + BOUNDS_ERROR, + UNKNOWN_SPECIES_ERROR, + UNKNOWN_ERROR + }; + + inline std::unordered_map InterpolationErrorTypeMap = { + {InterpolationErrorType::BOUNDS_ERROR, "Bounds Error"}, + {InterpolationErrorType::UNKNOWN_SPECIES_ERROR, "Unknown Species Error"}, + {InterpolationErrorType::UNKNOWN_ERROR, "Unknown Error"} + }; + + enum class TableAxes { + T9, + LOG_RHOYE, + MUE + }; + + struct BoundsErrorInfo { + TableAxes axis; + double axisMinValue; + double axisMaxValue; + double queryValue; + }; + + struct InterpolationError { + InterpolationErrorType type; + std::optional> boundsErrorInfo = std::nullopt; + }; + + struct IsotopeGrid { + std::vector t9_axis; + std::vector rhoYe_axis; + std::vector mue_axis; + + // index = (i_t9 * logRhoYe_axis.size() + j_rhoYe) + mue_axis.size() + k_mue + std::vector data; + }; + + constexpr std::string_view weak_reaction_type_name(const WeakReactionType t) noexcept { + switch (t) { + case WeakReactionType::BETA_PLUS_DECAY: return "bp"; + case WeakReactionType::BETA_MINUS_DECAY: return "bm"; + case WeakReactionType::ELECTRON_CAPTURE: return "ec"; + case WeakReactionType::POSITRON_CAPTURE: return "pc"; + } + return "Unknown"; + } + + struct WeakReactionEntry { + WeakReactionType type; + float T9; + float log_rhoYe; + float mu_e; + float log_rate; + float log_neutrino_loss; + + friend std::ostream& operator<<(std::ostream& os, const WeakReactionEntry& reaction) { + os << "WeakReactionEntry(type=" << weak_reaction_type_name(reaction.type) + << ", T9=" << reaction.T9 + << ", log_rhoYe=" << reaction.log_rhoYe + << ", mu_e=" << reaction.mu_e + << ", log_rate=" << reaction.log_rate + << ", log_neutrino_loss=" << reaction.log_neutrino_loss + << ")"; + return os; + } + }; +} \ No newline at end of file diff --git a/src/include/gridfire/screening/screening_abstract.h b/src/include/gridfire/screening/screening_abstract.h index f96038cf..1aafa874 100644 --- a/src/include/gridfire/screening/screening_abstract.h +++ b/src/include/gridfire/screening/screening_abstract.h @@ -44,7 +44,7 @@ namespace gridfire::screening { * * @param reactions The set of logical reactions in the network. * @param species A vector of all atomic species involved in the network. - * @param Y A vector of the molar abundances (mol/g) for each species. + * @param Y The current composition, providing molar abundances (mol/g) for each species. * @param T9 The temperature in units of 10^9 K. * @param rho The plasma density in g/cm^3. * @return A vector of screening factors (dimensionless), one for each reaction @@ -73,10 +73,10 @@ namespace gridfire::screening { [[nodiscard]] virtual std::vector calculateScreeningFactors( const reaction::ReactionSet& reactions, const std::vector& species, - const std::vector& Y, - const double T9, - const double rho - ) const = 0; + const std::vector &Y, + double T9, + double rho + ) const = 0; /** * @brief Calculates screening factors using CppAD types for automatic differentiation. @@ -88,7 +88,7 @@ namespace gridfire::screening { * * @param reactions The set of logical reactions in the network. * @param species A vector of all atomic species involved in the network. - * @param Y A vector of the molar abundances (mol/g) for each species, as AD types. + * @param Y The current composition, providing molar abundances (mol/g) for each species. * @param T9 The temperature in units of 10^9 K, as an AD type. * @param rho The plasma density in g/cm^3, as an AD type. * @return A vector of screening factors (dimensionless), as AD types. @@ -100,9 +100,9 @@ namespace gridfire::screening { [[nodiscard]] virtual std::vector calculateScreeningFactors( const reaction::ReactionSet& reactions, const std::vector& species, - const std::vector& Y, - const ADDouble T9, - const ADDouble rho + const std::vector> &Y, + ADDouble T9, + ADDouble rho ) const = 0; }; -} \ No newline at end of file +} diff --git a/src/include/gridfire/screening/screening_bare.h b/src/include/gridfire/screening/screening_bare.h index 5bc33d08..183be605 100644 --- a/src/include/gridfire/screening/screening_bare.h +++ b/src/include/gridfire/screening/screening_bare.h @@ -31,7 +31,7 @@ namespace gridfire::screening { * * @param reactions The set of logical reactions in the network. * @param species A vector of all atomic species (unused). - * @param Y A vector of the molar abundances (unused). + * @param Y A vector of the molar abundances. * @param T9 The temperature (unused). * @param rho The plasma density (unused). * @return A vector of doubles, with each element being 1.0, of the same @@ -54,9 +54,9 @@ namespace gridfire::screening { [[nodiscard]] std::vector calculateScreeningFactors( const reaction::ReactionSet& reactions, const std::vector& species, - const std::vector& Y, - const double T9, - const double rho + const std::vector &Y, + double T9, + double rho ) const override; /** @@ -68,7 +68,7 @@ namespace gridfire::screening { * * @param reactions The set of logical reactions in the network. * @param species A vector of all atomic species (unused). - * @param Y A vector of the molar abundances as AD types (unused). + * @param Y The current composition, providing molar abundances (mol/g) for each species (unused). * @param T9 The temperature as an AD type (unused). * @param rho The plasma density as an AD type (unused). * @return A vector of ADDouble, with each element being 1.0, of the same @@ -77,9 +77,9 @@ namespace gridfire::screening { [[nodiscard]] std::vector calculateScreeningFactors( const reaction::ReactionSet& reactions, const std::vector& species, - const std::vector& Y, - const ADDouble T9, - const ADDouble rho + const std::vector> &Y, + ADDouble T9, + ADDouble rho ) const override; private: /** @@ -92,7 +92,7 @@ namespace gridfire::screening { * @tparam T The numeric type, either `double` or `CppAD::AD`. * @param reactions The set of reactions for which to calculate factors. * @param species A vector of all atomic species (unused). - * @param Y A vector of molar abundances (unused). + * @param Y The current molar composition (unused). * @param T9 The temperature (unused). * @param rho The density (unused). * @return A vector of type `T` with all elements initialized to 1.0. @@ -102,8 +102,8 @@ namespace gridfire::screening { const reaction::ReactionSet& reactions, const std::vector& species, const std::vector& Y, - const T T9, - const T rho + const T& T9, + const T& rho ) const; }; @@ -117,7 +117,7 @@ namespace gridfire::screening { * @tparam T The numeric type, either `double` or `CppAD::AD`. * @param reactions The set of reactions, used to determine the size of the output vector. * @param species Unused parameter. - * @param Y Unused parameter. + * @param Y Unused parameter.` * @param T9 Unused parameter. * @param rho Unused parameter. * @return A `std::vector` of the same size as `reactions`, with all elements set to 1.0. @@ -127,8 +127,8 @@ namespace gridfire::screening { const reaction::ReactionSet &reactions, const std::vector &species, const std::vector &Y, - const T T9, - const T rho + const T& T9, + const T& rho ) const { return std::vector(reactions.size(), T(1.0)); // Bare screening returns 1.0 for all reactions } diff --git a/src/include/gridfire/screening/screening_intermediate.h b/src/include/gridfire/screening/screening_intermediate.h index 2f625456..c7f8da15 100644 --- a/src/include/gridfire/screening/screening_intermediate.h +++ b/src/include/gridfire/screening/screening_intermediate.h @@ -39,7 +39,7 @@ namespace gridfire::screening { std::vector IntermediateScreeningModel::calculateFactors_impl( const reaction::ReactionSet &reactions, const std::vector &species, - const std::vector &Y, + const std::vector& Y, const T T9, const T rho ) const { diff --git a/src/include/gridfire/screening/screening_weak.h b/src/include/gridfire/screening/screening_weak.h index ef8cd2ad..643a3ed8 100644 --- a/src/include/gridfire/screening/screening_weak.h +++ b/src/include/gridfire/screening/screening_weak.h @@ -33,7 +33,7 @@ namespace gridfire::screening { * * @param reactions The set of logical reactions in the network. * @param species A vector of all atomic species involved in the network. - * @param Y A vector of the molar abundances (mol/g) for each species. + * @param Y The composition object giving the current molar abundances. * @param T9 The temperature in units of 10^9 K. * @param rho The plasma density in g/cm^3. * @return A vector of screening factors (dimensionless), one for each reaction. @@ -50,9 +50,9 @@ namespace gridfire::screening { [[nodiscard]] std::vector calculateScreeningFactors( const reaction::ReactionSet& reactions, const std::vector& species, - const std::vector& Y, - const double T9, - const double rho + const std::vector &Y, + double T9, + double rho ) const override; /** @@ -64,7 +64,7 @@ namespace gridfire::screening { * * @param reactions The set of logical reactions in the network. * @param species A vector of all atomic species involved in the network. - * @param Y A vector of the molar abundances as AD types. + * @param Y The composition object giving the current molar abundances. * @param T9 The temperature as an AD type. * @param rho The plasma density as an AD type. * @return A vector of screening factors as AD types. @@ -72,9 +72,9 @@ namespace gridfire::screening { [[nodiscard]] std::vector> calculateScreeningFactors( const reaction::ReactionSet& reactions, const std::vector& species, - const std::vector>& Y, - const CppAD::AD T9, - const CppAD::AD rho + const std::vector> &Y, + CppAD::AD T9, + CppAD::AD rho ) const override; private: /// @brief Logger instance for recording trace and debug information. @@ -91,7 +91,7 @@ namespace gridfire::screening { * @tparam T The numeric type, either `double` or `CppAD::AD`. * @param reactions The set of reactions. * @param species A vector of all species in the network. - * @param Y A vector of molar abundances. + * @param Y The composition object with current molar abundances. * @param T9 The temperature in 10^9 K. * @param rho The density in g/cm^3. * @return A vector of screening factors of type `T`. @@ -101,8 +101,8 @@ namespace gridfire::screening { const reaction::ReactionSet& reactions, const std::vector& species, const std::vector& Y, - const T T9, - const T rho + T T9, + T rho ) const; }; @@ -115,7 +115,7 @@ namespace gridfire::screening { * @tparam T The numeric type (`double` or `CppAD::AD`). * @param reactions The set of reactions to be screened. * @param species The list of all species in the network. - * @param Y The molar abundances of the species. + * @param Y The composition object providing current molar abundances. * @param T9 The temperature in 10^9 K. * @param rho The density in g/cm^3. * @return A vector of screening factors, one for each reaction. diff --git a/src/include/gridfire/solver/solver.h b/src/include/gridfire/solver/solver.h index eb7ea88b..060a0448 100644 --- a/src/include/gridfire/solver/solver.h +++ b/src/include/gridfire/solver/solver.h @@ -1,14 +1,8 @@ #pragma once -#include "gridfire/engine/engine_graph.h" #include "gridfire/engine/engine_abstract.h" #include "gridfire/network.h" -#include "fourdst/logging/logging.h" -#include "fourdst/config/config.h" - -#include "quill/Logger.h" - #include #include #include @@ -101,249 +95,4 @@ namespace gridfire::solver { * @brief Type alias for a network solver strategy that uses a DynamicEngine. */ using DynamicNetworkSolverStrategy = NetworkSolverStrategy; - - /** - * @class DirectNetworkSolver - * @brief A network solver that directly integrates the reaction network ODEs. - * - * This solver uses a Runge-Kutta method to directly integrate the reaction network - * ODEs. It is simpler than the QSENetworkSolver, but it can be less efficient for - * stiff networks with disparate timescales. - * - * @implements DynamicNetworkSolverStrategy - */ - class DirectNetworkSolver final : public DynamicNetworkSolverStrategy { - public: - /** - * @brief Constructor for the DirectNetworkSolver. - * @param engine The dynamic engine to use for evaluating the network. - */ - using DynamicNetworkSolverStrategy::DynamicNetworkSolverStrategy; - - /** - * @struct TimestepContext - * @brief Context for the timestep callback function for the DirectNetworkSolver. - * - * This struct contains the context that will be passed to the callback function at the end of each timestep. - * It includes the current time, state, timestep size, cached results, and other relevant information. - * - * This type should be used when defining a callback function - * - * **Example:** - * @code - * #include "gridfire/solver/solver.h" - * - * #include - * #include - * - * static std::ofstream consumptionFile("consumption.txt"); - * void callback(const gridfire::solver::DirectNetworkSolver::TimestepContext& context) { - * int H1Index = context.engine.getSpeciesIndex(fourdst::atomic::H_1); - * int He4Index = context.engine.getSpeciesIndex(fourdst::atomic::He_4); - * - * consumptionFile << context.t << "," << context.state(H1Index) << "," << context.state(He4Index) << "\n"; - * } - * - * int main() { - * ... // Code to set up engine and solvers... - * solver.set_callback(callback); - * solver.evaluate(netIn); - * consumptionFile.close(); - * } - * @endcode - */ - struct TimestepContext final : public SolverContextBase { - const double t; ///< Current time. - const boost::numeric::ublas::vector& state; ///< Current state of the system. - const double dt; ///< Time step size. - const double cached_time; ///< Cached time for the last observed state. - const double last_observed_time; ///< Last time the state was observed. - const double last_step_time; ///< Last step time. - const double T9; ///< Temperature in units of 10^9 K. - const double rho; ///< Density in g/cm^3. - const std::optional>& cached_result; ///< Cached result of the step derivatives. - const int num_steps; ///< Total number of steps taken. - const DynamicEngine& engine; ///< Reference to the dynamic engine. - const std::vector& networkSpecies; - - TimestepContext( - const double t, - const boost::numeric::ublas::vector &state, - const double dt, - const double cached_time, - const double last_observed_time, - const double last_step_time, - const double t9, - const double rho, - const std::optional> &cached_result, - const int num_steps, - const DynamicEngine &engine, - const std::vector& networkSpecies - ); - - /** - * @brief Describe the context for callback functions. - * @return A vector of tuples, each containing a string for the parameter's name and a string for its type. - * - * This method provides a description of the context that will be passed to the callback function. - * The intent is that an end user can investigate the context and use this information to craft their own - * callback function. - * - * @implements SolverContextBase::describe - */ - [[nodiscard]] std::vector> describe() const override; - }; - - /** - * @brief Type alias for a timestep callback function. - * - * @brief The type alias for the callback function that will be called at the end of each timestep. - * - */ - using TimestepCallback = std::function; ///< Type alias for a timestep callback function. - - /** - * @brief Evaluates the network for a given timestep using direct integration. - * @param netIn The input conditions for the network. - * @return The output conditions after the timestep. - */ - NetOut evaluate(const NetIn& netIn) override; - - /** - * @brief Sets the callback function to be called at the end of each timestep. - * @param callback The callback function to be called at the end of each timestep. - * - * This function allows the user to set a callback function that will be called at the end of each timestep. - * The callback function will receive a gridfire::solver::DirectNetworkSolver::TimestepContext object. - */ - void set_callback(const std::any &callback) override; - - /** - * @brief Describe the context that will be passed to the callback function. - * @return A vector of tuples, each containing a string for the parameter's name and a string for its type. - * - * This method provides a description of the context that will be passed to the callback function. - * The intent is that an end user can investigate the context and use this information to craft their own - * callback function. - * - * @implements SolverContextBase::describe - */ - [[nodiscard]] std::vector> describe_callback_context() const override; - - - private: - /** - * @struct RHSManager - * @brief Functor for calculating the right-hand side of the ODEs. - * - * This functor is used by the ODE solver to calculate the time derivatives of the - * species abundances. It takes the current abundances as input and returns the - * time derivatives. - */ - struct RHSManager { - DynamicEngine& m_engine; ///< The engine used to evaluate the network. - const double m_T9; ///< Temperature in units of 10^9 K. - const double m_rho; ///< Density in g/cm^3. - - mutable double m_cached_time; - mutable std::optional> m_cached_result; - - mutable double m_last_observed_time = 0.0; ///< Last time the state was observed. - - - quill::Logger* m_logger = LogManager::getInstance().newFileLogger("integration.log", "GridFire"); ///< Logger instance. - mutable int m_num_steps = 0; - mutable double m_last_step_time = 1e-20; - - TimestepCallback& m_callback; - const std::vector& m_networkSpecies; - - /** - * @brief Constructor for the RHSFunctor. - * @param engine The engine used to evaluate the network. - * @param T9 Temperature in units of 10^9 K. - * @param rho Density in g/cm^3. - * @param callback callback function to be called at the end of each timestep. - * @param networkSpecies vector of species in the network in the correct order. - */ - RHSManager( - DynamicEngine& engine, - const double T9, - const double rho, - TimestepCallback& callback, - const std::vector& networkSpecies - ) : - m_engine(engine), - m_T9(T9), - m_rho(rho), - m_cached_time(0), - m_callback(callback), - m_networkSpecies(networkSpecies){} - - /** - * @brief Calculates the time derivatives of the species abundances. - * @param Y Vector of current abundances. - * @param dYdt Vector to store the time derivatives. - * @param t Current time. - */ - void operator()( - const boost::numeric::ublas::vector& Y, - boost::numeric::ublas::vector& dYdt, - double t - ) const; - - void observe(const boost::numeric::ublas::vector& state, double t) const; - void compute_and_cache(const boost::numeric::ublas::vector& state, double t) const; - - }; - - /** - * @struct JacobianFunctor - * @brief Functor for calculating the Jacobian matrix. - * - * This functor is used by the ODE solver to calculate the Jacobian matrix of the - * ODEs. It takes the current abundances as input and returns the Jacobian matrix. - */ - struct JacobianFunctor { - DynamicEngine& m_engine; ///< The engine used to evaluate the network. - const double m_T9; ///< Temperature in units of 10^9 K. - const double m_rho; ///< Density in g/cm^3. - - /** - * @brief Constructor for the JacobianFunctor. - * @param engine The engine used to evaluate the network. - * @param T9 Temperature in units of 10^9 K. - * @param rho Density in g/cm^3. - */ - JacobianFunctor( - DynamicEngine& engine, - const double T9, - const double rho - ) : - m_engine(engine), - m_T9(T9), - m_rho(rho) {} - - /** - * @brief Calculates the Jacobian matrix. - * @param Y Vector of current abundances. - * @param J Matrix to store the Jacobian matrix. - * @param t Current time. - * @param dfdt Vector to store the time derivatives (not used). - */ - void operator()( - const boost::numeric::ublas::vector& Y, - boost::numeric::ublas::matrix& J, - double t, - boost::numeric::ublas::vector& dfdt - ) const; - - }; - - private: - quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); ///< Logger instance. - Config& m_config = Config::getInstance(); ///< Configuration instance. - - TimestepCallback m_callback; - }; } \ No newline at end of file diff --git a/src/include/gridfire/solver/strategies/CVODE_solver_strategy.h b/src/include/gridfire/solver/strategies/CVODE_solver_strategy.h index ec73777c..7df1b7ca 100644 --- a/src/include/gridfire/solver/strategies/CVODE_solver_strategy.h +++ b/src/include/gridfire/solver/strategies/CVODE_solver_strategy.h @@ -55,7 +55,7 @@ namespace gridfire::solver { void set_callback(const std::any &callback) override; - bool get_stdout_logging_enabled() const; + [[nodiscard]] bool get_stdout_logging_enabled() const; void set_stdout_logging_enabled(const bool value); @@ -69,14 +69,14 @@ namespace gridfire::solver { const double last_step_time; const double T9; const double rho; - const int num_steps; + const size_t num_steps; const DynamicEngine& engine; const std::vector& networkSpecies; // Constructor TimestepContext( double t, const N_Vector& state, double dt, double last_step_time, - double t9, double rho, int num_steps, const DynamicEngine& engine, + double t9, double rho, size_t num_steps, const DynamicEngine& engine, const std::vector& networkSpecies ); @@ -104,8 +104,8 @@ namespace gridfire::solver { }; private: - Config& m_config = Config::getInstance(); - quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); + fourdst::config::Config& m_config = fourdst::config::Config::getInstance(); + quill::Logger* m_logger = fourdst::logging::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); diff --git a/src/include/gridfire/solver/strategies/triggers/engine_partitioning_trigger.h b/src/include/gridfire/solver/strategies/triggers/engine_partitioning_trigger.h index eb0d8fdb..c900605a 100644 --- a/src/include/gridfire/solver/strategies/triggers/engine_partitioning_trigger.h +++ b/src/include/gridfire/solver/strategies/triggers/engine_partitioning_trigger.h @@ -23,7 +23,7 @@ namespace gridfire::trigger::solver::CVODE { size_t numTriggers() const override; size_t numMisses() const override; private: - quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); + quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); mutable size_t m_hits = 0; mutable size_t m_misses = 0; mutable size_t m_updates = 0; @@ -48,7 +48,7 @@ namespace gridfire::trigger::solver::CVODE { size_t numTriggers() const override; size_t numMisses() const override; private: - quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); + quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); mutable size_t m_hits = 0; mutable size_t m_misses = 0; mutable size_t m_updates = 0; @@ -71,7 +71,7 @@ namespace gridfire::trigger::solver::CVODE { size_t numTriggers() const override; size_t numMisses() const override; private: - quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); + quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); mutable size_t m_hits = 0; mutable size_t m_misses = 0; mutable size_t m_updates = 0; diff --git a/src/include/gridfire/trigger/procedures/trigger_pprint.h b/src/include/gridfire/trigger/procedures/trigger_pprint.h index c1729462..f6721f89 100644 --- a/src/include/gridfire/trigger/procedures/trigger_pprint.h +++ b/src/include/gridfire/trigger/procedures/trigger_pprint.h @@ -5,7 +5,7 @@ #include namespace gridfire::trigger { - inline void printWhy(const TriggerResult& result, const int indent = 0) { + inline void printWhy(const TriggerResult& result, const int indent = 0) { // NOLINT(*-no-recursion) const std::string prefix(indent * 2, ' '); std::cout << prefix << "• [" << (result.value ? "TRUE" : "FALSE") << "] " << result.name << ": " << result.description << std::endl; diff --git a/src/include/gridfire/trigger/trigger_abstract.h b/src/include/gridfire/trigger/trigger_abstract.h index f2fe5ceb..20465de8 100644 --- a/src/include/gridfire/trigger/trigger_abstract.h +++ b/src/include/gridfire/trigger/trigger_abstract.h @@ -3,7 +3,6 @@ #include "gridfire/trigger/trigger_result.h" #include -#include namespace gridfire::trigger { template @@ -16,11 +15,11 @@ namespace gridfire::trigger { virtual void update(const TriggerContextStruct& ctx) = 0; virtual void reset() = 0; - virtual std::string name() const = 0; - virtual std::string describe() const = 0; - virtual TriggerResult why(const TriggerContextStruct& ctx) const = 0; + [[nodiscard]] virtual std::string name() const = 0; + [[nodiscard]] virtual std::string describe() const = 0; + [[nodiscard]] virtual TriggerResult why(const TriggerContextStruct& ctx) const = 0; - virtual size_t numTriggers() const = 0; - virtual size_t numMisses() const = 0; + [[nodiscard]] virtual size_t numTriggers() const = 0; + [[nodiscard]] virtual size_t numMisses() const = 0; }; } \ No newline at end of file diff --git a/src/include/gridfire/utils/logging.h b/src/include/gridfire/utils/logging.h index e08cf85a..e583dc71 100644 --- a/src/include/gridfire/utils/logging.h +++ b/src/include/gridfire/utils/logging.h @@ -4,11 +4,6 @@ #include #include -#include -#include -#include -#include -#include namespace gridfire::utils { @@ -64,9 +59,9 @@ namespace gridfire::utils { */ std::string formatNuclearTimescaleLogString( const DynamicEngine& engine, - const std::vector& Y, - const double T9, - const double rho + const fourdst::composition::Composition& composition, + double T9, + double rho ); diff --git a/src/include/gridfire/utils/table_format.h b/src/include/gridfire/utils/table_format.h index 0e8d0b00..7df92e9e 100644 --- a/src/include/gridfire/utils/table_format.h +++ b/src/include/gridfire/utils/table_format.h @@ -1,6 +1,7 @@ #pragma once #include +#include #include #include #include @@ -15,20 +16,20 @@ namespace gridfire::utils { public: virtual ~ColumnBase() = default; // Gets the string representation of the data at a given row - virtual std::string getCellData(size_t rowIndex) const = 0; + [[nodiscard]] virtual std::string getCellData(size_t rowIndex) const = 0; // Gets the header text for the column - virtual std::string getHeader() const = 0; + [[nodiscard]] virtual std::string getHeader() const = 0; // Gets the number of data rows in the column - virtual size_t getRowCount() const = 0; + [[nodiscard]] 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) {} + Column(std::string header, const std::vector& data) + : m_header(std::move(header)), m_data(data) {} - std::string getCellData(size_t rowIndex) const override { + [[nodiscard]] std::string getCellData(size_t rowIndex) const override { std::stringstream ss; if (rowIndex < m_data.size()) { ss << m_data[rowIndex]; @@ -36,11 +37,11 @@ namespace gridfire::utils { return ss.str(); } - std::string getHeader() const override { + [[nodiscard]] std::string getHeader() const override { return m_header; } - size_t getRowCount() const override { + [[nodiscard]] size_t getRowCount() const override { return m_data.size(); } private: @@ -74,8 +75,8 @@ namespace gridfire::utils { 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; + const size_t total_width = std::accumulate(col_widths.begin(), col_widths.end(), 0) + (num_cols * 3) + 1; // NOLINT(*-fold-init-type) + const 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 --- @@ -94,7 +95,7 @@ namespace gridfire::utils { // --- 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 << " " << std::left << std::setw(col_widths[j]) << columns[j]->getHeader() << " |"; // NOLINT(*-narrowing-conversions) } table_ss << "\n"; @@ -105,7 +106,7 @@ namespace gridfire::utils { 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 << " " << std::left << std::setw(col_widths[j]) << columns[j]->getCellData(i) << " |"; // NOLINT(*-narrowing-conversions) } table_ss << "\n"; } diff --git a/src/lib/engine/diagnostics/dynamic_engine_diagnostics.cpp b/src/lib/engine/diagnostics/dynamic_engine_diagnostics.cpp index bbe434c7..353f1d6d 100644 --- a/src/lib/engine/diagnostics/dynamic_engine_diagnostics.cpp +++ b/src/lib/engine/diagnostics/dynamic_engine_diagnostics.cpp @@ -5,7 +5,6 @@ #include #include #include -#include #include namespace gridfire::diagnostics { @@ -73,7 +72,7 @@ namespace gridfire::diagnostics { void inspect_species_balance( const DynamicEngine& engine, const std::string& species_name, - const std::vector& Y_full, + const fourdst::composition::Composition &comp, const double T9, const double rho ) { @@ -89,15 +88,15 @@ namespace gridfire::diagnostics { const int stoichiometry = reaction->stoichiometry(species_obj); if (stoichiometry == 0) continue; - const double flow = engine.calculateMolarReactionFlow(*reaction, Y_full, T9, rho); + const double flow = engine.calculateMolarReactionFlow(*reaction, comp, T9, rho); if (stoichiometry > 0) { - creation_ids.push_back(std::string(reaction->id())); + creation_ids.emplace_back(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_ids.emplace_back(reaction->id()); destruction_stoichiometry.push_back(stoichiometry); destruction_flows.push_back(flow); total_destruction_flow += std::abs(stoichiometry) * flow; @@ -129,38 +128,38 @@ namespace gridfire::diagnostics { void inspect_jacobian_stiffness( const DynamicEngine& engine, - const std::vector& Y_full, + const fourdst::composition::Composition &comp, const double T9, const double rho ) { - engine.generateJacobianMatrix(Y_full, T9, rho); + engine.generateJacobianMatrix(comp, 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; + std::optional max_diag_species = std::nullopt; + std::optional> max_off_diag_species = std::nullopt; - 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; } + for (const auto& rowSpecies : species_list) { + for (const auto& colSpecies : species_list) { + const double val = std::abs(engine.getJacobianMatrixEntry(rowSpecies, colSpecies)); + if (rowSpecies == colSpecies) { + if (val > max_diag) { max_diag = val; max_diag_species = colSpecies; } } else { - if (val > max_off_diag) { max_off_diag = val; max_off_diag_i = i; max_off_diag_j = j; } + if (val > max_off_diag) { max_off_diag = val; max_off_diag_species = {rowSpecies, colSpecies};} } } } std::cout << "\n--- Jacobian Stiffness Report ---" << std::endl; - if (max_diag_idx != -1) { + if (max_diag_species.has_value()) { std::cout << " Largest Diagonal Element (d(dYi/dt)/dYi): " << std::scientific << max_diag - << " for species " << species_list[max_diag_idx].name() << std::endl; + << " for species " << max_diag_species->name() << std::endl; } - if (max_off_diag_i != -1) { + if (max_off_diag_species.has_value()) { 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; + << " for d(" << max_off_diag_species->first.name() + << ")/d(" << max_off_diag_species->second.name() << ")" << std::endl; } std::cout << "---------------------------------" << std::endl; } diff --git a/src/lib/engine/engine_graph.cpp b/src/lib/engine/engine_graph.cpp index 352586db..c84f0c1a 100644 --- a/src/lib/engine/engine_graph.cpp +++ b/src/lib/engine/engine_graph.cpp @@ -12,7 +12,6 @@ #include "quill/LogMacros.h" #include -#include #include #include #include @@ -40,7 +39,8 @@ namespace gridfire { const fourdst::composition::Composition &composition, const partition::PartitionFunction& partitionFunction, const BuildDepthType buildDepth) : - m_reactions(build_reaclib_nuclear_network(composition, buildDepth, false)), + m_weakRateInterpolator(rates::weak::UNIFIED_WEAK_DATA), + m_reactions(build_reaclib_nuclear_network(composition, m_weakRateInterpolator, buildDepth, false)), m_depth(buildDepth), m_partitionFunction(partitionFunction.clone()) { @@ -50,15 +50,19 @@ namespace gridfire { GraphEngine::GraphEngine( const reaction::ReactionSet &reactions ) : - m_reactions(reactions) { + m_weakRateInterpolator(rates::weak::UNIFIED_WEAK_DATA), + m_reactions(reactions) + { syncInternalMaps(); } std::expected, expectations::StaleEngineError> GraphEngine::calculateRHSAndEnergy( - const std::vector &Y, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { + const double Ye = comp.getElectronAbundance(); + const double mue = 5.0e-3 * std::pow(rho * Ye, 1.0 / 3.0) + 0.5 * T9; if (m_usePrecomputation) { std::vector bare_rates; std::vector bare_reverse_rates; @@ -66,33 +70,35 @@ namespace gridfire { bare_reverse_rates.reserve(m_reactions.size()); // TODO: Add cache to this + for (const auto& reaction: m_reactions) { - bare_rates.push_back(reaction->calculate_rate(T9, rho, Y)); - bare_reverse_rates.push_back(calculateReverseRate(*reaction, T9, rho, Y)); + bare_rates.push_back(reaction->calculate_rate(T9, rho, Ye, mue, comp.getMolarAbundanceVector(), m_indexToSpeciesMap)); + bare_reverse_rates.push_back(calculateReverseRate(*reaction, T9, rho, comp)); } // --- The public facing interface can always use the precomputed version since taping is done internally --- - return calculateAllDerivativesUsingPrecomputation(Y, bare_rates, bare_reverse_rates, T9, rho); + return calculateAllDerivativesUsingPrecomputation(comp, bare_rates, bare_reverse_rates, T9, rho); } else { - return calculateAllDerivatives(Y, T9, rho); + return calculateAllDerivatives(comp.getMolarAbundanceVector(), T9, rho, Ye, mue); } } EnergyDerivatives GraphEngine::calculateEpsDerivatives( - const std::vector &Y, + const fourdst::composition::Composition &comp, 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) { + if (comp.getRegisteredSpecies().size() != numSpecies) { LOG_ERROR(m_logger, "Input abundance vector size ({}) does not match number of species in the network ({}).", - Y.size(), numSpecies); + comp.getRegisteredSpecies().size(), numSpecies); throw std::invalid_argument("Input abundance vector size does not match number of species in the network."); } std::vector x(numADInputs); + const std::vector Y = comp.getMolarAbundanceVector(); for (size_t i = 0; i < numSpecies; ++i) { x[i] = Y[i]; } @@ -148,6 +154,7 @@ namespace gridfire { // --- Network Graph Construction Methods --- void GraphEngine::collectNetworkSpecies() { + //TODO: Ensure consistent ordering in the m_networkSpecies vector so that it is ordered by species mass. m_networkSpecies.clear(); m_networkSpeciesMap.clear(); @@ -190,6 +197,10 @@ namespace gridfire { for (size_t i = 0; i < m_networkSpecies.size(); ++i) { m_speciesToIndexMap.insert({m_networkSpecies[i], i}); } + m_indexToSpeciesMap.clear(); + for (size_t i = 0; i < m_networkSpecies.size(); ++i) { + m_indexToSpeciesMap.insert({i, m_networkSpecies[i]}); + } } void GraphEngine::reserveJacobianMatrix() const { @@ -287,13 +298,18 @@ namespace gridfire { const reaction::Reaction &reaction, const double T9, const double rho, - const std::vector &Y + const fourdst::composition::Composition &comp ) const { if (!m_useReverseReactions) { LOG_TRACE_L3_LIMIT_EVERY_N(std::numeric_limits::max(), m_logger, "Reverse reactions are disabled. Returning 0.0 for reverse rate of reaction '{}'.", reaction.id()); return 0.0; // If reverse reactions are not used, return 0.0 } const double temp = T9 * 1e9; // Convert T9 to Kelvin + const double Ye = comp.getElectronAbundance(); + + // TODO: This is a dummy value for the electron chemical potential. We eventually need to replace this with an EOS call. + const double mue = 5.0e-3 * std::pow(rho * Ye, 1.0 / 3.0) + 0.5 * T9; + // In debug builds we check the units on kB to ensure it is in erg/K. This is removed in release builds to avoid overhead. (Note assert is a no-op in release builds) assert(Constants::getInstance().get("kB").unit == "erg / K"); @@ -301,7 +317,7 @@ namespace gridfire { const double kBMeV = m_constants.kB * 624151; // Convert kB to MeV/K NOTE: This relies on the fact that m_constants.kB is in erg/K! const double expFactor = std::exp(-reaction.qValue() / (kBMeV * temp)); double reverseRate = 0.0; - const double forwardRate = reaction.calculate_rate(T9, rho, Y); + const double forwardRate = reaction.calculate_rate(T9, rho, Ye, mue, comp.getMolarAbundanceVector(), m_indexToSpeciesMap); if (reaction.reactants().size() == 2 && reaction.products().size() == 2) { reverseRate = calculateReverseRateTwoBody(reaction, T9, forwardRate, expFactor); @@ -393,14 +409,17 @@ namespace gridfire { const reaction::Reaction &reaction, const double T9, const double rho, - const std::vector &Y, + const fourdst::composition::Composition& comp, const double reverseRate ) const { if (!m_useReverseReactions) { LOG_TRACE_L3_LIMIT_EVERY_N(std::numeric_limits::max(), m_logger, "Reverse reactions are disabled. Returning 0.0 for reverse rate of reaction '{}'.", reaction.id()); return 0.0; // If reverse reactions are not used, return 0.0 } - const double d_log_kFwd = reaction.calculate_forward_rate_log_derivative(T9, rho, Y); + double Ye = comp.getElectronAbundance(); + // TODO: This is a dummy value for the electron chemical potential. We eventually need to replace this with an EOS call. + double mue = 5.0e-3 * std::pow(rho * Ye, 1.0 / 3.0) + 0.5 * T9; + const double d_log_kFwd = reaction.calculate_forward_rate_log_derivative(T9, rho, Ye, mue, comp); auto log_deriv_pf_op = [&](double acc, const auto& species) { const double g = m_partitionFunction->evaluate(species.z(), species.a(), T9); @@ -486,7 +505,7 @@ namespace gridfire { void GraphEngine::rebuild(const fourdst::composition::Composition& comp, const BuildDepthType depth) { if (depth != m_depth) { m_depth = depth; - m_reactions = build_reaclib_nuclear_network(comp, m_depth, false); + m_reactions = build_reaclib_nuclear_network(comp, m_weakRateInterpolator, m_depth, false); 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."); @@ -494,7 +513,7 @@ namespace gridfire { } StepDerivatives GraphEngine::calculateAllDerivativesUsingPrecomputation( - const std::vector &Y_in, + const fourdst::composition::Composition& comp, const std::vector &bare_rates, const std::vector &bare_reverse_rates, const double T9, @@ -504,33 +523,27 @@ namespace gridfire { const std::vector screeningFactors = m_screeningModel->calculateScreeningFactors( m_reactions, m_networkSpecies, - Y_in, + comp.getMolarAbundanceVector(), T9, rho ); + // TODO: Fix up the precomputation to use the new comp in interface as opposed to a raw vector of molar abundances + // This will require carefully checking the way the precomputation is stashed. + // --- Optimized loop --- std::vector molarReactionFlows; molarReactionFlows.reserve(m_precomputedReactions.size()); for (const auto& precomp : m_precomputedReactions) { double forwardAbundanceProduct = 1.0; - // bool below_threshold = false; for (size_t i = 0; i < precomp.unique_reactant_indices.size(); ++i) { const size_t reactantIndex = precomp.unique_reactant_indices[i]; + const fourdst::atomic::Species& reactant = m_networkSpecies[reactantIndex]; const int power = precomp.reactant_powers[i]; - // const double abundance = Y_in[reactantIndex]; - // if (abundance < MIN_ABUNDANCE_THRESHOLD) { - // below_threshold = true; - // break; - // } - forwardAbundanceProduct *= std::pow(Y_in[reactantIndex], power); + forwardAbundanceProduct *= std::pow(comp.getMolarAbundance(reactant), power); } - // if (below_threshold) { - // molarReactionFlows.push_back(0.0); - // continue; // Skip this reaction if any reactant is below the abundance threshold - // } const double bare_rate = bare_rates[precomp.reaction_index]; const double screeningFactor = screeningFactors[precomp.reaction_index]; @@ -549,7 +562,9 @@ namespace gridfire { const double bare_reverse_rate = bare_reverse_rates[precomp.reaction_index]; double reverseAbundanceProduct = 1.0; for (size_t i = 0; i < precomp.unique_product_indices.size(); ++i) { - reverseAbundanceProduct *= std::pow(Y_in[precomp.unique_product_indices[i]], precomp.product_powers[i]); + const size_t productIndex = precomp.unique_product_indices[i]; + const fourdst::atomic::Species& product = m_networkSpecies[productIndex]; + reverseAbundanceProduct *= std::pow(comp.getMolarAbundance(product), precomp.product_powers[i]); } reverseMolarReactionFlow = screeningFactor * bare_reverse_rate * @@ -564,25 +579,28 @@ namespace gridfire { // --- Assemble molar abundance derivatives --- StepDerivatives result; - result.dydt.assign(m_networkSpecies.size(), 0.0); // Initialize derivatives to zero + for (const auto& species: m_networkSpecies) { + result.dydt[species] = 0.0; // Initialize the change in abundance for each network species to 0 + } for (size_t j = 0; j < m_precomputedReactions.size(); ++j) { const auto& precomp = m_precomputedReactions[j]; const double R_j = molarReactionFlows[j]; for (size_t i = 0; i < precomp.affected_species_indices.size(); ++i) { const size_t speciesIndex = precomp.affected_species_indices[i]; + const fourdst::atomic::Species& species = m_networkSpecies[speciesIndex]; + const int stoichiometricCoefficient = precomp.stoichiometric_coefficients[i]; // Update the derivative for this species - result.dydt[speciesIndex] += static_cast(stoichiometricCoefficient) * R_j; + result.dydt.at(species) += static_cast(stoichiometricCoefficient) * R_j; } } // --- Calculate the nuclear energy generation rate --- double massProductionRate = 0.0; // [mol][s^-1] - for (size_t i = 0; i < m_networkSpecies.size(); ++i) { - const auto& species = m_networkSpecies[i]; - massProductionRate += result.dydt[i] * species.mass() * m_constants.u; + for (const auto & species : m_networkSpecies) { + massProductionRate += result.dydt[species] * species.mass() * m_constants.u; } result.nuclearEnergyGenerationRate = -massProductionRate * m_constants.Na * m_constants.c * m_constants.c; // [erg][s^-1][g^-1] return result; @@ -631,21 +649,22 @@ namespace gridfire { m_stoichiometryMatrix.nnz()); // Assuming nnz() exists for compressed_matrix } - StepDerivatives GraphEngine::calculateAllDerivatives( - const std::vector &Y_in, - const double T9, - const double rho - ) const { - return calculateAllDerivatives(Y_in, T9, rho); - } - - StepDerivatives GraphEngine::calculateAllDerivatives( - const std::vector &Y_in, - const ADDouble &T9, - const ADDouble &rho - ) const { - return calculateAllDerivatives(Y_in, T9, rho); - } + // StepDerivatives GraphEngine::calculateAllDerivatives( + // const std::vector &Y, + // const double T9, + // const double rho + // ) const { + // return calculateAllDerivatives(Y, T9, rho); + // } + // + // StepDerivatives GraphEngine::calculateAllDerivatives( + // const std::vector &Y, + // ADDouble T9, + // ADDouble rho + // ) const { + // //TODO: Sort out why this template is not being found + // return calculateAllDerivatives(Y, T9, rho); + // } void GraphEngine::setScreeningModel(const screening::ScreeningType model) { m_screeningModel = screening::selectScreeningModel(model); @@ -670,15 +689,21 @@ namespace gridfire { double GraphEngine::calculateMolarReactionFlow( const reaction::Reaction &reaction, - const std::vector &Y, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { - return calculateMolarReactionFlow(reaction, Y, T9, rho); + + const double Ye = comp.getElectronAbundance(); + + // TODO: This is a dummy placeholder which must be replaced with an EOS call + const double mue = 5.0e-3 * std::pow(rho * Ye, 1.0 / 3.0) + 0.5 * T9; + + return calculateMolarReactionFlow(reaction, comp.getMolarAbundanceVector(), T9, rho, Ye, mue); } void GraphEngine::generateJacobianMatrix( - const std::vector &Y_dynamic, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { @@ -688,6 +713,7 @@ namespace gridfire { // 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(); for (size_t i = 0; i < numSpecies; ++i) { adInput[i] = std::max(Y_dynamic[i], 1e-99); // regularize the jacobian... } @@ -711,7 +737,7 @@ namespace gridfire { } void GraphEngine::generateJacobianMatrix( - const std::vector &Y_dynamic, + const fourdst::composition::Composition &comp, const double T9, const double rho, const SparsityPattern &sparsityPattern @@ -719,6 +745,7 @@ namespace gridfire { // --- Pack the input variables into a vector for CppAD --- const size_t numSpecies = m_networkSpecies.size(); std::vector x(numSpecies + 2, 0.0); + const std::vector& Y_dynamic = comp.getMolarAbundanceVector(); for (size_t i = 0; i < numSpecies; ++i) { x[i] = Y_dynamic[i]; } @@ -767,8 +794,13 @@ namespace gridfire { } } - double GraphEngine::getJacobianMatrixEntry(const int i, const int j) const { - // LOG_TRACE_L3(m_logger, "Getting jacobian matrix entry for {},{} = {}", i, j, m_jacobianMatrix(i, j)); + 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); } @@ -779,10 +811,10 @@ namespace gridfire { } int GraphEngine::getStoichiometryMatrixEntry( - const int speciesIndex, - const int reactionIndex + const fourdst::atomic::Species& species, + const reaction::Reaction &reaction ) const { - return m_stoichiometryMatrix(speciesIndex, reactionIndex); + return reaction.stoichiometry(species); } void GraphEngine::exportToDot(const std::string &filename) const { @@ -871,18 +903,21 @@ namespace gridfire { } std::expected, expectations::StaleEngineError> GraphEngine::getSpeciesTimescales( - const std::vector &Y, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { - auto [dydt, _] = calculateAllDerivatives(Y, T9, rho); + const double Ye = comp.getElectronAbundance(); + // TODO: This is a dummy placeholder which must be replaced with an EOS call + const double mue = 5.0e-3 * std::pow(rho * Ye, 1.0 / 3.0) + 0.5 * T9; + + auto [dydt, _] = calculateAllDerivatives(comp.getMolarAbundanceVector(), T9, rho, Ye, mue); std::unordered_map speciesTimescales; speciesTimescales.reserve(m_networkSpecies.size()); - for (size_t i = 0; i < m_networkSpecies.size(); ++i) { + for (const auto& species : m_networkSpecies) { double timescale = std::numeric_limits::infinity(); - const auto species = m_networkSpecies[i]; - if (std::abs(dydt[i]) > 0.0) { - timescale = std::abs(Y[i] / dydt[i]); + if (std::abs(dydt.at(species)) > 0.0) { + timescale = std::abs(comp.getMolarAbundance(species) / dydt.at(species)); } speciesTimescales.emplace(species, timescale); } @@ -890,24 +925,28 @@ namespace gridfire { } std::expected, expectations::StaleEngineError> GraphEngine::getSpeciesDestructionTimescales( - const std::vector &Y, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { - auto [dydt, _] = calculateAllDerivatives(Y, T9, rho); + const double Ye = comp.getElectronAbundance(); + // TODO: This is a dummy placeholder which must be replaced with an EOS call + const double mue = 5.0e-3 * std::pow(rho * Ye, 1.0 / 3.0) + 0.5 * T9; + + auto [dydt, _] = calculateAllDerivatives(comp.getMolarAbundanceVector(), T9, rho, Ye, mue); std::unordered_map speciesDestructionTimescales; speciesDestructionTimescales.reserve(m_networkSpecies.size()); for (const auto& species : m_networkSpecies) { double netDestructionFlow = 0.0; for (const auto& reaction : m_reactions) { if (reaction->stoichiometry(species) < 0) { - const auto flow = calculateMolarReactionFlow(*reaction, Y, T9, rho); + const auto flow = calculateMolarReactionFlow(*reaction, comp.getMolarAbundanceVector(), T9, rho, Ye, mue); netDestructionFlow += flow; } } double timescale = std::numeric_limits::infinity(); if (netDestructionFlow != 0.0) { - timescale = Y[getSpeciesIndex(species)] / netDestructionFlow; + timescale = comp.getMolarAbundance(species) / netDestructionFlow; } speciesDestructionTimescales.emplace(species, timescale); } @@ -964,12 +1003,21 @@ namespace gridfire { const CppAD::AD adT9 = adInput[numSpecies]; const CppAD::AD adRho = adInput[numSpecies + 1]; + // Dummy values for Ye and mue to let taping happen + const CppAD::AD adYe = 1.0; + const CppAD::AD adMue = 1.0; + // 5. Call the actual templated function // We let T9 and rho be constant, so we pass them as fixed values. - auto [dydt, nuclearEnergyGenerationRate] = calculateAllDerivatives>(adY, adT9, adRho); + auto [dydt, nuclearEnergyGenerationRate] = calculateAllDerivatives>(adY, adT9, adRho, adYe, adMue); - m_rhsADFun.Dependent(adInput, dydt); + // Extract the raw vector from the associative map + std::vector> dydt_vec; + dydt_vec.reserve(dydt.size()); + std::ranges::transform(dydt, std::back_inserter(dydt_vec),[](const auto& kv) { return kv.second; }); + + m_rhsADFun.Dependent(adInput, dydt_vec); LOG_TRACE_L1(m_logger, "AD tape recorded successfully for the RHS calculation. Number of independent variables: {}.", adInput.size()); @@ -996,16 +1044,19 @@ namespace gridfire { CppAD::Independent(adInput); - const std::vector> adY(adInput.begin(), adInput.begin() + numSpecies); + const std::vector> adY(adInput.begin(), adInput.begin() + static_cast(numSpecies)); const CppAD::AD adT9 = adInput[numSpecies]; const CppAD::AD adRho = adInput[numSpecies + 1]; - StepDerivatives> derivatives = calculateAllDerivatives>(adY, adT9, adRho); + // Dummy values for Ye and mue to let taping happen + const CppAD::AD adYe = 1.0; + const CppAD::AD adMue = 1.0; + + auto [dydt, nuclearEnergyGenerationRate] = calculateAllDerivatives>(adY, adT9, adRho, adYe, adMue); std::vector> adOutput(1); - adOutput[0] = derivatives.nuclearEnergyGenerationRate; + adOutput[0] = 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()); diff --git a/src/lib/engine/procedures/construction.cpp b/src/lib/engine/procedures/construction.cpp index 55ee8755..02535ba3 100644 --- a/src/lib/engine/procedures/construction.cpp +++ b/src/lib/engine/procedures/construction.cpp @@ -1,7 +1,11 @@ #include "gridfire/engine/procedures/construction.h" +#include "gridfire/reaction/weak/weak_interpolator.h" +#include "gridfire/reaction/weak/weak.h" +#include "gridfire/reaction/weak/weak_types.h" #include #include +#include #include "gridfire/reaction/reaction.h" #include "gridfire/reaction/reaclib.h" @@ -20,11 +24,11 @@ namespace gridfire { using fourdst::atomic::Species; - ReactionSet build_reaclib_nuclear_network( - const Composition &composition, - BuildDepthType maxLayers, - bool reverse - ) { + ReactionSet build_nuclear_network( + const Composition& composition, + const rates::weak::WeakRateInterpolator& weak_interpolator, + BuildDepthType maxLayers, + bool reverse_reaclib) { int depth; if (std::holds_alternative(maxLayers)) { depth = static_cast(std::get(maxLayers)); @@ -37,31 +41,71 @@ namespace gridfire { throw std::logic_error("Network build depth is set to 0. No reactions will be collected."); } - const auto& allReactions = reaclib::get_all_reaclib_reactions(); - std::vector remainingReactions; - for (const auto& reaction : allReactions) { - if (reaction->is_reverse() == reverse) { - remainingReactions.push_back(reaction.get()); + // --- Step 1: Create a single master pool that owns all possible reactions --- + ReactionSet master_reaction_pool; + + // Clone all relevant REACLIB reactions into the master pool + const auto& allReaclibReactions = reaclib::get_all_reaclib_reactions(); + for (const auto& reaction : allReaclibReactions) { + if (reaction->is_reverse() == reverse_reaclib) { + master_reaction_pool.add_reaction(reaction->clone()); } } - if (depth == static_cast(NetworkBuildDepth::Full)) { - LOG_INFO(logger, "Building full nuclear network with a total of {} reactions.", allReactions.size()); - const ReactionSet reactionSet(remainingReactions); - return reactionSet; + for (const auto& parent_species: weak_interpolator.available_isotopes()) { + master_reaction_pool.add_reaction( + std::make_unique( + parent_species, + rates::weak::WeakReactionType::BETA_PLUS_DECAY, + weak_interpolator + ) + ); + master_reaction_pool.add_reaction( + std::make_unique( + parent_species, + rates::weak::WeakReactionType::BETA_MINUS_DECAY, + weak_interpolator + ) + ); + master_reaction_pool.add_reaction( + std::make_unique( + parent_species, + rates::weak::WeakReactionType::ELECTRON_CAPTURE, + weak_interpolator + ) + ); + master_reaction_pool.add_reaction( + std::make_unique( + parent_species, + rates::weak::WeakReactionType::POSITRON_CAPTURE, + weak_interpolator + ) + ); } + // --- Step 2: Use non-owning raw pointers for the fast build algorithm --- + std::vector remainingReactions; + remainingReactions.reserve(master_reaction_pool.size()); + for(const auto& reaction : master_reaction_pool) { + remainingReactions.push_back(reaction.get()); + } + + if (depth == static_cast(NetworkBuildDepth::Full)) { + LOG_INFO(logger, "Building full nuclear network with a total of {} reactions.", remainingReactions.size()); + return master_reaction_pool; // No need to build layer by layer if we want the full network + } + + // --- Step 3: Execute the layered network build using observing pointers --- std::unordered_set availableSpecies; - for (const auto &entry: composition | std::views::values) { + for (const auto& entry : composition | std::views::values) { if (entry.mass_fraction() > 0.0) { availableSpecies.insert(entry.isotope()); } } - - std::vector collectedReactions; - + std::vector collectedReactionPtrs; LOG_INFO(logger, "Starting network construction with {} available species.", availableSpecies.size()); + for (int layer = 0; layer < depth && !remainingReactions.empty(); ++layer) { LOG_TRACE_L1(logger, "Collecting reactions for layer {} with {} remaining reactions. Currently there are {} available species", layer, remainingReactions.size(), availableSpecies.size()); std::vector reactionsForNextPass; @@ -70,7 +114,7 @@ namespace gridfire { reactionsForNextPass.reserve(remainingReactions.size()); - for (const auto &reaction : remainingReactions) { + for (Reaction* reaction : remainingReactions) { bool allReactantsAvailable = true; for (const auto& reactant : reaction->reactants()) { if (!availableSpecies.contains(reactant)) { @@ -80,7 +124,7 @@ namespace gridfire { } if (allReactantsAvailable) { - collectedReactions.push_back(reaction); + collectedReactionPtrs.push_back(reaction); newReactionsAdded = true; for (const auto& product : reaction->products()) { @@ -92,19 +136,29 @@ namespace gridfire { } if (!newReactionsAdded) { - LOG_INFO(logger, "No new reactions added in layer {}. Stopping network construction with {} reactions collected.", layer, collectedReactions.size()); + LOG_INFO(logger, "No new reactions added in layer {}. Stopping network construction.", layer); break; } - LOG_TRACE_L1(logger, "Layer {}: Collected {} reactions. New products this layer: {}", layer, collectedReactions.size(), newProductsThisLayer.size()); + LOG_TRACE_L1(logger, "Layer {}: Collected {} new reactions. New products this layer: {}", collectedReactionPtrs.size() - collectedReactionPtrs.size(), newProductsThisLayer.size()); availableSpecies.insert(newProductsThisLayer.begin(), newProductsThisLayer.end()); - remainingReactions = std::move(reactionsForNextPass); } - LOG_INFO(logger, "Network construction completed with {} reactions collected.", collectedReactions.size()); - const ReactionSet reactionSet(collectedReactions); - return reactionSet; + // --- Step 4: Construct the final ReactionSet by moving ownership --- + LOG_INFO(logger, "Network construction completed. Assembling final set of {} reactions.", collectedReactionPtrs.size()); + ReactionSet finalReactionSet; + + std::unordered_set collectedPtrSet(collectedReactionPtrs.begin(), collectedReactionPtrs.end()); + + for (auto& reaction_ptr : master_reaction_pool) { + if (reaction_ptr && collectedPtrSet.contains(reaction_ptr.get())) { + finalReactionSet.add_reaction(std::move(reaction_ptr)); + } + } + + return finalReactionSet; } + } \ No newline at end of file diff --git a/src/lib/engine/procedures/priming.cpp b/src/lib/engine/procedures/priming.cpp index 8e8600fc..52673b0b 100644 --- a/src/lib/engine/procedures/priming.cpp +++ b/src/lib/engine/procedures/priming.cpp @@ -17,7 +17,7 @@ namespace gridfire { const reaction::Reaction* findDominantCreationChannel ( const DynamicEngine& engine, const Species& species, - const std::vector& Y, + const fourdst::composition::Composition &comp, const double T9, const double rho ) { @@ -25,7 +25,7 @@ namespace gridfire { double maxFlow = -1.0; for (const auto& reaction : engine.getNetworkReactions()) { if (reaction->contains(species) && reaction->stoichiometry(species) > 0) { - const double flow = engine.calculateMolarReactionFlow(*reaction, Y, T9, rho); + const double flow = engine.calculateMolarReactionFlow(*reaction, comp, T9, rho); if (flow > maxFlow) { maxFlow = flow; dominateReaction = reaction.get(); @@ -67,7 +67,7 @@ namespace gridfire { * robustly primed composition. */ PrimingReport primeNetwork(const NetIn& netIn, DynamicEngine& engine) { - auto logger = LogManager::getInstance().getLogger("log"); + auto logger = fourdst::logging::LogManager::getInstance().getLogger("log"); // --- Initial Setup --- // Identify all species with zero initial mass fraction that need to be primed. @@ -127,9 +127,6 @@ namespace gridfire { } 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()); @@ -138,15 +135,14 @@ namespace gridfire { continue; } - const auto Y = primer.mapNetInToMolarAbundanceVector(tempNetIn); - const double destructionRateConstant = calculateDestructionRateConstant(primer, primingSpecies, Y, T9, rho); + const double destructionRateConstant = calculateDestructionRateConstant(primer, primingSpecies, tempComp, T9, rho); if (destructionRateConstant > 1e-99) { - const double creationRate = calculateCreationRate(primer, primingSpecies, Y, T9, rho); + const double creationRate = calculateCreationRate(primer, primingSpecies, tempComp, T9, rho); double equilibriumMassFraction = (creationRate / destructionRateConstant) * primingSpecies.mass(); if (std::isnan(equilibriumMassFraction)) equilibriumMassFraction = 0.0; - if (const reaction::Reaction* dominantChannel = findDominantCreationChannel(primer, primingSpecies, Y, T9, rho)) { + if (const reaction::Reaction* dominantChannel = findDominantCreationChannel(primer, primingSpecies, tempComp, T9, rho)) { // Store the request instead of applying it immediately. requests.push_back({primingSpecies, equilibriumMassFraction, dominantChannel->reactants()}); } else { @@ -403,19 +399,17 @@ namespace gridfire { double calculateDestructionRateConstant( const DynamicEngine& engine, - const fourdst::atomic::Species& species, - const std::vector& Y, + const Species& species, + const Composition& comp, const double T9, const double rho ) { - const size_t speciesIndex = engine.getSpeciesIndex(species); - std::vector Y_scaled(Y.begin(), Y.end()); - Y_scaled[speciesIndex] = 1.0; // Set the abundance of the species to 1.0 for rate constant calculation + //TODO: previously (when using raw vectors) I set y[speciesIndex] = 1.0 to let there be enough so that just the destruction rate could be found (without bottlenecks from abundance) we will need to do a similar thing here. double destructionRateConstant = 0.0; for (const auto& reaction: engine.getNetworkReactions()) { if (reaction->contains_reactant(species)) { const int stoichiometry = reaction->stoichiometry(species); - destructionRateConstant += std::abs(stoichiometry) * engine.calculateMolarReactionFlow(*reaction, Y_scaled, T9, rho); + destructionRateConstant += std::abs(stoichiometry) * engine.calculateMolarReactionFlow(*reaction, comp, T9, rho); } } return destructionRateConstant; @@ -424,7 +418,7 @@ namespace gridfire { double calculateCreationRate( const DynamicEngine& engine, const Species& species, - const std::vector& Y, + const Composition& comp, const double T9, const double rho ) { @@ -432,9 +426,9 @@ namespace gridfire { for (const auto& reaction: engine.getNetworkReactions()) { const int stoichiometry = reaction->stoichiometry(species); if (stoichiometry > 0) { - if (engine.calculateMolarReactionFlow(*reaction, Y, T9, rho) > 0.0) { + if (engine.calculateMolarReactionFlow(*reaction, comp, T9, rho) > 0.0) { } - creationRate += stoichiometry * engine.calculateMolarReactionFlow(*reaction, Y, T9, rho); + creationRate += stoichiometry * engine.calculateMolarReactionFlow(*reaction, comp, T9, rho); } } return creationRate; diff --git a/src/lib/engine/views/engine_adaptive.cpp b/src/lib/engine/views/engine_adaptive.cpp index 7b906088..684ce22a 100644 --- a/src/lib/engine/views/engine_adaptive.cpp +++ b/src/lib/engine/views/engine_adaptive.cpp @@ -95,8 +95,7 @@ namespace gridfire { LOG_TRACE_L1(m_logger, "Updating AdaptiveEngineView with new network input..."); - std::vector Y_Full; - std::vector allFlows = calculateAllReactionFlows(updatedNetIn, Y_Full); + auto [allFlows, composition] = calculateAllReactionFlows(updatedNetIn); double maxFlow = 0.0; @@ -110,11 +109,11 @@ namespace gridfire { const std::unordered_set reachableSpecies = findReachableSpecies(updatedNetIn); LOG_DEBUG(m_logger, "Found {} reachable species in adaptive engine view.", reachableSpecies.size()); - const std::vector finalReactions = cullReactionsByFlow(allFlows, reachableSpecies, Y_Full, maxFlow); + const std::vector finalReactions = cullReactionsByFlow(allFlows, reachableSpecies, composition, maxFlow); finalizeActiveSet(finalReactions); - auto [rescuedReactions, rescuedSpecies] = rescueEdgeSpeciesDestructionChannel(Y_Full, netIn.temperature/1e9, netIn.density, m_activeSpecies, m_activeReactions); + auto [rescuedReactions, rescuedSpecies] = rescueEdgeSpeciesDestructionChannel(composition, netIn.temperature/1e9, netIn.density, m_activeSpecies, m_activeReactions); for (const auto& reactionPtr : rescuedReactions) { m_activeReactions.add_reaction(*reactionPtr); @@ -145,59 +144,46 @@ namespace gridfire { } std::expected, expectations::StaleEngineError> AdaptiveEngineView::calculateRHSAndEnergy( - const std::vector &Y_culled, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { validateState(); - - const auto Y_full = mapCulledToFull(Y_culled); - - auto result = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho); + // TODO: Think about if I need to reach in and adjust the composition to zero out inactive species. + auto result = m_baseEngine.calculateRHSAndEnergy(comp, T9, rho); if (!result) { return std::unexpected{result.error()}; } - const auto [dydt, nuclearEnergyGenerationRate] = result.value(); - StepDerivatives culledResults; - culledResults.nuclearEnergyGenerationRate = nuclearEnergyGenerationRate; - culledResults.dydt = mapFullToCulled(dydt); - - return culledResults; + return result.value(); } EnergyDerivatives AdaptiveEngineView::calculateEpsDerivatives( - const std::vector &Y_culled, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { validateState(); - const auto Y_full = mapCulledToFull(Y_culled); - - return m_baseEngine.calculateEpsDerivatives(Y_full, T9, rho); + return m_baseEngine.calculateEpsDerivatives(comp, T9, rho); } void AdaptiveEngineView::generateJacobianMatrix( - const std::vector &Y_dynamic, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { validateState(); - const auto Y_full = mapCulledToFull(Y_dynamic); - - m_baseEngine.generateJacobianMatrix(Y_full, T9, rho); + m_baseEngine.generateJacobianMatrix(comp, T9, rho); } double AdaptiveEngineView::getJacobianMatrixEntry( - const int i_culled, - const int j_culled + const Species &rowSpecies, + const Species &colSpecies ) const { validateState(); - const size_t i_full = mapCulledToFullSpeciesIndex(i_culled); - const size_t j_full = mapCulledToFullSpeciesIndex(j_culled); - return m_baseEngine.getJacobianMatrixEntry(static_cast(i_full), static_cast(j_full)); + return m_baseEngine.getJacobianMatrixEntry(rowSpecies, colSpecies); } void AdaptiveEngineView::generateStoichiometryMatrix() { @@ -206,18 +192,16 @@ namespace gridfire { } int AdaptiveEngineView::getStoichiometryMatrixEntry( - const int speciesIndex_culled, - const int reactionIndex_culled + const Species &species, + const reaction::Reaction& reaction ) const { validateState(); - const size_t speciesIndex_full = mapCulledToFullSpeciesIndex(speciesIndex_culled); - const size_t reactionIndex_full = mapCulledToFullReactionIndex(reactionIndex_culled); - return m_baseEngine.getStoichiometryMatrixEntry(static_cast(speciesIndex_full), static_cast(reactionIndex_full)); + return m_baseEngine.getStoichiometryMatrixEntry(species, reaction); } double AdaptiveEngineView::calculateMolarReactionFlow( const reaction::Reaction &reaction, - const std::vector &Y_culled, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { @@ -227,9 +211,8 @@ namespace gridfire { m_logger -> flush_log(); throw std::runtime_error("Reaction not found in active reactions: " + std::string(reaction.id())); } - const auto Y = mapCulledToFull(Y_culled); - return m_baseEngine.calculateMolarReactionFlow(reaction, Y, T9, rho); + return m_baseEngine.calculateMolarReactionFlow(reaction, comp, T9, rho); } const reaction::ReactionSet & AdaptiveEngineView::getNetworkReactions() const { @@ -242,13 +225,12 @@ namespace gridfire { } std::expected, expectations::StaleEngineError> AdaptiveEngineView::getSpeciesTimescales( - const std::vector &Y_culled, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { validateState(); - const auto Y_full = mapCulledToFull(Y_culled); - const auto result = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho); + const auto result = m_baseEngine.getSpeciesTimescales(comp, T9, rho); if (!result) { return std::unexpected{result.error()}; @@ -270,15 +252,13 @@ namespace gridfire { std::expected, expectations::StaleEngineError> AdaptiveEngineView::getSpeciesDestructionTimescales( - const std::vector &Y, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { validateState(); - const std::vector Y_full = mapCulledToFull(Y); - - const auto result = m_baseEngine.getSpeciesDestructionTimescales(Y_full, T9, rho); + const auto result = m_baseEngine.getSpeciesDestructionTimescales(comp, T9, rho); if (!result) { return std::unexpected{result.error()}; } @@ -344,7 +324,7 @@ namespace gridfire { } size_t AdaptiveEngineView::mapCulledToFullSpeciesIndex(size_t culledSpeciesIndex) const { - if (culledSpeciesIndex < 0 || culledSpeciesIndex >= m_speciesIndexMap.size()) { + if (culledSpeciesIndex >= m_speciesIndexMap.size()) { LOG_ERROR(m_logger, "Culled index {} is out of bounds for species index map of size {}.", culledSpeciesIndex, m_speciesIndexMap.size()); m_logger->flush_log(); throw std::out_of_range("Culled index " + std::to_string(culledSpeciesIndex) + " is out of bounds for species index map of size " + std::to_string(m_speciesIndexMap.size()) + "."); @@ -353,7 +333,7 @@ namespace gridfire { } size_t AdaptiveEngineView::mapCulledToFullReactionIndex(size_t culledReactionIndex) const { - if (culledReactionIndex < 0 || culledReactionIndex >= m_reactionIndexMap.size()) { + if (culledReactionIndex >= m_reactionIndexMap.size()) { LOG_ERROR(m_logger, "Culled index {} is out of bounds for reaction index map of size {}.", culledReactionIndex, m_reactionIndexMap.size()); m_logger->flush_log(); throw std::out_of_range("Culled index " + std::to_string(culledReactionIndex) + " is out of bounds for reaction index map of size " + std::to_string(m_reactionIndexMap.size()) + "."); @@ -369,21 +349,17 @@ namespace gridfire { } } - // TODO: Change this to use a return value instead of an output parameter. - std::vector AdaptiveEngineView::calculateAllReactionFlows( - const NetIn &netIn, - std::vector &out_Y_Full + std::pair, fourdst::composition::Composition> AdaptiveEngineView::calculateAllReactionFlows( + const NetIn &netIn ) const { const auto& fullSpeciesList = m_baseEngine.getNetworkSpecies(); - out_Y_Full.clear(); - out_Y_Full.reserve(fullSpeciesList.size()); + fourdst::composition::Composition composition = netIn.composition; for (const auto& species: fullSpeciesList) { - if (netIn.composition.contains(species)) { - out_Y_Full.push_back(netIn.composition.getMolarAbundance(std::string(species.name()))); - } else { + if (!netIn.composition.contains(species)) { LOG_TRACE_L2(m_logger, "Species '{}' not found in composition. Setting abundance to 0.0.", species.name()); - out_Y_Full.push_back(0.0); + composition.registerSpecies(species); + composition.setMassFraction(species, 0.0); } } @@ -394,11 +370,11 @@ namespace gridfire { const auto& fullReactionSet = m_baseEngine.getNetworkReactions(); reactionFlows.reserve(fullReactionSet.size()); for (const auto& reaction : fullReactionSet) { - const double flow = m_baseEngine.calculateMolarReactionFlow(*reaction, out_Y_Full, T9, rho); + const double flow = m_baseEngine.calculateMolarReactionFlow(*reaction, composition, T9, rho); reactionFlows.push_back({reaction.get(), flow}); LOG_TRACE_L1(m_logger, "Reaction '{}' has flow rate: {:0.3E} [mol/s/g]", reaction->id(), flow); } - return reactionFlows; + return {reactionFlows, composition}; } std::unordered_set AdaptiveEngineView::findReachableSpecies( @@ -447,7 +423,7 @@ namespace gridfire { std::vector AdaptiveEngineView::cullReactionsByFlow( const std::vector &allFlows, const std::unordered_set &reachableSpecies, - const std::vector &Y_full, + const fourdst::composition::Composition &comp, const double maxFlow ) const { LOG_TRACE_L1(m_logger, "Culling reactions based on flow rates..."); @@ -464,9 +440,7 @@ namespace gridfire { bool zero_flow_due_to_reachable_reactants = false; if (flowRate < 1e-99 && flowRate > 0.0) { for (const auto& reactant: reactionPtr->reactants()) { - const auto it = std::ranges::find(m_baseEngine.getNetworkSpecies(), reactant); - const size_t index = std::distance(m_baseEngine.getNetworkSpecies().begin(), it); - if (Y_full[index] < 1e-99 && reachableSpecies.contains(reactant)) { + if (comp.getMolarAbundance(reactant) < 1e-99 && reachableSpecies.contains(reactant)) { LOG_TRACE_L1(m_logger, "Maintaining reaction '{}' with low flow ({:0.3E} [mol/s/g]) due to reachable reactant '{}'.", reactionPtr->id(), flowRate, reactant.name()); zero_flow_due_to_reachable_reactants = true; break; @@ -488,13 +462,13 @@ namespace gridfire { } AdaptiveEngineView::RescueSet AdaptiveEngineView::rescueEdgeSpeciesDestructionChannel( - const std::vector &Y_full, + const fourdst::composition::Composition &comp, const double T9, const double rho, const std::vector &activeSpecies, const reaction::ReactionSet &activeReactions ) const { - const auto result = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho); + const auto result = m_baseEngine.getSpeciesTimescales(comp, T9, rho); if (!result) { LOG_ERROR(m_logger, "Failed to get species timescales due to stale engine state."); throw exceptions::StaleEngineError("Failed to get species timescales"); @@ -565,8 +539,23 @@ namespace gridfire { allOtherReactantsAreAvailable = false; } } - if (allOtherReactantsAreAvailable && speciesToCheckIsConsumed) { - double rate = reaction->calculate_rate(T9, rho, Y_full); + if (allOtherReactantsAreAvailable) { + std::vector Y = comp.getMolarAbundanceVector(); + + const double Ye = comp.getElectronAbundance(); + // TODO: This is a dummy placeholder which must be replaced with an EOS call + const double mue = 5.0e-3 * std::pow(rho * Ye, 1.0 / 3.0) + 0.5 * T9; + + std::unordered_map speciesMassMap; + for (const auto &entry: comp | std::views::values) { + speciesMassMap[entry.isotope()] = entry.isotope().mass(); + } + std::unordered_map speciesIndexMap; + for (const auto& entry: comp | std::views::values) { + size_t distance = std::distance(speciesMassMap.begin(), speciesMassMap.find(entry.isotope())); + speciesIndexMap.emplace(distance, entry.isotope()); + } + double rate = reaction->calculate_rate(T9, rho, Ye, mue, Y, speciesIndexMap); if (rate > maxSpeciesConsumptionRate) { maxSpeciesConsumptionRate = rate; reactionsToRescue[species] = reaction.get(); diff --git a/src/lib/engine/views/engine_defined.cpp b/src/lib/engine/views/engine_defined.cpp index 5c30d628..f716980a 100644 --- a/src/lib/engine/views/engine_defined.cpp +++ b/src/lib/engine/views/engine_defined.cpp @@ -28,58 +28,48 @@ namespace gridfire { } std::expected, expectations::StaleEngineError> DefinedEngineView::calculateRHSAndEnergy( - const std::vector &Y_defined, + const fourdst::composition::Composition& comp, const double T9, const double rho ) const { validateNetworkState(); - const auto Y_full = mapViewToFull(Y_defined); - const auto result = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho); + const auto result = m_baseEngine.calculateRHSAndEnergy(comp, T9, rho); if (!result) { return std::unexpected{result.error()}; } - const auto [dydt, nuclearEnergyGenerationRate] = result.value(); - StepDerivatives definedResults; - definedResults.nuclearEnergyGenerationRate = nuclearEnergyGenerationRate; - definedResults.dydt = mapFullToView(dydt); - return definedResults; + return result.value(); } EnergyDerivatives DefinedEngineView::calculateEpsDerivatives( - const std::vector &Y_dynamic, + const fourdst::composition::Composition& comp, const double T9, const double rho ) const { validateNetworkState(); - const auto Y_full = mapViewToFull(Y_dynamic); - return m_baseEngine.calculateEpsDerivatives(Y_full, T9, rho); + return m_baseEngine.calculateEpsDerivatives(comp, T9, rho); } void DefinedEngineView::generateJacobianMatrix( - const std::vector &Y_dynamic, + const fourdst::composition::Composition& comp, const double T9, const double rho ) const { validateNetworkState(); - const auto Y_full = mapViewToFull(Y_dynamic); - m_baseEngine.generateJacobianMatrix(Y_full, T9, rho); + m_baseEngine.generateJacobianMatrix(comp, T9, rho); } double DefinedEngineView::getJacobianMatrixEntry( - const int i_defined, - const int j_defined + const Species& rowSpecies, + const Species& colSpecies ) const { validateNetworkState(); - const size_t i_full = mapViewToFullSpeciesIndex(i_defined); - const size_t j_full = mapViewToFullSpeciesIndex(j_defined); - - return m_baseEngine.getJacobianMatrixEntry(static_cast(i_full), static_cast(j_full)); + return m_baseEngine.getJacobianMatrixEntry(rowSpecies, colSpecies); } void DefinedEngineView::generateStoichiometryMatrix() { @@ -89,19 +79,17 @@ namespace gridfire { } int DefinedEngineView::getStoichiometryMatrixEntry( - const int speciesIndex_defined, - const int reactionIndex_defined + const Species& species, + const reaction::Reaction& reaction ) const { validateNetworkState(); - const size_t i_full = mapViewToFullSpeciesIndex(speciesIndex_defined); - const size_t j_full = mapViewToFullReactionIndex(reactionIndex_defined); - return m_baseEngine.getStoichiometryMatrixEntry(static_cast(i_full), static_cast(j_full)); + return m_baseEngine.getStoichiometryMatrixEntry(species, reaction); } double DefinedEngineView::calculateMolarReactionFlow( const reaction::Reaction &reaction, - const std::vector &Y_defined, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { @@ -112,8 +100,7 @@ namespace gridfire { m_logger -> flush_log(); throw std::runtime_error("Reaction not found in active reactions: " + std::string(reaction.id())); } - const auto Y_full = mapViewToFull(Y_defined); - return m_baseEngine.calculateMolarReactionFlow(reaction, Y_full, T9, rho); + return m_baseEngine.calculateMolarReactionFlow(reaction, comp, T9, rho); } const reaction::ReactionSet & DefinedEngineView::getNetworkReactions() const { @@ -131,14 +118,13 @@ namespace gridfire { } std::expected, expectations::StaleEngineError> DefinedEngineView::getSpeciesTimescales( - const std::vector &Y_defined, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { validateNetworkState(); - const auto Y_full = mapViewToFull(Y_defined); - const auto result = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho); + const auto result = m_baseEngine.getSpeciesTimescales(comp, T9, rho); if (!result) { return std::unexpected{result.error()}; } @@ -155,14 +141,13 @@ namespace gridfire { std::expected, expectations::StaleEngineError> DefinedEngineView::getSpeciesDestructionTimescales( - const std::vector &Y_defined, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { validateNetworkState(); - const auto Y_full = mapViewToFull(Y_defined); - const auto result = m_baseEngine.getSpeciesDestructionTimescales(Y_full, T9, rho); + const auto result = m_baseEngine.getSpeciesDestructionTimescales(comp, T9, rho); if (!result) { return std::unexpected{result.error()}; @@ -304,7 +289,7 @@ namespace gridfire { } size_t DefinedEngineView::mapViewToFullSpeciesIndex(size_t culledSpeciesIndex) const { - if (culledSpeciesIndex < 0 || culledSpeciesIndex >= m_speciesIndexMap.size()) { + if (culledSpeciesIndex >= m_speciesIndexMap.size()) { LOG_ERROR(m_logger, "Defined index {} is out of bounds for species index map of size {}.", culledSpeciesIndex, m_speciesIndexMap.size()); m_logger->flush_log(); throw std::out_of_range("Defined index " + std::to_string(culledSpeciesIndex) + " is out of bounds for species index map of size " + std::to_string(m_speciesIndexMap.size()) + "."); @@ -313,7 +298,7 @@ namespace gridfire { } size_t DefinedEngineView::mapViewToFullReactionIndex(size_t culledReactionIndex) const { - if (culledReactionIndex < 0 || culledReactionIndex >= m_reactionIndexMap.size()) { + if (culledReactionIndex >= m_reactionIndexMap.size()) { LOG_ERROR(m_logger, "Defined index {} is out of bounds for reaction index map of size {}.", culledReactionIndex, m_reactionIndexMap.size()); m_logger->flush_log(); throw std::out_of_range("Defined index " + std::to_string(culledReactionIndex) + " is out of bounds for reaction index map of size " + std::to_string(m_reactionIndexMap.size()) + "."); diff --git a/src/lib/engine/views/engine_multiscale.cpp b/src/lib/engine/views/engine_multiscale.cpp index 8f43c5a1..5dac52f2 100644 --- a/src/lib/engine/views/engine_multiscale.cpp +++ b/src/lib/engine/views/engine_multiscale.cpp @@ -4,15 +4,14 @@ #include #include +#include #include #include #include -#include #include -#include "gridfire/utils/logging.h" #include "quill/LogMacros.h" #include "quill/Logger.h" @@ -38,23 +37,23 @@ namespace { seed ^= hashed(v) + 0x9e3779b9 + (seed << 6) + (seed >> 2); } - std::vector> findConnectedComponentsBFS( - const std::unordered_map>& graph, - const std::vector& nodes + std::vector> findConnectedComponentsBFS( + const std::unordered_map>& graph, + const std::vector& nodes ) { - std::vector> components; - std::unordered_set visited; + std::vector> components; + std::unordered_set visited; - for (const size_t& start_node : nodes) { + for (const Species& start_node : nodes) { if (!visited.contains(start_node)) { - std::vector current_component; - std::queue q; + std::vector current_component; + std::queue q; q.push(start_node); visited.insert(start_node); while (!q.empty()) { - size_t u = q.front(); + Species u = q.front(); q.pop(); current_component.push_back(u); @@ -149,7 +148,6 @@ namespace { } namespace gridfire { - static int s_operator_parens_called = 0; using fourdst::atomic::Species; MultiscalePartitioningEngineView::MultiscalePartitioningEngineView( @@ -161,65 +159,53 @@ namespace gridfire { } std::expected, expectations::StaleEngineError> MultiscalePartitioningEngineView::calculateRHSAndEnergy( - const std::vector &Y_full, + const fourdst::composition::Composition& comp, const double T9, const double rho ) const { - if (Y_full.size() != getNetworkSpecies().size()) { - LOG_ERROR( - m_logger, - "Y_full size ({}) does not match the number of species in the network ({})", - Y_full.size(), - getNetworkSpecies().size() - ); - throw std::runtime_error("Y_full size does not match the number of species in the network. See logs for more details..."); - } - // Check the cache to see if the network needs to be repartitioned. Note that the QSECacheKey manages binning of T9, rho, and Y_full to ensure that small changes (which would likely not result in a repartitioning) do not trigger a cache miss. - const auto result = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho); + const auto result = m_baseEngine.calculateRHSAndEnergy(comp, T9, rho); if (!result) { return std::unexpected{result.error()}; } auto deriv = result.value(); - for (const size_t species_index : m_algebraic_species_indices) { - deriv.dydt[species_index] = 0.0; // Fix the algebraic species to the equilibrium abundances we calculate. + //TODO: Sort out how to deal with this. Need to return something like a step derivative but with the index consistent... + for (const auto& species : m_algebraic_species) { + deriv.dydt[species] = 0.0; // Fix the algebraic species to the equilibrium abundances we calculate. } return deriv; } EnergyDerivatives MultiscalePartitioningEngineView::calculateEpsDerivatives( - const std::vector &Y, + const fourdst::composition::Composition& comp, const double T9, const double rho ) const { - return m_baseEngine.calculateEpsDerivatives(Y, T9, rho); + return m_baseEngine.calculateEpsDerivatives(comp, T9, rho); } void MultiscalePartitioningEngineView::generateJacobianMatrix( - const std::vector &Y_full, + const fourdst::composition::Composition& comp, const double T9, const double rho ) const { // TODO: Add sparsity pattern to this to prevent base engine from doing unnecessary work. - m_baseEngine.generateJacobianMatrix(Y_full, T9, rho); + m_baseEngine.generateJacobianMatrix(comp, T9, rho); } double MultiscalePartitioningEngineView::getJacobianMatrixEntry( - const int i_full, - const int j_full + const Species& rowSpecies, + const Species& colSpecies ) 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. + if (std::ranges::contains(m_algebraic_species, colSpecies)) { return 0.0; } - if (std::ranges::contains(m_algebraic_species_indices, i_full)) { + if (std::ranges::contains(m_algebraic_species, rowSpecies)) { return 0.0; } - // Otherwise we need to query the full jacobian - return m_baseEngine.getJacobianMatrixEntry(i_full, j_full); + return m_baseEngine.getJacobianMatrixEntry(rowSpecies, colSpecies); } void MultiscalePartitioningEngineView::generateStoichiometryMatrix() { @@ -227,27 +213,31 @@ namespace gridfire { } int MultiscalePartitioningEngineView::getStoichiometryMatrixEntry( - const int speciesIndex, - const int reactionIndex + const Species& species, + const reaction::Reaction& reaction ) const { - return m_baseEngine.getStoichiometryMatrixEntry(speciesIndex, reactionIndex); + return m_baseEngine.getStoichiometryMatrixEntry(species, reaction); } double MultiscalePartitioningEngineView::calculateMolarReactionFlow( const reaction::Reaction &reaction, - const std::vector &Y_full, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { - assert(m_Y_algebraic.size() == m_algebraic_species_indices.size()); - // Fix the algebraic species to the equilibrium abundances we calculate. - std::vector Y_mutable = Y_full; - for (const auto& [index, Yi] : std::views::zip(m_algebraic_species_indices, m_Y_algebraic)) { - Y_mutable[index] = Yi; - + fourdst::composition::Composition comp_mutable = comp; + for (const auto& species : m_algebraic_species) { + // TODO: Check this conversion to mass fraction (also consider adding the ability to set molar abundance directly) + const double Yi = m_algebraic_abundances.at(species); + comp_mutable.setMassFraction(species, Yi * species.a() / (rho * 1e-3)); // Convert Yi (mol/g) to Xi (mass fraction) } - return m_baseEngine.calculateMolarReactionFlow(reaction, Y_mutable, T9, rho); + if (!comp_mutable.finalize()) { + LOG_ERROR(m_logger, "Failed to finalize composition after setting algebraic species abundances."); + m_logger->flush_log(); + throw std::runtime_error("Failed to finalize composition after setting algebraic species abundances."); + } + return m_baseEngine.calculateMolarReactionFlow(reaction, comp_mutable, T9, rho); } const reaction::ReactionSet & MultiscalePartitioningEngineView::getNetworkReactions() const { @@ -260,11 +250,11 @@ namespace gridfire { } std::expected, expectations::StaleEngineError> MultiscalePartitioningEngineView::getSpeciesTimescales( - const std::vector &Y, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { - const auto result = m_baseEngine.getSpeciesTimescales(Y, T9, rho); + const auto result = m_baseEngine.getSpeciesTimescales(comp, T9, rho); if (!result) { return std::unexpected{result.error()}; } @@ -277,11 +267,11 @@ namespace gridfire { std::expected, expectations::StaleEngineError> MultiscalePartitioningEngineView::getSpeciesDestructionTimescales( - const std::vector &Y, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { - const auto result = m_baseEngine.getSpeciesDestructionTimescales(Y, T9, rho); + const auto result = m_baseEngine.getSpeciesDestructionTimescales(comp, T9, rho); if (!result) { return std::unexpected{result.error()}; } @@ -298,13 +288,12 @@ namespace gridfire { NetIn baseUpdatedNetIn = netIn; baseUpdatedNetIn.composition = baseUpdatedComposition; const fourdst::composition::Composition equilibratedComposition = equilibrateNetwork(baseUpdatedNetIn); - std::vector Y_algebraic(m_algebraic_species_indices.size(), 0.0); - for (size_t i = 0; i < m_algebraic_species_indices.size(); ++i) { - const size_t species_index = m_algebraic_species_indices[i]; - Y_algebraic[i] = equilibratedComposition.getMolarAbundance(m_baseEngine.getNetworkSpecies()[species_index]); + std::unordered_map algebraicAbundances; + for (const auto& species : m_algebraic_species) { + algebraicAbundances[species] = equilibratedComposition.getMolarAbundance(species); } - m_Y_algebraic = std::move(Y_algebraic); + m_algebraic_abundances = std::move(algebraicAbundances); return equilibratedComposition; } @@ -335,13 +324,13 @@ namespace gridfire { return m_baseEngine; } - std::vector> MultiscalePartitioningEngineView::analyzeTimescalePoolConnectivity( - const std::vector> ×cale_pools, - const std::vector &Y, + std::vector> MultiscalePartitioningEngineView::analyzeTimescalePoolConnectivity( + const std::vector> ×cale_pools, + const fourdst::composition::Composition &comp, double T9, double rho ) const { - std::vector> final_connected_pools; + std::vector> final_connected_pools; for (const auto& pool : timescale_pools) { if (pool.empty()) { @@ -358,7 +347,7 @@ namespace gridfire { } void MultiscalePartitioningEngineView::partitionNetwork( - const std::vector &Y, + const fourdst::composition::Composition &comp, const double T9, const double rho ) { @@ -367,28 +356,25 @@ namespace gridfire { LOG_TRACE_L1(m_logger, "Clearing previous state..."); m_qse_groups.clear(); m_dynamic_species.clear(); - m_dynamic_species_indices.clear(); m_algebraic_species.clear(); - m_algebraic_species_indices.clear(); // --- Step 1. Identify distinct timescale regions --- LOG_TRACE_L1(m_logger, "Identifying fast reactions..."); - const std::vector> timescale_pools = partitionByTimescale(Y, T9, rho); + const std::vector> timescale_pools = partitionByTimescale(comp, T9, rho); LOG_TRACE_L1(m_logger, "Found {} timescale pools.", timescale_pools.size()); // --- Step 2. Select the mean slowest pool as the base dynamical group --- LOG_TRACE_L1(m_logger, "Identifying mean slowest pool..."); - const size_t mean_slowest_pool_index = identifyMeanSlowestPool(timescale_pools, Y, T9, rho); + const size_t mean_slowest_pool_index = identifyMeanSlowestPool(timescale_pools, comp, T9, rho); LOG_TRACE_L1(m_logger, "Mean slowest pool index: {}", mean_slowest_pool_index); // --- Step 3. Push the slowest pool into the dynamic species list --- - m_dynamic_species_indices = timescale_pools[mean_slowest_pool_index]; - for (const auto& index : m_dynamic_species_indices) { - m_dynamic_species.push_back(m_baseEngine.getNetworkSpecies()[index]); + for (const auto& slowSpecies : timescale_pools[mean_slowest_pool_index]) { + m_dynamic_species.push_back(slowSpecies); } // --- Step 4. Pack Candidate QSE Groups --- - std::vector> candidate_pools; + std::vector> candidate_pools; for (size_t i = 0; i < timescale_pools.size(); ++i) { if (i == mean_slowest_pool_index) continue; // Skip the slowest pool LOG_TRACE_L1(m_logger, "Group {} with {} species identified for potential QSE.", i, timescale_pools[i].size()); @@ -396,13 +382,13 @@ namespace gridfire { } LOG_TRACE_L1(m_logger, "Preforming connectivity analysis on timescale pools..."); - const std::vector> connected_pools = analyzeTimescalePoolConnectivity(candidate_pools, Y, T9, rho); + const std::vector> connected_pools = analyzeTimescalePoolConnectivity(candidate_pools, comp, T9, rho); LOG_TRACE_L1(m_logger, "Found {} connected pools (compared to {} timescale pools) for QSE analysis.", connected_pools.size(), timescale_pools.size()); // --- Step 5. Identify potential seed species for each candidate pool --- LOG_TRACE_L1(m_logger, "Identifying potential seed species for candidate pools..."); - const std::vector candidate_groups = constructCandidateGroups(connected_pools, Y, T9, rho); + const std::vector candidate_groups = constructCandidateGroups(connected_pools, comp, T9, rho); LOG_TRACE_L1(m_logger, "Found {} candidate QSE groups for further analysis", candidate_groups.size()); LOG_TRACE_L2( m_logger, @@ -413,17 +399,17 @@ namespace gridfire { 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) { + for (const auto& species : group.algebraic_species) { + ss << species.name(); + if (i < group.algebraic_species.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) { + for (const auto& species : group.seed_species) { + ss << species.name(); + if (i < group.seed_species.size() - 1) { ss << ", "; } i++; @@ -439,7 +425,7 @@ namespace gridfire { ); LOG_TRACE_L1(m_logger, "Validating candidate groups with flux analysis..."); - const auto [validated_groups, invalidate_groups] = validateGroupsWithFluxAnalysis(candidate_groups, Y, T9, rho); + const auto [validated_groups, invalidate_groups] = validateGroupsWithFluxAnalysis(candidate_groups, comp, T9, rho); LOG_TRACE_L1( m_logger, "Validated {} group(s) QSE groups. {}", @@ -465,8 +451,8 @@ 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]); + for (const auto& species : group.algebraic_species) { + m_dynamic_species.push_back(species); } } @@ -475,11 +461,9 @@ namespace gridfire { for (const auto& group : m_qse_groups) { // Add algebraic species to the algebraic set - for (const auto& index : group.algebraic_indices) { - if (std::ranges::find(m_algebraic_species_indices, index) == m_algebraic_species_indices.end()) { - m_algebraic_species.push_back(m_baseEngine.getNetworkSpecies()[index]); - m_algebraic_species_indices.push_back(index); - + for (const auto& species : group.algebraic_species) { + if (std::ranges::find(m_algebraic_species, species) == m_algebraic_species.end()) { + m_algebraic_species.push_back(species); } } } @@ -509,16 +493,15 @@ namespace gridfire { void MultiscalePartitioningEngineView::partitionNetwork( const NetIn &netIn ) { - const std::vector Y = packCompositionToVector(netIn.composition, m_baseEngine); const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9) const double rho = netIn.density; // Density in g/cm^3 - partitionNetwork(Y, T9, rho); + partitionNetwork(netIn.composition, T9, rho); } void MultiscalePartitioningEngineView::exportToDot( const std::string &filename, - const std::vector& Y, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { @@ -534,12 +517,12 @@ namespace gridfire { // --- 1. Pre-computation and Categorization --- // Categorize species into algebraic, seed, and core dynamic - std::unordered_set algebraic_indices; - std::unordered_set seed_indices; + std::unordered_set algebraic_species; + std::unordered_set seed_species; for (const auto& group : m_qse_groups) { if (group.is_in_equilibrium) { - algebraic_indices.insert(group.algebraic_indices.begin(), group.algebraic_indices.end()); - seed_indices.insert(group.seed_indices.begin(), group.seed_indices.end()); + algebraic_species.insert(group.algebraic_species.begin(), group.algebraic_species.end()); + seed_species.insert(group.seed_species.begin(), group.seed_species.end()); } } @@ -550,7 +533,7 @@ namespace gridfire { double max_log_flow = std::numeric_limits::lowest(); for (const auto& reaction : all_reactions) { - double flow = std::abs(m_baseEngine.calculateMolarReactionFlow(*reaction, Y, T9, rho)); + double flow = std::abs(m_baseEngine.calculateMolarReactionFlow(*reaction, comp, T9, rho)); reaction_flows.push_back(flow); if (flow > 1e-99) { // Avoid log(0) double log_flow = std::log10(flow); @@ -571,17 +554,16 @@ namespace gridfire { // Define all species nodes first, so they can be referenced by clusters and ranks later. dotFile << " // --- Species Nodes Definitions ---\n"; std::map> species_by_mass; - for (size_t i = 0; i < all_species.size(); ++i) { - const auto& species = all_species[i]; + for (const auto & species : all_species) { std::string fillcolor = "#f1f5f9"; // Default: Other/Uninvolved // Determine color based on category. A species can be a seed and also in the core dynamic group. // The more specific category (algebraic, then seed) takes precedence. - if (algebraic_indices.contains(i)) { + if (algebraic_species.contains(species)) { fillcolor = "#e0f2fe"; // Light Blue: Algebraic (in QSE) - } else if (seed_indices.contains(i)) { + } else if (seed_species.contains(species)) { fillcolor = "#a7f3d0"; // Light Green: Seed (Dynamic, feeds a QSE group) - } else if (std::ranges::contains(m_dynamic_species_indices, i)) { + } else if (std::ranges::contains(m_dynamic_species, species)) { fillcolor = "#dcfce7"; // Pale Green: Core Dynamic } dotFile << " \"" << species.name() << "\" [label=\"" << species.name() << "\", fillcolor=\"" << fillcolor << "\"];\n"; @@ -625,7 +607,7 @@ namespace gridfire { dotFile << " // --- QSE Group Clusters ---\n"; int group_counter = 0; for (const auto& group : m_qse_groups) { - if (!group.is_in_equilibrium || group.algebraic_indices.empty()) { + if (!group.is_in_equilibrium || group.algebraic_species.empty()) { continue; } dotFile << " subgraph cluster_qse_" << group_counter++ << " {\n"; @@ -640,11 +622,11 @@ namespace gridfire { dotFile << " color = \"#a7f3d0\";\n"; // Light green for seed species dotFile << " penwidth = 1.5;\n"; // Thinner border for seed cluster std::vector seed_node_ids; - seed_node_ids.reserve(group.seed_indices.size()); - for (const size_t species_idx : group.seed_indices) { + seed_node_ids.reserve(group.seed_species.size()); + for (const auto& species : group.seed_species) { std::stringstream ss; - ss << "node_" << group_counter << "_seed_" << species_idx; - dotFile << " " << ss.str() << " [label=\"" << all_species[species_idx].name() << "\"];\n"; + ss << "node_" << group_counter << "_seed_" << species.name(); + dotFile << " " << ss.str() << " [label=\"" << species.name() << "\"];\n"; seed_node_ids.push_back(ss.str()); } for (size_t i = 0; i < seed_node_ids.size() - 1; ++i) { @@ -657,11 +639,11 @@ namespace gridfire { dotFile << " color = \"#e0f2fe\";\n"; // Light blue for algebraic species dotFile << " penwidth = 1.5;\n"; // Thinner border for algebraic cluster std::vector algebraic_node_ids; - algebraic_node_ids.reserve(group.algebraic_indices.size()); - for (const size_t species_idx : group.algebraic_indices) { + algebraic_node_ids.reserve(group.algebraic_species.size()); + for (const Species& species : group.algebraic_species) { std::stringstream ss; - ss << "node_" << group_counter << "_algebraic_" << species_idx; - dotFile << " " << ss.str() << " [label=\"" << all_species[species_idx].name() << "\"];\n"; + ss << "node_" << group_counter << "_algebraic_" << species.name(); + dotFile << " " << ss.str() << " [label=\"" << species.name() << "\"];\n"; algebraic_node_ids.push_back(ss.str()); } // Make invisible edges between algebraic indices to keep them in top-down order @@ -756,69 +738,49 @@ namespace gridfire { } fourdst::composition::Composition MultiscalePartitioningEngineView::equilibrateNetwork( - const std::vector &Y, + const fourdst::composition::Composition &comp, const double T9, const double rho ) { - partitionNetwork(Y, T9, rho); - const std::vector Y_equilibrated = solveQSEAbundances(Y, T9, rho); - fourdst::composition::Composition composition; + partitionNetwork(comp, T9, rho); + fourdst::composition::Composition qseComposition = solveQSEAbundances(comp, T9, rho); - std::vector symbols; - symbols.reserve(m_baseEngine.getNetworkSpecies().size()); - for (const auto& species : m_baseEngine.getNetworkSpecies()) { - symbols.emplace_back(species.name()); - } - composition.registerSymbol(symbols); - - std::vector X; - X.reserve(Y_equilibrated.size()); - for (size_t i = 0; i < Y_equilibrated.size(); ++i) { - const double molarMass = m_baseEngine.getNetworkSpecies()[i].mass(); - X.push_back(Y_equilibrated[i] * molarMass); // Convert from molar abundance to mass fraction - } - - for (size_t i = 0; i < Y_equilibrated.size(); ++i) { - const auto& species = m_baseEngine.getNetworkSpecies()[i]; - if (X[i] < 0.0 && std::abs(X[i]) < 1e-20) { - composition.setMassFraction(std::string(species.name()), 0.0); // Avoid negative mass fractions - } else { - composition.setMassFraction(std::string(species.name()), X[i]); + for (const auto &symbol: qseComposition | std::views::keys) { + const double speciesMassFraction = qseComposition.getMassFraction(symbol); + if (speciesMassFraction < 0.0 && std::abs(speciesMassFraction) < 1e-20) { + qseComposition.setMassFraction(symbol, 0.0); // Avoid negative mass fractions } } - composition.finalize(true); + qseComposition.finalize(true); - return composition; + return qseComposition; } fourdst::composition::Composition MultiscalePartitioningEngineView::equilibrateNetwork( const NetIn &netIn ) { const PrimingReport primingReport = m_baseEngine.primeEngine(netIn); - const std::vector Y = packCompositionToVector(primingReport.primedComposition, m_baseEngine); const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9) const double rho = netIn.density; // Density in g/cm^3 - return equilibrateNetwork(Y, T9, rho); + return equilibrateNetwork(primingReport.primedComposition, T9, rho); } - size_t MultiscalePartitioningEngineView::getSpeciesIndex(const fourdst::atomic::Species &species) const { + size_t MultiscalePartitioningEngineView::getSpeciesIndex(const Species &species) const { return m_baseEngine.getSpeciesIndex(species); } - - std::vector> MultiscalePartitioningEngineView::partitionByTimescale( - const std::vector& Y_full, + std::vector> MultiscalePartitioningEngineView::partitionByTimescale( + const fourdst::composition::Composition &comp, const double T9, const double rho ) 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; + const auto result= m_baseEngine.getSpeciesDestructionTimescales(comp, T9, rho); + const auto netTimescale = m_baseEngine.getSpeciesTimescales(comp, T9, rho); if (!result) { LOG_ERROR(m_logger, "Failed to get species destruction timescales due to stale engine state"); m_logger->flush_log(); @@ -833,15 +795,15 @@ namespace gridfire { 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]); + std::vector> sorted_timescales; + for (const auto & all_specie : all_species) { + double timescale = all_timescales.at(all_specie); + double net_timescale = net_timescales.at(all_specie); 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); + LOG_TRACE_L3(m_logger, "Species {} has finite destruction timescale: destruction: {} s, net: {} s", all_specie.name(), timescale, net_timescale); + sorted_timescales.emplace_back(timescale, all_specie); } else { - LOG_TRACE_L3(m_logger, "Species {} has infinite or negative destruction timescale: destruction: {} s, net: {} s", all_species[i].name(), timescale, net_timescale); + LOG_TRACE_L3(m_logger, "Species {} has infinite or negative destruction timescale: destruction: {} s, net: {} s", all_specie.name(), timescale, net_timescale); } } @@ -853,7 +815,7 @@ namespace gridfire { } ); - std::vector> final_pools; + std::vector> final_pools; if (sorted_timescales.empty()) { return final_pools; } @@ -866,25 +828,25 @@ namespace gridfire { ABSOLUTE_QSE_TIMESCALE_THRESHOLD, ABSOLUTE_QSE_TIMESCALE_THRESHOLD / 3.156e7); LOG_TRACE_L1(m_logger, "Minimum gap threshold: {} orders of magnitude.", MIN_GAP_THRESHOLD); - std::vector dynamic_pool_indices; - std::vector> fast_candidates; + std::vector dynamic_pool_species; + std::vector> fast_candidates; // 1. First Pass: Absolute Timescale Cutoff - for (const auto& ts_pair : sorted_timescales) { - if (ts_pair.first > ABSOLUTE_QSE_TIMESCALE_THRESHOLD) { + for (const auto& [timescale, species] : sorted_timescales) { + if (timescale > ABSOLUTE_QSE_TIMESCALE_THRESHOLD) { LOG_TRACE_L3(m_logger, "Species {} with timescale {} is considered dynamic (slower than qse timescale threshold).", - all_species[ts_pair.second].name(), ts_pair.first); - dynamic_pool_indices.push_back(ts_pair.second); + species.name(), timescale); + dynamic_pool_species.push_back(species); } else { LOG_TRACE_L3(m_logger, "Species {} with timescale {} is a candidate fast species (faster than qse timescale threshold).", - all_species[ts_pair.second].name(), ts_pair.first); - fast_candidates.push_back(ts_pair); + species.name(), timescale); + fast_candidates.emplace_back(timescale, species); } } - if (!dynamic_pool_indices.empty()) { - LOG_TRACE_L1(m_logger, "Found {} dynamic species (slower than QSE timescale threshold).", dynamic_pool_indices.size()); - final_pools.push_back(dynamic_pool_indices); + if (!dynamic_pool_species.empty()) { + LOG_TRACE_L1(m_logger, "Found {} dynamic species (slower than QSE timescale threshold).", dynamic_pool_species.size()); + final_pools.push_back(dynamic_pool_species); } if (fast_candidates.empty()) { @@ -899,15 +861,15 @@ namespace gridfire { const double t2 = fast_candidates[i+1].first; if (std::log10(t1) - std::log10(t2) > MIN_GAP_THRESHOLD) { LOG_TRACE_L3(m_logger, "Detected gap between species {} (timescale {:0.2E}) and {} (timescale {:0.2E}).", - all_species[fast_candidates[i].second].name(), t1, - all_species[fast_candidates[i+1].second].name(), t2); + fast_candidates[i].second.name(), t1, + fast_candidates[i+1].second.name(), t2); split_points.push_back(i + 1); } } size_t last_split = 0; for (const size_t split : split_points) { - std::vector pool; + std::vector pool; for (size_t i = last_split; i < split; ++i) { pool.push_back(fast_candidates[i].second); } @@ -915,7 +877,7 @@ namespace gridfire { last_split = split; } - std::vector final_fast_pool; + std::vector final_fast_pool; for (size_t i = last_split; i < fast_candidates.size(); ++i) { final_fast_pool.push_back(fast_candidates[i].second); } @@ -928,8 +890,8 @@ namespace gridfire { for (const auto& pool : final_pools) { ss << "["; int ic = 0; - for (const auto& idx : pool) { - ss << all_species[idx].name(); + for (const auto& species : pool) { + ss << species.name(); if (ic < pool.size() - 1) { ss << ", "; } @@ -951,7 +913,7 @@ namespace gridfire { QSEGroup>> MultiscalePartitioningEngineView::validateGroupsWithFluxAnalysis( const std::vector &candidate_groups, - const std::vector &Y, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { constexpr double FLUX_RATIO_THRESHOLD = 5; @@ -959,28 +921,28 @@ namespace gridfire { 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 algebraic_group_members( + group.algebraic_species.begin(), + group.algebraic_species.end() ); - const std::unordered_set seed_group_members( - group.seed_indices.begin(), - group.seed_indices.end() + const std::unordered_set seed_group_members( + group.seed_species.begin(), + group.seed_species.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)); + const double flow = std::abs(m_baseEngine.calculateMolarReactionFlow(*reaction, comp, T9, rho)); if (flow == 0.0) { continue; // Skip reactions with zero flow } bool has_internal_algebraic_reactant = false; for (const auto& reactant : reaction->reactants()) { - if (algebraic_group_members.contains(m_baseEngine.getSpeciesIndex(reactant))) { + if (algebraic_group_members.contains(reactant)) { has_internal_algebraic_reactant = true; } } @@ -988,13 +950,13 @@ namespace gridfire { bool has_internal_algebraic_product = false; for (const auto& product : reaction->products()) { - if (algebraic_group_members.contains(m_baseEngine.getSpeciesIndex(product))) { + if (algebraic_group_members.contains(product)) { has_internal_algebraic_product = true; } } 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()); + LOG_TRACE_L3(m_logger, "{}: Skipping reaction {} as it has no internal algebraic species in reactants or products.", group.toString(), reaction->id()); continue; } @@ -1007,15 +969,14 @@ namespace gridfire { 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()); + if (algebraic_group_members.contains(species)) { + LOG_TRACE_L3(m_logger, "{}: Species {} is an algebraic participant in reaction {}.", group.toString(), 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()); + } else if (seed_group_members.contains(species)) { + LOG_TRACE_L3(m_logger, "{}: Species {} is a seed participant in reaction {}.", group.toString(), 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()); + LOG_TRACE_L3(m_logger, "{}: Species {} is an external participant in reaction {}.", group.toString(), species.name(), reaction->id()); external_participants++; } } @@ -1041,9 +1002,9 @@ namespace gridfire { [&]() -> 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) { + for (const auto& species: group.algebraic_species) { + ss << species.name(); + if (count < group.algebraic_species.size() - 1) { ss << ", "; } count++; @@ -1064,9 +1025,9 @@ namespace gridfire { [&]() -> 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) { + for (const auto& species : group.algebraic_species) { + ss << species.name(); + if (count < group.algebraic_species.size() - 1) { ss << ", "; } count++; @@ -1085,112 +1046,51 @@ namespace gridfire { return {validated_groups, invalidated_groups}; } - std::vector MultiscalePartitioningEngineView::solveQSEAbundances( - const std::vector &Y_full, + fourdst::composition::Composition MultiscalePartitioningEngineView::solveQSEAbundances( + const fourdst::composition::Composition &comp, const double T9, const double rho ) { LOG_TRACE_L1(m_logger, "Solving for QSE abundances..."); - auto Y_out = Y_full; - // Sort by timescale to solve fastest QSE groups first (these can seed slower groups) std::ranges::sort(m_qse_groups, [](const QSEGroup& a, const QSEGroup& b) { return a.mean_timescale < b.mean_timescale; }); + fourdst::composition::Composition outputComposition = comp; - for (const auto&[species_indices, is_in_equilibrium, algebraic_indices, seed_indices, mean_timescale] : m_qse_groups) { - if (!is_in_equilibrium || species_indices.empty()) { - LOG_TRACE_L1( - m_logger, - "Skipping QSE group with {} species ({} algebraic ({}), {} seeds ({})) as it is not in equilibrium.", - species_indices.size(), - algebraic_indices.size(), - [&]() -> std::string { - std::ostringstream os; - int count = 0; - for (const auto& idx : algebraic_indices) { - os << m_baseEngine.getNetworkSpecies()[idx].name(); - if (count < algebraic_indices.size() - 1) { - os << ", "; - } - count++; - } - return os.str(); - }(), - seed_indices.size(), - [&]() -> std::string { - std::ostringstream os; - int count = 0; - for (const auto& idx : seed_indices) { - os << m_baseEngine.getNetworkSpecies()[idx].name(); - if (count < seed_indices.size() - 1) { - os << ", "; - } - count++; - } - return os.str(); - }() - ); + for (const auto&[is_in_equilibrium, algebraic_species, seed_species, mean_timescale] : m_qse_groups) { + if (!is_in_equilibrium || (algebraic_species.empty() && seed_species.empty())) { continue; } - - LOG_TRACE_L1( - m_logger, - "Solving for QSE abundances for group with {} species ([{}] algebraic, [{}] seeds).", - species_indices.size(), - [&]() -> std::string { - std::stringstream ss; - int count = 0; - for (const auto& idx : algebraic_indices) { - ss << m_baseEngine.getNetworkSpecies()[idx].name(); - if (count < algebraic_indices.size() - 1) { - ss << ", "; - } - count++; - } - return ss.str(); - }(), - [&]() -> std::string { - std::stringstream ss; - int count = 0; - for (const auto& idx : seed_indices) { - ss << m_baseEngine.getNetworkSpecies()[idx].name(); - if (count < seed_indices.size() - 1) { - ss << ", "; - } - count++; - } - return ss.str(); - }() - ); - - std::vector qse_solve_indices; - std::vector seed_indices_vec; - - seed_indices_vec.reserve(species_indices.size()); - qse_solve_indices.reserve(species_indices.size()); - - for (size_t seed_idx : seed_indices) { - seed_indices_vec.emplace_back(seed_idx); + fourdst::composition::Composition normalized_composition = comp; + for (const auto& species: algebraic_species) { + if (!normalized_composition.contains(species)) { + normalized_composition.registerSpecies(species); + normalized_composition.setMassFraction(species, 0.0); + } + } + for (const auto& species: seed_species) { + if (!normalized_composition.contains(species)) { + normalized_composition.registerSpecies(species); + normalized_composition.setMassFraction(species, 0.0); + } } - for (size_t algebraic_idx : algebraic_indices) { - qse_solve_indices.emplace_back(algebraic_idx); - } - - if (qse_solve_indices.empty()) continue; - - Eigen::VectorXd Y_scale(qse_solve_indices.size()); - Eigen::VectorXd v_initial(qse_solve_indices.size()); - for (long i = 0; i < qse_solve_indices.size(); ++i) { + Eigen::VectorXd Y_scale(algebraic_species.size()); + Eigen::VectorXd v_initial(algebraic_species.size()); + long i = 0; + std::unordered_map species_to_index_map; + for (const auto& species : algebraic_species) { constexpr double abundance_floor = 1.0e-15; - const double initial_abundance = Y_full[qse_solve_indices[i]]; + const double initial_abundance = normalized_composition.getMolarAbundance(species); Y_scale(i) = std::max(initial_abundance, abundance_floor); v_initial(i) = std::asinh(initial_abundance / Y_scale(i)); // Scale the initial abundances using asinh + species_to_index_map.emplace(species, i); + i++; } - EigenFunctor functor(*this, qse_solve_indices, Y_full, T9, rho, Y_scale); + EigenFunctor functor(*this, algebraic_species, normalized_composition, T9, rho, Y_scale, species_to_index_map); Eigen::LevenbergMarquardt lm(functor); lm.parameters.ftol = 1.0e-10; lm.parameters.xtol = 1.0e-10; @@ -1200,53 +1100,40 @@ namespace gridfire { if (status <= 0 || status >= 4) { std::stringstream msg; - msg << "QSE solver failed with status: " << status << " for group with seed "; - if (seed_indices.size() == 1) { - msg << "nucleus " << m_baseEngine.getNetworkSpecies()[seed_indices_vec[0]].name(); - } else { - msg << "nuclei "; - // TODO: Refactor nice list printing into utility function somewhere - size_t counter = 0; - for (const auto& seed_idx : seed_indices) { - msg << m_baseEngine.getNetworkSpecies()[seed_idx].name(); - if (counter < seed_indices.size() - 2) { - msg << ", "; - } else if (counter == seed_indices.size() - 2) { - if (seed_indices.size() < 2) { - msg << " and "; - } else { - msg << ", and "; - } - } - ++counter; - } - } + //TODO: Add a better error message here and quill logging + msg << "QSE solver failed with status: " << status; throw std::runtime_error(msg.str()); } LOG_TRACE_L1(m_logger, "Minimization succeeded!"); Eigen::VectorXd Y_final_qse = Y_scale.array() * v_initial.array().sinh(); // Convert back to physical abundances using asinh scaling - for (long i = 0; i < qse_solve_indices.size(); ++i) { + i = 0; + for (const auto& species: algebraic_species) { LOG_TRACE_L1( m_logger, - "Species {} (index {}) started with abundance {} and ended with {}.", - m_baseEngine.getNetworkSpecies()[qse_solve_indices[i]].name(), - qse_solve_indices[i], - Y_full[qse_solve_indices[i]], + "During QSE solving species {} started with a molar abundance of {} and ended with an abundance of {}.", + species.name(), + normalized_composition.getMolarAbundance(species), Y_final_qse(i) ); - Y_out[qse_solve_indices[i]] = Y_final_qse(i); + //TODO: CHeck this conversion + double Xi = Y_final_qse(i) * species.mass(); // Convert from molar abundance to mass fraction + if (!outputComposition.contains(species)) { + outputComposition.registerSpecies(species); + } + outputComposition.setMassFraction(species, Xi); + i++; } } - return Y_out; + return outputComposition; } size_t MultiscalePartitioningEngineView::identifyMeanSlowestPool( - const std::vector> &pools, - const std::vector &Y, + const std::vector> &pools, + const fourdst::composition::Composition &comp, const double T9, const double rho ) const { - const auto& result = m_baseEngine.getSpeciesDestructionTimescales(Y, T9, rho); + const auto& result = m_baseEngine.getSpeciesDestructionTimescales(comp, T9, rho); if (!result) { LOG_ERROR(m_logger, "Failed to get species timescales due to stale engine state"); m_logger->flush_log(); @@ -1254,16 +1141,13 @@ namespace gridfire { } const std::unordered_map all_timescales = result.value(); - - const auto& all_species = m_baseEngine.getNetworkSpecies(); - size_t slowest_pool_index = 0; // Default to the first pool if no valid pool is found double slowest_mean_timescale = std::numeric_limits::min(); size_t count = 0; for (const auto& pool : pools) { double mean_timescale = 0.0; - for (const auto& species_idx : pool) { - const double timescale = all_timescales.at(all_species[species_idx]); + for (const auto& species : pool) { + const double timescale = all_timescales.at(species); mean_timescale += timescale; } mean_timescale = mean_timescale / static_cast(pool.size()); @@ -1272,8 +1156,8 @@ namespace gridfire { count, [&]() -> std::string { std::stringstream ss; size_t iCount = 0; - for (const auto& idx : pool) { - ss << all_species[idx].name() << ": " << all_timescales.at(all_species[idx]); + for (const auto& species : pool) { + ss << species.name() << ": " << all_timescales.at(species); if (iCount < pool.size() - 1) { ss << ", "; } @@ -1293,15 +1177,15 @@ namespace gridfire { return slowest_pool_index; } - std::unordered_map> MultiscalePartitioningEngineView::buildConnectivityGraph( - const std::vector &species_pool + std::unordered_map> MultiscalePartitioningEngineView::buildConnectivityGraph( + const std::vector &species_pool ) const { - std::unordered_map> connectivity_graph; - const std::set pool_set(species_pool.begin(), species_pool.end()); + std::unordered_map> connectivity_graph; + const std::set pool_set(species_pool.begin(), species_pool.end()); + const std::unordered_set pool_species = [&]() -> std::unordered_set { std::unordered_set result; - for (const auto& species_idx : species_pool) { - Species species = m_baseEngine.getNetworkSpecies()[species_idx]; + for (const auto& species : species_pool) { result.insert(species); } return result; @@ -1310,16 +1194,6 @@ namespace gridfire { std::map> speciesReactionMap; std::vector candidate_reactions; - auto getSpeciesIdx = [&](const std::vector &species) -> std::vector { - std::vector result; - result.reserve(species.size()); - for (const auto& s : species) { - size_t idx = m_baseEngine.getSpeciesIndex(s); - result.push_back(idx); - } - return result; - }; - for (const auto& reaction : m_baseEngine.getNetworkReactions()) { const std::vector &reactants = reaction->reactants(); const std::vector &products = reaction->products(); @@ -1329,19 +1203,22 @@ namespace gridfire { // Only consider reactions where at least one distinct reactant and product are in the pool if (has_distinct_reactant_and_product_species(pool_species, reactant_set, product_set)) { - std::vector involvedIDs = getSpeciesIdx(reactants); - std::vector involvedProducts = getSpeciesIdx(products); - involvedIDs.insert(involvedIDs.end(), involvedProducts.begin(), involvedProducts.end()); - std::set involvedSet(involvedIDs.begin(), involvedIDs.end()); + std::set involvedSet; + involvedSet.insert(reactants.begin(), reactants.end()); + involvedSet.insert(products.begin(), products.end()); - std::vector intersection; + std::vector intersection; intersection.reserve(involvedSet.size()); - std::ranges::set_intersection(pool_set, involvedSet, std::back_inserter(intersection)); + for (const auto& s : pool_species) { // Find intersection with pool species + if (involvedSet.contains(s)) { + intersection.push_back(s); + } + } // Add clique - for (const size_t& u : intersection) { - for (const size_t& v : intersection) { + for (const auto& u : intersection) { + for (const auto& v : intersection) { if (u != v) { // Avoid self-loops connectivity_graph[u].push_back(v); } @@ -1355,13 +1232,13 @@ namespace gridfire { } std::vector MultiscalePartitioningEngineView::constructCandidateGroups( - const std::vector> &candidate_pools, - const std::vector &Y, - const double T9, const double rho + const std::vector> &candidate_pools, + const fourdst::composition::Composition &comp, + const double T9, + const double rho ) const { - const auto& all_species = m_baseEngine.getNetworkSpecies(); const auto& all_reactions = m_baseEngine.getNetworkReactions(); - const auto& result = m_baseEngine.getSpeciesDestructionTimescales(Y, T9, rho); + const auto& result = m_baseEngine.getSpeciesDestructionTimescales(comp, T9, rho); if (!result) { LOG_ERROR(m_logger, "Failed to get species timescales due to stale engine state"); m_logger->flush_log(); @@ -1375,22 +1252,21 @@ namespace gridfire { // For each pool first identify all topological bridge connections std::vector> bridge_reactions; - for (const auto& species_idx : pool) { - Species ash = all_species[species_idx]; + for (const auto& ash: pool) { for (const auto& reaction : all_reactions) { if (reaction->contains(ash)) { // Check to make sure there is at least one reactant that is not in the pool // This lets seed nuclei bring mass into the QSE group. bool has_external_reactant = false; for (const auto& reactant : reaction->reactants()) { - if (std::ranges::find(pool, m_baseEngine.getSpeciesIndex(reactant)) == pool.end()) { + if (std::ranges::find(pool, reactant) == pool.end()) { has_external_reactant = true; LOG_TRACE_L3(m_logger, "Found external reactant {} in reaction {} for species {}.", reactant.name(), reaction->id(), ash.name()); break; // Found an external reactant, no need to check further } } if (has_external_reactant) { - double flow = std::abs(m_baseEngine.calculateMolarReactionFlow(*reaction, Y, T9, rho)); + double flow = std::abs(m_baseEngine.calculateMolarReactionFlow(*reaction, comp, T9, rho)); LOG_TRACE_L3(m_logger, "Found bridge reaction {} with flow {} for species {}.", reaction->id(), flow, ash.name()); bridge_reactions.emplace_back(reaction.get(), flow); } @@ -1418,26 +1294,24 @@ namespace gridfire { split_points.push_back(bridge_reactions.size() - 1); } - std::vector seed_indices; + std::vector seed_species; for (auto &reaction: bridge_reactions | std::views::keys) { for (const auto& fuel : reaction->reactants()) { - size_t fuel_idx = m_baseEngine.getSpeciesIndex(fuel); // Only add the fuel if it is not already in the pool - if (std::ranges::find(pool, fuel_idx) == pool.end()) { - seed_indices.push_back(fuel_idx); + if (std::ranges::find(pool, fuel) == pool.end()) { + seed_species.push_back(fuel); } } } - std::set all_indices(pool.begin(), pool.end()); - for (const auto& seed_idx : seed_indices) { - all_indices.insert(seed_idx); + std::set pool_species(pool.begin(), pool.end()); + for (const auto& species : seed_species) { + pool_species.insert(species); } - const std::set poolSet(pool.begin(), pool.end()); - const std::set seedSet(seed_indices.begin(), seed_indices.end()); + const std::set poolSet(pool.begin(), pool.end()); + const std::set seedSet(seed_species.begin(), seed_species.end()); double mean_timescale = 0.0; - for (const auto& pool_idx : poolSet) { - const auto& species = all_species[pool_idx]; + for (const auto& species : poolSet) { if (destruction_timescales.contains(species)) { mean_timescale += std::min(destruction_timescales.at(species), species.halfLife()); // Use the minimum of destruction timescale and half-life } else { @@ -1445,7 +1319,7 @@ namespace gridfire { } } mean_timescale /= static_cast(poolSet.size()); - QSEGroup qse_group(all_indices, false, poolSet, seedSet, mean_timescale); + QSEGroup qse_group(false, poolSet, seedSet, mean_timescale); candidate_groups.push_back(qse_group); } return candidate_groups; @@ -1453,43 +1327,54 @@ namespace gridfire { int MultiscalePartitioningEngineView::EigenFunctor::operator()(const InputType &v_qse, OutputType &f_qse) const { - s_operator_parens_called++; - std::vector y_trial = m_Y_full_initial; + fourdst::composition::Composition comp_trial = m_initial_comp; Eigen::VectorXd y_qse = m_Y_scale.array() * v_qse.array().sinh(); // Convert to physical abundances using asinh scaling - for (long i = 0; i < m_qse_solve_indices.size(); ++i) { - y_trial[m_qse_solve_indices[i]] = y_qse(i); + for (const auto& species: m_qse_solve_species) { + if (!comp_trial.contains(species)) { + comp_trial.registerSpecies(species); + } + comp_trial.setMassFraction(species, y_qse[static_cast(m_qse_solve_species_index_map.at(species))]); } - const auto result = m_view->getBaseEngine().calculateRHSAndEnergy(y_trial, m_T9, m_rho); + const auto result = m_view->getBaseEngine().calculateRHSAndEnergy(comp_trial, m_T9, m_rho); if (!result) { throw exceptions::StaleEngineError("Failed to calculate RHS and energy due to stale engine state"); } const auto&[dydt, nuclearEnergyGenerationRate] = result.value(); - f_qse.resize(static_cast(m_qse_solve_indices.size())); - for (long i = 0; i < m_qse_solve_indices.size(); ++i) { - f_qse(i) = dydt[m_qse_solve_indices[i]]; + f_qse.resize(static_cast(m_qse_solve_species.size())); + long i = 0; + // TODO: make sure that just counting up i is a valid approach, this is a possible place an indexing bug may have crept in + for (const auto& species: m_qse_solve_species) { + f_qse(i) = dydt.at(species); + i++; } return 0; // Success } int MultiscalePartitioningEngineView::EigenFunctor::df(const InputType &v_qse, JacobianType &J_qse) const { - std::vector y_trial = m_Y_full_initial; + fourdst::composition::Composition comp_trial = m_initial_comp; Eigen::VectorXd y_qse = m_Y_scale.array() * v_qse.array().sinh(); // Convert to physical abundances using asinh scaling - for (long i = 0; i < m_qse_solve_indices.size(); ++i) { - y_trial[m_qse_solve_indices[i]] = y_qse(i); + for (const auto& species: m_qse_solve_species) { + if (!comp_trial.contains(species)) { + comp_trial.registerSpecies(species); + } + comp_trial.setMassFraction(species, y_qse[static_cast(m_qse_solve_species_index_map.at(species))]); } - m_view->getBaseEngine().generateJacobianMatrix(y_trial, m_T9, m_rho); + m_view->getBaseEngine().generateJacobianMatrix(comp_trial, m_T9, m_rho); - J_qse.resize(static_cast(m_qse_solve_indices.size()), static_cast(m_qse_solve_indices.size())); - for (long i = 0; i < m_qse_solve_indices.size(); ++i) { - for (long j = 0; j < m_qse_solve_indices.size(); ++j) { - J_qse(i, j) = m_view->getBaseEngine().getJacobianMatrixEntry( - static_cast(m_qse_solve_indices[i]), - static_cast(m_qse_solve_indices[j]) + const long N = static_cast(m_qse_solve_species.size()); + J_qse.resize(N, N); + long rowID = 0; + long colID = 0; + for (const auto& rowSpecies : m_qse_solve_species) { + for (const auto& colSpecies: m_qse_solve_species) { + J_qse(rowID++, colID++) = m_view->getBaseEngine().getJacobianMatrixEntry( + rowSpecies, + colSpecies ); } } @@ -1564,18 +1449,18 @@ namespace gridfire { std::stringstream ss; ss << "QSEGroup(Algebraic: ["; size_t count = 0; - for (const auto& idx : algebraic_indices) { - ss << idx; - if (count < algebraic_indices.size() - 1) { + for (const auto& species : algebraic_species) { + ss << species.name(); + if (count < algebraic_species.size() - 1) { ss << ", "; } count++; } ss << "], Seed: ["; count = 0; - for (const auto& idx : seed_indices) { - ss << idx; - if (count < seed_indices.size() - 1) { + for (const auto& species : seed_species) { + ss << species.name(); + if (count < seed_species.size() - 1) { ss << ", "; } count++; @@ -1585,30 +1470,6 @@ namespace gridfire { } - std::string MultiscalePartitioningEngineView::QSEGroup::toString(const 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/io/network_file.cpp b/src/lib/io/network_file.cpp index 13081f30..d99a0c3a 100644 --- a/src/lib/io/network_file.cpp +++ b/src/lib/io/network_file.cpp @@ -39,7 +39,7 @@ namespace gridfire::io { } - SimpleReactionListFileParser::SimpleReactionListFileParser() {} + SimpleReactionListFileParser::SimpleReactionListFileParser() = default; ParsedNetworkData SimpleReactionListFileParser::parse(const std::string& filename) const { LOG_TRACE_L1(m_logger, "Parsing simple reaction list file: {}", filename); diff --git a/src/lib/network.cpp b/src/lib/network.cpp index 145e5a62..81eefbcd 100644 --- a/src/lib/network.cpp +++ b/src/lib/network.cpp @@ -19,7 +19,6 @@ // // *********************************************************************** */ #include "gridfire/network.h" -#include "gridfire/reaction/reaclib.h" #include "gridfire/reaction/reaction.h" #include @@ -74,7 +73,7 @@ namespace gridfire { } const auto ritr = std::find_if(str.rbegin(), std::string::const_reverse_iterator(startIt), [](const unsigned char ch){ return !std::isspace(ch); }); - return std::string(startIt, ritr.base()); + return {startIt, ritr.base()}; } } diff --git a/src/lib/partition/partition_rauscher_thielemann.cpp b/src/lib/partition/partition_rauscher_thielemann.cpp index b98bd20b..6c037a80 100644 --- a/src/lib/partition/partition_rauscher_thielemann.cpp +++ b/src/lib/partition/partition_rauscher_thielemann.cpp @@ -7,7 +7,6 @@ #include #include -#include #include #include diff --git a/src/lib/reaction/reaction.cpp b/src/lib/reaction/reaction.cpp index 8ef79aa1..d03a29db 100644 --- a/src/lib/reaction/reaction.cpp +++ b/src/lib/reaction/reaction.cpp @@ -45,15 +45,35 @@ namespace gridfire::reaction { m_rateCoefficients(sets), m_reverse(reverse) {} - double ReaclibReaction::calculate_rate(const double T9, const double rho, const std::vector& Y) const { + double ReaclibReaction::calculate_rate( + const double T9, + const double rho, + double Ye, + double mue, + const std::vector &Y, + const std::unordered_map& index_to_species_map + ) const { return calculate_rate(T9); } - CppAD::AD ReaclibReaction::calculate_rate(const CppAD::AD T9, const CppAD::AD rho, const std::vector>& Y) const { + CppAD::AD ReaclibReaction::calculate_rate( + const CppAD::AD T9, + const CppAD::AD rho, + CppAD::AD Ye, + CppAD::AD mue, + const std::vector>& Y, + const std::unordered_map& index_to_species_map + ) const { return calculate_rate>(T9); } - double ReaclibReaction::calculate_forward_rate_log_derivative(const double T9, const double rho, const std::vector& Y) const { + double ReaclibReaction::calculate_forward_rate_log_derivative( + const double T9, + const double rho, + double Ye, + double mue, + const fourdst::composition::Composition& comp + ) const { constexpr double r_p13 = 1.0 / 3.0; constexpr double r_p53 = 5.0 / 3.0; constexpr double r_p23 = 2.0 / 3.0; @@ -80,21 +100,11 @@ namespace gridfire::reaction { bool ReaclibReaction::contains_reactant(const Species& species) const { - for (const auto& reactant : m_reactants) { - if (reactant == species) { - return true; - } - } - return false; + return std::ranges::any_of(m_reactants, [&](const Species& r) { return r == species; }); } bool ReaclibReaction::contains_product(const Species& species) const { - for (const auto& product : m_products) { - if (product == species) { - return true; - } - } - return false; + return std::ranges::any_of(m_products, [&](const Species& p) { return p == species; }); } std::unordered_set ReaclibReaction::all_species() const { @@ -224,11 +234,20 @@ namespace gridfire::reaction { m_rates.push_back(reaction.rateCoefficients()); } - double LogicalReaclibReaction::calculate_rate(const double T9, const double rho, const std::vector& Y) const { + double LogicalReaclibReaction::calculate_rate( + const double T9, + const double rho, + double Ye, + double mue, const std::vector &Y, const std::unordered_map& index_to_species_map + ) const { return calculate_rate(T9); } - double LogicalReaclibReaction::calculate_forward_rate_log_derivative(const double T9, const double rho, const std::vector& Y) const { + double LogicalReaclibReaction::calculate_forward_rate_log_derivative( + const double T9, + const double rho, + double Ye, double mue, const fourdst::composition::Composition& comp + ) const { constexpr double r_p13 = 1.0 / 3.0; constexpr double r_p53 = 5.0 / 3.0; constexpr double r_p23 = 2.0 / 3.0; @@ -286,7 +305,8 @@ namespace gridfire::reaction { CppAD::AD LogicalReaclibReaction::calculate_rate( const CppAD::AD T9, const CppAD::AD rho, - const std::vector>& Y + CppAD::AD Ye, + CppAD::AD mue, const std::vector>& Y, const std::unordered_map& index_to_species_map ) const { return calculate_rate>(T9); } diff --git a/src/lib/reaction/weak/weak.cpp b/src/lib/reaction/weak/weak.cpp index 8c202be6..7f5c1085 100644 --- a/src/lib/reaction/weak/weak.cpp +++ b/src/lib/reaction/weak/weak.cpp @@ -7,20 +7,81 @@ #include #include #include +#include +#include -#define GRIDFIRE_WEAK_REACTION_LIB_SENTINEL -60.0 +#include "gridfire/reaction/weak/weak_interpolator.h" + +#include "xxhash64.h" + + +namespace { + fourdst::atomic::Species resolve_weak_product( + const gridfire::rates::weak::WeakReactionType type, + const fourdst::atomic::Species& reactant + ) { + using namespace fourdst::atomic; + using namespace gridfire::rates::weak; + + std::optional product; // Use optional so that we can start in a valid "null" state + switch (type) { + case WeakReactionType::BETA_MINUS_DECAY: + product = az_to_species(reactant.a(), reactant.z() + 1); + return product.value(); + case WeakReactionType::BETA_PLUS_DECAY: + product = az_to_species(reactant.a(), reactant.z() - 1); + return product.value(); + case WeakReactionType::ELECTRON_CAPTURE: + product = az_to_species(reactant.a(), reactant.z() - 1); + return product.value(); + case WeakReactionType::POSITRON_CAPTURE: + product = az_to_species(reactant.a(), reactant.z() + 1); + break; + } + if (!product.has_value()) { + throw std::runtime_error("Failed to resolve weak reaction product for reactant: " + std::string(reactant.name())); + } + return product.value(); + } + + std::string resolve_weak_id( + const gridfire::rates::weak::WeakReactionType type, + const fourdst::atomic::Species& reactant, + const fourdst::atomic::Species& product + ) { + using namespace gridfire::rates::weak; + + std::string id; + switch (type) { + case WeakReactionType::BETA_MINUS_DECAY: + id = std::format("{}(,ν|)e-,{}", reactant.name(), product.name()); + break; + case WeakReactionType::BETA_PLUS_DECAY: + id = std::format("{}(,ν)e+,{}", reactant.name(), product.name()); + break; + case WeakReactionType::ELECTRON_CAPTURE: + id = std::format("{}(e-,ν){}", reactant.name(), product.name()); + break; + case WeakReactionType::POSITRON_CAPTURE: + id = std::format("{}(e+,ν|){}", reactant.name(), product.name()); + break; + } + return id; + } +} namespace gridfire::rates::weak { WeakReactionMap::WeakReactionMap() { using namespace fourdst::atomic; + // ReSharper disable once CppUseStructuredBinding for (const auto& weak_reaction_record : UNIFIED_WEAK_DATA) { Species species = az_to_species(weak_reaction_record.A, weak_reaction_record.Z); if (weak_reaction_record.log_beta_minus > GRIDFIRE_WEAK_REACTION_LIB_SENTINEL) { m_weak_network[species].push_back( - WeakReaction{ + WeakReactionEntry{ WeakReactionType::BETA_MINUS_DECAY, weak_reaction_record.t9, weak_reaction_record.log_rhoye, @@ -32,7 +93,7 @@ namespace gridfire::rates::weak { } if (weak_reaction_record.log_beta_plus > GRIDFIRE_WEAK_REACTION_LIB_SENTINEL) { m_weak_network[species].push_back( - WeakReaction{ + WeakReactionEntry{ WeakReactionType::BETA_PLUS_DECAY, weak_reaction_record.t9, weak_reaction_record.log_rhoye, @@ -44,7 +105,7 @@ namespace gridfire::rates::weak { } if (weak_reaction_record.log_electron_capture > GRIDFIRE_WEAK_REACTION_LIB_SENTINEL) { m_weak_network[species].push_back( - WeakReaction{ + WeakReactionEntry{ WeakReactionType::ELECTRON_CAPTURE, weak_reaction_record.t9, weak_reaction_record.log_rhoye, @@ -56,7 +117,7 @@ namespace gridfire::rates::weak { } if (weak_reaction_record.log_positron_capture > GRIDFIRE_WEAK_REACTION_LIB_SENTINEL) { m_weak_network[species].push_back( - WeakReaction{ + WeakReactionEntry{ WeakReactionType::POSITRON_CAPTURE, weak_reaction_record.t9, weak_reaction_record.log_rhoye, @@ -69,26 +130,391 @@ namespace gridfire::rates::weak { } } - std::vector WeakReactionMap::get_all_reactions() const { - std::vector reactions; + std::vector WeakReactionMap::get_all_reactions() const { + std::vector reactions; for (const auto &species_reactions: m_weak_network | std::views::values) { reactions.insert(reactions.end(), species_reactions.begin(), species_reactions.end()); } return reactions; } - std::expected, bool> WeakReactionMap::get_species_reactions(const fourdst::atomic::Species &species) const { + std::expected, WeakMapError> WeakReactionMap::get_species_reactions( + const fourdst::atomic::Species &species) const { if (m_weak_network.contains(species)) { return m_weak_network.at(species); } - return std::unexpected(false); + return std::unexpected(WeakMapError::SPECIES_NOT_FOUND); } - std::expected, bool> WeakReactionMap::get_species_reactions(const std::string &species_name) const { - fourdst::atomic::Species species = fourdst::atomic::species.at(species_name); + std::expected, WeakMapError> WeakReactionMap::get_species_reactions( + const std::string &species_name) const { + const fourdst::atomic::Species species = fourdst::atomic::species.at(species_name); if (m_weak_network.contains(species)) { return m_weak_network.at(species); } - return std::unexpected(false); + return std::unexpected(WeakMapError::SPECIES_NOT_FOUND); } + + WeakReaction::WeakReaction( + const fourdst::atomic::Species &species, + const WeakReactionType type, + const WeakRateInterpolator &interpolator + ) : + m_reactant(species), + m_product(resolve_weak_product(type, species)), + m_reactant_a(species.a()), + m_reactant_z(species.z()), + m_product_a(m_product.a()), + m_product_z(m_product.z()), + m_id(resolve_weak_id(type, species, m_product)), + m_type(type), + m_interpolator(interpolator), + m_atomic(m_interpolator, m_reactant_a, m_reactant_z, m_type) {} + + double WeakReaction::calculate_rate( + const double T9, + const double rho, + const double Ye, + const double mue, + const std::vector &Y, + const std::unordered_map& index_to_species_map + ) const { + return calculate_rate(T9, rho, Ye, mue, Y, index_to_species_map); + } + + CppAD::AD WeakReaction::calculate_rate( + CppAD::AD T9, + CppAD::AD rho, + CppAD::AD Ye, + CppAD::AD mue, const std::vector> &Y, const std::unordered_map& index_to_species_map + ) const { + return static_cast>(0.0); + } + + bool WeakReaction::contains(const fourdst::atomic::Species &species) const { + return contains_reactant(species) || contains_product(species); + } + + bool WeakReaction::contains_reactant(const fourdst::atomic::Species& species) const { + if (m_reactant == species) { + return true; + } + return false; + } + + bool WeakReaction::contains_product(const fourdst::atomic::Species &species) const { + if (m_product == species) { + return true; + } + return false; + } + + std::unordered_set WeakReaction::all_species() const { + return {m_reactant, m_product}; + } + + std::unordered_set WeakReaction::reactant_species() const { + return {m_reactant}; + } + + std::unordered_set WeakReaction::product_species() const { + return {m_product}; + } + + int WeakReaction::stoichiometry(const fourdst::atomic::Species &species) const { + if (species == m_reactant) { + return -1; + } + if (species == m_product) { + return 1; + } + return 0; + } + + std::unordered_map WeakReaction::stoichiometry() const { + return { + {m_reactant, -1}, + {m_product, 1} + }; + } + + uint64_t WeakReaction::hash(const uint64_t seed) const { + const std::string reaction_string = std::format( + "{}:{}({})", + m_reactant.name(), + m_product.name(), + static_cast(m_type) + ); + return XXHash64::hash(reaction_string.data(), reaction_string.size(), seed); + } + + double WeakReaction::qValue() const { + // We ignore neutrino mass as it is negligible compared to other masses here. + double Q_MeV = 0.0; + + const double parentMass_u = m_reactant.mass(); + const double daughterMass_u = m_product.mass(); + const double electronMass_MeV = m_constants.electronMassMeV; + + const double nuclearMassDiff_MeV = (parentMass_u - daughterMass_u) * m_constants.u_to_MeV; + switch (m_type) { + case WeakReactionType::BETA_PLUS_DECAY: + Q_MeV = nuclearMassDiff_MeV - 2.0 * electronMass_MeV; + break; + case WeakReactionType::POSITRON_CAPTURE: + Q_MeV = nuclearMassDiff_MeV + 2.0 * electronMass_MeV; + break; + case WeakReactionType::BETA_MINUS_DECAY: + Q_MeV = nuclearMassDiff_MeV; + break; + case WeakReactionType::ELECTRON_CAPTURE: + Q_MeV = nuclearMassDiff_MeV; + break; + } + return Q_MeV; + } + + double WeakReaction::calculate_energy_generation_rate( + const double T9, + const double rho, + const double Ye, + const double mue, + const std::vector &Y, + const std::unordered_map &index_to_species_map + ) const { + std::expected rates = m_interpolator.get_rates( + static_cast(m_reactant_a), + static_cast(m_reactant_z), + T9, + std::log10(rho * Ye), + mue + ); + + if (!rates.has_value()) { + const InterpolationErrorType type = rates.error().type; + const std::string msg = std::format( + "Failed to interpolate weak rate for (A={}, Z={}) at T9={}, log10(rho*Ye)={}, mu_e={} with error: {}", + m_reactant.name(), m_reactant_a, m_reactant_z, T9, std::log10(rho * Ye), mue, InterpolationErrorTypeMap.at(type) + ); + throw std::runtime_error(msg); + } + + double logNeutrinoLossRate = 0.0; + if (m_type == WeakReactionType::BETA_PLUS_DECAY || m_type == WeakReactionType::ELECTRON_CAPTURE) { + logNeutrinoLossRate = rates->log_antineutrino_loss_bd; + } else if (m_type == WeakReactionType::BETA_MINUS_DECAY || m_type == WeakReactionType::POSITRON_CAPTURE) { + logNeutrinoLossRate = rates->log_neutrino_loss_ec; + } + + const double neutrinoLossRate = std::pow(10, logNeutrinoLossRate); + + const double EDeposited_MeV = qValue() - neutrinoLossRate; + + // We reimplement this logic here instead of calling calculate_rate() to avoid + // doing the interpolation twice (since the payload has already been interpolated). + const double logRate = get_log_rate_from_payload(rates.value()); + double lambda = 0.0; + if (logRate > GRIDFIRE_WEAK_REACTION_LIB_SENTINEL) { + lambda = std::pow(10, logRate); + } + return lambda * EDeposited_MeV; // returns in MeV / s + + } + + CppAD::AD WeakReaction::calculate_energy_generation_rate( + const CppAD::AD &T9, + const CppAD::AD &rho, + const CppAD::AD &Ye, + const CppAD::AD &mue, + const std::vector> &Y, + const std::unordered_map &index_to_species_map + ) const { + const CppAD::AD log_rhoYe = CppAD::log10(rho * Ye); + std::vector> ax = {T9, log_rhoYe, mue}; + std::vector> ay(1); + m_atomic(ax, ay); // TODO: Sort out why this isn't working and checkline 222 in weak.h where a similar line is + //TODO: think about how to get out neutrino loss in a autodiff safe way. This may mean I need to add an extra output to the atomic base + // so that I can get out both the rate and the neutrino loss rate. This will also mean that the sparsity pattern will need to + // be updated to account for the extra output. + CppAD::AD rateConstant = ay[0]; + + } + + std::unique_ptr WeakReaction::clone() const { + std::unique_ptr reaction_ptr = std::make_unique( + m_reactant, + m_type, + m_interpolator + ); + return reaction_ptr; + } + + double WeakReaction::get_log_rate_from_payload(const WeakRatePayload &payload) const { + double logRate = 0.0; + switch (m_type) { + case WeakReactionType::BETA_MINUS_DECAY: + logRate = payload.log_beta_minus; + break; + case WeakReactionType::BETA_PLUS_DECAY: + logRate = payload.log_beta_plus; + break; + case WeakReactionType::ELECTRON_CAPTURE: + logRate = payload.log_electron_capture; + break; + case WeakReactionType::POSITRON_CAPTURE: + logRate = payload.log_positron_capture; + break; + } + return logRate; + } + + bool WeakReaction::AtomicWeakRate::forward ( + const size_t p, + const size_t q, + const CppAD::vector &vx, + CppAD::vector &vy, + const CppAD::vector &tx, + CppAD::vector &ty + ) { + // Doing this explicitly (only allowing p == 0) makes forward mode AD impossible for now. + if (p != 0) { + return false; + } + const double T9 = tx[0]; + const double log10_rhoye = tx[1]; + const double mu_e = tx[2]; + + const std::expected result = m_interpolator.get_rates( + static_cast(m_a), + static_cast(m_z), + T9, + log10_rhoye, + mu_e + ); + if (!result.has_value()) { + const InterpolationErrorType type = result.error().type; + std::string msg = std::format( + "Failed to interpolate weak rate for (A={}, Z={}) at T9={}, log10(rho*Ye)={}, mu_e={} with error: {}", + m_a, m_z, T9, log10_rhoye, mu_e, InterpolationErrorTypeMap.at(type) + ); + } + switch (m_type) { + case WeakReactionType::BETA_MINUS_DECAY: + ty[0] = std::pow(10, result.value().log_beta_minus); + break; + case WeakReactionType::BETA_PLUS_DECAY: + ty[0] = std::pow(10, result.value().log_beta_plus); + break; + case WeakReactionType::ELECTRON_CAPTURE: + ty[0] = std::pow(10, result.value().log_electron_capture); + break; + case WeakReactionType::POSITRON_CAPTURE: + ty[0] = std::pow(10, result.value().log_positron_capture); + break; + } + + if (vx.size() > 0) { + vy[0] = vx[0] || vx[1] || vx[2]; // Sets the output sparsity pattern + } + return true; + } + + bool WeakReaction::AtomicWeakRate::reverse( + size_t q, + const CppAD::vector &tx, + const CppAD::vector &ty, + CppAD::vector &px, + const CppAD::vector &py + ) { + const double T9 = tx[0]; + const double log10_rhoye = tx[1]; + const double mu_e = tx[2]; + + const std::expected result = m_interpolator.get_rate_derivatives( + static_cast(m_a), + static_cast(m_z), + T9, + log10_rhoye, + mu_e + ); + + if (!result.has_value()) { + const InterpolationErrorType type = result.error().type; + const std::string msg = std::format( + "Failed to interpolate weak rate derivatives for (A={}, Z={}) at T9={}, log10(rho*Ye)={}, mu_e={} with error: {}", + m_a, m_z, T9, log10_rhoye, mu_e, InterpolationErrorTypeMap.at(type) + ); + throw std::runtime_error(msg); + } + + WeakRateDerivatives derivatives = result.value(); + + double dT9 = 0.0; + double dRho = 0.0; + double dMuE = 0.0; + switch (m_type) { + case WeakReactionType::BETA_MINUS_DECAY: + dT9 = py[0] * derivatives.d_log_beta_minus[0]; + dRho = py[0] * derivatives.d_log_beta_minus[1]; + dMuE = py[0] * derivatives.d_log_beta_minus[2]; + break; + case WeakReactionType::BETA_PLUS_DECAY: + dT9 = py[0] * derivatives.d_log_beta_plus[0]; + dRho = py[0] * derivatives.d_log_beta_plus[1]; + dMuE = py[0] * derivatives.d_log_beta_plus[2]; + break; + case WeakReactionType::ELECTRON_CAPTURE: + dT9 = py[0] * derivatives.d_log_electron_capture[0]; + dRho = py[0] * derivatives.d_log_electron_capture[1]; + dMuE = py[0] * derivatives.d_log_electron_capture[2]; + break; + case WeakReactionType::POSITRON_CAPTURE: + dT9 = py[0] * derivatives.d_log_positron_capture[0]; + dRho = py[0] * derivatives.d_log_positron_capture[1]; + dMuE = py[0] * derivatives.d_log_positron_capture[2]; + break; + } + + px[0] = py[0] * dT9; // d(rate)/dT9 + px[1] = py[0] * dRho; // d(rate)/dlogRhoYe + px[2] = py[0] * dMuE; // d(rate)/dMuE + + return true; + + } + + bool WeakReaction::AtomicWeakRate::for_sparse_jac( + size_t q, + const CppAD::vector > &r, + CppAD::vector > &s + ) { + s[0] = r[0]; + s[0].insert(r[1].begin(), r[1].end()); + s[0].insert(r[2].begin(), r[2].end()); + + return true; + } + + bool WeakReaction::AtomicWeakRate::rev_sparse_jac( + size_t q, + const CppAD::vector > &rt, + CppAD::vector > &st + ) { + // What this is saying is that each of the three input variables (T9, rho, Ye) + // all only affect the output variable (the rate) since there is only + // one output variable. + st[0] = rt[0]; + st[1] = rt[0]; + st[2] = rt[0]; + return true; + + } + + + + + + + + } diff --git a/src/lib/reaction/weak/weak_interpolator.cpp b/src/lib/reaction/weak/weak_interpolator.cpp new file mode 100644 index 00000000..a96c7a6f --- /dev/null +++ b/src/lib/reaction/weak/weak_interpolator.cpp @@ -0,0 +1,308 @@ +#include "gridfire/reaction/weak/weak_interpolator.h" +#include "gridfire/reaction/reaction.h" +#include "gridfire/reaction/weak/weak.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fourdst/composition/species.h" + +namespace gridfire::rates::weak { + + WeakRateInterpolator::WeakRateInterpolator(const RowDataTable &raw_data) { + std::map> grouped_rows; + for (const auto& row : raw_data) { + grouped_rows[pack_isotope_id(row.A, row.Z)].push_back(&row); + } + + for (auto const& [isotope_id, rows] : grouped_rows) { + IsotopeGrid grid; + + std::set unique_t9, unique_rhoYe, unique_mue; + for (const auto* row : rows) { + unique_t9.emplace(row->t9); + unique_rhoYe.emplace(row->log_rhoye); + unique_mue.emplace(row->mu_e); + } + + grid.t9_axis.reserve(unique_t9.size()); + grid.rhoYe_axis.reserve(unique_rhoYe.size()); + grid.mue_axis.reserve(unique_mue.size()); + + grid.t9_axis.insert(grid.t9_axis.begin(), unique_t9.begin(), unique_t9.end()); + grid.rhoYe_axis.insert(grid.rhoYe_axis.begin(), unique_rhoYe.begin(), unique_rhoYe.end()); + grid.mue_axis.insert(grid.mue_axis.begin(), unique_mue.begin(), unique_mue.end()); + + std::ranges::sort(grid.t9_axis); + std::ranges::sort(grid.rhoYe_axis); + std::ranges::sort(grid.mue_axis); + + const size_t nt9 = grid.t9_axis.size(); + const size_t nrhoYe = grid.rhoYe_axis.size(); + const size_t nmue = grid.mue_axis.size(); + + grid.data.resize(nt9 * nrhoYe * nmue); + + // Reverse map for quick index lookup + std::unordered_map t9_map, rhoYe_map, mue_map; + for (size_t i = 0; i < nt9; i++) { t9_map[grid.t9_axis[i]] = i; } + for (size_t j = 0; j < nrhoYe; j++) { rhoYe_map[grid.rhoYe_axis[j]] = j; } + for (size_t k = 0; k < nmue; k++) { mue_map[grid.mue_axis[k]] = k; } + + for (const auto* row: rows) { + size_t i_t9 = t9_map.at(row->t9); + size_t j_rhoYe = rhoYe_map.at(row->log_rhoye); + size_t k_mue = mue_map.at(row->mu_e); + + size_t index = (i_t9 * nrhoYe + j_rhoYe) * nmue + k_mue; + grid.data[index] = WeakRatePayload{ + row->log_beta_plus, + row->log_electron_capture, + row->log_neutrino_loss_ec, + row->log_beta_minus, + row->log_positron_capture, + row->log_antineutrino_loss_bd + }; + } + m_rate_table[isotope_id] = std::move(grid); + } + } + + std::vector WeakRateInterpolator::available_isotopes() const { + std::vector isotopes; + for (const auto &packed_id: m_rate_table | std::views::keys) { + const uint16_t A = static_cast(packed_id >> 8); + const uint8_t Z = static_cast(packed_id & 0xFF); + try { + fourdst::atomic::Species species = fourdst::atomic::az_to_species(A, Z); + isotopes.push_back(species); + } catch (const std::exception& e) { + throw std::runtime_error("Error converting A=" + std::to_string(A) + ", Z=" + std::to_string(Z) + " to Species: " + e.what()); + } + } + return isotopes; + } + + std::expected WeakRateInterpolator::get_rates( + const uint16_t A, + const uint8_t Z, + const double t9, + const double log_rhoYe, + const double mu_e + ) const { + const auto it = m_rate_table.find(pack_isotope_id(A, Z)); + if (it == m_rate_table.end()) { + return std::unexpected(InterpolationError{InterpolationErrorType::UNKNOWN_SPECIES_ERROR}); + } + const auto&[t9_axis, rhoYe_axis, mue_axis, data] = it->second; + + // Now find the bracketing indices for t9, log_rhoYe, and mu_e + auto find_lower_index = [](const std::vector& axis, const double value) -> std::optional { + const auto upperBoundIterator = std::ranges::upper_bound(axis, value); + if (upperBoundIterator == axis.begin() || upperBoundIterator == axis.end()) { + return std::nullopt; // Out of bounds + } + return std::distance(axis.begin(), upperBoundIterator) - 1; + }; + + const auto i_t9_opt = find_lower_index(t9_axis, t9); + const auto j_rhoYe_opt = find_lower_index(rhoYe_axis, log_rhoYe); + const auto k_mue_opt = find_lower_index(mue_axis, mu_e); + + if (!i_t9_opt || !j_rhoYe_opt || !k_mue_opt) { + std::unordered_map boundsInfo; + if (!i_t9_opt) { + boundsInfo[TableAxes::T9] = BoundsErrorInfo{ + TableAxes::T9, + t9_axis.front(), + t9_axis.back(), + t9 + }; + } + if (!j_rhoYe_opt) { + boundsInfo[TableAxes::LOG_RHOYE] = BoundsErrorInfo{ + TableAxes::LOG_RHOYE, + rhoYe_axis.front(), + rhoYe_axis.back(), + log_rhoYe + }; + } + if (!k_mue_opt) { + boundsInfo[TableAxes::MUE] = BoundsErrorInfo{ + TableAxes::MUE, + mue_axis.front(), + mue_axis.back(), + mu_e + }; + } + return std::unexpected( + InterpolationError{ + InterpolationErrorType::BOUNDS_ERROR, + boundsInfo + } + ); + } + + const size_t i = i_t9_opt.value(); + const size_t j = j_rhoYe_opt.value(); + const size_t k = k_mue_opt.value(); + + // Coordinates of the bounding cube + const double t1 = t9_axis[i]; + const double t2 = t9_axis[i + 1]; + const double r1 = rhoYe_axis[j]; + const double r2 = rhoYe_axis[j + 1]; + const double m1 = mue_axis[k]; + const double m2 = mue_axis[k + 1]; + + const double td = (t9 - t1) / (t2 - t1); + const double rd = (log_rhoYe - r1) / (r2 - r1); + const double md = (mu_e - m1) / (m2 - m1); + + auto lerp = [](const double v0, const double v1, const double t) { + return v0 * (1 - t) + v1 * t; + }; + + auto interpolationField = [&](auto field_accessor) { + const size_t nrhoYe = rhoYe_axis.size(); + const size_t nmue = mue_axis.size(); + + auto get_val = [&](const size_t i_t, const size_t j_r, const size_t k_m) { + return field_accessor(data[(i_t * nrhoYe + j_r) * nmue + k_m]); + }; + + const double c000 = get_val(i, j, k); + const double c001 = get_val(i, j, k + 1); + const double c010 = get_val(i, j + 1, k); + const double c011 = get_val(i, j + 1, k + 1); + const double c100 = get_val(i + 1, j, k); + const double c101 = get_val(i + 1, j, k + 1); + const double c110 = get_val(i + 1, j + 1, k); + const double c111 = get_val(i + 1, j + 1, k + 1); + + const double c00 = lerp(c000, c001, md); + const double c01 = lerp(c010, c011, md); + const double c10 = lerp(c100, c101, md); + const double c11 = lerp(c110, c111, md); + + const double c0 = lerp(c00, c01, rd); + const double c1 = lerp(c10, c11, rd); + + return lerp(c0, c1, td); + + }; + + WeakRatePayload result; + + result.log_beta_plus = interpolationField([](const WeakRatePayload& p) { return p.log_beta_plus; }); + result.log_electron_capture = interpolationField([](const WeakRatePayload& p) { return p.log_electron_capture; }); + result.log_neutrino_loss_ec = interpolationField([](const WeakRatePayload& p) { return p.log_neutrino_loss_ec; }); + result.log_beta_minus = interpolationField([](const WeakRatePayload& p) { return p.log_beta_minus; }); + result.log_positron_capture = interpolationField([](const WeakRatePayload& p) { return p.log_positron_capture; }); + result.log_antineutrino_loss_bd = interpolationField([](const WeakRatePayload& p) { return p.log_antineutrino_loss_bd; }); + return result; + } + + std::expected WeakRateInterpolator::get_rate_derivatives( + uint16_t A, + uint8_t Z, + double t9, + double log_rhoYe, + double mu_e + ) const { + WeakRateDerivatives result; + constexpr double eps = 1e-6; // Small perturbation for finite difference + + // Perturbations for finite difference + const double h_t9 = (t9 > 1e-9) ? t9 * eps : eps; + const auto payload_plus_t9 = get_rates(A, Z, t9 + h_t9, log_rhoYe, mu_e); + const auto payload_minus_t9 = get_rates(A, Z, t9 - h_t9, log_rhoYe, mu_e); + + const double h_rhoYe = (std::abs(log_rhoYe) > 1e-9) ? std::abs(log_rhoYe) * eps : eps; + const auto payload_plus_rhoYe = get_rates(A, Z, t9, log_rhoYe + h_rhoYe, mu_e); + const auto payload_minus_rhoYe = get_rates(A, Z, t9, log_rhoYe - h_rhoYe, mu_e); + + const double h_mue = (std::abs(mu_e) > 1e-9) ? std::abs(mu_e) * eps : eps; + const auto payload_plus_mue = get_rates(A, Z, t9, log_rhoYe, mu_e + h_mue); + const auto payload_minus_mue = get_rates(A, Z, t9, log_rhoYe, mu_e - h_mue); + + if (!payload_plus_t9 || !payload_minus_t9 || !payload_plus_rhoYe || !payload_minus_rhoYe || !payload_plus_mue || !payload_minus_mue) { + const auto it = m_rate_table.find(pack_isotope_id(A, Z)); + if (it == m_rate_table.end()) { + return std::unexpected(InterpolationError{InterpolationErrorType::UNKNOWN_SPECIES_ERROR}); + } + + const IsotopeGrid& grid = it->second; + InterpolationError error; + std::unordered_map boundsInfo; + if (!payload_minus_t9 || !payload_plus_t9) { + boundsInfo[TableAxes::T9] = BoundsErrorInfo{ + TableAxes::T9, + grid.t9_axis.front(), + grid.t9_axis.back(), + t9 + }; + } + if (!payload_minus_rhoYe || !payload_plus_rhoYe) { + boundsInfo[TableAxes::LOG_RHOYE] = BoundsErrorInfo{ + TableAxes::LOG_RHOYE, + grid.rhoYe_axis.front(), + grid.rhoYe_axis.back(), + log_rhoYe + }; + } + if (!payload_minus_mue || !payload_plus_mue) { + boundsInfo[TableAxes::MUE] = BoundsErrorInfo{ + TableAxes::MUE, + grid.mue_axis.front(), + grid.mue_axis.back(), + mu_e + }; + } + error.type = InterpolationErrorType::BOUNDS_ERROR; + error.boundsErrorInfo = boundsInfo; + return std::unexpected(error); + } + + // Derivatives wrt. T9 + const double t9_denominator = 2 * h_t9; + result.d_log_beta_plus[0] = (payload_plus_t9->log_beta_plus - payload_minus_t9->log_beta_plus) / t9_denominator; + result.d_log_beta_minus[0] = (payload_plus_t9->log_beta_minus - payload_minus_t9->log_beta_minus) / t9_denominator; + result.d_log_electron_capture[0] = (payload_plus_t9->log_electron_capture - payload_minus_t9->log_electron_capture) / t9_denominator; + result.d_log_neutrino_loss_ec[0] = (payload_plus_t9->log_neutrino_loss_ec - payload_minus_t9->log_neutrino_loss_ec) / t9_denominator; + result.d_log_positron_capture[0] = (payload_plus_t9->log_positron_capture - payload_minus_t9->log_positron_capture) / t9_denominator; + result.d_log_antineutrino_loss_bd[0] = (payload_plus_t9->log_antineutrino_loss_bd - payload_minus_t9->log_antineutrino_loss_bd) / t9_denominator; + + // Derivatives wrt. logRhoYe + const double rhoYe_denominator = 2 * h_rhoYe; + result.d_log_beta_plus[1] = (payload_plus_rhoYe->log_beta_plus - payload_minus_rhoYe->log_beta_plus) / rhoYe_denominator; + result.d_log_beta_minus[1] = (payload_plus_rhoYe->log_beta_minus - payload_minus_rhoYe->log_beta_minus) / rhoYe_denominator; + result.d_log_electron_capture[1] = (payload_plus_rhoYe->log_electron_capture - payload_minus_rhoYe->log_electron_capture) / rhoYe_denominator; + result.d_log_neutrino_loss_ec[1] = (payload_plus_rhoYe->log_neutrino_loss_ec - payload_minus_rhoYe->log_neutrino_loss_ec) / rhoYe_denominator; + result.d_log_positron_capture[1] = (payload_plus_rhoYe->log_positron_capture - payload_minus_rhoYe->log_positron_capture) / rhoYe_denominator; + result.d_log_antineutrino_loss_bd[1] = (payload_plus_rhoYe->log_antineutrino_loss_bd - payload_minus_rhoYe->log_antineutrino_loss_bd) / rhoYe_denominator; + + // Derivatives wrt. MuE + const double mue_denominator = 2 * h_mue; + result.d_log_beta_plus[2] = (payload_plus_mue->log_beta_plus - payload_minus_mue->log_beta_plus) / mue_denominator; + result.d_log_beta_minus[2] = (payload_plus_mue->log_beta_minus - payload_minus_mue->log_beta_minus) / mue_denominator; + result.d_log_electron_capture[2] = (payload_plus_mue->log_electron_capture - payload_minus_mue->log_electron_capture) / mue_denominator; + result.d_log_neutrino_loss_ec[2] = (payload_plus_mue->log_neutrino_loss_ec - payload_minus_mue->log_neutrino_loss_ec) / mue_denominator; + result.d_log_positron_capture[2] = (payload_plus_mue->log_positron_capture - payload_minus_mue->log_positron_capture) / mue_denominator; + result.d_log_antineutrino_loss_bd[2] = (payload_plus_mue->log_antineutrino_loss_bd - payload_minus_mue->log_antineutrino_loss_bd) / mue_denominator; + + return result; + } + + + uint32_t WeakRateInterpolator::pack_isotope_id(const uint16_t A, const uint8_t Z) { + return (static_cast(A) << 8) | static_cast(Z); + } + +} diff --git a/src/lib/screening/screening_bare.cpp b/src/lib/screening/screening_bare.cpp index 5d2d7e8b..838c93a1 100644 --- a/src/lib/screening/screening_bare.cpp +++ b/src/lib/screening/screening_bare.cpp @@ -12,7 +12,7 @@ namespace gridfire::screening { std::vector BareScreeningModel::calculateScreeningFactors( const reaction::ReactionSet &reactions, const std::vector& species, - const std::vector &Y, + const std::vector &Y, const ADDouble T9, const ADDouble rho ) const { diff --git a/src/lib/screening/screening_weak.cpp b/src/lib/screening/screening_weak.cpp index d3370b20..a866dcdd 100644 --- a/src/lib/screening/screening_weak.cpp +++ b/src/lib/screening/screening_weak.cpp @@ -12,7 +12,7 @@ namespace gridfire::screening { std::vector WeakScreeningModel::calculateScreeningFactors( const reaction::ReactionSet &reactions, const std::vector& species, - const std::vector &Y, + const std::vector &Y, const ADDouble T9, const ADDouble rho ) const { diff --git a/src/lib/solver/solver.cpp b/src/lib/solver/solver.cpp deleted file mode 100644 index 72d64c5f..00000000 --- a/src/lib/solver/solver.cpp +++ /dev/null @@ -1,356 +0,0 @@ -#include "gridfire/solver/solver.h" -#include "gridfire/engine/engine_graph.h" -#include "gridfire/network.h" -#include "gridfire/exceptions/error_engine.h" - -#include "fourdst/composition/atomicSpecies.h" -#include "fourdst/composition/composition.h" -#include "fourdst/config/config.h" - -#include "fourdst/plugin/plugin.h" - -#include "gridfire/interfaces/solver/solver_interfaces.h" - -#include "unsupported/Eigen/NonLinearOptimization" - -#include - -#include -#include -#include -#include - -#include "quill/LogMacros.h" - -namespace gridfire::solver { - NetOut DirectNetworkSolver::evaluate(const NetIn &netIn) { - namespace ublas = boost::numeric::ublas; - namespace odeint = boost::numeric::odeint; - using fourdst::composition::Composition; - - - const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9) - - const auto absTol = m_config.get("gridfire:solver:DirectNetworkSolver:absTol", 1.0e-8); - const auto relTol = m_config.get("gridfire:solver:DirectNetworkSolver:relTol", 1.0e-8); - - Composition equilibratedComposition = m_engine.update(netIn); - size_t numSpecies = m_engine.getNetworkSpecies().size(); - ublas::vector Y(numSpecies + 1); - - RHSManager manager(m_engine, T9, netIn.density, m_callback, m_engine.getNetworkSpecies()); - JacobianFunctor jacobianFunctor(m_engine, T9, netIn.density); - - auto populateY = [&](const Composition& comp) { - const size_t numSpeciesInternal = m_engine.getNetworkSpecies().size(); - Y.resize(numSpeciesInternal + 1); - for (size_t i = 0; i < numSpeciesInternal; i++) { - const auto& species = m_engine.getNetworkSpecies()[i]; - if (!comp.contains(species)) { - double lim = std::numeric_limits::min(); - LOG_DEBUG(m_logger, "Species '{}' not found in composition. Setting abundance to {:0.3E}.", species.name(), lim); - Y(i) = lim; // Species not in the composition, set to zero - } else { - Y(i) = comp.getMolarAbundance(species); - } - } - // TODO: a good starting point to make the temperature, density, and energy self consistent would be to turn this into an accumulator - Y(numSpeciesInternal) = 0.0; // Specific energy rate, initialized to zero - }; - - // This is a quick debug that can be turned on. For solar code input parameters (T~1.5e7K, ρ~1.5e3 g/cm^3) this should be near 8e-17 - // std::cout << "D/H: " << equilibratedComposition.getMolarAbundance("H-2") / equilibratedComposition.getMolarAbundance("H-1") << std::endl; - - populateY(equilibratedComposition); - const auto stepper = odeint::make_controlled>(absTol, relTol); - - double current_time = 0.0; - double current_initial_timestep = netIn.dt0; - double accumulated_energy = 0.0; - // size_t total_update_stages_triggered = 0; - - while (current_time < netIn.tMax) { - try { - odeint::integrate_adaptive( - stepper, - std::make_pair(manager, jacobianFunctor), - Y, - current_time, - netIn.tMax, - current_initial_timestep, - [&](const auto& state, double t) { - current_time = t; - manager.observe(state, t); - } - ); - current_time = netIn.tMax; - } catch (const exceptions::StaleEngineTrigger &e) { - LOG_INFO(m_logger, "Catching StaleEngineTrigger at t = {:0.3E} with T9 = {:0.3E}, rho = {:0.3E}. Triggering update stage (last stage took {} steps).", current_time, T9, netIn.density, e.totalSteps()); - exceptions::StaleEngineTrigger::state staleState = e.getState(); - accumulated_energy += e.energy(); // Add the specific energy rate to the accumulated energy - // total_update_stages_triggered++; - - 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) { - const auto& species = m_engine.getNetworkSpecies()[i]; - temp_comp.registerSpecies(species); - mass_fractions.push_back(e.getMolarAbundance(i) * species.mass()); // Convert from molar abundance to mass fraction - } - temp_comp.setMassFraction(m_engine.getNetworkSpecies(), mass_fractions); - temp_comp.finalize(true); - - NetIn netInTemp = netIn; - netInTemp.temperature = e.temperature(); - netInTemp.density = e.density(); - netInTemp.composition = temp_comp; - - Composition currentComposition = m_engine.update(netInTemp); - populateY(currentComposition); - Y(Y.size() - 1) = e.energy(); // Set the specific energy rate from the stale state - numSpecies = m_engine.getNetworkSpecies().size(); - - // current_initial_timestep = 0.001 * manager.m_last_step_time; // set the new timestep to the last successful timestep before repartitioning - } - } - - accumulated_energy += Y(Y.size() - 1); // Add the specific energy rate to the accumulated energy - - std::vector finalMassFractions(numSpecies); - for (size_t i = 0; i < numSpecies; ++i) { - const double molarMass = m_engine.getNetworkSpecies()[i].mass(); - finalMassFractions[i] = Y(i) * molarMass; // Convert from molar abundance to mass fraction - if (finalMassFractions[i] < MIN_ABUNDANCE_THRESHOLD) { - finalMassFractions[i] = 0.0; - } - } - - std::vector speciesNames; - speciesNames.reserve(numSpecies); - for (const auto& species : m_engine.getNetworkSpecies()) { - speciesNames.emplace_back(species.name()); - } - - Composition outputComposition(speciesNames); - outputComposition.setMassFraction(speciesNames, finalMassFractions); - outputComposition.finalize(true); - - NetOut netOut; - 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; - } - - void DirectNetworkSolver::set_callback(const std::any& callback) { - if (!callback.has_value()) { - m_callback = {}; - return; - } - - using FunctionPtrType = void (*)(const TimestepContext&); - - if (callback.type() == typeid(TimestepCallback)) { - m_callback = std::any_cast(callback); - } - else if (callback.type() == typeid(FunctionPtrType)) { - auto func_ptr = std::any_cast(callback); - m_callback = func_ptr; - } - else { - throw std::invalid_argument("Unsupported type passed to set_callback. " - "Provide a std::function or a matching function pointer."); - } - } - std::vector> DirectNetworkSolver::describe_callback_context() const { - const TimestepContext context( - 0.0, // time - boost::numeric::ublas::vector(), // state - 0.0, // dt - 0.0, // cached_time - 0.0, // last_observed_time - 0.0, // last_step_time - 0.0, // T9 - 0.0, // rho - std::nullopt, // cached_result - 0, // num_steps - m_engine, // engine, - {} - ); - return context.describe(); - } - - void DirectNetworkSolver::RHSManager::operator()( - const boost::numeric::ublas::vector &Y, - boost::numeric::ublas::vector &dYdt, - const double t - ) const { - const size_t numSpecies = m_engine.getNetworkSpecies().size(); - if (t != m_cached_time || !m_cached_result.has_value() || m_cached_result.value().dydt.size() != numSpecies + 1) { - compute_and_cache(Y, t); - } - const auto&[dydt, nuclearEnergyGenerationRate] = m_cached_result.value(); - dYdt.resize(numSpecies + 1); - std::ranges::copy(dydt, dYdt.begin()); - dYdt(numSpecies) = nuclearEnergyGenerationRate; // Set the last element to the specific energy rate - } - - void DirectNetworkSolver::RHSManager::observe( - const boost::numeric::ublas::vector &state, - const double t - ) const { - double dt = t - m_last_observed_time; - compute_and_cache(state, t); - LOG_INFO( - m_logger, - "(Step {}) Observed state at t = {:0.3E} (dt = {:0.3E})", - m_num_steps, - t, - dt - ); - std::ostringstream oss; - oss << std::scientific << std::setprecision(3); - oss << "(Step: " << std::setw(10) << m_num_steps << ") t = " << t << " (dt = " << dt << ", eps_nuc: " << state(state.size() - 1) << " [erg])\n"; - std::cout << oss.str(); - - fourdst::plugin::manager::PluginManager &pluginManager = fourdst::plugin::manager::PluginManager::getInstance(); - - if (pluginManager.has("gridfire/solver")) { - auto* plugin = pluginManager.get("gridfire/solver"); - plugin -> log_time(t, dt); - } - - // Callback logic - if (m_callback) { - LOG_TRACE_L1(m_logger, "Calling user callback function at t = {:0.3E} with dt = {:0.3E}", t, dt); - const TimestepContext context( - t, - state, - dt, - m_cached_time, - m_last_observed_time, - m_last_step_time, - m_T9, - m_rho, - m_cached_result, - m_num_steps, - m_engine, - m_networkSpecies - ); - - m_callback(context); - LOG_TRACE_L1(m_logger, "User callback function completed at t = {:0.3E} with dt = {:0.3E}", t, dt); - } - - m_last_observed_time = t; - m_last_step_time = dt; - - } - - void DirectNetworkSolver::RHSManager::compute_and_cache( - const boost::numeric::ublas::vector &state, - double t - ) const { - std::vector y_vec(state.begin(), state.end() - 1); - std::ranges::replace_if( - y_vec, - [](const double yi){ - return yi < 0.0; - }, - 0.0 // Avoid negative abundances - ); - - const auto result = m_engine.calculateRHSAndEnergy(y_vec, m_T9, m_rho); - if (!result) { - LOG_INFO(m_logger, - "Triggering update stage due to stale engine detected at t = {:0.3E} with T9 = {:0.3E}, rho = {:0.3E}, y_vec (size: {})", - t, m_T9, m_rho, y_vec.size()); - throw exceptions::StaleEngineTrigger({m_T9, m_rho, y_vec, t, m_num_steps, state(state.size() - 1)}); - } - m_cached_result = result.value(); - m_cached_time = t; - - m_num_steps++; - } - - void DirectNetworkSolver::JacobianFunctor::operator()( - const boost::numeric::ublas::vector &Y, - boost::numeric::ublas::matrix &J, - double t, - boost::numeric::ublas::vector &dfdt - ) const { - size_t numSpecies = m_engine.getNetworkSpecies().size(); - J.resize(numSpecies+1, numSpecies+1); - J.clear(); - for (size_t i = 0; i < numSpecies; ++i) { - for (size_t j = 0; j < numSpecies; ++j) { - J(i, j) = m_engine.getJacobianMatrixEntry(i, j); - } - } - } - - DirectNetworkSolver::TimestepContext::TimestepContext( - const double t, - const boost::numeric::ublas::vector &state, - const double dt, - const double cached_time, - const double last_observed_time, - const double last_step_time, - const double t9, - const double rho, - const std::optional> &cached_result, - const int num_steps, - const DynamicEngine &engine, - const std::vector &networkSpecies - ) - : t(t), - state(state), - dt(dt), - cached_time(cached_time), - last_observed_time(last_observed_time), - last_step_time(last_step_time), - T9(t9), - rho(rho), - cached_result(cached_result), - num_steps(num_steps), - engine(engine), - networkSpecies(networkSpecies) {} - - std::vector> DirectNetworkSolver::TimestepContext::describe() const { - return { - {"time", "double"}, - {"state", "boost::numeric::ublas::vector&"}, - {"dt", "double"}, - {"cached_time", "double"}, - {"last_observed_time", "double"}, - {"last_step_time", "double"}, - {"T9", "double"}, - {"rho", "double"}, - {"cached_result", "std::optional>&"}, - {"num_steps", "int"}, - {"engine", "DynamicEngine&"}, - {"networkSpecies", "std::vector&"} - }; - } - - -} \ No newline at end of file diff --git a/src/lib/solver/strategies/CVODE_solver_strategy.cpp b/src/lib/solver/strategies/CVODE_solver_strategy.cpp index 8777b2f4..6b3df074 100644 --- a/src/lib/solver/strategies/CVODE_solver_strategy.cpp +++ b/src/lib/solver/strategies/CVODE_solver_strategy.cpp @@ -17,6 +17,7 @@ #include #include "fourdst/composition/exceptions/exceptions_composition.h" +#include "gridfire/engine/engine_graph.h" #include "gridfire/solver/strategies/triggers/engine_partitioning_trigger.h" #include "gridfire/trigger/procedures/trigger_pprint.h" @@ -71,7 +72,7 @@ namespace { #elif SUNDIALS_HAVE_PTHREADS N_Vector vec = N_VNew_Pthreads(size, sun_ctx); #else - N_Vector vec = N_VNew_Serial(size, sun_ctx); + N_Vector vec = N_VNew_Serial(static_cast(size), sun_ctx); #endif check_cvode_flag(vec == nullptr ? -1 : 0, "N_VNew"); return vec; @@ -87,7 +88,7 @@ namespace gridfire::solver { const double last_step_time, const double t9, const double rho, - const int num_steps, + const size_t num_steps, const DynamicEngine &engine, const std::vector &networkSpecies ) : @@ -157,6 +158,7 @@ namespace gridfire::solver { user_data.engine = &m_engine; double current_time = 0; + // ReSharper disable once CppTooWideScope [[maybe_unused]] double last_callback_time = 0; m_num_steps = 0; double accumulated_energy = 0.0; @@ -169,7 +171,7 @@ namespace gridfire::solver { check_cvode_flag(CVodeSetUserData(m_cvode_mem, &user_data), "CVodeSetUserData"); - int flag = -1; + int flag{}; if (m_stdout_logging_enabled) { flag = CVode(m_cvode_mem, netIn.tMax, m_Y, ¤t_time, CV_ONE_STEP); } else { @@ -304,10 +306,8 @@ namespace gridfire::solver { 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 + outputComposition, T9, netIn.density ); @@ -374,9 +374,11 @@ namespace gridfire::solver { for (size_t j = 0; j < numSpecies; ++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(i, j); + J_data[j * N + i] = engine->getJacobianMatrixEntry(species_i, species_j); } } @@ -398,11 +400,30 @@ namespace gridfire::solver { const size_t numSpecies = m_engine.getNetworkSpecies().size(); sunrealtype* y_data = N_VGetArrayPointer(y); + // 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 + // the index in the vector being consistent. The trade off is that sometimes we need to construct a composition object + // which, at the moment, requires a somewhat expensive set of operations. Perhaps in the future we could enforce + // some consistent memory layout for the composition object to make this cheeper. That optimization would need to be + // done in the libcomposition library though... std::vector y_vec(y_data, y_data + numSpecies); + std::vector symbols; + symbols.reserve(numSpecies); + for (const auto& species : m_engine.getNetworkSpecies()) { + symbols.emplace_back(species.name()); + } + std::vector X; + X.reserve(numSpecies); + for (size_t i = 0; i < numSpecies; ++i) { + const double molarMass = m_engine.getNetworkSpecies()[i].mass(); + X.push_back(y_vec[i] * molarMass); // Convert from molar abundance to mass fraction + } + fourdst::composition::Composition composition(symbols, X); std::ranges::replace_if(y_vec, [](const double val) { return val < 0.0; }, 0.0); - const auto result = m_engine.calculateRHSAndEnergy(y_vec, data->T9, data->rho); + const auto result = m_engine.calculateRHSAndEnergy(composition, data->T9, data->rho); if (!result) { throw exceptions::StaleEngineTrigger({data->T9, data->rho, y_vec, t, m_num_steps, y_data[numSpecies]}); } @@ -411,7 +432,7 @@ namespace gridfire::solver { const auto& [dydt, nuclearEnergyGenerationRate] = result.value(); for (size_t i = 0; i < numSpecies; ++i) { - ydot_data[i] = dydt[i]; + ydot_data[i] = dydt.at(m_engine.getNetworkSpecies()[i]); } ydot_data[numSpecies] = nuclearEnergyGenerationRate; // Set the last element to the specific energy rate } @@ -513,6 +534,7 @@ namespace gridfire::solver { std::vector Y_full(y_data, y_data + num_components - 1); + std::ranges::replace_if( Y_full, [](const double val) { @@ -528,8 +550,9 @@ namespace gridfire::solver { 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())); + speciesNames.emplace_back(user_data.networkSpecies->at(i).name()); } + fourdst::composition::Composition composition(speciesNames, Y_full); if (err_ratios.empty()) { return; @@ -565,9 +588,9 @@ namespace gridfire::solver { 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); + diagnostics::inspect_jacobian_stiffness(*user_data.engine, composition, user_data.T9, user_data.rho); + diagnostics::inspect_species_balance(*user_data.engine, "N-14", composition, user_data.T9, user_data.rho); + diagnostics::inspect_species_balance(*user_data.engine, "n-1", composition, user_data.T9, user_data.rho); } } diff --git a/src/lib/solver/strategies/triggers/engine_partitioning_trigger.cpp b/src/lib/solver/strategies/triggers/engine_partitioning_trigger.cpp index 6d8789f4..0cb30a3b 100644 --- a/src/lib/solver/strategies/triggers/engine_partitioning_trigger.cpp +++ b/src/lib/solver/strategies/triggers/engine_partitioning_trigger.cpp @@ -94,13 +94,12 @@ namespace gridfire::trigger::solver::CVODE { } bool OffDiagonalTrigger::check(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) const { - const size_t numSpecies = ctx.engine.getNetworkSpecies().size(); - for (int row = 0; row < numSpecies; ++row) { - for (int col = 0; col < numSpecies; ++col) { - double DRowDCol = std::abs(ctx.engine.getJacobianMatrixEntry(row, col)); - if (row != col && DRowDCol > m_threshold) { + for (const auto& rowSpecies : ctx.engine.getNetworkSpecies()) { + for (const auto& colSpecies : ctx.engine.getNetworkSpecies()) { + double DRowDCol = std::abs(ctx.engine.getJacobianMatrixEntry(rowSpecies, colSpecies)); + if (rowSpecies != colSpecies && DRowDCol > m_threshold) { m_hits++; - LOG_TRACE_L2(m_logger, "OffDiagonalTrigger triggered at t = {} due to entry ({}, {}) = {}", ctx.t, row, col, DRowDCol); + LOG_TRACE_L2(m_logger, "OffDiagonalTrigger triggered at t = {} due to entry ({}, {}) = {}", ctx.t, rowSpecies.name(), colSpecies.name(), DRowDCol); return true; } } @@ -173,7 +172,7 @@ namespace gridfire::trigger::solver::CVODE { } bool TimestepCollapseTrigger::check(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) const { - if (m_timestep_window.size() < 1) { + if (m_timestep_window.empty()) { m_misses++; return false; } @@ -181,7 +180,7 @@ namespace gridfire::trigger::solver::CVODE { for (const auto& dt : m_timestep_window) { averageTimestep += dt; } - averageTimestep /= m_timestep_window.size(); + averageTimestep /= static_cast(m_timestep_window.size()); if (m_relative && (std::abs(ctx.dt - averageTimestep) / averageTimestep) >= m_threshold) { m_hits++; LOG_TRACE_L2(m_logger, "TimestepCollapseTrigger triggered at t = {} due to relative growth: dt = {}, average dt = {}, threshold = {}", ctx.t, ctx.dt, averageTimestep, m_threshold); @@ -250,6 +249,12 @@ namespace gridfire::trigger::solver::CVODE { using ctx_t = gridfire::solver::CVODESolverStrategy::TimestepContext; // Create the individual conditions that can trigger a repartitioning + // The current trigger logic is as follows + // 1. Trigger every 1000th time that the simulation time exceeds the simulationTimeInterval + // 2. OR if any off-diagonal Jacobian entry exceeds the offDiagonalThreshold + // 3. OR every 10th time that the timestep growth exceeds the timestepGrowthThreshold (relative or absolute) + + // TODO: This logic likely needs to be revisited; however, for now it is easy enough to change and test and it works reasonably well auto simulationTimeTrigger = std::make_unique>(std::make_unique(simulationTimeInterval), 1000); auto offDiagTrigger = std::make_unique(offDiagonalThreshold); auto timestepGrowthTrigger = std::make_unique>(std::make_unique(timestepGrowthThreshold, timestepGrowthRelative, timestepGrowthWindowSize), 10); diff --git a/src/lib/utils/logging.cpp b/src/lib/utils/logging.cpp index c05714b8..a764ef72 100644 --- a/src/lib/utils/logging.cpp +++ b/src/lib/utils/logging.cpp @@ -11,11 +11,11 @@ std::string gridfire::utils::formatNuclearTimescaleLogString( const DynamicEngine& engine, - std::vector const& Y, + const fourdst::composition::Composition& composition, const double T9, const double rho ) { - auto const& result = engine.getSpeciesTimescales(Y, T9, rho); + auto const& result = engine.getSpeciesTimescales(composition, T9, rho); if (!result) { std::ostringstream ss; ss << "Failed to get species timescales: " << result.error(); diff --git a/src/meson.build b/src/meson.build index c8ca70b3..99b1e43f 100644 --- a/src/meson.build +++ b/src/meson.build @@ -12,8 +12,9 @@ gridfire_sources = files( 'lib/engine/diagnostics/dynamic_engine_diagnostics.cpp', 'lib/reaction/reaction.cpp', 'lib/reaction/reaclib.cpp', + 'lib/reaction/weak/weak.cpp', + 'lib/reaction/weak/weak_interpolator.cpp', 'lib/io/network_file.cpp', - 'lib/solver/solver.cpp', 'lib/solver/strategies/CVODE_solver_strategy.cpp', 'lib/solver/strategies/triggers/engine_partitioning_trigger.cpp', 'lib/screening/screening_types.cpp', diff --git a/src/python/reaction/bindings.cpp b/src/python/reaction/bindings.cpp index 56c48a66..edecc3aa 100644 --- a/src/python/reaction/bindings.cpp +++ b/src/python/reaction/bindings.cpp @@ -48,7 +48,7 @@ void register_reaction_bindings(py::module &m) { .def( "calculate_rate", [](const gridfire::reaction::ReaclibReaction& self, const double T9, const double rho, const std::vector& Y) -> double { - return self.calculate_rate(T9, rho, Y); + return self.calculate_rate(T9, rho, 0, TODO, Y, TODO); }, py::arg("T9"), py::arg("rho"), @@ -211,7 +211,7 @@ void register_reaction_bindings(py::module &m) { .def( "calculate_rate", [](const gridfire::reaction::LogicalReaclibReaction& self, const double T9, const double rho, const std::vector& Y) -> double { - return self.calculate_rate(T9, rho, Y); + return self.calculate_rate(T9, rho, 0, TODO, Y, TODO); }, py::arg("T9"), "Calculate the reaction rate at a given temperature T9 (in units of 10^9 K)." diff --git a/subprojects/fourdst.wrap b/subprojects/fourdst.wrap index fc133ce0..d23a0c54 100644 --- a/subprojects/fourdst.wrap +++ b/subprojects/fourdst.wrap @@ -1,4 +1,4 @@ [wrap-git] url = https://github.com/4D-STAR/fourdst -revision = v0.8.0 +revision = v0.8.1 depth = 1 diff --git a/tests/graphnet_sandbox/main.cpp b/tests/graphnet_sandbox/main.cpp index a63cc543..513dd22a 100644 --- a/tests/graphnet_sandbox/main.cpp +++ b/tests/graphnet_sandbox/main.cpp @@ -5,12 +5,7 @@ #include "gridfire/engine/engine_approx8.h" #include "gridfire/engine/views/engine_adaptive.h" #include "gridfire/partition/partition_types.h" -#include "gridfire/engine/views/engine_defined.h" #include "gridfire/engine/views/engine_multiscale.h" -#include "gridfire/engine/procedures/priming.h" -#include "gridfire/io/network_file.h" - -#include "gridfire/solver/solver.h" #include "gridfire/solver/strategies/CVODE_solver_strategy.h" #include "gridfire/network.h" @@ -32,19 +27,6 @@ static std::terminate_handler g_previousHandler = nullptr; - - -void callback(const gridfire::solver::DirectNetworkSolver::TimestepContext& ctx) { - const auto H1IndexPtr = std::ranges::find(ctx.engine.getNetworkSpecies(), fourdst::atomic::H_1); - const auto He4IndexPtr = std::ranges::find(ctx.engine.getNetworkSpecies(), fourdst::atomic::He_4); - - const size_t H1Index = H1IndexPtr != ctx.engine.getNetworkSpecies().end() ? std::distance(ctx.engine.getNetworkSpecies().begin(), H1IndexPtr) : -1; - const size_t He4Index = He4IndexPtr != ctx.engine.getNetworkSpecies().end() ? std::distance(ctx.engine.getNetworkSpecies().begin(), He4IndexPtr) : -1; - - std::cout << "Time: " << ctx.t << ", H-1: " << ctx.state(H1Index) << ", He-4: " << ctx.state(He4Index) << "\n"; - -} - void measure_execution_time(const std::function& callback, const std::string& name) { const auto startTime = std::chrono::steady_clock::now();