From e98c9a4050eb728c96ff3316eca374b3db607ba8 Mon Sep 17 00:00:00 2001
From: Emily Boudreaux <thomas@boudreauxmail.com>
Date: Sat, 20 Dec 2025 16:08:33 -0500
Subject: [PATCH] feat(python): Python multi-test

Python now works with mulit-threading and zones
---
 subprojects/fourdst.wrap                    |  2 +-
 tests/python/py_test_multi.py               | 79 +++++++++++++++++++++
 tests/python/{test.py => py_test_single.py} |  0
 3 files changed, 80 insertions(+), 1 deletion(-)
 create mode 100644 tests/python/py_test_multi.py
 rename tests/python/{test.py => py_test_single.py} (100%)

diff --git a/subprojects/fourdst.wrap b/subprojects/fourdst.wrap
index 73a4d422..5ab46be9 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.17
+revision = v0.9.18
 depth = 1
diff --git a/tests/python/py_test_multi.py b/tests/python/py_test_multi.py
new file mode 100644
index 00000000..8b8188db
--- /dev/null
+++ b/tests/python/py_test_multi.py
@@ -0,0 +1,79 @@
+from gridfire.solver import GridSolver, PointSolver, GridSolverContext
+from gridfire.policy import MainSequencePolicy
+from fourdst.composition import Composition
+
+from fourdst.composition import CanonicalComposition
+from fourdst.atomic import Species
+from gridfire.type import NetIn
+
+import numpy as np
+
+def rescale_composition(comp_ref : Composition, ZZs : float, Y_primordial : float = 0.248) -> Composition:
+    CC : CanonicalComposition  = comp_ref.getCanonicalComposition()
+
+    dY_dZ = (CC.Y - Y_primordial) / CC.Z
+
+    Z_new = CC.Z * (10**ZZs)
+    Y_bulk_new = Y_primordial + (dY_dZ * Z_new)
+    X_new = 1.0 - Z_new - Y_bulk_new
+
+    if X_new < 0: raise ValueError(f"ZZs={ZZs} yields unphysical composition (X < 0)")
+
+    ratio_H  = X_new / CC.X if CC.X > 0 else 0
+    ratio_He = Y_bulk_new / CC.Y if CC.Y > 0 else 0
+    ratio_Z  = Z_new / CC.Z if CC.Z > 0 else 0
+
+    Y_new_list = []
+    newComp : Composition = Composition()
+    s: Species
+    for s  in comp_ref.getRegisteredSpecies():
+        Xi_ref = comp_ref.getMassFraction(s)
+
+        if s.el() == "H":
+            Xi_new = Xi_ref * ratio_H
+        elif s.el() == "He":
+            Xi_new = Xi_ref * ratio_He
+        else:
+            Xi_new = Xi_ref * ratio_Z
+
+        Y = Xi_new / s.mass()
+        newComp.registerSpecies(s)
+        newComp.setMolarAbundance(s, Y)
+
+    return newComp
+
+def init_composition(ZZs : float = 0) -> Composition:
+    Y_solar = [7.0262E-01, 9.7479E-06, 6.8955E-02, 2.5000E-04, 7.8554E-05, 6.0144E-04, 8.1031E-05, 2.1513E-05]
+    S = ["H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"]
+    return rescale_composition(Composition(S, Y_solar), ZZs)
+
+
+def init_netIn(temp: float, rho: float, time: float, comp: Composition) -> NetIn:
+    n : NetIn = NetIn()
+    n.temperature = temp
+    n.density = rho
+    n.tMax = time
+    n.dt0 = 1e-12
+    n.composition = comp
+    return n
+
+def years_to_seconds(years: float) -> float:
+    return years * 3.1536e7
+
+
+def main():
+    C = init_composition()
+    temps = np.linspace(1.5e7, 2e7, 100)
+    rhos = np.linspace(1.5e2, 1.5e2, 100)
+    netIns = []
+    for T, R in zip(temps, rhos):
+        netIns.append(init_netIn(T, R, years_to_seconds(100e6), C))
+    policy = MainSequencePolicy(C)
+    construct = policy.construct()
+    local_solver = PointSolver(construct.engine)
+    grid_solver = GridSolver(construct.engine, local_solver)
+    solver_ctx = GridSolverContext(construct.scratch_blob)
+    results = grid_solver.evaluate(solver_ctx, netIns)
+
+if __name__ == "__main__":
+    main()
diff --git a/tests/python/test.py b/tests/python/py_test_single.py
similarity index 100%
rename from tests/python/test.py
rename to tests/python/py_test_single.py