|
| 1 | +"""Thin wrapper around the ``pyvcell_mbsolver`` moving-boundary solver. |
| 2 | +
|
| 3 | +The solver reports its solution through observer callbacks rather than an output |
| 4 | +file: once per internal time step it hands back the moving front geometry and |
| 5 | +the per-element field values. :func:`solve_moving_boundary` runs the solver with |
| 6 | +a collecting observer that snapshots those callbacks at the requested output |
| 7 | +times into a :class:`~pyvcell.sim_results.moving_boundary_result.MovingBoundaryResult`. |
| 8 | +
|
| 9 | +Callers who need full control can use ``pyvcell_mbsolver`` directly; this wrapper |
| 10 | +covers the common "run it and give me the trajectory" case. |
| 11 | +""" |
| 12 | + |
| 13 | +from __future__ import annotations |
| 14 | + |
| 15 | +from collections.abc import Sequence |
| 16 | +from os import PathLike |
| 17 | +from typing import Any |
| 18 | + |
| 19 | +import numpy as np |
| 20 | +import pyvcell_mbsolver as mb |
| 21 | + |
| 22 | +from pyvcell.sim_results.moving_boundary_result import MovingBoundaryFrame, MovingBoundaryResult |
| 23 | + |
| 24 | +# Time tolerance (relative to the output step) for deciding a step has reached an |
| 25 | +# output time; the solver advances on its own internal step, not the output step. |
| 26 | +_OUTPUT_TIME_TOL = 1e-9 |
| 27 | + |
| 28 | + |
| 29 | +def solve_moving_boundary( |
| 30 | + setup_xml_file: PathLike[str] | str, |
| 31 | + species_names: Sequence[str], |
| 32 | + output_times: Sequence[float] | None = None, |
| 33 | +) -> MovingBoundaryResult: |
| 34 | + """Run the moving-boundary solver on a ``MovingBoundarySetup`` XML file. |
| 35 | +
|
| 36 | + Args: |
| 37 | + setup_xml_file: the ``*_mb.xml`` setup produced by |
| 38 | + ``libvcell.vcml_to_moving_boundary_input``. |
| 39 | + species_names: the volume species, in the order the solver reports |
| 40 | + concentrations (i.e. the order of the ``<species>`` entries in the |
| 41 | + setup file's ``<physiology>`` block). |
| 42 | + output_times: times at which to snapshot a frame. The first solver step |
| 43 | + at or past each requested time is kept; the final step is always |
| 44 | + kept. When None, every solver step is kept. |
| 45 | +
|
| 46 | + Returns: |
| 47 | + a :class:`MovingBoundaryResult` with one frame per kept output time. |
| 48 | + """ |
| 49 | + species = list(species_names) |
| 50 | + targets = sorted(float(t) for t in output_times) if output_times is not None else None |
| 51 | + result = MovingBoundaryResult(species_names=species) |
| 52 | + |
| 53 | + class _Collector(mb.SimulationObserver): |
| 54 | + def __init__(self) -> None: |
| 55 | + super().__init__() |
| 56 | + self._target_index = 0 |
| 57 | + self._keep = False |
| 58 | + self._time = 0.0 |
| 59 | + self._front: np.ndarray = np.empty((0, 2), dtype=float) |
| 60 | + self._buffer: list[tuple[float, float, int, int, list[float]]] = [] |
| 61 | + |
| 62 | + def on_time(self, t: float, generation: int, last: bool, geometry: Any) -> None: |
| 63 | + keep = bool(last) |
| 64 | + if targets is None: |
| 65 | + keep = True |
| 66 | + else: |
| 67 | + tol = _OUTPUT_TIME_TOL * (targets[-1] - targets[0] + 1.0) |
| 68 | + while self._target_index < len(targets) and t + tol >= targets[self._target_index]: |
| 69 | + keep = True |
| 70 | + self._target_index += 1 |
| 71 | + self._keep = keep |
| 72 | + if keep: |
| 73 | + self._time = t |
| 74 | + boundary = list(getattr(geometry, "boundary", []) or []) |
| 75 | + self._front = np.array(boundary, dtype=float).reshape(-1, 2) if boundary else np.empty((0, 2)) |
| 76 | + self._buffer = [] |
| 77 | + |
| 78 | + def on_element(self, node: Any) -> None: |
| 79 | + if not self._keep or node.is_outside: |
| 80 | + return |
| 81 | + self._buffer.append(( |
| 82 | + float(node.x), |
| 83 | + float(node.y), |
| 84 | + int(node.grid_i), |
| 85 | + int(node.grid_j), |
| 86 | + [float(node.concentration(i)) for i in range(len(species))], |
| 87 | + )) |
| 88 | + |
| 89 | + def on_iteration_complete(self) -> None: |
| 90 | + if not self._keep: |
| 91 | + return |
| 92 | + self._keep = False |
| 93 | + x = np.array([row[0] for row in self._buffer], dtype=float) |
| 94 | + y = np.array([row[1] for row in self._buffer], dtype=float) |
| 95 | + grid_i = np.array([row[2] for row in self._buffer], dtype=int) |
| 96 | + grid_j = np.array([row[3] for row in self._buffer], dtype=int) |
| 97 | + concentrations = { |
| 98 | + name: np.array([row[4][i] for row in self._buffer], dtype=float) for i, name in enumerate(species) |
| 99 | + } |
| 100 | + result.frames.append( |
| 101 | + MovingBoundaryFrame( |
| 102 | + time=self._time, |
| 103 | + front=self._front, |
| 104 | + x=x, |
| 105 | + y=y, |
| 106 | + grid_i=grid_i, |
| 107 | + grid_j=grid_j, |
| 108 | + concentrations=concentrations, |
| 109 | + ) |
| 110 | + ) |
| 111 | + |
| 112 | + def on_complete(self) -> None: |
| 113 | + pass |
| 114 | + |
| 115 | + solver = mb.MovingBoundarySolver.from_xml(str(setup_xml_file)) |
| 116 | + solver.add_element_observer(_Collector(), name="pyvcell_collector") |
| 117 | + solver.run() |
| 118 | + return result |
0 commit comments