tboudreaux/GridFire
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

feat(GridFire): added weak electron screening

Browse Source
This commit is contained in:
Emily Boudreaux
2025-07-01 11:40:03 -04:00
parent 40b28477ed
commit 0c16e81e98
26 changed files with 1408 additions and 444 deletions
Download Patch File Download Diff File Expand all files Collapse all files

23
src/network/include/gridfire/engine/engine_abstract.h
View File

@@ -1,6 +1,9 @@
#pragma once
#include "gridfire/reaction/reaction.h"
#include "gridfire/network.h"
#include "gridfire/screening/screening_abstract.h"
#include "gridfire/screening/screening_types.h"
#include <vector>
#include <unordered_map>
@@ -80,7 +83,7 @@ namespace gridfire {
* @brief Get the list of species in the network.
* @return Vector of Species objects representing all network species.
*/
virtual const std::vector<fourdst::atomic::Species>& getNetworkSpecies() const = 0;
[[nodiscard]] virtual const std::vector<fourdst::atomic::Species>& getNetworkSpecies() const = 0;
/**
* @brief Calculate the right-hand side (dY/dt) and energy generation.
@@ -94,7 +97,7 @@ namespace gridfire {
* time derivatives of all species and the specific nuclear energy generation
* rate for the current state.
*/
virtual StepDerivatives<double> calculateRHSAndEnergy(
[[nodiscard]] virtual StepDerivatives<double> calculateRHSAndEnergy(
const std::vector<double>& Y,
double T9,
double rho
@@ -141,7 +144,7 @@ namespace gridfire {
*
* The Jacobian must have been generated by generateJacobianMatrix() before calling this.
*/
virtual double getJacobianMatrixEntry(
[[nodiscard]] virtual double getJacobianMatrixEntry(
int i,
int j
) const = 0;
@@ -163,7 +166,7 @@ namespace gridfire {
*
* The stoichiometry matrix must have been generated by generateStoichiometryMatrix().
*/
virtual int getStoichiometryMatrixEntry(
[[nodiscard]] virtual int getStoichiometryMatrixEntry(
int speciesIndex,
int reactionIndex
) const = 0;
@@ -180,7 +183,7 @@ namespace gridfire {
* This method computes the net rate at which the given reaction proceeds
* under the current state.
*/
virtual double calculateMolarReactionFlow(
[[nodiscard]] virtual double calculateMolarReactionFlow(
const reaction::Reaction& reaction,
const std::vector<double>& Y,
double T9,
@@ -192,7 +195,7 @@ namespace gridfire {
*
* @return Reference to the LogicalReactionSet containing all reactions.
*/
virtual const reaction::LogicalReactionSet& getNetworkReactions() const = 0;
[[nodiscard]] virtual const reaction::LogicalReactionSet& getNetworkReactions() const = 0;
/**
* @brief Compute timescales for all species in the network.
@@ -205,10 +208,16 @@ namespace gridfire {
* This method estimates the timescale for abundance change of each species,
* which can be used for timestep control, diagnostics, and reaction network culling.
*/
virtual std::unordered_map<fourdst::atomic::Species, double> getSpeciesTimescales(
[[nodiscard]] virtual std::unordered_map<fourdst::atomic::Species, double> getSpeciesTimescales(
const std::vector<double>& Y,
double T9,
double rho
) const = 0;
virtual void update(const NetIn& netIn) = 0;
virtual void setScreeningModel(screening::ScreeningType model) = 0;
[[nodiscard]] virtual screening::ScreeningType getScreeningModel() const = 0;
};
}

31
src/network/include/gridfire/engine/engine_graph.h
View File

@@ -8,10 +8,13 @@
#include "gridfire/network.h"
#include "gridfire/reaction/reaction.h"
#include "gridfire/engine/engine_abstract.h"
#include "gridfire/screening/screening_abstract.h"
#include "gridfire/screening/screening_types.h"
#include <string>
#include <unordered_map>
#include <vector>
#include <memory>
#include <boost/numeric/ublas/matrix_sparse.hpp>
@@ -120,7 +123,7 @@ namespace gridfire {
*
* @see StepDerivatives
*/
StepDerivatives<double> calculateRHSAndEnergy(
[[nodiscard]] StepDerivatives<double> calculateRHSAndEnergy(
const std::vector<double>& Y,
const double T9,
const double rho
@@ -165,7 +168,7 @@ namespace gridfire {
* This method computes the net rate at which the given reaction proceeds
* under the current state.
*/
double calculateMolarReactionFlow(
[[nodiscard]] double calculateMolarReactionFlow(
const reaction::Reaction& reaction,
const std::vector<double>&Y,
const double T9,
@@ -243,6 +246,8 @@ namespace gridfire {
double rho
) const override;
void update(const NetIn& netIn) override;
/**
* @brief Checks if a given species is involved in the network.
*
@@ -293,6 +298,10 @@ namespace gridfire {
const std::string& filename
) const;
void setScreeningModel(screening::ScreeningType) override;
[[nodiscard]] screening::ScreeningType getScreeningModel() const override;
private:
reaction::LogicalReactionSet m_reactions; ///< Set of REACLIB reactions in the network.
@@ -307,6 +316,9 @@ namespace gridfire {
CppAD::ADFun<double> m_rhsADFun; ///< CppAD function for the right-hand side of the ODE.
screening::ScreeningType m_screeningType = screening::ScreeningType::BARE; ///< Screening type for the reaction network. Default to no screening.
std::unique_ptr<screening::ScreeningModel> m_screeningModel = screening::selectScreeningModel(m_screeningType);
Config& m_config = Config::getInstance();
Constants& m_constants = Constants::getInstance(); ///< Access to physical constants.
quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
@@ -427,7 +439,7 @@ namespace gridfire {
* specific nuclear energy generation rate for the current state.
*/
template<IsArithmeticOrAD T>
StepDerivatives<T> calculateAllDerivatives(
[[nodiscard]] StepDerivatives<T> calculateAllDerivatives(
const std::vector<T> &Y_in,
T T9,
T rho
@@ -445,7 +457,7 @@ namespace gridfire {
* specific nuclear energy generation rate for the current state using
* double precision arithmetic.
*/
StepDerivatives<double> calculateAllDerivatives(
[[nodiscard]] StepDerivatives<double> calculateAllDerivatives(
const std::vector<double>& Y_in,
const double T9,
const double rho
@@ -463,7 +475,7 @@ namespace gridfire {
* specific nuclear energy generation rate for the current state using
* automatic differentiation.
*/
StepDerivatives<ADDouble> calculateAllDerivatives(
[[nodiscard]] StepDerivatives<ADDouble> calculateAllDerivatives(
const std::vector<ADDouble>& Y_in,
const ADDouble &T9,
const ADDouble &rho
@@ -474,6 +486,13 @@ namespace gridfire {
template<IsArithmeticOrAD T>
StepDerivatives<T> GraphEngine::calculateAllDerivatives(
const std::vector<T> &Y_in, T T9, T rho) const {
std::vector<T> screeningFactors = m_screeningModel->calculateScreeningFactors(
m_reactions,
m_networkSpecies,
Y_in,
T9,
rho
);
// --- Setup output derivatives structure ---
StepDerivatives<T> result;
@@ -512,7 +531,7 @@ namespace gridfire {
const auto& reaction = m_reactions[reactionIndex];
// 1. Calculate reaction rate
const T molarReactionFlow = calculateMolarReactionFlow<T>(reaction, Y, T9, rho);
const T molarReactionFlow = screeningFactors[reactionIndex] * calculateMolarReactionFlow<T>(reaction, Y, T9, rho);
// 2. Use the rate to update all relevant species derivatives (dY/dt)
for (size_t speciesIndex = 0; speciesIndex < m_networkSpecies.size(); ++speciesIndex) {

42
src/network/include/gridfire/engine/views/engine_adaptive.h
View File

@@ -1,6 +1,8 @@
#pragma once
#include "gridfire/engine/engine_abstract.h"
#include "engine_view_abstract.h"
#include "gridfire/engine/views/engine_view_abstract.h"
#include "gridfire/screening/screening_abstract.h"
#include "gridfire/screening/screening_types.h"
#include "gridfire/network.h"
#include "fourdst/composition/atomicSpecies.h"
@@ -71,13 +73,13 @@ namespace gridfire {
* @see AdaptiveEngineView::constructSpeciesIndexMap()
* @see AdaptiveEngineView::constructReactionIndexMap()
*/
void update(const NetIn& netIn);
void update(const NetIn& netIn) override;
/**
* @brief Gets the list of active species in the network.
* @return A const reference to the vector of active species.
*/
const std::vector<fourdst::atomic::Species>& getNetworkSpecies() const override;
[[nodiscard]] const std::vector<fourdst::atomic::Species>& getNetworkSpecies() const override;
/**
* @brief Calculates the right-hand side (dY/dt) and energy generation for the active species.
@@ -95,7 +97,7 @@ namespace gridfire {
* @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called).
* @see AdaptiveEngineView::update()
*/
StepDerivatives<double> calculateRHSAndEnergy(
[[nodiscard]] StepDerivatives<double> calculateRHSAndEnergy(
const std::vector<double> &Y_culled,
const double T9,
const double rho
@@ -134,7 +136,7 @@ namespace gridfire {
* @throws std::out_of_range If the culled index is out of bounds for the species index map.
* @see AdaptiveEngineView::update()
*/
double getJacobianMatrixEntry(
[[nodiscard]] double getJacobianMatrixEntry(
const int i_culled,
const int j_culled
) const override;
@@ -164,7 +166,7 @@ namespace gridfire {
* @throws std::out_of_range If the culled index is out of bounds for the species or reaction index map.
* @see AdaptiveEngineView::update()
*/
int getStoichiometryMatrixEntry(
[[nodiscard]] int getStoichiometryMatrixEntry(
const int speciesIndex_culled,
const int reactionIndex_culled
) const override;
@@ -184,7 +186,7 @@ namespace gridfire {
* @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called).
* @throws std::runtime_error If the reaction is not part of the active reactions in the adaptive engine view.
*/
double calculateMolarReactionFlow(
[[nodiscard]] double calculateMolarReactionFlow(
const reaction::Reaction &reaction,
const std::vector<double> &Y_culled,
double T9,
@@ -196,7 +198,7 @@ namespace gridfire {
*
* @return Reference to the LogicalReactionSet containing all active reactions.
*/
const reaction::LogicalReactionSet& getNetworkReactions() const override;
[[nodiscard]] const reaction::LogicalReactionSet& getNetworkReactions() const override;
/**
* @brief Computes timescales for all active species in the network.
@@ -211,7 +213,7 @@ namespace gridfire {
*
* @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called).
*/
std::unordered_map<fourdst::atomic::Species, double> getSpeciesTimescales(
[[nodiscard]] std::unordered_map<fourdst::atomic::Species, double> getSpeciesTimescales(
const std::vector<double> &Y_culled,
double T9,
double rho
@@ -221,7 +223,11 @@ namespace gridfire {
* @brief Gets the base engine.
* @return A const reference to the base engine.
*/
const DynamicEngine& getBaseEngine() const override { return m_baseEngine; }
[[nodiscard]] const DynamicEngine& getBaseEngine() const override { return m_baseEngine; }
void setScreeningModel(screening::ScreeningType model) override;
[[nodiscard]] screening::ScreeningType getScreeningModel() const override;
private:
using Config = fourdst::config::Config;
using LogManager = fourdst::logging::LogManager;
@@ -257,7 +263,7 @@ namespace gridfire {
*
* @see AdaptiveEngineView::update()
*/
std::vector<size_t> constructSpeciesIndexMap() const;
[[nodiscard]] std::vector<size_t> constructSpeciesIndexMap() const;
/**
* @brief Constructs the reaction index map.
@@ -269,7 +275,7 @@ namespace gridfire {
*
* @see AdaptiveEngineView::update()
*/
std::vector<size_t> constructReactionIndexMap() const;
[[nodiscard]] std::vector<size_t> constructReactionIndexMap() const;
/**
* @brief Maps a vector of culled abundances to a vector of full abundances.
@@ -278,7 +284,7 @@ namespace gridfire {
* @return A vector of abundances for the full network, with the abundances of the active
* species copied from the culled vector.
*/
std::vector<double> mapCulledToFull(const std::vector<double>& culled) const;
[[nodiscard]] std::vector<double> mapCulledToFull(const std::vector<double>& culled) const;
/**
* @brief Maps a vector of full abundances to a vector of culled abundances.
@@ -287,7 +293,7 @@ namespace gridfire {
* @return A vector of abundances for the active species, with the abundances of the active
* species copied from the full vector.
*/
std::vector<double> mapFullToCulled(const std::vector<double>& full) const;
[[nodiscard]] std::vector<double> mapFullToCulled(const std::vector<double>& full) const;
/**
* @brief Maps a culled species index to a full species index.
@@ -297,7 +303,7 @@ namespace gridfire {
*
* @throws std::out_of_range If the culled index is out of bounds for the species index map.
*/
size_t mapCulledToFullSpeciesIndex(size_t culledSpeciesIndex) const;
[[nodiscard]] size_t mapCulledToFullSpeciesIndex(size_t culledSpeciesIndex) const;
/**
* @brief Maps a culled reaction index to a full reaction index.
@@ -307,7 +313,7 @@ namespace gridfire {
*
* @throws std::out_of_range If the culled index is out of bounds for the reaction index map.
*/
size_t mapCulledToFullReactionIndex(size_t culledReactionIndex) const;
[[nodiscard]] size_t mapCulledToFullReactionIndex(size_t culledReactionIndex) const;
/**
* @brief Validates that the AdaptiveEngineView is not stale.
@@ -320,10 +326,10 @@ namespace gridfire {
const NetIn& netIn,
std::vector<double>& out_Y_Full
) const;
std::unordered_set<fourdst::atomic::Species> findReachableSpecies(
[[nodiscard]] std::unordered_set<fourdst::atomic::Species> findReachableSpecies(
const NetIn& netIn
) const;
std::vector<const reaction::LogicalReaction*> cullReactionsByFlow(
[[nodiscard]] std::vector<const reaction::LogicalReaction*> cullReactionsByFlow(
const std::vector<ReactionFlow>& allFlows,
const std::unordered_set<fourdst::atomic::Species>& reachableSpecies,
const std::vector<double>& Y_full,

34
src/network/include/gridfire/engine/views/engine_defined.h
View File

@@ -1,7 +1,9 @@
#pragma once
#include "engine_view_abstract.h"
#include "../engine_abstract.h"
#include "gridfire/engine/views/engine_view_abstract.h"
#include "gridfire/engine/engine_abstract.h"
#include "gridfire/io/network_file.h"
#include "gridfire/network.h"
#include "fourdst/config/config.h"
#include "fourdst/logging/logging.h"
@@ -13,7 +15,11 @@
namespace gridfire{
class FileDefinedEngineView final: public DynamicEngine, public EngineView<DynamicEngine> {
public:
explicit FileDefinedEngineView(DynamicEngine& baseEngine, const std::string& fileName);
explicit FileDefinedEngineView(
DynamicEngine& baseEngine,
const std::string& fileName,
const io::NetworkFileParser& parser
);
// --- EngineView Interface ---
const DynamicEngine& getBaseEngine() const override;
@@ -53,6 +59,14 @@ namespace gridfire{
const double T9,
const double rho
) const override;
void update(const NetIn &netIn) override;
void setNetworkFile(const std::string& fileName);
void setScreeningModel(screening::ScreeningType model) override;
[[nodiscard]] screening::ScreeningType getScreeningModel() const override;
private:
using Config = fourdst::config::Config;
using LogManager = fourdst::logging::LogManager;
@@ -60,12 +74,17 @@ namespace gridfire{
quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
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.
std::vector<fourdst::atomic::Species> m_activeSpecies; ///< Active species in the defined engine.
reaction::LogicalReactionSet m_activeReactions; ///< Active reactions in the defined engine.
std::vector<size_t> m_speciesIndexMap; ///< Maps indices of active species to indices in the full network.
std::vector<size_t> m_reactionIndexMap; ///< Maps indices of active reactions to indices in the full network.
bool m_isStale = true;
private:
void buildFromFile(const std::string& fileName);
@@ -100,7 +119,7 @@ namespace gridfire{
* @return A vector of abundances for the full network, with the abundances of the active
* species copied from the culled vector.
*/
std::vector<double> mapCulledToFull(const std::vector<double>& culled) const;
std::vector<double> mapViewToFull(const std::vector<double>& culled) const;
/**
* @brief Maps a vector of full abundances to a vector of culled abundances.
@@ -109,7 +128,7 @@ namespace gridfire{
* @return A vector of abundances for the active species, with the abundances of the active
* species copied from the full vector.
*/
std::vector<double> mapFullToCulled(const std::vector<double>& full) const;
std::vector<double> mapFullToView(const std::vector<double>& full) const;
/**
* @brief Maps a culled species index to a full species index.
@@ -119,7 +138,7 @@ namespace gridfire{
*
* @throws std::out_of_range If the culled index is out of bounds for the species index map.
*/
size_t mapCulledToFullSpeciesIndex(size_t culledSpeciesIndex) const;
size_t mapViewToFullSpeciesIndex(size_t culledSpeciesIndex) const;
/**
* @brief Maps a culled reaction index to a full reaction index.
@@ -129,7 +148,8 @@ namespace gridfire{
*
* @throws std::out_of_range If the culled index is out of bounds for the reaction index map.
*/
size_t mapCulledToFullReactionIndex(size_t culledReactionIndex) const;
size_t mapViewToFullReactionIndex(size_t culledReactionIndex) const;
void validateNetworkState() const;
};
}

48
src/network/include/gridfire/io/network_file.h
View File

@@ -0,0 +1,48 @@
#pragma once
#include "fourdst/config/config.h"
#include "fourdst/logging/logging.h"
#include "quill/Logger.h"
#include <string>
#include <vector>
namespace gridfire::io {
struct ParsedNetworkData {
std::vector<std::string> reactionPENames;
};
class NetworkFileParser {
public:
virtual ~NetworkFileParser() = default;
[[nodiscard]] virtual ParsedNetworkData parse(const std::string& filename) const = 0;
};
class SimpleReactionListFileParser final : public NetworkFileParser {
public:
explicit SimpleReactionListFileParser();
ParsedNetworkData parse(const std::string& filename) const override;
private:
using Config = fourdst::config::Config;
using LogManager = fourdst::logging::LogManager;
Config& m_config = Config::getInstance();
quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
};
class MESANetworkFileParser final : public NetworkFileParser {
public:
explicit MESANetworkFileParser(const std::string& filename);
ParsedNetworkData parse(const std::string& filename) const override;
private:
using Config = fourdst::config::Config;
using LogManager = fourdst::logging::LogManager;
Config& m_config = Config::getInstance();
quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
std::string m_filename;
};
}

509
src/network/include/gridfire/reaction/reaction.h
View File

@@ -11,6 +11,7 @@
#include "cppad/cppad.hpp"
#include "xxhash64.h"
/**
* @file reaction.h
@@ -252,6 +253,10 @@ namespace gridfire::reaction {
*/
[[nodiscard]] uint64_t hash(uint64_t seed = 0) const;
friend std::ostream& operator<<(std::ostream& os, const Reaction& r) {
return os << "(Reaction:" << r.m_id << ")";
}
protected:
quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
std::string m_id; ///< Unique identifier for the reaction (e.g., "h1+h1=>h2+e+nu").
@@ -288,162 +293,6 @@ namespace gridfire::reaction {
}
};
/**
* @class ReactionSet
* @brief A collection of Reaction objects.
*
* This class manages a set of individual `Reaction` objects, providing
* efficient lookup by ID and functionality to query the entire set.
*
* Example:
* @code
* ReactionSet my_set({reaction1, reaction2});
* my_set.add_reaction(reaction3);
* if (my_set.contains("h1(p,g)h2")) {
* const Reaction& r = my_set["h1(p,g)h2"];
* }
* @endcode
*/
class ReactionSet {
public:
/**
* @brief Constructs a ReactionSet from a vector of reactions.
* @param reactions The initial vector of Reaction objects.
*/
explicit ReactionSet(std::vector<Reaction> reactions);
/**
* @brief Copy constructor.
* @param other The ReactionSet to copy.
*/
ReactionSet(const ReactionSet& other);
/**
* @brief Copy assignment operator.
* @param other The ReactionSet to assign from.
* @return A reference to this ReactionSet.
*/
ReactionSet& operator=(const ReactionSet& other);
/**
* @brief Virtual destructor.
*/
virtual ~ReactionSet() = default;
/**
* @brief Adds a reaction to the set.
* @param reaction The Reaction to add.
*/
virtual void add_reaction(Reaction reaction);
/**
* @brief Removes a reaction from the set.
* @param reaction The Reaction to remove.
*/
virtual void remove_reaction(const Reaction& reaction);
/**
* @brief Checks if the set contains a reaction with the given ID.
* @param id The ID of the reaction to find.
* @return True if the reaction is in the set, false otherwise.
*/
[[nodiscard]] bool contains(const std::string_view& id) const;
/**
* @brief Checks if the set contains the given reaction.
* @param reaction The Reaction to find.
* @return True if the reaction is in the set, false otherwise.
*/
[[nodiscard]] bool contains(const Reaction& reaction) const;
/**
* @brief Gets the number of reactions in the set.
* @return The size of the set.
*/
[[nodiscard]] virtual size_t size() const { return m_reactions.size(); }
/**
* @brief Removes all reactions from the set.
*/
void clear();
/**
* @brief Checks if any reaction in the set involves the given species.
* @param species The species to check for.
* @return True if the species is involved in any reaction.
*/
[[nodiscard]] bool contains_species(const fourdst::atomic::Species& species) const;
/**
* @brief Checks if any reaction in the set contains the given species as a reactant.
* @param species The species to check for.
* @return True if the species is a reactant in any reaction.
*/
[[nodiscard]] bool contains_reactant(const fourdst::atomic::Species& species) const;
/**
* @brief Checks if any reaction in the set contains the given species as a product.
* @param species The species to check for.
* @return True if the species is a product in any reaction.
*/
[[nodiscard]] bool contains_product(const fourdst::atomic::Species& species) const;
/**
* @brief Accesses a reaction by its index.
* @param index The index of the reaction to access.
* @return A const reference to the Reaction.
* @throws std::out_of_range if the index is out of bounds.
*/
[[nodiscard]] virtual const Reaction& operator[](size_t index) const;
/**
* @brief Accesses a reaction by its ID.
* @param id The ID of the reaction to access.
* @return A const reference to the Reaction.
* @throws std::out_of_range if no reaction with the given ID exists.
*/
[[nodiscard]] const Reaction& operator[](const std::string_view& id) const;
/**
* @brief Compares this set with another for equality.
* @param other The other ReactionSet to compare with.
* @return True if the sets are equal (same size and hash).
*/
bool operator==(const ReactionSet& other) const;
/**
* @brief Compares this set with another for inequality.
* @param other The other ReactionSet to compare with.
* @return True if the sets are not equal.
*/
bool operator!=(const ReactionSet& other) const;
/**
* @brief Computes a hash for the entire set.
* @param seed The seed for the hash function.
* @return A 64-bit hash value.
* @details The algorithm computes the hash of each individual reaction,
* sorts the hashes, and then computes a final hash over the sorted list
* of hashes. This ensures the hash is order-independent.
*/
[[nodiscard]] uint64_t hash(uint64_t seed = 0) const;
/** @name Iterators
* Provides iterators to loop over the reactions in the set.
*/
///@{
auto begin() { return m_reactions.begin(); }
[[nodiscard]] auto begin() const { return m_reactions.cbegin(); }
auto end() { return m_reactions.end(); }
[[nodiscard]] auto end() const { return m_reactions.cend(); }
///@}
private:
quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
std::vector<Reaction> m_reactions;
std::string m_id;
std::unordered_map<std::string, Reaction> m_reactionNameMap; ///< Maps reaction IDs to Reaction objects for quick lookup.
};
/**
@@ -508,6 +357,12 @@ namespace gridfire::reaction {
auto end() { return m_rates.end(); }
[[nodiscard]] auto end() const { return m_rates.cend(); }
///@}
///
friend std::ostream& operator<<(std::ostream& os, const LogicalReaction& r) {
os << "(LogicalReaction: " << r.id() << ", reverse: " << r.is_reverse() << ")";
return os;
}
private:
std::vector<std::string> m_sources; ///< List of source labels.
@@ -543,31 +398,128 @@ namespace gridfire::reaction {
}
};
/**
* @class LogicalReactionSet
* @brief A collection of LogicalReaction objects.
*
* This class takes a `ReactionSet` and groups individual `Reaction` objects
* into `LogicalReaction` objects based on their `peName`. This provides a
* view of the network where all rates for the same physical process are combined.
*/
class LogicalReactionSet final : public ReactionSet {
template <typename ReactionT>
class TemplatedReactionSet final {
public:
/**
* @brief Deleted default constructor.
* @brief Constructs a ReactionSet from a vector of reactions.
* @param reactions The initial vector of Reaction objects.
*/
LogicalReactionSet() = delete;
explicit TemplatedReactionSet(std::vector<ReactionT> reactions);
/**
* @brief Constructs a LogicalReactionSet from a ReactionSet.
* @param reactionSet The set of individual reactions to group.
* @details This constructor iterates through the provided `ReactionSet`,
* groups reactions by their `peName`, and creates a `LogicalReaction` for each group.
* @brief Copy constructor.
* @param other The ReactionSet to copy.
*/
explicit LogicalReactionSet(const ReactionSet& reactionSet);
TemplatedReactionSet(const TemplatedReactionSet<ReactionT>& other);
/**
* @brief Copy assignment operator.
* @param other The ReactionSet to assign from.
* @return A reference to this ReactionSet.
*/
TemplatedReactionSet<ReactionT>& operator=(const TemplatedReactionSet<ReactionT>& other);
/**
* @brief Adds a reaction to the set.
* @param reaction The Reaction to add.
*/
void add_reaction(ReactionT reaction);
/**
* @brief Removes a reaction from the set.
* @param reaction The Reaction to remove.
*/
void remove_reaction(const ReactionT& reaction);
/**
* @brief Checks if the set contains a reaction with the given ID.
* @param id The ID of the reaction to find.
* @return True if the reaction is in the set, false otherwise.
*/
[[nodiscard]] bool contains(const std::string_view& id) const;
/**
* @brief Checks if the set contains the given reaction.
* @param reaction The Reaction to find.
* @return True if the reaction is in the set, false otherwise.
*/
[[nodiscard]] bool contains(const Reaction& reaction) const;
/**
* @brief Gets the number of reactions in the set.
* @return The size of the set.
*/
[[nodiscard]] size_t size() const { return m_reactions.size(); }
/**
* @brief Removes all reactions from the set.
*/
void clear();
/**
* @brief Checks if any reaction in the set involves the given species.
* @param species The species to check for.
* @return True if the species is involved in any reaction.
*/
[[nodiscard]] bool contains_species(const fourdst::atomic::Species& species) const;
/**
* @brief Checks if any reaction in the set contains the given species as a reactant.
* @param species The species to check for.
* @return True if the species is a reactant in any reaction.
*/
[[nodiscard]] bool contains_reactant(const fourdst::atomic::Species& species) const;
/**
* @brief Checks if any reaction in the set contains the given species as a product.
* @param species The species to check for.
* @return True if the species is a product in any reaction.
*/
[[nodiscard]] bool contains_product(const fourdst::atomic::Species& species) const;
/**
* @brief Accesses a reaction by its index.
* @param index The index of the reaction to access.
* @return A const reference to the Reaction.
* @throws std::out_of_range if the index is out of bounds.
*/
[[nodiscard]] const ReactionT& operator[](size_t index) const;
/**
* @brief Accesses a reaction by its ID.
* @param id The ID of the reaction to access.
* @return A const reference to the Reaction.
* @throws std::out_of_range if no reaction with the given ID exists.
*/
[[nodiscard]] const ReactionT& operator[](const std::string_view& id) const;
/**
* @brief Compares this set with another for equality.
* @param other The other ReactionSet to compare with.
* @return True if the sets are equal (same size and hash).
*/
bool operator==(const TemplatedReactionSet& other) const;
/**
* @brief Compares this set with another for inequality.
* @param other The other ReactionSet to compare with.
* @return True if the sets are not equal.
*/
bool operator!=(const TemplatedReactionSet& other) const;
/**
* @brief Computes a hash for the entire set.
* @param seed The seed for the hash function.
* @return A 64-bit hash value.
* @details The algorithm computes the hash of each individual reaction,
* sorts the hashes, and then computes a final hash over the sorted list
* of hashes. This ensures the hash is order-independent.
*/
[[nodiscard]] uint64_t hash(uint64_t seed = 0) const;
/** @name Iterators
* Provides iterators to loop over the logical reactions in the set.
* Provides iterators to loop over the reactions in the set.
*/
///@{
auto begin() { return m_reactions.begin(); }
@@ -575,25 +527,208 @@ namespace gridfire::reaction {
auto end() { return m_reactions.end(); }
[[nodiscard]] auto end() const { return m_reactions.cend(); }
///@}
///
friend std::ostream& operator<<(std::ostream& os, const TemplatedReactionSet<ReactionT>& r) {
os << "(ReactionSet: [";
int counter = 0;
for (const auto& reaction : r.m_reactions) {
os << reaction;
if (counter < r.m_reactions.size() - 2) {
os << ", ";
} else if (counter == r.m_reactions.size() - 2) {
os << " and ";
}
++counter;
}
os << "])";
return os;
}
/**
* @brief Gets the number of logical reactions in the set.
* @return The size of the set.
*/
[[nodiscard]] size_t size() const { return m_reactions.size(); }
/**
* @brief Accesses a logical reaction by its index.
* @param index The index of the logical reaction.
* @return A const reference to the LogicalReaction.
*/
[[nodiscard]] const LogicalReaction& operator[](size_t index) const { return m_reactions[index]; }
[[nodiscard]] std::unordered_set<fourdst::atomic::Species> getReactionSetSpecies() const;
private:
quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
std::vector<LogicalReaction> m_reactions;
std::vector<ReactionT> m_reactions;
std::string m_id;
std::unordered_map<std::string, LogicalReaction> m_reactionNameMap; ///< Maps reaction IDs to LogicalReaction objects for quick lookup.
std::unordered_map<std::string, ReactionT> m_reactionNameMap; ///< Maps reaction IDs to Reaction objects for quick lookup.
};
using ReactionSet = TemplatedReactionSet<Reaction>; ///< A set of reactions, typically from a single source like REACLIB.
using LogicalReactionSet = TemplatedReactionSet<LogicalReaction>; ///< A set of logical reactions.
LogicalReactionSet packReactionSetToLogicalReactionSet(const ReactionSet& reactionSet);
template <typename ReactionT>
TemplatedReactionSet<ReactionT>::TemplatedReactionSet(
std::vector<ReactionT> reactions
) :
m_reactions(std::move(reactions)) {
if (m_reactions.empty()) {
return; // Case where the reactions will be added later.
}
m_reactionNameMap.reserve(reactions.size());
for (const auto& reaction : m_reactions) {
m_id += reaction.id();
m_reactionNameMap.emplace(reaction.id(), reaction);
}
}
template <typename ReactionT>
TemplatedReactionSet<ReactionT>::TemplatedReactionSet(const TemplatedReactionSet<ReactionT> &other) {
m_reactions.reserve(other.m_reactions.size());
for (const auto& reaction_ptr: other.m_reactions) {
m_reactions.push_back(reaction_ptr);
}
m_reactionNameMap.reserve(other.m_reactionNameMap.size());
for (const auto& reaction_ptr : m_reactions) {
m_reactionNameMap.emplace(reaction_ptr.id(), reaction_ptr);
}
}
template <typename ReactionT>
TemplatedReactionSet<ReactionT>& TemplatedReactionSet<ReactionT>::operator=(const TemplatedReactionSet<ReactionT> &other) {
if (this != &other) {
TemplatedReactionSet temp(other);
std::swap(m_reactions, temp.m_reactions);
std::swap(m_reactionNameMap, temp.m_reactionNameMap);
}
return *this;
}
template <typename ReactionT>
void TemplatedReactionSet<ReactionT>::add_reaction(ReactionT reaction) {
m_reactions.emplace_back(reaction);
m_id += m_reactions.back().id();
m_reactionNameMap.emplace(m_reactions.back().id(), m_reactions.back());
}
template <typename ReactionT>
void TemplatedReactionSet<ReactionT>::remove_reaction(const ReactionT& reaction) {
if (!m_reactionNameMap.contains(std::string(reaction.id()))) {
return;
}
m_reactionNameMap.erase(std::string(reaction.id()));
std::erase_if(m_reactions, [&reaction](const Reaction& r) {
return r == reaction;
});
}
template <typename ReactionT>
bool TemplatedReactionSet<ReactionT>::contains(const std::string_view& id) const {
for (const auto& reaction : m_reactions) {
if (reaction.id() == id) {
return true;
}
}
return false;
}
template <typename ReactionT>
bool TemplatedReactionSet<ReactionT>::contains(const Reaction& reaction) const {
for (const auto& r : m_reactions) {
if (r == reaction) {
return true;
}
}
return false;
}
template <typename ReactionT>
void TemplatedReactionSet<ReactionT>::clear() {
m_reactions.clear();
m_reactionNameMap.clear();
}
template <typename ReactionT>
bool TemplatedReactionSet<ReactionT>::contains_species(const fourdst::atomic::Species& species) const {
for (const auto& reaction : m_reactions) {
if (reaction.contains(species)) {
return true;
}
}
return false;
}
template <typename ReactionT>
bool TemplatedReactionSet<ReactionT>::contains_reactant(const fourdst::atomic::Species& species) const {
for (const auto& r : m_reactions) {
if (r.contains_reactant(species)) {
return true;
}
}
return false;
}
template <typename ReactionT>
bool TemplatedReactionSet<ReactionT>::contains_product(const fourdst::atomic::Species& species) const {
for (const auto& r : m_reactions) {
if (r.contains_product(species)) {
return true;
}
}
return false;
}
template <typename ReactionT>
const ReactionT& TemplatedReactionSet<ReactionT>::operator[](const size_t index) const {
if (index >= m_reactions.size()) {
m_logger -> flush_log();
throw std::out_of_range("Index" + std::to_string(index) + " out of range for ReactionSet of size " + std::to_string(m_reactions.size()) + ".");
}
return m_reactions[index];
}
template <typename ReactionT>
const ReactionT& TemplatedReactionSet<ReactionT>::operator[](const std::string_view& id) const {
if (auto it = m_reactionNameMap.find(std::string(id)); it != m_reactionNameMap.end()) {
return it->second;
}
m_logger -> flush_log();
throw std::out_of_range("Species " + std::string(id) + " does not exist in ReactionSet.");
}
template <typename ReactionT>
bool TemplatedReactionSet<ReactionT>::operator==(const TemplatedReactionSet<ReactionT>& other) const {
if (size() != other.size()) {
return false;
}
return hash() == other.hash();
}
template <typename ReactionT>
bool TemplatedReactionSet<ReactionT>::operator!=(const TemplatedReactionSet<ReactionT>& other) const {
return !(*this == other);
}
template <typename ReactionT>
uint64_t TemplatedReactionSet<ReactionT>::hash(uint64_t seed) const {
if (m_reactions.empty()) {
return XXHash64::hash(nullptr, 0, seed);
}
std::vector<uint64_t> individualReactionHashes;
individualReactionHashes.reserve(m_reactions.size());
for (const auto& reaction : m_reactions) {
individualReactionHashes.push_back(reaction.hash(seed));
}
std::ranges::sort(individualReactionHashes);
const auto data = static_cast<const void*>(individualReactionHashes.data());
const size_t sizeInBytes = individualReactionHashes.size() * sizeof(uint64_t);
return XXHash64::hash(data, sizeInBytes, seed);
}
template<typename ReactionT>
std::unordered_set<fourdst::atomic::Species> TemplatedReactionSet<ReactionT>::getReactionSetSpecies() const {
std::unordered_set<fourdst::atomic::Species> species;
for (const auto& reaction : m_reactions) {
const auto reactionSpecies = reaction.all_species();
species.insert(reactionSpecies.begin(), reactionSpecies.end());
}
return species;
}
}

33
src/network/include/gridfire/screening/screening_abstract.h Normal file
View File

@@ -0,0 +1,33 @@
#pragma once
#include "gridfire/reaction/reaction.h"
#include "fourdst/composition/atomicSpecies.h"
#include "cppad/cppad.hpp"
#include <vector>
namespace gridfire::screening {
class ScreeningModel {
public:
using ADDouble = CppAD::AD<double>;
virtual ~ScreeningModel() = default;
virtual std::vector<double> calculateScreeningFactors(
const reaction::LogicalReactionSet& reactions,
const std::vector<fourdst::atomic::Species>& species,
const std::vector<double>& Y,
const double T9,
const double rho
) const = 0;
virtual std::vector<ADDouble> calculateScreeningFactors(
const reaction::LogicalReactionSet& reactions,
const std::vector<fourdst::atomic::Species>& species,
const std::vector<ADDouble>& Y,
const ADDouble T9,
const ADDouble rho
) const = 0;
};
}

48
src/network/include/gridfire/screening/screening_bare.h Normal file
View File

@@ -0,0 +1,48 @@
#pragma once
#include "gridfire/screening/screening_abstract.h"
#include "gridfire/reaction/reaction.h"
#include "cppad/cppad.hpp"
namespace gridfire::screening {
class BareScreeningModel final : public ScreeningModel {
using ADDouble = CppAD::AD<double>;
public:
[[nodiscard]] std::vector<double> calculateScreeningFactors(
const reaction::LogicalReactionSet& reactions,
const std::vector<fourdst::atomic::Species>& species,
const std::vector<double>& Y,
const double T9,
const double rho
) const override;
[[nodiscard]] std::vector<ADDouble> calculateScreeningFactors(
const reaction::LogicalReactionSet& reactions,
const std::vector<fourdst::atomic::Species>& species,
const std::vector<ADDouble>& Y,
const ADDouble T9,
const ADDouble rho
) const override;
private:
template <typename T>
[[nodiscard]] std::vector<T> calculateFactors_impl(
const reaction::LogicalReactionSet& reactions,
const std::vector<fourdst::atomic::Species>& species,
const std::vector<T>& Y,
const T T9,
const T rho
) const;
};
template<typename T>
std::vector<T> BareScreeningModel::calculateFactors_impl(
const reaction::LogicalReactionSet &reactions,
const std::vector<fourdst::atomic::Species> &species,
const std::vector<T> &Y,
const T T9,
const T rho
) const {
return std::vector<T>(reactions.size(), T(1.0)); // Bare screening returns 1.0 for all reactions
}
}

14
src/network/include/gridfire/screening/screening_types.h Normal file
View File

@@ -0,0 +1,14 @@
#pragma once
#include "gridfire/screening/screening_abstract.h"
#include <memory>
namespace gridfire::screening {
enum class ScreeningType {
BARE, ///< No screening applied
WEAK, ///< Weak screening model
};
std::unique_ptr<ScreeningModel> selectScreeningModel(ScreeningType type);
}

117
src/network/include/gridfire/screening/screening_weak.h Normal file
View File

@@ -0,0 +1,117 @@
#pragma once
#include "gridfire/screening/screening_abstract.h"
#include "gridfire/reaction/reaction.h"
#include "fourdst/logging/logging.h"
#include "quill/Logger.h"
#include "quill/LogMacros.h"
#include "cppad/cppad.hpp"
namespace gridfire::screening {
class WeakScreeningModel final : public ScreeningModel {
public:
[[nodiscard]] std::vector<double> calculateScreeningFactors(
const reaction::LogicalReactionSet& reactions,
const std::vector<fourdst::atomic::Species>& species,
const std::vector<double>& Y,
const double T9,
const double rho
) const override;
[[nodiscard]] std::vector<CppAD::AD<double>> calculateScreeningFactors(
const reaction::LogicalReactionSet& reactions,
const std::vector<fourdst::atomic::Species>& species,
const std::vector<CppAD::AD<double>>& Y,
const CppAD::AD<double> T9,
const CppAD::AD<double> rho
) const override;
private:
quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
private:
template <typename T>
[[nodiscard]] std::vector<T> calculateFactors_impl(
const reaction::LogicalReactionSet& reactions,
const std::vector<fourdst::atomic::Species>& species,
const std::vector<T>& Y,
const T T9,
const T rho
) const;
};
template <typename T>
std::vector<T> WeakScreeningModel::calculateFactors_impl(
const reaction::LogicalReactionSet& reactions,
const std::vector<fourdst::atomic::Species>& species,
const std::vector<T>& Y,
const T T9,
const T rho
) const {
LOG_TRACE_L1(
m_logger,
"Calculating weak screening factors for {} reactions...",
reactions.size()
);
// --- CppAD Safe low temp checking ---
const T zero(0.0);
const T one(1.0);
const T low_temp_threshold(1e-9);
const T low_T_flag = CppAD::CondExpLt(T9, low_temp_threshold, zero, one);
// --- Calculate composition-dependent terms ---
// ζ = ∑(Z_i^2 + Z_i) * X_i / A_i
// This simplifies somewhat to ζ = ∑ (Z_i^2 + Z_i) * Y_i
// Where Y_i is the molar abundance (mol/g)
T zeta(0.0);
for (size_t i = 0; i < species.size(); ++i) {
const T Z = species[i].m_z;
zeta += (Z * Z + Z) * Y[i];
}
// --- Constant prefactors ---
const T T7 = T9 * static_cast<T>(100.00);
const T T7_safe = CppAD::CondExpLe(T7, low_temp_threshold, low_temp_threshold, T7);
const T prefactor = static_cast<T>(0.188) * CppAD::sqrt(rho / (T7_safe * T7_safe * T7_safe)) * CppAD::sqrt(zeta);
// --- Loop through reactions and calculate screening factors for each ---
std::vector<T> factors;
factors.reserve(reactions.size());
for (const auto& reaction : reactions) {
T H_12(0.0); // screening abundance term
const auto& reactants = reaction.reactants();
const bool isTripleAlpha = (
reactants.size() == 3 &&
reactants[0].m_z == 2 &&
reactants[1].m_z == 2 &&
reactants[2].m_z == 2 &&
reactants[0] == reactants[1] &&
reactants[1] == reactants[2]
);
if (reactants.size() == 2) {
LOG_TRACE_L3(m_logger, "Calculating screening factor for reaction: {}", reaction.peName());
const T Z1 = static_cast<T>(reactants[0].m_z);
const T Z2 = static_cast<T>(reactants[1].m_z);
H_12 = prefactor * Z1 * Z2;
}
else if (isTripleAlpha) {
LOG_TRACE_L3(m_logger, "Special case for triple alpha process in reaction: {}", reaction.peName());
// Special case for triple alpha process
const T Z_alpha = static_cast<T>(2.0);
const T H_alpha_alpha = prefactor * Z_alpha * Z_alpha;
H_12 = static_cast<T>(3.0) * H_alpha_alpha; // Triple alpha process
}
// For 1 body reactions H_12 remains 0 so e^H_12 will be 1.0 (screening does not apply)
// Aside from triple alpha, all other astrophysically relevant reactions are 2-body in the weak screening regime
H_12 *= low_T_flag; // Apply low temperature flag to screening factor
H_12 = CppAD::CondExpGe(H_12, static_cast<T>(2.0), static_cast<T>(2.0), H_12); // Caps the screening factor at 10 to avoid numerical issues
factors.push_back(CppAD::exp(H_12));
// std::cout << "Screening factor: " << reaction.peName() << " : " << factors.back() << "(" << H_12 << ")" << std::endl;
}
return factors;
}
}

78
src/network/include/gridfire/solver/solver.h
View File

@@ -2,7 +2,7 @@
#include "gridfire/engine/engine_graph.h"
#include "gridfire/engine/engine_abstract.h"
#include "gridfire/engine/engine_adaptive.h"
#include "../engine/views/engine_adaptive.h"
#include "gridfire/network.h"
#include "fourdst/logging/logging.h"
@@ -10,7 +10,7 @@
#include "quill/Logger.h"
#include "Eigen/Dense"
#include "unsupported/Eigen/NonLinearOptimization" // Required for LevenbergMarquardt
#include <vector>
@@ -95,13 +95,13 @@ namespace gridfire::solver {
* @see AdaptiveEngineView
* @see DynamicEngine::getSpeciesTimescales()
*/
class QSENetworkSolver final : public AdaptiveNetworkSolverStrategy {
class QSENetworkSolver final : public DynamicNetworkSolverStrategy {
public:
/**
* @brief Constructor for the QSENetworkSolver.
* @param engine The adaptive engine view to use for evaluating the network.
*/
using AdaptiveNetworkSolverStrategy::AdaptiveNetworkSolverStrategy;
using DynamicNetworkSolverStrategy::DynamicNetworkSolverStrategy;
/**
* @brief Evaluates the network for a given timestep using the QSE approach.
@@ -310,9 +310,13 @@ namespace gridfire::solver {
using InputType = Eigen::Matrix<T, Eigen::Dynamic, 1>;
using OutputType = Eigen::Matrix<T, Eigen::Dynamic, 1>;
using JacobianType = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;
enum {
InputsAtCompileTime = Eigen::Dynamic,
ValuesAtCompileTime = Eigen::Dynamic
};
DynamicEngine& m_engine; ///< The engine used to evaluate the network.
const std::vector<double>& m_YDynamic; ///< Abundances of the dynamic species.
const std::vector<double>& m_YFull; ///< The full, initial abundance vector
const std::vector<size_t>& m_dynamicSpeciesIndices; ///< Indices of the dynamic species.
const std::vector<size_t>& m_QSESpeciesIndices; ///< Indices of the QSE species.
const double m_T9; ///< Temperature in units of 10^9 K.
@@ -321,7 +325,7 @@ namespace gridfire::solver {
/**
* @brief Constructor for the EigenFunctor.
* @param engine The engine used to evaluate the network.
* @param YDynamic Abundances of the dynamic species.
* @param YFull Abundances of the dynamic species.
* @param dynamicSpeciesIndices Indices of the dynamic species.
* @param QSESpeciesIndices Indices of the QSE species.
* @param T9 Temperature in units of 10^9 K.
@@ -329,34 +333,37 @@ namespace gridfire::solver {
*/
EigenFunctor(
DynamicEngine& engine,
const std::vector<double>& YDynamic,
const std::vector<double>& YFull,
const std::vector<size_t>& dynamicSpeciesIndices,
const std::vector<size_t>& QSESpeciesIndices,
const double T9,
const double rho
) :
m_engine(engine),
m_YDynamic(YDynamic),
m_YFull(YFull),
m_dynamicSpeciesIndices(dynamicSpeciesIndices),
m_QSESpeciesIndices(QSESpeciesIndices),
m_T9(T9),
m_rho(rho) {}
int values() const { return m_QSESpeciesIndices.size(); }
int inputs() const { return m_QSESpeciesIndices.size(); }
/**
* @brief Calculates the residual vector for the QSE species.
* @param v_QSE Input vector of QSE species abundances (logarithmic).
* @param v_QSE_log Input vector of QSE species abundances (logarithmic).
* @param f_QSE Output vector of residuals.
* @return 0 for success.
*/
int operator()(const InputType& v_QSE, OutputType& f_QSE) const;
int operator()(const InputType& v_QSE_log, OutputType& f_QSE) const;
/**
* @brief Calculates the Jacobian matrix for the QSE species.
* @param v_QSE Input vector of QSE species abundances (logarithmic).
* @param v_QSE_log Input vector of QSE species abundances (logarithmic).
* @param J_QSE Output Jacobian matrix.
* @return 0 for success.
*/
int df(const InputType& v_QSE, JacobianType& J_QSE) const;
int df(const InputType& v_QSE_log, JacobianType& J_QSE) const;
};
private:
quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); ///< Logger instance.
@@ -404,6 +411,7 @@ namespace gridfire::solver {
const double m_T9; ///< Temperature in units of 10^9 K.
const double m_rho; ///< Density in g/cm^3.
const size_t m_numSpecies; ///< The number of species in the network.
quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); ///< Logger instance.
/**
* @brief Constructor for the RHSFunctor.
@@ -485,49 +493,45 @@ namespace gridfire::solver {
};
template<typename T>
int QSENetworkSolver::EigenFunctor<T>::operator()(const InputType &v_QSE, OutputType &f_QSE) const {
std::vector<double> YFull(m_engine.getNetworkSpecies().size(), 0.0);
Eigen::VectorXd Y_QSE(v_QSE.array().exp());
for (size_t i = 0; i < m_dynamicSpeciesIndices.size(); ++i) {
YFull[m_dynamicSpeciesIndices[i]] = m_YDynamic[i];
}
int QSENetworkSolver::EigenFunctor<T>::operator()(const InputType &v_QSE_log, OutputType &f_QSE) const {
std::vector<double> y = m_YFull; // Full vector of species abundances
Eigen::VectorXd Y_QSE = v_QSE_log.array().exp();
for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) {
YFull[m_QSESpeciesIndices[i]] = Y_QSE(i);
}
const auto [full_dYdt, specificEnergyGenerationRate] = m_engine.calculateRHSAndEnergy(YFull, m_T9, m_rho);
f_QSE.resize(m_QSESpeciesIndices.size());
for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) {
f_QSE(i) = full_dYdt[m_QSESpeciesIndices[i]];
y[m_QSESpeciesIndices[i]] = Y_QSE(i);
}
const auto [dydt, specificEnergyRate] = m_engine.calculateRHSAndEnergy(y, m_T9, m_rho);
f_QSE.resize(m_QSESpeciesIndices.size());
for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) {
f_QSE(i) = dydt[m_QSESpeciesIndices[i]];
}
return 0; // Success
}
template <typename T>
int QSENetworkSolver::EigenFunctor<T>::df(const InputType& v_QSE, JacobianType& J_QSE) const {
std::vector<double> YFull(m_engine.getNetworkSpecies().size(), 0.0);
Eigen::VectorXd Y_QSE(v_QSE.array().exp());
for (size_t i = 0; i < m_dynamicSpeciesIndices.size(); ++i) {
YFull[m_dynamicSpeciesIndices[i]] = m_YDynamic[i];
}
int QSENetworkSolver::EigenFunctor<T>::df(const InputType& v_QSE_log, JacobianType& J_QSE) const {
std::vector<double> y = m_YFull;
Eigen::VectorXd Y_QSE = v_QSE_log.array().exp();
for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) {
YFull[m_QSESpeciesIndices[i]] = Y_QSE(i);
y[m_QSESpeciesIndices[i]] = Y_QSE(i);
}
m_engine.generateJacobianMatrix(YFull, m_T9, m_rho);
m_engine.generateJacobianMatrix(y, m_T9, m_rho);
Eigen::MatrixXd J_orig(m_QSESpeciesIndices.size(), m_QSESpeciesIndices.size());
J_QSE.resize(m_QSESpeciesIndices.size(), m_QSESpeciesIndices.size());
for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) {
for (size_t j = 0; j < m_QSESpeciesIndices.size(); ++j) {
J_orig(i, j) = m_engine.getJacobianMatrixEntry(m_QSESpeciesIndices[i], m_QSESpeciesIndices[j]);
J_QSE(i, j) = m_engine.getJacobianMatrixEntry(m_QSESpeciesIndices[i], m_QSESpeciesIndices[j]);
}
}
J_QSE = J_orig;
for (long j = 0; j < J_QSE.cols(); ++j) {
J_QSE.col(j) *= Y_QSE(j); // Chain rule for log space
J_QSE(j, j) *= Y_QSE(j); // chain rule for log space transformation
}
return 0; // Success
return 0;
}

15
src/network/include/gridfire/utils/logging.h Normal file
View File

@@ -0,0 +1,15 @@
#pragma once
#include "gridfire/engine/engine_abstract.h"
#include <string>
#include <vector>
namespace gridfire::utils {
std::string formatNuclearTimescaleLogString(
const DynamicEngine& engine,
const std::vector<double>& Y,
const double T9,
const double rho
);
}

1
src/network/lib/engine/engine_approx8.cpp
View File

@@ -509,7 +509,6 @@ namespace gridfire::approx8{
vector_type Approx8Network::convert_netIn(const NetIn &netIn) {
vector_type y(Approx8Net::nVar, 0.0);
y[Approx8Net::ih1] = netIn.composition.getNumberFraction("H-1");
std::cout << "Approx8::convert_netIn -> H-1 fraction: " << y[Approx8Net::ih1] << std::endl;
y[Approx8Net::ihe3] = netIn.composition.getNumberFraction("He-3");
y[Approx8Net::ihe4] = netIn.composition.getNumberFraction("He-4");
y[Approx8Net::ic12] = netIn.composition.getNumberFraction("C-12");

14
src/network/lib/engine/engine_graph.cpp
View File

@@ -1,6 +1,7 @@
#include "gridfire/engine/engine_graph.h"
#include "gridfire/reaction/reaction.h"
#include "gridfire/network.h"
#include "gridfire/screening/screening_types.h"
#include "fourdst/composition/species.h"
#include "fourdst/composition/atomicSpecies.h"
@@ -262,6 +263,15 @@ namespace gridfire {
return calculateAllDerivatives<ADDouble>(Y_in, T9, rho);
}
void GraphEngine::setScreeningModel(const screening::ScreeningType model) {
m_screeningModel = screening::selectScreeningModel(model);
m_screeningType = model;
}
screening::ScreeningType GraphEngine::getScreeningModel() const {
return m_screeningType;
}
double GraphEngine::calculateMolarReactionFlow(
const reaction::Reaction &reaction,
const std::vector<double> &Y,
@@ -439,6 +449,10 @@ namespace gridfire {
return speciesTimescales;
}
void GraphEngine::update(const NetIn &netIn) {
return; // No-op for GraphEngine, as it does not support manually triggering updates.
}
void GraphEngine::recordADTape() {
LOG_TRACE_L1(m_logger, "Recording AD tape for the RHS calculation...");

10
src/network/lib/engine/views/engine_adaptive.cpp
View File

@@ -213,6 +213,14 @@ namespace gridfire {
}
void AdaptiveEngineView::setScreeningModel(const screening::ScreeningType model) {
m_baseEngine.setScreeningModel(model);
}
screening::ScreeningType AdaptiveEngineView::getScreeningModel() const {
return m_baseEngine.getScreeningModel();
}
std::vector<double> AdaptiveEngineView::mapCulledToFull(const std::vector<double>& culled) const {
std::vector<double> full(m_baseEngine.getNetworkSpecies().size(), 0.0);
for (size_t i_culled = 0; i_culled < culled.size(); ++i_culled) {
@@ -338,7 +346,7 @@ namespace gridfire {
const double maxFlow
) const {
LOG_TRACE_L1(m_logger, "Culling reactions based on flow rates...");
const double relative_culling_threshold = m_config.get<double>("gridfire:AdaptiveEngineView:RelativeCullingThreshold", 1e-75);
const auto relative_culling_threshold = m_config.get<double>("gridfire:AdaptiveEngineView:RelativeCullingThreshold", 1e-75);
double absoluteCullingThreshold = relative_culling_threshold * maxFlow;
LOG_DEBUG(m_logger, "Relative culling threshold: {:0.3E} ({})", relative_culling_threshold, absoluteCullingThreshold);
std::vector<const reaction::LogicalReaction*> culledReactions;

313
src/network/lib/engine/views/engine_defined.cpp
View File

@@ -1,3 +1,310 @@
//
// Created by Emily Boudreaux on 6/30/25.
//
#include "gridfire/engine/views/engine_defined.h"
#include "quill/LogMacros.h"
namespace gridfire {
using fourdst::atomic::Species;
FileDefinedEngineView::FileDefinedEngineView(
DynamicEngine &baseEngine,
const std::string &fileName,
const io::NetworkFileParser &parser
):
m_baseEngine(baseEngine),
m_fileName(fileName),
m_parser(parser),
m_activeSpecies(baseEngine.getNetworkSpecies()),
m_activeReactions(baseEngine.getNetworkReactions()) {
buildFromFile(fileName);
}
const DynamicEngine & FileDefinedEngineView::getBaseEngine() const {
return m_baseEngine;
}
const std::vector<Species> & FileDefinedEngineView::getNetworkSpecies() const {
return m_activeSpecies;
}
StepDerivatives<double> FileDefinedEngineView::calculateRHSAndEnergy(
const std::vector<double> &Y_defined,
const double T9,
const double rho
) const {
validateNetworkState();
const auto Y_full = mapViewToFull(Y_defined);
const auto [dydt, nuclearEnergyGenerationRate] = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho);
StepDerivatives<double> definedResults;
definedResults.nuclearEnergyGenerationRate = nuclearEnergyGenerationRate;
definedResults.dydt = mapFullToView(dydt);
return definedResults;
}
void FileDefinedEngineView::generateJacobianMatrix(
const std::vector<double> &Y_defined,
const double T9,
const double rho
) {
validateNetworkState();
const auto Y_full = mapViewToFull(Y_defined);
m_baseEngine.generateJacobianMatrix(Y_full, T9, rho);
}
double FileDefinedEngineView::getJacobianMatrixEntry(
const int i_defined,
const int j_defined
) const {
validateNetworkState();
const size_t i_full = mapViewToFullSpeciesIndex(i_defined);
const size_t j_full = mapViewToFullSpeciesIndex(j_defined);
return m_baseEngine.getJacobianMatrixEntry(i_full, j_full);
}
void FileDefinedEngineView::generateStoichiometryMatrix() {
validateNetworkState();
m_baseEngine.generateStoichiometryMatrix();
}
int FileDefinedEngineView::getStoichiometryMatrixEntry(
const int speciesIndex_defined,
const int reactionIndex_defined
) const {
validateNetworkState();
const size_t i_full = mapViewToFullSpeciesIndex(speciesIndex_defined);
const size_t j_full = mapViewToFullReactionIndex(reactionIndex_defined);
return m_baseEngine.getStoichiometryMatrixEntry(i_full, j_full);
}
double FileDefinedEngineView::calculateMolarReactionFlow(
const reaction::Reaction &reaction,
const std::vector<double> &Y_defined,
const double T9,
const double rho
) const {
validateNetworkState();
if (!m_activeReactions.contains(reaction)) {
LOG_ERROR(m_logger, "Reaction '{}' is not part of the active reactions in the file defined engine view.", reaction.id());
m_logger -> flush_log();
throw std::runtime_error("Reaction not found in active reactions: " + std::string(reaction.id()));
}
const auto Y_full = mapViewToFull(Y_defined);
return m_baseEngine.calculateMolarReactionFlow(reaction, Y_full, T9, rho);
}
const reaction::LogicalReactionSet & FileDefinedEngineView::getNetworkReactions() const {
validateNetworkState();
return m_activeReactions;
}
std::unordered_map<Species, double> FileDefinedEngineView::getSpeciesTimescales(
const std::vector<double> &Y_defined,
const double T9,
const double rho
) const {
validateNetworkState();
const auto Y_full = mapViewToFull(Y_defined);
const auto fullTimescales = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho);
std::unordered_map<Species, double> definedTimescales;
for (const auto& active_species : m_activeSpecies) {
if (fullTimescales.contains(active_species)) {
definedTimescales[active_species] = fullTimescales.at(active_species);
}
}
return definedTimescales;
}
void FileDefinedEngineView::update(const NetIn &netIn) {
if (m_isStale) {
buildFromFile(m_fileName);
}
}
void FileDefinedEngineView::setNetworkFile(const std::string &fileName) {
m_isStale = true;
m_fileName = fileName;
LOG_DEBUG(m_logger, "File '{}' set to '{}'. FileDefinedNetworkView is now stale! You MUST call update() before you use it!", m_fileName, fileName);
}
void FileDefinedEngineView::setScreeningModel(const screening::ScreeningType model) {
m_baseEngine.setScreeningModel(model);
}
screening::ScreeningType FileDefinedEngineView::getScreeningModel() const {
return m_baseEngine.getScreeningModel();
}
std::vector<size_t> FileDefinedEngineView::constructSpeciesIndexMap() const {
LOG_TRACE_L1(m_logger, "Constructing species index map for file defined engine view...");
std::unordered_map<Species, size_t> fullSpeciesReverseMap;
const auto& fullSpeciesList = m_baseEngine.getNetworkSpecies();
fullSpeciesReverseMap.reserve(fullSpeciesList.size());
for (size_t i = 0; i < fullSpeciesList.size(); ++i) {
fullSpeciesReverseMap[fullSpeciesList[i]] = i;
}
std::vector<size_t> speciesIndexMap;
speciesIndexMap.reserve(m_activeSpecies.size());
for (const auto& active_species : m_activeSpecies) {
auto it = fullSpeciesReverseMap.find(active_species);
if (it != fullSpeciesReverseMap.end()) {
speciesIndexMap.push_back(it->second);
} else {
LOG_ERROR(m_logger, "Species '{}' not found in full species map.", active_species.name());
m_logger -> flush_log();
throw std::runtime_error("Species not found in full species map: " + std::string(active_species.name()));
}
}
LOG_TRACE_L1(m_logger, "Species index map constructed with {} entries.", speciesIndexMap.size());
return speciesIndexMap;
}
std::vector<size_t> FileDefinedEngineView::constructReactionIndexMap() const {
LOG_TRACE_L1(m_logger, "Constructing reaction index map for file defined engine view...");
// --- Step 1: Create a reverse map using the reaction's unique ID as the key. ---
std::unordered_map<std::string_view, size_t> fullReactionReverseMap;
const auto& fullReactionSet = m_baseEngine.getNetworkReactions();
fullReactionReverseMap.reserve(fullReactionSet.size());
for (size_t i_full = 0; i_full < fullReactionSet.size(); ++i_full) {
fullReactionReverseMap[fullReactionSet[i_full].id()] = i_full;
}
// --- Step 2: Build the final index map using the active reaction set. ---
std::vector<size_t> reactionIndexMap;
reactionIndexMap.reserve(m_activeReactions.size());
for (const auto& active_reaction_ptr : m_activeReactions) {
auto it = fullReactionReverseMap.find(active_reaction_ptr.id());
if (it != fullReactionReverseMap.end()) {
reactionIndexMap.push_back(it->second);
} else {
LOG_ERROR(m_logger, "Active reaction '{}' not found in base engine during reaction index map construction.", active_reaction_ptr.id());
m_logger->flush_log();
throw std::runtime_error("Mismatch between active reactions and base engine.");
}
}
LOG_TRACE_L1(m_logger, "Reaction index map constructed with {} entries.", reactionIndexMap.size());
return reactionIndexMap;
}
void FileDefinedEngineView::buildFromFile(const std::string &fileName) {
LOG_TRACE_L1(m_logger, "Building file defined engine view from {}...", fileName);
auto [reactionPENames] = m_parser.parse(fileName);
m_activeReactions.clear();
m_activeSpecies.clear();
std::unordered_set<Species> seenSpecies;
const auto& fullNetworkReactionSet = m_baseEngine.getNetworkReactions();
for (const auto& peName : reactionPENames) {
if (!fullNetworkReactionSet.contains(peName)) {
LOG_ERROR(m_logger, "Reaction with name '{}' not found in the base engine's network reactions. Aborting...", peName);
m_logger->flush_log();
throw std::runtime_error("Reaction with name '" + std::string(peName) + "' not found in the base engine's network reactions.");
}
auto reaction = fullNetworkReactionSet[peName];
for (const auto& reactant : reaction.reactants()) {
if (!seenSpecies.contains(reactant)) {
seenSpecies.insert(reactant);
m_activeSpecies.push_back(reactant);
}
}
for (const auto& product : reaction.products()) {
if (!seenSpecies.contains(product)) {
seenSpecies.insert(product);
m_activeSpecies.push_back(product);
}
}
m_activeReactions.add_reaction(reaction);
}
LOG_TRACE_L1(m_logger, "File defined engine view built with {} active species and {} active reactions.", m_activeSpecies.size(), m_activeReactions.size());
LOG_DEBUG(m_logger, "Active species: {}", [this]() -> std::string {
std::string result;
for (const auto& species : m_activeSpecies) {
result += std::string(species.name()) + ", ";
}
if (!result.empty()) {
result.pop_back(); // Remove last space
result.pop_back(); // Remove last comma
}
return result;
}());
LOG_DEBUG(m_logger, "Active reactions: {}", [this]() -> std::string {
std::string result;
for (const auto& reaction : m_activeReactions) {
result += std::string(reaction.id()) + ", ";
}
if (!result.empty()) {
result.pop_back(); // Remove last space
result.pop_back(); // Remove last comma
}
return result;
}());
m_speciesIndexMap = constructSpeciesIndexMap();
m_reactionIndexMap = constructReactionIndexMap();
m_isStale = false;
}
std::vector<double> FileDefinedEngineView::mapViewToFull(const std::vector<double>& culled) const {
std::vector<double> full(m_baseEngine.getNetworkSpecies().size(), 0.0);
for (size_t i_culled = 0; i_culled < culled.size(); ++i_culled) {
const size_t i_full = m_speciesIndexMap[i_culled];
full[i_full] += culled[i_culled];
}
return full;
}
std::vector<double> FileDefinedEngineView::mapFullToView(const std::vector<double>& full) const {
std::vector<double> culled(m_activeSpecies.size(), 0.0);
for (size_t i_culled = 0; i_culled < m_activeSpecies.size(); ++i_culled) {
const size_t i_full = m_speciesIndexMap[i_culled];
culled[i_culled] = full[i_full];
}
return culled;
}
size_t FileDefinedEngineView::mapViewToFullSpeciesIndex(size_t culledSpeciesIndex) const {
if (culledSpeciesIndex < 0 || culledSpeciesIndex >= static_cast<int>(m_speciesIndexMap.size())) {
LOG_ERROR(m_logger, "Defined index {} is out of bounds for species index map of size {}.", culledSpeciesIndex, m_speciesIndexMap.size());
m_logger->flush_log();
throw std::out_of_range("Defined index " + std::to_string(culledSpeciesIndex) + " is out of bounds for species index map of size " + std::to_string(m_speciesIndexMap.size()) + ".");
}
return m_speciesIndexMap[culledSpeciesIndex];
}
size_t FileDefinedEngineView::mapViewToFullReactionIndex(size_t culledReactionIndex) const {
if (culledReactionIndex < 0 || culledReactionIndex >= static_cast<int>(m_reactionIndexMap.size())) {
LOG_ERROR(m_logger, "Defined index {} is out of bounds for reaction index map of size {}.", culledReactionIndex, m_reactionIndexMap.size());
m_logger->flush_log();
throw std::out_of_range("Defined index " + std::to_string(culledReactionIndex) + " is out of bounds for reaction index map of size " + std::to_string(m_reactionIndexMap.size()) + ".");
}
return m_reactionIndexMap[culledReactionIndex];
}
void FileDefinedEngineView::validateNetworkState() const {
if (m_isStale) {
LOG_ERROR(m_logger, "File defined engine view is stale. Please call update() with a valid NetIn object.");
m_logger->flush_log();
throw std::runtime_error("File defined engine view is stale. Please call update() with a valid NetIn object.");
}
}
}

80
src/network/lib/io/network_file.cpp
View File

@@ -1,3 +1,77 @@
//
// Created by Emily Boudreaux on 6/30/25.
//
#include "gridfire/io/network_file.h"
#include <string>
#include <vector>
#include <algorithm>
#include <fstream>
#include <stdexcept>
#include "quill/LogMacros.h"
namespace gridfire::io {
namespace {
inline void ltrim(std::string &s) {
s.erase(
s.begin(),
std::ranges::find_if(s,
[](const unsigned char ch) {
return !std::isspace(ch);
})
);
}
inline void rtrim(std::string &s) {
s.erase(
std::find_if(
s.rbegin(),
s.rend(),
[](const unsigned char ch) {
return !std::isspace(ch);
}).base(),
s.end()
);
}
inline void trim(std::string &s) {
ltrim(s);
rtrim(s);
}
}
SimpleReactionListFileParser::SimpleReactionListFileParser() {}
ParsedNetworkData SimpleReactionListFileParser::parse(const std::string& filename) const {
LOG_TRACE_L1(m_logger, "Parsing simple reaction list file: {}", filename);
std::ifstream file(filename);
if (!file.is_open()) {
LOG_ERROR(m_logger, "Failed to open file: {}", filename);
m_logger -> flush_log();
throw std::runtime_error("Could not open file: " + filename);
}
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);
const size_t comment_pos = line.find('#');
if (comment_pos != std::string::npos) {
line = line.substr(0, comment_pos);
}
trim(line);
if (line.empty()) {
continue; // Skip empty lines
}
parsed.reactionPENames.push_back(line);
}
LOG_TRACE_L1(m_logger, "Parsed {} reactions from file: {}", parsed.reactionPENames.size(), filename);
return parsed;
}
}

2
src/network/lib/network.cpp
View File

@@ -84,7 +84,7 @@ namespace gridfire {
}
}
const ReactionSet reactionSet(reaclibReactions);
return LogicalReactionSet(reactionSet);
return packReactionSetToLogicalReactionSet(reactionSet);
}
// Trim whitespace from both ends of a string

6
src/network/lib/reaction/reaclib.cpp
View File

@@ -122,10 +122,12 @@ namespace gridfire::reaclib {
}
// The ReactionSet takes the vector of all individual reactions.
reaction::ReactionSet reaction_set(std::move(reaction_list));
const reaction::ReactionSet reaction_set(std::move(reaction_list));
// The LogicalReactionSet groups reactions by their peName, which is what we want.
s_all_reaclib_reactions_ptr = new reaction::LogicalReactionSet(reaction_set);
s_all_reaclib_reactions_ptr = new reaction::LogicalReactionSet(
reaction::packReactionSetToLogicalReactionSet(reaction_set)
);
s_initialized = true;
}

177
src/network/lib/reaction/reaction.cpp
View File

@@ -137,152 +137,6 @@ namespace gridfire::reaction {
return XXHash64::hash(m_id.data(), m_id.size(), seed);
}
ReactionSet::ReactionSet(
std::vector<Reaction> reactions
) :
m_reactions(std::move(reactions)) {
if (m_reactions.empty()) {
return; // Case where the reactions will be added later.
}
m_reactionNameMap.reserve(reactions.size());
for (const auto& reaction : m_reactions) {
m_id += reaction.id();
m_reactionNameMap.emplace(reaction.id(), reaction);
}
}
ReactionSet::ReactionSet(const ReactionSet &other) {
m_reactions.reserve(other.m_reactions.size());
for (const auto& reaction_ptr: other.m_reactions) {
m_reactions.push_back(reaction_ptr);
}
m_reactionNameMap.reserve(other.m_reactionNameMap.size());
for (const auto& reaction_ptr : m_reactions) {
m_reactionNameMap.emplace(reaction_ptr.id(), reaction_ptr);
}
}
ReactionSet & ReactionSet::operator=(const ReactionSet &other) {
if (this != &other) {
ReactionSet temp(other);
std::swap(m_reactions, temp.m_reactions);
std::swap(m_reactionNameMap, temp.m_reactionNameMap);
}
return *this;
}
void ReactionSet::add_reaction(Reaction reaction) {
m_reactions.emplace_back(reaction);
m_id += m_reactions.back().id();
m_reactionNameMap.emplace(m_reactions.back().id(), m_reactions.back());
}
void ReactionSet::remove_reaction(const Reaction& reaction) {
if (!m_reactionNameMap.contains(std::string(reaction.id()))) {
return;
}
m_reactionNameMap.erase(std::string(reaction.id()));
std::erase_if(m_reactions, [&reaction](const Reaction& r) {
return r == reaction;
});
}
bool ReactionSet::contains(const std::string_view& id) const {
for (const auto& reaction : m_reactions) {
if (reaction.id() == id) {
return true;
}
}
return false;
}
bool ReactionSet::contains(const Reaction& reaction) const {
for (const auto& r : m_reactions) {
if (r == reaction) {
return true;
}
}
return false;
}
void ReactionSet::clear() {
m_reactions.clear();
m_reactionNameMap.clear();
}
bool ReactionSet::contains_species(const Species& species) const {
for (const auto& reaction : m_reactions) {
if (reaction.contains(species)) {
return true;
}
}
return false;
}
bool ReactionSet::contains_reactant(const Species& species) const {
for (const auto& r : m_reactions) {
if (r.contains_reactant(species)) {
return true;
}
}
return false;
}
bool ReactionSet::contains_product(const Species& species) const {
for (const auto& r : m_reactions) {
if (r.contains_product(species)) {
return true;
}
}
return false;
}
const Reaction& ReactionSet::operator[](const size_t index) const {
if (index >= m_reactions.size()) {
m_logger -> flush_log();
throw std::out_of_range("Index" + std::to_string(index) + " out of range for ReactionSet of size " + std::to_string(m_reactions.size()) + ".");
}
return m_reactions[index];
}
const Reaction& ReactionSet::operator[](const std::string_view& id) const {
if (auto it = m_reactionNameMap.find(std::string(id)); it != m_reactionNameMap.end()) {
return it->second;
}
m_logger -> flush_log();
throw std::out_of_range("Species " + std::string(id) + " does not exist in ReactionSet.");
}
bool ReactionSet::operator==(const ReactionSet& other) const {
if (size() != other.size()) {
return false;
}
return hash() == other.hash();
}
bool ReactionSet::operator!=(const ReactionSet& other) const {
return !(*this == other);
}
uint64_t ReactionSet::hash(uint64_t seed) const {
if (m_reactions.empty()) {
return XXHash64::hash(nullptr, 0, seed);
}
std::vector<uint64_t> individualReactionHashes;
individualReactionHashes.reserve(m_reactions.size());
for (const auto& reaction : m_reactions) {
individualReactionHashes.push_back(reaction.hash(seed));
}
std::ranges::sort(individualReactionHashes);
const void* data = static_cast<const void*>(individualReactionHashes.data());
size_t sizeInBytes = individualReactionHashes.size() * sizeof(uint64_t);
return XXHash64::hash(data, sizeInBytes, seed);
}
LogicalReaction::LogicalReaction(const std::vector<Reaction>& reactants) :
@@ -344,21 +198,21 @@ namespace gridfire::reaction {
return calculate_rate<CppAD::AD<double>>(T9);
}
LogicalReactionSet::LogicalReactionSet(const ReactionSet &reactionSet) :
ReactionSet(std::vector<Reaction>()) {
LogicalReactionSet packReactionSetToLogicalReactionSet(const ReactionSet& reactionSet) {
std::unordered_map<std::string_view, std::vector<Reaction>> groupedReactions;
std::unordered_map<std::string_view, std::vector<Reaction>> grouped_reactions;
for (const auto& reaction: reactionSet) {
groupedReactions[reaction.peName()].push_back(reaction);
}
for (const auto& reaction : reactionSet) {
grouped_reactions[reaction.peName()].push_back(reaction);
}
m_reactions.reserve(grouped_reactions.size());
m_reactionNameMap.reserve(grouped_reactions.size());
for (const auto &reactions_for_peName: grouped_reactions | std::views::values) {
LogicalReaction logical_reaction(reactions_for_peName);
m_reactionNameMap.emplace(logical_reaction.id(), logical_reaction);
m_reactions.push_back(std::move(logical_reaction));
std::vector<LogicalReaction> reactions;
reactions.reserve(groupedReactions.size());
for (const auto &reactionsGroup: groupedReactions | std::views::values) {
LogicalReaction logicalReaction(reactionsGroup);
reactions.push_back(logicalReaction);
}
return LogicalReactionSet(std::move(reactions));
}
}
@@ -376,4 +230,11 @@ namespace std {
return s.hash(0);
}
};
template<>
struct hash<gridfire::reaction::LogicalReactionSet> {
size_t operator()(const gridfire::reaction::LogicalReactionSet& s) const noexcept {
return s.hash(0);
}
};
} // namespace std

31
src/network/lib/screening/screening_bare.cpp Normal file
View File

@@ -0,0 +1,31 @@
#include "gridfire/screening/screening_bare.h"
#include "fourdst/composition/atomicSpecies.h"
#include "cppad/cppad.hpp"
#include <vector>
namespace gridfire::screening {
using ADDouble = CppAD::AD<double>;
std::vector<ADDouble> BareScreeningModel::calculateScreeningFactors(
const reaction::LogicalReactionSet &reactions,
const std::vector<fourdst::atomic::Species>& species,
const std::vector<ADDouble> &Y,
const ADDouble T9,
const ADDouble rho
) const {
return calculateFactors_impl<ADDouble>(reactions, species, Y, T9, rho);
}
std::vector<double> BareScreeningModel::calculateScreeningFactors(
const reaction::LogicalReactionSet &reactions,
const std::vector<fourdst::atomic::Species>& species,
const std::vector<double> &Y,
const double T9,
const double rho
) const {
return calculateFactors_impl<double>(reactions, species, Y, T9, rho);
}
}

19
src/network/lib/screening/screening_types.cpp Normal file
View File

@@ -0,0 +1,19 @@
#include "gridfire/screening/screening_abstract.h"
#include "gridfire/screening/screening_types.h"
#include "gridfire/screening/screening_weak.h"
#include "gridfire/screening/screening_bare.h"
#include <memory>
namespace gridfire::screening {
std::unique_ptr<ScreeningModel> selectScreeningModel(const ScreeningType type) {
switch (type) {
case ScreeningType::WEAK:
return std::make_unique<WeakScreeningModel>();
case ScreeningType::BARE:
return std::make_unique<BareScreeningModel>();
default:
return std::make_unique<BareScreeningModel>();
}
}
}

31
src/network/lib/screening/screening_weak.cpp Normal file
View File

@@ -0,0 +1,31 @@
#include "gridfire/screening/screening_weak.h"
#include "fourdst/composition/atomicSpecies.h"
#include "cppad/cppad.hpp"
#include <vector>
namespace gridfire::screening {
using ADDouble = CppAD::AD<double>;
std::vector<ADDouble> WeakScreeningModel::calculateScreeningFactors(
const reaction::LogicalReactionSet &reactions,
const std::vector<fourdst::atomic::Species>& species,
const std::vector<ADDouble> &Y,
const ADDouble T9,
const ADDouble rho
) const {
return calculateFactors_impl<ADDouble>(reactions, species, Y, T9, rho);
}
std::vector<double> WeakScreeningModel::calculateScreeningFactors(
const reaction::LogicalReactionSet &reactions,
const std::vector<fourdst::atomic::Species>& species,
const std::vector<double> &Y,
const double T9,
const double rho
) const {
return calculateFactors_impl<double>(reactions, species, Y, T9, rho);
}
}

89
src/network/lib/solver/solver.cpp
View File

@@ -2,6 +2,8 @@
#include "gridfire/engine/engine_graph.h"
#include "gridfire/network.h"
#include "gridfire/utils/logging.h"
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/composition/composition.h"
#include "fourdst/config/config.h"
@@ -15,6 +17,7 @@
#include <unordered_map>
#include <string>
#include <stdexcept>
#include <iomanip>
#include "quill/LogMacros.h"
@@ -59,6 +62,21 @@ namespace gridfire::solver {
Eigen::VectorXd Y_QSE;
try {
Y_QSE = calculateSteadyStateAbundances(Y_sanitized_initial, T9, rho, indices);
LOG_DEBUG(m_logger, "QSE Abundances: {}", [*this](const dynamicQSESpeciesIndices& indices, const Eigen::VectorXd& Y_QSE) -> std::string {
std::stringstream ss;
ss << std::scientific << std::setprecision(5);
for (size_t i = 0; i < indices.QSESpeciesIndices.size(); ++i) {
ss << std::string(m_engine.getNetworkSpecies()[indices.QSESpeciesIndices[i]].name()) + ": ";
ss << Y_QSE(i);
if (i < indices.QSESpeciesIndices.size() - 2) {
ss << ", ";
} else if (i == indices.QSESpeciesIndices.size() - 2) {
ss << ", and ";
}
}
return ss.str();
}(indices, Y_QSE));
} catch (const std::runtime_error& e) {
LOG_ERROR(m_logger, "Failed to calculate steady state abundances. Aborting QSE evaluation.");
m_logger->flush_log();
@@ -191,12 +209,56 @@ namespace gridfire::solver {
const dynamicQSESpeciesIndices &indices
) const {
LOG_TRACE_L1(m_logger, "Calculating steady state abundances for QSE species...");
LOG_WARNING(m_logger, "QSE solver logic not yet implemented, assuming all QSE species have 0 abundance.");
// --- Prepare the QSE species vector ---
Eigen::VectorXd v_qse(indices.QSESpeciesIndices.size());
v_qse.setZero();
return v_qse.array();
if (indices.QSESpeciesIndices.empty()) {
LOG_DEBUG(m_logger, "No QSE species to solve for.");
return Eigen::VectorXd(0);
}
// Use the EigenFunctor with Eigen's nonlinear solver
EigenFunctor<double> functor(
m_engine,
Y,
indices.dynamicSpeciesIndices,
indices.QSESpeciesIndices,
T9,
rho
);
Eigen::VectorXd v_qse_log_initial(indices.QSESpeciesIndices.size());
for (size_t i = 0; i < indices.QSESpeciesIndices.size(); ++i) {
v_qse_log_initial(i) = std::log(std::max(Y[indices.QSESpeciesIndices[i]], 1e-99));
}
const static std::unordered_map<Eigen::LevenbergMarquardtSpace::Status, const char*> statusMessages = {
{Eigen::LevenbergMarquardtSpace::NotStarted, "Not started"},
{Eigen::LevenbergMarquardtSpace::Running, "Running"},
{Eigen::LevenbergMarquardtSpace::ImproperInputParameters, "Improper input parameters"},
{Eigen::LevenbergMarquardtSpace::RelativeReductionTooSmall, "Relative reduction too small"},
{Eigen::LevenbergMarquardtSpace::RelativeErrorTooSmall, "Relative error too small"},
{Eigen::LevenbergMarquardtSpace::RelativeErrorAndReductionTooSmall, "Relative error and reduction too small"},
{Eigen::LevenbergMarquardtSpace::CosinusTooSmall, "Cosine too small"},
{Eigen::LevenbergMarquardtSpace::TooManyFunctionEvaluation, "Too many function evaluations"},
{Eigen::LevenbergMarquardtSpace::FtolTooSmall, "Function tolerance too small"},
{Eigen::LevenbergMarquardtSpace::XtolTooSmall, "X tolerance too small"},
{Eigen::LevenbergMarquardtSpace::GtolTooSmall, "Gradient tolerance too small"},
{Eigen::LevenbergMarquardtSpace::UserAsked, "User asked to stop"}
};
Eigen::LevenbergMarquardt lm(functor);
const Eigen::LevenbergMarquardtSpace::Status info = lm.minimize(v_qse_log_initial);
if (info <= 0 || info >= 4) {
LOG_ERROR(m_logger, "QSE species minimization failed with status: {} ({})",
static_cast<int>(info), statusMessages.at(info));
throw std::runtime_error(
"QSE species minimization failed with status: " + std::to_string(static_cast<int>(info)) +
" (" + std::string(statusMessages.at(info)) + ")"
);
}
LOG_DEBUG(m_logger, "QSE species minimization completed successfully with status: {} ({})",
static_cast<int>(info), statusMessages.at(info));
return v_qse_log_initial.array().exp();
}
NetOut QSENetworkSolver::initializeNetworkWithShortIgnition(const NetIn &netIn) const {
@@ -232,9 +294,14 @@ namespace gridfire::solver {
preIgnition.tMax = ignitionTime;
preIgnition.dt0 = ignitionStepSize;
const auto prevScreeningModel = m_engine.getScreeningModel();
LOG_DEBUG(m_logger, "Setting screening model to BARE for high temperature and density ignition.");
m_engine.setScreeningModel(screening::ScreeningType::BARE);
DirectNetworkSolver ignitionSolver(m_engine);
NetOut postIgnition = ignitionSolver.evaluate(preIgnition);
LOG_INFO(m_logger, "Network ignition completed in {} steps.", postIgnition.num_steps);
m_engine.setScreeningModel(prevScreeningModel);
LOG_DEBUG(m_logger, "Restoring previous screening model: {}", static_cast<int>(prevScreeningModel));
return postIgnition;
}
@@ -360,7 +427,8 @@ namespace gridfire::solver {
speciesNames.push_back(std::string(species.name()));
}
Composition outputComposition(speciesNames, finalMassFractions);
Composition outputComposition(speciesNames);
outputComposition.setMassFraction(speciesNames, finalMassFractions);
outputComposition.finalize(true);
NetOut netOut;
@@ -377,6 +445,15 @@ namespace gridfire::solver {
double t
) const {
const std::vector<double> y(Y.begin(), m_numSpecies + Y.begin());
// std::string timescales = utils::formatNuclearTimescaleLogString(
// m_engine,
// y,
// m_T9,
// m_rho
// );
// LOG_TRACE_L2(m_logger, "{}", timescales);
auto [dydt, eps] = m_engine.calculateRHSAndEnergy(y, m_T9, m_rho);
dYdt.resize(m_numSpecies + 1);
std::ranges::copy(dydt, dYdt.begin());

60
src/network/lib/utils/logging.cpp Normal file
View File

@@ -0,0 +1,60 @@
#include "gridfire/utils/logging.h"
#include "gridfire/engine/engine_abstract.h"
#include <sstream>
#include <iomanip>
#include <algorithm>
#include <ranges>
#include <string_view>
#include <string>
#include <iostream>
#include <vector>
std::string gridfire::utils::formatNuclearTimescaleLogString(
const DynamicEngine& engine,
std::vector<double> const& Y,
const double T9,
const double rho
) {
auto const& timescales = engine.getSpeciesTimescales(Y, T9, rho);
// Figure out how wide the "Species" column needs to be:
std::size_t maxNameLen = std::string_view("Species").size();
for (const auto &key: timescales | std::views::keys) {
std::string_view name = key.name();
maxNameLen = std::max(maxNameLen, name.size());
}
// Pick a fixed width for the timescale column:
constexpr int timescaleWidth = 12;
std::ostringstream ss;
ss << "== Timescales (s) ==\n";
// Header row
ss << std::left << std::setw(static_cast<int>(maxNameLen) + 2) << "Species"
<< std::right << std::setw(timescaleWidth) << "Timescale (s)" << "\n";
// Underline
ss << std::string(static_cast<int>(maxNameLen) + 2 + timescaleWidth, '=') << "\n";
ss << std::scientific;
// Data rows
for (auto const& [species, timescale] : timescales) {
const std::string_view name = species.name();
ss << std::left << std::setw(static_cast<int>(maxNameLen) + 2) << name;
if (std::isinf(timescale)) {
ss << std::right << std::setw(timescaleWidth) << "inf" << "\n";
} else {
ss << std::right << std::setw(timescaleWidth)
<< std::scientific << std::setprecision(3) << timescale << "\n";
}
}
// Footer underline
ss << std::string(static_cast<int>(maxNameLen) + 2 + timescaleWidth, '=') << "\n";
return ss.str();
}

19
src/network/meson.build
View File

@@ -3,10 +3,16 @@ network_sources = files(
'lib/network.cpp',
'lib/engine/engine_approx8.cpp',
'lib/engine/engine_graph.cpp',
'lib/engine/engine_adaptive.cpp',
'lib/engine/views/engine_adaptive.cpp',
'lib/engine/views/engine_defined.cpp',
'lib/reaction/reaction.cpp',
'lib/reaction/reaclib.cpp',
'lib/io/network_file.cpp',
'lib/solver/solver.cpp',
'lib/screening/screening_types.cpp',
'lib/screening/screening_weak.cpp',
'lib/screening/screening_bare.cpp',
'lib/utils/logging.cpp',
)
@@ -39,12 +45,19 @@ network_dep = declare_dependency(
network_headers = files(
'include/gridfire/network.h',
'include/gridfire/engine/engine_abstract.h',
'include/gridfire/engine/engine_view_abstract.h',
'include/gridfire/engine/views/engine_view_abstract.h',
'include/gridfire/engine/engine_approx8.h',
'include/gridfire/engine/engine_graph.h',
'include/gridfire/engine/engine_adaptive.h',
'include/gridfire/engine/views/engine_adaptive.h',
'include/gridfire/engine/views/engine_defined.h',
'include/gridfire/reaction/reaction.h',
'include/gridfire/reaction/reaclib.h',
'include/gridfire/io/network_file.h',
'include/gridfire/solver/solver.h',
'include/gridfire/screening/screening_abstract.h',
'include/gridfire/screening/screening_bare.h',
'include/gridfire/screening/screening_weak.h',
'include/gridfire/screening/screening_types.h',
'include/gridfire/utils/logging.h',
)
install_headers(network_headers, subdir : 'gridfire')