Add moving-boundary solver support

docs/getting-started/installation.md

@@ -15,9 +15,13 @@ features live behind optional extras — install the ones you need, or everythin
 
 ```bash
 pip install "pyvcell[viz]"     # plotting / VTK / PyVista
-pip install "pyvcell[all]"     # full feature set (solver, viz, remote, io, convert, native)
+pip install "pyvcell[mb]"      # moving-boundary solver (pyvcell-mbsolver)
+pip install "pyvcell[all]"     # full feature set (solver, viz, remote, io, convert, native, mb)
 ```
 
+> The moving-boundary solver (`mb`) pulls `pyvcell-mbsolver` plus
+> `libvcell >= 0.0.17` (the `native` extra) for moving-boundary input generation.
+
 ## Install for development (uv)
 
 ```bash

pyproject.toml

@@ -59,10 +59,17 @@ convert = [
     "sympy>=1.13.1,<2",
 ]
 native = [
-    "libvcell (>=0.0.15.3)",
+    # 0.0.17 ships universal py3-none wheels and includes moving-boundary support.
+    "libvcell (>=0.0.17)",
+]
+# Moving-boundary solver. libvcell (the `native` extra) generates the solver
+# input; pyvcell-mbsolver runs it.
+mb = [
+    "pyvcell-mbsolver>=1.0.3,<2",
+    "pyvcell[native]",
 ]
 all = [
-    "pyvcell[solver,viz,remote,io,convert,native]",
+    "pyvcell[solver,viz,remote,io,convert,native,mb]",
 ]
 
 [dependency-groups]

pyvcell/_internal/solvers/mbsolver.py

@@ -0,0 +1,118 @@
+"""Thin wrapper around the ``pyvcell_mbsolver`` moving-boundary solver.
+
+The solver reports its solution through observer callbacks rather than an output
+file: once per internal time step it hands back the moving front geometry and
+the per-element field values. :func:`solve_moving_boundary` runs the solver with
+a collecting observer that snapshots those callbacks at the requested output
+times into a :class:`~pyvcell.sim_results.moving_boundary_result.MovingBoundaryResult`.
+
+Callers who need full control can use ``pyvcell_mbsolver`` directly; this wrapper
+covers the common "run it and give me the trajectory" case.
+"""
+
+from __future__ import annotations
+
+from collections.abc import Sequence
+from os import PathLike
+from typing import Any
+
+import numpy as np
+import pyvcell_mbsolver as mb
+
+from pyvcell.sim_results.moving_boundary_result import MovingBoundaryFrame, MovingBoundaryResult
+
+# Time tolerance (relative to the output step) for deciding a step has reached an
+# output time; the solver advances on its own internal step, not the output step.
+_OUTPUT_TIME_TOL = 1e-9
+
+
+def solve_moving_boundary(
+    setup_xml_file: PathLike[str] | str,
+    species_names: Sequence[str],
+    output_times: Sequence[float] | None = None,
+) -> MovingBoundaryResult:
+    """Run the moving-boundary solver on a ``MovingBoundarySetup`` XML file.
+
+    Args:
+        setup_xml_file: the ``*_mb.xml`` setup produced by
+            ``libvcell.vcml_to_moving_boundary_input``.
+        species_names: the volume species, in the order the solver reports
+            concentrations (i.e. the order of the ``<species>`` entries in the
+            setup file's ``<physiology>`` block).
+        output_times: times at which to snapshot a frame. The first solver step
+            at or past each requested time is kept; the final step is always
+            kept. When None, every solver step is kept.
+
+    Returns:
+        a :class:`MovingBoundaryResult` with one frame per kept output time.
+    """
+    species = list(species_names)
+    targets = sorted(float(t) for t in output_times) if output_times is not None else None
+    result = MovingBoundaryResult(species_names=species)
+
+    class _Collector(mb.SimulationObserver):
+        def __init__(self) -> None:
+            super().__init__()
+            self._target_index = 0
+            self._keep = False
+            self._time = 0.0
+            self._front: np.ndarray = np.empty((0, 2), dtype=float)
+            self._buffer: list[tuple[float, float, int, int, list[float]]] = []
+
+        def on_time(self, t: float, generation: int, last: bool, geometry: Any) -> None:
+            keep = bool(last)
+            if targets is None:
+                keep = True
+            else:
+                tol = _OUTPUT_TIME_TOL * (targets[-1] - targets[0] + 1.0)
+                while self._target_index < len(targets) and t + tol >= targets[self._target_index]:
+                    keep = True
+                    self._target_index += 1
+            self._keep = keep
+            if keep:
+                self._time = t
+                boundary = list(getattr(geometry, "boundary", []) or [])
+                self._front = np.array(boundary, dtype=float).reshape(-1, 2) if boundary else np.empty((0, 2))
+                self._buffer = []
+
+        def on_element(self, node: Any) -> None:
+            if not self._keep or node.is_outside:
+                return
+            self._buffer.append((
+                float(node.x),
+                float(node.y),
+                int(node.grid_i),
+                int(node.grid_j),
+                [float(node.concentration(i)) for i in range(len(species))],
+            ))
+
+        def on_iteration_complete(self) -> None:
+            if not self._keep:
+                return
+            self._keep = False
+            x = np.array([row[0] for row in self._buffer], dtype=float)
+            y = np.array([row[1] for row in self._buffer], dtype=float)
+            grid_i = np.array([row[2] for row in self._buffer], dtype=int)
+            grid_j = np.array([row[3] for row in self._buffer], dtype=int)
+            concentrations = {
+                name: np.array([row[4][i] for row in self._buffer], dtype=float) for i, name in enumerate(species)
+            }
+            result.frames.append(
+                MovingBoundaryFrame(
+                    time=self._time,
+                    front=self._front,
+                    x=x,
+                    y=y,
+                    grid_i=grid_i,
+                    grid_j=grid_j,
+                    concentrations=concentrations,
+                )
+            )
+
+        def on_complete(self) -> None:
+            pass
+
+    solver = mb.MovingBoundarySolver.from_xml(str(setup_xml_file))
+    solver.add_element_observer(_Collector(), name="pyvcell_collector")
+    solver.run()
+    return result

pyvcell/sim_results/moving_boundary_result.py

@@ -0,0 +1,73 @@
+"""Result of a moving-boundary simulation.
+
+Unlike the finite-volume solver, the moving-boundary solver does not write a
+fixed Cartesian grid. Its solution lives on a mesh that moves with the boundary,
+so a result is a time series of frames, where each frame carries the moving
+front (a polygon of ``(x, y)`` points) and the field values of the "inside"
+mesh elements at that time. This mirrors VCell's own moving-boundary data model
+(elements + per-species values on the moving mesh), rather than the
+zarr/4-D Cartesian :class:`~pyvcell.sim_results.result.Result`.
+
+These objects hold only numpy arrays and plain Python, so they can be inspected
+and post-processed without any heavy optional dependency.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass, field
+
+import numpy as np
+
+
+@dataclass
+class MovingBoundaryFrame:
+    """A single output time step of a moving-boundary simulation.
+
+    Attributes:
+        time: the simulation time of this frame.
+        front: ``(N, 2)`` array of the moving-front boundary polygon ``(x, y)``.
+        x, y: ``(M,)`` coordinates of the inside mesh elements.
+        grid_i, grid_j: ``(M,)`` integer grid indices of the inside elements.
+        concentrations: species name -> ``(M,)`` array of element concentrations,
+            aligned with ``x``/``y``/``grid_i``/``grid_j``.
+    """
+
+    time: float
+    front: np.ndarray
+    x: np.ndarray
+    y: np.ndarray
+    grid_i: np.ndarray
+    grid_j: np.ndarray
+    concentrations: dict[str, np.ndarray] = field(default_factory=dict)
+
+
+@dataclass
+class MovingBoundaryResult:
+    """The full time series of a moving-boundary simulation.
+
+    Attributes:
+        species_names: the volume species, in the same order the solver reports
+            concentrations.
+        frames: the output-time frames, in increasing time order.
+    """
+
+    species_names: list[str]
+    frames: list[MovingBoundaryFrame] = field(default_factory=list)
+
+    @property
+    def times(self) -> np.ndarray:
+        """The output times as a 1-D array."""
+        return np.array([frame.time for frame in self.frames], dtype=float)
+
+    def front(self, time_index: int) -> np.ndarray:
+        """The moving-front polygon ``(N, 2)`` at the given output index."""
+        return self.frames[time_index].front
+
+    def concentrations(self, species_name: str, time_index: int) -> np.ndarray:
+        """Inside-element concentrations of ``species_name`` at the given output index."""
+        return self.frames[time_index].concentrations[species_name]
+
+    def __repr__(self) -> str:
+        n = len(self.frames)
+        span = f"{self.frames[0].time:g}..{self.frames[-1].time:g}" if n else "empty"
+        return f"MovingBoundaryResult(species={self.species_names}, frames={n}, t={span})"

pyvcell/vcml/__init__.py

@@ -21,10 +21,12 @@
     Biomodel,
     BoundaryType,
     Compartment,
+    FrontVelocity,
     Kinetics,
     KineticsParameter,
     Model,
     ModelParameter,
+    MovingBoundarySolverOptions,
     Reaction,
     Simulation,
     Species,
@@ -68,6 +70,7 @@
 if TYPE_CHECKING:
     # Heavy names — eager only for type checkers / IDE autocomplete.
     from pyvcell._internal.geometry import SegmentedImageGeometry
+    from pyvcell.sim_results.moving_boundary_result import MovingBoundaryResult
     from pyvcell.vcml.field import Field
     from pyvcell.vcml.session import SimulationJob, VCellSession
     from pyvcell.vcml.utils import (
@@ -90,6 +93,7 @@
         write_sbml_file,
         write_vcml_file,
     )
+    from pyvcell.vcml.vcml_mb_simulation import simulate_moving_boundary
     from pyvcell.vcml.vcml_remote import connect, logout
     from pyvcell.vcml.vcml_simulation import cartesian_mesh_from_geometry, simulate
     from pyvcell.vcml.vcml_writer import VcmlWriter
@@ -104,6 +108,8 @@
     **dict.fromkeys(["SimulationJob", "VCellSession"], "pyvcell.vcml.session"),
     **dict.fromkeys(["connect", "logout"], "pyvcell.vcml.vcml_remote"),
     **dict.fromkeys(["simulate", "cartesian_mesh_from_geometry"], "pyvcell.vcml.vcml_simulation"),
+    "simulate_moving_boundary": "pyvcell.vcml.vcml_mb_simulation",
+    "MovingBoundaryResult": "pyvcell.sim_results.moving_boundary_result",
     **dict.fromkeys(["get_workspace_dir", "set_workspace_dir"], "pyvcell.vcml.workspace"),
     **dict.fromkeys(
         [
@@ -135,6 +141,7 @@
     "libvcell": "native",
     "pyvcell_fvsolver": "solver",
     "fvsolver": "solver",
+    "pyvcell_mbsolver": "mb",
     "vtk": "viz",
     "pyvista": "viz",
     "matplotlib": "viz",
@@ -192,6 +199,7 @@ def __dir__() -> list[str]:
     "Constant",
     "Effect",
     "Field",
+    "FrontVelocity",
     "Geometry",
     "Image",
     "JumpCondition",
@@ -206,6 +214,8 @@ def __dir__() -> list[str]:
     "MembraneSubDomain",
     "Model",
     "ModelParameter",
+    "MovingBoundaryResult",
+    "MovingBoundarySolverOptions",
     "OdeEquation",
     "ParticleInitialCount",
     "ParticleJumpProcess",
@@ -246,6 +256,7 @@ def __dir__() -> list[str]:
     "restore_stdout",
     "set_workspace_dir",
     "simulate",
+    "simulate_moving_boundary",
     "suppress_stdout",
     "to_antimony_str",
     "to_sbml_str",

pyvcell/vcml/models.py

@@ -14,6 +14,12 @@
 from pyvcell.vcml.models_app import (
     CompartmentMapping as CompartmentMapping,
 )
+from pyvcell.vcml.models_app import (
+    FrontVelocity as FrontVelocity,
+)
+from pyvcell.vcml.models_app import (
+    MovingBoundarySolverOptions as MovingBoundarySolverOptions,
+)
 from pyvcell.vcml.models_app import (
     ReactionMapping as ReactionMapping,
 )