From 4ccef4f3a95f02b810fa18eb76b75e0e9ff7a80e Mon Sep 17 00:00:00 2001 From: Mike Arpaia Date: Thu, 24 Sep 2026 15:48:31 -0600 Subject: [PATCH] Specify conservative cell occupancy transport with executable reference cases --- .../architecture/0025-cell-occupied-volume.md | 104 ++++++ docs/architecture/README.md | 1 + docs/microfluidics.md | 4 + .../src/microsimulator/occupancy_reference.py | 314 ++++++++++++++++++ python/tests/test_occupancy_reference.py | 178 ++++++++++ 5 files changed, 601 insertions(+) create mode 100644 docs/architecture/0025-cell-occupied-volume.md create mode 100644 python/src/microsimulator/occupancy_reference.py create mode 100644 python/tests/test_occupancy_reference.py diff --git a/docs/architecture/0025-cell-occupied-volume.md b/docs/architecture/0025-cell-occupied-volume.md new file mode 100644 index 0000000..e55fe66 --- /dev/null +++ b/docs/architecture/0025-cell-occupied-volume.md @@ -0,0 +1,104 @@ +# ADR 0025: conservative extracellular storage with coarse geometric porosity + +- Status: selected numerical design; executable CPU reference only +- Date: 2026-09-24 +- Scope: issue #13. Native implementation and enabling this model in simulations are separate work. + +## Decision and current behavior + +Select **coarse geometric porosity** as an opt-in transport model. Estimate the union of cell capsules inside each non-wall voxel, store solute amount per voxel, and use a declared porosity closure for face conductance. This approximates cell exclusion without claiming to resolve sub-voxel fluid passages, membrane boundary layers, displacement flow, or hydrodynamic forces. It is preferable here to treating smoothed biochemical biomass density as an exact solid fraction. + +Current production transport continues to use full non-wall voxel volume: `SignalGridSpec.voxel_volume()` is `hx*hy*hz`, and `cpu_coupled.cpp` divides scattered cell amount rates by that volume. `colony_volume_fraction()` deposits conserved biochemical biomass for empirical resistance; it may exceed one and is not this occupancy representation. No production defaults or checkpoint formats change in this contribution. + +The reference is `python/src/microsimulator/occupancy_reference.py`. It uses float64 midpoint quadrature and a dense backward-Euler solve, intentionally unsuitable for large simulations. Its interfaces carry explicit amounts, volumes, faces, and ledgers so later backends can be compared without inheriting their implementation. + +## State, geometry, and units + +For voxel i, geometric volume `V_i = hx*hy*hz`, accessible fraction `epsilon_i`, accessible storage `W_i = epsilon_i V_i`, and extracellular concentration `c_i`, define the authoritative solute amount `N_i = W_i c_i`. Concentration is amount per accessible fluid volume, including in partially occupied voxels. Length has units L, time T, concentration A/L³, and amount A. + +Occupancy depends only on the current cell centers, normalized directions, nonnegative centerline lengths, positive radii, lattice centers/spacing, and the existing binary transport-wall mask. It does not depend on species, growth-rate attributes, cell type, stationary attachment, or the biomass-resistance averaging radius. All cells, including mobile ones, exclude storage. Mechanical wall primitives do not define a second transport mask: only the declared voxel obstacle mask clips storage. + +Each voxel uses m³ midpoint samples on its physical box centered at `origin + index*spacing`. A sample is occupied when its distance to any capsule centerline segment is at most the corresponding radius. Count the union, never summed capsule fractions, so overlaps cannot make epsilon negative. Outside-domain capsule pieces are ignored. A wall voxel has epsilon zero regardless of cells. Default reference m is 8; m is a reproducible numerical parameter requiring convergence checks. Even a degenerate lattice axis retains its physical voxel thickness for occupancy and storage; it does not turn capsules into disks. + +Biochemical biomass remains `B = pi*r²*(length + 2*r)`. Geometric capsule volume is `V_geom = pi*r²*length + 4*pi*r³/3`. Native division preserves B but, for equal radii, reduces total capsule volume by `2*pi*r³/3`. This closure deliberately exposes the resulting extra fluid storage. It represents coarse division geometry, not physical septum formation or calibrated cell volume. Never transfer extracellular solute into daughters merely because their geometric occupancy changed. + +## Conservative transport balance + +For an internal oriented face i→j, use harmonic closure `a_f = 2 epsilon_i epsilon_j / (epsilon_i + epsilon_j)` when both voxels are accessible, otherwise zero. This is a coarse permeability/aperture approximation, not geometric face intersection. Let `K_f = D a_f A_f / d_f` and `Q_f = a_f A_f u_f`. D has units L²/T, K and Q have units L³/T, and u is intrinsic accessible-fluid velocity in L/T. Positive Q flows from i to j. The outward amount rate is + +```text +F_ij = K_f (c_i - c_j) + max(Q_f, 0) c_i + min(Q_f, 0) c_j +``` + +The identical face value enters the neighbor with the opposite sign. With amount source S_i, accessible-fluid reaction source b_i, and first-order loss lambda_i: + +```text +dN_i/dt = -sum_faces F_ij + S_i + W_i b_i - lambda_i N_i +``` + +After the geometric remap described below, freeze W and face coefficients over one transport step. The first native implementation should offer backward Euler: + +```text +(W + dt L) c_new = N_remapped + dt (S + W b + reservoir_inflow) +N_new = W c_new +``` + +L has diffusion conductances, upwind outgoing fluxes, and `lambda_i W_i` on its diagonal, with negative incoming neighbor coefficients. Zero-storage rows are isolated identity rows with zero right-hand side. The reference reports before/after total amount plus separately integrated source, reaction, and boundary amounts; internal faces cancel. Do not infer conservation from changes in concentration or from an unweighted sum of concentrations. + +No-flux boundaries have K=Q=0. Periodic pairs contribute one shared internal face. Fixed reservoirs retain the existing **exterior lattice-center** convention: distance d is one spacing, not a half spacing. For this closure, use exterior epsilon=1 and the same harmonic rule. Reservoir exchange is `K(c_res-c_i) - Q*c_upwind` with Q positive outward, and is included in the boundary ledger. In a singleton axis production transport has no face operator, matching the existing degenerate-axis rule. Affine b is concentration/time per accessible fluid; a physical source specified per whole voxel instead becomes an explicit amount rate, never an implicit second volume conversion. + +## Velocity convention and flow coupling + +The current face field is a velocity used directly by transport and also sampled by cell drift. With the existing binary wall mask and epsilon=1 on fluid sites, intrinsic velocity and whole-open-face volumetric velocity coincide. Partial porosity introduces a distinction that cannot be inferred from old arrays. + +In the new opt-in mode, retain **intrinsic velocity** in `SignalGridVelocityField` and derive Q by multiplying `a_f*A_f` once. If a flow solver produces integrated flux Q, convert to intrinsic u by dividing by `a_f*A_f` once on open faces and require Q=0 on closed faces. Never multiply an already aperture-weighted flux by epsilon again. Constant-advection inputs use the same explicit convention. + +Existing flow solutions generally satisfy continuity for their current binary fluid geometry, not for these new weighted faces. They must not be advertised as occupancy-consistent flow without a weighted projection/re-solve. For fixed occupancy, validate `sum Q_out = 0` in interior voxels. Moving occupancy would require `dW/dt + sum Q_out = 0` for incompressible displaced fluid; this initial model instead uses the explicit conservative remap below. That remap is a solute bookkeeping closure and does not solve fluid displacement. Prescribed velocity experiments must declare this limitation; coupling a resolved displacement-flow model is subsequent work. Transport remains amount-conservative even for a prescribed divergent field, but uniform concentration need not remain uniform. + +## Occupancy changes, closed storage, and transaction order + +Fractions smaller than `epsilon_cutoff = 1e-8` are treated as zero, including for face closure. The cutoff is part of the model configuration and checkpoint state, not a hidden denominator floor. Keep N fixed in each voxel whose new W remains positive, so concentration changes to N/W. Newly accessible voxels start with their existing amount, normally zero; removal does not invent extracellular solute. + +For voxels changing to zero storage, redistribute their entire amount over newly accessible recipient voxels in the same face-connected component of the **union of old and new accessible voxels**, weighted by recipients' new W. Connectivity uses regular voxel neighbors including declared periodic pairs and excludes persistent wall/closed voxels. Use sorted voxel order for deterministic accumulation. Existing recipient amounts are retained. If a component closes completely while holding any positive amount, reject the entire geometry/transport transaction. If it contains exactly zero amount, closure is valid. No epsilon floor, disappearing amount, cross-wall transfer, or silent clipping is allowed. A near-zero positive W can still create high concentrations; finite/positivity and solver checks may reject the step rather than alter its mass balance. + +This remap is intentionally nonlocal within its transition component and can instantaneously mix expelled solute. It does not reconstruct a membrane trajectory or predict a swept-volume velocity. Its domain, cutoff, and redistribution rule must therefore be stated in scientific use. A donor is emptied once even when several cells overlap it. Reject invalid input before changing simulation state. + +The proposed opt-in controller/native stage order is: + +1. Snapshot complete state and authoritative N at the last committed geometry. Apply validated regulation, divisions, removals, and their geometry callbacks. Recompute occupancy and conservatively remap if geometry changed. +2. Advance biological growth and intracellular dilution using B. Recompute occupancy for post-growth geometry and remap N again. Construct one set of accessible exchange weights, sample this remapped pre-transport concentration, evaluate biology, then solve transport/reactions and scatter cell amount rates using the same weights. This explicitly changes the sampling point from the legacy stage; it needs a separately named opt-in contract and native split-stage work. +3. Apply configured flow drift and mechanical relaxation, then recompute/remap once at final geometry. Any direct geometry edit or removal outside the controller must invoke the same barrier before the next transport/export/checkpoint operation. No geometry change may leave old W paired with new cells. +4. Commit geometry, biomass/species, N, W, time, and all ledgers together only after every stage validates. A failure restores the complete transaction, including controller/RNG state. Zero-time topology edits still require remapping. Unchanged geometry reuses occupancy without resampling. + +This ordering follows the existing regulation/topology → biology → drift/mechanics outline while making its new geometry barriers explicit. It is not implemented by the reference module. Production work must add atomic staging; installing only a new denominator inside `cpu_coupled.cpp` would be incorrect. + +## Cell sampling and scattering + +Choose a declared physical support that reaches accessible extracellular sites, restricted to the connected fluid component containing a deterministic nearest-accessible anchor. Resolve equal-distance anchors by flattened index and explicitly validate the maximum search radius; no cross-wall search. For nonnegative geometric kernel values phi_i in that support, use `w_i = phi_i W_i / sum(phi_j W_j)`. Sampling is `sum(w_i c_i)`; a cellular extracellular amount rate J scatters `w_i J`. The concentration-rate conversion is `(w_i J)/W_i` exactly once, for accessible sites only. The same weights give a partition of unity and the sampling/scattering adjoint relation. The reference `exchange_weights` receives an already connected support; topology construction remains native implementation work. + +With epsilon=1 and equal voxel volumes this reduces to the current normalized trilinear weights and J/V conversion. With resolved excluded cell centers, the old eight center-neighbor stencil may have no accessible sites; expanding or surface-based physical support is required, with a declared radius and refinement study. Reject zero accessible support. Uptake must be limited by an explicitly conservative coupled solve or reject an unaffordable step; never clamp negative extracellular concentrations. Opposite intracellular amount uses B, not geometric occupancy. Washout of a cell exports its intracellular amount through a separate biological ledger; it does not remove a voxel's extracellular solute. + +## Executable reference cases and tolerances + +Run `uv run python -m pytest python/tests/test_occupancy_reference.py -v`. The cases are deterministic CPU tests: + +| Case | Amount evidence | +| --- | --- | +| Empty grid | Float64 reference equals existing native CPU backward Euler for `[0,1,0]`, within rtol 2e-6/atol 2e-7; total amount remains 1. Empty geometry gives epsilon=1. | +| Partially occupied closed domain and unequal storage | W=[0.5,1.5], initial c=[8,0], total N=4; diffusion tends to equal c=2 with unequal amounts [1,3]. Every-step ledger residual is below 1e-12. | +| Changing occupancy, closure, reopening | N=[2,3,1], W changes [1,1,1]→[0,0.25,0.75]; remap gives [0,3.5,2.5], total 6. Reopening keeps new-site amount zero. Entire-component closure rejects without mutation; a persistent wall prevents redistribution to an unrelated recipient. | +| Division/removal, overlap, wall-adjacent geometry | One native-shaped parent becomes two daughters with unchanged B and smaller geometric volume. Union quadrature does not double-count duplicate cells, wall voxels remain inaccessible, and remapping through division/removal preserves total N. | +| Boundaries, reactions, exchange and advection | Unequal storage, an internal advective face, inflow/outflow reservoirs, decay and cellular source all contribute to the signed ledger; residual below 1e-12. A face test verifies aperture enters Q once. | +| Refinement | Backward-Euler timestep error approximately halves for 10/20/40 steps; the empty-limit centered operator error approximately quarters for 10/20/40 voxels. Sphere quadrature at m=8/16/32 improves the finest volume error to below 2%. | + +The 1e-12 reference balance tolerance applies to these order-one float64 cases. General checks scale by `max(1, |N_before|, |N_after|, sum(abs(external terms)))`; scientific units must be normalized explicitly. Initial native float32 gates: per-step relative ledger residual <=5e-6, 1000-step closed-case drift <=5e-5, and concentrations versus float64 reference within rtol 2e-4/atol 2e-6 in normalized test units. These are acceptance targets to measure, not verified GPU results. + +Subsequent native validation must run h, h/2, h/4 at fixed physical cell/device dimensions and exchange support, with m, 2m, 4m independently; also dt, dt/2, dt/4. Track occupied volume, total amount, concentration L1/Linf errors, boundary/reaction/cell ledgers, cutoff crossings, and solver residuals. Require convergence of the scientific observable, not only conservation. The midpoint geometry estimate can oscillate across resolutions; never assert a universal smooth-interface order from one placement. Test rotated/translated rods, overlapping capsules, nearly blocked passages, all-zero storage, wall contacts, division/removal, and restart at a geometry barrier. Vary epsilon_cutoff by factors of ten. Native backends must implement their own kernels and pass the same cases without CPU fallback. + +## Compatibility and implementation work + +Old checkpoints and all ordinary runs retain occupancy-disabled full-voxel semantics. Do not reinterpret old concentration arrays as fluid-volume concentrations. Enabling occupancy on an existing state is an explicit conversion: preserve legacy amount `N=V*c`, compute W, then apply closure remapping and derive c. Such a conversion can reject a sealed component and must be recorded in provenance. + +A future checkpoint version must record model kind/version, lattice and wall geometry, quadrature resolution, cutoff, exchange-support/anchor rule, remap rule, velocity convention, authoritative per-species N and committed W/geometry revision, plus any required solver history. N and W participate in integrity checks. On restore, authenticate before migration, validate N>=0 and zero amount at W=0, verify geometry/occupancy consistency, and resume at a committed barrier without repeating remapping. Record precision/backend provenance; recomputing occupancy with another quadrature algorithm may change results and requires an explicit conversion. The current reference adds no fields to checkpoints. + +Required follow-up contributions are (1) native state/configuration and checkpoint conversion, (2) conservative geometry rasterization/transition connectivity, (3) CPU weighted-storage operator and atomic coupled staging, (4) independent Metal and CUDA kernels/reductions/solvers, (5) accessible cell exchange and uptake budgets, (6) weighted-flow interface/projection with explicit drift semantics, and (7) end-to-end geometry/removal/restart and refinement validation. Existing flow resistance may coexist as a calibrated closure; it must not be relabeled as geometric exclusion or counted a second time in transport porosity. diff --git a/docs/architecture/README.md b/docs/architecture/README.md index a823b8c..643f936 100644 --- a/docs/architecture/README.md +++ b/docs/architecture/README.md @@ -24,6 +24,7 @@ The [numerical contract](numerical-contract.md) is the best starting point for w - [Persistent fixed rod cells](0009-fixed-cells.md) - [Biomass, growth, division, and uptake](0024-biomass-accounting.md) - [Typed species rate plans](0003-species-rates.md) +- [Cell-occupied extracellular volume: design and CPU reference](0025-cell-occupied-volume.md) - [Grid signaling and cell coupling](0006-grid-signaling.md) - [Crank-Nicolson signal transport](0008-crank-nicolson-signals.md) - [Neighbor diffusion](0010-neighbor-diffusion.md) diff --git a/docs/microfluidics.md b/docs/microfluidics.md index 833dc95..286f1b8 100644 --- a/docs/microfluidics.md +++ b/docs/microfluidics.md @@ -31,6 +31,10 @@ Attached biomass can change flow resistance through a conservatively smoothed de Free-cell motion uses the local velocity and a finite-aspect Jeffery orientation approximation, followed by contact relaxation. This kinematic coupling approximates rods as equivalent spheroids for rotation. Cell-scale hydrodynamic forces, lubrication, and predictive adhesion or detachment are outside its scope. The [flow-drift design](architecture/0021-flow-drift.md) specifies the approximation and integration limits. +## Cell-occupied extracellular volume + +Current native transport stores concentration per full non-wall voxel; cells do not yet exclude extracellular storage. The smoothed biochemical biomass density used for flow resistance is neither bounded geometric occupancy nor a resolved fluid fraction. [ADR 0025](architecture/0025-cell-occupied-volume.md) selects an opt-in coarse geometric-porosity model and provides executable CPU reference cases, including conservative amount remapping as geometry changes. It is a numerical design, not an enabled production feature. Its harmonic face closure and component-level redistribution do not resolve fluid passages around individual cells, membrane transport layers, or displacement flow; geometric, spatial, and timestep refinement remain required. + ## Interpreting results The [analytic flow benchmarks](tutorials/flow-solvers.md#numerical-evidence) test profile convergence, flux routing, and agreement between the solvers in a shared thin-gap regime. The [controlled nutrient study](tutorials/nutrient-validation.md) measures spatial growth, nutrient balance, and sensitivity to grid spacing, timestep, and flow-refresh interval. It isolates attached-population growth and transport; the interactive tutorials exercise division, mechanics, and washout separately. diff --git a/python/src/microsimulator/occupancy_reference.py b/python/src/microsimulator/occupancy_reference.py new file mode 100644 index 0000000..7b0544d --- /dev/null +++ b/python/src/microsimulator/occupancy_reference.py @@ -0,0 +1,314 @@ +"""Float64 design reference for ADR 0025; not a native transport implementation. + +Amounts are authoritative. Coarse geometric porosity is distinct from the +biochemical biomass density used by the flow-resistance closure. +""" + +from __future__ import annotations + +import math +from collections.abc import Sequence +from dataclasses import dataclass +from typing import cast + +import numpy as np +from numpy.typing import NDArray + +Array = NDArray[np.float64] +EPSILON_CUTOFF = 1.0e-8 + + +@dataclass(frozen=True) +class Capsule: + center: tuple[float, float, float] + direction: tuple[float, float, float] + length: float + radius: float + + def __post_init__(self) -> None: + values = (*self.center, *self.direction, self.length, self.radius) + if not all(math.isfinite(value) for value in values): + raise ValueError("capsule geometry must be finite") + if len(self.center) != 3 or len(self.direction) != 3: + raise ValueError("capsule vectors must have three coordinates") + if self.length < 0 or self.radius <= 0 or math.hypot(*self.direction) == 0: + raise ValueError("capsule requires nonnegative length, positive radius and direction") + + +def geometric_porosity( + centers: Array, + spacing: tuple[float, float, float], + cells: Sequence[Capsule], + *, + subdivisions: int = 8, + walls: Sequence[bool] | None = None, +) -> Array: + """Midpoint quadrature of the union of capsules, clipped to fluid voxels. + + centers use native lattice-center convention. Walls use the existing binary + voxel mask, not a second interpretation of mechanical constraint surfaces. + Overlaps count once. The result is a coarse storage fraction, not a resolved + aperture or a sub-voxel connectivity claim. + """ + points = np.asarray(centers, dtype=np.float64) + if points.ndim != 2 or points.shape[1] != 3 or not np.isfinite(points).all(): + raise ValueError("centers must be finite N by 3 coordinates") + if len(spacing) != 3 or any(not math.isfinite(h) or h <= 0 for h in spacing): + raise ValueError("spacing must contain three positive finite lengths") + if ( + isinstance(subdivisions, bool) + or not isinstance(cast(object, subdivisions), int) + or subdivisions < 1 + ): + raise ValueError("subdivisions must be a positive integer") + solid = np.zeros(len(points), dtype=np.bool_) if walls is None else np.asarray(walls) + if solid.shape != (len(points),) or solid.dtype != np.bool_: + raise ValueError("walls must contain one Boolean per voxel") + samples = (np.arange(subdivisions, dtype=np.float64) + 0.5) / subdivisions - 0.5 + offsets = np.stack(np.meshgrid(samples, samples, samples, indexing="ij"), axis=-1) + offsets = offsets.reshape(-1, 3) * np.asarray(spacing) + result = np.zeros(len(points)) + for index, center in enumerate(points): + if solid[index]: + continue + coordinates = center + offsets + occupied = np.zeros(len(coordinates), dtype=np.bool_) + for cell in cells: + direction = np.asarray(cell.direction, dtype=np.float64) + direction = direction / np.linalg.norm(direction) + relative = coordinates - np.asarray(cell.center) + axial = np.clip(relative @ direction, -cell.length / 2, cell.length / 2) + distance = relative - axial[:, None] * direction + squared = cast(Array, np.sum(distance * distance, axis=1)) + within = squared <= cell.radius * cell.radius + occupied = np.logical_or(occupied, within) + result[index] = 1.0 - np.mean(occupied) + result[result < EPSILON_CUTOFF] = 0 + return result + + +def _vector(values: Sequence[float] | Array, name: str, count: int | None = None) -> Array: + result = np.asarray(values, dtype=np.float64) + if result.ndim != 1 or not np.isfinite(result).all(): + raise ValueError(f"{name} must be a finite vector") + if count is not None and result.shape != (count,): + raise ValueError(f"{name} size mismatch") + return result + + +def accessible_volumes(porosity: Sequence[float] | Array, voxel_volume: float) -> Array: + epsilon = _vector(porosity, "porosity") + if np.any(epsilon < 0) or np.any(epsilon > 1): + raise ValueError("porosity must lie in [0, 1]") + if not math.isfinite(voxel_volume) or voxel_volume <= 0: + raise ValueError("voxel volume must be finite and positive") + return np.where(epsilon < EPSILON_CUTOFF, 0.0, epsilon) * voxel_volume + + +def concentration(amount: Sequence[float] | Array, volume: Sequence[float] | Array) -> Array: + n = _vector(amount, "amount") + w = _vector(volume, "accessible volume", len(n)) + if np.any(n < 0) or np.any(w < 0) or np.any((w == 0) & (n != 0)): + raise ValueError("nonnegative amounts require accessible storage") + return np.divide(n, w, out=np.zeros_like(n), where=w > 0) + + +def remap_amounts( + amount: Sequence[float] | Array, + old_volume: Sequence[float] | Array, + new_volume: Sequence[float] | Array, + neighbors: Sequence[tuple[int, int]], +) -> Array: + """Keep surviving voxel amounts; expel closing storage conservatively. + + Recipients are all newly accessible voxels in the old-or-new accessible + face-connected component, weighted by new accessible volume. A closing + component with nonzero amount fails atomically. Inputs are never mutated. + """ + n = _vector(amount, "amount") + old = _vector(old_volume, "old volume", len(n)) + new = _vector(new_volume, "new volume", len(n)) + concentration(n, old) + if np.any(new < 0): + raise ValueError("new volume must be nonnegative") + result = n.copy() + adjacency: list[list[int]] = [[] for _ in n] + for first, second in neighbors: + if first == second or not 0 <= first < len(n) or not 0 <= second < len(n): + raise ValueError("invalid neighbor edge") + adjacency[first].append(second) + adjacency[second].append(first) + active = (old > 0) | (new > 0) + visited: set[int] = set() + for start in range(len(n)): + if not active[start] or start in visited: + continue + pending = [start] + component: list[int] = [] + while pending: + index = pending.pop() + if index in visited or not active[index]: + continue + visited.add(index) + component.append(index) + pending.extend(adjacency[index]) + component.sort() + donors = [index for index in component if new[index] == 0] + recipients = [index for index in component if new[index] > 0] + expelled = math.fsum(float(n[index]) for index in donors) + if expelled > 0 and not recipients: + raise ValueError( + f"closing component at voxel {start} has solute but no accessible recipient" + ) + result[donors] = 0 + if expelled: + capacity = math.fsum(float(new[index]) for index in recipients) + for index in recipients: + result[index] += expelled * new[index] / capacity + concentration(result, new) + return result + + +@dataclass(frozen=True) +class Face: + """Internal oriented face: diffusive conductance L^3/T, fluid flux L^3/T.""" + + first: int + second: int + conductance: float + volume_flux: float = 0.0 + + +def porosity_face( + first: int, + second: int, + epsilon_first: float, + epsilon_second: float, + *, + diffusion: float, + area: float, + distance: float, + intrinsic_velocity: float = 0.0, +) -> Face: + """Harmonic porosity closure; aperture is applied exactly once to flux.""" + values = (epsilon_first, epsilon_second, diffusion, area, distance, intrinsic_velocity) + if not all(math.isfinite(value) for value in values): + raise ValueError("face data must be finite") + if not 0 <= epsilon_first <= 1 or not 0 <= epsilon_second <= 1: + raise ValueError("face porosities must lie in [0, 1]") + if diffusion < 0 or area <= 0 or distance <= 0: + raise ValueError("invalid face geometry or diffusion") + aperture = ( + 0.0 + if min(epsilon_first, epsilon_second) < EPSILON_CUTOFF + else 2 * epsilon_first * epsilon_second / (epsilon_first + epsilon_second) + ) + return Face( + first, second, diffusion * aperture * area / distance, aperture * area * intrinsic_velocity + ) + + +@dataclass(frozen=True) +class ReservoirFace: + """Exterior reservoir: positive flux leaves domain; concentration is amount/L^3.""" + + site: int + concentration: float + conductance: float = 0.0 + volume_flux: float = 0.0 + + +@dataclass(frozen=True) +class Balance: + before: float + after: float + source: float + reaction: float + boundary: float + + @property + def residual(self) -> float: + return self.after - self.before - self.source - self.reaction - self.boundary + + +def backward_euler( + amount: Sequence[float] | Array, + volume: Sequence[float] | Array, + faces: Sequence[Face], + dt: float, + *, + source: Sequence[float] | Array | None = None, + loss: Sequence[float] | Array | None = None, + reservoirs: Sequence[ReservoirFace] = (), +) -> tuple[Array, Balance]: + """Dense float64 finite-volume reference with an explicit amount ledger. + + source is amount/time, loss is 1/time. Interior faces are equal/opposite; + first-order advection and reservoir/loss terms are implicit. Not scalable. + """ + n = _vector(amount, "amount") + w = _vector(volume, "volume", len(n)) + concentration(n, w) + if not math.isfinite(dt) or dt < 0: + raise ValueError("dt must be finite and nonnegative") + s = np.zeros_like(n) if source is None else _vector(source, "source", len(n)) + k = np.zeros_like(n) if loss is None else _vector(loss, "loss", len(n)) + if np.any(k < 0) or np.any((w == 0) & (s != 0)): + raise ValueError("loss must be nonnegative; sources require accessible storage") + operator = np.diag(k * w) + rhs = n + dt * s + for face in faces: + i, j, g, q = face.first, face.second, face.conductance, face.volume_flux + if i == j or not 0 <= i < len(n) or not 0 <= j < len(n): + raise ValueError("invalid transport face indices") + if not math.isfinite(g) or g < 0 or not math.isfinite(q): + raise ValueError("invalid transport coefficients") + if (w[i] == 0 or w[j] == 0) and (g != 0 or q != 0): + raise ValueError("closed storage cannot have an open face") + operator[i, i] += g + max(q, 0) + operator[j, j] += g + max(-q, 0) + operator[i, j] -= g + max(-q, 0) + operator[j, i] -= g + max(q, 0) + for face in reservoirs: + i, c, g, q = face.site, face.concentration, face.conductance, face.volume_flux + if not 0 <= i < len(n) or w[i] == 0: + raise ValueError("reservoir must connect accessible storage") + if not all(math.isfinite(value) for value in (c, g, q)) or c < 0 or g < 0: + raise ValueError("invalid reservoir coefficients") + operator[i, i] += g + max(q, 0) + rhs[i] += dt * (g + max(-q, 0)) * c + matrix = np.diag(w) + dt * operator + for i in range(len(n)): + if w[i] == 0: + matrix[i, i] = 1.0 + c = np.linalg.solve(matrix, rhs) + updated = c * w + if not np.isfinite(updated).all() or np.any(updated < 0): + raise ValueError("step produced invalid amount; no clipping is permitted") + boundary = dt * math.fsum( + f.conductance * (f.concentration - c[f.site]) + - f.volume_flux * (c[f.site] if f.volume_flux >= 0 else f.concentration) + for f in reservoirs + ) + return updated, Balance( + float(n.sum()), + float(updated.sum()), + float(dt * s.sum()), + float(-dt * np.dot(k, updated)), + float(boundary), + ) + + +def exchange_weights( + kernel: Sequence[float] | Array, + volume: Sequence[float] | Array, +) -> Array: + """Accessible-volume weighted partition of unity in a declared connected support.""" + base = _vector(kernel, "kernel") + accessible = _vector(volume, "volume", len(base)) + if np.any(base < 0) or np.any(accessible < 0): + raise ValueError("exchange kernel and volume must be nonnegative") + weights = base * accessible + if np.any(weights < 0) or not np.isfinite(weights).all() or weights.sum() <= 0: + raise ValueError("cell has no valid accessible exchange support") + return weights / weights.sum() diff --git a/python/tests/test_occupancy_reference.py b/python/tests/test_occupancy_reference.py new file mode 100644 index 0000000..880c2ad --- /dev/null +++ b/python/tests/test_occupancy_reference.py @@ -0,0 +1,178 @@ +from __future__ import annotations + +import math + +import numpy as np +import pytest +from microsimulator.occupancy_reference import ( + Capsule, + Face, + ReservoirFace, + accessible_volumes, + backward_euler, + concentration, + exchange_weights, + geometric_porosity, + porosity_face, + remap_amounts, +) + + +def test_empty_grid_matches_existing_native_backward_euler() -> None: + from microsimulator import GridShape, SignalGridSpec, SignalIntegrationKind, Simulation, Vec3 + + shape = GridShape() + shape.x, shape.y, shape.z = 3, 1, 1 + spec = SignalGridSpec() + spec.shape, spec.signal_count = shape, 1 + spec.diffusion, spec.advection = [1.0], [Vec3()] + spec.integration = SignalIntegrationKind.BACKWARD_EULER + simulation = Simulation() + simulation.configure_signal_grid(spec, [0.0, 1.0, 0.0]) + simulation.step(0.1) + result, balance = backward_euler([0, 1, 0], [1, 1, 1], [Face(0, 1, 1), Face(1, 2, 1)], 0.1) + np.testing.assert_allclose(result, simulation.signal_levels, rtol=2e-6, atol=2e-7) + assert abs(balance.residual) < 1e-12 + centers = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]]) + np.testing.assert_array_equal(geometric_porosity(centers, (1, 1, 1), []), [1, 1]) + + +def test_partial_closed_grid_uses_amount_and_unequal_storage() -> None: + volume = accessible_volumes([0.25, 0.75], 2.0) + amount = volume * np.array([8.0, 0.0]) + face = porosity_face(0, 1, 0.25, 0.75, diffusion=1, area=1, distance=1) + initial = amount.sum() + for _ in range(100): + amount, ledger = backward_euler(amount, volume, [face], 1.0) + assert abs(ledger.residual) < 1e-12 + np.testing.assert_allclose(concentration(amount, volume), [2, 2], atol=1e-12) + assert abs(amount.sum() - initial) < 1e-12 + assert not math.isclose(float(amount[0]), float(amount[1])) + + +def test_changing_occupancy_and_full_closure_conserve_without_division_by_zero() -> None: + amount = np.array([2.0, 3.0, 1.0]) + old = np.ones(3) + new = accessible_volumes([0, 0.25, 0.75], 1) + result = remap_amounts(amount, old, new, [(0, 1), (1, 2)]) + np.testing.assert_allclose(result, [0, 3.5, 2.5], atol=1e-15) + assert result.sum() == amount.sum() + assert np.isfinite(concentration(result, new)).all() + reopened = remap_amounts(result, new, [1, 1, 1], [(0, 1), (1, 2)]) + assert reopened[0] == 0 # no invented solute in newly exposed storage + assert reopened.sum() == amount.sum() + with pytest.raises(ValueError, match="no accessible recipient"): + remap_amounts(amount, old, [0, 0, 0], [(0, 1), (1, 2)]) + np.testing.assert_array_equal(amount, [2, 3, 1]) # rejection is atomic + # A persistent wall separates the only potential recipient. + with pytest.raises(ValueError, match="no accessible recipient"): + remap_amounts([1, 0, 0], [1, 0, 1], [0, 0, 1], [(0, 1), (1, 2)]) + np.testing.assert_array_equal(accessible_volumes([1e-12, 0], 1), [0, 0]) + np.testing.assert_array_equal(remap_amounts([0], [1], [0], []), [0]) + + +def test_geometry_union_wall_clipping_division_and_removal() -> None: + centers = np.array( + [(x, y, z) for x in range(-3, 4) for y in range(-1, 2) for z in range(-1, 2)], + dtype=np.float64, + ) + parent = Capsule((0, 0, 0), (1, 0, 0), 4, 0.4) + daughters = [ + Capsule((-1.2, 0, 0), (1, 0, 0), 1.6, 0.4), + Capsule((1.2, 0, 0), (1, 0, 0), 1.6, 0.4), + ] + old = geometric_porosity(centers, (1, 1, 1), [parent], subdivisions=16) + np.testing.assert_array_equal( + old, geometric_porosity(centers, (1, 1, 1), [parent, parent], subdivisions=16) + ) + divided = geometric_porosity(centers, (1, 1, 1), daughters, subdivisions=16) + assert divided.sum() > old.sum() # native division reduces geometric solid volume + b_parent = math.pi * 0.4**2 * (4 + 0.8) + b_daughters = 2 * math.pi * 0.4**2 * (1.6 + 0.8) + assert math.isclose(b_parent, b_daughters) + assert math.isclose( + (math.pi * 0.4**2 * 4 + 4 / 3 * math.pi * 0.4**3) + - 2 * (math.pi * 0.4**2 * 1.6 + 4 / 3 * math.pi * 0.4**3), + 2 / 3 * math.pi * 0.4**3, + ) + # No voxel closes in this geometry transition; amount remains voxel-local. + amount = old * 2 + updated = remap_amounts(amount, old, divided, []) + removed = remap_amounts(updated, divided, np.ones_like(old), []) + assert abs(removed.sum() - amount.sum()) < 1e-12 + wall_mask = [bool(index == 31) for index in range(len(centers))] + with_wall = geometric_porosity(centers, (1, 1, 1), [parent], walls=wall_mask) + assert with_wall[31] == 0 + assert np.all((with_wall >= 0) & (with_wall <= 1)) + + +def test_boundary_reaction_and_cell_exchange_have_explicit_amount_ledgers() -> None: + volume = np.array([0.25, 0.75]) + weights = exchange_weights([0.5, 0.5], volume) + c = np.array([2.0, 4.0]) + source = weights * 3.0 + assert float(weights @ c) == 3.5 + assert float(source.sum()) == 3.0 + updated, ledger = backward_euler( + c * volume, + volume, + [Face(0, 1, 0.1, 0.05)], + 0.2, + source=source, + loss=[0.3, 0.7], + reservoirs=[ReservoirFace(0, 5, 0.2, -0.1), ReservoirFace(1, 0, 0, 0.1)], + ) + assert np.all(updated >= 0) + assert ledger.boundary > 0 and ledger.reaction < 0 and ledger.source > 0 + assert abs(ledger.residual) < 1e-12 + face = porosity_face(0, 1, 0.25, 0.75, diffusion=2, area=3, distance=4, intrinsic_velocity=5) + assert face.volume_flux == 0.375 * 3 * 5 # porosity appears exactly once + with pytest.raises(ValueError, match="accessible exchange"): + exchange_weights([1, 0], [0, 1]) + + +def test_diffusion_timestep_refinement_and_geometric_quadrature_refinement() -> None: + # Two-cell antisymmetric diffusion mode has eigenvalue -2 for epsilon=1. + errors: list[float] = [] + for steps in (10, 20, 40): + amount = np.array([1.5, 0.5]) + for _ in range(steps): + amount, _ = backward_euler(amount, [1, 1], [Face(0, 1, 1)], 1 / steps) + errors.append(abs(float(amount[0]) - (1 + 0.5 * math.exp(-2)))) + assert errors[2] < 0.55 * errors[1] < 0.31 * errors[0] + sphere = Capsule((0, 0, 0), (1, 0, 0), 0, 0.5) + exact = 4 / 3 * math.pi * 0.5**3 + geometric_errors: list[float] = [] + for resolution in (8, 16, 32): + epsilon = geometric_porosity(np.zeros((1, 3)), (2, 2, 2), [sphere], subdivisions=resolution) + geometric_errors.append(abs(float((1 - epsilon[0]) * 8) - exact)) + assert geometric_errors[-1] < geometric_errors[0] + assert geometric_errors[-1] < 0.02 * exact + + +def test_invalid_reference_inputs_fail_without_silent_clipping() -> None: + with pytest.raises(ValueError): + accessible_volumes([1.1], 1) + with pytest.raises(ValueError): + concentration([1], [0]) + with pytest.raises(ValueError): + backward_euler([0], [1], [], 1, source=[-1]) + with pytest.raises(ValueError): + backward_euler([0, 1], [0, 1], [Face(0, 1, 1)], 0.1) + + +def test_empty_limit_spatial_operator_is_second_order() -> None: + errors: list[float] = [] + for count in (10, 20, 40): + h = 1 / count + centers = (np.arange(count, dtype=np.float64) + 0.5) * h + values = 2 + np.cos(math.pi * centers) + rate = np.zeros(count) + for i in range(count - 1): + face = porosity_face(i, i + 1, 1, 1, diffusion=1, area=1, distance=h) + amount_flux = face.conductance * (values[i + 1] - values[i]) + rate[i] += amount_flux / h + rate[i + 1] -= amount_flux / h + exact = -(math.pi**2) * np.cos(math.pi * centers) + errors.append(float(np.max(np.abs(rate - exact)))) + assert errors[2] < 0.26 * errors[1] < 0.07 * errors[0]