0D validity study

Author: hmonroedd

examples/sparging_standard.py

@@ -23,18 +23,23 @@
 
 # standard_input = get_sim_input_standard()
 my_input = get_sim_input_LIBRA_Pi()
+# my_func = my_input.h_l
+# my_input.h_l = lambda z: my_func(z) * 9.4
 
 
 print(f"Pi = {my_input.get_Pi_number():.2f}")
 print(f"steady state c_T2 = {my_input.get_c_T2_SS():.2e}")
 print(f"Bo = {my_input.get_Bo():.2e}")
 
+tau_ave = my_input.get_tau_ave()
+print(f"tau_ave = {tau_ave:~.2e}")
+
 my_input.c_T2_init = 3e-11 * ureg.molT2 / ureg.m**3
 
 
 my_simulation = Simulation(
     my_input,
-    t_final=3 * ureg.days,
+    t_final=5 * tau_ave,
     dispersion_on=True,
     constant_profiles=False,
 )
@@ -55,15 +60,15 @@
         "h_l",
     ]
 
-    tau_pred = my_input.get_tau()
-    dt = (tau_pred * 0.02).to(
+    tau_pred = my_input.get_tau_ave()
+    dt = (tau_pred * 0.005).to(
         "s"
     )  # 2% of tau_pred gives 1% error on tau compared to fine mesh
-    dx = my_input.dx_from_Pe(2)  # grid Pe = 2
-    my_input.K_s *= 1e-2
+    dx = my_input.dx_from_Pe(0.1)  # grid Pe = 2
+    my_input.K_s *= 0.36
     print(f"dx={dx:~.2e}, dt={dt:~.2e}")
     t_start = time.perf_counter()
-    output = my_simulation.solve(dt=dt, dx=dx, verbose=True)
+    output = my_simulation.solve(dt=dt, dx=dx, verbose=False)
     elapsed = time.perf_counter() - t_start
     print(f"solve() took {elapsed:.1f} s")
 

paper/generate_paper_data.py

@@ -21,23 +21,29 @@
 spatial profiles), plus, for a manual post-run sanity check, the inventory decay
 series of the finest temporal run of each regime.
 
+Also exposes `generate_validity_data`, which produces the data for the 0D
+approximation validity study (see its docstring below).
+
 Run with::
 
     python paper/generate_paper_data.py
 
