From c3bc75a7f4c74253d8f33691cad87910b615df5b Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 24 Jul 2025 08:37:52 -0400 Subject: [PATCH] docs(GridFire): added loads of docs and supressed yaml-cpp shadow warnings --- .gitignore | 2 + src/include/gridfire/engine/engine.h | 174 ++++ src/include/gridfire/engine/engine_graph.h | 13 +- .../gridfire/engine/views/engine_multiscale.h | 878 ++++++++++++++++-- .../gridfire/screening/screening_weak.h | 2 +- src/lib/engine/engine_graph.cpp | 4 +- src/lib/engine/views/engine_multiscale.cpp | 47 +- src/lib/io/network_file.cpp | 4 +- src/python/engine/trampoline/py_engine.cpp | 46 +- src/python/engine/trampoline/py_engine.h | 7 + subprojects/fourdst.wrap | 4 +- tests/python/test.py | 2 +- 12 files changed, 1061 insertions(+), 122 deletions(-) diff --git a/.gitignore b/.gitignore index 75a5fe4d..55fafd44 100644 --- a/.gitignore +++ b/.gitignore @@ -57,6 +57,8 @@ tags *.private *.private/ +*.bin + subprojects/mfem/ subprojects/tetgen/ subprojects/yaml-cpp/ diff --git a/src/include/gridfire/engine/engine.h b/src/include/gridfire/engine/engine.h index cb74acca..028315a9 100644 --- a/src/include/gridfire/engine/engine.h +++ b/src/include/gridfire/engine/engine.h @@ -1,3 +1,177 @@ +/** + * @file engine.h + * @brief Core header for the GridFire reaction network engine module. + * + * This module defines the core interfaces and classes for reaction network + * engines in GridFire. It provides abstract base classes for engines, + * dynamic engines, and engine views, as well as concrete engine + * implementations and view implementations. + * + * The engine module is designed to support a wide range of reaction network + * simulations, from simple single-zone calculations to complex multi-zone + * simulations with adaptive network topologies. + * + * @section EngineDesign Engine Design + * + * The engine module is built around the following key concepts: + * + * - **Engine:** The base class for all reaction network engines. It defines + * the minimal interface for evaluating the right-hand side (dY/dt) and + * energy generation rate for a given set of abundances, temperature, and + * density. + * + * - **DynamicEngine:** An extension of the Engine class that supports + * Jacobian and stoichiometry operations, as well as the ability to + * dynamically modify the reaction network. + * + * - **EngineView:** An abstract base class for "views" of reaction network + * engines. Engine views provide a way to dynamically or adaptively + * modify the network topology without modifying the underlying physics + * engine. + * + * @section EngineComposition Engine Composition + * + * Engines and engine views can be composed to create complex reaction network + * simulations. For example, an AdaptiveEngineView can be used to dynamically + * cull species and reactions from a GraphEngine, reducing the computational + * cost of the simulation. + * + * The order in which engines and engine views are composed is important. The + * base engine should always be the innermost engine, and the engine views + * should be layered on top of the base engine. + * + * @section AvailableEngines Available Engines + * + * The engine module provides the following concrete engine implementations: + * + * - **GraphEngine:** A reaction network engine that uses a graph-based + * representation of the reaction network. It uses sparse matrices for + * efficient storage and computation of the stoichiometry and Jacobian + * matrices. + * + * @section AvailableViews Available Views + * + * The engine module provides the following engine view implementations: + * + * - **AdaptiveEngineView:** An engine view that dynamically adapts the + * reaction network based on runtime conditions. It culls species and + * reactions with low reaction flow rates, reducing the computational + * cost of the simulation. + * + * - **DefinedEngineView:** An engine view that restricts the reaction + * network to a predefined set of species and reactions. This can be + * useful for simulating specific reaction pathways or for comparing + * results with other codes. + * + * - **MultiscalePartitioningEngineView:** An engine view that partitions the + * reaction network into multiple groups based on timescales. This can be + * useful for simulating stiff reaction networks, where some reactions + * occur much faster than others. + * + * - **NetworkPrimingEngineView:** An engine view that primes the reaction + * network with a specific species or set of species. This can be useful + * for igniting a reaction network or for studying the effects of specific + * species on the network. + * + * @section UsageExamples Usage Examples + * + * @subsection GraphEngineExample GraphEngine Example + * + * The following code shows how to create a GraphEngine from a composition: + * + * @code + * #include "gridfire/engine/engine_graph.h" + * #include "fourdst/composition/composition.h" + * + * // Create a composition + * fourdst::composition::Composition composition; + * + * // Create a GraphEngine + * gridfire::GraphEngine engine(composition); + * @endcode + * + * @subsection AdaptiveEngineViewExample AdaptiveEngineView Example + * + * The following code shows how to create an AdaptiveEngineView from a + * GraphEngine: + * + * @code + * #include "gridfire/engine/views/engine_adaptive.h" + * #include "gridfire/engine/engine_graph.h" + * #include "fourdst/composition/composition.h" + * + * // Create a composition + * fourdst::composition::Composition composition; + * + * // Create a GraphEngine + * gridfire::GraphEngine baseEngine(composition); + * + * // Create an AdaptiveEngineView + * gridfire::AdaptiveEngineView engine(baseEngine); + * @endcode + * + * @subsection DefinedEngineViewExample DefinedEngineView Example + * + * The following code shows how to create a DefinedEngineView from a + * GraphEngine: + * + * @code + * #include "gridfire/engine/views/engine_defined.h" + * #include "gridfire/engine/engine_graph.h" + * #include "fourdst/composition/composition.h" + * + * // Create a composition + * fourdst::composition::Composition composition; + * + * // Create a GraphEngine + * gridfire::GraphEngine baseEngine(composition); + * + * // Create a DefinedEngineView + * std::vector peNames = {"p(p,e+)d", "he4(a,g)be8"}; + * gridfire::DefinedEngineView engine(peNames, baseEngine); + * @endcode + * + * @subsection MultiscalePartitioningEngineViewExample MultiscalePartitioningEngineView Example + * + * The following code shows how to create a MultiscalePartitioningEngineView from a + * GraphEngine: + * + * @code + * #include "gridfire/engine/views/engine_multiscale.h" + * #include "gridfire/engine/engine_graph.h" + * #include "fourdst/composition/composition.h" + * + * // Create a composition + * fourdst::composition::Composition composition; + * + * // Create a GraphEngine + * gridfire::GraphEngine baseEngine(composition); + * + * // Create a MultiscalePartitioningEngineView + * gridfire::MultiscalePartitioningEngineView engine(baseEngine); + * @endcode + * + * @subsection NetworkPrimingEngineViewExample NetworkPrimingEngineView Example + * + * The following code shows how to create a NetworkPrimingEngineView from a + * GraphEngine: + * + * @code + * #include "gridfire/engine/views/engine_priming.h" + * #include "gridfire/engine/engine_graph.h" + * #include "fourdst/composition/composition.h" + * + * // Create a composition + * fourdst::composition::Composition composition; + * + * // Create a GraphEngine + * gridfire::GraphEngine baseEngine(composition); + * + * // Create a NetworkPrimingEngineView + * std::string primingSymbol = "p"; + * gridfire::NetworkPrimingEngineView engine(primingSymbol, baseEngine); + * @endcode + */ #pragma once #include "gridfire/engine/engine_abstract.h" diff --git a/src/include/gridfire/engine/engine_graph.h b/src/include/gridfire/engine/engine_graph.h index f48bc310..8f77bb6a 100644 --- a/src/include/gridfire/engine/engine_graph.h +++ b/src/include/gridfire/engine/engine_graph.h @@ -833,14 +833,6 @@ namespace gridfire { // --- Pre-setup (flags to control conditionals in an AD safe / branch aware manner) --- // ----- Constants for AD safe calculations --- const T zero = static_cast(0.0); - const T one = static_cast(1.0); - - // ----- Initialize variables for molar concentration product and thresholds --- - // Note: the logic here is that we use CppAD::CondExprLt to test thresholds and if they are less we set the flag - // to zero so that the final returned reaction flow is 0. This is as opposed to standard if statements - // which create branches that break the AD tape. - const T Y_threshold = static_cast(MIN_ABUNDANCE_THRESHOLD); - T threshold_flag = one; // --- 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); @@ -864,9 +856,6 @@ namespace gridfire { const size_t species_index = species_it->second; const T Yi = Y[species_index]; - // --- Check if the species abundance is below the threshold where we ignore reactions --- - // threshold_flag *= CppAD::CondExpLt(Yi, Y_threshold, zero, one); - // --- If count is > 1 , we need to raise the molar concentration to the power of count since there are really count bodies in that reaction --- molar_concentration_product *= CppAD::pow(Yi, static_cast(count)); // ni^count @@ -881,7 +870,7 @@ namespace gridfire { // the tape more expensive to record, but it will also mean that we only need to record it once for // the entire network. const T densityTerm = CppAD::pow(rho, totalReactants > 1 ? static_cast(totalReactants - 1) : zero); // Density raised to the power of (N-1) for N reactants - return molar_concentration_product * k_reaction * threshold_flag * densityTerm; + return molar_concentration_product * k_reaction * densityTerm; } }; \ No newline at end of file diff --git a/src/include/gridfire/engine/views/engine_multiscale.h b/src/include/gridfire/engine/views/engine_multiscale.h index 631029ea..afb5e83d 100644 --- a/src/include/gridfire/engine/views/engine_multiscale.h +++ b/src/include/gridfire/engine/views/engine_multiscale.h @@ -7,12 +7,43 @@ #include "unsupported/Eigen/NonLinearOptimization" namespace gridfire { + /** + * @brief Configuration struct for the QSE cache. + * + * @purpose This struct defines the tolerances used to determine if a QSE cache key + * is considered a hit. It allows for tuning the sensitivity of the cache. + * + * @how It works by providing binning widths for temperature, density, and abundances. + * When a `QSECacheKey` is created, it uses these tolerances to discretize the + * continuous physical values into bins. If two sets of conditions fall into the + * same bins, they will produce the same hash and be considered a cache hit. + * + * @par Usage Example: + * Although not typically set by the user directly, the `QSECacheKey` uses this + * internally. A smaller tolerance (e.g., `T9_tol = 1e-4`) makes the cache more + * sensitive, leading to more frequent re-partitions, while a larger tolerance + * (`T9_tol = 1e-2`) makes it less sensitive. + */ struct QSECacheConfig { double T9_tol; ///< Absolute tolerance to produce the same hash for T9. double rho_tol; ///< Absolute tolerance to produce the same hash for rho. double Yi_tol; ///< Absolute tolerance to produce the same hash for species abundances. }; + /** + * @brief Key struct for the QSE abundance cache. + * + * @purpose This struct is used as the key for the QSE abundance cache (`m_qse_abundance_cache`) + * within the `MultiscalePartitioningEngineView`. Its primary goal is to avoid + * expensive re-partitioning and QSE solves for thermodynamic conditions that are + * "close enough" to previously computed ones. + * + * @how It works by storing the temperature (`m_T9`), density (`m_rho`), and species + * abundances (`m_Y`). A pre-computed hash is generated in the constructor by + * calling the `hash()` method. This method discretizes the continuous physical + * values into bins using the tolerances defined in `QSECacheConfig`. The `operator==` + * simply compares the pre-computed hash values for fast lookups in the `std::unordered_map`. + */ struct QSECacheKey { double m_T9; double m_rho; @@ -27,16 +58,46 @@ namespace gridfire { 1e-3 // Default tolerance for species abundances }; + /** + * @brief Constructs a QSECacheKey. + * + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @param Y Species molar abundances. + * + * @post The `m_hash` member is computed and stored. + */ QSECacheKey( const double T9, const double rho, const std::vector& Y ); + /** + * @brief Computes the hash value for this key. + * + * @return The computed hash value. + * + * @how This method combines the hashes of the binned temperature, density, and + * each species abundance. The `bin()` static method is used for discretization. + */ size_t hash() const; + /** + * @brief Converts a value to a discrete bin based on a tolerance. + * @param value The value to bin. + * @param tol The tolerance (bin width) to use for binning. + * @return The bin number as a long integer. + * + * @how The algorithm is `floor(value / tol)`. + */ static long bin(double value, double tol); + /** + * @brief Equality operator for QSECacheKey. + * @param other The other QSECacheKey to compare to. + * @return True if the pre-computed hashes are equal, false otherwise. + */ bool operator==(const QSECacheKey& other) const; }; @@ -47,7 +108,7 @@ namespace std { template <> struct hash { /** - * @brief Computes the hash of a QSECacheKey. + * @brief Computes the hash of a QSECacheKey for use in `std::unordered_map`. * @param key The QSECacheKey to hash. * @return The pre-computed hash value of the key. */ @@ -59,37 +120,209 @@ namespace std { } // namespace std namespace gridfire { + /** + * @class MultiscalePartitioningEngineView + * @brief An engine view that partitions the reaction network into multiple groups based on timescales. + * + * @purpose This class is designed to accelerate the integration of stiff nuclear reaction networks. + * It identifies species that react on very short timescales ("fast" species) and treats them + * as being in Quasi-Steady-State Equilibrium (QSE). Their abundances are solved for algebraically, + * removing their stiff differential equations from the system. The remaining "slow" or "dynamic" + * species are integrated normally. This significantly improves the stability and performance of + * the solver. + * + * @how The core logic resides in the `partitionNetwork()` and `equilibrateNetwork()` methods. + * The partitioning process involves: + * 1. **Timescale Analysis:** Using `getSpeciesDestructionTimescales` from the base engine, + * all species are sorted by their characteristic timescales. + * 2. **Gap Detection:** The sorted list of timescales is scanned for large gaps (e.g., several + * orders of magnitude) to create distinct "timescale pools". + * 3. **Connectivity Analysis:** Each pool is analyzed for internal reaction connectivity to + * form cohesive groups. + * 4. **Flux Validation:** Candidate QSE groups are validated by comparing the total reaction + * flux *within* the group to the flux *leaving* the group. A high internal-to-external + * flux ratio indicates a valid QSE group. + * 5. **QSE Solve:** For valid QSE groups, `solveQSEAbundances` uses a Levenberg-Marquardt + * nonlinear solver (`Eigen::LevenbergMarquardt`) to find the equilibrium abundances of the + * "algebraic" species, holding the "seed" species constant. + * + * All calculations are cached using `QSECacheKey` to avoid re-partitioning and re-solving for + * similar thermodynamic conditions. + * + * @par Usage Example: + * @code + * // 1. Create a base engine (e.g., GraphEngine) + * gridfire::GraphEngine baseEngine(composition); + * + * // 2. Wrap it with the MultiscalePartitioningEngineView + * gridfire::MultiscalePartitioningEngineView multiscaleEngine(baseEngine); + * + * // 3. Before integration, update the view to partition the network + * // and find the initial equilibrium state. + * NetIn initialConditions = { .composition = composition, .temperature = 1e8, .density = 1e3 }; + * fourdst::composition::Composition equilibratedComp = multiscaleEngine.update(initialConditions); + * + * // 4. Use the multiscaleEngine for integration. It will use the cached QSE solution. + * // The integrator will call calculateRHSAndEnergy, etc. on the multiscaleEngine. + * auto Y_initial = multiscaleEngine.mapNetInToMolarAbundanceVector({equilibratedComp, ...}); + * auto derivatives = multiscaleEngine.calculateRHSAndEnergy(Y_initial, T9, rho); + * @endcode + * + * @implements DynamicEngine + * @implements EngineView + */ class MultiscalePartitioningEngineView final: public DynamicEngine, public EngineView { + /** + * @brief Type alias for a QSE partition. + * + * A QSE partition is a tuple containing the fast species, their indices, + * the slow species, and their indices. + */ typedef std::tuple, std::vector, std::vector, std::vector> QSEPartition; public: + /** + * @brief Constructs a MultiscalePartitioningEngineView. + * + * @param baseEngine The underlying GraphEngine to which this view delegates calculations. + * It must be a `GraphEngine` and not a more general `DynamicEngine` + * because this view relies on its specific implementation details. + */ explicit MultiscalePartitioningEngineView(GraphEngine& baseEngine); + /** + * @brief Gets the list of species in the network. + * @return A const reference to the vector of `Species` objects representing all species + * in the underlying base engine. This view does not alter the species list itself, + * only how their abundances are evolved. + */ [[nodiscard]] const std::vector & getNetworkSpecies() const override; + /** + * @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 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 + * `StaleEngineError` if the engine's QSE cache does not contain a solution + * for the given state. + * + * @purpose To compute the time derivatives for the ODE solver. This implementation + * modifies the derivatives from the base engine to enforce the QSE condition. + * + * @how It first performs a lookup in the QSE abundance cache (`m_qse_abundance_cache`). + * If a cache hit occurs, it calls the base engine's `calculateRHSAndEnergy`. It then + * manually sets the time derivatives (`dydt`) of all identified algebraic species to zero, + * effectively removing their differential equations from the system being solved. + * + * @pre The engine must have been updated via `update()` or `equilibrateNetwork()` for the + * current thermodynamic conditions, so that a valid entry exists in the QSE cache. + * @post The returned derivatives will have `dydt=0` for all algebraic species. + * + * @throws StaleEngineError If the QSE cache does not contain an entry for the given + * (T9, rho, Y_full). This indicates `update()` was not called recently enough. + */ [[nodiscard]] std::expected, expectations::StaleEngineError> calculateRHSAndEnergy( const std::vector &Y_full, double T9, double rho ) const override; + /** + * @brief Generates the Jacobian matrix for the current state. + * + * @param Y_full Vector of current molar abundances. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * + * @purpose To compute the Jacobian matrix required by implicit ODE solvers. + * + * @how It first performs a QSE cache lookup. On a hit, it delegates the full Jacobian + * calculation to the base engine. While this view could theoretically return a + * modified, sparser Jacobian reflecting the QSE constraints, the current implementation + * returns the full Jacobian from the base engine. The solver is expected to handle the + * algebraic constraints (e.g., via `dydt=0` from `calculateRHSAndEnergy`). + * + * @pre The engine must have a valid QSE cache entry for the given state. + * @post The base engine's internal Jacobian is updated. + * + * @throws exceptions::StaleEngineError If the QSE cache misses, as it cannot proceed + * without a valid partition. + */ void generateJacobianMatrix( const std::vector &Y_full, double T9, double rho ) const override; + /** + * @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. + * @return Value of the Jacobian matrix at (i_full, j_full). + * + * @purpose To provide Jacobian entries to an implicit solver. + * + * @how This method directly delegates to the base engine's `getJacobianMatrixEntry`. + * It does not currently modify the Jacobian to reflect the QSE algebraic constraints, + * as these are handled by setting `dY/dt = 0` in `calculateRHSAndEnergy`. + * + * @pre `generateJacobianMatrix()` must have been called for the current state. + */ [[nodiscard]] double getJacobianMatrixEntry( int i_full, int j_full ) const override; + /** + * @brief Generates the stoichiometry matrix for the network. + * + * @purpose To prepare the stoichiometry matrix for later queries. + * + * @how This method delegates directly to the base engine's `generateStoichiometryMatrix()`. + * The stoichiometry is based on the full, unpartitioned network. + */ void generateStoichiometryMatrix() override; + /** + * @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. + * @return Stoichiometric coefficient for the species in the reaction. + * + * @purpose To query the stoichiometric relationship between a species and a reaction. + * + * @how This method delegates directly to the base engine's `getStoichiometryMatrixEntry()`. + * + * @pre `generateStoichiometryMatrix()` must have been called. + */ [[nodiscard]] int getStoichiometryMatrixEntry( int speciesIndex, int reactionIndex ) 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 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). + * + * @purpose To compute the net rate of a single reaction. + * + * @how It first checks the QSE cache. On a hit, it retrieves the cached equilibrium + * abundances for the algebraic species. It creates a mutable copy of `Y_full`, + * overwrites the algebraic species abundances with the cached equilibrium values, + * and then calls the base engine's `calculateMolarReactionFlow` with this modified + * abundance vector. + * + * @pre The engine must have a valid QSE cache entry for the given state. + * @throws StaleEngineError If the QSE cache misses. + */ [[nodiscard]] double calculateMolarReactionFlow( const reaction::Reaction &reaction, const std::vector &Y_full, @@ -97,38 +330,170 @@ namespace gridfire { double rho ) const override; + /** + * @brief Gets the set of logical reactions in the network. + * + * @return A const reference to the `LogicalReactionSet` from the base engine, + * containing all reactions in the full network. + */ [[nodiscard]] const reaction::LogicalReactionSet & getNetworkReactions() const override; + /** + * @brief Sets the set of logical reactions in the network. + * + * @param reactions The set of logical reactions to use. + * + * @purpose To modify the reaction network. + * + * @how This operation is not supported by the `MultiscalePartitioningEngineView` as it + * would invalidate the partitioning logic. It logs a critical error and throws an + * exception. Network modifications should be done on the base engine before it is + * wrapped by this view. + * + * @throws exceptions::UnableToSetNetworkReactionsError Always. + */ void setNetworkReactions( const reaction::LogicalReactionSet &reactions ) override; + /** + * @brief Computes timescales for all species in the network. + * + * @param Y 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 `std::expected` containing a map from `Species` to their characteristic + * timescales (s) on success, or a `StaleEngineError` on failure. + * + * @purpose To get the characteristic timescale `Y / (dY/dt)` for each species. + * + * @how It delegates the calculation to the base engine. For any species identified + * as algebraic (in QSE), it manually sets their timescale to 0.0 to signify + * that they equilibrate instantaneously on the timescale of the solver. + * + * @pre The engine must have a valid QSE cache entry for the given state. + * @throws StaleEngineError If the QSE cache misses. + */ [[nodiscard]] std::expected, expectations::StaleEngineError> getSpeciesTimescales( const std::vector &Y, double T9, double rho ) const override; + /** + * @brief Computes destruction timescales for all species in the network. + * + * @param Y 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 `std::expected` containing a map from `Species` to their characteristic + * destruction timescales (s) on success, or a `StaleEngineError` on failure. + * + * @purpose To get the timescale for species destruction, which is used as the primary + * metric for network partitioning. + * + * @how It delegates the calculation to the base engine. For any species identified + * as algebraic (in QSE), it manually sets their timescale to 0.0. + * + * @pre The engine must have a valid QSE cache entry for the given state. + * @throws StaleEngineError If the QSE cache misses. + */ [[nodiscard]] std::expected, expectations::StaleEngineError> getSpeciesDestructionTimescales( const std::vector &Y, double T9, double rho ) const override; + /** + * @brief Updates the internal state of the engine, performing partitioning and QSE equilibration. + * + * @param netIn A struct containing the current network input: temperature, density, and composition. + * @return The new composition after QSE species have been brought to equilibrium. + * + * @purpose This is the main entry point for preparing the multiscale engine for use. It + * triggers the network partitioning and solves for the initial QSE abundances, caching the result. + * + * @how + * 1. It first checks the QSE cache. If a valid entry already exists for the input state, + * it returns the input composition, as no work is needed. + * 2. If the cache misses, it calls `equilibrateNetwork()`. + * 3. `equilibrateNetwork()` in turn calls `partitionNetwork()` to define the dynamic and + * algebraic species sets. + * 4. It then calls `solveQSEAbundances()` to compute the equilibrium abundances. + * 5. The resulting equilibrium abundances for the algebraic species are stored in the + * `m_qse_abundance_cache`. + * 6. A new `fourdst::composition::Composition` object reflecting the equilibrated state + * is created and returned. + * + * @pre The `netIn` struct should contain a valid physical state. + * @post The engine is partitioned (`m_dynamic_species`, `m_algebraic_species`, etc. are populated). + * The `m_qse_abundance_cache` is populated with the QSE solution for the given state. + * The returned composition reflects the new equilibrium. + */ fourdst::composition::Composition update( const NetIn &netIn ) override; + /** + * @brief Checks if the engine's internal state is stale relative to the provided conditions. + * + * @param netIn A struct containing the current network input. + * @return `true` if the engine is stale, `false` otherwise. + * + * @purpose To determine if `update()` needs to be called. + * + * @how It creates a `QSECacheKey` from the `netIn` data and checks for its + * existence in the `m_qse_abundance_cache`. A cache miss indicates the engine is + * stale because it does not have a valid QSE partition for the current conditions. + * It also queries the base engine's `isStale()` method. + */ bool isStale(const NetIn& netIn) override; + /** + * @brief Sets the electron screening model. + * + * @param model The type of screening model to use for reaction rate calculations. + * + * @how This method delegates directly to the base engine's `setScreeningModel()`. + */ void setScreeningModel( screening::ScreeningType model ) override; + /** + * @brief Gets the current electron screening model. + * + * @return The currently active screening model type. + * + * @how This method delegates directly to the base engine's `getScreeningModel()`. + */ [[nodiscard]] screening::ScreeningType getScreeningModel() const override; + /** + * @brief Gets the base engine. + * + * @return A const reference to the base engine. + */ const DynamicEngine & getBaseEngine() const override; + /** + * @brief Analyzes the connectivity of timescale pools. + * + * @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 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 + * single connected component. + * + * @purpose To merge timescale pools that are strongly connected by reactions, forming + * cohesive groups for QSE analysis. + * + * @how For each pool, it builds a reaction connectivity graph using `buildConnectivityGraph`. + * 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, @@ -136,16 +501,66 @@ namespace gridfire { double rho ) const; + /** + * @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 T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * + * @purpose To perform the core partitioning logic that identifies which species are "fast" + * (and can be treated algebraically) and which are "slow" (and must be integrated dynamically). + * + * @how + * 1. **`partitionByTimescale`**: Gets species destruction timescales from the base engine, + * sorts them, and looks for large gaps to create timescale "pools". + * 2. **`identifyMeanSlowestPool`**: The pool with the slowest average timescale is designated + * as the core set of dynamic species. + * 3. **`analyzeTimescalePoolConnectivity`**: The other (faster) pools are analyzed for + * reaction connectivity to form cohesive groups. + * 4. **`constructCandidateGroups`**: These connected groups are processed to identify "seed" + * species (dynamic species that feed the group) and "algebraic" species (the rest). + * 5. **`validateGroupsWithFluxAnalysis`**: The groups are validated by ensuring their internal + * reaction flux is much larger than the flux connecting them to the outside network. + * + * @pre The input state (Y, T9, rho) must be a valid physical state. + * @post The internal member variables `m_qse_groups`, `m_dynamic_species`, and + * `m_algebraic_species` (and their index maps) are populated with the results of the + * partitioning. + */ void partitionNetwork( const std::vector& Y, double T9, double rho ); + /** + * @brief Partitions the network based on timescales from a `NetIn` struct. + * + * @param netIn A struct containing the current network input. + * + * @purpose A convenience overload for `partitionNetwork`. + * + * @how It unpacks the `netIn` struct into `Y`, `T9`, and `rho` and then calls the + * primary `partitionNetwork` method. + */ void partitionNetwork( const NetIn& netIn ); + /** + * @brief Exports the network to a DOT file for visualization. + * + * @param filename The name of the DOT file to create. + * @param Y 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. + * + * @purpose To visualize the partitioned network graph. + * + * @how This method delegates the DOT file export to the base engine. It does not + * currently add any partitioning information to the output graph. + */ void exportToDot( const std::string& filename, const std::vector& Y, @@ -153,27 +568,112 @@ namespace gridfire { const double rho ) const; + /** + * @brief Gets the index of a species in the full network. + * + * @param species The species to get the index of. + * @return The index of the species in the base engine's network. + * + * @how This method delegates directly to the base engine's `getSpeciesIndex()`. + */ [[nodiscard]] int getSpeciesIndex(const fourdst::atomic::Species &species) const override; + /** + * @brief Maps a `NetIn` struct to a molar abundance vector for the full network. + * + * @param netIn A struct containing the current network input. + * @return A vector of molar abundances corresponding to the species order in the base engine. + * + * @how This method delegates directly to the base engine's `mapNetInToMolarAbundanceVector()`. + */ [[nodiscard]] std::vector mapNetInToMolarAbundanceVector(const NetIn &netIn) const override; + /** + * @brief Primes the engine with a specific species. + * + * @param netIn A struct containing the current network input. + * @return A `PrimingReport` struct containing information about the priming process. + * + * @purpose To prepare the network for ignition or specific pathway studies. + * + * @how This method delegates directly to the base engine's `primeEngine()`. The + * multiscale view does not currently interact with the priming process. + */ [[nodiscard]] PrimingReport primeEngine(const NetIn &netIn) override; + /** + * @brief Gets the fast species in the network. + * + * @return A vector of species identified as "fast" or "algebraic" by the partitioning. + * + * @purpose To allow external queries of the partitioning results. + * + * @how It returns a copy of the `m_algebraic_species` member vector. + * + * @pre `partitionNetwork()` must have been called. + */ [[nodiscard]] std::vector getFastSpecies() const; + /** + * @brief Gets the dynamic species in the network. + * + * @return A const reference to the vector of species identified as "dynamic" or "slow". + * + * @purpose To allow external queries of the partitioning results. + * + * @how It returns a const reference to the `m_dynamic_species` member vector. + * + * @pre `partitionNetwork()` must have been called. + */ [[nodiscard]] const std::vector& getDynamicSpecies() const; + /** + * @brief Equilibrates the network by partitioning and solving for QSE abundances. + * + * @param Y 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. + * + * @purpose A convenience method to run the full QSE analysis and get an equilibrated + * composition object as a result. + * + * @how It first calls `partitionNetwork()` with the given state to define the QSE groups. + * Then, it calls `solveQSEAbundances()` to compute the new equilibrium abundances for the + * algebraic species. Finally, it packs the resulting full abundance vector into a new + * `fourdst::composition::Composition` object and returns it. + * + * @pre The input state (Y, T9, rho) must be a valid physical state. + * @post The engine's internal partition is updated. A new composition object is returned. + */ fourdst::composition::Composition equilibrateNetwork( const std::vector &Y, double T9, double rho ); + /** + * @brief Equilibrates the network using QSE from a `NetIn` struct. + * + * @param netIn A struct containing the current network input. + * @return The equilibrated composition. + * + * @purpose A convenience overload for `equilibrateNetwork`. + * + * @how It unpacks the `netIn` struct into `Y`, `T9`, and `rho` and then calls the + * primary `equilibrateNetwork` method. + */ fourdst::composition::Composition equilibrateNetwork( const NetIn &netIn ); private: + /** + * @brief Struct representing a QSE group. + * + * @purpose A container to hold all information about a set of species that are potentially + * 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. @@ -181,46 +681,51 @@ namespace gridfire { std::set seed_indices; ///< Indices of dynamic species in this group. double mean_timescale; ///< Mean timescale of the group. + /** + * @brief Less-than operator for QSEGroup, used for sorting. + * @param other The other QSEGroup to compare to. + * @return True if this group's mean timescale is less than the other's. + */ bool operator<(const QSEGroup& other) const; + /** + * @brief Greater-than operator for QSEGroup. + * @param other The other QSEGroup to compare to. + * @return True if this group's mean timescale is greater than the other's. + */ bool operator>(const QSEGroup& other) const; + /** + * @brief Equality operator for QSEGroup. + * @param other The other QSEGroup to compare to. + * @return True if the sets of species indices are identical. + */ bool operator==(const QSEGroup& other) const; + /** + * @brief Inequality operator for QSEGroup. + * @param other The other QSEGroup to compare to. + * @return True if the sets of species indices are not identical. + */ bool operator!=(const QSEGroup& other) const; - - friend std::ostream& operator<<(std::ostream& os, const QSEGroup& group) { - os << "QSEGroup(species_indices: ["; - int count = 0; - for (const auto& idx : group.species_indices) { - os << idx; - if (count < group.species_indices.size() - 1) { - os << ", "; - } - count++; - } - count = 0; - os << "], is_in_equilibrium: " << group.is_in_equilibrium - << ", mean_timescale: " << group.mean_timescale - << " s, algebraic_indices: ["; - for (const auto& idx : group.algebraic_indices) { - os << idx; - if (count < group.algebraic_indices.size() - 1) { - os << ", "; - } - count++; - } - count = 0; - os << "], seed_indices: ["; - for (const auto& idx : group.seed_indices) { - os << idx; - if (count < group.seed_indices.size() - 1) { - os << ", "; - } - count++; - } - os << "])"; - return os; - } }; + /** + * @brief Functor for solving QSE abundances using Eigen's nonlinear optimization. + * + * @purpose This struct provides the objective function (`operator()`) and its Jacobian + * (`df`) to Eigen's Levenberg-Marquardt solver. The goal is to find the abundances + * of algebraic species that make their time derivatives (`dY/dt`) equal to zero. + * + * @how + * - **`operator()`**: Takes a vector `v_qse` (scaled abundances of algebraic species) as input. + * It constructs a full trial abundance vector `y_trial`, calls the base engine's + * `calculateRHSAndEnergy`, and returns the `dY/dt` values for the algebraic species. + * The solver attempts to drive this return vector to zero. + * - **`df`**: Computes the Jacobian of the objective function. It calls the base engine's + * `generateJacobianMatrix` and extracts the sub-matrix corresponding to the algebraic + * species. It applies the chain rule to account for the `asinh` scaling used on the + * abundances. + * + * The abundances are scaled using `asinh` to handle the large dynamic range and ensure positivity. + */ struct EigenFunctor { using InputType = Eigen::Matrix; using OutputType = Eigen::Matrix; @@ -230,13 +735,41 @@ namespace gridfire { ValuesAtCompileTime = Eigen::Dynamic }; + /** + * @brief Pointer to the MultiscalePartitioningEngineView instance. + */ MultiscalePartitioningEngineView* m_view; + /** + * @brief Indices of the species to solve for in the QSE group. + */ const std::vector& m_qse_solve_indices; + /** + * @brief Initial abundances of all species in the full network. + */ const std::vector& m_Y_full_initial; + /** + * @brief Temperature in units of 10^9 K. + */ const double m_T9; + /** + * @brief Density in g/cm^3. + */ const double m_rho; + /** + * @brief Scaling factors for the species abundances, used to improve solver stability. + */ const Eigen::VectorXd& m_Y_scale; + /** + * @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 T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @param Y_scale Scaling factors for the species abundances. + */ EigenFunctor( MultiscalePartitioningEngineView& view, const std::vector& qse_solve_indices, @@ -252,14 +785,40 @@ namespace gridfire { m_rho(rho), m_Y_scale(Y_scale) {} + /** + * @brief Gets the number of output values from the functor (size of the residual vector). + * @return The number of algebraic species being solved. + */ [[nodiscard]] int values() const { return m_qse_solve_indices.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]] int inputs() const { return m_qse_solve_indices.size(); } + /** + * @brief Evaluates the functor's residual vector `f_qse = dY_alg/dt`. + * @param v_qse The input vector of scaled algebraic abundances. + * @param f_qse The output residual vector. + * @return 0 on success. + */ int operator()(const InputType& v_qse, OutputType& f_qse) const; + /** + * @brief Evaluates the Jacobian of the functor, `J_qse = d(f_qse)/d(v_qse)`. + * @param v_qse The input vector of scaled algebraic abundances. + * @param J_qse The output Jacobian matrix. + * @return 0 on success. + */ int df(const InputType& v_qse, JacobianType& J_qse) const; }; + /** + * @brief Struct for tracking cache statistics. + * + * @purpose A simple utility to monitor the performance of the QSE cache by counting + * hits and misses for various engine operations. + */ struct CacheStats { enum class operators { CalculateRHSAndEnergy, @@ -271,6 +830,9 @@ namespace gridfire { All }; + /** + * @brief Map from operators to their string names for logging. + */ std::map operatorsNameMap = { {operators::CalculateRHSAndEnergy, "calculateRHSAndEnergy"}, {operators::GenerateJacobianMatrix, "generateJacobianMatrix"}, @@ -280,9 +842,18 @@ namespace gridfire { {operators::Other, "other"} }; + /** + * @brief Total number of cache hits. + */ size_t m_hit = 0; + /** + * @brief Total number of cache misses. + */ size_t m_miss = 0; + /** + * @brief Map from operators to the number of cache hits for that operator. + */ std::map m_operatorHits = { {operators::CalculateRHSAndEnergy, 0}, {operators::GenerateJacobianMatrix, 0}, @@ -292,72 +863,160 @@ namespace gridfire { {operators::Other, 0} }; + /** + * @brief Map from operators to the number of cache misses for that operator. + */ + std::map m_operatorMisses = { + {operators::CalculateRHSAndEnergy, 0}, + {operators::GenerateJacobianMatrix, 0}, + {operators::CalculateMolarReactionFlow, 0}, + {operators::GetSpeciesTimescales, 0}, + {operators::GetSpeciesDestructionTimescales, 0}, + {operators::Other, 0} + }; - friend std::ostream& operator<<(std::ostream& os, const CacheStats& stats) { - os << "CacheStats(hit: " << stats.m_hit << ", miss: " << stats.m_miss << ")"; - return os; - } + /** + * @brief Increments the hit counter for a given operator. + * @param op The operator that resulted in a cache hit. + * @throws std::invalid_argument if `op` is `All`. + */ + void hit(const operators op=operators::Other); + /** + * @brief Increments the miss counter for a given operator. + * @param op The operator that resulted in a cache miss. + * @throws std::invalid_argument if `op` is `All`. + */ + void miss(const operators op=operators::Other); - void hit(const operators op=operators::Other) { - if (op==operators::All) { - throw std::invalid_argument("Cannot use 'All' as an operator for hit/miss."); - } - m_hit++; - m_operatorHits[op]++; - } - void miss(const operators op=operators::Other) { - if (op==operators::All) { - throw std::invalid_argument("Cannot use 'All' as an operator for hit/miss."); - } - m_miss++; - m_operatorHits[op]++; - } - - [[nodiscard]] size_t hits(const operators op=operators::All) const { - if (op==operators::All) { - return m_hit; - } - return m_operatorHits.at(op); - } - [[nodiscard]] size_t misses(const operators op=operators::All) const { - if (op==operators::All) { - return m_miss; - } - return m_operatorHits.at(op); - } + /** + * @brief Gets the number of hits for a specific operator or all operators. + * @param op The operator to get the number of hits for. Defaults to `All`. + * @return The number of hits. + */ + [[nodiscard]] size_t hits(const operators op=operators::All) const; + /** + * @brief Gets the number of misses for a specific operator or all operators. + * @param op The operator to get the number of misses for. Defaults to `All`. + * @return The number of misses. + */ + [[nodiscard]] size_t misses(const operators op=operators::All) const; }; private: + /** + * @brief Logger instance for logging messages. + */ quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); - GraphEngine& m_baseEngine; ///< The base engine to which this view delegates calculations. - std::vector m_qse_groups; ///< The list of identified equilibrium groups. - std::vector m_dynamic_species; ///< The simplified set of species presented to the solver. - std::vector m_dynamic_species_indices; ///< Indices mapping the dynamic species back to the master engine's list. - std::vector m_algebraic_species; ///< Species that are algebraic in the QSE groups. - std::vector m_algebraic_species_indices; ///< Indices of algebraic species in the full network. + /** + * @brief The base engine to which this view delegates calculations. + */ + GraphEngine& m_baseEngine; + /** + * @brief The list of identified equilibrium groups. + */ + std::vector m_qse_groups; + /** + * @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; + /** + * @brief Indices of algebraic species in the full network. + */ + std::vector m_algebraic_species_indices; - std::vector m_activeSpeciesIndices; ///< Indices of active species in the full network. - std::vector m_activeReactionIndices; ///< Indices of active reactions in the full network. + /** + * @brief Indices of all species considered active in the current partition (dynamic + algebraic). + */ + std::vector m_activeSpeciesIndices; + /** + * @brief Indices of all reactions involving only active species. + */ + std::vector m_activeReactionIndices; // TODO: Enhance the hashing for the cache to consider not just T and rho but also the current abundance in some careful way that automatically ignores small changes (i.e. network should only be repartitioned sometimes) - std::unordered_map> m_qse_abundance_cache; ///< Cache for QSE abundances based on T9 and rho. - mutable CacheStats m_cacheStats; ///< Statistics for the QSE abundance cache. + /** + * @brief Cache for QSE abundances based on T9, rho, and Y. + * + * @purpose This is the core of the caching mechanism. It stores the results of QSE solves + * to avoid re-computation. The key is a `QSECacheKey` which hashes the thermodynamic + * state, and the value is the vector of solved molar abundances for the algebraic species. + */ + mutable std::unordered_map> m_qse_abundance_cache; + /** + * @brief Statistics for the QSE abundance cache. + */ + mutable CacheStats m_cacheStats; private: + /** + * @brief Partitions the network by timescale. + * + * @param Y_full 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 + * timescale pool. + * + * @purpose To group species into "pools" based on their destruction timescales. + * + * @how It retrieves all species destruction timescales from the base engine, sorts them, + * and then iterates through the sorted list, creating a new pool whenever it detects + * 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, double T9, double rho ) const; + /** + * @brief Builds a connectivity graph from a set of fast reaction indices. + * + * @param fast_reaction_indices A set of indices for reactions considered "fast". + * @return An unordered map representing the adjacency list of the connectivity graph, + * where keys are species indices and values are vectors of connected species indices. + * + * @purpose To represent the reaction pathways among a subset of reactions. + * + * @how It iterates through the specified fast reactions. For each reaction, it creates + * a two-way edge in the graph between every reactant and every product, signifying + * that mass can flow between them. + */ std::unordered_map> buildConnectivityGraph( const std::unordered_set &fast_reaction_indices ) const; + /** + * @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 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. + * + * @purpose To ensure that a candidate QSE group is truly in equilibrium by checking that + * the reaction fluxes *within* the group are much larger than the fluxes + * *leaving* the group. + * + * @how For each candidate group, it calculates the sum of all internal reaction fluxes and + * the sum of all external (bridge) reaction fluxes. If the ratio of internal to external + * flux exceeds a configurable threshold, the group is considered valid and is added + * to the returned vector. + */ std::vector validateGroupsWithFluxAnalysis( const std::vector &candidate_groups, const std::vector& Y, @@ -365,12 +1024,45 @@ namespace gridfire { double rho ) const; + /** + * @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 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. + * + * @purpose To find the equilibrium abundances of the algebraic species that satisfy + * the QSE conditions. + * + * @how It uses the Levenberg-Marquardt algorithm via Eigen's `LevenbergMarquardt` class. + * The problem is defined by the `EigenFunctor` which computes the residuals and + * Jacobian for the QSE equations. + * + * @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, double T9, double rho ); + /** + * @brief Identifies the pool with the slowest mean timescale. + * + * @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 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. + * + * @purpose To identify the core set of dynamic species that will not be part of any QSE group. + * + * @how It calculates the geometric mean of the destruction timescales for all species in each + * pool and returns the index of the pool with the maximum mean timescale. + */ size_t identifyMeanSlowestPool( const std::vector>& pools, const std::vector &Y, @@ -378,14 +1070,46 @@ namespace gridfire { double rho ) const; + /** + * @brief Builds a connectivity graph from a species pool. + * + * @param species_pool A vector of species indices representing a species pool. + * @return An unordered map representing the adjacency list of the connectivity graph. + * + * @purpose To find reaction connections within a specific group of species. + * + * @how It iterates through all reactions in the base engine. If a reaction involves + * at least two distinct species from the input `species_pool` (one as a reactant + * 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 ) const; + /** + * @brief Constructs candidate QSE groups from connected timescale pools. + * + * @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 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. + * + * @how For each input pool, it identifies "bridge" reactions that connect the pool to + * species outside the pool. The reactants of these bridge reactions that are *not* in the + * pool are identified as "seed" species. The original pool members are the "algebraic" + * species. It then bundles the seed and algebraic species into a `QSEGroup` struct. + * + * @pre The `candidate_pools` should be connected components from `analyzeTimescalePoolConnectivity`. + * @post A list of candidate `QSEGroup` objects is returned. + */ std::vector constructCandidateGroups( const std::vector>& candidate_pools, const std::vector& Y, - double T9, double rho + double T9, + double rho ) const; }; } diff --git a/src/include/gridfire/screening/screening_weak.h b/src/include/gridfire/screening/screening_weak.h index f5e25dea..a81ea63a 100644 --- a/src/include/gridfire/screening/screening_weak.h +++ b/src/include/gridfire/screening/screening_weak.h @@ -78,7 +78,7 @@ namespace gridfire::screening { ) const override; private: /// @brief Logger instance for recording trace and debug information. - quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); + [[maybe_unused]] quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); private: /** diff --git a/src/lib/engine/engine_graph.cpp b/src/lib/engine/engine_graph.cpp index 6f9ad5da..6f97e11b 100644 --- a/src/lib/engine/engine_graph.cpp +++ b/src/lib/engine/engine_graph.cpp @@ -41,8 +41,8 @@ namespace gridfire { const partition::PartitionFunction& partitionFunction, const BuildDepthType buildDepth) : m_reactions(build_reaclib_nuclear_network(composition, buildDepth, false)), - m_partitionFunction(partitionFunction.clone()), - m_depth(buildDepth) + m_depth(buildDepth), + m_partitionFunction(partitionFunction.clone()) { syncInternalMaps(); } diff --git a/src/lib/engine/views/engine_multiscale.cpp b/src/lib/engine/views/engine_multiscale.cpp index 7612c50e..8440b610 100644 --- a/src/lib/engine/views/engine_multiscale.cpp +++ b/src/lib/engine/views/engine_multiscale.cpp @@ -5,6 +5,7 @@ #include "fourdst/logging/logging.h" +#include #include #include #include @@ -482,18 +483,7 @@ namespace gridfire { // --- 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); - LOG_TRACE_L1( - m_logger, - "Found {} candidate QSE groups for further analysis: {}", - candidate_groups.size(), - [&]() -> std::string { - std::stringstream ss; - for (const auto& group : candidate_groups) { - ss << group << " "; - } - return ss.str(); - }() - ); + LOG_TRACE_L1(m_logger, "Found {} candidate QSE groups for further analysis", candidate_groups.size()); LOG_TRACE_L1(m_logger, "Validating candidate groups with flux analysis..."); const std::vector validated_groups = validateGroupsWithFluxAnalysis(candidate_groups, Y, T9, rho); @@ -1401,9 +1391,7 @@ namespace gridfire { // Add clique for (const size_t& u : intersection) { - const auto& uSpecies = m_baseEngine.getNetworkSpecies()[u]; for (const size_t& v : intersection) { - const auto & vSpecies = m_baseEngine.getNetworkSpecies()[v]; if (u != v) { // Avoid self-loops connectivity_graph[u].push_back(v); } @@ -1623,4 +1611,35 @@ namespace gridfire { return !(*this == other); } + void MultiscalePartitioningEngineView::CacheStats::hit(const operators op) { + if (op == operators::All) { + throw std::invalid_argument("Cannot use 'ALL' as an operator for a hit"); + } + + m_hit ++; + m_operatorHits[op]++; + } + void MultiscalePartitioningEngineView::CacheStats::miss(const operators op) { + if (op == operators::All) { + throw std::invalid_argument("Cannot use 'ALL' as an operator for a miss"); + } + + m_miss ++; + m_operatorMisses[op]++; + } + + size_t MultiscalePartitioningEngineView::CacheStats::hits(const operators op) const { + if (op == operators::All) { + return m_hit; + } + return m_operatorHits.at(op); + } + + size_t MultiscalePartitioningEngineView::CacheStats::misses(const operators op) const { + if (op == operators::All) { + return m_miss; + } + return m_operatorMisses.at(op); + } + } diff --git a/src/lib/io/network_file.cpp b/src/lib/io/network_file.cpp index 7e3e460f..13081f30 100644 --- a/src/lib/io/network_file.cpp +++ b/src/lib/io/network_file.cpp @@ -53,10 +53,8 @@ namespace gridfire::io { ParsedNetworkData parsed; std::string line; - int line_number = 0; while (std::getline(file, line)) { - line_number++; - LOG_TRACE_L3(m_logger, "Parsing reaction list file {}, line {}: {}", filename, line_number, line); + LOG_TRACE_L3(m_logger, "Parsing reaction list file {}, line: {}", filename, line); const size_t comment_pos = line.find('#'); if (comment_pos != std::string::npos) { diff --git a/src/python/engine/trampoline/py_engine.cpp b/src/python/engine/trampoline/py_engine.cpp index 41082fc0..0c8fb827 100644 --- a/src/python/engine/trampoline/py_engine.cpp +++ b/src/python/engine/trampoline/py_engine.cpp @@ -14,11 +14,24 @@ namespace py = pybind11; const std::vector& PyEngine::getNetworkSpecies() const { - PYBIND11_OVERRIDE_PURE( - std::vector, - gridfire::Engine, /* Base class */ - getNetworkSpecies - ); + /* + * Acquire the GIL (Global Interpreter Lock) for thread safety + * with the Python interpreter. + */ + py::gil_scoped_acquire gil; + + /* + * get_override() looks for a Python method that overrides this C++ one. + */ + py::function override = py::get_override(this, "getNetworkSpecies"); + + if (override) { + py::object result = override(); + m_species_cache = result.cast>(); + return m_species_cache; + } + + py::pybind11_fail("Tried to call pure virtual function \"DynamicEngine::getNetworkSpecies\""); } std::expected, gridfire::expectations::StaleEngineError> PyEngine::calculateRHSAndEnergy(const std::vector &Y, double T9, double rho) const { @@ -35,11 +48,24 @@ std::expected, gridfire::expectations::StaleEn ///////////////////////////////////// const std::vector& PyDynamicEngine::getNetworkSpecies() const { - PYBIND11_OVERRIDE_PURE( - std::vector, - gridfire::DynamicEngine, /* Base class */ - getNetworkSpecies - ); + /* + * Acquire the GIL (Global Interpreter Lock) for thread safety + * with the Python interpreter. + */ + py::gil_scoped_acquire gil; + + /* + * get_override() looks for a Python method that overrides this C++ one. + */ + py::function override = py::get_override(this, "getNetworkSpecies"); + + if (override) { + py::object result = override(); + m_species_cache = result.cast>(); + return m_species_cache; + } + + py::pybind11_fail("Tried to call pure virtual function \"DynamicEngine::getNetworkSpecies\""); } std::expected, gridfire::expectations::StaleEngineError> PyDynamicEngine::calculateRHSAndEnergy(const std::vector &Y, double T9, double rho) const { PYBIND11_OVERRIDE_PURE( diff --git a/src/python/engine/trampoline/py_engine.h b/src/python/engine/trampoline/py_engine.h index 2a303125..4fe6d697 100644 --- a/src/python/engine/trampoline/py_engine.h +++ b/src/python/engine/trampoline/py_engine.h @@ -13,9 +13,12 @@ class PyEngine final : public gridfire::Engine { public: const std::vector& getNetworkSpecies() const override; std::expected,gridfire::expectations::StaleEngineError> calculateRHSAndEnergy(const std::vector &Y, double T9, double rho) const override; +private: + mutable std::vector m_species_cache; }; class PyDynamicEngine final : public gridfire::DynamicEngine { +public: const std::vector& getNetworkSpecies() const override; std::expected,gridfire::expectations::StaleEngineError> calculateRHSAndEnergy(const std::vector &Y, double T9, double rho) const override; void generateJacobianMatrix(const std::vector &Y_dynamic, double T9, double rho) const override; @@ -41,6 +44,10 @@ class PyDynamicEngine final : public gridfire::DynamicEngine { void rebuild(const fourdst::composition::Composition& comp, gridfire::BuildDepthType depth) override { throw std::logic_error("Setting network depth not supported by this engine."); } +private: + mutable std::vector m_species_cache; + + }; class PyEngineView final : public gridfire::EngineView { diff --git a/subprojects/fourdst.wrap b/subprojects/fourdst.wrap index 28348f18..328f8539 100644 --- a/subprojects/fourdst.wrap +++ b/subprojects/fourdst.wrap @@ -1,4 +1,4 @@ [wrap-git] url = https://github.com/4D-STAR/fourdst -revision = v0.5.0 -depth = 1 \ No newline at end of file +revision = v0.5.1 +depth = 1 diff --git a/tests/python/test.py b/tests/python/test.py index a3a1cc24..286fda0f 100644 --- a/tests/python/test.py +++ b/tests/python/test.py @@ -16,7 +16,7 @@ netIn = NetIn() netIn.composition = comp netIn.temperature = 1.5e7 netIn.density = 1.5e2 -netIn.tMax = 3e14 +netIn.tMax = 3e17 netIn.dt0 = 1e-12 baseEngine = GraphEngine(comp, 2)