diff --git a/src/include/gridfire/engine/types/graph.h b/src/include/gridfire/engine/types/graph.h
deleted file mode 100644
index e69de29b..00000000
diff --git a/src/lib/engine/types/graph.cpp b/src/lib/engine/types/graph.cpp
deleted file mode 100644
index e69de29b..00000000
diff --git a/tests/graphnet_sandbox/WeakReactionsGraph.svg b/tests/graphnet_sandbox/WeakReactionsGraph.svg
deleted file mode 100644
index b70e338d..00000000
--- a/tests/graphnet_sandbox/WeakReactionsGraph.svg
+++ /dev/null
@@ -1,17411 +0,0 @@
-
-
-
-
-
diff --git a/tests/graphnet_sandbox/main.cpp b/tests/graphnet_sandbox/main.cpp
index f4dd0d94..ef2ca0a3 100644
--- a/tests/graphnet_sandbox/main.cpp
+++ b/tests/graphnet_sandbox/main.cpp
@@ -23,89 +23,122 @@
#include "fourdst/atomic/species.h"
#include "fourdst/composition/utils.h"
+#include "gridfire/exceptions/error_solver.h"
+
+#include
+
+#include "gridfire/reaction/reaclib.h"
+#include "gridfire/solver/strategies/CVODE_solver_strategy.h"
static std::terminate_handler g_previousHandler = nullptr;
+boost::json::object g_reaction_contribution_history;
+static std::vector>>> g_callbackHistory;
+static bool s_wrote_abundance_history = false;
+static bool s_wrote_reaction_history = false;
+void quill_terminate_handler();
-void quill_terminate_handler()
-{
- quill::Backend::stop();
- if (g_previousHandler)
- g_previousHandler();
- else
- std::abort();
+
+inline std::unique_ptr build_partition_function() {
+ using gridfire::partition::BasePartitionType;
+ const auto partitionFunction = gridfire::partition::CompositePartitionFunction({
+ BasePartitionType::RauscherThielemann,
+ BasePartitionType::GroundState
+ });
+ return std::make_unique(partitionFunction);
}
-gridfire::NetIn init() {
+gridfire::NetIn init(const double temp) {
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::TraceL1);
+ logger->set_log_level(quill::LogLevel::TraceL2);
using namespace gridfire;
const std::vector X = {0.7081145999999999, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4};
const std::vector symbols = {"H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"};
- fourdst::composition::Composition composition = fourdst::composition::buildCompositionFromMassFractions(symbols, X);
+ const fourdst::composition::Composition composition = fourdst::composition::buildCompositionFromMassFractions(symbols, X);
NetIn netIn;
netIn.composition = composition;
- netIn.temperature = 1.5e7;
- netIn.density = 1.6e2;
+ netIn.temperature = temp;
+ netIn.density = 1.5e2;
netIn.energy = 0;
- netIn.tMax = 3e16;
+ netIn.tMax = 3e17;
netIn.dt0 = 1e-12;
return netIn;
}
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;
+ std::vector logSpecies = {
+ 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
+ };
- 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;
+ std::vector initial;
+ std::vector final;
+ std::vector delta;
+ std::vector fractional;
+ for (const auto& species : logSpecies) {
+ double initial_X = netIn.composition.getMassFraction(species);
+ double final_X = netOut.composition.getMassFraction(species);
+ double delta_X = final_X - initial_X;
+ double fractionalChange = (delta_X) / initial_X * 100.0;
- 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;
+ initial.push_back(initial_X);
+ final.push_back(final_X);
+ delta.push_back(delta_X);
+ fractional.push_back(fractionalChange);
+ }
- 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;
+ initial.push_back(0.0); // Placeholder for energy
+ final.push_back(netOut.energy);
+ delta.push_back(netOut.energy);
+ fractional.push_back(0.0); // Placeholder for energy
- double finalEnergy = netOut.energy;
- double dEpsdT = netOut.dEps_dT;
- double dEpsdRho = netOut.dEps_dRho;
+ initial.push_back(0.0);
+ final.push_back(netOut.dEps_dT);
+ delta.push_back(netOut.dEps_dT);
+ fractional.push_back(0.0);
- double initialMeanMolecularMass = netIn.composition.getMeanParticleMass();
- double finalMeanMolecularMass = netOut.composition.getMeanParticleMass();
+ initial.push_back(0.0);
+ final.push_back(netOut.dEps_dRho);
+ delta.push_back(netOut.dEps_dRho);
+ fractional.push_back(0.0);
- std::vector rowLabels = {"H-1", "He-4", "Be-7", "C-12", "ε", "dε/dT", "dε/dρ", "<μ>"};
- std::vector header = {"Parameter", "Initial", "Final", "Δ", "% Change"};
- std::vector H1Data = {initialH1_X, finalH1_X, deltaH1_X, fractionalChangeH1};
- std::vector He4Data = {initialHe4_X, finalHe4_X, deltaHe4_X, fractionalChangeHe4};
- std::vector Be7Data = {initialBe7_X, finalBe7_X, deltaBe7_X, fractionalChangeBe7};
- std::vector C12Data = {initialC12_X, finalC12_X, deltaC12_X, fractionalChangeC12};
- std::vector energyData = {0.0, finalEnergy, finalEnergy, 0.0};
- std::vector dEpsdTData = {0.0, dEpsdT, dEpsdT, 0.0};
- std::vector dEpsdRhoData = {0.0, dEpsdRho, dEpsdRho, 0.0};
- std::vector meanMolecularMassData = {initialMeanMolecularMass, finalMeanMolecularMass, finalMeanMolecularMass - initialMeanMolecularMass,
- (finalMeanMolecularMass - initialMeanMolecularMass) / initialMeanMolecularMass * 100.0};
+ initial.push_back(netIn.composition.getMeanParticleMass());
+ final.push_back(netOut.composition.getMeanParticleMass());
+ delta.push_back(final.back() - initial.back());
+ fractional.push_back((final.back() - initial.back()) / initial.back() * 100.0);
- gridfire::utils::Column paramCol(header[0], rowLabels);
- gridfire::utils::Column initialCol(header[1], {H1Data[0], He4Data[0], Be7Data[0], C12Data[0], energyData[0], dEpsdTData[0], dEpsdRhoData[0], meanMolecularMassData[0]});
- gridfire::utils::Column finalCol (header[2], {H1Data[1], He4Data[1], Be7Data[1], C12Data[1], energyData[1], dEpsdTData[1], dEpsdRhoData[1], meanMolecularMassData[1]});
- gridfire::utils::Column deltaCol (header[3], {H1Data[2], He4Data[2], Be7Data[2], C12Data[2], energyData[2], dEpsdTData[2], dEpsdRhoData[2], meanMolecularMassData[2]});
- gridfire::utils::Column percentCol(header[4], {H1Data[3], He4Data[3], Be7Data[3], C12Data[3], energyData[3], dEpsdTData[3], dEpsdRhoData[3], meanMolecularMassData[3]});
+ std::vector rowLabels = [&]() -> std::vector {
+ std::vector labels;
+ for (const auto& species : logSpecies) {
+ labels.push_back(std::string(species.name()));
+ }
+ labels.push_back("ε");
+ labels.push_back("dε/dT");
+ labels.push_back("dε/dρ");
+ labels.push_back("<μ>");
+ return labels;
+ }();
+
+
+ gridfire::utils::Column paramCol("Parameter", rowLabels);
+ gridfire::utils::Column initialCol("Initial", initial);
+ gridfire::utils::Column finalCol ("Final", final);
+ gridfire::utils::Column deltaCol ("δ", delta);
+ gridfire::utils::Column percentCol("% Change", fractional);
std::vector> columns;
columns.push_back(std::make_unique>(paramCol));
@@ -119,17 +152,59 @@ void log_results(const gridfire::NetOut& netOut, const gridfire::NetIn& netIn) {
}
-static std::vector>>> g_callbackHistory;
void record_abundance_history_callback(const gridfire::solver::CVODESolverStrategy::TimestepContext& ctx) {
+ s_wrote_abundance_history = true;
const auto& engine = ctx.engine;
- std::unordered_map> abundances;
+ // std::unordered_map> abundances;
+ std::vector Y;
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]);
+ double y = N_VGetArrayPointer(ctx.state)[sid];
+ Y.push_back(y > 0.0 ? y : 0.0); // Regularize tiny negative abundances to zero
+ }
+
+ fourdst::composition::Composition comp(engine.getNetworkSpecies(), Y);
+
+
+ std::unordered_map> abundances;
+ for (const auto& sp : comp | std::views::keys) {
+ abundances.emplace(std::string(sp.name()), std::make_pair(sp.mass(), comp.getMolarAbundance(sp)));
}
g_callbackHistory.emplace_back(ctx.t, abundances);
}
+size_t g_iters = 0;
+void record_contribution_callback(const gridfire::solver::CVODESolverStrategy::TimestepContext& ctx) {
+ s_wrote_reaction_history = true;
+ boost::json::object timestep;
+ boost::json::object reaction_contribution;
+ boost::json::object species_abundance;
+ std::set activeSpecies(ctx.engine.getNetworkSpecies().begin(), ctx.engine.getNetworkSpecies().end());
+ for (const auto& [species, contributions] : ctx.reactionContributionMap) {
+ boost::json::object species_obj;
+ for (const auto& [reaction_id, contribution] : contributions) {
+ species_obj[reaction_id] = contribution;
+ }
+ reaction_contribution[std::string(species.name())] = species_obj;
+
+ double y;
+ if (activeSpecies.contains(species)) {
+ const size_t sid = ctx.engine.getSpeciesIndex(species);
+ y = N_VGetArrayPointer(ctx.state)[sid];
+ } else {
+ y = 0.0;
+ }
+
+ species_abundance[std::string(species.name())] = y;
+ }
+ timestep["t"] = ctx.t;
+ timestep["dt"] = ctx.dt;
+ timestep["reaction_contribution"] = reaction_contribution;
+ timestep["species_abundance"] = species_abundance;
+ g_reaction_contribution_history[std::to_string(g_iters)] = timestep;
+ g_iters++;
+}
+
void save_callback_data(const std::string_view filename) {
std::set unique_species;
for (const auto &abundances: g_callbackHistory | std::views::values) {
@@ -171,18 +246,46 @@ void save_callback_data(const std::string_view filename) {
csvFile.close();
}
+void log_callback_data() {
+ if (s_wrote_abundance_history) {
+ std::cout << "Saving abundance history to abundance_history.csv" << std::endl;
+ save_callback_data("abundance_history.csv");
+ }
+
+ if (s_wrote_reaction_history) {
+ std::cout << "Saving reaction history to reaction_contribution_history.json" << std::endl;
+ std::ofstream jsonFile("reaction_contribution_history.json", std::ios::out);
+ jsonFile << boost::json::serialize(g_reaction_contribution_history);
+ jsonFile.close();
+ }
+}
+
+void quill_terminate_handler()
+{
+ log_callback_data();
+ quill::Backend::stop();
+ if (g_previousHandler)
+ g_previousHandler();
+ else
+ std::abort();
+}
+
+void callback_main(const gridfire::solver::CVODESolverStrategy::TimestepContext& ctx) {
+ record_abundance_history_callback(ctx);
+ record_contribution_callback(ctx);
+}
+
int main() {
using namespace gridfire;
- const NetIn netIn = init();
+ const NetIn netIn = init(1.5e7);
policy::MainSequencePolicy stellarPolicy(netIn.composition);
+ stellarPolicy.construct();
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);
-
+ solver.set_callback(solver::CVODESolverStrategy::TimestepCallback(callback_main));
+ const NetOut netOut = solver.evaluate(netIn);
log_results(netOut, netIn);
- save_callback_data("abundance_history.csv");
+ log_callback_data();
}
\ No newline at end of file
diff --git a/tests/graphnet_sandbox/post.dat b/tests/graphnet_sandbox/post.dat
deleted file mode 100644
index e69de29b..00000000
diff --git a/tests/graphnet_sandbox/prior.dat b/tests/graphnet_sandbox/prior.dat
deleted file mode 100644
index e69de29b..00000000