Generic salt/feed profiles, Jacobian sparsity, unpack helper

Author: mironovb

chromosim/models/params_ev_aex.py

@@ -0,0 +1,64 @@
+from __future__ import annotations
+from types import SimpleNamespace as NS
+from ..profiles.salt_and_feed import salt_profile_clamped_linear, feed_profile_by_load_volume
+
+def params_ev_aex():
+    """
+    Multi-species parameter set matching your MATLAB 'params_ev_aex.m'.
+    Returns a SimpleNamespace so attributes can be attached.
+    """
+    par = NS()
+    # Geometry / grid
+    par.L    = 0.025         # m
+    par.dID  = 0.007         # m
+    par.eps  = 0.40
+    par.N    = 80
+    par.dz   = par.L / par.N
+
+    # Hydro
+    par.u    = 1.0e-4        # m s^-1
+    par.Dax  = 8e-11         # m^2 s^-1
+
+    # EV classes
+    par.zeta_mV = [-12, -25, -45]
+    par.frac    = [0.10, 0.15, 0.75]
+    par.Nsp     = len(par.zeta_mV)
+
+    # Scaling / capacity
+    par.Cfeed_total = 7.5e9        # particles m^-3
+    par.qmax        = 4e13         # particles m^-3 wall
+    par.Cscale      = par.Cfeed_total
+    par.Qscale      = par.qmax
+    par.CplotScale  = 1e10
+    par.CplotLabel  = r"10^{10} per m^{-3}"
+
+    # Kinetics & zeta mapping
+    par.kd0   = 2e-3
+    par.K0    = 1e-22
+    par.gamma = 0.53
+    par.Imax  = 1.2
+    par.fI    = lambda I: max(0.0, 1.0 - float(I)/par.Imax)
+
+    # Gradient
+    par.grad_start = 460.0
+    par.grad_end   = par.grad_start + 1200.0
+    par.I_load     = 0.05
+    par.I_elute    = 1.0
+    par.salt_profile = lambda t: salt_profile_clamped_linear(t, par)
+
+    # Feed (physical units)
+    # flow & load (optional but enables volume-based load windows)
+    par.flow_mL_min = 1.0
+    par.load_mL     = 5.0
+    par.feed_profile = lambda t: feed_profile_by_load_volume(t, par)
+
+    # Time / solver tolerances (used by runner)
+    par.tspan = (0.0, 3000.0)       # s
+    par.rtol  = 1e-6
+    par.atol  = 1e-6
+    par.max_step = 8.0
+
+    # Plot helper
+    par.B_start = 10.0
+    par.B_end   = 100.0
+    return par

chromosim/models/params_ev_aex_single.py

@@ -0,0 +1,61 @@
+from __future__ import annotations
+from types import SimpleNamespace as NS
+from ..profiles.salt_and_feed import salt_profile_clamped_linear, feed_profile_by_load_volume
+
+def params_ev_aex_single():
+    """Single-species parameter set matching your MATLAB 'params_ev_aex_single.m'."""
+    par = NS()
+    # Geometry / grid
+    par.L    = 0.025
+    par.dID  = 0.007
+    par.eps  = 0.40
+    par.N    = 80
+    par.dz   = par.L / par.N
+
+    # Hydro
+    par.u    = 1.0e-4
+    par.Dax  = 9e-11
+
+    # EV class (single)
+    par.zeta_mV = [-25]
+    par.frac    = [1.0]
+    par.Nsp     = 1
+
+    # Scaling / capacity
+    par.Cfeed_total = 5e9
+    par.qmax        = 4e14
+    par.Cscale      = par.Cfeed_total
+    par.Qscale      = par.qmax
+    par.CplotScale  = 1e10
+    par.CplotLabel  = r"10^{10} per m^{-3}"
+
+    # Kinetics & zeta mapping
+    par.kd0   = 1e-2
+    par.K0    = 1e-22
+    par.gamma = 0.750
+    par.Imax  = 1.2
+    par.fI    = lambda I: max(0.0, 1.0 - float(I)/par.Imax)
+
+    # Gradient
+    par.grad_start = 12.0*60.0
+    par.grad_end   = par.grad_start + 20.0*60.0
+    par.I_load     = 0.01
+    par.I_elute    = 1.0
+    par.salt_profile = lambda t: salt_profile_clamped_linear(t, par)
+
+    # Feed (physical units)
+    # flow & load for a clear 5 mL at 1 mL/min → 300 s load
+    par.flow_mL_min = 1.0
+    par.load_mL     = 5.0
+    par.feed_profile = lambda t: feed_profile_by_load_volume(t, par)
+
+    # Time / solver tolerances
+    par.tspan = (0.0, 45.0*60.0)
+    par.rtol  = 1e-6
+    par.atol  = 1e-6
+    par.max_step = 8.0
+
+    # Plot helper
+    par.B_start = 10.0
+    par.B_end   = 100.0
+    return par

