From e3d7bed39c582c8d8e86360071998c7db5b9fbb1 Mon Sep 17 00:00:00 2001
From: Emily Boudreaux <thomas@boudreauxmail.com>
Date: Tue, 29 Jul 2025 07:49:17 -0400
Subject: [PATCH] test(reverse-rates): testing reverse rates

---
 tests/graphnet_sandbox/main.cpp | 28 ++++++++++++++--------------
 tests/python/test.py            | 21 ++++++++++++++-------
 2 files changed, 28 insertions(+), 21 deletions(-)

diff --git a/tests/graphnet_sandbox/main.cpp b/tests/graphnet_sandbox/main.cpp
index 1f93c2c6..4f7f54cd 100644
--- a/tests/graphnet_sandbox/main.cpp
+++ b/tests/graphnet_sandbox/main.cpp
@@ -71,26 +71,26 @@ int main() {
 
     NetIn netIn;
     netIn.composition = composition;
-    netIn.temperature = 1.5e7;
-    netIn.density = 1.6e2;
+    netIn.temperature = 5e9;
+    netIn.density = 1.6e6;
     netIn.energy = 0;
-    netIn.tMax = 3.1536e17; // ~ 10Gyr
+    // netIn.tMax = 3.1536e17; // ~ 10Gyr
+    netIn.tMax = 1e-14;
     netIn.dt0 = 1e-12;
-    for (const auto& [symbol, entry] : netIn.composition) {
-        std::cout << symbol << ": " << entry.mass_fraction() << "\n";
-    }
 
-    // GraphEngine ReaclibEngine(composition, partitionFunction, NetworkBuildDepth::SecondOrder);
-    //
-    // ReaclibEngine.setPrecomputation(true);
-    // ReaclibEngine.setUseReverseReactions(false);
+    GraphEngine ReaclibEngine(composition, partitionFunction, NetworkBuildDepth::SecondOrder);
+    ReaclibEngine.setUseReverseReactions(true);
     // ReaclibEngine.setScreeningModel(screening::ScreeningType::WEAK);
     //
-    // MultiscalePartitioningEngineView partitioningView(ReaclibEngine);
-    // AdaptiveEngineView adaptiveView(partitioningView);
+    MultiscalePartitioningEngineView partitioningView(ReaclibEngine);
+    AdaptiveEngineView adaptiveView(partitioningView);
     //
-    // solver::DirectNetworkSolver solver(adaptiveView);
-    // NetOut netOut;
+    solver::DirectNetworkSolver solver(adaptiveView);
+    NetOut netOut;
+    netOut = solver.evaluate(netIn);
+    std::cout << "Initial H-1: " << netIn.composition.getMassFraction("H-1") << std::endl;
+    std::cout << "NetOut H-1: " << netOut.composition.getMassFraction("H-1") << std::endl;
+    std::cout << "Consumed " << (netIn.composition.getMassFraction("H-1") - netOut.composition.getMassFraction("H-1")) * 100 << " % H-1 by mass" << std::endl;
     // measure_execution_time([&](){netOut = solver.evaluate(netIn);}, "DirectNetworkSolver Evaluation");
     // std::cout << "DirectNetworkSolver completed in " << netOut.num_steps << " steps.\n";
     // std::cout << "Final composition:\n";
diff --git a/tests/python/test.py b/tests/python/test.py
index 286fda0f..4923f097 100644
--- a/tests/python/test.py
+++ b/tests/python/test.py
@@ -7,28 +7,35 @@ from fourdst.composition import Composition
 symbols : list[str] = ["H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"]
 X : list[float]     = [0.708, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4]
 
-comp = Composition();
+
+comp = Composition()
 comp.registerSymbol(symbols)
 comp.setMassFraction(symbols, X)
 comp.finalize(True)
 
+print(f"Initial H-1 mass fraction {comp.getMassFraction("H-1")}")
+
 netIn = NetIn()
 netIn.composition = comp
 netIn.temperature = 1.5e7
-netIn.density = 1.5e2
-netIn.tMax = 3e17
+netIn.density = 1.6e2
+netIn.tMax = 1e-9
 netIn.dt0 = 1e-12
 
-baseEngine = GraphEngine(comp, 2)
-
+baseEngine = GraphEngine(netIn.composition, 2)
 baseEngine.setUseReverseReactions(False)
 
 qseEngine = MultiscalePartitioningEngineView(baseEngine)
+
 adaptiveEngine = AdaptiveEngineView(qseEngine)
 
 solver = DirectNetworkSolver(adaptiveEngine)
-results = solver.evaluate(netIn);
-print(results.composition.getMassFraction("H-1"))
+
+results = solver.evaluate(netIn)
+
+print(f"Final H-1 mass fraction {results.composition.getMassFraction("H-1")}")
+
+