Generic profiles (#2)

* Python: add finite-volume derivative matrix builder (build_matrices)

* Scaffold ChromoSim Python package: matrices, sparsity, multi-species RHS

* Generic salt/feed profiles, Jacobian sparsity, unpack helper

Author: mironovb

chromosim/models/__init__.py

@@ -0,0 +1,80 @@
+from __future__ import annotations
+import numpy as np
+
+def _get(par, key):
+    try:
+        return getattr(par, key)
+    except AttributeError:
+        return par[key]
+
+def ev_column_ode_multi(t: float, y: np.ndarray, par) -> np.ndarray:
+    """
+    Multi-species 1D convection–dispersion + competitive Langmuir adsorption
+    with zeta- and ionic-strength–dependent kinetics, integrated in *scaled* vars.
+
+    y = [ c(:) ; qhat(:) ], with c,Q scaled by Cscale,Qscale respectively.
+    """
+    N   = int(_get(par, "N"))
+    Ns  = int(_get(par, "Nsp"))
+    dz  = float(_get(par, "dz"))
+    u   = float(_get(par, "u"))
+    Dax = float(_get(par, "Dax"))
+    eps = float(_get(par, "eps"))
+
+    Cscale = float(_get(par, "Cscale"))
+    Qscale = float(_get(par, "Qscale"))
+    qmax   = float(_get(par, "qmax"))
+
+    zeta_mV = np.asarray(_get(par, "zeta_mV"), dtype=float).reshape(-1)  # (Ns,)
+    NN = N * Ns
+
+    c    = y[0:NN].reshape(N, Ns)            # dimensionless
+    qhat = y[NN:2*NN].reshape(N, Ns)         # dimensionless
+
+    C = c    * Cscale                         # physical
+    Q = qhat * Qscale
+
+    I  = float(_get(par, "salt_profile")(t))  # ionic strength (scalar)
+    fI = float(_get(par, "fI")(I))
+
+    K0    = float(_get(par, "K0"))
+    gamma = float(_get(par, "gamma"))
+    Ki_vec = K0 * np.exp(gamma * np.abs(zeta_mV) * fI)  # (Ns,)
+
+    # desorption (constant or function of I)
+    kd_fun = getattr(par, "kd_fun", None) if hasattr(par, "kd_fun") else _get(par, "kd0")
+    if callable(kd_fun):
+        kd_vec = np.asarray(kd_fun(I), dtype=float).reshape(-1)
+        if kd_vec.size == 1:
+            kd_vec = np.repeat(kd_vec, Ns)
+    else:
+        kd_vec = np.ones(Ns) * float(_get(par, "kd0"))
+
+    ka_vec = Ki_vec * kd_vec  # (Ns,)
+
+    ka = np.broadcast_to(ka_vec, (N, Ns))
+    kd = np.broadcast_to(kd_vec, (N, Ns))
+
+    Qsum      = np.sum(Q, axis=1)             # (N,)
+    Qfree     = np.maximum(qmax - Qsum, 0.0)  # (N,)
+    Qfree_mat = Qfree[:, None]                # (N,1) -> broadcast
+
+    dQdt    = ka * C * Qfree_mat - kd * Q
+    dqhatdt = dQdt / Qscale
+
+    D1 = _get(par, "D1")
+    D2 = _get(par, "D2")
+    dcdz  = D1.dot(c)
+    d2cdz = D2.dot(c)
+
+    beta = ((1.0 - eps) / eps) * (Qscale / Cscale)
+    dcdt = -u * dcdz + Dax * d2cdz - beta * dqhatdt
+
+    # inlet boundary: Dirichlet via ghost/upwind (scaled)
+    Cin_phys = np.asarray(_get(par, "feed_profile")(t), dtype=float).reshape(-1)
+    if Cin_phys.size != Ns:
+        raise ValueError(f"feed_profile(t) must return length {Ns}, got {Cin_phys.size}")
+    cin = Cin_phys / Cscale
+    dcdt[0, :] += (-u/dz) * (c[0, :] - cin) + (Dax/dz**2) * (cin - c[0, :])
+
+    return np.concatenate([dcdt.reshape(-1), dqhatdt.reshape(-1)])

chromosim/models/ev_column_ode_multi.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.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/models/params_ev_aex_single.py

@@ -0,0 +1,28 @@
+from __future__ import annotations
+import numpy as np
+from scipy.sparse import spdiags, kron, eye, bmat, csr_matrix
+
+def sparsity_pattern_multi(N: int, Ns: int):
+    """
+    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 = 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
+
+    Nc = N * Ns
+    B  = eye(Nc, format='csr', dtype=bool)    # diagonal
+    Iq = eye(Nc, format='csr', dtype=bool)
+
+    # assemble block matrix
+    S = bmat([[Ac.astype(bool), B],
+              [B,               Iq]], format='csr')
+    return S

chromosim/numerics/sparsity.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/__init__.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