From 5a1a904e71f440ff8c2aee4589e851aec12028df Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 13 Apr 2026 07:19:18 -0400 Subject: [PATCH] refactor(GridFire): updated outputs --- .../SingleZoneSolver/gf_wall_vs_temp.cpp | 134 ++++++++ benchmarks/SingleZoneSolver/main.cpp | 232 +++++++++----- benchmarks/SingleZoneSolver/meson.build | 6 + subprojects/fourdst.wrap | 2 +- tests/graphnet_sandbox/main.cpp | 30 +- tests/python/py_test_single.py | 12 +- tools/gf_bbn/main.cpp | 298 ++++++++++++++++++ tools/gf_bbn/meson.build | 1 + tools/gf_quick/main.cpp | 277 ++++++---------- tools/meson.build | 3 +- 10 files changed, 704 insertions(+), 291 deletions(-) create mode 100644 benchmarks/SingleZoneSolver/gf_wall_vs_temp.cpp create mode 100644 tools/gf_bbn/main.cpp create mode 100644 tools/gf_bbn/meson.build diff --git a/benchmarks/SingleZoneSolver/gf_wall_vs_temp.cpp b/benchmarks/SingleZoneSolver/gf_wall_vs_temp.cpp new file mode 100644 index 00000000..d014687e --- /dev/null +++ b/benchmarks/SingleZoneSolver/gf_wall_vs_temp.cpp @@ -0,0 +1,134 @@ +// ReSharper disable CppUnusedIncludeDirective +#include +#include +#include +#include +#include + +#include "gridfire/gridfire.h" +#include // Required for parallel_setup + +#include "fourdst/composition/composition.h" +#include "fourdst/logging/logging.h" +#include "fourdst/atomic/species.h" +#include "fourdst/composition/utils.h" + +#include "quill/Logger.h" +#include "quill/Backend.h" + +#include + +#include "gridfire/reaction/reaclib.h" +#include "gridfire/utils/gf_omp.h" + +template +[[nodiscard]] constexpr auto linspace(T start, T end, std::size_t num_points) -> std::vector { + if (num_points == 0) { + return {}; + } + + if (num_points == 1) { + return {start}; + } + + return std::views::iota(0uz, num_points) + | std::views::transform([=](std::size_t i) -> T { + const T t = static_cast(i) / static_cast(num_points - 1); + + return std::lerp(start, end, t); + }) + | std::ranges::to>(); +} + +gridfire::NetIn init(const double temp, const double rho, const double tMax) { + std::setlocale(LC_ALL, ""); + quill::Logger* logger = fourdst::logging::LogManager::getInstance().getLogger("log"); + 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"}; + + + const fourdst::composition::Composition composition = fourdst::composition::buildCompositionFromMassFractions(symbols, X); + + NetIn netIn; + netIn.composition = composition; + netIn.temperature = temp; + netIn.density = rho; + netIn.energy = 0; + + netIn.tMax = tMax; + netIn.dt0 = 1e-12; + + return netIn; +} + + +int main() { + GF_PAR_INIT() + using namespace gridfire; + + constexpr double temp_init = 1.5e7; + constexpr double rho_init = 1.5e2; + constexpr double tMax = 3.1536e+12; + + NetIn netIn = init(temp_init, rho_init, tMax); + + policy::MainSequencePolicy stellarPolicy(netIn.composition); + const policy::ConstructionResults construct = stellarPolicy.construct(); + std::println("Sandbox Engine Stack: {}", stellarPolicy); + std::println("Scratch Blob State: {}", *construct.scratch_blob); + + + // arrays to store timings + + // Total number of interpolated data points + constexpr size_t N = 20; + std::array eval_times{}; + auto density = linspace(10.0, 5.0e2, N); + auto temperature = linspace(4e6,3e7, N); + + solver::PointSolverContext solverCtx(*construct.scratch_blob); + solverCtx.set_stdout_logging(false); + solver::PointSolver solver(construct.engine); + + auto startTime = std::chrono::high_resolution_clock::now(); + + + size_t i = 0; + for (const auto temp : temperature) { + for (const auto dens : density) { + std::println("Evaluation {:3}/{:5} ({:3.0f}%): ρ = {:10.4E}, T = {:10.4E}", i + 1, N*N, 100.0*((static_cast(i)+1.0)/(N*N)), dens, temp); + netIn.temperature = temp; + netIn.density = dens; + try { + auto start_eval_time = std::chrono::high_resolution_clock::now(); + const NetOut netOut = solver.evaluate(solverCtx, netIn); + auto end_eval_time = std::chrono::high_resolution_clock::now(); + std::chrono::duration eval_elapsed = end_eval_time - start_eval_time; + eval_times[i] = eval_elapsed.count(); + } catch (const gridfire::exceptions::GridFireError& e) { + std::cerr << "Error during evaluation " << (i + 1) << ": " << e.what() << std::endl; + eval_times[i] = std::numeric_limits::quiet_NaN(); + } + i++; + } + } + auto endTime = std::chrono::high_resolution_clock::now(); + + std::println("Total time for {} evaluations: {} seconds", N, (endTime - startTime).count()); + for (size_t j = 0; j < static_cast(N*N); ++j) { + std::println("Evaluation {}: {} seconds", j + 1, eval_times[j]); + } + + std::ofstream outfile("gf_wall_vs_temp_results.csv"); + outfile << "Evaluation,Density,Temperature,TimeSeconds\n"; + size_t j = 0; + for (const auto temp: temperature) { + for (const auto dens: density ) { + outfile << (j + 1) << "," << dens << ","<< temp << "," << eval_times[j] << "\n"; + j++; + } + } +} \ No newline at end of file diff --git a/benchmarks/SingleZoneSolver/main.cpp b/benchmarks/SingleZoneSolver/main.cpp index 72edb819..d450a66b 100644 --- a/benchmarks/SingleZoneSolver/main.cpp +++ b/benchmarks/SingleZoneSolver/main.cpp @@ -15,11 +15,15 @@ #include "quill/Logger.h" #include "quill/Backend.h" +#include "nlohmann/json.hpp" #include +#include + #include "gridfire/reaction/reaclib.h" #include "gridfire/utils/gf_omp.h" +#include "gridfire/utils/config.h" gridfire::NetIn init(const double temp, const double rho, const double tMax) { std::setlocale(LC_ALL, ""); @@ -63,97 +67,157 @@ int main() { std::println("Scratch Blob State: {}", *construct.scratch_blob); - constexpr size_t runs = 10; - auto startTime = std::chrono::high_resolution_clock::now(); + constexpr size_t runs = 100; + nlohmann::json results; + nlohmann::json metadata; - // arrays to store timings - std::array, runs> setup_times; - std::array, runs> eval_times; - std::array serial_results; - for (size_t i = 0; i < runs; ++i) { - auto start_setup_time = std::chrono::high_resolution_clock::now(); - solver::PointSolverContext solverCtx(*construct.scratch_blob); - solverCtx.set_stdout_logging(false); - solver::PointSolver solver(construct.engine); - auto end_setup_time = std::chrono::high_resolution_clock::now(); - std::chrono::duration setup_elapsed = end_setup_time - start_setup_time; - setup_times[i] = setup_elapsed; + const auto now = std::chrono::system_clock::now(); + std::string now_str = std::format("{:%Y-%m-%d %H:%M:%S}", now); - auto start_eval_time = std::chrono::high_resolution_clock::now(); - const NetOut netOut = solver.evaluate(solverCtx, netIn); - auto end_eval_time = std::chrono::high_resolution_clock::now(); - serial_results[i] = netOut; - std::chrono::duration eval_elapsed = end_eval_time - start_eval_time; - eval_times[i] = eval_elapsed; - } - auto endTime = std::chrono::high_resolution_clock::now(); - std::chrono::duration elapsed = endTime - startTime; - std::println(""); + metadata["Datetime"] = now_str; + metadata["GF_Version"] = version::toString(); - // Summarize serial timings - double total_setup_time = 0.0; - double total_eval_time = 0.0; - for (size_t i = 0; i < runs; ++i) { - total_setup_time += setup_times[i].count(); - total_eval_time += eval_times[i].count(); - } - std::println("Average Setup Time over {} runs: {:.6f} seconds", runs, total_setup_time / runs); - std::println("Average Evaluation Time over {} runs: {:.6f} seconds", runs, total_eval_time / runs); - std::println("Total Time for {} runs: {:.6f} seconds", runs, elapsed.count()); - - - std::array parallelResults; - std::array, runs> setupTimes; - std::array, runs> evalTimes; - std::array, runs> workspaces; - for (size_t i = 0; i < runs; ++i) { - workspaces[i] = construct.scratch_blob->clone_structure(); + utsname buffer{}; + if (uname(&buffer) == 0) { + std::string osName = buffer.sysname; +#ifdef __APPLE__ + if (osName == "Darwin") osName = "macOS"; +#endif + metadata["OS"] = osName; + metadata["OS Version"] = buffer.release; + metadata["Architecture"] = buffer.machine; + } else { + metadata["OS"] = "Unknown"; } +#if defined(__clang__) + metadata["Compiler"] = "Clang " __clang_version__; +#elif defined(__GNUC__) + metadata["Compiler"] = "GCC " __VERSION__; +#else + metadata["Compiler"] = "Unknown"; +#endif - // Parallel runs - startTime = std::chrono::high_resolution_clock::now(); + metadata["Threads"] = omp_get_max_threads(); + metadata["Runs"] = runs; + metadata["Temperature"] = temp; + metadata["Density"] = rho; + metadata["tMax_per_run_s"] = tMax; - GF_OMP(parallel for, for (size_t i = 0; i < runs; ++i)) { - auto start_setup_time = std::chrono::high_resolution_clock::now(); - solver::PointSolverContext solverCtx(*construct.scratch_blob); - solverCtx.set_stdout_logging(false); - solver::PointSolver solver(construct.engine); - auto end_setup_time = std::chrono::high_resolution_clock::now(); - std::chrono::duration setup_elapsed = end_setup_time - start_setup_time; - setupTimes[i] = setup_elapsed; - auto start_eval_time = std::chrono::high_resolution_clock::now(); - parallelResults[i] = solver.evaluate(solverCtx, netIn); - auto end_eval_time = std::chrono::high_resolution_clock::now(); - std::chrono::duration eval_elapsed = end_eval_time - start_eval_time; - evalTimes[i] = eval_elapsed; - } - endTime = std::chrono::high_resolution_clock::now(); - elapsed = endTime - startTime; - std::println(""); + results["Metadata"] = metadata; - // Summarize parallel timings - total_setup_time = 0.0; - total_eval_time = 0.0; - for (size_t i = 0; i < runs; ++i) { - total_setup_time += setupTimes[i].count(); - total_eval_time += evalTimes[i].count(); + for (size_t rID = 0; rID < runs; rID++) { + nlohmann::json run_result; + nlohmann::json run_metadata; + run_metadata["num_zones"] = rID; + run_result["metadata"] = run_metadata; + + + auto startTime = std::chrono::high_resolution_clock::now(); + + // arrays to store timings + std::array, runs> setup_times{}; + std::array, runs> eval_times{}; + std::array serial_results; + for (size_t i = 0; i < rID; ++i) { + auto start_setup_time = std::chrono::high_resolution_clock::now(); + solver::PointSolverContext solverCtx(*construct.scratch_blob); + solverCtx.set_stdout_logging(false); + solver::PointSolver solver(construct.engine); + auto end_setup_time = std::chrono::high_resolution_clock::now(); + std::chrono::duration setup_elapsed = end_setup_time - start_setup_time; + setup_times[i] = setup_elapsed; + + auto start_eval_time = std::chrono::high_resolution_clock::now(); + const NetOut netOut = solver.evaluate(solverCtx, netIn); + auto end_eval_time = std::chrono::high_resolution_clock::now(); + serial_results[i] = netOut; + std::chrono::duration eval_elapsed = end_eval_time - start_eval_time; + eval_times[i] = eval_elapsed; + } + auto endTime = std::chrono::high_resolution_clock::now(); + std::chrono::duration elapsed = endTime - startTime; + std::println(""); + + nlohmann::json point_solver_time_results; + point_solver_time_results["total_time_s"] = elapsed.count(); + run_result["Serial"] = point_solver_time_results; + + // Summarize serial timings + double total_setup_time = 0.0; + double total_eval_time = 0.0; + for (size_t i = 0; i < rID; ++i) { + total_setup_time += setup_times[i].count(); + total_eval_time += eval_times[i].count(); + } + std::println("Average Setup Time over {} runs: {:.6f} seconds", runs, total_setup_time / runs); + std::println("Average Evaluation Time over {} runs: {:.6f} seconds", runs, total_eval_time / runs); + std::println("Total Time for {} runs: {:.6f} seconds", runs, elapsed.count()); + + + std::array parallelResults; + std::array, runs> setupTimes; + std::array, runs> evalTimes; + std::array, runs> workspaces; + for (size_t i = 0; i < rID; ++i) { + workspaces[i] = construct.scratch_blob->clone_structure(); + } + + + // Parallel runs + startTime = std::chrono::high_resolution_clock::now(); + + GF_OMP(parallel for, for (size_t i = 0; i < rID; ++i)) { + auto start_setup_time = std::chrono::high_resolution_clock::now(); + solver::PointSolverContext solverCtx(*construct.scratch_blob); + solverCtx.set_stdout_logging(false); + solver::PointSolver solver(construct.engine); + auto end_setup_time = std::chrono::high_resolution_clock::now(); + std::chrono::duration setup_elapsed = end_setup_time - start_setup_time; + setupTimes[i] = setup_elapsed; + auto start_eval_time = std::chrono::high_resolution_clock::now(); + parallelResults[i] = solver.evaluate(solverCtx, netIn); + auto end_eval_time = std::chrono::high_resolution_clock::now(); + std::chrono::duration eval_elapsed = end_eval_time - start_eval_time; + evalTimes[i] = eval_elapsed; + } + endTime = std::chrono::high_resolution_clock::now(); + elapsed = endTime - startTime; + std::println(""); + + nlohmann::json grid_solver_results; + grid_solver_results["total_time_s"] = elapsed.count(); + run_result["Parallel"] = grid_solver_results; + + // Summarize parallel timings + total_setup_time = 0.0; + total_eval_time = 0.0; + for (size_t i = 0; i < runs; ++i) { + total_setup_time += setupTimes[i].count(); + total_eval_time += evalTimes[i].count(); + } + + std::println("Average Parallel Setup Time over {} runs: {:.6f} seconds", runs, total_setup_time / runs); + std::println("Average Parallel Evaluation Time over {} runs: {:.6f} seconds", runs, total_eval_time / runs); + std::println("Total Parallel Time for {} runs: {:.6f} seconds", runs, elapsed.count()); + + std::println("========== Summary =========="); + std::println("Serial Runs:"); + std::println(" Average Setup Time: {:.6f} seconds", total_setup_time / runs); + std::println(" Average Evaluation Time: {:.6f} seconds", total_eval_time / runs); + std::println("Parallel Runs:"); + std::println(" Average Setup Time: {:.6f} seconds", total_setup_time / runs); + std::println(" Average Evaluation Time: {:.6f} seconds", total_eval_time / runs); + std::println("Difference:"); + std::println(" Setup Time Difference: {:.6f} seconds", (total_setup_time / runs) - (total_setup_time / runs)); + std::println(" Evaluation Time Difference: {:.6f} seconds", (total_eval_time / runs) - (total_eval_time / runs)); + std::println(" Setup Time Fractional Difference: {:.2f}%", ((total_setup_time / runs) - (total_setup_time / runs)) / (total_setup_time / runs) * 100.0); + std::println(" Evaluation Time Fractional Difference: {:.2f}%", ((total_eval_time / runs) - (total_eval_time / runs)) / (total_eval_time / runs) * 100.0); + + results[std::format("Run_{}", rID)] = run_result; } - std::println("Average Parallel Setup Time over {} runs: {:.6f} seconds", runs, total_setup_time / runs); - std::println("Average Parallel Evaluation Time over {} runs: {:.6f} seconds", runs, total_eval_time / runs); - std::println("Total Parallel Time for {} runs: {:.6f} seconds", runs, elapsed.count()); - - std::println("========== Summary =========="); - std::println("Serial Runs:"); - std::println(" Average Setup Time: {:.6f} seconds", total_setup_time / runs); - std::println(" Average Evaluation Time: {:.6f} seconds", total_eval_time / runs); - std::println("Parallel Runs:"); - std::println(" Average Setup Time: {:.6f} seconds", total_setup_time / runs); - std::println(" Average Evaluation Time: {:.6f} seconds", total_eval_time / runs); - std::println("Difference:"); - std::println(" Setup Time Difference: {:.6f} seconds", (total_setup_time / runs) - (total_setup_time / runs)); - std::println(" Evaluation Time Difference: {:.6f} seconds", (total_eval_time / runs) - (total_eval_time / runs)); - std::println(" Setup Time Fractional Difference: {:.2f}%", ((total_setup_time / runs) - (total_setup_time / runs)) / (total_setup_time / runs) * 100.0); - std::println(" Evaluation Time Fractional Difference: {:.2f}%", ((total_eval_time / runs) - (total_eval_time / runs)) / (total_eval_time / runs) * 100.0); -} \ No newline at end of file + std::ofstream o("gf_single_zone_solver_benchmark_results.json"); + o << std::setw(4) << results << std::endl; + o.close(); +} diff --git a/benchmarks/SingleZoneSolver/meson.build b/benchmarks/SingleZoneSolver/meson.build index 68c2aa8e..6d110b53 100644 --- a/benchmarks/SingleZoneSolver/meson.build +++ b/benchmarks/SingleZoneSolver/meson.build @@ -3,3 +3,9 @@ executable( 'main.cpp', dependencies: [gridfire_dep], ) + +executable( + 'gf_wall_vs_temp', + 'gf_wall_vs_temp.cpp', + dependencies: [gridfire_dep] +) diff --git a/subprojects/fourdst.wrap b/subprojects/fourdst.wrap index a857546c..f1f5f5b8 100644 --- a/subprojects/fourdst.wrap +++ b/subprojects/fourdst.wrap @@ -1,4 +1,4 @@ [wrap-git] url = https://github.com/4D-STAR/fourdst -revision = v0.9.19 +revision = v0.9.21 depth = 1 diff --git a/tests/graphnet_sandbox/main.cpp b/tests/graphnet_sandbox/main.cpp index 795b212a..3de8db86 100644 --- a/tests/graphnet_sandbox/main.cpp +++ b/tests/graphnet_sandbox/main.cpp @@ -5,10 +5,8 @@ #include #include -#include "gridfire/gridfire.h" #include // Required for parallel_setup -#include "fourdst/composition/composition.h" #include "fourdst/logging/logging.h" #include "fourdst/atomic/species.h" #include "fourdst/composition/utils.h" @@ -19,6 +17,8 @@ #include +#include "gridfire/gridfire.h" +#include "fourdst/composition/composition.h" #include "gridfire/utils/gf_omp.h" @@ -229,31 +229,21 @@ void callback_main(const gridfire::solver::PointSolverTimestepContext& ctx) { } int main() { - GF_PAR_INIT(); using namespace gridfire; - constexpr size_t breaks = 1; - double temp = 1.5e7; - double rho = 1.5e2; - double tMax = 3.1536e+16/breaks; + constexpr double temp = 1.5e7; + constexpr double rho = 1.6e2; + constexpr double tMax = 3e17; const NetIn netIn = init(temp, rho, tMax); - policy::MainSequencePolicy stellarPolicy(netIn.composition); - auto [engine, ctx_template] = stellarPolicy.construct(); - std::println("Sandbox Engine Stack: {}", stellarPolicy); - std::println("Scratch Blob State: {}", *ctx_template); - constexpr size_t nZones = 100; - std::array netIns; - for (size_t zone = 0; zone < nZones; ++zone) { - netIns[zone] = netIn; - netIns[zone].temperature = 1.5e7; - } + auto engine = engine::GraphEngine(netIn.composition, engine::NetworkBuildDepth::Full); + auto blob = engine.constructStateBlob(); const solver::PointSolver localSolver(engine); - solver::GridSolverContext solverCtx(*ctx_template); - const solver::GridSolver gridSolver(engine, localSolver); + solver::PointSolverContext solverCtx(*blob); - std::vector netOuts = gridSolver.evaluate(solverCtx, netIns | std::ranges::to()); + auto result = localSolver.evaluate(solverCtx, netIn, false, false); + std::cout << result << std::endl; } diff --git a/tests/python/py_test_single.py b/tests/python/py_test_single.py index 7e5f4276..e7274ad1 100644 --- a/tests/python/py_test_single.py +++ b/tests/python/py_test_single.py @@ -6,6 +6,8 @@ from fourdst.composition import CanonicalComposition from fourdst.atomic import Species from gridfire.type import NetIn +from logger import StepLogger + def rescale_composition(comp_ref : Composition, ZZs : float, Y_primordial : float = 0.248) -> Composition: CC : CanonicalComposition = comp_ref.getCanonicalComposition() @@ -61,13 +63,17 @@ def years_to_seconds(years: float) -> float: def main(): C = init_composition() - netIn = init_netIn(2.75e6, 1.5e1, years_to_seconds(10e9), C) + netIn = init_netIn(1.5e7, 1.6e2, years_to_seconds(10e9), C) policy = MainSequencePolicy(C) construct = policy.construct() solver = PointSolver(construct.engine) solver_ctx = PointSolverContext(construct.scratch_blob) - results = solver.evaluate(solver_ctx, netIn, False, False) - print(results) + + stepLogger = StepLogger() + solver_ctx.callback = lambda ctx: stepLogger.log_step(ctx); + solver.evaluate(solver_ctx, netIn, False, False) + stepLogger.to_json("test_single.json", TestName="test_single"); + if __name__ == "__main__": main() diff --git a/tools/gf_bbn/main.cpp b/tools/gf_bbn/main.cpp new file mode 100644 index 00000000..d2e2669a --- /dev/null +++ b/tools/gf_bbn/main.cpp @@ -0,0 +1,298 @@ +// ReSharper disable CppUnusedIncludeDirective +#include +#include +#include +#include +#include + +#include "gridfire/gridfire.h" + +#include "fourdst/composition/composition.h" +#include "fourdst/logging/logging.h" +#include "fourdst/atomic/species.h" +#include "fourdst/composition/utils.h" + +#include "quill/Logger.h" +#include "quill/Backend.h" +#include "CLI/CLI.hpp" + +#include +#include + +#include "gridfire/utils/gf_omp.h" + +#include "nlohmann/json.hpp" + +struct IntermediateResult { + double time; + fourdst::composition::Composition comp; + double current_energy; + double current_neutrino_loss_rate; +}; + +static std::vector g_callbackHistory; + +gridfire::NetIn init(const double temp, const double rho, const double tMax) { + std::setlocale(LC_ALL, ""); + quill::Logger* logger = fourdst::logging::LogManager::getInstance().getLogger("log"); + logger->set_log_level(quill::LogLevel::Info); + + using namespace gridfire; + constexpr double XpXn = 7.17; + constexpr double Xn = 1.0 / (1.0 + XpXn); + constexpr double Xp = 1 - Xn; + + const std::vector X = {Xp, Xn}; + const std::vector symbols = {"H-1", "n-1"}; + + + const fourdst::composition::Composition composition = fourdst::composition::buildCompositionFromMassFractions(symbols, X); + + NetIn netIn; + netIn.composition = composition; + netIn.temperature = temp; + netIn.density = rho; + netIn.energy = 0; + + netIn.tMax = tMax; + netIn.dt0 = 1e-12; + + return netIn; +} + +void log_results(const gridfire::NetOut& netOut, const gridfire::NetIn& netIn) { + 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 + }; + + 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; + + initial.push_back(initial_X); + final.push_back(final_X); + delta.push_back(delta_X); + fractional.push_back(fractionalChange); + } + + 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 + + initial.push_back(0.0); + final.push_back(netOut.dEps_dT); + delta.push_back(netOut.dEps_dT); + fractional.push_back(0.0); + + initial.push_back(0.0); + final.push_back(netOut.dEps_dRho); + delta.push_back(netOut.dEps_dRho); + fractional.push_back(0.0); + + initial.push_back(0.0); + final.push_back(netOut.specific_neutrino_energy_loss); + delta.push_back(netOut.specific_neutrino_energy_loss); + fractional.push_back(0.0); + + initial.push_back(0.0); + final.push_back(netOut.specific_neutrino_flux); + delta.push_back(netOut.specific_neutrino_flux); + fractional.push_back(0.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); + + std::vector rowLabels = [&]() -> std::vector { + std::vector labels; + for (const auto& species : logSpecies) { + labels.emplace_back(species.name()); + } + labels.emplace_back("ε"); + labels.emplace_back("dε/dT"); + labels.emplace_back("dε/dρ"); + labels.emplace_back("Eν"); + labels.emplace_back("Fν"); + labels.emplace_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)); + columns.push_back(std::make_unique>(initialCol)); + columns.push_back(std::make_unique>(finalCol)); + columns.push_back(std::make_unique>(deltaCol)); + columns.push_back(std::make_unique>(percentCol)); + + + gridfire::utils::print_table("Simulation Results", columns); +} + + +void record_abundance_history_callback(const gridfire::solver::PointSolverTimestepContext& ctx) { + const auto& engine = ctx.engine; + std::vector Y; + for (const auto& species : engine.getNetworkSpecies(ctx.state_ctx)) { + const size_t sid = engine.getSpeciesIndex(ctx.state_ctx, species); + double y = N_VGetArrayPointer(ctx.state)[sid]; + Y.push_back(y > 0.0 ? y : 0.0); // Regularize tiny negative abundances to zero + } + + const fourdst::composition::Composition comp(engine.getNetworkSpecies(ctx.state_ctx), Y); + IntermediateResult stepResult; + stepResult.comp = comp; + stepResult.time = ctx.t; + stepResult.current_energy = ctx.current_total_energy; + stepResult.current_neutrino_loss_rate = ctx.current_neutrino_energy_loss_rate; + g_callbackHistory.push_back(stepResult); +} + + + +void callback_main(const gridfire::solver::PointSolverTimestepContext& ctx) { + record_abundance_history_callback(ctx); +} + +void save_callback(const std::string& filename) { + // Save to JSON + nlohmann::json j; + for (const auto& record : g_callbackHistory) { + nlohmann::json entry; + entry["time"] = record.time; + entry["current_energy"] = record.current_energy; + entry["current_neutrino_loss_rate"] = record.current_neutrino_loss_rate; + // make a sub-json for composition + nlohmann::json comp_json; + for (const auto& [species, abundance] : record.comp) { + comp_json[species.name()] = abundance; + } + entry["composition"] = comp_json; + j.push_back(entry); + } + std::ofstream ofs(filename); + ofs << j.dump(4); +} + +double T9(const double age) { + return 10.0/std::sqrt(age); +} + +double density(const double age) { + return 4e-5 * std::pow(T9(age), 3); +} + +int main(int argc, char** argv) { + GF_PAR_INIT(); + using namespace gridfire; + + double tMax = 3600; + double h = 0.1; + + CLI::App app("GridFire Quick BBN Test"); + app.add_option("--tmax", tMax, "Maximum Time in seconds")->default_val(std::format("{:5.2E}", tMax)); + app.add_option("--s", h, "Geometric Timestep Scale Factor")->default_val(std::format("{:5.2f}", h)); + + CLI11_PARSE(app, argc, argv); + + NetIn netIn = init(0, 0, tMax); + const engine::GraphEngine engine(netIn.composition); + auto blob = std::make_unique(); + blob->enroll(); + auto* graph_engine_state = engine::scratch::get_state(*blob); + graph_engine_state->initialize(engine); + + + solver::PointSolverContext solver_ctx(*blob); + solver::PointSolver solver(engine); + solver_ctx.stdout_logging=false; + + double current_time = 180; + nlohmann::json j; + nlohmann::json meta; + + meta["tMax"] = tMax; + meta["tStart"] = current_time; + meta["h"] = h; + j.push_back(meta); + + nlohmann::json steps; + + while (current_time < tMax) { + nlohmann::json entry; + double current_dt = h * current_time; + double next_time = current_time + current_dt; + + netIn.tMax = current_dt; + double current_temp = T9(current_time) * 1e9; + double next_temp = T9(next_time) * 1e9; + double burn_temp = (current_temp + next_temp)/2.0; + + double current_density = density(current_time); + double next_density = density(next_time); + double burn_density = (current_density + next_density)/2.0; + + netIn.temperature = burn_temp; + netIn.density = burn_density; + + + fourdst::composition::Composition initial_comp = netIn.composition; + + NetOut result = solver.evaluate(solver_ctx, netIn); + netIn.composition = result.composition; + std::println("Time: {:5.2E} (+{:5.2E}), Burn Temp: {:5.2E}, Burn Density: {:5.2E}", current_time, current_dt, burn_temp, burn_density); + + const fourdst::composition::Composition& comp = result.composition; + + auto Xi = [&](const std::string& symbol) -> double { + if (!initial_comp.contains(symbol)) { + return 0; + } + return initial_comp.getMassFraction(symbol); + }; + + entry["t"] = current_time; + entry["T"] = burn_temp; + entry["D"] = burn_density; + entry["mu"] = comp.getMeanParticleMass(); + auto lifetimes = engine.getSpeciesTimescales(*blob, comp, burn_temp, burn_density); + if (lifetimes.has_value()) { + entry["tau_n"] = lifetimes.value().at(fourdst::atomic::n_1); + } + + for (const auto& [sp, finalY] : comp) { + double initialY = Xi(std::string(sp.name())); + entry[std::format("{}_f", sp.name())] = finalY; + entry[std::format("{}_i", sp.name())] = initialY; + } + + steps.push_back(entry); + + + current_time += current_dt; + } + j.push_back(steps); + std::ofstream out("BBNResults.json"); + out << j.dump(4); +} diff --git a/tools/gf_bbn/meson.build b/tools/gf_bbn/meson.build new file mode 100644 index 00000000..1fae70d2 --- /dev/null +++ b/tools/gf_bbn/meson.build @@ -0,0 +1 @@ +executable('gf_bbn', 'main.cpp', dependencies: [gridfire_dep, cli11_dep]) diff --git a/tools/gf_quick/main.cpp b/tools/gf_quick/main.cpp index 0c1c7038..ae5d9b58 100644 --- a/tools/gf_quick/main.cpp +++ b/tools/gf_quick/main.cpp @@ -21,128 +21,26 @@ #include "gridfire/utils/gf_omp.h" +#include "nlohmann/json.hpp" -static std::terminate_handler g_previousHandler = nullptr; -static std::vector>>> g_callbackHistory; -static bool s_wrote_abundance_history = false; -void quill_terminate_handler(); +struct IntermediateResult { + double time{}; + fourdst::composition::Composition comp; + gridfire::reaction::ReactionSet reactions; + std::vector reaction_flows; + double current_energy{}; + double current_neutrino_loss_rate{}; -using namespace fourdst::composition; -Composition rescale(const Composition& comp, double target_X, double target_Z) { - // 1. Validate inputs - if (target_X < 0.0 || target_Z < 0.0 || (target_X + target_Z) > 1.0 + 1e-14) { - throw std::invalid_argument("Target mass fractions X and Z must be non-negative and sum to <= 1.0"); - } + gridfire::reaction::ReactionSet inactive_reactions; + std::vector inactive_reaction_flows; +}; - // Force high precision for the target Y to ensure X+Y+Z = 1.0 exactly in our logic - long double ld_target_X = static_cast(target_X); - long double ld_target_Z = static_cast(target_Z); - long double ld_target_Y = 1.0L - ld_target_X - ld_target_Z; - - // Clamp Y to 0 if it dipped slightly below due to precision (e.g. X+Z=1.0000000001) - if (ld_target_Y < 0.0L) ld_target_Y = 0.0L; - - // 2. Manually calculate current Mass Totals (bypass getCanonicalComposition to avoid crashes) - long double total_mass_H = 0.0L; - long double total_mass_He = 0.0L; - long double total_mass_Z = 0.0L; - - // We need to iterate and identify species types manually - // Standard definition: H (z=1), He (z=2), Metals (z>2) - // Note: We use long double accumulators to prevent summation drift - for (const auto& [spec, molar_abundance] : comp) { - // Retrieve atomic properties. - // Note: usage assumes fourdst::atomic::Species has .z() and .mass() - // consistent with the provided composition.cpp - int z = spec.z(); - double a = spec.mass(); - - long double mass_contribution = static_cast(molar_abundance) * static_cast(a); - - if (z == 1) { - total_mass_H += mass_contribution; - } else if (z == 2) { - total_mass_He += mass_contribution; - } else { - total_mass_Z += mass_contribution; - } - } - - long double total_mass_current = total_mass_H + total_mass_He + total_mass_Z; - - // Edge case: Empty composition - if (total_mass_current <= 0.0L) { - // Return empty or throw? If input was empty, return empty. - if (comp.size() == 0) return comp; - throw std::runtime_error("Input composition has zero total mass."); - } - - // 3. Calculate Scaling Factors - // Factor = (Target_Mass_Fraction / Old_Mass_Fraction) - // = (Target_Mass_Fraction) / (Old_Group_Mass / Total_Mass) - // = (Target_Mass_Fraction * Total_Mass) / Old_Group_Mass - - long double scale_H = 0.0L; - long double scale_He = 0.0L; - long double scale_Z = 0.0L; - - if (ld_target_X > 1e-16L) { - if (total_mass_H <= 1e-19L) { - throw std::runtime_error("Cannot rescale Hydrogen to " + std::to_string(target_X) + - " because input has no Hydrogen."); - } - scale_H = (ld_target_X * total_mass_current) / total_mass_H; - } - - if (ld_target_Y > 1e-16L) { - if (total_mass_He <= 1e-19L) { - throw std::runtime_error("Cannot rescale Helium to " + std::to_string((double)ld_target_Y) + - " because input has no Helium."); - } - scale_He = (ld_target_Y * total_mass_current) / total_mass_He; - } - - if (ld_target_Z > 1e-16L) { - if (total_mass_Z <= 1e-19L) { - throw std::runtime_error("Cannot rescale Metals to " + std::to_string(target_Z) + - " because input has no Metals."); - } - scale_Z = (ld_target_Z * total_mass_current) / total_mass_Z; - } - - // 4. Apply Scaling and Construct New Vectors - std::vector new_species; - std::vector new_abundances; - new_species.reserve(comp.size()); - new_abundances.reserve(comp.size()); - - for (const auto& [spec, abundance] : comp) { - new_species.push_back(spec); - - long double factor = 0.0L; - int z = spec.z(); - - if (z == 1) { - factor = scale_H; - } else if (z == 2) { - factor = scale_He; - } else { - factor = scale_Z; - } - - // Calculate new abundance in long double then cast back - long double new_val_ld = static_cast(abundance) * factor; - new_abundances.push_back(static_cast(new_val_ld)); - } - - return Composition(new_species, new_abundances); -} +static std::vector g_callbackHistory; gridfire::NetIn init(const double temp, const double rho, const double tMax) { 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::Info); + 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}; @@ -255,9 +153,7 @@ void log_results(const gridfire::NetOut& netOut, const gridfire::NetIn& netIn) { void record_abundance_history_callback(const gridfire::solver::PointSolverTimestepContext& ctx) { - s_wrote_abundance_history = true; const auto& engine = ctx.engine; - // std::unordered_map> abundances; std::vector Y; for (const auto& species : engine.getNetworkSpecies(ctx.state_ctx)) { const size_t sid = engine.getSpeciesIndex(ctx.state_ctx, species); @@ -265,80 +161,84 @@ void record_abundance_history_callback(const gridfire::solver::PointSolverTimest Y.push_back(y > 0.0 ? y : 0.0); // Regularize tiny negative abundances to zero } - fourdst::composition::Composition comp(engine.getNetworkSpecies(ctx.state_ctx), Y); + const fourdst::composition::Composition comp(engine.getNetworkSpecies(ctx.state_ctx), Y); + IntermediateResult stepResult; + stepResult.comp = comp; + stepResult.time = ctx.t; + stepResult.current_energy = ctx.current_total_energy; + stepResult.current_neutrino_loss_rate = ctx.current_neutrino_energy_loss_rate; + stepResult.reactions = engine.getNetworkReactions(ctx.state_ctx); - - 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))); + for (const auto& reactionPtr : stepResult.reactions) { + double flow = engine.calculateMolarReactionFlow(ctx.state_ctx, *reactionPtr, comp, ctx.T9, ctx.rho); + stepResult.reaction_flows.push_back(flow); } - g_callbackHistory.emplace_back(ctx.t, abundances); + + stepResult.inactive_reactions = engine.getInactiveNetworkReactions(ctx.state_ctx); + for (const auto& reactionPtr : stepResult.inactive_reactions) { + double flow = engine.getInactiveReactionMolarReactionFlow(ctx.state_ctx, *reactionPtr, comp, ctx.T9, ctx.rho); + stepResult.inactive_reaction_flows.push_back(flow); + } + g_callbackHistory.push_back(stepResult); } -void save_callback_data(const std::string_view filename) { - std::set 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(); -} - -void log_callback_data(const double temp) { - if (s_wrote_abundance_history) { - std::cout << "Saving abundance history to abundance_history.csv" << std::endl; - save_callback_data("abundance_history_" + std::to_string(temp) + ".csv"); - } - -} - -void quill_terminate_handler() -{ - log_callback_data(1.5e7); - quill::Backend::stop(); - if (g_previousHandler) - g_previousHandler(); - else - std::abort(); -} void callback_main(const gridfire::solver::PointSolverTimestepContext& ctx) { record_abundance_history_callback(ctx); } +void save_callback(const std::string& filename) { + // Save to JSON + nlohmann::json j; + for (const auto& record : g_callbackHistory) { + nlohmann::json entry; + entry["time"] = record.time; + entry["current_energy"] = record.current_energy; + entry["current_neutrino_loss_rate"] = record.current_neutrino_loss_rate; + // make a sub-json for composition + nlohmann::json comp_json; + for (const auto& [species, abundance] : record.comp) { + comp_json[species.name()] = abundance; + } + entry["composition"] = comp_json; + entry["reactions"] = nlohmann::json::array(); + for (const auto& [reaction, flow] : std::views::zip(record.reactions, record.reaction_flows)) { + nlohmann::json reaction_info; + reaction_info["id"] = reaction->id(); + reaction_info["flow"] = flow; + reaction_info["species"] = nlohmann::json::array(); + reaction_info["Q"] = reaction->qValue(); + for (const auto& sp : reaction->all_species()) { + nlohmann::json species_info; + species_info["name"] = sp.name(); + species_info["stoichiometry"] = reaction->stoichiometry(sp); + reaction_info["species"].push_back(species_info); + } + entry["reactions"].push_back(reaction_info); + } + + entry["inactive_reactions"] = nlohmann::json::array(); + for (const auto& [reaction, flow] : std::views::zip(record.inactive_reactions, record.inactive_reaction_flows)) { + nlohmann::json reaction_info; + reaction_info["id"] = reaction->id(); + reaction_info["flow"] = flow; + reaction_info["species"] = nlohmann::json::array(); + reaction_info["Q"] = reaction->qValue(); + for (const auto& sp : reaction->all_species()) { + nlohmann::json species_info; + species_info["name"] = sp.name(); + species_info["stoichiometry"] = reaction->stoichiometry(sp); + reaction_info["species"].push_back(species_info); + } + entry["inactive_reactions"].push_back(reaction_info); + } + j.push_back(entry); + } + std::ofstream ofs(filename); + ofs << j.dump(4); +} + int main(int argc, char** argv) { GF_PAR_INIT(); using namespace gridfire; @@ -346,12 +246,18 @@ int main(int argc, char** argv) { double temp = 1.5e7; double rho = 1.5e2; double tMax = 3.1536e+16; + bool save_intermediate_results = false; + bool display_trigger = false; + std::string output_filename = "abundance_history.json"; CLI::App app("GridFire Quick CLI Test"); app.add_option("--temp", temp, "Initial Temperature")->default_val(std::format("{:5.2E}", temp)); app.add_option("--rho", rho, "Initial Density")->default_val(std::format("{:5.2E}", rho)); app.add_option("--tmax", tMax, "Maximum Time")->default_val(std::format("{:5.2E}", tMax)); + app.add_option("--save_intermediate_results", save_intermediate_results, "Save Intermediate Results")->default_val("false"); + app.add_option("--output", output_filename, "Output filename for intermediate results")->default_val("abundance_history.json"); + app.add_option("--display_trigger_explanations", display_trigger, "Display trigger explanations during run")->default_val("false"); CLI11_PARSE(app, argc, argv); NetIn netIn = init(temp, rho, tMax); @@ -361,7 +267,14 @@ int main(int argc, char** argv) { solver::PointSolverContext solver_context(*ctx_template); solver::PointSolver solver(engine); + if (save_intermediate_results) { + solver_context.callback = solver::TimestepCallback(callback_main); + } - NetOut result = solver.evaluate(solver_context, netIn); + NetOut result = solver.evaluate(solver_context, netIn, display_trigger); log_results(result, netIn); + + if (save_intermediate_results) { + save_callback(output_filename); + } } diff --git a/tools/meson.build b/tools/meson.build index d4fae098..301ab23a 100644 --- a/tools/meson.build +++ b/tools/meson.build @@ -2,4 +2,5 @@ if get_option('build_tools') subdir('gf_config') subdir('gf_quick') subdir('gf_multi') -endif \ No newline at end of file + subdir('gf_bbn') +endif