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

feat(Comoposition-Tracking): updated GridFire to use new, molar-abundance based, version of libcomposition (v2.0.6)

Browse Source 
This entailed a major rewrite of the composition handling from each engine and engine view along with the solver and primer. The intent here is to let Compositions be constructed from the same extensive property which the solver tracks internally. This addressed C0 discontinuity issues in the tracked molar abundances of species which were introduced by repeadidly swaping from molar abundance space to mass fraction space and back. This also allowed for a simplification of the primeNetwork method. Specifically the mass borrowing system was dramatically simplified as molar abundances are extensive.
This commit is contained in:
Emily Boudreaux
2025-11-10 10:40:03 -05:00
parent 534a44448b
commit a7a4a30028
57 changed files with 1878 additions and 2823 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
src/include/gridfire/engine/engine_abstract.h
View File

@@ -10,6 +10,8 @@
#include "gridfire/expectations/expected_engine.h"
#include "fourdst/composition/composition_abstract.h"
#include <vector>
#include <unordered_map>
#include <utility>
@@ -58,7 +60,7 @@ namespace gridfire {
*/
template <IsArithmeticOrAD T>
struct StepDerivatives {
std::map<fourdst::atomic::Species, T> dydt; ///< Derivatives of abundances (dY/dt for each species).
std::map<fourdst::atomic::Species, T> dydt{}; ///< Derivatives of abundances (dY/dt for each species).
T nuclearEnergyGenerationRate = T(0.0); ///< Specific energy generation rate (e.g., erg/g/s).
StepDerivatives() : dydt(), nuclearEnergyGenerationRate(T(0.0)) {}
@@ -120,7 +122,7 @@ namespace gridfire {
* rate for the current state.
*/
[[nodiscard]] virtual std::expected<StepDerivatives<double>, expectations::StaleEngineError> calculateRHSAndEnergy(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const = 0;
@@ -153,20 +155,20 @@ namespace gridfire {
* for the current state. The matrix can then be accessed via getJacobianMatrixEntry().
*/
virtual void generateJacobianMatrix(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const = 0;
virtual void generateJacobianMatrix(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho,
const std::vector<fourdst::atomic::Species>& activeSpecies
) const = 0;
virtual void generateJacobianMatrix(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho,
const SparsityPattern& sparsityPattern
@@ -223,7 +225,7 @@ namespace gridfire {
*/
[[nodiscard]] virtual double calculateMolarReactionFlow(
const reaction::Reaction& reaction,
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const = 0;
@@ -240,7 +242,7 @@ namespace gridfire {
* generation rate with respect to temperature and density for the current state.
*/
[[nodiscard]] virtual EnergyDerivatives calculateEpsDerivatives(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const = 0;
@@ -266,13 +268,13 @@ namespace gridfire {
* which can be used for timestep control, diagnostics, and reaction network culling.
*/
[[nodiscard]] virtual std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesTimescales(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const = 0;
[[nodiscard]] virtual std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesDestructionTimescales(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const = 0;
@@ -390,7 +392,7 @@ namespace gridfire {
* which may involve adding or removing species and reactions based on the
* specified depth. However, not all engines support this operation.
*/
virtual void rebuild(const fourdst::composition::Composition& comp, BuildDepthType depth) {
virtual void rebuild(const fourdst::composition::CompositionAbstract &comp, BuildDepthType depth) {
throw std::logic_error("Setting network depth not supported by this engine.");
// ReSharper disable once CppDFAUnreachableCode
}
@@ -408,7 +410,7 @@ namespace gridfire {
* example, by either QSE partitioning or reaction flow rate culling
*/
virtual fourdst::composition::Composition collectComposition(
fourdst::composition::Composition& comp
fourdst::composition::CompositionAbstract &comp
) const = 0;
};

331
src/include/gridfire/engine/engine_approx8.h
View File

@@ -1,331 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
#include <array>
#include <boost/numeric/odeint.hpp>
#include "gridfire/network.h"
/**
* @file approx8.h
* @brief Header file for the Approx8 nuclear reaction network.
*
* This file contains the definitions and declarations for the Approx8 nuclear reaction network.
* The network is based on Frank Timmes' "approx8" and includes 8 isotopes and various nuclear reactions.
* The rates are evaluated using a fitting function with coefficients from reaclib.jinaweb.org.
*/
namespace gridfire::approx8{
/**
* @typedef vector_type
* @brief Alias for a vector of doubles using Boost uBLAS.
*/
typedef boost::numeric::ublas::vector< double > vector_type;
/**
* @typedef matrix_type
* @brief Alias for a matrix of doubles using Boost uBLAS.
*/
typedef boost::numeric::ublas::matrix< double > matrix_type;
/**
* @typedef vec7
* @brief Alias for a std::array of 7 doubles.
*/
typedef std::array<double,7> vec7;
/**
* @struct Approx8Net
* @brief Contains constants and arrays related to the nuclear network.
*/
struct Approx8Net{
static constexpr int ih1=0;
static constexpr int ihe3=1;
static constexpr int ihe4=2;
static constexpr int ic12=3;
static constexpr int in14=4;
static constexpr int io16=5;
static constexpr int ine20=6;
static constexpr int img24=7;
static constexpr int iTemp=img24+1;
static constexpr int iDensity =iTemp+1;
static constexpr int iEnergy=iDensity+1;
static constexpr int nIso=img24+1; // number of isotopes
static constexpr int nVar=iEnergy+1; // number of variables
static constexpr std::array<int,nIso> aIon = {
1,
3,
4,
12,
14,
16,
20,
24
};
static constexpr std::array<double,nIso> mIon = {
1.67262164e-24,
5.00641157e-24,
6.64465545e-24,
1.99209977e-23,
2.32462686e-23,
2.65528858e-23,
3.31891077e-23,
3.98171594e-23
};
};
/**
* @brief Multiplies two arrays and sums the resulting elements.
* @param a First array.
* @param b Second array.
* @return Sum of the product of the arrays.
* @example
* @code
* vec7 a = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0};
* vec7 b = {0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5};
* double result = sum_product(a, b);
* @endcode
*/
double sum_product( const vec7 &a, const vec7 &b);
/**
* @brief Returns an array of T9 terms for the nuclear reaction rate fit.
* @param T Temperature in GigaKelvin.
* @return Array of T9 terms.
* @example
* @code
* double T = 1.5;
* vec7 T9_array = get_T9_array(T);
* @endcode
*/
vec7 get_T9_array(const double &T);
/**
* @brief Evaluates the nuclear reaction rate given the T9 array and coefficients.
* @param T9 Array of T9 terms.
* @param coef Array of coefficients.
* @return Evaluated rate.
* @example
* @code
* vec7 T9 = get_T9_array(1.5);
* vec7 coef = {1.0, 0.1, 0.01, 0.001, 0.0001, 0.00001, 0.000001};
* double rate = rate_fit(T9, coef);
* @endcode
*/
double rate_fit(const vec7 &T9, const vec7 &coef);
/**
* @brief Calculates the rate for the reaction p + p -> d.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double pp_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction p + d -> he3.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double dp_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction he3 + he3 -> he4 + 2p.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double he3he3_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction he3(he3,2p)he4.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double he3he4_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction he4 + he4 + he4 -> c12.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double triple_alpha_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction c12 + p -> n13.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double c12p_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction c12 + he4 -> o16.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double c12a_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction n14(p,g)o15 - o15 + p -> c12 + he4.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double n14p_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction n14(a,g)f18 assumed to go on to ne20.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double n14a_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction n15(p,a)c12 (CNO I).
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double n15pa_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction n15(p,g)o16 (CNO II).
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double n15pg_rate(const vec7 &T9);
/**
* @brief Calculates the fraction for the reaction n15(p,g)o16.
* @param T9 Array of T9 terms.
* @return Fraction of the reaction.
*/
double n15pg_frac(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction o16(p,g)f17 then f17 -> o17(p,a)n14.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double o16p_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction o16(a,g)ne20.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double o16a_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction ne20(a,g)mg24.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double ne20a_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction c12(c12,a)ne20.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double c12c12_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction c12(o16,a)mg24.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double c12o16_rate(const vec7 &T9);
/**
* @struct Jacobian
* @brief Functor to calculate the Jacobian matrix for implicit solvers.
*/
struct Jacobian {
/**
* @brief Calculates the Jacobian matrix.
* @param y State vector.
* @param J Jacobian matrix.
* @param dfdt Derivative of the state vector.
*/
void operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) const;
};
/**
* @struct ODE
* @brief Functor to calculate the derivatives for the ODE solver.
*/
struct ODE {
/**
* @brief Calculates the derivatives of the state vector.
* @param y State vector.
* @param dydt Derivative of the state vector.
*/
void operator() ( const vector_type &y, vector_type &dydt, double /* t */) const;
};
/**
* @class Approx8Network
* @brief Class for the Approx8 nuclear reaction network.
*/
class Approx8Network final : public Network {
public:
Approx8Network();
/**
* @brief Evaluates the nuclear network.
* @param netIn Input parameters for the network.
* @return Output results from the network.
*/
NetOut evaluate(const NetIn &netIn) override;
/**
* @brief Sets whether the solver should use a stiff method.
* @param stiff Boolean indicating if a stiff method should be used.
*/
void setStiff(bool stiff) override;
/**
* @brief Checks if the solver is using a stiff method.
* @return Boolean indicating if a stiff method is being used.
*/
bool isStiff() const override { return m_stiff; }
private:
vector_type m_y;
double m_tMax = 0;
double m_dt0 = 0;
bool m_stiff = false;
/**
* @brief Converts the input parameters to the internal state vector.
* @param netIn Input parameters for the network.
* @return Internal state vector.
*/
static vector_type convert_netIn(const NetIn &netIn);
};
} // namespace nnApprox8

56
src/include/gridfire/engine/engine_graph.h
View File

@@ -1,6 +1,6 @@
#pragma once
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/atomic/atomicSpecies.h"
#include "fourdst/composition/composition.h"
#include "fourdst/logging/logging.h"
#include "fourdst/config/config.h"
@@ -150,7 +150,7 @@ namespace gridfire {
* @see StepDerivatives
*/
[[nodiscard]] std::expected<StepDerivatives<double>, expectations::StaleEngineError> calculateRHSAndEnergy(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -172,7 +172,7 @@ namespace gridfire {
* @see StepDerivatives
*/
[[nodiscard]] std::expected<StepDerivatives<double>, expectations::StaleEngineError> calculateRHSAndEnergy(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract& comp,
double T9,
double rho,
const reaction::ReactionSet &activeReactions
@@ -193,7 +193,7 @@ namespace gridfire {
* @see EnergyDerivatives
*/
[[nodiscard]] EnergyDerivatives calculateEpsDerivatives(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -217,7 +217,7 @@ namespace gridfire {
* @see EnergyDerivatives
*/
[[nodiscard]] EnergyDerivatives calculateEpsDerivatives(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho,
const reaction::ReactionSet &activeReactions
@@ -237,7 +237,7 @@ namespace gridfire {
* @see getJacobianMatrixEntry()
*/
void generateJacobianMatrix(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -255,7 +255,7 @@ namespace gridfire {
* @see generateJacobianMatrix()
*/
void generateJacobianMatrix(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho,
const std::vector<fourdst::atomic::Species>& activeSpecies
@@ -277,7 +277,7 @@ namespace gridfire {
* @see getJacobianMatrixEntry()
*/
void generateJacobianMatrix(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho,
const SparsityPattern &sparsityPattern
@@ -306,7 +306,7 @@ namespace gridfire {
*/
[[nodiscard]] double calculateMolarReactionFlow(
const reaction::Reaction& reaction,
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -388,7 +388,7 @@ namespace gridfire {
* which can be used for timestep control or diagnostics.
*/
[[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesTimescales(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -407,7 +407,7 @@ namespace gridfire {
* calculations with different reaction sets without modifying the engine's internal state.
*/
[[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesTimescales(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho,
const reaction::ReactionSet &activeReactions
@@ -426,7 +426,7 @@ namespace gridfire {
* which can be useful for understanding reaction flows and equilibrium states.
*/
[[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesDestructionTimescales(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -445,7 +445,7 @@ namespace gridfire {
* calculations with different reaction sets without modifying the engine's internal state.
*/
[[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesDestructionTimescales(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho,
const reaction::ReactionSet &activeReactions
@@ -606,7 +606,7 @@ namespace gridfire {
const reaction::Reaction &reaction,
double T9,
double rho,
const fourdst::composition::Composition& comp
const fourdst::composition::CompositionAbstract &comp
) const;
/**
@@ -733,7 +733,7 @@ namespace gridfire {
* and build depth. It updates all internal data structures accordingly.
*/
void rebuild(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
BuildDepthType depth
) override;
@@ -748,25 +748,25 @@ namespace gridfire {
* have a molar abundance set to 0.
* @throws BadCollectionError If the input composition contains species not present in the network species set
*/
fourdst::composition::Composition collectComposition(fourdst::composition::Composition &comp) const override;
fourdst::composition::Composition collectComposition(fourdst::composition::CompositionAbstract &comp) const override;
private:
struct PrecomputedReaction {
// Forward cacheing
size_t reaction_index;
reaction::ReactionType reaction_type;
uint64_t reaction_hash;
std::vector<size_t> unique_reactant_indices;
std::vector<int> reactant_powers;
double symmetry_factor;
std::vector<size_t> affected_species_indices;
std::vector<int> stoichiometric_coefficients;
size_t reaction_index{};
reaction::ReactionType reaction_type{};
uint64_t reaction_hash{};
std::vector<size_t> unique_reactant_indices{};
std::vector<int> reactant_powers{};
double symmetry_factor{};
std::vector<size_t> affected_species_indices{};
std::vector<int> stoichiometric_coefficients{};
// Reverse cacheing
std::vector<size_t> unique_product_indices; ///< Unique product indices for reverse reactions.
std::vector<int> product_powers; ///< Powers of each unique product in the reverse reaction.
double reverse_symmetry_factor; ///< Symmetry factor for reverse reactions.
std::vector<size_t> unique_product_indices{}; ///< Unique product indices for reverse reactions.
std::vector<int> product_powers{}; ///< Powers of each unique product in the reverse reaction.
double reverse_symmetry_factor{}; ///< Symmetry factor for reverse reactions.
};
struct constants {
@@ -955,7 +955,7 @@ namespace gridfire {
[[nodiscard]] StepDerivatives<double> calculateAllDerivativesUsingPrecomputation(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const std::vector<double> &bare_rates,
const std::vector<double> &bare_reverse_rates,
double T9,

4
src/include/gridfire/engine/procedures/construction.h
View File

@@ -3,7 +3,7 @@
#include "gridfire/reaction/reaction.h"
#include "gridfire/engine/types/building.h"
#include "fourdst/composition/composition.h"
#include "fourdst/composition/composition_abstract.h"
#include <variant>
@@ -113,7 +113,7 @@ namespace gridfire {
* @throws std::logic_error If the resolved network depth is zero (no reactions can be collected).
*/
reaction::ReactionSet build_nuclear_network(
const fourdst::composition::Composition &composition,
const fourdst::composition::CompositionAbstract &composition,
const rates::weak::WeakRateInterpolator &weakInterpolator,
BuildDepthType maxLayers = NetworkBuildDepth::Full,
NetworkConstructionFlags ReactionTypes = NetworkConstructionFlags::DEFAULT

2
src/include/gridfire/engine/procedures/priming.h
View File

@@ -4,7 +4,7 @@
#include "gridfire/engine/engine_graph.h"
#include "gridfire/network.h"
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/atomic/atomicSpecies.h"
namespace gridfire {

4
src/include/gridfire/engine/types/reporting.h
View File

@@ -2,14 +2,12 @@
#include <map>
#include <string>
#include <ranges>
// Required for PrimingReport fields and streaming
#include <vector>
#include <utility>
#include <ostream>
#include <sstream>
#include "fourdst/composition/composition.h"
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/atomic/atomicSpecies.h"
namespace gridfire {

20
src/include/gridfire/engine/views/engine_adaptive.h
View File

@@ -5,7 +5,7 @@
#include "gridfire/screening/screening_types.h"
#include "gridfire/network.h"
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/atomic/atomicSpecies.h"
#include "fourdst/config/config.h"
#include "fourdst/logging/logging.h"
@@ -102,7 +102,7 @@ namespace gridfire {
* @see AdaptiveEngineView::update()
*/
[[nodiscard]] std::expected<StepDerivatives<double>, expectations::StaleEngineError> calculateRHSAndEnergy(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -116,7 +116,7 @@ namespace gridfire {
* @return A struct containing the derivatives of the energy generation rate with respect to temperature and density.
*/
[[nodiscard]] EnergyDerivatives calculateEpsDerivatives(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -135,20 +135,20 @@ namespace gridfire {
* @see AdaptiveEngineView::update()
*/
void generateJacobianMatrix(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
void generateJacobianMatrix(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho,
const std::vector<fourdst::atomic::Species> &activeSpecies
) const override;
void generateJacobianMatrix(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho,
const SparsityPattern &sparsityPattern
@@ -220,7 +220,7 @@ namespace gridfire {
*/
[[nodiscard]] double calculateMolarReactionFlow(
const reaction::Reaction &reaction,
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -248,13 +248,13 @@ namespace gridfire {
* @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called).
*/
[[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesTimescales(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
[[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesDestructionTimescales(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -303,7 +303,7 @@ namespace gridfire {
[[nodiscard]] PrimingReport primeEngine(const NetIn &netIn) override;
fourdst::composition::Composition collectComposition(fourdst::composition::Composition &comp) const override;
fourdst::composition::Composition collectComposition(fourdst::composition::CompositionAbstract &comp) const override;
private:
using Config = fourdst::config::Config;
using LogManager = fourdst::logging::LogManager;

18
src/include/gridfire/engine/views/engine_defined.h
View File

@@ -43,13 +43,13 @@ namespace gridfire{
* @throws std::runtime_error If the view is stale (i.e., `update()` has not been called after `setNetworkFile()`).
*/
[[nodiscard]] std::expected<StepDerivatives<double>, expectations::StaleEngineError> calculateRHSAndEnergy(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
[[nodiscard]] EnergyDerivatives calculateEpsDerivatives(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -64,7 +64,7 @@ namespace gridfire{
* @throws std::runtime_error If the view is stale.
*/
void generateJacobianMatrix(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -80,7 +80,7 @@ namespace gridfire{
* @throws std::runtime_error If the view is stale.
*/
void generateJacobianMatrix(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho,
const std::vector<fourdst::atomic::Species> &activeSpecies
@@ -97,7 +97,7 @@ namespace gridfire{
* @throws std::runtime_error If the view is stale.
*/
void generateJacobianMatrix(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho,
const SparsityPattern &sparsityPattern
@@ -150,7 +150,7 @@ namespace gridfire{
*/
[[nodiscard]] double calculateMolarReactionFlow(
const reaction::Reaction& reaction,
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -175,13 +175,13 @@ namespace gridfire{
* @throws std::runtime_error If the view is stale.
*/
[[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesTimescales(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
[[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesDestructionTimescales(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -221,7 +221,7 @@ namespace gridfire{
[[nodiscard]] PrimingReport primeEngine(const NetIn &netIn) override;
fourdst::composition::Composition collectComposition(fourdst::composition::Composition &comp) const override;
fourdst::composition::Composition collectComposition(fourdst::composition::CompositionAbstract &comp) const override;
protected:
bool m_isStale = true;
GraphEngine& m_baseEngine;

18
src/include/gridfire/engine/views/engine_multiscale.h
View File

@@ -231,13 +231,13 @@ namespace gridfire {
* (T9, rho, Y_full). This indicates `update()` was not called recently enough.
*/
[[nodiscard]] std::expected<StepDerivatives<double>, expectations::StaleEngineError> calculateRHSAndEnergy(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
[[nodiscard]] EnergyDerivatives calculateEpsDerivatives(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -266,7 +266,7 @@ namespace gridfire {
* without a valid partition.
*/
void generateJacobianMatrix(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -295,7 +295,7 @@ namespace gridfire {
* @throws exceptions::StaleEngineError If the QSE cache misses.
*/
void generateJacobianMatrix(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho,
const std::vector<fourdst::atomic::Species> &activeSpecies
@@ -323,7 +323,7 @@ namespace gridfire {
* @throws exceptions::StaleEngineError If the QSE cache misses.
*/
void generateJacobianMatrix(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho,
const SparsityPattern &sparsityPattern
@@ -407,7 +407,7 @@ namespace gridfire {
*/
[[nodiscard]] double calculateMolarReactionFlow(
const reaction::Reaction &reaction,
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -461,7 +461,7 @@ namespace gridfire {
* @throws StaleEngineError If the QSE cache misses.
*/
[[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesTimescales(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -487,7 +487,7 @@ namespace gridfire {
* @throws StaleEngineError If the QSE cache misses.
*/
[[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesDestructionTimescales(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -792,7 +792,7 @@ namespace gridfire {
* @return New composition which is comp + any edits from lower levels + the equilibrium abundances of all algebraic species.
* @throws BadCollectionError: if there is a species in the algebraic species set which does not show up in the reported composition from the base engine.:w
*/
fourdst::composition::Composition collectComposition(fourdst::composition::Composition &comp) const override;
fourdst::composition::Composition collectComposition(fourdst::composition::CompositionAbstract &comp) const override;
private:

2
src/include/gridfire/engine/views/engine_priming.h
View File

@@ -3,7 +3,7 @@
#include "gridfire/engine/views/engine_defined.h"
#include "fourdst/logging/logging.h"
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/atomic/atomicSpecies.h"
#include "quill/Logger.h"

3
src/include/gridfire/exceptions/error_policy.h
View File

@@ -13,6 +13,7 @@
#include <exception>
#include <string>
#include <utility>
namespace gridfire::exceptions {
/**
@@ -28,7 +29,7 @@ namespace gridfire::exceptions {
* @brief Constructs a PolicyError with a descriptive message.
* @param msg The error message.
*/
explicit PolicyError(const std::string& msg) : m_message(msg) {};
explicit PolicyError(std::string msg) : m_message(std::move(msg)) {};
/**
* @brief Returns the explanatory string.

16
src/include/gridfire/exceptions/error_utils.h
View File

@@ -5,14 +5,9 @@
#include <utility>
namespace gridfire::exceptions {
class UtilityError : public std::exception {};
class HashingError final : public UtilityError {
class UtilityError : public std::exception {
public:
explicit HashingError() = default;
explicit HashingError(std::string message)
: m_message(std::move(message)) {}
explicit UtilityError(std::string message) : m_message(std::move(message)) {}
[[nodiscard]] const char* what() const noexcept override {
return m_message.c_str();
@@ -20,4 +15,11 @@ namespace gridfire::exceptions {
private:
std::string m_message;
};
class HashingError final : public UtilityError {
public:
explicit HashingError(const std::string &message) : UtilityError(message) {}
};
}

40
src/include/gridfire/network.h
View File

@@ -22,13 +22,7 @@
#include <vector>
#include "fourdst/logging/logging.h"
#include "fourdst/config/config.h"
#include "fourdst/composition/species.h"
#include "fourdst/composition/composition.h"
#include "fourdst/constants/const.h"
#include "quill/Logger.h"
#include <unordered_map>
@@ -56,8 +50,6 @@ namespace gridfire {
double density; ///< Density in g/cm^3
double energy; ///< Energy in ergs
double culling = 0.0; ///< Culling threshold for reactions (default is 0.0, meaning no culling)
[[nodiscard]] std::vector<double> MolarAbundance() const;
};
struct NetOut {
@@ -73,36 +65,4 @@ namespace gridfire {
}
};
class Network {
public:
explicit Network(const NetworkFormat format = NetworkFormat::APPROX8);
virtual ~Network() = default;
[[nodiscard]] NetworkFormat getFormat() const;
NetworkFormat setFormat(const NetworkFormat format);
/**
* @brief Evaluate the network based on the input parameters.
*
* @param netIn Input parameters for the network evaluation.
* @return NetOut Output results from the network evaluation.
*/
virtual NetOut evaluate(const NetIn &netIn) = 0;
[[nodiscard]] virtual bool isStiff() const { return m_stiff; }
virtual void setStiff(const bool stiff) { m_stiff = stiff; }
protected:
fourdst::config::Config& m_config; ///< Configuration instance
fourdst::logging::LogManager& m_logManager; ///< Log manager instance
quill::Logger* m_logger; ///< Logger instance
NetworkFormat m_format; ///< Format of the network
fourdst::constant::Constants& m_constants;
bool m_stiff = false; ///< Flag indicating if the network is stiff
};
} // namespace nuclearNetwork

2
src/include/gridfire/policy/policy_abstract.h
View File

@@ -21,7 +21,7 @@
#pragma once
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/atomic/atomicSpecies.h"
#include "gridfire/reaction/reaction.h"
#include "gridfire/engine/engine_abstract.h"

2
src/include/gridfire/policy/policy_logical.h
View File

@@ -45,7 +45,7 @@ namespace gridfire::policy {
auto end() { return m_chain_policies.end(); }
[[nodiscard]] auto end() const { return m_chain_policies.cend(); }
protected:
std::vector<std::unique_ptr<ReactionChainPolicy>> m_chain_policies;
std::vector<std::unique_ptr<ReactionChainPolicy>> m_chain_policies{};
reaction::ReactionSet m_reactions;
};
}

17
src/include/gridfire/policy/stellar_policy.h
View File

@@ -25,7 +25,7 @@
#include "fourdst/composition/composition.h"
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/atomic/atomicSpecies.h"
#include "gridfire/partition/composite/partition_composite.h"
#include "gridfire/policy/chains.h"
@@ -91,7 +91,7 @@ namespace gridfire::policy {
* LowMassMainSequencePolicy policy(species, mass_fractions);
* @endcode
*/
explicit MainSequencePolicy(std::vector<fourdst::atomic::Species> seed_species, std::vector<double> mass_fractions);
explicit MainSequencePolicy(std::vector<fourdst::atomic::Species> seed_species, const std::vector<double> &mass_fractions);
/**
* @brief Returns the name of the policy.
@@ -143,16 +143,7 @@ namespace gridfire::policy {
*/
[[nodiscard]] NetworkPolicyStatus getStatus() const override;
private:
std::set<fourdst::atomic::Species> m_seed_species = {
fourdst::atomic::H_1,
fourdst::atomic::He_3,
fourdst::atomic::He_4,
fourdst::atomic::C_12,
fourdst::atomic::N_14,
fourdst::atomic::O_16,
fourdst::atomic::Ne_20,
fourdst::atomic::Mg_24
};
std::set<fourdst::atomic::Species> m_seed_species;
std::unique_ptr<ReactionChainPolicy> m_reaction_policy = std::make_unique<MainSequenceReactionChainPolicy>();
fourdst::composition::Composition m_initializing_composition;
@@ -162,7 +153,7 @@ namespace gridfire::policy {
NetworkPolicyStatus m_status = NetworkPolicyStatus::UNINITIALIZED;
private:
static std::unique_ptr<partition::PartitionFunction> build_partition_function();
NetworkPolicyStatus check_status() const;
[[nodiscard]] NetworkPolicyStatus check_status() const;
};

2
src/include/gridfire/reaction/reaction.h
View File

@@ -3,7 +3,7 @@
#include <ranges>
#include <string_view>
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/atomic/atomicSpecies.h"
#include "fourdst/logging/logging.h"
#include "quill/Logger.h"
#include <unordered_map>

2198
src/include/gridfire/reaction/reactions_data.h
View File

File diff suppressed because it is too large Load Diff

2
src/include/gridfire/reaction/weak/weak.h
View File

@@ -8,7 +8,7 @@
#include "gridfire/engine/engine_abstract.h"
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/atomic/atomicSpecies.h"
#include "fourdst/constants/const.h"
#include "cppad/cppad.hpp"

2
src/include/gridfire/reaction/weak/weak_interpolator.h
View File

@@ -1,7 +1,7 @@
#pragma once
#include "gridfire/reaction/weak/weak_types.h"
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/atomic/atomicSpecies.h"
#include "fourdst/logging/logging.h"
#include <unordered_map>

2
src/include/gridfire/screening/screening_abstract.h
View File

@@ -2,7 +2,7 @@
#include "gridfire/reaction/reaction.h"
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/atomic/atomicSpecies.h"
#include "cppad/cppad.hpp"

2
src/include/gridfire/solver/solver.h
View File

@@ -14,7 +14,7 @@ namespace gridfire::solver {
* @struct SolverContextBase
* @brief Base class for solver callback contexts.
*
* This struct serves as a base class for contexts that can be papubl;ssed to solver callbacks, it enforces
* This struct serves as a base class for contexts that can be passed to solver callbacks, it enforces
* that derived classes implement a `describe` method that returns a vector of tuples describing
* the context that a callback will receive when called.
*/

16
src/include/gridfire/solver/strategies/CVODE_solver_strategy.h
View File

@@ -5,7 +5,7 @@
#include "gridfire/network.h"
#include "gridfire/exceptions/exceptions.h"
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/atomic/atomicSpecies.h"
#include "fourdst/config/config.h"
@@ -176,7 +176,7 @@ namespace gridfire::solver {
const DynamicEngine& engine; ///< Reference to the engine.
const std::vector<fourdst::atomic::Species>& networkSpecies; ///< Species layout.
const size_t currentConvergenceFailures; ///< Total number of convergence failures
const size_t currentNonlinearIterations; ///< Total number of non linear iterations
const size_t currentNonlinearIterations; ///< Total number of non-linear iterations
/**
* @brief Construct a context snapshot.
@@ -216,12 +216,12 @@ namespace gridfire::solver {
* to CVODE, then the driver loop inspects and rethrows.
*/
struct CVODEUserData {
CVODESolverStrategy* solver_instance; // Pointer back to the class instance
DynamicEngine* engine;
double T9;
double rho;
double energy;
const std::vector<fourdst::atomic::Species>* networkSpecies;
CVODESolverStrategy* solver_instance{}; // Pointer back to the class instance
DynamicEngine* engine{};
double T9{};
double rho{};
double energy{};
const std::vector<fourdst::atomic::Species>* networkSpecies{};
std::unique_ptr<exceptions::StaleEngineTrigger> captured_exception = nullptr;
};

81
src/include/gridfire/utils/general_composition.h
View File

@@ -1,81 +0,0 @@
#pragma once
#include "fourdst/composition/composition.h"
#include "fourdst/composition/atomicSpecies.h"
#include <ranges>
namespace gridfire::utils {
inline double massFractionFromMolarAbundanceAndComposition (
const fourdst::composition::Composition& composition,
const fourdst::atomic::Species& species,
const double Yi
) {
double sum = 0;
for (const auto& [symbol, entry] : composition) {
if (entry.isotope() == species) {
sum += species.mass() * Yi;
} else {
sum += entry.isotope().mass() * composition.getMolarAbundance(symbol);
}
}
return (species.mass() * Yi) / sum;
};
/**
* @brief Convert a vector of molar abundances into a vector of mass fractions
* @param molarAbundances Vector of molar abundances
* @param molarMasses Vector of molar masses
*
* @note The vectors molarAbundances and molarMasses must be parallel. This function does not provide any checks
* to ensure that the correct molar mass is being used with the correct molar abundance.
* @return A vector of molar masses such that each molar mass < 1 and the sum of all is = 1
*/
inline std::vector<double> massFractionFromMolarAbundanceAndMolarMass (
const std::vector<double>& molarAbundances,
const std::vector<double>& molarMasses
) noexcept {
assert(molarMasses.size() == molarAbundances.size());
assert(!molarMasses.empty());
double totalMass = 0;
std::vector<double> masses;
masses.reserve(molarMasses.size());
for (const auto [m, Y] : std::views::zip(molarMasses, molarAbundances)) {
const double mass = m * Y;
totalMass += mass;
masses.push_back(mass);
}
assert(totalMass > 0);
std::vector<double> massFractions;
massFractions.reserve(masses.size());
std::ranges::transform(
masses,
std::back_inserter(massFractions),
[&totalMass](const double speciesMass) {
const double Xi = speciesMass / totalMass;
if (std::abs(Xi) < 1e-16 && Xi < 0) {
return 0.0;
}
return Xi;
});
return massFractions;
}
inline std::vector<double> molarMassVectorFromComposition(
const fourdst::composition::Composition& composition
) {
std::map<fourdst::atomic::Species, double> molarMassMap;
for (const auto &entry: composition | std::views::values) {
molarMassMap.emplace(entry.isotope(), entry.isotope().mass());
}
std::vector<double> molarMassVector;
molarMassVector.reserve(molarMassMap.size());
for (const auto molarMass : molarMassMap | std::views::values) {
molarMassVector.push_back(molarMass);
}
return molarMassVector;
}
}

122
src/include/gridfire/utils/table_format.h
View File

@@ -9,9 +9,50 @@
#include <algorithm>
#include <numeric>
#include <memory>
#include <format>
#include <print>
#include <cstdlib>
#include <cwchar>
#include <clocale>
namespace gridfire::utils {
inline size_t visual_width(const std::string& s) {
// IMPORTANT: std::setlocale(LC_ALL, "") must be called once in main()
// for mbtowc and wcwidth to function correctly with the system's locale.
size_t width = 0;
std::mbtowc(nullptr, 0, 0); // Reset multi-byte state
const char* p = s.c_str();
const char* end = s.c_str() + s.length();
while (p < end) {
wchar_t wc;
// Convert the next multi-byte char to a wide char
int byte_len = std::mbtowc(&wc, p, end - p);
if (byte_len <= 0) {
// Invalid byte sequence or null char.
// Treat as a 1-width '?' and advance by 1 byte to avoid infinite loop.
width++;
p++;
continue;
}
// Get the visual width of the wide char
int char_width = wcwidth(wc);
if (char_width != -1) {
width += char_width;
}
// else: char_width == -1 means non-printable/control char; treat as 0 width
p += byte_len; // Advance by the number of bytes consumed
}
return width;
}
class ColumnBase {
public:
virtual ~ColumnBase() = default;
@@ -117,6 +158,87 @@ namespace gridfire::utils {
return table_ss.str();
}
inline void print_table(const std::string& tableName, const std::vector<std::unique_ptr<ColumnBase>>& columns) {
// --- 1. Handle Empty Table ---
if (columns.empty()) {
std::println("{} \n(Table has no columns)\n", tableName);
return;
}
// --- 2. Determine dimensions and calculate column widths (using visual_width) ---
size_t num_cols = columns.size();
size_t num_rows = 0;
for (const auto& col : columns) {
num_rows = std::max(num_rows, col->getRowCount());
}
std::vector<size_t> col_widths(num_cols);
for (size_t j = 0; j < num_cols; ++j) {
// Start with header width
col_widths[j] = visual_width(columns[j]->getHeader());
// Find max width in all data cells
for (size_t i = 0; i < num_rows; ++i) {
col_widths[j] = std::max(col_widths[j], visual_width(columns[j]->getCellData(i)));
}
}
// --- 3. Print the table using std::print / std::println ---
// --- Table Title ---
// NOLINTNEXTLINE(*-fold-init-type)
const size_t total_width = std::accumulate(col_widths.begin(), col_widths.end(), 0UL) + (num_cols * 3) + 1;
const size_t title_padding_len = (total_width > visual_width(tableName)) ? (total_width - visual_width(tableName)) / 2 : 0;
// Print padding, then title
std::print("{: <{}}", "", title_padding_len); // Left-aligned empty string "" of width title_padding_len
std::println("{}", tableName);
// --- Helper to draw horizontal border ---
auto draw_border = [&]() {
std::print("+");
for (const size_t width : col_widths) {
// std::string(width + 2, '-') is still the easiest way to repeat a char
std::print("{:-<{}}+", "", width + 2); // Prints '-' repeated (width + 2) times
}
std::println("");
};
// --- Draw Top Border ---
draw_border();
// --- Helper to print a cell with correct padding ---
auto print_cell = [&](const std::string& text, size_t width) {
size_t text_width = visual_width(text);
size_t padding = (width >= text_width) ? (width - text_width) : 0;
// Print text and then the manual padding
std::print(" {}{: <{}} |", text, "", padding);
};
// --- Draw Header Row ---
std::print("|");
for (size_t j = 0; j < num_cols; ++j) {
print_cell(columns[j]->getHeader(), col_widths[j]);
}
std::println("");
// --- Draw Separator ---
draw_border();
// --- Draw Data Rows ---
for (size_t i = 0; i < num_rows; ++i) {
std::print("|");
for (size_t j = 0; j < num_cols; ++j) {
print_cell(columns[j]->getCellData(i), col_widths[j]);
}
std::println("");
}
// --- Draw Bottom Border ---
draw_border();
}
}

1
src/lib/engine/diagnostics/dynamic_engine_diagnostics.cpp
View File

@@ -1,6 +1,7 @@
#include "gridfire/engine/diagnostics/dynamic_engine_diagnostics.h"
#include "gridfire/engine/engine_abstract.h"
#include "gridfire/utils/table_format.h"
#include "fourdst/atomic/species.h"
#include <vector>
#include <string>

529
src/lib/engine/engine_approx8.cpp
View File

@@ -1,529 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include <cmath>
#include <stdexcept>
#include <array>
#include <boost/numeric/odeint.hpp>
#include "fourdst/constants/const.h"
#include "fourdst/config/config.h"
#include "quill/LogMacros.h"
#include "gridfire/engine/engine_approx8.h"
#include "gridfire/network.h"
/* Nuclear reaction network in cgs units based on Frank Timmes' "approx8".
At this time it does neither screening nor neutrino losses. It includes
the following 8 isotopes:
h1
he3
he4
c12
n14
o16
ne20
mg24
and the following nuclear reactions:
---pp chain---
p(p,e+)d
d(p,g)he3
he3(he3,2p)he4
---CNO cycle---
c12(p,g)n13 - n13 -> c13 + p -> n14
n14(p,g)o15 - o15 + p -> c12 + he4
n14(a,g)f18 - proceeds to ne20
n15(p,a)c12 - / these two n15 reactions are \ CNO I
n15(p,g)o16 - \ used to calculate a fraction / CNO II
o16(p,g)f17 - f17 + e -> o17 and then o17 + p -> n14 + he4
---alpha captures---
c12(a,g)o16
triple alpha
o16(a,g)ne20
ne20(a,g)mg24
c12(c12,a)ne20
c12(o16,a)mg24
At present the rates are all evaluated using a fitting function.
The coefficients to the fit are from reaclib.jinaweb.org .
*/
namespace gridfire::approx8{
// using namespace std;
using namespace boost::numeric::odeint;
//helper functions
// a function to multiply two arrays and then sum the resulting elements: sum(a*b)
double sum_product( const vec7 &a, const vec7 &b){
double sum=0;
for (size_t i=0; i < a.size(); i++) {
sum += a[i] * b[i];
}
return sum;
}
// the fit to the nuclear reaction rates is of the form:
// exp( a0 + a1/T9 + a2/T9^(1/3) + a3*T9^(1/3) + a4*T9 + a5*T9^(5/3) + log(T9) )
// this function returns an array of the T9 terms in that order, where T9 is the temperatures in GigaKelvin
vec7 get_T9_array(const double &T) {
vec7 arr;
const double T9=1e-9*T;
const double T913=pow(T9,1./3.);
arr[0]=1;
arr[1]=1/T9;
arr[2]=1/T913;
arr[3]=T913;
arr[4]=T9;
arr[5]=pow(T9,5./3.);
arr[6]=log(T9);
return arr;
}
// this function uses the two preceding functions to evaluate the rate given the T9 array and the coefficients
double rate_fit(const vec7 &T9, const vec7 &coef){
return exp(sum_product(T9,coef));
}
// p + p -> d; this, like some of the other rates, this is a composite of multiple fits
double pp_rate(const vec7 &T9) {
constexpr vec7 a1 = {-34.78630, 0,-3.511930, 3.100860, -0.1983140, 1.262510e-2, -1.025170};
constexpr vec7 a2 = { -4.364990e+1,-2.460640e-3,-2.750700,-4.248770e-1,1.598700e-2,-6.908750e-4,-2.076250e-1};
return rate_fit(T9,a1) + rate_fit(T9,a2);
}
// p + d -> he3
double dp_rate(const vec7 &T9) {
constexpr vec7 a1 = {7.528980, 0, -3.720800, 0.8717820, 0, 0,-0.6666670};
constexpr vec7 a2 = {8.935250, 0, -3.720800, 0.1986540, 0, 0, 0.3333330};
return rate_fit(T9,a1) + rate_fit(T9,a2);
}
// he3 + he3 -> he4 + 2p
double he3he3_rate(const vec7 &T9){
constexpr vec7 a = {2.477880e+01,0,-12.27700,-0.1036990,-6.499670e-02,1.681910e-02,-6.666670e-01};
return rate_fit(T9,a);
}
// he3(he3,2p)he4 ** (missing both coefficients but have a reaction)
double he3he4_rate(const vec7 &T9){
constexpr vec7 a1 = {1.560990e+01,0.000000e+00,-1.282710e+01,-3.082250e-02,-6.546850e-01,8.963310e-02,-6.666670e-01};
constexpr vec7 a2 = {1.770750e+01,0.000000e+00,-1.282710e+01,-3.812600e+00,9.422850e-02,-3.010180e-03,1.333330e+00};
return rate_fit(T9,a1) + rate_fit(T9,a2);
}
// he4 + he4 + he4 -> c12 ** (missing middle coefficient but have other two)
double triple_alpha_rate(const vec7 &T9){
constexpr vec7 a1 = {-9.710520e-01,0.000000e+00,-3.706000e+01,2.934930e+01,-1.155070e+02,-1.000000e+01,-1.333330e+00};
constexpr vec7 a2 = {-1.178840e+01,-1.024460e+00,-2.357000e+01,2.048860e+01,-1.298820e+01,-2.000000e+01,-2.166670e+00};
constexpr vec7 a3 = {-2.435050e+01,-4.126560e+00,-1.349000e+01,2.142590e+01,-1.347690e+00,8.798160e-02,-1.316530e+01};
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
}
// c12 + p -> n13
double c12p_rate(const vec7 &T9){
constexpr vec7 a1={1.714820e+01,0.000000e+00,-1.369200e+01,-2.308810e-01,4.443620e+00,-3.158980e+00,-6.666670e-01};
constexpr vec7 a2={1.754280e+01,-3.778490e+00,-5.107350e+00,-2.241110e+00,1.488830e-01,0.000000e+00,-1.500000e+00};
return rate_fit(T9,a1) + rate_fit(T9,a2);
}
// c12 + he4 -> o16 ** (missing first coefficient but have the second)
double c12a_rate(const vec7 &T9){
constexpr vec7 a1={6.965260e+01,-1.392540e+00,5.891280e+01,-1.482730e+02,9.083240e+00,-5.410410e-01,7.035540e+01};
constexpr vec7 a2={2.546340e+02,-1.840970e+00,1.034110e+02,-4.205670e+02,6.408740e+01,-1.246240e+01,1.373030e+02};
return rate_fit(T9,a1) + rate_fit(T9,a2);
}
// n14(p,g)o15 - o15 + p -> c12 + he4
double n14p_rate(const vec7 &T9){
constexpr vec7 a1={1.701000e+01,0.000000e+00,-1.519300e+01,-1.619540e-01,-7.521230e+00,-9.875650e-01,-6.666670e-01};
constexpr vec7 a2={2.011690e+01,0.000000e+00,-1.519300e+01,-4.639750e+00,9.734580e+00,-9.550510e+00,3.333330e-01};
constexpr vec7 a3={7.654440e+00,-2.998000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
constexpr vec7 a4={6.735780e+00,-4.891000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,6.820000e-02};
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3) + rate_fit(T9,a4);
}
// n14(a,g)f18 assumed to go on to ne20
double n14a_rate(const vec7 &T9){
constexpr vec7 a1={2.153390e+01,0.000000e+00,-3.625040e+01,0.000000e+00,0.000000e+00,-5.000000e+00,-6.666670e-01};
constexpr vec7 a2={1.968380e-01,-5.160340e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
constexpr vec7 a3={1.389950e+01,-1.096560e+01,-5.622700e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
}
// n15(p,a)c12 (CNO I)
double n15pa_rate(const vec7 &T9){
constexpr vec7 a1 = {2.747640e+01,0.000000e+00,-1.525300e+01,1.593180e+00,2.447900e+00,-2.197080e+00,-6.666670e-01};
constexpr vec7 a2 = {-4.873470e+00,-2.021170e+00,0.000000e+00,3.084970e+01,-8.504330e+00,-1.544260e+00,-1.500000e+00};
constexpr vec7 a3 = {2.089720e+01,-7.406000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
constexpr vec7 a4 = {-6.575220e+00,-1.163800e+00,0.000000e+00,2.271050e+01,-2.907070e+00,2.057540e-01,-1.500000e+00};
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3) + rate_fit(T9,a4);
}
// n15(p,g)o16 (CNO II)
double n15pg_rate(const vec7 &T9){
constexpr vec7 a1 = {2.001760e+01,0.000000e+00,-1.524000e+01,3.349260e-01,4.590880e+00,-4.784680e+00,-6.666670e-01};
constexpr vec7 a2 = {6.590560e+00,-2.923150e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
constexpr vec7 a3 = {1.454440e+01,-1.022950e+01,0.000000e+00,0.000000e+00,4.590370e-02,0.000000e+00,-1.500000e+00};
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
}
double n15pg_frac(const vec7 &T9){
const double f1=n15pg_rate(T9);
const double f2=n15pa_rate(T9);
return f1/(f1+f2);
}
// o16(p,g)f17 then f17 -> o17(p,a)n14
double o16p_rate(const vec7 &T9){
constexpr vec7 a={1.909040e+01,0.000000e+00,-1.669600e+01,-1.162520e+00,2.677030e-01,-3.384110e-02,-6.666670e-01};
return rate_fit(T9,a);
}
// o16(a,g)ne20
double o16a_rate(const vec7 &T9){
constexpr vec7 a1={2.390300e+01,0.000000e+00,-3.972620e+01,-2.107990e-01,4.428790e-01,-7.977530e-02,-6.666670e-01};
constexpr vec7 a2={3.885710e+00,-1.035850e+01,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
constexpr vec7 a3={9.508480e+00,-1.276430e+01,0.000000e+00,-3.659250e+00,7.142240e-01,-1.075080e-03,-1.500000e+00};
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
}
// ne20(a,g)mg24
double ne20a_rate(const vec7 &T9){
constexpr vec7 a1={2.450580e+01,0.000000e+00,-4.625250e+01,5.589010e+00,7.618430e+00,-3.683000e+00,-6.666670e-01};
constexpr vec7 a2={-3.870550e+01,-2.506050e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
constexpr vec7 a3={1.983070e+00,-9.220260e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
constexpr vec7 a4={-8.798270e+00,-1.278090e+01,0.000000e+00,1.692290e+01,-2.573250e+00,2.089970e-01,-1.500000e+00};
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3) + rate_fit(T9,a4);
}
// c12(c12,a)ne20
double c12c12_rate(const vec7 &T9){
constexpr vec7 a={6.128630e+01,0.000000e+00,-8.416500e+01,-1.566270e+00,-7.360840e-02,-7.279700e-02,-6.666670e-01};
return rate_fit(T9,a);
}
// c12(o16,a)mg24
double c12o16_rate(const vec7 &T9){
constexpr vec7 a={4.853410e+01,3.720400e-01,-1.334130e+02,5.015720e+01,-3.159870e+00,1.782510e-02,-2.370270e+01};
return rate_fit(T9,a);
}
// for Boost ODE solvers either a struct or a class is required
// a Jacobian matrix for implicit solvers
void Jacobian::operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) const {
fourdst::constant::Constants& constants = fourdst::constant::Constants::getInstance();
const double avo = constants.get("N_a").value;
const double clight = constants.get("c").value;
// EOS
const vec7 T9=get_T9_array(y[Approx8Net::iTemp]);
// evaluate rates once per call
const double rpp=pp_rate(T9);
const double r33=he3he3_rate(T9);
const double r34=he3he4_rate(T9);
const double r3a=triple_alpha_rate(T9);
const double rc12p=c12p_rate(T9);
const double rc12a=c12a_rate(T9);
const double rn14p=n14p_rate(T9);
const double rn14a=n14a_rate(T9);
const double ro16p=o16p_rate(T9);
const double ro16a=o16a_rate(T9);
const double rne20a=ne20a_rate(T9);
const double r1212=c12c12_rate(T9);
const double r1216=c12o16_rate(T9);
const double pFrac=n15pg_frac(T9);
const double aFrac=1-pFrac;
const double yh1 = y[Approx8Net::ih1];
const double yhe3 = y[Approx8Net::ihe3];
const double yhe4 = y[Approx8Net::ihe4];
const double yc12 = y[Approx8Net::ic12];
const double yn14 = y[Approx8Net::in14];
const double yo16 = y[Approx8Net::io16];
const double yne20 = y[Approx8Net::ine20];
// zero all elements to begin
for (int i=0; i < Approx8Net::nVar; i++) {
for (int j=0; j < Approx8Net::nVar; j++) {
J(i,j)=0.0;
}
}
// h1 jacobian elements
J(Approx8Net::ih1,Approx8Net::ih1) = -3*yh1*rpp - 2*yc12*rc12p -2*yn14*rn14p -2*yo16*ro16p;
J(Approx8Net::ih1,Approx8Net::ihe3) = 2*yhe3*r33 - yhe4*r34;
J(Approx8Net::ih1,Approx8Net::ihe4) = -yhe3*r34;
J(Approx8Net::ih1,Approx8Net::ic12) = -2*yh1*rc12p;
J(Approx8Net::ih1,Approx8Net::in14) = -2*yh1*rn14p;
J(Approx8Net::ih1,Approx8Net::io16) = -2*yh1*ro16p;
// he3 jacobian elements
J(Approx8Net::ihe3,Approx8Net::ih1) = yh1*rpp;
J(Approx8Net::ihe3,Approx8Net::ihe3) = -2*yhe3*r33 - yhe4*r34;
J(Approx8Net::ihe3,Approx8Net::ihe4) = -yhe3*r34;
// he4 jacobian elements
J(Approx8Net::ihe4,Approx8Net::ih1) = yn14*aFrac*rn14p + yo16*ro16p;
J(Approx8Net::ihe4,Approx8Net::ihe3) = yhe3*r33 - yhe4*r34;
J(Approx8Net::ihe4,Approx8Net::ihe4) = yhe3*r34 - 1.5*yhe4*yhe4*r3a - yc12*rc12a - 1.5*yn14*rn14a - yo16*ro16a - yne20*rne20a;
J(Approx8Net::ihe4,Approx8Net::ic12) = -yhe4*rc12a + yc12*r1212 + yo16*r1216;
J(Approx8Net::ihe4,Approx8Net::in14) = yh1*aFrac*rn14p - 1.5*yhe4*rn14a;
J(Approx8Net::ihe4,Approx8Net::io16) = yh1*ro16p - yhe4*ro16a + yc12*r1216;
J(Approx8Net::ihe4,Approx8Net::ine20) = -yhe4*rne20a;
// c12 jacobian elements
J(Approx8Net::ic12,Approx8Net::ih1) = -yc12*rc12p + yn14*aFrac*rn14p;
J(Approx8Net::ic12,Approx8Net::ihe4) = 0.5*yhe4*yhe4*r3a - yhe4*rc12a;
J(Approx8Net::ic12,Approx8Net::ic12) = -yh1*rc12p - yhe4*rc12a - yo16*r1216 - 2*yc12*r1212;
J(Approx8Net::ic12,Approx8Net::in14) = yh1*yn14*aFrac*rn14p;
J(Approx8Net::ic12,Approx8Net::io16) = -yc12*r1216;
// n14 jacobian elements
J(Approx8Net::in14,Approx8Net::ih1) = yc12*rc12p - yn14*rn14p + yo16*ro16p;
J(Approx8Net::in14,Approx8Net::ihe4) = -yn14*rn14a;
J(Approx8Net::in14,Approx8Net::ic12) = yh1*rc12p;
J(Approx8Net::in14,Approx8Net::in14) = -yh1*rn14p - yhe4*rn14a;
J(Approx8Net::in14,Approx8Net::io16) = yo16*ro16p;
// o16 jacobian elements
J(Approx8Net::io16,Approx8Net::ih1) = yn14*pFrac*rn14p - yo16*ro16p;
J(Approx8Net::io16,Approx8Net::ihe4) = yc12*rc12a - yo16*ro16a;
J(Approx8Net::io16,Approx8Net::ic12) = yhe4*rc12a - yo16*r1216;
J(Approx8Net::io16,Approx8Net::in14) = yh1*pFrac*rn14p;
J(Approx8Net::io16,Approx8Net::io16) = yh1*ro16p - yc12*r1216 -yhe4*ro16a;
// ne20 jacobian elements
J(Approx8Net::ine20,Approx8Net::ihe4) = yn14*rn14a + yo16*ro16a - yne20*rne20a;
J(Approx8Net::ine20,Approx8Net::ic12) = yc12*r1212;
J(Approx8Net::ine20,Approx8Net::in14) = yhe4*rn14a;
J(Approx8Net::ine20,Approx8Net::io16) = yo16*ro16a;
J(Approx8Net::ine20,Approx8Net::ine20) = -yhe4*rne20a;
// mg24 jacobian elements
J(Approx8Net::img24,Approx8Net::ihe4) = yne20*rne20a;
J(Approx8Net::img24,Approx8Net::ic12) = yo16*r1216;
J(Approx8Net::img24,Approx8Net::io16) = yc12*r1216;
J(Approx8Net::img24,Approx8Net::ine20) = yhe4*rne20a;
// energy accounting
for (int j=0; j<Approx8Net::nIso; j++) {
for (int i=0; i<Approx8Net::nIso; i++) {
J(Approx8Net::iEnergy,j) += J(i,j)*Approx8Net::mIon[i];
}
J(Approx8Net::iEnergy,j) *= -avo*clight*clight;
}
}
void ODE::operator() ( const vector_type &y, vector_type &dydt, double /* t */) const {
const fourdst::constant::Constants& constants = fourdst::constant::Constants::getInstance();
const double avo = constants.get("N_a").value;
const double clight = constants.get("c").value;
// EOS
const double T = y[Approx8Net::iTemp];
const double den = y[Approx8Net::iDensity];
const vec7 T9=get_T9_array(T);
// rates
const double rpp=den*pp_rate(T9);
const double r33=den*he3he3_rate(T9);
const double r34=den*he3he4_rate(T9);
const double r3a=den*den*triple_alpha_rate(T9);
const double rc12p=den*c12p_rate(T9);
const double rc12a=den*c12a_rate(T9);
const double rn14p=den*n14p_rate(T9);
const double rn14a=n14a_rate(T9);
const double ro16p=den*o16p_rate(T9);
const double ro16a=den*o16a_rate(T9);
const double rne20a=den*ne20a_rate(T9);
const double r1212=den*c12c12_rate(T9);
const double r1216=den*c12o16_rate(T9);
const double pFrac=n15pg_frac(T9);
const double aFrac=1-pFrac;
const double yh1 = y[Approx8Net:: ih1];
const double yhe3 = y[Approx8Net:: ihe3];
const double yhe4 = y[Approx8Net:: ihe4];
const double yc12 = y[Approx8Net:: ic12];
const double yn14 = y[Approx8Net:: in14];
const double yo16 = y[Approx8Net:: io16];
const double yne20 = y[Approx8Net::ine20];
dydt[Approx8Net::ih1] = -1.5*yh1*yh1*rpp;
dydt[Approx8Net::ih1] += yhe3*yhe3*r33;
dydt[Approx8Net::ih1] += -yhe3*yhe4*r34;
dydt[Approx8Net::ih1] += -2*yh1*yc12*rc12p;
dydt[Approx8Net::ih1] += -2*yh1*yn14*rn14p;
dydt[Approx8Net::ih1] += -2*yh1*yo16*ro16p;
dydt[Approx8Net::ihe3] = 0.5*yh1*yh1*rpp;
dydt[Approx8Net::ihe3] += -yhe3*yhe3*r33;
dydt[Approx8Net::ihe3] += -yhe3*yhe4*r34;
dydt[Approx8Net::ihe4] = 0.5*yhe3*yhe3*r33;
dydt[Approx8Net::ihe4] += yhe3*yhe4*r34;
dydt[Approx8Net::ihe4] += -yhe4*yc12*rc12a;
dydt[Approx8Net::ihe4] += yh1*yn14*aFrac*rn14p;
dydt[Approx8Net::ihe4] += yh1*yo16*ro16p;
dydt[Approx8Net::ihe4] += -0.5*yhe4*yhe4*yhe4*r3a;
dydt[Approx8Net::ihe4] += -yhe4*yo16*ro16a;
dydt[Approx8Net::ihe4] += 0.5*yc12*yc12*r1212;
dydt[Approx8Net::ihe4] += yc12*yo16*r1216;
dydt[Approx8Net::ihe4] += -yhe4*yne20*rne20a;
dydt[Approx8Net::ic12] = (1./6.)*yhe4*yhe4*yhe4*r3a;
dydt[Approx8Net::ic12] += -yhe4*yc12*rc12a;
dydt[Approx8Net::ic12] += -yh1*yc12*rc12p;
dydt[Approx8Net::ic12] += yh1*yn14*aFrac*rn14p;
dydt[Approx8Net::ic12] += -yc12*yc12*r1212;
dydt[Approx8Net::ic12] += -yc12*yo16*r1216;
dydt[Approx8Net::in14] = yh1*yc12*rc12p;
dydt[Approx8Net::in14] += -yh1*yn14*rn14p;
dydt[Approx8Net::in14] += yh1*yo16*ro16p;
dydt[Approx8Net::in14] += -yhe4*yn14*rn14a;
dydt[Approx8Net::io16] = yhe4*yc12*rc12a;
dydt[Approx8Net::io16] += yh1*yn14*pFrac*rn14p;
dydt[Approx8Net::io16] += -yh1*yo16*ro16p;
dydt[Approx8Net::io16] += -yc12*yo16*r1216;
dydt[Approx8Net::io16] += -yhe4*yo16*ro16a;
dydt[Approx8Net::ine20] = 0.5*yc12*yc12*r1212;
dydt[Approx8Net::ine20] += yhe4*yn14*rn14a;
dydt[Approx8Net::ine20] += yhe4*yo16*ro16a;
dydt[Approx8Net::ine20] += -yhe4*yne20*rne20a;
dydt[Approx8Net::img24] = yc12*yo16*r1216;
dydt[Approx8Net::img24] += yhe4*yne20*rne20a;
dydt[Approx8Net::iTemp] = 0.;
dydt[Approx8Net::iDensity] = 0.;
// energy accounting
double eNuc = 0.;
for (int i=0; i<Approx8Net::nIso; i++) {
eNuc += Approx8Net::mIon[i]*dydt[i];
}
dydt[Approx8Net::iEnergy] = -eNuc*avo*clight*clight;
}
Approx8Network::Approx8Network() : Network(APPROX8) {}
NetOut Approx8Network::evaluate(const NetIn &netIn) {
m_y = convert_netIn(netIn);
m_tMax = netIn.tMax;
m_dt0 = netIn.dt0;
const auto stiff_abs_tol = m_config.get<double>("Network:Approx8:Stiff:AbsTol", 1.0e-6);
const auto stiff_rel_tol = m_config.get<double>("Network:Approx8:Stiff:RelTol", 1.0e-6);
const auto nonstiff_abs_tol = m_config.get<double>("Network:Approx8:NonStiff:AbsTol", 1.0e-6);
const auto nonstiff_rel_tol = m_config.get<double>("Network:Approx8:NonStiff:RelTol", 1.0e-6);
int num_steps = -1;
if (m_stiff) {
LOG_DEBUG(m_logger, "Using stiff solver for Approx8Network");
num_steps = integrate_const(
make_dense_output<rosenbrock4<double>>(stiff_abs_tol, stiff_rel_tol),
std::make_pair(ODE(), Jacobian()),
m_y,
0.0,
m_tMax,
m_dt0
);
} else {
LOG_DEBUG(m_logger, "Using non stiff solver for Approx8Network");
num_steps = integrate_const (
make_dense_output<runge_kutta_dopri5<vector_type>>(nonstiff_abs_tol, nonstiff_rel_tol),
ODE(),
m_y,
0.0,
m_tMax,
m_dt0
);
}
double ySum = 0.0;
for (int i = 0; i < Approx8Net::nIso; i++) {
m_y[i] *= Approx8Net::aIon[i];
ySum += m_y[i];
}
for (int i = 0; i < Approx8Net::nIso; i++) {
m_y[i] /= ySum;
}
NetOut netOut;
std::vector<double> outComposition;
outComposition.reserve(Approx8Net::nVar);
for (int i = 0; i < Approx8Net::nIso; i++) {
outComposition.push_back(m_y[i]);
}
netOut.energy = m_y[Approx8Net::iEnergy];
netOut.num_steps = num_steps;
const std::vector<std::string> symbols = {"H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"};
netOut.composition = fourdst::composition::Composition(symbols, outComposition);
return netOut;
}
void Approx8Network::setStiff(bool stiff) {
m_stiff = stiff;
}
vector_type Approx8Network::convert_netIn(const NetIn &netIn) {
vector_type y(Approx8Net::nVar, 0.0);
y[Approx8Net::ih1] = netIn.composition.getNumberFraction("H-1");
y[Approx8Net::ihe3] = netIn.composition.getNumberFraction("He-3");
y[Approx8Net::ihe4] = netIn.composition.getNumberFraction("He-4");
y[Approx8Net::ic12] = netIn.composition.getNumberFraction("C-12");
y[Approx8Net::in14] = netIn.composition.getNumberFraction("N-14");
y[Approx8Net::io16] = netIn.composition.getNumberFraction("O-16");
y[Approx8Net::ine20] = netIn.composition.getNumberFraction("Ne-20");
y[Approx8Net::img24] = netIn.composition.getNumberFraction("Mg-24");
y[Approx8Net::iTemp] = netIn.temperature;
y[Approx8Net::iDensity] = netIn.density;
y[Approx8Net::iEnergy] = netIn.energy;
return y;
}
};
// main program

99
src/lib/engine/engine_graph.cpp
View File

@@ -8,8 +8,8 @@
#include "gridfire/utils/hashing.h"
#include "gridfire/utils/table_format.h"
#include "fourdst/composition/species.h"
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/atomic/species.h"
#include "fourdst/atomic/atomicSpecies.h"
#include "quill/LogMacros.h"
@@ -66,7 +66,7 @@ namespace gridfire {
}
std::expected<StepDerivatives<double>, expectations::StaleEngineError> GraphEngine::calculateRHSAndEnergy(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
@@ -74,7 +74,7 @@ namespace gridfire {
}
std::expected<StepDerivatives<double>, expectations::StaleEngineError> GraphEngine::calculateRHSAndEnergy(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho,
const reaction::ReactionSet &activeReactions
@@ -119,7 +119,7 @@ namespace gridfire {
}
EnergyDerivatives GraphEngine::calculateEpsDerivatives(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
@@ -127,7 +127,7 @@ namespace gridfire {
}
EnergyDerivatives GraphEngine::calculateEpsDerivatives(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho,
const reaction::ReactionSet &activeReactions
@@ -357,7 +357,7 @@ namespace gridfire {
const reaction::Reaction &reaction,
const double T9,
const double rho,
const fourdst::composition::Composition &comp
const fourdst::composition::CompositionAbstract &comp
) const {
if (!m_useReverseReactions) {
LOG_TRACE_L3_LIMIT_EVERY_N(std::numeric_limits<int>::max(), m_logger, "Reverse reactions are disabled. Returning 0.0 for reverse rate of reaction '{}'.", reaction.id());
@@ -534,8 +534,8 @@ namespace gridfire {
std::vector<double> GraphEngine::mapNetInToMolarAbundanceVector(const NetIn &netIn) const {
std::vector<double> Y(m_networkSpecies.size(), 0.0); // Initialize with zeros
for (const auto& [symbol, entry] : netIn.composition) {
Y[getSpeciesIndex(entry.isotope())] = netIn.composition.getMolarAbundance(symbol); // Map species to their molar abundance
for (const auto& [sp, y] : netIn.composition) {
Y[getSpeciesIndex(sp)] = y; // Map species to their molar abundance
}
return Y; // Return the vector of molar abundances
}
@@ -544,24 +544,13 @@ namespace gridfire {
NetIn fullNetIn;
fourdst::composition::Composition composition;
std::vector<std::string> symbols;
symbols.reserve(m_networkSpecies.size());
for (const auto &symbol: m_networkSpecies) {
symbols.emplace_back(symbol.name());
}
composition.registerSymbol(symbols);
for (const auto& [symbol, entry] : netIn.composition) {
if (m_networkSpeciesMap.contains(symbol)) {
composition.setMassFraction(symbol, entry.mass_fraction());
} else {
composition.setMassFraction(symbol, 0.0);
for (const auto& sp : m_networkSpecies) {
composition.registerSpecies(sp);
if (netIn.composition.contains(sp)) {
composition.setMolarAbundance(sp, netIn.composition.getMolarAbundance(sp));
}
}
const bool didFinalize = composition.finalize(true);
if (!didFinalize) {
LOG_ERROR(m_logger, "Failed to finalize composition during priming. Check input mass fractions for validity.");
throw std::runtime_error("Failed to finalize composition during network priming.");
}
fullNetIn.composition = composition;
fullNetIn.temperature = netIn.temperature;
fullNetIn.density = netIn.density;
@@ -580,7 +569,7 @@ namespace gridfire {
return m_depth;
}
void GraphEngine::rebuild(const fourdst::composition::Composition& comp, const BuildDepthType depth) {
void GraphEngine::rebuild(const fourdst::composition::CompositionAbstract &comp, const BuildDepthType depth) {
if (depth != m_depth) {
m_depth = depth;
m_reactions = build_nuclear_network(comp, m_weakRateInterpolator, m_depth);
@@ -592,27 +581,25 @@ namespace gridfire {
}
fourdst::composition::Composition GraphEngine::collectComposition(
fourdst::composition::Composition &comp
fourdst::composition::CompositionAbstract &comp
) const {
for (const auto &speciesName: comp | std::views::keys) {
if (!m_networkSpeciesMap.contains(speciesName)) {
throw exceptions::BadCollectionError("Cannot collect composition from GraphEngine as " + speciesName + " present in input composition does not exist in the network species map");
for (const auto &species: comp.getRegisteredSpecies()) {
if (!m_networkSpeciesMap.contains(species.name())) {
throw exceptions::BadCollectionError("Cannot collect composition from GraphEngine as " + std::string(species.name()) + " present in input composition does not exist in the network species map");
}
}
fourdst::composition::Composition result;
for (const auto& species : m_networkSpecies ) {
result.registerSpecies(species);
if (comp.hasSpecies(species)) {
result.setMassFraction(species, comp.getMassFraction(species));
} else {
result.setMassFraction(species, 0.0);
if (comp.contains(species)) {
result.setMolarAbundance(species, comp.getMolarAbundance(species));
}
}
return result;
}
StepDerivatives<double> GraphEngine::calculateAllDerivativesUsingPrecomputation(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
const std::vector<double> &bare_rates,
const std::vector<double> &bare_reverse_rates,
const double T9,
@@ -646,6 +633,10 @@ namespace gridfire {
const fourdst::atomic::Species& reactant = m_networkSpecies[reactantIndex];
const int power = precomputedReaction.reactant_powers[i];
if (!comp.contains(reactant)) {
forwardAbundanceProduct = 0.0;
break; // No need to continue if one of the reactants has zero abundance
}
forwardAbundanceProduct *= std::pow(comp.getMolarAbundance(reactant), power);
}
@@ -787,7 +778,7 @@ namespace gridfire {
double GraphEngine::calculateMolarReactionFlow(
const reaction::Reaction &reaction,
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
@@ -811,22 +802,17 @@ namespace gridfire {
}
void GraphEngine::generateJacobianMatrix(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
fourdst::composition::Composition mutableComp = comp;
fourdst::composition::Composition mutableComp;
for (const auto& species : m_networkSpecies) {
if (!comp.hasSpecies(species)) {
mutableComp.registerSpecies(species);
mutableComp.setMassFraction(species, 0.0);
mutableComp.registerSpecies(species);
if (comp.contains(species)) {
mutableComp.setMolarAbundance(species, comp.getMolarAbundance(species));
}
}
const bool didFinalize = mutableComp.finalize(false);
if (!didFinalize) {
LOG_CRITICAL(m_logger, "Could not finalize the composition used to generate the jacobian matrix!");
throw std::runtime_error("Could not finalize the composition used to generate the jacobian matrix");
}
LOG_TRACE_L1_LIMIT_EVERY_N(1000, m_logger, "Generating jacobian matrix for T9={}, rho={}..", T9, rho);
const size_t numSpecies = m_networkSpecies.size();
@@ -863,7 +849,7 @@ namespace gridfire {
}
void GraphEngine::generateJacobianMatrix(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho,
const std::vector<fourdst::atomic::Species> &activeSpecies
@@ -895,7 +881,7 @@ namespace gridfire {
}
void GraphEngine::generateJacobianMatrix(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho,
const SparsityPattern &sparsityPattern
@@ -1102,7 +1088,7 @@ namespace gridfire {
}
std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> GraphEngine::getSpeciesTimescales(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
@@ -1110,7 +1096,7 @@ namespace gridfire {
}
std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> GraphEngine::getSpeciesTimescales(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho,
const reaction::ReactionSet &activeReactions
@@ -1146,7 +1132,7 @@ namespace gridfire {
}
std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> GraphEngine::getSpeciesDestructionTimescales(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
@@ -1154,7 +1140,7 @@ namespace gridfire {
}
std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> GraphEngine::getSpeciesDestructionTimescales(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho,
const reaction::ReactionSet &activeReactions
@@ -1213,14 +1199,8 @@ namespace gridfire {
for (const auto& species : m_networkSpecies) {
if (!netIn.composition.contains(species)) {
baseUpdatedComposition.registerSpecies(species);
baseUpdatedComposition.setMassFraction(species, 0.0);
}
}
const bool didFinalize = baseUpdatedComposition.finalize(false);
if (!didFinalize) {
LOG_ERROR(m_logger, "Failed to finalize composition during update. Check input mass fractions for validity.");
throw std::runtime_error("Failed to finalize composition during network update.");
}
return baseUpdatedComposition;
}
@@ -1403,7 +1383,8 @@ namespace gridfire {
// We can pass a dummy comp and rho because reverse rates should only be calculated for strong reactions whose
// rates of progression do not depend on composition or density.
const double reverseRate = m_engine.calculateReverseRate(m_reaction, T9, 0.0, {});
const fourdst::composition::Composition dummyComp;
const double reverseRate = m_engine.calculateReverseRate(m_reaction, T9, 0.0, dummyComp);
ty[0] = reverseRate; // Store the reverse rate in the output vector
if (vx.size() > 0) {

13
src/lib/engine/procedures/construction.cpp
View File

@@ -6,12 +6,12 @@
#include <ranges>
#include <stdexcept>
#include <memory>
#include <cmath>
#include "fourdst/atomic/species.h"
#include "gridfire/reaction/reaction.h"
#include "gridfire/reaction/reaclib.h"
#include "fourdst/composition/composition.h"
#include "fourdst/composition/composition_abstract.h"
#include "fourdst/logging/logging.h"
@@ -146,11 +146,10 @@ namespace {
namespace gridfire {
using reaction::ReactionSet;
using reaction::Reaction;
using fourdst::composition::Composition;
using fourdst::atomic::Species;
ReactionSet build_nuclear_network(
const Composition& composition,
const fourdst::composition::CompositionAbstract &composition,
const rates::weak::WeakRateInterpolator &weakInterpolator,
BuildDepthType maxLayers,
NetworkConstructionFlags ReactionTypes
@@ -202,9 +201,9 @@ namespace gridfire {
// --- Step 3: Execute the layered network build using observing pointers ---
std::unordered_set<Species> availableSpecies;
for (const auto& entry : composition | std::views::values) {
if (entry.mass_fraction() > 0.0) {
availableSpecies.insert(entry.isotope());
for (const auto& sp : composition.getRegisteredSpecies()) {
if (composition.getMolarAbundance(sp) > 0.0) {
availableSpecies.insert(sp);
}
}

211
src/lib/engine/procedures/priming.cpp
View File

@@ -1,4 +1,7 @@
#include "gridfire/engine/procedures/priming.h"
#include "fourdst/atomic/species.h"
#include "fourdst/composition/utils.h"
#include "gridfire/engine/views/engine_priming.h"
#include "gridfire/engine/procedures/construction.h"
#include "gridfire/solver/solver.h"
@@ -11,30 +14,6 @@
#include "quill/LogMacros.h"
namespace {
// Create a dummy wrapper Composition to measure the unrestricted flow rate of species
class UnrestrictedComposition final : public fourdst::composition::Composition {
private:
const fourdst::atomic::Species& m_unrestrictedSpecies;
public:
explicit UnrestrictedComposition(const Composition& baseComposition, const fourdst::atomic::Species& unrestrictedSpecies) :
Composition(baseComposition),
m_unrestrictedSpecies(unrestrictedSpecies) {}
double getMolarAbundance(const fourdst::atomic::Species &species) const override {
if (species == m_unrestrictedSpecies) {
return 1.0; // Set to a high value to simulate unrestricted abundance
}
return Composition::getMolarAbundance(species);
}
double getMolarAbundance(const std::string &symbol) const override {
if (symbol == m_unrestrictedSpecies.name()) {
return 1.0; // Set to a high value to simulate unrestricted abundance
}
return Composition::getMolarAbundance(symbol);
}
};
bool isReactionIgnorable(
const gridfire::reaction::Reaction& reaction,
const std::optional<std::vector<gridfire::reaction::ReactionType>>& reactionsTypesToIgnore
@@ -117,18 +96,18 @@ namespace gridfire {
GraphEngine& engine,
const std::optional<std::vector<reaction::ReactionType>>& ignoredReactionTypes
) {
auto logger = fourdst::logging::LogManager::getInstance().getLogger("log");
auto logger = LogManager::getInstance().getLogger("log");
// --- Initial Setup ---
// Identify all species with zero initial mass fraction that need to be primed.
// Identify all species with zero initial abundance that need to be primed.
std::vector<Species> speciesToPrime;
for (const auto &entry: netIn.composition | std::views::values) {
if (entry.mass_fraction() == 0.0) {
speciesToPrime.push_back(entry.isotope());
for (const auto &[sp, y]: netIn.composition) {
if (y == 0.0) {
speciesToPrime.push_back(sp);
}
}
// sort primingSpecies by mass number, lightest to heaviest. This ensures we prime in a physically sensible order.
// Sort priming species by mass number, lightest to heaviest.
std::ranges::sort(speciesToPrime, [](const Species& a, const Species& b) {
return a.mass() < b.mass();
});
@@ -147,46 +126,37 @@ namespace gridfire {
const double T9 = netIn.temperature / 1e9;
const double rho = netIn.density;
const auto initialReactionSet = engine.getNetworkReactions();
report.status = PrimingReportStatus::FULL_SUCCESS;
report.success = true;
// Create a mutable map of the mass fractions that we will modify.
std::unordered_map<Species, double> currentMassFractions;
for (const auto& entry : netIn.composition | std::views::values) {
currentMassFractions[entry.isotope()] = entry.mass_fraction();
// Build full set of species
std::set<Species> allSpecies;
for (const auto &sp: netIn.composition | std::views::keys) {
allSpecies.insert(sp);
}
// Ensure all species to be primed exist in the map, initialized to zero.
for (const auto& entry : speciesToPrime) {
currentMassFractions[entry] = 0.0;
for (const auto& sp : speciesToPrime) {
allSpecies.insert(sp);
}
// Rebuild the engine with the full network to ensure all possible creation channels are available.
engine.rebuild(netIn.composition, NetworkBuildDepth::Full);
// --- STAGE 1: Calculation and Bookkeeping (No Modifications) ---
// In this stage, we calculate all required mass transfers but do not apply them yet.
// Initialize mutable molar abundances for all species
std::map<Species, double> molarAbundances;
for (const auto& sp : allSpecies) {
molarAbundances[sp] = netIn.composition.contains(sp) ? netIn.composition.getMolarAbundance(sp) : 0.0;
}
// A struct to hold the result of each individual priming calculation.
struct MassTransferRequest {
Species species_to_prime;
double equilibrium_mass_fraction;
std::vector<Species> reactants;
};
std::vector<MassTransferRequest> requests;
// --- Prime Each Species ---
// Since molar abundances are independent, we can directly calculate and apply changes
std::unordered_map<Species, double> totalMolarAbundanceChanges;
for (const auto& primingSpecies : speciesToPrime) {
// Create a temporary composition reflecting the current state for rate calculations.
Composition tempComp;
for(const auto& [sp, mf] : currentMassFractions) {
tempComp.registerSymbol(std::string(sp.name()));
tempComp.setMassFraction(sp, std::max(0.0, mf));
}
bool didFinalize = tempComp.finalize(true);
if (!didFinalize) {
LOG_ERROR(logger, "Failed to finalize temporary composition during priming.");
report.success = false;
report.status = PrimingReportStatus::FAILED_TO_FINALIZE_COMPOSITION;
continue;
Composition tempComp(allSpecies);
for (const auto& [sp, abundance] : molarAbundances) {
tempComp.setMolarAbundance(sp, abundance);
}
NetworkPrimingEngineView primer(primingSpecies, engine);
@@ -206,6 +176,9 @@ namespace gridfire {
ignoredReactionTypes
);
double equilibriumMolarAbundance = 0.0;
std::vector<Species> reactants;
if (destructionRateConstant > 1e-99) {
const double creationRate = calculateCreationRate(
primer,
@@ -215,10 +188,13 @@ namespace gridfire {
rho,
ignoredReactionTypes
);
double equilibriumMassFraction = (creationRate / destructionRateConstant) * primingSpecies.mass();
// ReSharper disable once CppDFAUnusedValue
if (std::isnan(equilibriumMassFraction)) equilibriumMassFraction = 0.0;
// Equilibrium: creation rate = destruction rate constant * molar abundance
equilibriumMolarAbundance = creationRate / destructionRateConstant;
if (std::isnan(equilibriumMolarAbundance)) {
equilibriumMolarAbundance = 0.0;
}
if (const reaction::Reaction* dominantChannel = findDominantCreationChannel(
primer,
@@ -228,102 +204,74 @@ namespace gridfire {
rho,
ignoredReactionTypes)
) {
// Store the request instead of applying it immediately.
requests.push_back({primingSpecies, equilibriumMassFraction, dominantChannel->reactants()});
reactants = dominantChannel->reactants();
} else {
LOG_TRACE_L1(logger, "Failed to find dominant creation channel for {}.", primingSpecies.name());
report.status = PrimingReportStatus::FAILED_TO_FIND_CREATION_CHANNEL;
reactants.clear(); // Use fallback
}
} else {
LOG_TRACE_L2(logger, "No destruction channel found for {}. Using fallback abundance.", primingSpecies.name());
// For species with no destruction, we can't calculate an equilibrium.
// We add a request with a tiny fallback abundance to ensure it exists in the network.
requests.push_back({primingSpecies, 1e-40, {}});
// For species with no destruction, use a tiny fallback abundance
equilibriumMolarAbundance = 1e-40;
}
}
// --- STAGE 2: Collective Scaling Based on Reactant Availability ---
// Now, we determine the total demand for each reactant and find a global scaling factor.
// Add the equilibrium molar abundance to the primed species
molarAbundances.at(primingSpecies) += equilibriumMolarAbundance;
totalMolarAbundanceChanges[primingSpecies] += equilibriumMolarAbundance;
std::unordered_map<Species, double> total_mass_debits;
for (const auto& req : requests) {
if (req.reactants.empty()) continue; // Skip fallbacks which don't consume reactants.
double totalReactantMass = 0.0;
for (const auto& reactant : req.reactants) {
totalReactantMass += reactant.mass();
}
if (totalReactantMass == 0.0) continue;
for (const auto& reactant : req.reactants) {
const double massToSubtract = req.equilibrium_mass_fraction * (reactant.mass() / totalReactantMass);
total_mass_debits[reactant] += massToSubtract;
}
}
double globalScalingFactor = 1.0;
for (const auto& [reactant, total_debit] : total_mass_debits) {
double availableMass;
if (currentMassFractions.contains(reactant)) {
availableMass = currentMassFractions.at(reactant);
} else {
availableMass = 0.0;
}
if (total_debit > availableMass && availableMass > 0) {
globalScalingFactor = std::min(globalScalingFactor, availableMass / total_debit);
}
}
if (globalScalingFactor < 1.0) {
LOG_WARNING(logger, "Priming was limited by reactant availability. All transfers will be scaled by {:.4e}", globalScalingFactor);
}
// --- STAGE 3: Application of Scaled Mass Transfers ---
// Finally, apply all the transfers, scaled by our global factor.
std::unordered_map<Species, double> totalMassFractionChanges;
for (const auto&[species_to_prime, equilibrium_mass_fraction, reactants] : requests) {
const double scaled_equilibrium_mf = equilibrium_mass_fraction * globalScalingFactor;
// Add the scaled mass to the primed species.
currentMassFractions.at(species_to_prime) += scaled_equilibrium_mf;
totalMassFractionChanges[species_to_prime] += scaled_equilibrium_mf;
// Subtract the scaled mass from the reactants.
// Subtract from reactants proportionally to their stoichiometry
if (!reactants.empty()) {
double totalReactantMass = 0.0;
for (const auto& reactant : reactants) {
totalReactantMass += reactant.mass();
}
if (totalReactantMass == 0.0) continue;
const double totalStoichiometry = static_cast<double>(reactants.size());
const double abundancePerReactant = equilibriumMolarAbundance / totalStoichiometry;
for (const auto& reactant : reactants) {
const double massToSubtract = scaled_equilibrium_mf * (reactant.mass() / totalReactantMass);
if (massToSubtract != 0) {
currentMassFractions.at(reactant) -= massToSubtract;
totalMassFractionChanges[reactant] -= massToSubtract;
LOG_TRACE_L1(logger, "Transferring {:.4e} molar abundance from {} to {} during priming.",
abundancePerReactant, reactant.name(), primingSpecies.name());
if (!molarAbundances.contains(reactant)) {
continue;
}
molarAbundances.at(reactant) -= abundancePerReactant;
totalMolarAbundanceChanges[reactant] -= abundancePerReactant;
// Ensure non-negative abundances
if (molarAbundances.at(reactant) < 0.0) {
LOG_WARNING(logger, "Species {} went negative during priming. Clamping to zero.", reactant.name());
totalMolarAbundanceChanges[reactant] += molarAbundances.at(reactant); // Adjust change to reflect clamp
molarAbundances.at(reactant) = 0.0;
}
}
}
}
// --- Final Composition Construction ---
std::vector<std::string> final_symbols;
std::vector<double> final_mass_fractions;
for(const auto& [species, mass_fraction] : currentMassFractions) {
final_symbols.emplace_back(species.name());
final_mass_fractions.push_back(std::max(0.0, mass_fraction)); // Ensure no negative mass fractions.
std::vector<Species> final_species;
std::vector<double> final_molar_abundances;
for (const auto& [species, abundance] : molarAbundances) {
final_species.emplace_back(species);
final_molar_abundances.push_back(std::max(0.0, abundance));
LOG_TRACE_L1(logger, "After priming, species {} has molar abundance {:.4e} (had {:0.4e} prior to priming).",
species.name(),
std::max(0.0, abundance),
netIn.composition.getMolarAbundance(species));
}
Composition primedComposition(final_symbols, final_mass_fractions, true);
Composition primedComposition(final_species, final_molar_abundances);
report.primedComposition = primedComposition;
for (const auto& [species, change] : totalMassFractionChanges) {
report.massFractionChanges.emplace_back(species, change);
// Convert molar abundance changes to mass fraction changes for reporting
for (const auto& [species, molarChange] : totalMolarAbundanceChanges) {
double massFractionChange = molarChange * species.mass();
report.massFractionChanges.emplace_back(species, massFractionChange);
}
// Restore the engine to its original, smaller network state.
engine.setNetworkReactions(initialReactionSet);
return report;
}
@@ -335,7 +283,8 @@ namespace gridfire {
const double rho,
const std::optional<std::vector<reaction::ReactionType>> &reactionTypesToIgnore
) {
const UnrestrictedComposition unrestrictedComp(composition, species); // Create a composition that simulates an enormous source abundance of the target species (getMolarAbundance(species) always returns 1.0)
Composition unrestrictedComp(composition);
unrestrictedComp.setMolarAbundance(species, 1.0);
double destructionRateConstant = 0.0;
for (const auto& reaction: engine.getNetworkReactions()) {
if (isReactionIgnorable(*reaction, reactionTypesToIgnore)) continue;

46
src/lib/engine/views/engine_adaptive.cpp
View File

@@ -91,12 +91,6 @@ namespace gridfire {
updatedNetIn.composition = baseUpdatedComposition;
bool didFinalize = updatedNetIn.composition.finalize(false);
if (!didFinalize) {
LOG_ERROR(m_logger, "Failed to finalize composition during adaptive engine view update. Check input mass fractions for validity.");
throw std::runtime_error("Failed to finalize composition during adaptive engine view update.");
}
LOG_TRACE_L1(m_logger, "Updating AdaptiveEngineView with new network input...");
auto [allFlows, composition] = calculateAllReactionFlows(updatedNetIn);
@@ -148,7 +142,7 @@ namespace gridfire {
}
std::expected<StepDerivatives<double>, expectations::StaleEngineError> AdaptiveEngineView::calculateRHSAndEnergy(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
@@ -164,7 +158,7 @@ namespace gridfire {
}
EnergyDerivatives AdaptiveEngineView::calculateEpsDerivatives(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
@@ -173,7 +167,7 @@ namespace gridfire {
}
void AdaptiveEngineView::generateJacobianMatrix(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
@@ -181,7 +175,7 @@ namespace gridfire {
}
void AdaptiveEngineView::generateJacobianMatrix(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho,
const std::vector<Species> &activeSpecies
@@ -192,7 +186,7 @@ namespace gridfire {
}
void AdaptiveEngineView::generateJacobianMatrix(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho,
const SparsityPattern &sparsityPattern
@@ -225,7 +219,7 @@ namespace gridfire {
double AdaptiveEngineView::calculateMolarReactionFlow(
const reaction::Reaction &reaction,
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
@@ -249,7 +243,7 @@ namespace gridfire {
}
std::expected<std::unordered_map<Species, double>, expectations::StaleEngineError> AdaptiveEngineView::getSpeciesTimescales(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
@@ -276,7 +270,7 @@ namespace gridfire {
std::expected<std::unordered_map<Species, double>, expectations::StaleEngineError>
AdaptiveEngineView::getSpeciesDestructionTimescales(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
@@ -308,8 +302,8 @@ namespace gridfire {
std::vector<double> AdaptiveEngineView::mapNetInToMolarAbundanceVector(const NetIn &netIn) const {
std::vector<double> Y(m_activeSpecies.size(), 0.0); // Initialize with zeros
for (const auto& [symbol, entry] : netIn.composition) {
Y[getSpeciesIndex(entry.isotope())] = netIn.composition.getMolarAbundance(symbol); // Map species to their molar abundance
for (const auto& [species, y] : netIn.composition) {
Y[getSpeciesIndex(species)] = y; // Map species to their molar abundance
}
return Y; // Return the vector of molar abundances
}
@@ -319,14 +313,13 @@ namespace gridfire {
}
fourdst::composition::Composition AdaptiveEngineView::collectComposition(
fourdst::composition::Composition &comp
fourdst::composition::CompositionAbstract &comp
) const {
fourdst::composition::Composition result = m_baseEngine.collectComposition(comp); // Step one is to bubble the results from lower levels of the engine chain up
for (const auto& species : m_activeSpecies) {
if (!result.hasSpecies(species)) {
if (!result.contains(species)) {
result.registerSpecies(species);
result.setMassFraction(species, 0.0);
}
}
@@ -396,9 +389,8 @@ namespace gridfire {
for (const auto& species: fullSpeciesList) {
if (!netIn.composition.contains(species)) {
LOG_TRACE_L2(m_logger, "Species '{}' not found in composition. Setting abundance to 0.0.", species.name());
LOG_TRACE_L2(m_logger, "Species '{}' not found in composition. Registering", species.name());
composition.registerSpecies(species);
composition.setMassFraction(species, 0.0);
}
}
@@ -411,7 +403,7 @@ namespace gridfire {
for (const auto& reaction : fullReactionSet) {
const double flow = m_baseEngine.calculateMolarReactionFlow(*reaction, composition, T9, rho);
reactionFlows.push_back({reaction.get(), flow});
LOG_TRACE_L1(m_logger, "Reaction '{}' has flow rate: {:0.3E} [mol/s/g]", reaction->id(), flow);
LOG_TRACE_L3(m_logger, "Reaction '{}' has flow rate: {:0.3E} [mol/s/g]", reaction->id(), flow);
}
return {reactionFlows, composition};
}
@@ -586,13 +578,13 @@ namespace gridfire {
const double mue = 5.0e-3 * std::pow(rho * Ye, 1.0 / 3.0) + 0.5 * T9;
std::unordered_map<Species, double> speciesMassMap;
for (const auto &entry: comp | std::views::values) {
speciesMassMap[entry.isotope()] = entry.isotope().mass();
for (const auto &sp: comp | std::views::keys) {
speciesMassMap[sp] = sp.mass();
}
std::unordered_map<size_t, Species> speciesIndexMap;
for (const auto& entry: comp | std::views::values) {
size_t distance = std::distance(speciesMassMap.begin(), speciesMassMap.find(entry.isotope()));
speciesIndexMap.emplace(distance, entry.isotope());
for (const auto& sp: comp | std::views::keys) {
size_t distance = std::distance(speciesMassMap.begin(), speciesMassMap.find(sp));
speciesIndexMap.emplace(distance, sp);
}
double rate = reaction->calculate_rate(T9, rho, Ye, mue, Y, speciesIndexMap);
if (rate > maxSpeciesConsumptionRate) {

63
src/lib/engine/views/engine_defined.cpp
View File

@@ -1,7 +1,9 @@
#include "gridfire/engine/views/engine_defined.h"
#include "gridfire/engine/engine_graph.h"
#include <ranges>
#include "fourdst/atomic/species.h"
#include "fourdst/atomic/atomicSpecies.h"
#include "fourdst/composition/decorators/composition_masked.h"
#include "quill/LogMacros.h"
@@ -13,20 +15,7 @@
#include <unordered_map>
#include <utility>
namespace {
class MaskedComposition final : public fourdst::composition::Composition {
private:
std::set<fourdst::atomic::Species> m_activeSpecies;
public:
MaskedComposition(const Composition& baseComposition, const std::set<fourdst::atomic::Species>& activeSpecies) :
Composition(baseComposition),
m_activeSpecies(activeSpecies) {}
bool contains(const fourdst::atomic::Species& species) const override {
return Composition::contains(species) && m_activeSpecies.contains(species);
}
};
}
#include "fourdst/composition/exceptions/exceptions_composition.h"
namespace gridfire {
using fourdst::atomic::Species;
@@ -52,13 +41,13 @@ namespace gridfire {
}
std::expected<StepDerivatives<double>, expectations::StaleEngineError> DefinedEngineView::calculateRHSAndEnergy(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
validateNetworkState();
const MaskedComposition masked(comp, m_activeSpecies);
const fourdst::composition::MaskedComposition masked(comp, m_activeSpecies);
const auto result = m_baseEngine.calculateRHSAndEnergy(masked, T9, rho, m_activeReactions);
if (!result) {
@@ -69,19 +58,19 @@ namespace gridfire {
}
EnergyDerivatives DefinedEngineView::calculateEpsDerivatives(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
validateNetworkState();
const MaskedComposition masked(comp, m_activeSpecies);
const fourdst::composition::MaskedComposition masked(comp, m_activeSpecies);
return m_baseEngine.calculateEpsDerivatives(masked, T9, rho, m_activeReactions);
}
void DefinedEngineView::generateJacobianMatrix(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
@@ -89,12 +78,12 @@ namespace gridfire {
if (!m_activeSpeciesVectorCache.has_value()) {
m_activeSpeciesVectorCache = std::vector<Species>(m_activeSpecies.begin(), m_activeSpecies.end());
}
const MaskedComposition masked(comp, m_activeSpecies);
const fourdst::composition::MaskedComposition masked(comp, m_activeSpecies);
m_baseEngine.generateJacobianMatrix(masked, T9, rho, m_activeSpeciesVectorCache.value());
}
void DefinedEngineView::generateJacobianMatrix(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho,
const std::vector<fourdst::atomic::Species> &activeSpecies
@@ -106,18 +95,18 @@ namespace gridfire {
activeSpecies.end()
);
const MaskedComposition masked(comp, activeSpeciesSet);
const fourdst::composition::MaskedComposition masked(comp, activeSpeciesSet);
m_baseEngine.generateJacobianMatrix(masked, T9, rho, activeSpecies);
}
void DefinedEngineView::generateJacobianMatrix(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho,
const SparsityPattern &sparsityPattern
) const {
validateNetworkState();
const MaskedComposition masked(comp, m_activeSpecies);
const fourdst::composition::MaskedComposition masked(comp, m_activeSpecies);
m_baseEngine.generateJacobianMatrix(masked, T9, rho, sparsityPattern);
}
@@ -170,7 +159,7 @@ namespace gridfire {
double DefinedEngineView::calculateMolarReactionFlow(
const reaction::Reaction &reaction,
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
@@ -182,7 +171,7 @@ namespace gridfire {
throw std::runtime_error("Reaction not found in active reactions: " + std::string(reaction.id()));
}
const MaskedComposition masked(comp, m_activeSpecies);
const fourdst::composition::MaskedComposition masked(comp, m_activeSpecies);
return m_baseEngine.calculateMolarReactionFlow(reaction, masked, T9, rho);
}
@@ -202,12 +191,12 @@ namespace gridfire {
}
std::expected<std::unordered_map<Species, double>, expectations::StaleEngineError> DefinedEngineView::getSpeciesTimescales(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
validateNetworkState();
const MaskedComposition masked(comp, m_activeSpecies);
const fourdst::composition::MaskedComposition masked(comp, m_activeSpecies);
const auto result = m_baseEngine.getSpeciesTimescales(masked, T9, rho, m_activeReactions);
if (!result) {
@@ -225,12 +214,12 @@ namespace gridfire {
}
std::expected<std::unordered_map<Species, double>, expectations::StaleEngineError> DefinedEngineView::getSpeciesDestructionTimescales(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
validateNetworkState();
const MaskedComposition masked(comp, m_activeSpecies);
const fourdst::composition::MaskedComposition masked(comp, m_activeSpecies);
const auto result = m_baseEngine.getSpeciesDestructionTimescales(masked, T9, rho, m_activeReactions);
@@ -257,7 +246,6 @@ namespace gridfire {
return m_baseEngine.isStale(netIn);
}
void DefinedEngineView::setScreeningModel(const screening::ScreeningType model) {
m_baseEngine.setScreeningModel(model);
}
@@ -282,10 +270,10 @@ namespace gridfire {
std::vector<double> DefinedEngineView::mapNetInToMolarAbundanceVector(const NetIn &netIn) const {
std::vector<double> Y(m_activeSpecies.size(), 0.0); // Initialize with zeros
for (const auto& [symbol, entry] : netIn.composition) {
auto it = std::ranges::find(m_activeSpecies, entry.isotope());
for (const auto& [sp, y] : netIn.composition) {
auto it = std::ranges::find(m_activeSpecies, sp);
if (it != m_activeSpecies.end()) {
Y[getSpeciesIndex(entry.isotope())] = netIn.composition.getMolarAbundance(symbol); // Map species to their molar abundance
Y[getSpeciesIndex(sp)] = y; // Map species to their molar abundance
}
}
return Y; // Return the vector of molar abundances
@@ -296,14 +284,13 @@ namespace gridfire {
}
fourdst::composition::Composition DefinedEngineView::collectComposition(
fourdst::composition::Composition &comp
fourdst::composition::CompositionAbstract &comp
) const {
fourdst::composition::Composition result = m_baseEngine.collectComposition(comp);
for (const auto& species : m_activeSpecies) {
if (!result.hasSpecies(species)) {
if (!result.contains(species)) {
result.registerSpecies(species);
result.setMassFraction(species, 0.0);
}
}
return result;

127
src/lib/engine/views/engine_multiscale.cpp
View File

@@ -1,7 +1,6 @@
#include "gridfire/engine/views/engine_multiscale.h"
#include "gridfire/exceptions/error_engine.h"
#include "gridfire/engine/procedures/priming.h"
#include "gridfire/utils/general_composition.h"
#include <stdexcept>
#include <vector>
@@ -14,6 +13,7 @@
#include <algorithm>
#include "fourdst/atomic/species.h"
#include "quill/LogMacros.h"
#include "quill/Logger.h"
@@ -182,7 +182,7 @@ namespace gridfire {
}
std::expected<StepDerivatives<double>, expectations::StaleEngineError> MultiscalePartitioningEngineView::calculateRHSAndEnergy(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
@@ -200,7 +200,7 @@ namespace gridfire {
}
EnergyDerivatives MultiscalePartitioningEngineView::calculateEpsDerivatives(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
@@ -208,7 +208,7 @@ namespace gridfire {
}
void MultiscalePartitioningEngineView::generateJacobianMatrix(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
@@ -217,7 +217,7 @@ namespace gridfire {
}
void MultiscalePartitioningEngineView::generateJacobianMatrix(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho,
const std::vector<Species> &activeSpecies
@@ -253,7 +253,7 @@ namespace gridfire {
}
void MultiscalePartitioningEngineView::generateJacobianMatrix(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho,
const SparsityPattern &sparsityPattern
@@ -288,34 +288,20 @@ namespace gridfire {
double MultiscalePartitioningEngineView::calculateMolarReactionFlow(
const reaction::Reaction &reaction,
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
// Fix the algebraic species to the equilibrium abundances we calculate.
fourdst::composition::Composition comp_mutable = comp;
const bool didFinalize = comp_mutable.finalize(false);
if (!didFinalize) {
LOG_ERROR(m_logger, "Failed to finalize composition before setting algebraic species abundances.");
m_logger->flush_log();
throw std::runtime_error("Failed to finalize composition before setting algebraic species abundances.");
fourdst::composition::Composition comp_mutable;
for (const auto& species : comp.getRegisteredSpecies()) {
comp_mutable.registerSpecies(species);
comp_mutable.setMolarAbundance(species, comp.getMolarAbundance(species));
}
for (const auto& species : m_algebraic_species) {
const double Yi = m_algebraic_abundances.at(species);
double Xi = utils::massFractionFromMolarAbundanceAndComposition(comp_mutable, species, Yi);
comp_mutable.setMassFraction(species, Xi); // Convert Yi (mol/g) to Xi (mass fraction)
if (!comp_mutable.finalize(false)) {
LOG_ERROR(m_logger, "Failed to finalize composition after setting algebraic species abundance for species '{}'.", species.name());
m_logger->flush_log();
throw std::runtime_error("Failed to finalize composition after setting algebraic species abundance for species: " + std::string(species.name()));
}
}
if (!comp_mutable.finalize()) {
LOG_ERROR(m_logger, "Failed to finalize composition after setting algebraic species abundances.");
m_logger->flush_log();
throw std::runtime_error("Failed to finalize composition after setting algebraic species abundances.");
comp_mutable.setMolarAbundance(species, Yi);
}
return m_baseEngine.calculateMolarReactionFlow(reaction, comp_mutable, T9, rho);
@@ -331,7 +317,7 @@ namespace gridfire {
}
std::expected<std::unordered_map<Species, double>, expectations::StaleEngineError> MultiscalePartitioningEngineView::getSpeciesTimescales(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
@@ -348,7 +334,7 @@ namespace gridfire {
std::expected<std::unordered_map<Species, double>, expectations::StaleEngineError>
MultiscalePartitioningEngineView::getSpeciesDestructionTimescales(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
const double T9,
const double rho
) const {
@@ -371,7 +357,7 @@ namespace gridfire {
const fourdst::composition::Composition equilibratedComposition = equilibrateNetwork(baseUpdatedNetIn);
std::unordered_map<Species, double> algebraicAbundances;
for (const auto& species : m_algebraic_species) {
algebraicAbundances[species] = equilibratedComposition.getMolarAbundance(species);
algebraicAbundances.emplace(species, equilibratedComposition.getMolarAbundance(species));
}
m_algebraic_abundances = std::move(algebraicAbundances);
@@ -791,8 +777,8 @@ namespace gridfire {
std::vector<double> MultiscalePartitioningEngineView::mapNetInToMolarAbundanceVector(const NetIn &netIn) const {
std::vector<double> Y(m_dynamic_species.size(), 0.0); // Initialize with zeros
for (const auto& [symbol, entry] : netIn.composition) {
Y[getSpeciesIndex(entry.isotope())] = netIn.composition.getMolarAbundance(symbol); // Map species to their molar abundance
for (const auto& [sp, y] : netIn.composition) {
Y[getSpeciesIndex(sp)] = y; // Map species to their molar abundance
}
return Y; // Return the vector of molar abundances
}
@@ -826,20 +812,12 @@ namespace gridfire {
partitionNetwork(comp, T9, rho);
fourdst::composition::Composition qseComposition = solveQSEAbundances(comp, T9, rho);
for (const auto &symbol: qseComposition | std::views::keys) {
const double speciesMassFraction = qseComposition.getMassFraction(symbol);
if (speciesMassFraction < 0.0 && std::abs(speciesMassFraction) < 1e-20) {
qseComposition.setMassFraction(symbol, 0.0); // Avoid negative mass fractions
for (const auto &[sp, y]: qseComposition) {
if (y < 0.0 && std::abs(y) < 1e-20) {
qseComposition.setMolarAbundance(sp, 0.0); // Avoid negative mass fractions
}
}
bool didFinalize = qseComposition.finalize(true);
if (!didFinalize) {
LOG_ERROR(m_logger, "Failed to finalize composition after solving QSE abundances.");
m_logger->flush_log();
throw std::runtime_error("Failed to finalize composition after solving QSE abundances.");
}
return qseComposition;
}
@@ -875,18 +853,12 @@ namespace gridfire {
}
fourdst::composition::Composition MultiscalePartitioningEngineView::collectComposition(
fourdst::composition::Composition &comp
fourdst::composition::CompositionAbstract &comp
) const {
fourdst::composition::Composition result = m_baseEngine.collectComposition(comp);
bool didFinalize = result.finalize(false);
if (!didFinalize) {
std::string msg = "Failed to finalize collected composition from MultiscalePartitioningEngine view after calling base engines collectComposition method.";
LOG_ERROR(m_logger, "{}", msg);
throw exceptions::BadCollectionError(msg);
}
std::map<Species, double> Ym; // Use an ordered map here so that this is ordered by atomic mass (which is the </> comparator for Species)
for (const auto& [speciesName, entry] : result) {
Ym.emplace(entry.isotope(), result.getMolarAbundance(speciesName));
for (const auto& [sp, y] : result) {
Ym.emplace(sp, y);
}
for (const auto& [species, Yi] : m_algebraic_abundances) {
if (!Ym.contains(species)) {
@@ -894,21 +866,16 @@ namespace gridfire {
}
Ym.at(species) = Yi;
}
std::vector<double> M;
std::vector<double> Y;
std::vector<std::string> speciesNames;
M.reserve(Ym.size());
Y.reserve(Ym.size());
for (const auto& [species, Yi] : Ym) {
M.emplace_back(species.mass());
Y.emplace_back(Yi);
speciesNames.emplace_back(species.name());
}
std::vector<double> X = utils::massFractionFromMolarAbundanceAndMolarMass(Y, M);
return fourdst::composition::Composition(speciesNames, X);
return {speciesNames, Y};
}
size_t MultiscalePartitioningEngineView::getSpeciesIndex(const Species &species) const {
@@ -1254,20 +1221,13 @@ namespace gridfire {
for (const auto& species: algebraic_species) {
if (!normalized_composition.contains(species)) {
normalized_composition.registerSpecies(species);
normalized_composition.setMassFraction(species, 0.0);
}
}
for (const auto& species: seed_species) {
if (!normalized_composition.contains(species)) {
normalized_composition.registerSpecies(species);
normalized_composition.setMassFraction(species, 0.0);
}
}
bool normCompFinalizedOkay = normalized_composition.finalize(true);
if (!normCompFinalizedOkay) {
LOG_ERROR(m_logger, "Failed to finalize composition before QSE solve.");
throw std::runtime_error("Failed to finalize composition before QSE solve.");
}
Eigen::VectorXd Y_scale(algebraic_species.size());
Eigen::VectorXd v_initial(algebraic_species.size());
@@ -1319,20 +1279,13 @@ namespace gridfire {
Y_final_qse(i)
);
// double Xi = Y_final_qse(i) * species.mass(); // Convert from molar abundance to mass fraction
double Xi = utils::massFractionFromMolarAbundanceAndComposition(normalized_composition, species, Y_final_qse(i));
if (!outputComposition.hasSpecies(species)) {
if (!outputComposition.contains (species)) {
outputComposition.registerSpecies(species);
}
outputComposition.setMassFraction(species, Xi);
outputComposition.setMolarAbundance(species, Y_final_qse(i));
i++;
}
}
bool didFinalize = outputComposition.finalize(false);
if (!didFinalize) {
LOG_ERROR(m_logger, "Failed to finalize composition after solving QSE abundances.");
m_logger->flush_log();
throw std::runtime_error("Failed to finalize composition after solving QSE abundances.");
}
return outputComposition;
}
@@ -1540,22 +1493,11 @@ namespace gridfire {
Eigen::VectorXd y_qse = v_qse.array().exp(); // Convert to physical abundances using exponential scaling
for (const auto& species: m_qse_solve_species) {
if (!comp_trial.hasSymbol(std::string(species.name()))) {
if (!comp_trial.contains(std::string(species.name()))) {
comp_trial.registerSpecies(species);
}
auto index = static_cast<long>(m_qse_solve_species_index_map.at(species));
const double molarAbundance = y_qse[index];
double massFraction = utils::massFractionFromMolarAbundanceAndComposition(m_initial_comp, species, molarAbundance);
comp_trial.setMassFraction(species, massFraction);
}
const bool didFinalize = comp_trial.finalize(false);
if (!didFinalize) {
LOG_TRACE_L1(m_view->m_logger, "While evaluating the functor, failed to finalize composition. This is likely because the solver took a step outside of physical abundances. This is not an error; rather, the solver will be told to take a different step.");
f_qse.resize(static_cast<long>(m_qse_solve_species.size()));
f_qse.setConstant(1.0e20); // Return a large residual to indicate failure
return 0;
comp_trial.setMolarAbundance(species, y_qse[index]);
}
const auto result = m_view->getBaseEngine().calculateRHSAndEnergy(comp_trial, m_T9, m_rho);
@@ -1608,20 +1550,11 @@ namespace gridfire {
Eigen::VectorXd y_qse = v_qse.array().exp(); // Convert to physical abundances using exponential scaling
for (const auto& species: m_qse_solve_species) {
if (!comp_trial.hasSymbol(std::string(species.name()))) {
if (!comp_trial.contains(std::string(species.name()))) {
comp_trial.registerSpecies(species);
}
const double molarAbundance = y_qse[static_cast<long>(m_qse_solve_species_index_map.at(species))];
double massFraction = utils::massFractionFromMolarAbundanceAndComposition(m_initial_comp, species, molarAbundance);
comp_trial.setMassFraction(species, massFraction);
}
const bool didFinalize = comp_trial.finalize(false);
if (!didFinalize) {
LOG_TRACE_L1(m_view->m_logger, "While evaluating the Jacobian, failed to finalize composition. This is likely because the solver took a step outside of physical abundances. This is not an error; rather, the solver will be told to take a different step. Returning Identity");
J_qse.resize(static_cast<long>(m_qse_solve_species.size()), static_cast<long>(m_qse_solve_species.size()));
J_qse.setIdentity();
return 0;
comp_trial.setMolarAbundance(species, molarAbundance);
}
std::vector<Species> qse_species_vector(m_qse_solve_species.begin(), m_qse_solve_species.end());

3
src/lib/engine/views/engine_priming.cpp
View File

@@ -1,7 +1,7 @@
#include "gridfire/engine/views/engine_priming.h"
#include "gridfire/solver/solver.h"
#include "fourdst/composition/species.h"
#include "fourdst/atomic/species.h"
#include "quill/LogMacros.h"
#include "quill/Logger.h"
@@ -9,7 +9,6 @@
#include <vector>
#include <string>
#include <unordered_set>
#include <stdexcept>
#include <unordered_map>

8
src/lib/io/generative/python.cpp
View File

@@ -1,11 +1,9 @@
#include "gridfire/io/generative/python.h"
#include "fourdst/composition/species.h"
#include "fourdst/atomic/atomicSpecies.h"
#include <string>
#include <vector>
#include <format>
#include <algorithm>
#include <ranges>
#include <fstream>
#include <optional>
@@ -14,7 +12,7 @@
namespace {
template <typename T>
std::string join(std::vector<T> arr, std::string delim) {
std::string join(std::vector<T> arr, const std::string& delim) {
if (arr.empty()) return {};
size_t total = delim.size() * (arr.size() - 1);
@@ -142,7 +140,7 @@ namespace gridfire::io::gen {
std::string exportEngineToPy(const gridfire::DynamicEngine& engine) {
auto reactions = engine.getNetworkReactions();
std::vector<std::string> functions;
functions.push_back(R"(import numpy as np
functions.emplace_back(R"(import numpy as np
from typing import Dict, List, Tuple, Callable)");
functions.push_back(tfFunc);
for (const auto& reaction : reactions) {

36
src/lib/network.cpp
View File

@@ -23,43 +23,7 @@
#include <ranges>
#include "quill/LogMacros.h"
namespace gridfire {
std::vector<double> NetIn::MolarAbundance() const {
std::vector <double> y;
y.reserve(composition.getRegisteredSymbols().size());
const auto [fst, snd] = composition.getComposition();
for (const auto &name: fst | std::views::keys) {
y.push_back(composition.getMolarAbundance(name));
}
return y;
}
Network::Network(const NetworkFormat format) :
m_config(fourdst::config::Config::getInstance()),
m_logManager(fourdst::logging::LogManager::getInstance()),
m_logger(m_logManager.getLogger("log")),
m_format(format),
m_constants(fourdst::constant::Constants::getInstance()){
if (format == NetworkFormat::UNKNOWN) {
LOG_ERROR(m_logger, "nuclearNetwork::Network::Network() called with UNKNOWN format");
m_logger->flush_log();
throw std::runtime_error("nuclearNetwork::Network::Network() called with UNKNOWN format");
}
}
NetworkFormat Network::getFormat() const {
return m_format;
}
NetworkFormat Network::setFormat(const NetworkFormat format) {
const NetworkFormat oldFormat = m_format;
m_format = format;
return oldFormat;
}
// Trim whitespace from both ends of a string
std::string trim_whitespace(const std::string& str) {
auto startIt = str.begin();

4
src/lib/partition/partition_ground.cpp
View File

@@ -2,8 +2,8 @@
#include <ranges>
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/composition/species.h"
#include "fourdst/atomic/atomicSpecies.h"
#include "fourdst/atomic/species.h"
#include "quill/LogMacros.h"

3
src/lib/policy/chains.cpp
View File

@@ -101,8 +101,7 @@ namespace gridfire::policy {
"d(p,g)he3",
"he4(he3,g)be7",
"be7(p,g)b8",
"b8(,e+)be8",
"be8(,a)he4"
"b8(,e+ a)he4"
}, 0.001) {}
std::unique_ptr<ReactionChainPolicy> ProtonProtonIChainPolicy::clone() const {

1
src/lib/policy/policy_logical.cpp
View File

@@ -33,6 +33,7 @@ namespace gridfire::policy {
return std::make_unique<MultiReactionChainPolicy>(
[this]() {
std::vector<std::unique_ptr<ReactionChainPolicy>> chain_policies;
chain_policies.reserve(m_chain_policies.size());
for (const auto &ch : m_chain_policies) {
chain_policies.push_back(ch->clone());
}

41
src/lib/policy/stellar_policy.cpp
View File

@@ -6,12 +6,28 @@
#include "gridfire/engine/engine_graph.h"
#include "gridfire/engine/views/engine_views.h"
#include "fourdst/composition/exceptions/exceptions_composition.h"
#include "fourdst/atomic/species.h"
#include "fourdst/composition/utils.h"
namespace {
std::set<fourdst::atomic::Species> initialize_seed_species() {
return {
fourdst::atomic::H_1,
fourdst::atomic::He_3,
fourdst::atomic::He_4,
fourdst::atomic::C_12,
fourdst::atomic::N_14,
fourdst::atomic::O_16,
fourdst::atomic::Ne_20,
fourdst::atomic::Mg_24
};
}
}
namespace gridfire::policy {
MainSequencePolicy::MainSequencePolicy(const fourdst::composition::Composition& composition) {
MainSequencePolicy::MainSequencePolicy(const fourdst::composition::Composition& composition) : m_seed_species(initialize_seed_species()) {
for (const auto& species : m_seed_species) {
if (!composition.hasSpecies(species)) {
if (!composition.contains(species)) {
throw exceptions::MissingSeedSpeciesError("Cannot initialize MainSequencePolicy: Required Seed species " + std::string(species.name()) + " is missing from the provided composition.");
}
}
@@ -19,22 +35,14 @@ namespace gridfire::policy {
m_partition_function = build_partition_function();
}
MainSequencePolicy::MainSequencePolicy(std::vector<fourdst::atomic::Species> seed_species, std::vector<double> mass_fractions) {
MainSequencePolicy::MainSequencePolicy(std::vector<fourdst::atomic::Species> seed_species, const std::vector<double> &mass_fractions) {
for (const auto& species : m_seed_species) {
if (std::ranges::find(seed_species, species) == seed_species.end()) {
throw exceptions::MissingSeedSpeciesError("Cannot initialize MainSequencePolicy: Required Seed species " + std::string(species.name()) + " is missing from the provided composition.");
}
}
for (const auto& [species, x] : std::views::zip(seed_species, mass_fractions)) {
m_initializing_composition.registerSpecies(species);
m_initializing_composition.setMassFraction(species, x);
}
const bool didFinalize = m_initializing_composition.finalize(true);
if (!didFinalize) {
throw fourdst::composition::exceptions::CompositionNotFinalizedError("Failed to finalize initial composition for MainSequencePolicy.");
}
m_initializing_composition = fourdst::composition::buildCompositionFromMassFractions(seed_species, mass_fractions);
m_partition_function = build_partition_function();
}
@@ -45,6 +53,11 @@ namespace gridfire::policy {
m_network_stack.emplace_back(
std::make_unique<GraphEngine>(m_initializing_composition, *m_partition_function, NetworkBuildDepth::ThirdOrder, NetworkConstructionFlags::DEFAULT)
);
auto& graphRepr = dynamic_cast<GraphEngine&>(*m_network_stack.back().get());
graphRepr.setUseReverseReactions(false);
m_network_stack.emplace_back(
std::make_unique<MultiscalePartitioningEngineView>(*m_network_stack.back().get())
);
@@ -85,7 +98,7 @@ namespace gridfire::policy {
inline NetworkPolicyStatus MainSequencePolicy::check_status() const {
for (const auto& species : m_seed_species) {
if (!m_initializing_composition.hasSpecies(species)) {
if (!m_initializing_composition.contains(species)) {
return NetworkPolicyStatus::MISSING_KEY_SPECIES;
}
}

4
src/lib/reaction/reaclib.cpp
View File

@@ -1,5 +1,5 @@
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/composition/species.h"
#include "fourdst/atomic/atomicSpecies.h"
#include "fourdst/atomic/species.h"
#include "gridfire/reaction/reaclib.h"
#include "gridfire/reaction/reactions_data.h"

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

@@ -10,7 +10,7 @@
#include "quill/LogMacros.h"
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/atomic/atomicSpecies.h"
#include "xxhash64.h"

3
src/lib/reaction/weak/weak.cpp
View File

@@ -1,8 +1,6 @@
#include "gridfire/reaction/weak/weak_rate_library.h"
#include "gridfire/reaction/weak/weak.h"
#include "fourdst/composition/species.h"
#include <array>
#include <ranges>
#include <unordered_map>
@@ -13,6 +11,7 @@
#include "gridfire/reaction/weak/weak_interpolator.h"
#include "xxhash64.h"
#include "fourdst/atomic/species.h"
namespace {

3
src/lib/reaction/weak/weak_interpolator.cpp
View File

@@ -4,7 +4,6 @@
#include "gridfire/utils/hashing.h"
#include <algorithm>
#include <map>
#include <set>
#include <unordered_map>
#include <vector>
@@ -12,7 +11,7 @@
#include <expected>
#include <ranges>
#include "fourdst/composition/species.h"
#include "fourdst/atomic/species.h"
#include "quill/LogMacros.h"

2
src/lib/screening/screening_bare.cpp
View File

@@ -1,6 +1,6 @@
#include "gridfire/screening/screening_bare.h"
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/atomic/atomicSpecies.h"
#include "cppad/cppad.hpp"

2
src/lib/screening/screening_weak.cpp
View File

@@ -1,6 +1,6 @@
#include "gridfire/screening/screening_weak.h"
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/atomic/atomicSpecies.h"
#include "cppad/cppad.hpp"

207
src/lib/solver/strategies/CVODE_solver_strategy.cpp
View File

@@ -20,7 +20,6 @@
#include "gridfire/engine/engine_graph.h"
#include "gridfire/solver/strategies/triggers/engine_partitioning_trigger.h"
#include "gridfire/trigger/procedures/trigger_pprint.h"
#include "gridfire/utils/general_composition.h"
namespace {
std::unordered_map<int, std::string> cvode_ret_code_map {
@@ -68,9 +67,9 @@ namespace {
N_Vector init_sun_vector(uint64_t size, SUNContext sun_ctx) {
#ifdef SUNDIALS_HAVE_OPENMP
N_Vector vec = N_VNew_OpenMP(size, 0, sun_ctx);
const N_Vector vec = N_VNew_OpenMP(size, 0, sun_ctx);
#elif SUNDIALS_HAVE_PTHREADS
N_Vector vec = N_VNew_Pthreads(size, sun_ctx);
const N_Vector vec = N_VNew_Pthreads(size, sun_ctx);
#else
const N_Vector vec = N_VNew_Serial(static_cast<long long>(size), sun_ctx);
#endif
@@ -273,6 +272,8 @@ namespace gridfire::solver {
total_convergence_failures += nlcfails;
total_steps += n_steps;
sunrealtype* end_of_step_abundances = N_VGetArrayPointer(ctx.state);
LOG_INFO(
m_logger,
"Engine Update Triggered at time {} ({} update{} triggered). Current total specific energy {} [erg/g]",
@@ -284,7 +285,7 @@ namespace gridfire::solver {
fourdst::composition::Composition temp_comp;
std::vector<double> mass_fractions;
long int num_species_at_stop = m_engine.getNetworkSpecies().size();
auto num_species_at_stop = static_cast<long int>(m_engine.getNetworkSpecies().size());
if (num_species_at_stop > m_Y->ops->nvgetlength(m_Y) - 1) {
LOG_ERROR(
@@ -296,29 +297,126 @@ namespace gridfire::solver {
throw std::runtime_error("Number of species at engine update exceeds the number of species in the CVODE solver. This should never happen.");
}
mass_fractions.reserve(num_species_at_stop);
for (const auto& species: m_engine.getNetworkSpecies()) {
const size_t sid = m_engine.getSpeciesIndex(species);
temp_comp.registerSpecies(species);
double y = end_of_step_abundances[sid];
if (y > 0.0) {
temp_comp.setMolarAbundance(species, y);
}
}
#ifndef NDEBUG
for (long int i = 0; i < num_species_at_stop; ++i) {
const auto& species = m_engine.getNetworkSpecies()[i];
temp_comp.registerSpecies(species);
mass_fractions.push_back(y_data[i] * species.mass()); // Convert from molar abundance to mass fraction
}
temp_comp.setMassFraction(m_engine.getNetworkSpecies(), mass_fractions);
bool didFinalize = temp_comp.finalize(true);
if (!didFinalize) {
LOG_ERROR(m_logger, "Failed to finalize composition during engine update. Check input mass fractions for validity.");
throw std::runtime_error("Failed to finalize composition during engine update.");
if (std::abs(temp_comp.getMolarAbundance(species) - y_data[i]) > 1e-12) {
throw exceptions::UtilityError("Conversion from solver state to composition molar abundance failed verification.");
}
}
#endif
NetIn netInTemp = netIn;
NetIn netInTemp;
netInTemp.temperature = T9 * 1e9; // Convert back to Kelvin
netInTemp.density = netIn.density;
netInTemp.composition = temp_comp;
LOG_DEBUG(
m_logger,
"Prior to Engine update composition is (molar abundance) {}",
[&]() -> std::string {
std::stringstream ss;
const size_t ns = temp_comp.size();
size_t count = 0;
for (const auto &symbol : temp_comp | std::views::keys) {
ss << symbol << ": " << temp_comp.getMolarAbundance(symbol);
if (count < ns - 1) {
ss << ", ";
}
count++;
}
return ss.str();
}()
);
LOG_DEBUG(
m_logger,
"Prior to Engine Update active reactions are: {}",
[&]() -> std::string {
std::stringstream ss;
const gridfire::reaction::ReactionSet& reactions = m_engine.getNetworkReactions();
size_t count = 0;
for (const auto& reaction : reactions) {
ss << reaction -> id();
if (count < reactions.size() - 1) {
ss << ", ";
}
count++;
}
return ss.str();
}()
);
fourdst::composition::Composition currentComposition = m_engine.update(netInTemp);
LOG_DEBUG(
m_logger,
"After to Engine update composition is (molar abundance) {}",
[&]() -> std::string {
std::stringstream ss;
const size_t ns = currentComposition.size();
size_t count = 0;
for (const auto &symbol : currentComposition | std::views::keys) {
ss << symbol << ": " << currentComposition.getMolarAbundance(symbol);
if (count < ns - 1) {
ss << ", ";
}
count++;
}
return ss.str();
}()
);
LOG_DEBUG(
m_logger,
"Fractional Abundance Changes: {}",
[&]() -> std::string {
std::stringstream ss;
const size_t ns = currentComposition.size();
size_t count = 0;
for (const auto &symbol : currentComposition | std::views::keys) {
if (!temp_comp.contains(symbol)) {
ss << symbol << ": New Species";
} else {
const double old_X = temp_comp.getMolarAbundance(symbol);
const double new_X = currentComposition.getMolarAbundance(symbol);
const double frac_change = (new_X - old_X) / (old_X + std::numeric_limits<double>::epsilon());
ss << symbol << ": " << frac_change;
}
if (count < ns - 1) {
ss << ", ";
}
count++;
}
return ss.str();
}()
);
LOG_DEBUG(
m_logger,
"After Engine Update active reactions are: {}",
[&]() -> std::string {
std::stringstream ss;
const gridfire::reaction::ReactionSet& reactions = m_engine.getNetworkReactions();
size_t count = 0;
for (const auto& reaction : reactions) {
ss << reaction -> id();
if (count < reactions.size() - 1) {
ss << ", ";
}
count++;
}
return ss.str();
}()
);
LOG_INFO(
m_logger,
"Due to a triggered stale engine the composition was updated from size {} to {} species.",
"Due to a triggered engine update the composition was updated from size {} to {} species.",
num_species_at_stop,
m_engine.getNetworkSpecies().size()
);
@@ -350,40 +448,21 @@ namespace gridfire::solver {
sunrealtype* y_data = N_VGetArrayPointer(m_Y);
accumulated_energy += y_data[numSpecies];
std::vector<double> y_vec(y_data, y_data + numSpecies);
std::vector<double> finalMassFractions(numSpecies);
for (size_t i = 0; i < numSpecies; ++i) {
const double molarMass = m_engine.getNetworkSpecies()[i].mass();
finalMassFractions[i] = y_data[i] * molarMass; // Convert from molar abundance to mass fraction
if (finalMassFractions[i] < MIN_ABUNDANCE_THRESHOLD) {
finalMassFractions[i] = 0.0;
for (size_t i = 0; i < y_vec.size(); ++i) {
if (y_vec[i] < 0 && std::abs(y_vec[i]) < 1e-16) {
y_vec[i] = 0.0; // Regularize tiny negative abundances to zero
}
}
std::vector<std::string> speciesNames;
speciesNames.reserve(numSpecies);
for (const auto& species : m_engine.getNetworkSpecies()) {
speciesNames.emplace_back(species.name());
}
LOG_TRACE_L2(m_logger, "Constructing final composition= with {} species", numSpecies);
LOG_TRACE_L2(m_logger, "Constructing final composition= with {} species", speciesNames.size());
fourdst::composition::Composition topLevelComposition(speciesNames);
topLevelComposition.setMassFraction(speciesNames, finalMassFractions);
bool didFinalizeTopLevel = topLevelComposition.finalize(true);
if (!didFinalizeTopLevel) {
LOG_ERROR(m_logger, "Failed to finalize top level reconstructed composition after CVODE integration. Check output mass fractions for validity.");
throw std::runtime_error("Failed to finalize output composition after CVODE integration.");
}
fourdst::composition::Composition topLevelComposition(m_engine.getNetworkSpecies(), y_vec);
fourdst::composition::Composition outputComposition = m_engine.collectComposition(topLevelComposition);
assert(outputComposition.getRegisteredSymbols().size() == equilibratedComposition.getRegisteredSymbols().size());
bool didFinalizeOutput = outputComposition.finalize(false);
if (!didFinalizeOutput) {
LOG_ERROR(m_logger, "Failed to finalize output composition after CVODE integration.");
throw std::runtime_error("Failed to finalize output composition after CVODE integration.");
}
LOG_TRACE_L2(m_logger, "Final composition constructed successfully!");
LOG_TRACE_L2(m_logger, "Constructing output data...");
@@ -511,19 +590,7 @@ namespace gridfire::solver {
// some consistent memory layout for the composition object to make this cheeper. That optimization would need to be
// done in the libcomposition library though...
std::vector<double> y_vec(y_data, y_data + numSpecies);
std::vector<std::string> symbols;
symbols.reserve(numSpecies);
for (const auto& species : m_engine.getNetworkSpecies()) {
symbols.emplace_back(species.name());
}
std::vector<double> M;
M.reserve(numSpecies);
for (size_t i = 0; i < numSpecies; ++i) {
const double molarMass = m_engine.getNetworkSpecies()[i].mass();
M.push_back(molarMass);
}
std::vector<double> X = utils::massFractionFromMolarAbundanceAndMolarMass(y_vec, M);
fourdst::composition::Composition composition(symbols, X);
fourdst::composition::Composition composition(m_engine.getNetworkSpecies(), y_vec);
const auto result = m_engine.calculateRHSAndEnergy(composition, data->T9, data->rho);
if (!result) {
@@ -623,14 +690,9 @@ namespace gridfire::solver {
check_cvode_flag(CVodeGetEstLocalErrors(m_cvode_mem, m_YErr), "CVodeGetEstLocalErrors");
sunrealtype *y_data = N_VGetArrayPointer(m_Y);
sunrealtype *y_err_data = N_VGetArrayPointer(m_YErr);
std::vector<double> err_ratios;
std::vector<std::string> speciesNames;
const auto absTol = m_config.get<double>("gridfire:solver:CVODESolverStrategy:absTol", 1.0e-8);
const auto relTol = m_config.get<double>("gridfire:solver:CVODESolverStrategy:relTol", 1.0e-8);
const size_t num_components = N_VGetLength(m_Y);
err_ratios.resize(num_components - 1);
@@ -646,26 +708,13 @@ namespace gridfire::solver {
0.0
);
std::vector<double> M;
M.reserve(num_components);
for (size_t i = 0; i < num_components - 1; i++) {
const double weight = relTol * std::abs(y_data[i]) + absTol;
if (weight == 0.0) continue; // Skip components with zero weight
const double err_ratio = std::abs(y_err_data[i]) / weight;
err_ratios[i] = err_ratio;
speciesNames.emplace_back(user_data.networkSpecies->at(i).name());
M.push_back(user_data.networkSpecies->at(i).mass());
}
std::vector<double> X = utils::massFractionFromMolarAbundanceAndMolarMass(Y_full, M);
fourdst::composition::Composition composition(speciesNames, X);
fourdst::composition::Composition composition(user_data.engine->getNetworkSpecies(), Y_full);
if (err_ratios.empty()) {
return;
}
std::vector<size_t> indices(speciesNames.size());
std::vector<size_t> indices(composition.size());
for (size_t i = 0; i < indices.size(); ++i) {
indices[i] = i;
}
@@ -677,29 +726,29 @@ namespace gridfire::solver {
}
);
std::vector<std::string> sorted_speciesNames;
std::vector<fourdst::atomic::Species> sorted_species;
std::vector<double> sorted_err_ratios;
sorted_speciesNames.reserve(indices.size());
sorted_species.reserve(indices.size());
sorted_err_ratios.reserve(indices.size());
for (const auto idx: indices) {
sorted_speciesNames.push_back(speciesNames[idx]);
sorted_species.push_back(composition.getSpeciesAtIndex(idx));
sorted_err_ratios.push_back(err_ratios[idx]);
}
std::vector<std::unique_ptr<utils::ColumnBase>> columns;
columns.push_back(std::make_unique<utils::Column<std::string>>("Species", sorted_speciesNames));
columns.push_back(std::make_unique<utils::Column<fourdst::atomic::Species>>("Species", sorted_species));
columns.push_back(std::make_unique<utils::Column<double>>("Error Ratio", sorted_err_ratios));
std::cout << utils::format_table("Species Error Ratios", columns) << std::endl;
if (displayJacobianStiffness) {
diagnostics::inspect_jacobian_stiffness(*user_data.engine, composition, user_data.T9, user_data.rho);
for (const auto& species : sorted_speciesNames) {
diagnostics::inspect_species_balance(*user_data.engine, species, composition, user_data.T9, user_data.rho);
for (const auto& species : sorted_species) {
diagnostics::inspect_species_balance(*user_data.engine, std::string(species.name()), composition, user_data.T9, user_data.rho);
}
}

10
src/lib/solver/strategies/triggers/engine_partitioning_trigger.cpp
View File

@@ -344,11 +344,11 @@ namespace gridfire::trigger::solver::CVODE {
}
float ConvergenceFailureTrigger::current_mean() const {
float acc = 0;
size_t acc = 0;
for (const auto nlcfails: m_window) {
acc += nlcfails;
}
return acc / m_windowSize;
return static_cast<float>(acc) / static_cast<float>(m_windowSize);
}
bool ConvergenceFailureTrigger::abs_failure(
@@ -364,7 +364,7 @@ namespace gridfire::trigger::solver::CVODE {
const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx
) const {
const float mean = current_mean();
if (ctx.currentConvergenceFailures - mean > m_relativeFailureRate * mean) {
if (static_cast<float>(ctx.currentConvergenceFailures) - mean > m_relativeFailureRate * mean) {
return true;
}
return false;
@@ -395,8 +395,8 @@ namespace gridfire::trigger::solver::CVODE {
// Combine the triggers using logical OR
auto orTriggerA = std::make_unique<OrTrigger<ctx_t>>(std::move(simulationTimeTrigger), std::move(offDiagTrigger));
auto orTriggerB = std::make_unique<OrTrigger<ctx_t>>(std::move(orTriggerA), std::move(timestepGrowthTrigger));
auto orTriggerC = std::make_unique<OrTrigger<ctx_t>>(std::move(orTriggerB), std::move(convergenceFailureTrigger));
// auto orTriggerC = std::make_unique<OrTrigger<ctx_t>>(std::move(orTriggerB), std::move(convergenceFailureTrigger));
return orTriggerC;
return convergenceFailureTrigger;
}
}

1
src/lib/utils/logging.cpp
View File

@@ -7,7 +7,6 @@
#include <ranges>
#include <string_view>
#include <string>
#include <vector>
std::string gridfire::utils::formatNuclearTimescaleLogString(
const DynamicEngine& engine,

1
src/meson.build
View File

@@ -1,7 +1,6 @@
# Define the library
gridfire_sources = files(
'lib/network.cpp',
'lib/engine/engine_approx8.cpp',
'lib/engine/engine_graph.cpp',
'lib/engine/views/engine_adaptive.cpp',
'lib/engine/views/engine_defined.cpp',

12
src/python/engine/trampoline/py_engine.h
View File

@@ -14,7 +14,7 @@ public:
const std::vector<fourdst::atomic::Species>& getNetworkSpecies() const override;
std::expected<gridfire::StepDerivatives<double>,gridfire::expectations::StaleEngineError> calculateRHSAndEnergy(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -27,7 +27,7 @@ public:
const std::vector<fourdst::atomic::Species>& getNetworkSpecies() const override;
std::expected<gridfire::StepDerivatives<double>, gridfire::expectations::StaleEngineError> calculateRHSAndEnergy(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -39,7 +39,7 @@ public:
) const override;
void generateJacobianMatrix(
const fourdst::composition::Composition &comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho,
const std::vector<fourdst::atomic::Species> &activeSpecies
@@ -66,7 +66,7 @@ public:
double calculateMolarReactionFlow(
const gridfire::reaction::Reaction &reaction,
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -78,7 +78,7 @@ public:
) override;
std::expected<std::unordered_map<fourdst::atomic::Species, double>, gridfire::expectations::StaleEngineError> getSpeciesTimescales(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
double T9,
double rho
) const override;
@@ -119,7 +119,7 @@ public:
throw std::logic_error("Network depth not supported by this engine.");
}
void rebuild(
const fourdst::composition::Composition& comp,
const fourdst::composition::CompositionAbstract &comp,
gridfire::BuildDepthType depth
) override {
throw std::logic_error("Setting network depth not supported by this engine.");

2
src/python/reaction/bindings.cpp
View File

@@ -333,12 +333,14 @@ void register_reaction_bindings(py::module &m) {
.def(
"__getitem__",
py::overload_cast<size_t>(&gridfire::reaction::ReactionSet::operator[], py::const_),
py::return_value_policy::reference,
py::arg("index"),
"Get a LogicalReaclibReaction by index."
)
.def(
"__getitem___",
py::overload_cast<const std::string_view&>(&gridfire::reaction::ReactionSet::operator[], py::const_),
py::return_value_policy::reference,
py::arg("id"),
"Get a LogicalReaclibReaction by its ID."
)

2
subprojects/fourdst.wrap
View File

@@ -1,4 +1,4 @@
[wrap-git]
url = https://github.com/4D-STAR/fourdst
revision = v0.8.5
revision = v0.9.4
depth = 1

218
tests/graphnet_sandbox/main.cpp
View File

@@ -1,13 +1,11 @@
#include <iostream>
#include <fstream>
#include <format>
#include "gridfire/engine/engine_graph.h"
#include "gridfire/engine/engine_approx8.h"
#include "gridfire/engine/views/engine_adaptive.h"
#include "gridfire/partition/partition_types.h"
#include "gridfire/engine/views/engine_multiscale.h"
#include "gridfire/solver/strategies/CVODE_solver_strategy.h"
#include "gridfire/policy/stellar_policy.h"
#include "gridfire/utils/table_format.h"
#include "gridfire/network.h"
@@ -19,25 +17,16 @@
#include "quill/Logger.h"
#include "quill/LogMacros.h"
#include "quill/Backend.h"
#include "quill/Frontend.h"
#include <chrono>
#include <clocale>
#include <functional>
#include "gridfire/engine/views/engine_defined.h"
#include "gridfire/partition/composite/partition_composite.h"
#include "fourdst/atomic/species.h"
#include "fourdst/composition/utils.h"
static std::terminate_handler g_previousHandler = nullptr;
void measure_execution_time(const std::function<void()>& callback, const std::string& name)
{
const auto startTime = std::chrono::steady_clock::now();
callback();
const auto endTime = std::chrono::steady_clock::now();
const auto duration = std::chrono::duration_cast<std::chrono::nanoseconds>(endTime - startTime);
std::cout << "Execution time for " << name << ": " << duration.count()/1e9 << " s\n";
}
void quill_terminate_handler()
{
quill::Backend::stop();
@@ -47,62 +36,18 @@ void quill_terminate_handler()
std::abort();
}
int main(int argc, char* argv[]){
// Valid usages are either
// ./graphnet_sandbox
//or
// ./graphnet_sandbox --plug <plugin_bundle_path>
if (argc == 3 && std::string(argv[1]) == "--plug") {
std::filesystem::path pluginBundlePath(argv[2]);
if (!std::filesystem::exists(pluginBundlePath)) {
std::cerr << "Error: Plugin bundle path does not exist: " << pluginBundlePath << "\n";
std::cerr << "Usage: " << argv[0] << " [--plug <plugin_bundle_path>]\n";
return 1;
}
std::cout << "Loading plugin bundle from: " << pluginBundlePath << "\n";
fourdst::plugin::bundle::PluginBundle pluginBundle(pluginBundlePath);
}
if (argc == 2 && std::string(argv[1]) != "--plug") {
std::cerr << "Invalid argument: " << argv[1] << "\n";
std::cerr << "Usage: " << argv[0] << " [--plug <plugin_bundle_path>]\n";
return 1;
}
if (argc == 2 && std::string(argv[1]) == "--plug") {
std::cerr << "Error: No plugin bundle path provided.\n";
std::cerr << "Usage: " << argv[0] << " [--plug <plugin_bundle_path>]\n";
return 1;
}
if (argc > 3) {
std::cerr << "Too many arguments provided.\n";
std::cerr << "Usage: " << argv[0] << " [--plug <plugin_bundle_path>]\n";
return 1;
}
gridfire::NetIn init() {
std::setlocale(LC_ALL, "");
g_previousHandler = std::set_terminate(quill_terminate_handler);
quill::Logger* logger = fourdst::logging::LogManager::getInstance().getLogger("log");
logger->set_log_level(quill::LogLevel::TraceL2);
LOG_INFO(logger, "Starting Adaptive Engine View Example...");
logger->set_log_level(quill::LogLevel::TraceL1);
using namespace gridfire;
const std::vector<double> comp = {0.708, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4};
const std::vector<double> X = {0.7081145999999999, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4};
const std::vector<std::string> symbols = {"H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"};
fourdst::composition::Composition composition;
composition.registerSymbol(symbols, true);
composition.setMassFraction(symbols, comp);
bool didFinalize = composition.finalize(true);
if (!didFinalize) {
std::cerr << "Failed to finalize initial composition." << std::endl;
return 1;
}
using partition::BasePartitionType;
const auto partitionFunction = partition::CompositePartitionFunction({
BasePartitionType::RauscherThielemann,
BasePartitionType::GroundState
});
fourdst::composition::Composition composition = fourdst::composition::buildCompositionFromMassFractions(symbols, X);
NetIn netIn;
netIn.composition = composition;
@@ -110,37 +55,134 @@ int main(int argc, char* argv[]){
netIn.density = 1.6e2;
netIn.energy = 0;
// TODO: There is a bug when I get to very low concentrations of hydrogen where the solver will crash. I suspect this can be resolved with triggers
netIn.tMax = 3e17;
// netIn.tMax = 1e-14;
netIn.tMax = 3e16;
netIn.dt0 = 1e-12;
GraphEngine ReaclibEngine(composition, partitionFunction, NetworkBuildDepth::SecondOrder, NetworkConstructionFlags::DEFAULT);
return netIn;
}
ReaclibEngine.setPrecomputation(true);
ReaclibEngine.setUseReverseReactions(false);
void log_results(const gridfire::NetOut& netOut, const gridfire::NetIn& netIn) {
double initialH1_X = netIn.composition.getMassFraction("H-1");
double finalH1_X = netOut.composition.getMassFraction("H-1");
double deltaH1_X = finalH1_X - initialH1_X;
double fractionalChangeH1 = (deltaH1_X) / initialH1_X * 100.0;
DefinedEngineView subsetView({"p(p,e+)d", "d(p,g)he3", "n(p,g)d", "d(n,g)t"}, ReaclibEngine);
MultiscalePartitioningEngineView partitioningView(ReaclibEngine);
double initialHe4_X = netIn.composition.getMassFraction("He-4");
double finalHe4_X = netOut.composition.getMassFraction("He-4");
double deltaHe4_X = finalHe4_X - initialHe4_X;
double fractionalChangeHe4 = (deltaHe4_X) / initialHe4_X * 100.0;
AdaptiveEngineView adaptiveView(partitioningView);
double initialBe7_X = 0.0;
double finalBe7_X = netOut.composition.getMassFraction("Be-7");
double deltaBe7_X = finalBe7_X - initialBe7_X;
double fractionalChangeBe7 = (deltaBe7_X) / initialBe7_X * 100.0;
double initialC12_X = netIn.composition.getMassFraction("C-12");
double finalC12_X = netOut.composition.getMassFraction("C-12");
double deltaC12_X = finalC12_X - initialC12_X;
double fractionalChangeC12 = (deltaC12_X) / initialC12_X * 100.0;
double finalEnergy = netOut.energy;
double dEpsdT = netOut.dEps_dT;
double dEpsdRho = netOut.dEps_dRho;
double initialMeanMolecularMass = netIn.composition.getMeanParticleMass();
double finalMeanMolecularMass = netOut.composition.getMeanParticleMass();
std::vector<std::string> rowLabels = {"H-1", "He-4", "Be-7", "C-12", "ε", "dε/dT", "dε/dρ", "<μ>"};
std::vector<std::string> header = {"Parameter", "Initial", "Final", "Δ", "% Change"};
std::vector<double> H1Data = {initialH1_X, finalH1_X, deltaH1_X, fractionalChangeH1};
std::vector<double> He4Data = {initialHe4_X, finalHe4_X, deltaHe4_X, fractionalChangeHe4};
std::vector<double> Be7Data = {initialBe7_X, finalBe7_X, deltaBe7_X, fractionalChangeBe7};
std::vector<double> C12Data = {initialC12_X, finalC12_X, deltaC12_X, fractionalChangeC12};
std::vector<double> energyData = {0.0, finalEnergy, finalEnergy, 0.0};
std::vector<double> dEpsdTData = {0.0, dEpsdT, dEpsdT, 0.0};
std::vector<double> dEpsdRhoData = {0.0, dEpsdRho, dEpsdRho, 0.0};
std::vector<double> meanMolecularMassData = {initialMeanMolecularMass, finalMeanMolecularMass, finalMeanMolecularMass - initialMeanMolecularMass,
(finalMeanMolecularMass - initialMeanMolecularMass) / initialMeanMolecularMass * 100.0};
gridfire::utils::Column<std::string> paramCol(header[0], rowLabels);
gridfire::utils::Column<double> initialCol(header[1], {H1Data[0], He4Data[0], Be7Data[0], C12Data[0], energyData[0], dEpsdTData[0], dEpsdRhoData[0], meanMolecularMassData[0]});
gridfire::utils::Column<double> finalCol (header[2], {H1Data[1], He4Data[1], Be7Data[1], C12Data[1], energyData[1], dEpsdTData[1], dEpsdRhoData[1], meanMolecularMassData[1]});
gridfire::utils::Column<double> deltaCol (header[3], {H1Data[2], He4Data[2], Be7Data[2], C12Data[2], energyData[2], dEpsdTData[2], dEpsdRhoData[2], meanMolecularMassData[2]});
gridfire::utils::Column<double> percentCol(header[4], {H1Data[3], He4Data[3], Be7Data[3], C12Data[3], energyData[3], dEpsdTData[3], dEpsdRhoData[3], meanMolecularMassData[3]});
std::vector<std::unique_ptr<gridfire::utils::ColumnBase>> columns;
columns.push_back(std::make_unique<gridfire::utils::Column<std::string>>(paramCol));
columns.push_back(std::make_unique<gridfire::utils::Column<double>>(initialCol));
columns.push_back(std::make_unique<gridfire::utils::Column<double>>(finalCol));
columns.push_back(std::make_unique<gridfire::utils::Column<double>>(deltaCol));
columns.push_back(std::make_unique<gridfire::utils::Column<double>>(percentCol));
fourdst::composition::Composition updatedComposition = adaptiveView.update(netIn);
gridfire::utils::print_table("Simulation Results", columns);
}
adaptiveView.generateJacobianMatrix(updatedComposition, netIn.temperature/1e9, netIn.density);
solver::CVODESolverStrategy solver(adaptiveView);
NetOut netOut;
netOut = solver.evaluate(netIn, true);
static std::vector<std::pair<double, std::unordered_map<std::string, std::pair<double, double>>>> g_callbackHistory;
void record_abundance_history_callback(const gridfire::solver::CVODESolverStrategy::TimestepContext& ctx) {
const auto& engine = ctx.engine;
std::unordered_map<std::string, std::pair<double, double>> abundances;
for (const auto& species : engine.getNetworkSpecies()) {
const size_t sid = engine.getSpeciesIndex(species);
abundances[std::string(species.name())] = std::make_pair(species.mass(), N_VGetArrayPointer(ctx.state)[sid]);
}
g_callbackHistory.emplace_back(ctx.t, abundances);
}
std::cout << "Initial H-1: " << netIn.composition.getMassFraction("H-1") << std::endl;
std::cout << "NetOut H-1: " << netOut.composition.getMassFraction("H-1") << std::endl;
double initialHydrogen = netIn.composition.getMassFraction("H-1");
double finalHydrogen = netOut.composition.getMassFraction("H-1");
double fractionalConsumedHydrogen = (initialHydrogen - finalHydrogen) / initialHydrogen * 100.0;
std::cout << "Fractional consumed hydrogen: " << fractionalConsumedHydrogen << "%" << std::endl;
std::cout << "Final D abundance " << netOut.composition.getMolarAbundance("H-2") << std::endl;
std::cout << netOut << std::endl;
void save_callback_data(const std::string_view filename) {
std::set<std::string> unique_species;
for (const auto &abundances: g_callbackHistory | std::views::values) {
for (const auto &species_name: abundances | std::views::keys) {
unique_species.insert(species_name);
}
}
std::ofstream csvFile(filename.data(), std::ios::out);
csvFile << "t,";
size_t i = 0;
for (const auto& species_name : unique_species) {
csvFile << species_name;
if (i < unique_species.size() - 1) {
csvFile << ",";
}
i++;
}
csvFile << "\n";
for (const auto& [time, data] : g_callbackHistory) {
csvFile << time << ",";
size_t j = 0;
for (const auto& species_name : unique_species) {
if (!data.contains(species_name)) {
csvFile << "0.0";
} else {
csvFile << data.at(species_name).second;
}
if (j < unique_species.size() - 1) {
csvFile << ",";
}
++j;
}
csvFile << "\n";
}
csvFile.close();
}
int main() {
using namespace gridfire;
const NetIn netIn = init();
policy::MainSequencePolicy stellarPolicy(netIn.composition);
DynamicEngine& engine = stellarPolicy.construct();
solver::CVODESolverStrategy solver(engine);
solver.set_callback(solver::CVODESolverStrategy::TimestepCallback(record_abundance_history_callback));
const NetOut netOut = solver.evaluate(netIn, false);
log_results(netOut, netIn);
save_callback_data("abundance_history.csv");
}

0
tests/graphnet_sandbox/post.dat Normal file
View File

0
tests/graphnet_sandbox/prior.dat Normal file
View File

61
utils/reaclib/format.py
View File

@@ -69,9 +69,11 @@ from collections import defaultdict
from re import Match
from typing import List, Tuple, Any, LiteralString
import numpy as np
from serif.atomic import species
from serif.atomic import Species
from serif.constants import Constants
from fourdst.atomic import species
from fourdst.atomic import Species
from fourdst.constants import Constants
import hashlib
from collections import Counter
import math
@@ -281,13 +283,13 @@ def translate_names_to_species(names: List[str]) -> List[Species]:
try:
sp.append(species[name])
except Exception as e:
print("Error: Species not found in database:", name, e)
raise ReaclibParseError(f"Species '{name}' not found in species database.", line_content=name)
return sp
def determine_reaction_type(reactants: List[str],
products: List[str],
qValue: float
qValue: float,
chapter: int
) -> Tuple[str, List[str], List[str], str, str, str]:
"""
Analyze a reaction for quantum number conservation and classify projectiles/ejectiles.
@@ -323,8 +325,6 @@ def determine_reaction_type(reactants: List[str],
- Identifies weak (leptonic) and photonic processes.
- Determines projectiles/ejectiles based on mass and reaction type.
"""
if abs(qValue - 4.621) < 1e-6:
print("Looking at he3(he3, 2p)he4")
# --- helper look-ups ----------------------------------------------------
reactantSpecies = translate_names_to_species(reactants)
productSpecies = translate_names_to_species(products)
@@ -360,24 +360,41 @@ def determine_reaction_type(reactants: List[str],
projectiles: List[str] = []
ejectiles: List[str] = []
debug = False
if reactants == ['b8'] and products == ['be8']:
debug = True
BETA_PLUS_THRESHOLD_ENERGY = 1.022 # MeV
# -----------------------------------------------------------------------
# 1. Charged-lepton bookkeeping (|ΔZ| = 1) ------------------------------
# -----------------------------------------------------------------------
if abs(dZ) == 1:
# Proton → neutron (β⁻ / e- capture)
if dZ == -1:
# Electron capture when (i) exo-thermic and (ii) nucleus count unchanged
if qValue > 0 and dN == 0:
projectiles.append("e-") # write e- as projectile
else:
ejectiles.append("e-") # β⁻ decay: e- is an ejectile
if debug:
print("============")
print("Reactant Species: ", reactantSpecies)
print("Product Species: ", productSpecies)
print("Target Species: ", targetSpecies)
print("Residual Species: ", residualSpecies)
print("Nuclear Projectiles: ", nuclearProjectiles)
print("Nuclear Ejectiles: ", nuclearEjectiles)
print("aReact, aProd: ", aReact, aProd)
print("zReact, zProd: ", zReact, zProd)
print("nReact, nProd: ", nReact, nProd)
print("dA, dZ, dN: ", dA, dZ, dN)
print("qValue: ", qValue)
# Neutron → proton (β⁺ / positron capture – capture is vanishingly rare)
elif dZ == 1:
ejectiles.append("e+") # β⁺ / weak-proton capture
if dZ == -1 and chapter == 1:
ejectiles.append("e-") # β- decay
if dZ == -1 and chapter != 1:
projectiles.append("e+") # positron capture
if dZ == 1 and chapter == 1 and qValue >= BETA_PLUS_THRESHOLD_ENERGY:
ejectiles.append("e+") # β+ Decay
if dZ == 1 and chapter == 1 and qValue < BETA_PLUS_THRESHOLD_ENERGY:
projectiles.append("e-") # electron capture
if dZ == 1 and chapter != 1:
ejectiles.append("e+") # Positron as byproduct of two body reaction
# Neutrino companion is implicit – never written
# (dL is automatically fixed by hiding ν or ν̄)
# -----------------------------------------------------------------------
# 2. Photon bookkeeping (ΔZ = 0) ----------------------------------------
@@ -452,12 +469,10 @@ def extract_groups(match: re.Match, reverse: bool) -> Reaction:
chapter = int(groups[0].strip())
rawGroup = groups[1].strip()
rList, pList = get_rp(rawGroup, chapter)
if 'c12' in rList and 'mg24' in pList:
print("Found it!")
if reverse:
rList, pList = pList, rList
qValue = float(groups[3].strip())
target, proj, ejec, residual, key, rType = determine_reaction_type(rList, pList, qValue)
target, proj, ejec, residual, key, rType = determine_reaction_type(rList, pList, qValue, chapter)
reaction = Reaction(
reactants=rList,
products=pList,