aquakin.BiofilmReactor#

class aquakin.BiofilmReactor(model, conditions, *, n_layers, thickness, area_per_volume, diffusivity, boundary_layer, boundary_diffusivity=None, soluble_mask=None, fixed_mask=None, max_density=None, packing_fraction=0.8, k_att=0.0, attach_mask=None, k_det=0.0, detach_mask=None, clamp_bulk=False, feed=None, dilution_rate=0.0, biofilm_reactions=None, rtol=1e-06, atol=None, integrator=IntegratorConfig(order=3, factormax=3.0, colored_jacobian='auto', dtmax=None, max_steps=100000, solver=None), diff=DifferentiationConfig(mode='reverse', method='stable', check_finite=True, adjoint_max_steps=100000, adjoint_low_memory=False))[source]#

Bases: GradientCheckMixin

Stateless 1-D layered biofilm reactor (diffusion–reaction over depth).

Parameters:

  • model (CompiledModel) – Compiled reaction model, run in every compartment.

  • conditions (SpatialConditions) – Condition fields. The same conditions are used in every compartment (location loc_idx=0).

  • n_layers (int) – Number of biofilm layers between the bulk and the wall.

  • thickness (float) – Total biofilm thickness L_f (metres). Layer thickness is thickness / n_layers.

  • area_per_volume (float) – Biofilm-area-to-bulk-volume ratio A_V (m^2 biofilm per m^3 bulk, i.e. 1/m). Sets how strongly the (small) biofilm exchange moves the (large) bulk pool.

  • diffusivity (float or jnp.ndarray) – Effective diffusivity D_eff of solubles inside the biofilm (m^2/day). A scalar applies to every soluble; an array of shape (n_species,) gives a per-species value (entries for particulates are ignored). Particulates never diffuse.

  • boundary_layer (float) – External boundary-layer thickness L_bl (metres) across which the bulk exchanges with the surface layer. The per-species mass-transfer coefficient is k_L = D_bl / L_bl.

  • boundary_diffusivity (float or jnp.ndarray, optional) – Diffusivity used in the external boundary layer (m^2/day). The boundary layer is liquid (outside the biofilm matrix), so it should carry the free-water diffusivity D_w, not the density-reduced in-biofilm D_eff. Same shape rules as diffusivity. None (default) reuses diffusivity for backward compatibility, which understates the bulk<->film transfer by the biofilm reduction factor; pass the water values to model the boundary layer correctly.

  • soluble_mask (jnp.ndarray, optional) – Boolean (n_species,) mask: True for solubles (which diffuse), False for particulates (which do not). Defaults to the S/X name heuristic. This controls diffusion only — whether a species is held fixed is a separate question governed by fixed_mask.

  • fixed_mask (jnp.ndarray, optional) – Boolean (n_species,) mask: True for species held fixed (net rate zeroed everywhere — the “mature biofilm” sustained source/sink), False for species that evolve. Defaults to ~soluble_mask (every particulate held fixed), which is correct for inert particulates (biomass, inert solids) but wrong for reactive particulates that do not diffuse yet must still react — e.g. elemental sulfur or precipitated FeS, whose inventory genuinely drains and fills. For such species pass a fixed_mask that holds only the inert biomass/solids fixed and lets the reactive particulates evolve; otherwise a held-fixed reactive particulate becomes a non-depleting source/sink and silently breaks mass balance.

  • biofilm_reactions (list of str or jnp.ndarray, optional) – Which reactions are biofilm processes — they run in the biofilm layers only, never in the well-mixed bulk. Given as a list of reaction names or a boolean (n_reactions,) mask. The remaining reactions (bulk-suspended and bulk-chemical) run in the bulk only. None (default) runs every reaction in every compartment, appropriate for a single-phase model. For a WATS-style model the biofilm reactions are those carrying the A_V area factor; bulk reactions carry the suspended biomass [X_BH] or are abiotic. This separation, not a zeroed biomass state, is what keeps the two phases from leaking into each other.

  • rtol (float) – Relative solver tolerance, applied across the whole layered state.

  • atol (float or jnp.ndarray, optional) – Absolute solver tolerance. None (default) uses the per-component default_atol() noise floor scaled off the model’s reference concentrations, tiled across every compartment (so g/m3 ASM/ADM states get a sensible floor instead of a fixed scalar ~9 orders too tight). A scalar broadcasts to the whole state; a (n_species,) array is tiled across the n_layers+1 compartments.

  • integrator (IntegratorConfig, optional) – Integrator / step-size configuration (ESDIRK order, factormax, dtmax, max_steps, an explicit solver); see BatchReactor. A stiff per-layer-biomass model with a tight dtmax can exceed the default max_steps; raise it if the solve raises a max-steps error.

  • diff (DifferentiationConfig, optional) – Autodiff configuration (mode, method); see BatchReactor.

  • packing_fraction (float)

  • k_att (float)

  • attach_mask (jnp.ndarray | None)

  • k_det (float)

  • detach_mask (jnp.ndarray | None)

  • clamp_bulk (bool)

  • dilution_rate (float)

