added postprocessing.py + implemented pint in SimulationResults

Author: hmonroedd

src/sparging/animation.py

@@ -36,17 +36,21 @@ def __init__(
         hspace : float, optional
             Vertical spacing between subplots
         """
-        self.times_hr = np.array(results.times) * sec_to_hour
-        self.c_T2_solutions = np.array(results.c_T2_solutions)
-        self.y_T2_solutions = np.array(results.y_T2_solutions)
-        self.x_ct = results.x_ct
-        self.x_y = results.x_y
-        self.inventories_T2_salt = results.inventories_T2_salt
+        self.times_hr = results.times.to("hour").magnitude
+        self.c_T2_solutions = results.c_T2_solutions.to("molT2/m^3").magnitude
+        self.y_T2_solutions = results.y_T2_solutions.to("dimensionless").magnitude
+        self.x_ct = results.x_ct.to("m").magnitude
+        self.x_y = results.x_y.to("m").magnitude
+        self.inventories_T2_salt = results.inventories_T2_salt.to("molT2").magnitude
         self.source_T2 = (
-            None if results.sources_T2 is None else np.array(results.sources_T2)
+            None
+            if results.sources_T2 is None
+            else results.sources_T2.to("molT2/s").magnitude
         )
         self.fluxes_T2 = (
-            None if results.fluxes_T2 is None else np.array(results.fluxes_T2)
+            None
+            if results.fluxes_T2 is None
+            else results.fluxes_T2.to("molT2/s").magnitude
         )
         self.show_activity = show_activity
         self.figsize = figsize

src/sparging/inputs.py

@@ -385,7 +385,7 @@ def get_sim_input_standard() -> SimulationInput:
     "height": 1 * ureg.m,
     "nozzle_diameter": 1.4 * ureg.mm,
     "nb_nozzle": 5 * ureg.dimensionless,
-    "temperature": 500 * ureg.celsius,
+    "temperature": 600 * ureg.celsius,
     "P_top": 1 * ureg.atm,
     "flow_g_mol": 100 * ureg.sccm,
     "tbr": 0.1 * ureg("triton / neutron"),

src/sparging/model.py

@@ -39,19 +39,19 @@
 
 
 @dataclass
-class SimulationResults:  # TODO implement pint in this class # TODO change list to np.array
-    times: list
-    c_T2_solutions: list
-    y_T2_solutions: list
-    J_T2_solutions: list
-    x_ct: np.ndarray
-    x_y: np.ndarray
-    inventories_T2_salt: np.ndarray
-    sources_T2: list
-    fluxes_T2: list
+class SimulationResults:
+    times: np.ndarray[pint.Quantity]
+    c_T2_solutions: np.ndarray[pint.Quantity]
+    y_T2_solutions: np.ndarray[pint.Quantity]
+    J_T2_solutions: np.ndarray[pint.Quantity]
+    x_ct: np.ndarray[pint.Quantity]
+    x_y: np.ndarray[pint.Quantity]
+    inventories_T2_salt: np.ndarray[pint.Quantity]
+    sources_T2: np.ndarray[pint.Quantity]
+    fluxes_T2: np.ndarray[pint.Quantity]
     sim_input: SimulationInput
-    dt: int | None = None
-    dx: int | None = None
+    dt: pint.Quantity
+    dx: pint.Quantity
 
     keys_to_ignore_results = [  # TODO do it the other way: keys_to_include_results
         "c_T2_solutions",
@@ -428,22 +428,20 @@ def post_process(t):
             # advance time
             t += dt
 
-        inventories_T2_salt = np.array(inventories_T2_salt)
-
         # TODO reattach units using wrapping
         # https://pint.readthedocs.io/en/stable/advanced/performance.html#a-safer-method-wrapping
         results = SimulationResults(
-            times=times,
-            c_T2_solutions=c_T2_solutions,
-            y_T2_solutions=y_T2_solutions,
-            J_T2_solutions=J_T2_solutions,
-            x_ct=x_ct,
-            x_y=x_y,
-            inventories_T2_salt=inventories_T2_salt,
-            sources_T2=sources_T2,
-            fluxes_T2=fluxes_T2,
+            times=np.array(times) * ureg("s"),
+            c_T2_solutions=np.array(c_T2_solutions) * ureg("molT2/m^3"),
+            y_T2_solutions=np.array(y_T2_solutions) * ureg("dimensionless"),
+            J_T2_solutions=np.array(J_T2_solutions) * ureg("molT2/m^3/s"),
+            x_ct=x_ct * ureg("m"),
+            x_y=x_y * ureg("m"),
+            inventories_T2_salt=np.array(inventories_T2_salt) * ureg("molT2"),
+            sources_T2=np.array(sources_T2) * ureg("molT2/s"),
+            fluxes_T2=np.array(fluxes_T2) * ureg("molT2/s"),
             sim_input=self.sim_input,
-            dt=dt,
-            dx=dx,
+            dt=dt * ureg("s"),
+            dx=dx * ureg("m"),
         )
         return results

src/sparging/postprocess.py

@@ -0,0 +1,50 @@
+from __future__ import annotations
+from sparging.config import ureg
+from typing import TYPE_CHECKING
+import numpy as np
+
+if TYPE_CHECKING:
+    import pint
+    from sparging.model import SimulationResults
+    import matplotlib.pyplot as plt
+
+
+def indices_from_times(sim_output: SimulationResults, times_ref: list[pint.Quantity]):
+    i_matching = [np.argmin(np.abs(sim_output.times - t)) for t in times_ref]
+    times_matching = [sim_output.times[i] for i in i_matching]
+    return i_matching, times_matching
+
+
+def plot_profile(
+    sim_output: SimulationResults,
+    ax: plt.Axes,
+    var_name: str,
+    times: list[pint.Quantity],
+):
+    i_to_plot, _ = indices_from_times(sim_output, times)
+    y_to_plot = getattr(sim_output, var_name)
+    for i in i_to_plot:
+        ax.plot(
+            sim_output.x_ct,
+            y_to_plot[i],
+            label=f"t = {sim_output.times[i] / 3600:.2f} h",
+        )
+    ax.legend()
+    ax.set_xlabel(r"z [m]")
+    ax.set_ylabel(var_name)  # TODO leverage pint for units
+    ax.set_title(var_name + " profile")
+    ax.grid()
+
+
+def plot_signal(sim_output: SimulationResults, ax: plt.Axes, var_name: str):
+    y_to_plot = getattr(sim_output, var_name)
+    ax.plot(
+        sim_output.times,
+        y_to_plot,
+        label=var_name,
+    )
+    ax.legend()
+    ax.set_xlabel(r"t [s]")
+    ax.set_ylabel(var_name)  # TODO leverage pint for units
+    ax.set_title(var_name + " profile")
+    ax.grid()