chromosim/numerics/build_matrices.py

@@ -1,42 +1,58 @@
 from __future__ import annotations
+from collections.abc import Mapping
 import numpy as np
-from scipy.sparse import diags, csc_matrix
+from scipy.sparse import diags
+from .sparsity import sparsity_pattern_multi
+
+def _get(par, key):
+    """Get a parameter from either a Mapping (dict) or an object with attrs."""
+    if isinstance(par, Mapping):
+        return par[key]
+    return getattr(par, key)
+
+def _set(par, key, value):
+    """Set a parameter on either a Mapping (dict) or an object with attrs."""
+    if isinstance(par, Mapping):
+        par[key] = value
+    else:
+        setattr(par, key, value)
 
 def build_matrices(par):
     """
     Assemble finite-volume derivative matrices in z and attach to `par`.
-      par.N, par.dz must exist.
+    Requires: par.N (int), par.dz (float).
     Adds:
       par.D1  -- first derivative (backward/upwind)  [N x N] sparse
       par.D2  -- second derivative (central)         [N x N] sparse
     """
-    N  = int(getattr(par, "N", par["N"]))
-    dz = float(getattr(par, "dz", par["dz"]))
+    N  = int(_get(par, "N"))
+    dz = float(_get(par, "dz"))
 
     e = np.ones(N)
-    # Backward/upwind first derivative
-    D1 = diags([-e, e], offsets=[-1, 0], shape=(N, N), format="csc") / dz
-    D1 = D1.tolil()
+
+    # 1) Backward/upwind first derivative
+    D1 = diags([-e, e], offsets=[-1, 0], shape=(N, N), format="lil") / dz
     D1[0, :] = 0  # inlet row handled via boundary condition
     D1 = D1.tocsc()
 
-    # Central second derivative
-    D2 = diags([e, -2*e, e], offsets=[-1, 0, 1], shape=(N, N), format="csc") / (dz**2)
-    D2 = D2.tolil()
-    D2[0, 0] = -2.0 / (dz**2)  # ghost for inlet Dirichlet
-    D2[-1, -1] = 1.0 / (dz**2) # Danckwerts outlet
+    # 2) Central second derivative
+    D2 = diags([e, -2.0 * e, e], offsets=[-1, 0, 1], shape=(N, N), format="lil") / (dz**2)
+    D2[0, 0]   = -2.0 / (dz**2)  # ghost for inlet Dirichlet
+    D2[-1, -1] =  1.0 / (dz**2)  # Danckwerts outlet
     D2 = D2.tocsc()
 
-    setattr(par, "D1", D1)
-    setattr(par, "D2", D2)
+    _set(par, "D1", D1)
+    _set(par, "D2", D2)
     return par
 
-# Optional helper to attach Jacobian sparsity later
-from .sparsity import jacobian_sparsity_multi
-
 def attach_jacobian_sparsity(par):
-    N  = int(getattr(par, "N", par["N"]))
-    Ns = int(getattr(par, "Nsp", par["Nsp"]))
-    J  = jacobian_sparsity_multi(N, Ns)
-    setattr(par, "Jpattern", J)
+    """
+    Attach Jacobian sparsity pattern for multi-species model.
+    Requires: par.N, par.Nsp.
+    Adds: par.Jpattern
+    """
+    N   = int(_get(par, "N"))
+    Nsp = int(_get(par, "Nsp"))
+    Jpat = sparsity_pattern_multi(N, Nsp)
+    _set(par, "Jpattern", Jpat)
     return par

chromosim/numerics/sparsity.py

@@ -1,25 +1,28 @@
 from __future__ import annotations
 import numpy as np
-from scipy.sparse import diags, eye, bmat, csc_matrix
+from scipy.sparse import spdiags, kron, eye, bmat, csr_matrix
 