-Output is written to ``paper/runs/convergence_study/``.
+Output is written to ``paper/runs/convergence_study/`` and
+``paper/runs/0D_validity_study/``.
 """
 
 from __future__ import annotations
 
 import json
 import logging
+import multiprocessing as mp
 import time
 import warnings
 from datetime import datetime
 from pathlib import Path
 
 import numpy as np
+import pandas as pd
 
 from sparging import get_sim_input_LIBRA_Pi, Simulation, ureg
 from sparging import helpers
@@ -231,6 +237,162 @@ def convergence_study(out_dir: Path = OUT_DIR) -> Path:
     return out_path
 
 
+# ---------------------------------------------------------------------------
+# 0D approximation validity study
+# ---------------------------------------------------------------------------
+"""
+Quantifies when the 0D analytical approximation of the extraction time
+(SimulationInput.get_tau / get_tau_ave) is valid, by comparing it against the
+full 1D ARD model over a sampled space, and correlating the disagreement with
+the dimensionless groups governing each simplifying assumption: Pi (negligible
+interfacial partial pressure / solubility), G_mix_pred (perfectly mixed
+liquid) and G_P (constant hydrodynamic parameters along z).
+
+Two things are varied, log-uniformly over +-DECADES around nominal LIBRA-Pi
+values: the product h_l*K_s (via a multiplier applied to a coin-flip choice of
+h_l or K_s -- Pi depends on the product, tau_pred depends only on h_l, so
+comparing the two populations tells apart Pi-controlled behaviour from a pure
+h_l effect), and the tank height H (area held fixed). This is a numerical-
+validity study, not a LIBRA-Pi operating study: sampled values may be
+unphysical for the real experiment, that is fine and intended.
+
+Run with::
+
+    python -c "from paper.generate_paper_data import generate_validity_data; generate_validity_data()"
+
+Output is written to ``paper/runs/0D_validity_study/``.
+"""
+
+OUT_DIR_VALIDITY = Path("paper/runs/0D_validity_study")
+
+N_SAMPLES = 300
+SEED = 42
+DECADES = 2.0  # log-uniform sampling range (+-) for both the h_l*K_s multiplier and H
+T_FINAL_IN_TAU = 2  # simulate 2x tau_pred_ave -- enough to see the decay
+DT_FRACTION_OF_TAU = 0.01  # dt = tau_pred_ave * this
+MESH_PE = 2  # mesh Peclet number for dx (matches examples/sparging_standard.py)
+MIN_CELLS = 20  # floor on n_cells: dx_from_Pe(MESH_PE) does not scale with H, so
+# small-H samples would otherwise plan for very few mesh cells
+OTHER_GROUP_THRESHOLD = 0.1  # "other groups < 0.1" highlighting rule used in the plots
+RMSE_FLAG_THRESHOLD = 1e-5  # normalized-RMSE above this -> non-exponential decay
+
+
+def _make_validity_input(height: "ureg.Quantity", factor: str, multiplier: float):
+    """LIBRA-Pi-like input with height `height` (area fixed at the nominal
+    value) and h_l or K_s scaled by `multiplier`. Both the profile (h_l) and
+    scalar (K_s) overrides are read fresh by every downstream getter and by
+    Simulation.solve(), so the change is automatically consistent everywhere
+    (ARD coefficients, tau_pred, Pi, ...)."""
+    from sparging import (
+        SimulationInput,
+        LIBRA_PI_GEOM,
+        LIBRA_PI_MAT,
+        LIBRA_PI_OPERATING_PARAMS,
+        LIBRA_PI_SPARGING_PARAMS,
+    )
+
+    geom = LIBRA_PI_GEOM.copy()
+    geom.height = height
+    inp = SimulationInput.from_parameters(
+        geom,
+        LIBRA_PI_MAT.copy(),
+        LIBRA_PI_OPERATING_PARAMS.copy(),
+        LIBRA_PI_SPARGING_PARAMS.copy(),
+    )
+    inp.c_T2_init = C_T2_INIT
+    if factor == "h_l":
+        base_h_l = inp.h_l
+        inp.h_l = lambda z, _f=base_h_l, _m=multiplier: _f(z) * _m
+    else:
+        inp.K_s = inp.K_s * multiplier
+    return inp
+
+
+def _draw_sample_specs(n_samples: int, seed: int) -> list[dict]:
+    """Draw the (multiplier, factor, height) sample space. Plain-dict specs
+    (only picklable primitives) so they can be sent to a worker process."""
+    rng = np.random.default_rng(seed)
+    H_nominal = get_sim_input_LIBRA_Pi().height.to("m").magnitude
+    return [
+        {
+            "sample_id": i,
+            "multiplier": float(10 ** rng.uniform(-DECADES, DECADES)),
+            "factor": str(rng.choice(["h_l", "K_s"])),
+            "height_m": float(H_nominal * 10 ** rng.uniform(-DECADES, DECADES)),
+        }
+        for i in range(n_samples)
+    ]
+
+
+def _run_one_sample(spec: dict) -> dict:
+    """Build the input from `spec` (not a pre-built SimulationInput: profile
+    closures aren't picklable, so a spawned worker process must rebuild the
+    input itself), run the 1D ARD model once, and return one row of the
+    validity dataset (see SimulationResults.validity_record)."""
+    warnings.filterwarnings("ignore")  # silence curve_fit / pint fit warnings
+
+    inp = _make_validity_input(
+        spec["height_m"] * ureg.m, spec["factor"], spec["multiplier"]
+    )
+
+    tau_pred_ave = inp.get_tau_ave()
+    dt = (tau_pred_ave * DT_FRACTION_OF_TAU).to("s")
+    t_final = (T_FINAL_IN_TAU * tau_pred_ave).to("s")
+    n_cells = max(
+        int(round((inp.height / inp.dx_from_Pe(MESH_PE)).to("dimensionless").magnitude)),
+        MIN_CELLS,
+    )
+    dx = (inp.height / n_cells).to("m")
+
+    sim = Simulation(inp, t_final=t_final, dispersion_on=True, constant_profiles=False)
+    out = sim.solve(dt=dt, dx=dx)
+    return out.validity_record(sample=spec, t_0=0 * ureg.s, t_end=t_final)
+
+
+def generate_validity_data(
+    n_samples: int = N_SAMPLES,
+    seed: int = SEED,
+    n_workers: int = 4,
+    out_dir: Path = OUT_DIR_VALIDITY,
+) -> Path:
+    """Run the 0D-approximation validity study: one 1D ARD solve per sampled
+    (h_l*K_s multiplier, height) point, in parallel, and write the results to
+    a single CSV. Spawn-based multiprocessing: each worker rebuilds its own
+    SimulationInput from a picklable spec, which sidesteps both the
+    unpicklable profile closures and dolfinx/mpi4py process-fork hazards.
+    """
+    out_dir.mkdir(parents=True, exist_ok=True)
+    specs = _draw_sample_specs(n_samples, seed)
+
+    t_start = time.time()
+    with mp.get_context("spawn").Pool(n_workers) as pool:
+        rows = list(pool.imap_unordered(_run_one_sample, specs))
+    elapsed = time.time() - t_start
+    logger.info("0D validity study: %d samples in %.1f s", len(rows), elapsed)
+
+    df = pd.DataFrame(rows).sort_values("sample_id").reset_index(drop=True)
+    csv_path = out_dir / "validity_data.csv"
+    df.to_csv(csv_path, index=False)
+
+    with open(out_dir / "metadata.json", "w") as f:
+        json.dump(
+            {
+                "git_commit": helpers.get_git_hash(),
+                "date": datetime.now().isoformat(),
+                "n_samples": len(df),
+                "decades": DECADES,
+                "t_final_in_tau": T_FINAL_IN_TAU,
+                "dt_fraction_of_tau": DT_FRACTION_OF_TAU,
+                "mesh_Pe": MESH_PE,
+                "other_group_threshold": OTHER_GROUP_THRESHOLD,
+                "rmse_flag_threshold": RMSE_FLAG_THRESHOLD,
+            },
+            f,
+            indent=2,
+        )
+    return csv_path
+
+
 if __name__ == "__main__":
     logging.basicConfig(
         level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s"

src/sparging/input_examples.py

@@ -110,7 +110,8 @@ def get_sim_input_LIBRA_Pi() -> SimulationInput:
     """
     geom = ColumnGeometry(
         area=0.13 * ureg.m**2,
-        height=0.93 * ureg.m,
+        # height=0.93 * ureg.m,
+        height=0.43 * ureg.m,
         nozzle_diameter=0.002 * ureg.m,
         nb_nozzle=4 * ureg.dimensionless,
     )

