post processing functions

Co-authored-by: Copilot <copilot@github.com>

Author: hmonroedd

sparging101.py

@@ -4,10 +4,9 @@
 from sparging import animation
 from sparging.model import Simulation
 from sparging.inputs import (
-    get_sim_input,
-    librapi_input_dict,
     get_sim_input_standard,
 )
+import sparging.postprocess as pp
 import logging
 from typing import TYPE_CHECKING
 import numpy as np
@@ -30,8 +29,8 @@
         my_input.get_c_T2_SS().to('molT2/m^3')
     }"
 )
-T_99 = (np.log(100) * my_input.get_tau()).to("seconds")
-print(f"T_99% = {T_99.to('hours')}")
+tau = my_input.get_tau().to("seconds")
+print(f"tau = {tau.to('hours')}")
 print(f"Partial pressure number PP = {my_input.get_PP_number()}")
 
 
@@ -48,32 +47,30 @@ def profile_source_T(z: pint.Quantity | list[float], height: pint.Quantity = Non
         return np.pi / 2 * np.sin(np.pi / height * z)  # normalized
     else:
         raise NotImplementedError("z must be either a float or a pint.Quantity")
-    # return 0.5 * (1 + np.cos(0.5 * np.pi / (1 * ureg.m) * z))
 
 
-my_input.profile_source_T = profile_source_T
-my_input.signal_irr = lambda t: 1 if t < T_99 else 0
+t_final = 50 * ureg.hour
+t_irr = 20 * ureg.hour
+my_input.signal_irr = lambda t: 1 if t <= t_irr else 0
+my_input.signal_sparging = lambda t: 1
 
 my_simulation = Simulation(
     my_input,
-    t_final=2 * T_99,
+    t_final=t_final,
     profile_pressure_hydrostatic=True,
 )
 
 if __name__ == "__main__":
     # my_simulation.sim_input.E_g *= 1e5
     # my_simulation.sim_input.E_l *= 1e-5
     output = my_simulation.solve(fast_solve=True)
-
-    # # save output to file
-    # output.profiles_to_csv(f"output_{tank_height}m.csv")
-
-    # # plot results
-    # from sparging import plotting
-    # plotting.plot_animation(output)
-
-    # import matplotlib.pyplot as plt
-
-    # plt.plot(output.times, output.inventories_T2_salt)
-    # plt.show()
+    output.to_json("temp.json")
+    popt, pcov = pp.fit_exp(
+        output.inventories_T2_salt, output.times, 0 * ureg.s, t_irr, phase="rampup"
+    )
+    print(f"Fitted parameters: n0 = {popt[1]}, tau = {popt[0].to('h')}")
+    popt, pcov = pp.fit_exp(
+        output.inventories_T2_salt, output.times, t_irr, t_final, phase="decay"
+    )
+    print(f"Fitted parameters: n0 = {popt[1]}, tau = {popt[0].to('h')}")
     animation.create_animation(output, show_activity=False)

src/sparging/config.py

@@ -5,6 +5,7 @@
 ureg = UnitRegistry(
     autoconvert_offset_to_baseunit=True
 )  # to deal with offset units (eg: degree celsius)
+ureg.setup_matplotlib(True)
 ureg.formatter.default_format = ".3e~D"
 ureg.define("triton = [tritium] = T")
 ureg.define(f"molT = {const.N_A} * triton")

src/sparging/model.py

@@ -226,7 +226,7 @@ def solve(
         tank_volume = self.sim_input.volume.to("m**3").magnitude
         a = self.sim_input.a.to("1/m").magnitude
         h_l = self.sim_input.h_l.to("m/s").magnitude
-        K_s = self.sim_input.K_s.to("mol/m**3/Pa").magnitude
+        K_s = self.sim_input.K_s.to("mol/m**3/Pa").magnitude  # convert to molT2 ?
         P_0 = self.sim_input.P_bottom.to("Pa").magnitude
         T = self.sim_input.temperature.to("K").magnitude
         eps_g = self.sim_input.eps_g.to("dimensionless").magnitude
@@ -474,6 +474,8 @@ def post_process(t):
 
         # SOLVE
         while t < t_final:
+            # advance time
+            t += dt
             # update time-dependent terms
             gen_T2_ave.value = (
                 Q_T2 / tank_volume * self.sim_input.signal_irr(t * ureg.s)
@@ -487,9 +489,6 @@ def post_process(t):
 
             post_process(t)
 
-            # advance time
-            t += dt
-
         # TODO reattach units using wrapping
         # https://pint.readthedocs.io/en/stable/advanced/performance.html#a-safer-method-wrapping
         results = SimulationResults(

src/sparging/postprocess.py

@@ -2,19 +2,23 @@
 from sparging.config import ureg
 from typing import TYPE_CHECKING
 import numpy as np
+import warnings
+from scipy.optimize import curve_fit
 
 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]
-
-    assert len(i_matching) == len(times_ref)
-    return i_matching, times_matching
+def idx_from_t(times: list[pint.Quantity], timestamp: pint.Quantity) -> pint.Quantity:
+    idx = np.argmin(np.abs(times - timestamp))
+    time = times[idx].to(timestamp.units)
+    if not np.isclose(times[idx], timestamp):
+        warnings.warn(
+            f"Requested time {timestamp} does not exactly match discrete times. Using closest time {time} for index calculation."
+        )
+    return idx
 
 
 def plot_profile(
@@ -25,26 +29,76 @@ def plot_profile(
     colors: list = None,
     **kwargs,
 ):
-    i_to_plot, _ = indices_from_times(sim_output, times)
+    idx_to_plot = [idx_from_t(sim_output.times, t) for t in times]
     y_to_plot = getattr(sim_output, var_name)
-    for j, idx in enumerate(i_to_plot):
+    for j, idx in enumerate(idx_to_plot):
         ax.plot(
             sim_output.x_ct,
             y_to_plot[idx],
             label=f"t = {sim_output.times[idx]:.2f}",
             c=colors[j] if colors is not None else None,
             **kwargs,
         )
-    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, **kwargs):
-    ax.plot(sim_output.times, getattr(sim_output, var_name), **kwargs)
-    # ax.set_xlabel(r"t [s]")
-    # ax.set_ylabel(var_name)  # TODO leverage pint for units
-    # ax.set_title(var_name + " profile")
-    # ax.grid()
+
+
+# TODO remove
+# def plot_signal(sim_output: SimulationResults, ax: plt.Axes, var_name: str, **kwargs):
+#     ax.plot(sim_output.times, getattr(sim_output, var_name), **kwargs)
+#     # ax.set_xlabel(r"t [s]")
+#     # ax.set_ylabel(var_name)  # TODO leverage pint for units
+#     # ax.set_title(var_name + " profile")
+#     # ax.grid()
+
+
+def get_residual_fraction(
+    T2_inventories: np.ndarray[pint.Quantity],
+    times: np.ndarray[pint.Quantity],
+    t_0: pint.Quantity,
+    t_end: pint.Quantity,
+) -> pint.Quantity:
+    i_0 = idx_from_t(times, t_0)
+    i_end = idx_from_t(times, t_end)
+    return T2_inventories[i_end] / T2_inventories[i_0]
+
+
+def fit_exp(
+    vec: np.ndarray[pint.Quantity],
+    times: np.ndarray[pint.Quantity],
+    t_0: pint.Quantity,
+    t_end: pint.Quantity,
+    phase: str,
+):  # TODO add warning if error or std is too large
+    """
+    - phase = 'decay' or 'rampup'
+    """
+
+    def fitting_func(t, tau, n0):
+        match phase:
+            case "decay":
+                return n0 * np.exp(-t / tau)
+            case "rampup":
+                return n0 * (1 - np.exp(-t / tau))
+            case _:
+                raise ValueError("Invalid phase. Must be 'decay' or 'rampup'.")
+
+    idx_0 = idx_from_t(times, t_0)
+    idx_end = idx_from_t(times, t_end)
+    t_0 = times[idx_0]
+    t_end = times[idx_end]
+    print(
+        f"Fitting from t={t_0.to('hour')} to t={t_end.to('hour')} (indices {idx_0} to {idx_end})"
+    )
+    tau_guess = 10000 * ureg.s
+    n0_guess = vec[idx_0] if phase == "decay" else vec[idx_end]
+
+    # wrapped_fitting_func = ureg.wraps(vec.units, (None, vec.units, "s"))(fitting_func)
+
+    popt, pcov = curve_fit(
+        fitting_func,
+        (times[idx_0 : idx_end + 1] - t_0).to("s").magnitude,
+        vec[idx_0 : idx_end + 1].magnitude,
+        p0=[tau_guess.to("s").magnitude, n0_guess.magnitude],
+    )
+    print(f"std:{np.sqrt(np.diag(pcov))}")
+
+    return (popt[0] * ureg.s, popt[1] * vec.units), pcov