-def jacobian_sparsity_multi(N: int, Ns: int) -> csc_matrix:
+def sparsity_pattern_multi(N: int, Ns: int):
     """
-    2*N*Ns x 2*N*Ns logical Jacobian sparsity for multi-species column model.
-
-    State layout (row-major):
-      y = [ c(0,0..Ns-1), ..., c(N-1,0..Ns-1),
-            qhat(0,0..Ns-1), ..., qhat(N-1,0..Ns-1) ]
+    Boolean sparse pattern for the Jacobian of size 2*N*Ns × 2*N*Ns.
+    State ordering:
+        [ c(1..N, species 1..Ns) ; qhat(1..N, species 1..Ns) ]
+    Blocks:
+        dc/dc  ~ tri-diagonal in z for each species
+        dc/dq  ~ diagonal (local kinetics)
+        dq/dc  ~ diagonal (local kinetics)
+        dq/dq  ~ diagonal
     """
     e = np.ones(N)
-    T = diags([e, e, e], offsets=[-1, 0, 1], shape=(N, N), format="csc")
+    T = spdiags([e, e, e], [-1, 0, 1], N, N, format='csr')  # tri-diagonal
+    # replicate per species
+    Ac = kron(eye(Ns, format='csr'), (T != 0), format='csr')  # Nc×Nc, Nc=N*Ns
 
-    # replicate tri-diagonal axial stencil per species
-    Ac = bmat([[T if s == r else None for s in range(Ns)]
-               for r in range(Ns)], format="csc")
     Nc = N * Ns
-    I  = eye(Nc, format="csc")
+    B  = eye(Nc, format='csr', dtype=bool)    # diagonal
+    Iq = eye(Nc, format='csr', dtype=bool)
 
-    S = bmat([[Ac, I],
-              [I,  I]], format="csc")
-    S.data[:] = 1
+    # assemble block matrix
+    S = bmat([[Ac.astype(bool), B],
+              [B,               Iq]], format='csr')
     return S

chromosim/profiles/__init__.py

@@ -0,0 +1,42 @@
+from __future__ import annotations
+import numpy as np
+
+def salt_profile_clamped_linear(t, par):
+    """
+    Clamped linear ramp:
+        I(t) = I_load + (I_elute - I_load) * clamp((t - grad_start)/(grad_end - grad_start), 0, 1)
+    Accepts scalar or array t. Returns same shape as t.
+    """
+    tarr = np.atleast_1d(np.asarray(t, dtype=float))
+    frac = (tarr - par.grad_start) / (par.grad_end - par.grad_start)
+    frac = np.clip(frac, 0.0, 1.0)
+    I = par.I_load + (par.I_elute - par.I_load) * frac
+    return I if np.ndim(t) else float(I.item())
+
+def feed_profile_by_load_volume(t, par):
+    """
+    Physical inlet concentration vector for Ns species (length Ns).
+    During the 'load' window, Cin = Cfeed_total * frac; otherwise zeros.
+    The load window is computed as:
+        load_duration_s = 60 * load_mL / flow_mL_min
+    If those fields are missing, defaults to 300 s.
+    Works with scalar or array t. For scalar t returns 1×Ns array.
+    """
+    Ns = par.Nsp
+    tarr = np.atleast_1d(np.asarray(t, dtype=float))
+    Cin = np.zeros((tarr.size, Ns), dtype=float)
+
+    # derive load duration if not given
+    if not hasattr(par, "load_duration_s"):
+        if hasattr(par, "load_mL") and hasattr(par, "flow_mL_min"):
+            par.load_duration_s = 60.0 * float(par.load_mL) / float(par.flow_mL_min)
+        else:
+            par.load_duration_s = 300.0  # sensible default (5 min)
+
+    mask_load = tarr < par.load_duration_s
+    if np.any(mask_load):
+        Cin[mask_load, :] = (par.Cfeed_total * np.asarray(par.frac)).reshape(1, Ns)
+
+    if np.ndim(t) == 0:
+        return Cin[0, :]
+    return Cin

chromosim/profiles/salt_and_feed.py

@@ -0,0 +1,15 @@
+from __future__ import annotations
+import numpy as np
+
+def unpack_state(Y: np.ndarray, par):
+    """
+    Reshape solution (Nt × (2*N*Ns)) into:
+        C : Nt × N × Ns   (mobile phase, physical or scaled depending on caller)
+        q : Nt × N × Ns   (bound phase)
+    """
+    Nt, total = Y.shape
+    N, Ns = par.N, par.Nsp
+    NN = N * Ns
+    C = Y[:, :NN].reshape(Nt, N, Ns)
+    q = Y[:, NN:].reshape(Nt, N, Ns)
+    return C, q