src/sparging/inputs.py

@@ -172,22 +172,23 @@ def get_tau(self) -> pint.Quantity:
         """characteristic time of the sparger under the small partial pressure (SPP) approximation"""
         return (self.eps_l0 / (self.h_l0 * self.a_0)).to("seconds")
 
-    def get_tau_ave(self) -> pint.Quantity:
-        """characteristic time computed from a, h_l and eps_l profiles integrated
-        over the tank height, instead of evaluated at the bottom (z=0)."""
+    def _height_integral(
+        self, profile: Callable[[pint.Quantity], pint.Quantity], unit: str
+    ) -> pint.Quantity:
+        """integral of a z-dependent profile (evaluated in `unit`) over [0, height]."""
         from scipy.integrate import quad
 
-        def eps_l_integrand(z_m: float) -> float:
-            return (1 - self.eps_g(z_m * ureg.m)).to("dimensionless").magnitude
-
-        def a_h_l_integrand(z_m: float) -> float:
-            z = z_m * ureg.m
-            return (self.a(z) * self.h_l(z)).to("1/s").magnitude
+        def integrand(z_m: float) -> float:
+            return profile(z_m * ureg.m).to(unit).magnitude
 
         height_m = self.height.to("m").magnitude
-        int_eps_l = quad(eps_l_integrand, 0, height_m)[0] * ureg.m
-        int_a_h_l = quad(a_h_l_integrand, 0, height_m)[0] * ureg("m/s")
+        return quad(integrand, 0, height_m)[0] * ureg(unit) * ureg.m
 
+    def get_tau_ave(self) -> pint.Quantity:
+        """characteristic time computed from a, h_l and eps_l profiles integrated
+        over the tank height, instead of evaluated at the bottom (z=0)."""
+        int_eps_l = self._height_integral(lambda z: 1 - self.eps_g(z), "dimensionless")
+        int_a_h_l = self._height_integral(lambda z: self.a(z) * self.h_l(z), "1/s")
         return (int_eps_l / int_a_h_l).to("seconds")
 
     def get_c_T2_SS(self) -> pint.Quantity:
@@ -207,6 +208,36 @@ def get_Pi_number(self) -> pint.Quantity:
             / (self.eps_g0 * self.graph.nodes["v_g0"]["value"])
         ).to("dimensionless")
 
+    def get_Pi_ave(self) -> pint.Quantity:
+        """Partial pressure number computed from a, h_l and u_g profiles integrated
+        over the tank height, instead of evaluated at the bottom (z=0). Consistent
+        with get_tau_ave() (same height-averaged a*h_l)."""
+        int_a_h_l = self._height_integral(lambda z: self.a(z) * self.h_l(z), "1/s")
+        int_u_g = self._height_integral(self.u_g, "m/s")
+        a_h_l_ave = int_a_h_l / self.height
+        u_g_ave = int_u_g / self.height
+        return (
+            self.K_s
+            * (const_R * self.temperature)
+            * self.height
+            * a_h_l_ave
+            / u_g_ave
+        ).to("dimensionless")
+
+    def get_G_mix_pred(self) -> pint.Quantity:
+        """Predicted mixing number: ratio of the liquid dispersive mixing time
+        scale (H^2/E_l) to the height-averaged extraction time tau_pred_ave.
+        Governs the "perfectly mixed liquid" assumption. "Predicted" since it
+        uses the analytical tau_pred_ave rather than a simulated/fitted tau."""
+        return ((self.height**2 / self.E_l) / self.get_tau_ave()).to("dimensionless")
+
+    def get_G_P(self) -> pint.Quantity:
+        """Relative hydrostatic pressure variation along the tank height
+        (P_bottom - P_top) / P_top. Governs the "constant hydrodynamic
+        parameters along z" assumption."""
+        P_top = self.P_l(self.height)
+        return ((self.P_l0 - P_top) / P_top).to("dimensionless")
+
     def get_dP_dx(self) -> pint.Quantity:
         """
         returns 1/c_T * dP_T/dx in the SPP approximation (linearized around P_T = 0))