Notes

The state is a (n_layers + 1, n_species) array: row 0 is the bulk, rows 1..``n_layers`` are the biofilm layers from surface to wall. The wall is a no-flux boundary. Species in fixed_mask have their net rate zeroed everywhere (held fixed); all others evolve. Diffusion is governed separately by soluble_mask.

Examples

Sulfur biofilm models carry reactive non-diffusing particulates — elemental sulfur X_S0 and precipitated X_FeS — whose inventory genuinely drains and fills. The default fixed_mask would freeze them and silently break mass balance (the reactor warns when it would). Build a mask that holds only the inert solids fixed and lets every reactive pool evolve:

>>> net = aquakin.load_model("wats_sewer_khalil_paper_balanced_biofilm_multispecies")
>>> inert = {"X_I"}   # the only genuinely inert, non-depleting solid
>>> fixed_mask = jnp.array([s in inert for s in net.species])
>>> reactor = aquakin.BiofilmReactor(
...     net, conditions, n_layers=6, thickness=8e-4, area_per_volume=50.0,
...     diffusivity=1e-4, boundary_layer=1e-4, fixed_mask=fixed_mask)

Everything not named in inert — the heterotrophs and functional-group biomass, the stored-substrate reservoirs X_S1/X_S2, and the reactive sulfur pools X_S0/X_FeS — then evolves and conserves mass. (The areal *_biofilm variant, by contrast, deliberately freezes its biomass as a sustained “mature biofilm” reservoir, so its mask holds the biomass fixed too; the rule is always “freeze only what is genuinely non-depleting”.)

__init__(model, conditions, *, n_layers, thickness, area_per_volume, diffusivity, boundary_layer, boundary_diffusivity=None, soluble_mask=None, fixed_mask=None, max_density=None, packing_fraction=0.8, k_att=0.0, attach_mask=None, k_det=0.0, detach_mask=None, clamp_bulk=False, feed=None, dilution_rate=0.0, biofilm_reactions=None, rtol=1e-06, atol=None, integrator=IntegratorConfig(order=3, factormax=3.0, colored_jacobian='auto', dtmax=None, max_steps=100000, solver=None), diff=DifferentiationConfig(mode='reverse', method='stable', check_finite=True, adjoint_max_steps=100000, adjoint_low_memory=False))[source]#
Parameters:

  • model (CompiledModel)

  • conditions (SpatialConditions)

  • n_layers (int)

  • thickness (float)

  • area_per_volume (float)

  • boundary_layer (float)

  • soluble_mask (Array | None)

  • fixed_mask (Array | None)

  • packing_fraction (float)

  • k_att (float)

  • attach_mask (Array | None)

  • k_det (float)

  • detach_mask (Array | None)

  • clamp_bulk (bool)

  • dilution_rate (float)

  • rtol (float)

  • integrator (IntegratorConfig)

  • diff (DifferentiationConfig)

Return type:

None

Methods

__init__(model, conditions, *, n_layers, ...)

check_gradient_finite(grad_value, *[, what])

Raise a friendly error if a reverse-mode gradient is non-finite.

solve(C0, t_span[, t_eval, params, ...])

Integrate the layered biofilm over a time span.

solve_sensitivity(C0, t_span[, t_eval, ...])

Solve and return the forward sensitivity of the bulk dC/dtheta.

steady_state(C0, params, *[, conditions, ...])

Solve for the steady-state profile by pseudo-transient continuation.
solve(C0, t_span, t_eval=None, *, params=None, conditions=None, time_unit=None)[source]#

Integrate the layered biofilm over a time span.

Parameters:

  • C0 (jnp.ndarray) – Initial state. Either (n_species,) — the same composition in the bulk and every layer — or (n_layers + 1, n_species) to set the bulk and each layer explicitly (row 0 bulk, rows 1.. surface to wall). The latter sets the stratified particulate (biomass) profile.

  • t_span (tuple of float) – (t_start, t_end) integration interval, in the model’s time unit unless time_unit is given. The required second positional argument.

  • t_eval (jnp.ndarray, optional) – Times at which to record the solution. If None only the endpoint.

  • params (jnp.ndarray, optional, keyword-only) – Rate constant vector, shape (n_params,). Defaults to model.default_parameters(). Keyword-only so a positional t_span can never land in it.

  • conditions (SpatialConditions, optional) – Override the reactor conditions for this call.

  • time_unit (str, optional) – The time unit t_span / t_eval are in ("s"/"min"/ "h"/"d"); see BatchReactor.solve(). Default None uses the model’s native unit.

