Calculates the molar reaction flow rate for all reactions in the full network.
+
This method iterates through all reactions in the base engine's network and calculates their molar flow rates based on the provided network input conditions (temperature, density, and composition). It also constructs a vector of molar abundances for all species in the full network.
+
Parameters
+
+
netIn
The current network input, containing temperature, density, and composition.
+
out_Y_Full
A vector that will be populated with the molar abundances of all species in the full network.
+
+
+
+
Returns
A vector of ReactionFlow structs, each containing a pointer to a reaction and its calculated flow rate.
+
Algorithm:
+
Clears and reserves space in out_Y_Full.
+
Iterates through all species in the base engine's network.
+
For each species, it retrieves the molar abundance from netIn.composition. If the species is not found, its abundance is set to 0.0.
+
Converts the temperature from Kelvin to T9.
+
Iterates through all reactions in the base engine's network.
+
For each reaction, it calls the base engine's calculateMolarReactionFlow to get the flow rate.
+
Stores the reaction pointer and its flow rate in a ReactionFlow struct and adds it to the returned vector.
Culls reactions from the network based on their flow rates.
+
This method filters the list of all reactions, keeping only those with a flow rate above an absolute culling threshold. The threshold is calculated by multiplying the maximum flow rate by a relative culling threshold read from the configuration.
+
Parameters
+
+
allFlows
A vector of all reactions and their flow rates.
+
reachableSpecies
A set of all species reachable from the initial fuel.
+
Y_full
A vector of molar abundances for all species in the full network.
+
maxFlow
The maximum reaction flow rate in the network.
+
+
+
+
Returns
A vector of pointers to the reactions that have been kept after culling.
+
Algorithm:
+
Retrieves the RelativeCullingThreshold from the configuration.
+
Calculates the absoluteCullingThreshold by multiplying maxFlow with the relative threshold.
+
Iterates through allFlows.
+
A reaction is kept if its flowRate is greater than the absoluteCullingThreshold.
+
The pointers to the kept reactions are stored in a vector and returned.
Finalizes the set of active species and reactions.
+
This method takes the final list of culled reactions and populates the m_activeReactions and m_activeSpecies members. The active species are determined by collecting all reactants and products from the final reactions. The active species list is then sorted by mass.
+
Parameters
+
+
finalReactions
A vector of pointers to the reactions to be included in the active set.
+
+
+
+
Postcondition
+
m_activeReactions is cleared and populated with the reactions from finalReactions.
+
m_activeSpecies is cleared and populated with all unique species present in finalReactions.
Finds all species that are reachable from the initial fuel through the reaction network.
+
This method performs a connectivity analysis to identify all species that can be produced starting from the initial fuel species. A species is considered part of the initial fuel if its mass fraction is above a certain threshold (ABUNDANCE_FLOOR).
+
Parameters
+
+
netIn
The current network input, containing the initial composition.
+
+
+
+
Returns
An unordered set of all reachable species.
+
Algorithm:
+
Initializes a set reachable and a queue to_visit with the initial fuel species.
+
Iteratively processes the reaction network until no new species can be reached.
+
In each pass, it iterates through all reactions in the base engine's network.
+
If all reactants of a reaction are in the reachable set, all products of that reaction are added to the reachable set.
+
The process continues until a full pass over all reactions does not add any new species to the reachable set.
The type of screening model to use for reaction rate calculations.
+
+
+
+
This method allows changing the screening model at runtime. Screening corrections account for the electrostatic shielding of nuclei by electrons, which affects reaction rates in dense stellar plasmas.
A struct containing the current network input, such as temperature, density, and composition.
+
+
+
+
This method is intended to be implemented by derived classes to update their internal state based on the provided network conditions. For example, an adaptive engine might use this to re-evaluate which reactions and species are active. For other engines that do not support manually updating, this method might do nothing.
An engine view that uses a user-defined reaction network from a file.
+
This class implements an EngineView that restricts the reaction network to a specific set of reactions defined in an external file. It acts as a filter or a view on a larger, more comprehensive base engine. The file provides a list of reaction identifiers, and this view will only consider those reactions and the species involved in them.
+
This is useful for focusing on a specific sub-network for analysis, debugging, or performance reasons, without modifying the underlying full network.
+
The view maintains mappings between the indices of its active (defined) species and reactions and the corresponding indices in the full network of the base engine. All calculations are delegated to the base engine after mapping the inputs from the view's context to the full network context, and the results are mapped back.
Builds the active species and reaction sets from a file.
+
This method uses the provided parser to read reaction names from the given file. It then finds these reactions in the base engine's full network and populates the m_activeReactions and m_activeSpecies members. Finally, it constructs the index maps for the active sets.
+
Parameters
+
+
fileName
The path to the network definition file.
+
+
+
+
Postcondition
+
m_activeReactions and m_activeSpecies are populated.
+
m_speciesIndexMap and m_reactionIndexMap are constructed.
+
m_isStale is set to false.
+
+
+
Exceptions
+
+
std::runtime_error
If a reaction from the file is not found in the base engine.
This function must be implemented by derived classes to compute the time derivatives of all species and the specific nuclear energy generation rate for the current state.
Generate the Jacobian matrix for the current state.
+
Generates the Jacobian matrix for the active species.
Parameters
-
Y
Vector of current abundances.
-
T9
Temperature in units of 10^9 K.
-
rho
Density in g/cm^3.
+
Y_defined
A vector of abundances for the active species.
+
T9
The temperature in units of 10^9 K.
+
rho
The density in g/cm^3.
+
+
+
+
Exceptions
+
+
std::runtime_error
If the view is stale.
-
This method must compute and store the Jacobian matrix (∂(dY/dt)_i/∂Y_j) for the current state. The matrix can then be accessed via getJacobianMatrixEntry().
This method must be implemented by derived classes to provide access to the base engine. The returned reference should remain valid for the lifetime of the EngineView.
Compute timescales for all species in the network.
+
Computes timescales for all active species in the network.
Parameters
-
Y
Vector of current abundances.
+
Y_defined
Vector of current abundances for the active species.
T9
Temperature in units of 10^9 K.
rho
Density in g/cm^3.
Returns
Map from Species to their characteristic timescales (s).
-
This method estimates the timescale for abundance change of each species, which can be used for timestep control, diagnostics, and reaction network culling.
This method checks if the view is stale (e.g., after setNetworkFile was called). If it is, it rebuilds the active network from the currently set file. The netIn parameter is not used by this implementation but is required by the interface.
+
Parameters
+
+
netIn
The current network input (unused).
+
+
+
+
Postcondition
If the view was stale, it is rebuilt and is no longer stale.
Flag to enable or disable using precomputed reactions for efficiency. Mathematically, this should not change the results. Generally end users should not need to change this.
This method calculates the time derivatives of all species and the specific nuclear energy generation rate for the current state using automatic differentiation.
This method calculates the time derivatives of all species and the specific nuclear energy generation rate for the current state using double precision arithmetic.
This method reserves space for the Jacobian matrix, which is used to store the partial derivatives of the right-hand side of the ODE with respect to the species abundances.
The type of screening model to use for reaction rate calculations.
+
+
+
+
This method allows changing the screening model at runtime. Screening corrections account for the electrostatic shielding of nuclei by electrons, which affects reaction rates in dense stellar plasmas.
This method synchronizes the internal maps used by the engine, including the species map, reaction ID map, and species-to-index map. It also generates the stoichiometry matrix and records the AD tape.
A struct containing the current network input, such as temperature, density, and composition.
+
+
+
+
This method is intended to be implemented by derived classes to update their internal state based on the provided network conditions. For example, an adaptive engine might use this to re-evaluate which reactions and species are active. For other engines that do not support manually updating, this method might do nothing.
This method validates the composition against the current reaction set. If the composition is not compatible with the reaction set, the reaction set is rebuilt from the composition.
True if all reactions conserve mass and charge, false otherwise.
This method checks that all reactions in the network conserve mass and charge. If any reaction does not conserve mass or charge, an error message is logged and false is returned.
Flag to enable or disable using precomputed reactions for efficiency. Mathematically, this should not change the results. Generally end users should not need to change this.
This class defines the interface for parsing files that contain reaction network definitions. Derived classes must implement the parse method to handle specific file formats.
Parses a network file and returns the parsed data.
+
This is a pure virtual function that must be implemented by derived classes. It takes a filename as input and returns a ParsedNetworkData struct containing the information extracted from the file.
This method reads the specified file line by line. It trims whitespace from each line, ignores lines that are empty or start with a '#' comment character, and stores the remaining lines as reaction names.
A screening model that applies no screening effect.
+
This class implements the ScreeningModel interface but returns a screening factor of 1.0 for all reactions, regardless of the plasma conditions. It represents the case of bare, unscreened nuclei and serves as a baseline or can be used when screening effects are negligible or intentionally ignored.
Template implementation for calculating screening factors.
+
Template implementation for the bare screening model.
+
This private helper function contains the core logic for both the double and ADDouble versions of calculateScreeningFactors. It is templated to handle both numeric types seamlessly.
+
Template Parameters
+
+
T
The numeric type, either double or CppAD::AD<double>.
+
+
+
+
Parameters
+
+
reactions
The set of reactions for which to calculate factors.
+
species
A vector of all atomic species (unused).
+
Y
A vector of molar abundances (unused).
+
T9
The temperature (unused).
+
rho
The density (unused).
+
+
+
+
Returns
A vector of type T with all elements initialized to 1.0.
+
This function provides the actual implementation for calculateFactors_impl. It creates a vector of the appropriate numeric type (T) and size, and initializes all its elements to 1.0, representing no screening.
+
Template Parameters
+
+
T
The numeric type, either double or CppAD::AD<double>.
+
+
+
+
Parameters
+
+
reactions
The set of reactions, used to determine the size of the output vector.
+
species
Unused parameter.
+
Y
Unused parameter.
+
T9
Unused parameter.
+
rho
Unused parameter.
+
+
+
+
Returns
A std::vector<T> of the same size as reactions, with all elements set to 1.0.
Calculates screening factors for AD types, which are always 1.0.
+
This implementation returns a vector of AD-typed screening factors where every element is 1.0. This is the automatic differentiation-compatible version.
+
Parameters
+
+
reactions
The set of logical reactions in the network.
+
species
A vector of all atomic species (unused).
+
Y
A vector of the molar abundances as AD types (unused).
+
T9
The temperature as an AD type (unused).
+
rho
The plasma density as an AD type (unused).
+
+
+
+
Returns
A vector of ADDouble, with each element being 1.0, of the same size as the reactions set.
Calculates screening factors, which are always 1.0.
+
This implementation returns a vector of screening factors where every element is 1.0, effectively applying no screening correction to the reaction rates.
+
Parameters
+
+
reactions
The set of logical reactions in the network.
+
species
A vector of all atomic species (unused).
+
Y
A vector of the molar abundances (unused).
+
T9
The temperature (unused).
+
rho
The plasma density (unused).
+
+
+
+
Returns
A vector of doubles, with each element being 1.0, of the same size as the reactions set.
+
Algorithm The function simply creates and returns a std::vector<double> of the same size as the input reactions set, with all elements initialized to 1.0.
An abstract base class for plasma screening models.
+
This class defines the interface for models that calculate the enhancement factor for nuclear reaction rates due to the electrostatic screening of interacting nuclei by the surrounding plasma. Concrete implementations of this class will provide specific screening prescriptions (e.g., WEAK, BARE, STRONG, etc.).
+
The interface provides methods for calculating screening factors for both standard double-precision inputs and for CppAD's automatic differentiation types, allowing the screening contributions to be included in Jacobian calculations.
Ensures that derived class destructors are called correctly.
+
Member Function Documentation
@@ -222,6 +237,21 @@ Public Member Functions
+
Calculates screening factors using CppAD types for automatic differentiation.
+
This is a pure virtual function that provides an overload of calculateScreeningFactors for use with CppAD. It allows the derivatives of the screening factors with respect to abundances, temperature, and density to be computed automatically.
+
Parameters
+
+
reactions
The set of logical reactions in the network.
+
species
A vector of all atomic species involved in the network.
+
Y
A vector of the molar abundances (mol/g) for each species, as AD types.
+
T9
The temperature in units of 10^9 K, as an AD type.
+
rho
The plasma density in g/cm^3, as an AD type.
+
+
+
+
Returns
A vector of screening factors (dimensionless), as AD types.
+
Note This method is essential for including the effects of screening in the Jacobian matrix of the reaction network.
Calculates screening factors for a set of reactions.
+
This is a pure virtual function that must be implemented by derived classes. It computes the screening enhancement factor for each reaction in the provided set based on the given plasma conditions.
+
Parameters
+
+
reactions
The set of logical reactions in the network.
+
species
A vector of all atomic species involved in the network.
+
Y
A vector of the molar abundances (mol/g) for each species.
+
T9
The temperature in units of 10^9 K.
+
rho
The plasma density in g/cm^3.
+
+
+
+
Returns
A vector of screening factors (dimensionless), one for each reaction in the reactions set, in the same order.
+
Pre-conditions
+
The size of the Y vector must match the size of the species vector.
+
T9 and rho must be positive.
+
+
Post-conditions
+
The returned vector will have the same size as the reactions set.
+
Each element in the returned vector will be >= 1.0.
+
+
Usage
// Assume 'model' is a std::unique_ptr<ScreeningModel> to a concrete implementation
+
// and other parameters (reactions, species, Y, T9, rho) are initialized.
Implements the weak screening model based on the Debye-Hückel approximation.
+
This class provides a concrete implementation of the ScreeningModel interface for the weak screening regime, following the formulation of Salpeter (1954). This approach applies the Debye-Hückel theory to model the electrostatic shielding of nuclei in a plasma. It is applicable to non-degenerate, non-relativistic plasmas where thermal energy dominates the electrostatic potential energy.
Template implementation for calculating weak screening factors.
+
Core implementation of the weak screening calculation (Debye-Hückel model).
+
This private helper function contains the core logic for calculating weak screening factors. It is templated to handle both double and CppAD::AD<double> numeric types, avoiding code duplication.
+
Template Parameters
+
+
T
The numeric type, either double or CppAD::AD<double>.
+
+
+
+
Parameters
+
+
reactions
The set of reactions.
+
species
A vector of all species in the network.
+
Y
A vector of molar abundances.
+
T9
The temperature in 10^9 K.
+
rho
The density in g/cm^3.
+
+
+
+
Returns
A vector of screening factors of type T.
+
This function calculates the screening factor exp(H_12) for each reaction, based on the Debye-Hückel approximation as formulated by Salpeter (1954).
+
Template Parameters
+
+
T
The numeric type (double or CppAD::AD<double>).
+
+
+
+
Parameters
+
+
reactions
The set of reactions to be screened.
+
species
The list of all species in the network.
+
Y
The molar abundances of the species.
+
T9
The temperature in 10^9 K.
+
rho
The density in g/cm^3.
+
+
+
+
Returns
A vector of screening factors, one for each reaction.
+
Algorithm
+
Low-Temperature Cutoff: If T9 is below a small threshold (1e-9), screening is effectively turned off to prevent numerical instability.
+
Zeta Factor (ζ): A composition-dependent term is calculated: ζ = ∑(Z_i² + Z_i) * Y_i, where Z_i is the charge and Y_i is the molar abundance of species i.
+
Prefactor: A key prefactor is computed: prefactor = 0.188 * sqrt(ρ / T₇³) * sqrt(ζ), where T₇ is the temperature in units of 10^7 K.
+
Screening Term (H_12): For each reaction, the term H_12 is calculated:
+
For a two-body reaction (reactants Z₁ and Z₂): H_12 = prefactor * Z₁ * Z₂.
+
For the triple-alpha reaction (3 * He4): H_12 = 3 * (prefactor * Z_α * Z_α).
+
For one-body reactions (decays), H_12 is 0, so the factor is 1.
+
+
+
Capping: The value of H_12 is capped at 2.0 to prevent excessively large and unphysical screening factors (exp(2) ≈ 7.4).
+
Final Factor: The screening factor for the reaction is exp(H_12).
Calculates weak screening factors using CppAD types.
+
This is the automatic differentiation-compatible version of the method. It allows the derivatives of the screening factors to be computed with respect to plasma conditions.
+
Parameters
+
+
reactions
The set of logical reactions in the network.
+
species
A vector of all atomic species involved in the network.
193// A. Build the basic network from the composition's species with non zero mass fractions.
-
194// B. rebuild a new composition from the reaction set's reactants + products (with the mass fractions from the things that are only products set to 0)
-
195// C. Rebuild the reaction set from the new composition
-
196// D. Cull reactions where all reactants have mass fractions below the culling threshold.
-
197// E. Be careful about maintaining caching through all of this
+
164// This scenario indicates a severe data integrity issue:
+
165// a reactant is part of a reaction but not in the network's species map.
+
166 LOG_ERROR(m_logger, "CRITICAL ERROR: Reactant species '{}' in reaction '{}' not found in network species map during conservation validation.",
208// A. Build the basic network from the composition's species with non zero mass fractions.
+
209// B. rebuild a new composition from the reaction set's reactants + products (with the mass fractions from the things that are only products set to 0)
+
210// C. Rebuild the reaction set from the new composition
+
211// D. Cull reactions where all reactants have mass fractions below the culling threshold.
+
212// E. Be careful about maintaining caching through all of this
+
213
+
214
+
215// This allows for dynamic network modification while retaining caching for networks which are very similar.
464throw std::runtime_error("Cannot record AD tape: No species in the network.");
-
465 }
-
466constsize_t numADInputs = numSpecies + 2; // Note here that by not letting T9 and rho be independent variables, we are constraining the network to a constant temperature and density during each evaluation.
-
467
-
468// --- CppAD Tape Recording ---
-
469// 1. Declare independent variable (adY)
-
470// We also initialize the dummy variable for tape recording (these tell CppAD what the derivative chain looks like).
-
471// Their numeric values are irrelevant except for in so far as they avoid numerical instabilities.
-
472
-
473// Distribute total mass fraction uniformly between species in the dummy variable space
566throw std::runtime_error("Cannot record AD tape: No species in the network.");
+
567 }
+
568constsize_t numADInputs = numSpecies + 2; // Note here that by not letting T9 and rho be independent variables, we are constraining the network to a constant temperature and density during each evaluation.
+
569
+
570// --- CppAD Tape Recording ---
+
571// 1. Declare independent variable (adY)
+
572// We also initialize the dummy variable for tape recording (these tell CppAD what the derivative chain looks like).
+
573// Their numeric values are irrelevant except for in so far as they avoid numerical instabilities.
+
574
+
575// Distribute total mass fraction uniformly between species in the dummy variable space
A factory function to select and create a screening model.
Typedef Documentation
@@ -171,14 +176,24 @@ Functions
+
+
Enumerates the available plasma screening models.
+
This enum provides a set of identifiers for the different screening prescriptions that can be used in the reaction rate calculations.
-
Enumerator
BARE
No screening applied.
+
Enumerator
BARE
No screening applied. The screening factor is always 1.0.
-
WEAK
Weak screening model.
+
WEAK
Weak screening model (Salpeter, 1954).
+
This model is suitable for non-degenerate, non-relativistic plasmas where the electrostatic potential energy between ions is small compared to their thermal kinetic energy. The screening enhancement factor is calculated as exp(H_12).
+
Algorithm
+
A composition-dependent term, ζ = ∑(Z_i^2 + Z_i) * Y_i, is calculated, where Z_i is the charge and Y_i is the molar abundance of each species.
+
A prefactor is computed: prefactor = 0.188 * sqrt(ρ / T₇³) * sqrt(ζ), where ρ is the density and T₇ is the temperature in 10^7 K.
+
For a reaction between two nuclei with charges Z₁ and Z₂, the enhancement term is H_12 = prefactor * Z₁ * Z₂.
+
The final screening factor is exp(H_12). A special calculation is performed for the triple-alpha reaction.
A factory function to select and create a screening model.
+
This function returns a std::unique_ptr to a concrete implementation of the ScreeningModel abstract base class, based on the specified ScreeningType. This allows for easy switching between different screening prescriptions at runtime.
+
Parameters
+
+
type
The ScreeningType enum value specifying which model to create.
+
+
+
+
Returns
A std::unique_ptr<ScreeningModel> holding an instance of the requested screening model.
+
Algorithm The function uses a switch statement to determine which concrete model to instantiate. If the provided type does not match a known case, it defaults to creating a BareScreeningModel to ensure safe behavior.
+
Post-conditions
+
A non-null std::unique_ptr<ScreeningModel> is always returned.
Formats a map of nuclear species timescales into a human-readable string.
Function Documentation
@@ -140,6 +141,46 @@ Functions
+
Formats a map of nuclear species timescales into a human-readable string.
+
This function takes a reaction network engine and the current plasma conditions to calculate the characteristic timescales for each species. It then formats this information into a neatly aligned ASCII table, which is suitable for logging or printing to the console.
+
Parameters
+
+
engine
A constant reference to a DynamicEngine object, used to calculate the species timescales.
+
Y
A vector of the molar abundances (mol/g) for each species.
+
T9
The temperature in units of 10^9 K.
+
rho
The plasma density in g/cm^3.
+
+
+
+
Returns
A std::string containing the formatted table of species and their timescales.
+
Pre-conditions
+
The engine must be in a valid state.
+
The size of the Y vector must be consistent with the number of species expected by the engine.
+
+
Algorithm
+
Calls the getSpeciesTimescales method on the provided engine to get the timescale for each species under the given conditions.
+
Determines the maximum length of the species names to dynamically set the width of the "Species" column for proper alignment.
+
Uses a std::ostringstream to build the output string.
+
Constructs a header for the table with titles "Species" and "Timescale (s)".
+
Iterates through the map of timescales, adding a row to the table for each species.
+
Timescales are formatted in scientific notation with 3 digits of precision.
+
Special handling is included to print "inf" for infinite timescales.
+
The final string, including header and footer lines, is returned.
+
+
Usage
// Assume 'my_engine' is a valid DynamicEngine object and Y, T9, rho are initialized.
diff --git a/src/network/include/gridfire/engine/engine_abstract.h b/src/network/include/gridfire/engine/engine_abstract.h
index a99ad8b2..69abdafe 100644
--- a/src/network/include/gridfire/engine/engine_abstract.h
+++ b/src/network/include/gridfire/engine/engine_abstract.h
@@ -214,10 +214,56 @@ namespace gridfire {
double rho
) const = 0;
+ /**
+ * @brief Update the internal state of the engine.
+ *
+ * @param netIn A struct containing the current network input, such as
+ * temperature, density, and composition.
+ *
+ * This method is intended to be implemented by derived classes to update
+ * their internal state based on the provided network conditions. For example,
+ * an adaptive engine might use this to re-evaluate which reactions and species
+ * are active. For other engines that do not support manually updating, this
+ * method might do nothing.
+ *
+ * @par Usage Example:
+ * @code
+ * NetIn input = { ... };
+ * myEngine.update(input);
+ * @endcode
+ *
+ * @post The internal state of the engine is updated to reflect the new conditions.
+ */
virtual void update(const NetIn& netIn) = 0;
+ /**
+ * @brief Set the electron screening model.
+ *
+ * @param model The type of screening model to use for reaction rate calculations.
+ *
+ * This method allows changing the screening model at runtime. Screening corrections
+ * account for the electrostatic shielding of nuclei by electrons, which affects
+ * reaction rates in dense stellar plasmas.
+ *
+ * @par Usage Example:
+ * @code
+ * myEngine.setScreeningModel(screening::ScreeningType::WEAK);
+ * @endcode
+ *
+ * @post The engine will use the specified screening model for subsequent rate calculations.
+ */
virtual void setScreeningModel(screening::ScreeningType model) = 0;
+ /**
+ * @brief Get the current electron screening model.
+ *
+ * @return The currently active screening model type.
+ *
+ * @par Usage Example:
+ * @code
+ * screening::ScreeningType currentModel = myEngine.getScreeningModel();
+ * @endcode
+ */
[[nodiscard]] virtual screening::ScreeningType getScreeningModel() const = 0;
};
}
\ No newline at end of file
diff --git a/src/network/include/gridfire/engine/views/engine_adaptive.h b/src/network/include/gridfire/engine/views/engine_adaptive.h
index 1cc8cb2f..da0d76ed 100644
--- a/src/network/include/gridfire/engine/views/engine_adaptive.h
+++ b/src/network/include/gridfire/engine/views/engine_adaptive.h
@@ -225,23 +225,59 @@ namespace gridfire {
*/
[[nodiscard]] const DynamicEngine& getBaseEngine() const override { return m_baseEngine; }
+ /**
+ * @brief Sets the screening model for the base engine.
+ *
+ * This method delegates the call to the base engine to set the electron screening model.
+ *
+ * @param model The electron screening model to set.
+ *
+ * @par Usage Example:
+ * @code
+ * AdaptiveEngineView engineView(...);
+ * engineView.setScreeningModel(screening::ScreeningType::WEAK);
+ * @endcode
+ *
+ * @post The screening model of the base engine is updated.
+ */
void setScreeningModel(screening::ScreeningType model) override;
+ /**
+ * @brief Gets the screening model from the base engine.
+ *
+ * This method delegates the call to the base engine to get the screening model.
+ *
+ * @return The current screening model type.
+ *
+ * @par Usage Example:
+ * @code
+ * AdaptiveEngineView engineView(...);
+ * screening::ScreeningType model = engineView.getScreeningModel();
+ * @endcode
+ */
[[nodiscard]] screening::ScreeningType getScreeningModel() const override;
private:
using Config = fourdst::config::Config;
using LogManager = fourdst::logging::LogManager;
+ /** @brief A reference to the singleton Config instance, used for retrieving configuration parameters. */
Config& m_config = Config::getInstance();
+ /** @brief A pointer to the logger instance, used for logging messages. */
quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
+ /** @brief The underlying engine to which this view delegates calculations. */
DynamicEngine& m_baseEngine;
+ /** @brief The set of species that are currently active in the network. */
std::vector m_activeSpecies;
+ /** @brief The set of reactions that are currently active in the network. */
reaction::LogicalReactionSet m_activeReactions;
+ /** @brief A map from the indices of the active species to the indices of the corresponding species in the full network. */
std::vector m_speciesIndexMap;
+ /** @brief A map from the indices of the active reactions to the indices of the corresponding reactions in the full network. */
std::vector m_reactionIndexMap;
+ /** @brief A flag indicating whether the view is stale and needs to be updated. */
bool m_isStale = true;
private:
@@ -322,19 +358,92 @@ namespace gridfire {
*/
void validateState() const;
+ /**
+ * @brief Calculates the molar reaction flow rate for all reactions in the full network.
+ *
+ * This method iterates through all reactions in the base engine's network and calculates
+ * their molar flow rates based on the provided network input conditions (temperature, density,
+ * and composition). It also constructs a vector of molar abundances for all species in the
+ * 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.
+ *
+ * @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.
+ */
std::vector calculateAllReactionFlows(
const NetIn& netIn,
std::vector& out_Y_Full
) const;
+ /**
+ * @brief Finds all species that are reachable from the initial fuel through the reaction network.
+ *
+ * This method performs a connectivity analysis to identify all species that can be produced
+ * starting from the initial fuel species. A species is considered part of the initial fuel if its
+ * mass fraction is above a certain threshold (`ABUNDANCE_FLOOR`).
+ *
+ * @param netIn The current network input, containing the initial composition.
+ * @return An unordered set of all reachable species.
+ *
+ * @par Algorithm:
+ * 1. Initializes a set `reachable` and a queue `to_visit` with the initial fuel species.
+ * 2. Iteratively processes the reaction network until no new species can be reached.
+ * 3. In each pass, it iterates through all reactions in the base engine's network.
+ * 4. If all reactants of a reaction are in the `reachable` set, all products of that reaction are added to the `reachable` set.
+ * 5. The process continues until a full pass over all reactions does not add any new species to the `reachable` set.
+ */
[[nodiscard]] std::unordered_set findReachableSpecies(
const NetIn& netIn
) const;
+ /**
+ * @brief Culls reactions from the network based on their flow rates.
+ *
+ * This method filters the list of all reactions, keeping only those with a flow rate
+ * above an absolute culling threshold. The threshold is calculated by multiplying the
+ * maximum flow rate by a relative culling threshold read from the configuration.
+ *
+ * @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 maxFlow The maximum reaction flow rate in the network.
+ * @return A vector of pointers to the reactions that have been kept after culling.
+ *
+ * @par Algorithm:
+ * 1. Retrieves the `RelativeCullingThreshold` from the configuration.
+ * 2. Calculates the `absoluteCullingThreshold` by multiplying `maxFlow` with the relative threshold.
+ * 3. Iterates through `allFlows`.
+ * 4. A reaction is kept if its `flowRate` is greater than the `absoluteCullingThreshold`.
+ * 5. The pointers to the kept reactions are stored in a vector and returned.
+ */
[[nodiscard]] std::vector cullReactionsByFlow(
const std::vector& allFlows,
const std::unordered_set& reachableSpecies,
const std::vector& Y_full,
double maxFlow
) const;
+ /**
+ * @brief Finalizes the set of active species and reactions.
+ *
+ * This method takes the final list of culled reactions and populates the
+ * `m_activeReactions` and `m_activeSpecies` members. The active species are
+ * determined by collecting all reactants and products from the final reactions.
+ * The active species list is then sorted by mass.
+ *
+ * @param finalReactions A vector of pointers to the reactions to be included in the active set.
+ *
+ * @post
+ * - `m_activeReactions` is cleared and populated with the reactions from `finalReactions`.
+ * - `m_activeSpecies` is cleared and populated with all unique species present in `finalReactions`.
+ * - `m_activeSpecies` is sorted by atomic mass.
+ */
void finalizeActiveSet(
const std::vector& finalReactions
);
diff --git a/src/network/include/gridfire/engine/views/engine_defined.h b/src/network/include/gridfire/engine/views/engine_defined.h
index 90d61ffc..5ec910f4 100644
--- a/src/network/include/gridfire/engine/views/engine_defined.h
+++ b/src/network/include/gridfire/engine/views/engine_defined.h
@@ -13,8 +13,45 @@
#include
namespace gridfire{
+ /**
+ * @class FileDefinedEngineView
+ * @brief An engine view that uses a user-defined reaction network from a file.
+ *
+ * This class implements an EngineView that restricts the reaction network to a specific set of
+ * reactions defined in an external file. It acts as a filter or a view on a larger, more
+ * comprehensive base engine. The file provides a list of reaction identifiers, and this view
+ * will only consider those reactions and the species involved in them.
+ *
+ * This is useful for focusing on a specific sub-network for analysis, debugging, or performance
+ * reasons, without modifying the underlying full network.
+ *
+ * The view maintains mappings between the indices of its active (defined) species and reactions
+ * and the corresponding indices in the full network of the base engine. All calculations are
+ * delegated to the base engine after mapping the inputs from the view's context to the full
+ * network context, and the results are mapped back.
+ *
+ * @implements DynamicEngine
+ * @implements EngineView
+ */
class FileDefinedEngineView final: public DynamicEngine, public EngineView {
public:
+ /**
+ * @brief Constructs a FileDefinedEngineView.
+ *
+ * @param baseEngine The underlying DynamicEngine to which this view delegates calculations.
+ * @param fileName The path to the file that defines the reaction network for this view.
+ * @param parser A reference to a parser object capable of parsing the network file.
+ *
+ * @par Usage Example:
+ * @code
+ * MyParser parser;
+ * DynamicEngine baseEngine(...);
+ * FileDefinedEngineView view(baseEngine, "my_network.net", parser);
+ * @endcode
+ *
+ * @post The view is initialized with the reactions and species from the specified file.
+ * @throws std::runtime_error If a reaction from the file is not found in the base engine.
+ */
explicit FileDefinedEngineView(
DynamicEngine& baseEngine,
const std::string& fileName,
@@ -22,104 +59,246 @@ namespace gridfire{
);
// --- EngineView Interface ---
+ /**
+ * @brief Gets the base engine.
+ * @return A const reference to the base engine.
+ */
const DynamicEngine& getBaseEngine() const override;
// --- Engine Interface ---
+ /**
+ * @brief Gets the list of active species in the network defined by the file.
+ * @return A const reference to the vector of active species.
+ */
const std::vector& getNetworkSpecies() const override;
// --- DynamicEngine Interface ---
+ /**
+ * @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 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
+ * nuclear energy generation rate.
+ *
+ * @throws std::runtime_error If the view is stale (i.e., `update()` has not been called after `setNetworkFile()`).
+ */
StepDerivatives calculateRHSAndEnergy(
const std::vector& Y_defined,
const double T9,
const double rho
) const override;
+ /**
+ * @brief Generates the Jacobian matrix for the active species.
+ *
+ * @param Y_defined A vector of abundances for the active species.
+ * @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_defined,
const double T9,
const double rho
) 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).
+ *
+ * @throws std::runtime_error If the view is stale.
+ * @throws std::out_of_range If an index is out of bounds.
+ */
double getJacobianMatrixEntry(
const int i_defined,
const int j_defined
) const override;
+ /**
+ * @brief Generates the stoichiometry matrix for the active reactions and species.
+ *
+ * @throws std::runtime_error If the view is stale.
+ */
void generateStoichiometryMatrix() override;
+ /**
+ * @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.
+ * @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.
+ */
int getStoichiometryMatrixEntry(
const int speciesIndex_defined,
const int reactionIndex_defined
) 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 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).
+ *
+ * @throws std::runtime_error If the view is stale or if the reaction is not in the active set.
+ */
double calculateMolarReactionFlow(
const reaction::Reaction& reaction,
const std::vector& Y_defined,
const double T9,
const double rho
) const override;
+ /**
+ * @brief Gets the set of active logical reactions in the network.
+ *
+ * @return Reference to the LogicalReactionSet containing all active reactions.
+ *
+ * @throws std::runtime_error If the view is stale.
+ */
const reaction::LogicalReactionSet& getNetworkReactions() const override;
+ /**
+ * @brief Computes timescales for all active species in the network.
+ *
+ * @param Y_defined Vector of current abundances for the active species.
+ * @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).
+ *
+ * @throws std::runtime_error If the view is stale.
+ */
std::unordered_map getSpeciesTimescales(
const std::vector& Y_defined,
const double T9,
const double rho
) const override;
+ /**
+ * @brief Updates the engine view if it is marked as stale.
+ *
+ * This method checks if the view is stale (e.g., after `setNetworkFile` was called).
+ * If it is, it rebuilds the active network from the currently set file.
+ * The `netIn` parameter is not used by this implementation but is required by the interface.
+ *
+ * @param netIn The current network input (unused).
+ *
+ * @post If the view was stale, it is rebuilt and is no longer stale.
+ */
void update(const NetIn &netIn) override;
+ /**
+ * @brief Sets a new network file to define the active reactions.
+ *
+ * @param fileName The path to the new network definition file.
+ *
+ * @par Usage Example:
+ * @code
+ * view.setNetworkFile("another_network.net");
+ * view.update(netIn); // Must be called before using the view again
+ * @endcode
+ *
+ * @post The view is marked as stale. `update()` must be called before further use.
+ */
void setNetworkFile(const std::string& fileName);
+ /**
+ * @brief Sets the screening model for the base engine.
+ *
+ * @param model The screening model to set.
+ */
void setScreeningModel(screening::ScreeningType model) override;
+ /**
+ * @brief Gets the screening model from the base engine.
+ *
+ * @return The current screening model type.
+ */
[[nodiscard]] screening::ScreeningType getScreeningModel() const override;
private:
using Config = fourdst::config::Config;
using LogManager = fourdst::logging::LogManager;
+ /** @brief A reference to the singleton Config instance. */
Config& m_config = Config::getInstance();
+ /** @brief A pointer to the logger instance. */
quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
+ /** @brief The underlying engine to which this view delegates calculations. */
DynamicEngine& m_baseEngine;
- std::string m_fileName; ///< Name of the file defining the reaction set considered by the engine view.
- const io::NetworkFileParser& m_parser; ///< Parser for the network file.
+ ///< Name of the file defining the reaction set considered by the engine view.
+ std::string m_fileName;
+ ///< Parser for the network file.
+ const io::NetworkFileParser& m_parser;
- std::vector m_activeSpecies; ///< Active species in the defined engine.
- reaction::LogicalReactionSet m_activeReactions; ///< Active reactions in the defined engine.
+ ///< Active species in the defined engine.
+ std::vector m_activeSpecies;
+ ///< Active reactions in the defined engine.
+ reaction::LogicalReactionSet m_activeReactions;
- std::vector m_speciesIndexMap; ///< Maps indices of active species to indices in the full network.
- std::vector m_reactionIndexMap; ///< Maps indices of active reactions to indices in the full network.
+ ///< Maps indices of active species to indices in the full network.
+ std::vector m_speciesIndexMap;
+ ///< Maps indices of active reactions to indices in the full network.
+ std::vector m_reactionIndexMap;
+ /** @brief A flag indicating whether the view is stale and needs to be updated. */
bool m_isStale = true;
private:
+ /**
+ * @brief Builds the active species and reaction sets from a file.
+ *
+ * This method uses the provided parser to read reaction names from the given file.
+ * It then finds these reactions in the base engine's full network and populates
+ * the `m_activeReactions` and `m_activeSpecies` members. Finally, it constructs
+ * the index maps for the active sets.
+ *
+ * @param fileName The path to the network definition file.
+ *
+ * @post
+ * - `m_activeReactions` and `m_activeSpecies` are populated.
+ * - `m_speciesIndexMap` and `m_reactionIndexMap` are constructed.
+ * - `m_isStale` is set to false.
+ *
+ * @throws std::runtime_error If a reaction from the file is not found in the base engine.
+ */
void buildFromFile(const std::string& fileName);
/**
* @brief Constructs the species index map.
*
- * @return A vector mapping culled species indices to full species indices.
+ * @return A vector mapping defined species indices to full species indices.
*
* This method creates a map from the indices of the active species to the indices of the
* corresponding species in the full network.
*
- * @see AdaptiveEngineView::update()
+ * @throws std::runtime_error If an active species is not found in the base engine's species list.
*/
std::vector constructSpeciesIndexMap() const;
/**
* @brief Constructs the reaction index map.
*
- * @return A vector mapping culled reaction indices to full reaction indices.
+ * @return A vector mapping defined reaction indices to full reaction indices.
*
* This method creates a map from the indices of the active reactions to the indices of the
* corresponding reactions in the full network.
*
- * @see AdaptiveEngineView::update()
+ * @throws std::runtime_error If an active reaction is not found in the base engine's reaction list.
*/
std::vector constructReactionIndexMap() const;
/**
* @brief Maps a vector of culled abundances to a vector of full abundances.
*
- * @param culled A vector of abundances for the active species.
+ * @param defined A vector of abundances for the active species.
* @return A vector of abundances for the full network, with the abundances of the active
- * species copied from the culled vector.
+ * species copied from the defined vector.
*/
- std::vector mapViewToFull(const std::vector& culled) const;
+ std::vector mapViewToFull(const std::vector& defined) const;
/**
* @brief Maps a vector of full abundances to a vector of culled abundances.
@@ -133,23 +312,28 @@ namespace gridfire{
/**
* @brief Maps a culled species index to a full species index.
*
- * @param culledSpeciesIndex The index of the species in the culled species list.
+ * @param definedSpeciesIndex The index of the species in the defined species list.
* @return The index of the corresponding species in the full network.
*
- * @throws std::out_of_range If the culled index is out of bounds for the species index map.
+ * @throws std::out_of_range If the defined index is out of bounds for the species index map.
*/
- size_t mapViewToFullSpeciesIndex(size_t culledSpeciesIndex) const;
+ size_t mapViewToFullSpeciesIndex(size_t definedSpeciesIndex) const;
/**
* @brief Maps a culled reaction index to a full reaction index.
*
- * @param culledReactionIndex The index of the reaction in the culled reaction list.
+ * @param definedReactionIndex The index of the reaction in the defined reaction list.
* @return The index of the corresponding reaction in the full network.
*
- * @throws std::out_of_range If the culled index is out of bounds for the reaction index map.
+ * @throws std::out_of_range If the defined index is out of bounds for the reaction index map.
*/
- size_t mapViewToFullReactionIndex(size_t culledReactionIndex) const;
+ size_t mapViewToFullReactionIndex(size_t definedReactionIndex) const;
+ /**
+ * @brief Validates that the FileDefinedEngineView is not stale.
+ *
+ * @throws std::runtime_error If the view is stale (i.e., `update()` has not been called after the view was made stale).
+ */
void validateNetworkState() const;
};
}
\ No newline at end of file
diff --git a/src/network/include/gridfire/io/network_file.h b/src/network/include/gridfire/io/network_file.h
index 4ed81f00..d742ec7b 100644
--- a/src/network/include/gridfire/io/network_file.h
+++ b/src/network/include/gridfire/io/network_file.h
@@ -10,21 +10,111 @@
namespace gridfire::io {
+ /**
+ * @struct ParsedNetworkData
+ * @brief Holds the data parsed from a network file.
+ *
+ * This struct is used to return the results of parsing a reaction network
+ * file. It contains the list of reaction names that define the network.
+ */
struct ParsedNetworkData {
+ /**
+ * @brief A vector of reaction names in their PEN-style format.
+ *
+ * Projectile, Ejectile style names p(p,e+)d is a common format for representing
+ * nuclear reactions as strings.
+ */
std::vector reactionPENames;
};
+ /**
+ * @class NetworkFileParser
+ * @brief An abstract base class for network file parsers.
+ *
+ * This class defines the interface for parsing files that contain
+ * reaction network definitions. Derived classes must implement the `parse`
+ * method to handle specific file formats.
+ */
class NetworkFileParser {
public:
+ /**
+ * @brief Virtual destructor for the base class.
+ */
virtual ~NetworkFileParser() = default;
+ /**
+ * @brief Parses a network file and returns the parsed data.
+ *
+ * This is a pure virtual function that must be implemented by derived
+ * classes. It takes a filename as input and returns a `ParsedNetworkData`
+ * struct containing the information extracted from the file.
+ *
+ * @param filename The path to the network file to parse.
+ * @return A `ParsedNetworkData` struct containing the parsed reaction data.
+ *
+ * @throws std::runtime_error If the file cannot be opened or a parsing
+ * error occurs.
+ *
+ * @b Usage
+ * @code
+ * std::unique_ptr parser = std::make_unique();
+ * try {
+ * ParsedNetworkData data = parser->parse("my_reactions.txt");
+ * for (const auto& reaction_name : data.reactionPENames) {
+ * // ... process reaction name
+ * }
+ * } catch (const std::runtime_error& e) {
+ * // ... handle error
+ * }
+ * @endcode
+ */
[[nodiscard]] virtual ParsedNetworkData parse(const std::string& filename) const = 0;
};
+ /**
+ * @class SimpleReactionListFileParser
+ * @brief A parser for simple text files containing a list of reactions.
+ *
+ * This parser reads a file where each line contains a single reaction name.
+ * It supports comments (lines starting with '#') and ignores empty lines.
+ *
+ * @implements NetworkFileParser
+ */
class SimpleReactionListFileParser final : public NetworkFileParser {
public:
+ /**
+ * @brief Constructs a SimpleReactionListFileParser.
+ *
+ * @post The parser is initialized and ready to parse files.
+ */
explicit SimpleReactionListFileParser();
+ /**
+ * @brief Parses a simple reaction list file.
+ *
+ * This method reads the specified file line by line. It trims whitespace
+ * from each line, ignores lines that are empty or start with a '#'
+ * comment character, and stores the remaining lines as reaction names.
+ *
+ * @param filename The path to the simple reaction list file.
+ * @return A `ParsedNetworkData` struct containing the list of reaction names.
+ *
+ * @throws std::runtime_error If the file cannot be opened for reading.
+ *
+ * @b Algorithm
+ * 1. Opens the specified file.
+ * 2. Reads the file line by line.
+ * 3. For each line, it removes any trailing comments (starting with '#').
+ * 4. Trims leading and trailing whitespace.
+ * 5. If the line is not empty, it is added to the list of reaction names.
+ * 6. Returns the populated `ParsedNetworkData` struct.
+ *
+ * @b Usage
+ * @code
+ * SimpleReactionListFileParser parser;
+ * ParsedNetworkData data = parser.parse("reactions.txt");
+ * @endcode
+ */
ParsedNetworkData parse(const std::string& filename) const override;
private:
using Config = fourdst::config::Config;
@@ -33,9 +123,38 @@ namespace gridfire::io {
quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
};
+ /**
+ * @class MESANetworkFileParser
+ * @brief A parser for MESA-format network files.
+ *
+ * This class is designed to parse reaction network files that follow the
+ * format used by the MESA stellar evolution code.
+ *
+ * @implements NetworkFileParser
+ */
class MESANetworkFileParser final : public NetworkFileParser {
public:
+ /**
+ * @brief Constructs a MESANetworkFileParser.
+ *
+ * @param filename The path to the MESA network file. This may be used
+ * to pre-configure the parser.
+ *
+ * @post The parser is initialized with the context of the given file.
+ */
explicit MESANetworkFileParser(const std::string& filename);
+ /**
+ * @brief Parses a MESA-format network file.
+ *
+ * This method will read and interpret the structure of a MESA network
+ * file to extract the list of reactions.
+ *
+ * @param filename The path to the MESA network file.
+ * @return A `ParsedNetworkData` struct containing the list of reaction names.
+ *
+ * @throws std::runtime_error If the file cannot be opened or if it
+ * contains formatting errors.
+ */
ParsedNetworkData parse(const std::string& filename) const override;
private:
using Config = fourdst::config::Config;
diff --git a/src/network/include/gridfire/screening/screening_abstract.h b/src/network/include/gridfire/screening/screening_abstract.h
index 0529fcda..0a03ce89 100644
--- a/src/network/include/gridfire/screening/screening_abstract.h
+++ b/src/network/include/gridfire/screening/screening_abstract.h
@@ -9,11 +9,67 @@
#include
namespace gridfire::screening {
+ /**
+ * @class ScreeningModel
+ * @brief An abstract base class for plasma screening models.
+ *
+ * This class defines the interface for models that calculate the enhancement
+ * factor for nuclear reaction rates due to the electrostatic screening of
+ * interacting nuclei by the surrounding plasma. Concrete implementations of
+ * this class will provide specific screening prescriptions (e.g., WEAK,
+ * BARE, STRONG, etc.).
+ *
+ * The interface provides methods for calculating screening factors for both
+ * standard double-precision inputs and for CppAD's automatic differentiation
+ * types, allowing the screening contributions to be included in Jacobian
+ * calculations.
+ */
class ScreeningModel {
public:
+ /// @brief Alias for CppAD Automatic Differentiation type for double precision.
using ADDouble = CppAD::AD;
+ /**
+ * @brief Virtual destructor.
+ *
+ * Ensures that derived class destructors are called correctly.
+ */
virtual ~ScreeningModel() = default;
+ /**
+ * @brief Calculates screening factors for a set of reactions.
+ *
+ * This is a pure virtual function that must be implemented by derived
+ * classes. It computes the screening enhancement factor for each reaction
+ * in the provided set based on the given plasma conditions.
+ *
+ * @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 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
+ * in the `reactions` set, in the same order.
+ *
+ * @b Pre-conditions
+ * - The size of the `Y` vector must match the size of the `species` vector.
+ * - `T9` and `rho` must be positive.
+ *
+ * @b Post-conditions
+ * - The returned vector will have the same size as the `reactions` set.
+ * - Each element in the returned vector will be >= 1.0.
+ *
+ * @b Usage
+ * @code
+ * // Assume 'model' is a std::unique_ptr to a concrete implementation
+ * // and other parameters (reactions, species, Y, T9, rho) are initialized.
+ * std::vector screening_factors = model->calculateScreeningFactors(
+ * reactions, species, Y, T9, rho
+ * );
+ * for (size_t i = 0; i < reactions.size(); ++i) {
+ * // ... use screening_factors[i] ...
+ * }
+ * @endcode
+ */
virtual std::vector calculateScreeningFactors(
const reaction::LogicalReactionSet& reactions,
const std::vector& species,
@@ -22,6 +78,25 @@ namespace gridfire::screening {
const double rho
) const = 0;
+ /**
+ * @brief Calculates screening factors using CppAD types for automatic differentiation.
+ *
+ * This is a pure virtual function that provides an overload of
+ * `calculateScreeningFactors` for use with CppAD. It allows the derivatives
+ * of the screening factors with respect to abundances, temperature, and
+ * density to be computed automatically.
+ *
+ * @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 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.
+ *
+ * @b Note
+ * This method is essential for including the effects of screening in the
+ * Jacobian matrix of the reaction network.
+ */
virtual std::vector calculateScreeningFactors(
const reaction::LogicalReactionSet& reactions,
const std::vector& species,
diff --git a/src/network/include/gridfire/screening/screening_bare.h b/src/network/include/gridfire/screening/screening_bare.h
index 52bff705..f1a66682 100644
--- a/src/network/include/gridfire/screening/screening_bare.h
+++ b/src/network/include/gridfire/screening/screening_bare.h
@@ -6,9 +6,51 @@
#include "cppad/cppad.hpp"
namespace gridfire::screening {
+ /**
+ * @class BareScreeningModel
+ * @brief A screening model that applies no screening effect.
+ *
+ * This class implements the `ScreeningModel` interface but returns a
+ * screening factor of 1.0 for all reactions, regardless of the plasma
+ * conditions. It represents the case of bare, unscreened nuclei and serves
+ * as a baseline or can be used when screening effects are negligible or
+ * intentionally ignored.
+ *
+ * @implements ScreeningModel
+ */
class BareScreeningModel final : public ScreeningModel {
+ /// @brief Alias for CppAD Automatic Differentiation type for double precision.
using ADDouble = CppAD::AD;
public:
+ /**
+ * @brief Calculates screening factors, which are always 1.0.
+ *
+ * This implementation returns a vector of screening factors where every
+ * element is 1.0, effectively applying no screening correction to the
+ * reaction rates.
+ *
+ * @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 T9 The temperature (unused).
+ * @param rho The plasma density (unused).
+ * @return A vector of doubles, with each element being 1.0, of the same
+ * size as the `reactions` set.
+ *
+ * @b Algorithm
+ * The function simply creates and returns a `std::vector` of the
+ * same size as the input `reactions` set, with all elements initialized to 1.0.
+ *
+ * @b Usage
+ * @code
+ * BareScreeningModel bare_model;
+ * // ... (initialize reactions, species, Y, T9, rho)
+ * std::vector factors = bare_model.calculateScreeningFactors(
+ * reactions, species, Y, T9, rho
+ * );
+ * // 'factors' will contain [1.0, 1.0, ...]
+ * @endcode
+ */
[[nodiscard]] std::vector calculateScreeningFactors(
const reaction::LogicalReactionSet& reactions,
const std::vector& species,
@@ -17,6 +59,21 @@ namespace gridfire::screening {
const double rho
) const override;
+ /**
+ * @brief Calculates screening factors for AD types, which are always 1.0.
+ *
+ * This implementation returns a vector of AD-typed screening factors where
+ * every element is 1.0. This is the automatic differentiation-compatible
+ * version.
+ *
+ * @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 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
+ * size as the `reactions` set.
+ */
[[nodiscard]] std::vector calculateScreeningFactors(
const reaction::LogicalReactionSet& reactions,
const std::vector& species,
@@ -25,6 +82,21 @@ namespace gridfire::screening {
const ADDouble rho
) const override;
private:
+ /**
+ * @brief Template implementation for calculating screening factors.
+ *
+ * This private helper function contains the core logic for both the `double`
+ * and `ADDouble` versions of `calculateScreeningFactors`. It is templated
+ * to handle both numeric types seamlessly.
+ *
+ * @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 T9 The temperature (unused).
+ * @param rho The density (unused).
+ * @return A vector of type `T` with all elements initialized to 1.0.
+ */
template
[[nodiscard]] std::vector calculateFactors_impl(
const reaction::LogicalReactionSet& reactions,
@@ -35,6 +107,21 @@ namespace gridfire::screening {
) const;
};
+ /**
+ * @brief Template implementation for the bare screening model.
+ *
+ * This function provides the actual implementation for `calculateFactors_impl`.
+ * It creates a vector of the appropriate numeric type (`T`) and size, and
+ * initializes all its elements to 1.0, representing no 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 T9 Unused parameter.
+ * @param rho Unused parameter.
+ * @return A `std::vector` of the same size as `reactions`, with all elements set to 1.0.
+ */
template
std::vector BareScreeningModel::calculateFactors_impl(
const reaction::LogicalReactionSet &reactions,
diff --git a/src/network/include/gridfire/screening/screening_types.h b/src/network/include/gridfire/screening/screening_types.h
index 77c2168e..1e850790 100644
--- a/src/network/include/gridfire/screening/screening_types.h
+++ b/src/network/include/gridfire/screening/screening_types.h
@@ -5,10 +5,66 @@
#include
namespace gridfire::screening {
+ /**
+ * @enum ScreeningType
+ * @brief Enumerates the available plasma screening models.
+ *
+ * This enum provides a set of identifiers for the different screening
+ * prescriptions that can be used in the reaction rate calculations.
+ */
enum class ScreeningType {
- BARE, ///< No screening applied
- WEAK, ///< Weak screening model
+ BARE, ///< No screening applied. The screening factor is always 1.0.
+ /**
+ * @brief Weak screening model (Salpeter, 1954).
+ *
+ * This model is suitable for non-degenerate, non-relativistic plasmas
+ * where the electrostatic potential energy between ions is small compared
+ * to their thermal kinetic energy. The screening enhancement factor is
+ * calculated as `exp(H_12)`.
+ *
+ * @b Algorithm
+ * 1. A composition-dependent term, `ζ = ∑(Z_i^2 + Z_i) * Y_i`, is calculated,
+ * where Z_i is the charge and Y_i is the molar abundance of each species.
+ * 2. A prefactor is computed: `prefactor = 0.188 * sqrt(ρ / T₇³) * sqrt(ζ)`,
+ * where ρ is the density and T₇ is the temperature in 10^7 K.
+ * 3. For a reaction between two nuclei with charges Z₁ and Z₂, the enhancement
+ * term is `H_12 = prefactor * Z₁ * Z₂`.
+ * 4. The final screening factor is `exp(H_12)`.
+ * A special calculation is performed for the triple-alpha reaction.
+ */
+ WEAK,
};
+ /**
+ * @brief A factory function to select and create a screening model.
+ *
+ * This function returns a `std::unique_ptr` to a concrete implementation of
+ * the `ScreeningModel` abstract base class, based on the specified `ScreeningType`.
+ * This allows for easy switching between different screening prescriptions at runtime.
+ *
+ * @param type The `ScreeningType` enum value specifying which model to create.
+ * @return A `std::unique_ptr` holding an instance of the
+ * requested screening model.
+ *
+ * @b Algorithm
+ * The function uses a `switch` statement to determine which concrete model to
+ * instantiate. If the provided `type` does not match a known case, it defaults
+ * to creating a `BareScreeningModel` to ensure safe behavior.
+ *
+ * @b Post-conditions
+ * - A non-null `std::unique_ptr` is always returned.
+ *
+ * @b Usage
+ * @code
+ * // Select the weak screening model
+ * auto screening_model = gridfire::screening::selectScreeningModel(gridfire::screening::ScreeningType::WEAK);
+ *
+ * // Use the model to calculate screening factors
+ * // (assuming other parameters are initialized)
+ * std::vector factors = screening_model->calculateScreeningFactors(
+ * reactions, species, Y, T9, rho
+ * );
+ * @endcode
+ */
std::unique_ptr selectScreeningModel(ScreeningType type);
}
\ No newline at end of file
diff --git a/src/network/include/gridfire/screening/screening_weak.h b/src/network/include/gridfire/screening/screening_weak.h
index a1c1286a..94c82e4f 100644
--- a/src/network/include/gridfire/screening/screening_weak.h
+++ b/src/network/include/gridfire/screening/screening_weak.h
@@ -10,8 +10,43 @@
#include "cppad/cppad.hpp"
namespace gridfire::screening {
+ /**
+ * @class WeakScreeningModel
+ * @brief Implements the weak screening model based on the Debye-Hückel approximation.
+ *
+ * This class provides a concrete implementation of the `ScreeningModel`
+ * interface for the weak screening regime, following the formulation of
+ * Salpeter (1954). This approach applies the Debye-Hückel theory to model the
+ * electrostatic shielding of nuclei in a plasma. It is applicable to
+ * non-degenerate, non-relativistic plasmas where thermal energy dominates
+ * the electrostatic potential energy.
+ *
+ * @implements ScreeningModel
+ */
class WeakScreeningModel final : public ScreeningModel {
public:
+ /**
+ * @brief Calculates weak screening factors for a set of reactions.
+ *
+ * This method computes the screening enhancement factor for each reaction
+ * based on the Salpeter (1954) formula.
+ *
+ * @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 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.
+ *
+ * @b Usage
+ * @code
+ * WeakScreeningModel weak_model;
+ * // ... (initialize reactions, species, Y, T9, rho)
+ * std::vector factors = weak_model.calculateScreeningFactors(
+ * reactions, species, Y, T9, rho
+ * );
+ * @endcode
+ */
[[nodiscard]] std::vector calculateScreeningFactors(
const reaction::LogicalReactionSet& reactions,
const std::vector& species,
@@ -20,6 +55,20 @@ namespace gridfire::screening {
const double rho
) const override;
+ /**
+ * @brief Calculates weak screening factors using CppAD types.
+ *
+ * This is the automatic differentiation-compatible version of the method.
+ * It allows the derivatives of the screening factors to be computed with
+ * respect to plasma conditions.
+ *
+ * @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 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.
+ */
[[nodiscard]] std::vector> calculateScreeningFactors(
const reaction::LogicalReactionSet& reactions,
const std::vector& species,
@@ -28,9 +77,25 @@ namespace gridfire::screening {
const CppAD::AD rho
) const override;
private:
+ /// @brief Logger instance for recording trace and debug information.
quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
private:
+ /**
+ * @brief Template implementation for calculating weak screening factors.
+ *
+ * This private helper function contains the core logic for calculating
+ * weak screening factors. It is templated to handle both `double` and
+ * `CppAD::AD` numeric types, avoiding code duplication.
+ *
+ * @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 T9 The temperature in 10^9 K.
+ * @param rho The density in g/cm^3.
+ * @return A vector of screening factors of type `T`.
+ */
template
[[nodiscard]] std::vector calculateFactors_impl(
const reaction::LogicalReactionSet& reactions,
@@ -41,6 +106,37 @@ namespace gridfire::screening {
) const;
};
+ /**
+ * @brief Core implementation of the weak screening calculation (Debye-Hückel model).
+ *
+ * This function calculates the screening factor `exp(H_12)` for each reaction,
+ * based on the Debye-Hückel approximation as formulated by Salpeter (1954).
+ *
+ * @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 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.
+ *
+ * @b Algorithm
+ * 1. **Low-Temperature Cutoff**: If T9 is below a small threshold (1e-9),
+ * screening is effectively turned off to prevent numerical instability.
+ * 2. **Zeta Factor (ζ)**: A composition-dependent term is calculated:
+ * `ζ = ∑(Z_i² + Z_i) * Y_i`, where Z_i is the charge and Y_i is the
+ * molar abundance of species i.
+ * 3. **Prefactor**: A key prefactor is computed:
+ * `prefactor = 0.188 * sqrt(ρ / T₇³) * sqrt(ζ)`,
+ * where T₇ is the temperature in units of 10^7 K.
+ * 4. **Screening Term (H_12)**: For each reaction, the term H_12 is calculated:
+ * - For a two-body reaction (reactants Z₁ and Z₂): `H_12 = prefactor * Z₁ * Z₂`.
+ * - For the triple-alpha reaction (3 * He4): `H_12 = 3 * (prefactor * Z_α * Z_α)`.
+ * - For one-body reactions (decays), H_12 is 0, so the factor is 1.
+ * 5. **Capping**: The value of H_12 is capped at 2.0 to prevent excessively large
+ * and unphysical screening factors (exp(2) ≈ 7.4).
+ * 6. **Final Factor**: The screening factor for the reaction is `exp(H_12)`.
+ */
template
std::vector WeakScreeningModel::calculateFactors_impl(
const reaction::LogicalReactionSet& reactions,
diff --git a/src/network/include/gridfire/utils/logging.h b/src/network/include/gridfire/utils/logging.h
index e3e6f96f..c1b8d0b4 100644
--- a/src/network/include/gridfire/utils/logging.h
+++ b/src/network/include/gridfire/utils/logging.h
@@ -6,6 +6,56 @@
#include
namespace gridfire::utils {
+ /**
+ * @brief Formats a map of nuclear species timescales into a human-readable string.
+ *
+ * This function takes a reaction network engine and the current plasma
+ * conditions to calculate the characteristic timescales for each species.
+ * It then formats this information into a neatly aligned ASCII table, which
+ * is suitable for logging or printing to the console.
+ *
+ * @param engine A constant reference to a `DynamicEngine` object, used to
+ * calculate the species timescales.
+ * @param Y A vector of the 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 std::string containing the formatted table of species and their
+ * timescales.
+ *
+ * @b Pre-conditions
+ * - The `engine` must be in a valid state.
+ * - The size of the `Y` vector must be consistent with the number of species
+ * expected by the `engine`.
+ *
+ * @b Algorithm
+ * 1. Calls the `getSpeciesTimescales` method on the provided `engine` to get
+ * the timescale for each species under the given conditions.
+ * 2. Determines the maximum length of the species names to dynamically set the
+ * width of the "Species" column for proper alignment.
+ * 3. Uses a `std::ostringstream` to build the output string.
+ * 4. Constructs a header for the table with titles "Species" and "Timescale (s)".
+ * 5. Iterates through the map of timescales, adding a row to the table for
+ * each species.
+ * 6. Timescales are formatted in scientific notation with 3 digits of precision.
+ * 7. Special handling is included to print "inf" for infinite timescales.
+ * 8. The final string, including header and footer lines, is returned.
+ *
+ * @b Usage
+ * @code
+ * // Assume 'my_engine' is a valid DynamicEngine object and Y, T9, rho are initialized.
+ * std::string log_output = gridfire::utils::formatNuclearTimescaleLogString(my_engine, Y, T9, rho);
+ * std::cout << log_output;
+ *
+ * // Example Output:
+ * // == Timescales (s) ==
+ * // Species Timescale (s)
+ * // ==========================
+ * // h1 1.234e+05
+ * // he4 inf
+ * // c12 8.765e-02
+ * // ==========================
+ * @endcode
+ */
std::string formatNuclearTimescaleLogString(
const DynamicEngine& engine,
const std::vector& Y,
Public Member Functions inherited from gridfire::Engine