Return type:

BiofilmSolution

solve_sensitivity(C0, t_span, t_eval=None, *, params, sens_params, conditions=None, sens_rtol=None, sens_atol=None, param_scale=None, shared_factor=None)[source]#

Solve and return the forward sensitivity of the bulk dC/dtheta.

Integrates the augmented [y; S] system over the full layered state with adaptive control over both, so the sensitivity is exact and finite without the dtmax cap that ordinary AD through this stiff diffusion–reaction solve needs (see aquakin.integrate.forward_sensitivity). This is the canonical use case: the biofilm models are stiff enough that capping dtmax for a reverse-mode gradient is an ~10x penalty, and this removes it.

Parameters:

  • C0 (Array) – As for solve() (C0 may be (n_species,) or (n_layers + 1, n_species)).

  • t_span (tuple[float, float]) – As for solve() (C0 may be (n_species,) or (n_layers + 1, n_species)).

  • t_eval (Array | None) – As for solve() (C0 may be (n_species,) or (n_layers + 1, n_species)).

  • conditions (SpatialConditions | None) – As for solve() (C0 may be (n_species,) or (n_layers + 1, n_species)).

  • params (jnp.ndarray) – Full parameter vector; keyword-only, matching solve().

  • sens_params (list of str or int) – Namespaced parameter names (or integer indices into params).

  • sens_rtol (float | None) – Sensitivity error-control tolerances (CVODES defaults; see augmented_forward_sensitivity()).

  • sens_atol – Sensitivity error-control tolerances (CVODES defaults; see augmented_forward_sensitivity()).

  • param_scale – Sensitivity error-control tolerances (CVODES defaults; see augmented_forward_sensitivity()).

  • shared_factor (bool, optional) – CVODES simultaneous-corrector linear solve. None (default) auto-selects True for more than one sensitivity parameter – the regime where factorising the shared diagonal block once is markedly cheaper than the dense augmented solve – else False.

Returns:

  • sol (BiofilmSolution) – The usual solution (bulk C and full profile).

  • S (jnp.ndarray) – Sensitivity of the bulk (measurable) concentration, dC_bulk/dtheta, shape (n_t, n_species, n_sens_params) – aligned with sol.C.

Return type:

tuple[BiofilmSolution, Array]

steady_state(C0, params, *, conditions=None, warmup=20.0, rtol=1e-06, atol=1e-08, newton_steps=200)[source]#

Solve for the steady-state profile by pseudo-transient continuation.

Instead of integrating to steady state (slow for a maturing biofilm), this finds y* such that f(y*, params) = 0 directly, by pseudo-transient continuation (PTC) – damped-Newton steps that ramp from a stable time-stepping move to a full Newton step as the residual falls, robust on the stiff/slow biofilm where a plain Newton / Levenberg–Marquardt root-find stalls (see aquakin.plant.steady.solve_steady_state()). A short forward integration to warmup seeds the iteration (the seed is detached from the gradient). The result is differentiable w.r.t. params via the implicit function theorem, so it composes with calibrate() – the intended use is the continuous-feed maturation (feed=..., dilution_rate>0) whose steady bulk is the predicted effluent and whose biofilm profile is a downstream batch IC.

Parameters:

  • C0 (jnp.ndarray) – Initial guess, (n_species,) (broadcast) or (n_layers+1, n_species).

  • params (jnp.ndarray) – Rate constant vector.

  • conditions (SpatialConditions, optional) – Override the reactor conditions for this call.

  • warmup (float) – Forward-integration time used to seed the iteration. 0 uses C0 directly.

  • rtol (float) – Convergence tolerance on the scaled steady-state residual max_i |f_i| / max(|y_i|, 1).

  • atol (float) – Retained for backward compatibility; unused (PTC converges on the relative residual rtol).

  • newton_steps (int) – Maximum PTC iterations.

Returns:

With t = [inf] and a single time slice holding the steady profile.

Return type:

BiofilmSolution

Notes

Not compatible with clamp_bulk or held-fixed species (their RHS rows are identically zero, making the residual Jacobian singular). Use the continuous feed to drive the bulk instead.