Examples

The voids.examples sub-package provides deterministic synthetic networks and images intended for testing, documentation, and reproducible demonstrations.

---

Demo Networks

voids.examples.demo

make_linear_chain_network

make_linear_chain_network(
    num_pores=3,
    *,
    axis="x",
    length=1.0,
    cross_section=1.0,
    bulk_volume=10.0,
    pore_volume=1.0,
    throat_volume=0.5,
    throat_length=1.0,
    hydraulic_conductance=1.0,
)

Build a deterministic one-dimensional pore-throat chain.

Parameters:

Name	Type	Description	Default
num_pores	int	Number of pores in the chain. The number of throats is num_pores - 1.	3
axis	str	Axis along which the chain is embedded.	'x'
length	float	Sample length along the chosen axis.	1.0
cross_section	float	Cross-sectional area normal to the flow axis.	1.0
bulk_volume	float	Bulk sample volume associated with the toy problem.	10.0
pore_volume	float	Synthetic pore and throat void volumes.	1.0
throat_volume	float	Synthetic pore and throat void volumes.	1.0
throat_length	float	Length assigned to each throat.	1.0
hydraulic_conductance	float	Precomputed throat hydraulic conductance.	1.0

Returns:

Type	Description
Network	Synthetic line network with canonical inlet and outlet labels.

Raises:

Type	Description
ValueError	If the number of pores, axis, or geometric parameters are invalid.

Notes

The pore coordinates are uniformly spaced so that the pore positions satisfy

x_k = k * length / (num_pores - 1)

along the selected axis. The function is intended for solver smoke tests, tutorials, and regression examples rather than realistic porous-media reconstruction.

Source code in src/voids/examples/demo.py

def make_linear_chain_network(
    num_pores: int = 3,
    *,
    axis: str = "x",
    length: float = 1.0,
    cross_section: float = 1.0,
    bulk_volume: float = 10.0,
    pore_volume: float = 1.0,
    throat_volume: float = 0.5,
    throat_length: float = 1.0,
    hydraulic_conductance: float = 1.0,
) -> Network:
    """Build a deterministic one-dimensional pore-throat chain.

    Parameters
    ----------
    num_pores :
        Number of pores in the chain. The number of throats is
        ``num_pores - 1``.
    axis :
        Axis along which the chain is embedded.
    length :
        Sample length along the chosen axis.
    cross_section :
        Cross-sectional area normal to the flow axis.
    bulk_volume :
        Bulk sample volume associated with the toy problem.
    pore_volume, throat_volume :
        Synthetic pore and throat void volumes.
    throat_length :
        Length assigned to each throat.
    hydraulic_conductance :
        Precomputed throat hydraulic conductance.

    Returns
    -------
    Network
        Synthetic line network with canonical inlet and outlet labels.

    Raises
    ------
    ValueError
        If the number of pores, axis, or geometric parameters are invalid.

    Notes
    -----
    The pore coordinates are uniformly spaced so that the pore positions satisfy

    ``x_k = k * length / (num_pores - 1)``

    along the selected axis. The function is intended for solver smoke tests,
    tutorials, and regression examples rather than realistic porous-media
    reconstruction.
    """

    if num_pores < 2:
        raise ValueError("num_pores must be >= 2")
    if axis not in _AXIS_TO_INDEX:
        raise ValueError("axis must be one of 'x', 'y', or 'z'")
    if length <= 0 or cross_section <= 0 or bulk_volume <= 0:
        raise ValueError("length, cross_section, and bulk_volume must be positive")
    if pore_volume < 0 or throat_volume < 0 or throat_length < 0 or hydraulic_conductance < 0:
        raise ValueError("pore/throat properties must be nonnegative")

    coords = np.zeros((num_pores, 3), dtype=float)
    coords[:, _AXIS_TO_INDEX[axis]] = np.linspace(0.0, float(length), num_pores)
    throat_conns = np.column_stack(
        [np.arange(num_pores - 1, dtype=np.int64), np.arange(1, num_pores, dtype=np.int64)]
    )

    pore_labels: dict[str, np.ndarray] = {
        f"inlet_{axis}min": np.zeros(num_pores, dtype=bool),
        f"outlet_{axis}max": np.zeros(num_pores, dtype=bool),
        "boundary": np.zeros(num_pores, dtype=bool),
    }
    pore_labels[f"inlet_{axis}min"][0] = True
    pore_labels[f"outlet_{axis}max"][-1] = True
    pore_labels["boundary"][[0, -1]] = True

    sample = SampleGeometry(
        bulk_volume=float(bulk_volume),
        lengths={axis: float(length)},
        cross_sections={axis: float(cross_section)},
    )
    provenance = Provenance(
        source_kind="synthetic_demo",
        extraction_method="linear_chain",
        user_notes={"num_pores": int(num_pores), "axis": axis},
    )

    return Network(
        throat_conns=throat_conns,
        pore_coords=coords,
        sample=sample,
        provenance=provenance,
        pore={"volume": np.full(num_pores, float(pore_volume), dtype=float)},
        throat={
            "volume": np.full(num_pores - 1, float(throat_volume), dtype=float),
            "length": np.full(num_pores - 1, float(throat_length), dtype=float),
            "hydraulic_conductance": np.full(
                num_pores - 1, float(hydraulic_conductance), dtype=float
            ),
        },
        pore_labels=pore_labels,
    )

---

Manufactured Void Images

voids.examples.manufactured

make_manufactured_void_image

make_manufactured_void_image(shape=(48, 48, 48))

Create a deterministic synthetic 3-D void-space image.

Parameters:

Name	Type	Description	Default
shape	tuple[int, int, int]	Output image shape in voxels.	(48, 48, 48)

Returns:

Type	Description
ndarray	Boolean array with shape shape where True denotes void space.

Notes

The construction is intentionally simple: a chain of overlapping spheres spans the x-direction, while a few side branches create off-axis connectivity. The result is not intended as a geological model. It is a manufactured test image for extraction workflows such as porespy.snow2.

Source code in src/voids/examples/manufactured.py

def make_manufactured_void_image(shape: tuple[int, int, int] = (48, 48, 48)) -> np.ndarray:
    """Create a deterministic synthetic 3-D void-space image.

    Parameters
    ----------
    shape :
        Output image shape in voxels.

    Returns
    -------
    numpy.ndarray
        Boolean array with shape ``shape`` where ``True`` denotes void space.

    Notes
    -----
    The construction is intentionally simple: a chain of overlapping spheres
    spans the x-direction, while a few side branches create off-axis
    connectivity. The result is not intended as a geological model. It is a
    manufactured test image for extraction workflows such as ``porespy.snow2``.
    """

    nx, ny, nz = shape
    X, Y, Z = np.indices(shape)
    im = np.zeros(shape, dtype=bool)

    chain = [
        (6, ny // 2, nz // 2, 7),
        (14, ny // 2 + 1, nz // 2, 7),
        (22, ny // 2 - 1, nz // 2 + 1, 7),
        (30, ny // 2, nz // 2 - 1, 7),
        (38, ny // 2 + 1, nz // 2, 7),
    ]
    branches = [
        (20, ny // 2 + 10, nz // 2, 5),
        (28, ny // 2 - 10, nz // 2 + 2, 5),
        (34, ny // 2 + 6, nz // 2 + 8, 4),
    ]
    for cx, cy, cz, r in chain + branches:
        mask = (X - cx) ** 2 + (Y - cy) ** 2 + (Z - cz) ** 2 <= r**2
        im |= mask

    y0 = ny // 2
    z0 = nz // 2
    im[12:17, y0 - 1 : y0 + 2, z0 - 1 : z0 + 2] = True

    return im

save_default_manufactured_void_image

save_default_manufactured_void_image(path)

Write the manufactured void image to a NumPy .npy file.

Parameters:

Name	Type	Description	Default
path	str | Path	Destination file path.	required

Returns:

Type	Description
Path	Resolved path that was written.

Source code in src/voids/examples/manufactured.py

def save_default_manufactured_void_image(path: str | Path) -> Path:
    """Write the manufactured void image to a NumPy ``.npy`` file.

    Parameters
    ----------
    path :
        Destination file path.

    Returns
    -------
    pathlib.Path
        Resolved path that was written.
    """

    path = Path(path)
    path.parent.mkdir(parents=True, exist_ok=True)
    np.save(path, make_manufactured_void_image())
    return path

---

Mesh Networks

voids.examples.mesh

make_cartesian_mesh_network

make_cartesian_mesh_network(
    shape,
    *,
    spacing=1.0,
    pore_radius=None,
    throat_radius=None,
    thickness=None,
    units=None,
)

Build a regular mesh-like pore network with one pore per mesh node.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	Number of pores along each active axis. Typical examples are (20, 20) and (20, 20, 20).	required
spacing	float	Center-to-center pore spacing.	1.0
pore_radius	float | None	Synthetic geometric radii used to construct pore and throat attributes.	None
throat_radius	float | None	Synthetic geometric radii used to construct pore and throat attributes.	None
thickness	float | None	Extrusion thickness for 2-D meshes. Ignored for 3-D meshes.	None
units	dict[str, str] | None	Optional unit metadata stored in :class:SampleGeometry.	None

Returns:

Type	Description
Network	Synthetic Cartesian lattice network with geometry, labels, and sample metadata.

Raises:

Type	Description
ValueError	If the shape, spacing, or geometric radii are invalid.

Notes

Each mesh node becomes one pore, and each nearest-neighbor pair becomes one throat. The resulting graph is a regular square or cubic lattice. For the current synthetic geometry model, the throat core length is

L_core = spacing - 2 * pore_radius

and the throat volume is approximated as

V_throat = A_throat * L_core.

This makes the example useful for solver verification and scaling studies, while remaining intentionally simpler than an image-derived pore network.

Source code in src/voids/examples/mesh.py

def make_cartesian_mesh_network(
    shape: Sequence[int],
    *,
    spacing: float = 1.0,
    pore_radius: float | None = None,
    throat_radius: float | None = None,
    thickness: float | None = None,
    units: dict[str, str] | None = None,
) -> Network:
    """Build a regular mesh-like pore network with one pore per mesh node.

    Parameters
    ----------
    shape :
        Number of pores along each active axis. Typical examples are ``(20, 20)``
        and ``(20, 20, 20)``.
    spacing :
        Center-to-center pore spacing.
    pore_radius, throat_radius :
        Synthetic geometric radii used to construct pore and throat attributes.
    thickness :
        Extrusion thickness for 2-D meshes. Ignored for 3-D meshes.
    units :
        Optional unit metadata stored in :class:`SampleGeometry`.

    Returns
    -------
    Network
        Synthetic Cartesian lattice network with geometry, labels, and sample
        metadata.

    Raises
    ------
    ValueError
        If the shape, spacing, or geometric radii are invalid.

    Notes
    -----
    Each mesh node becomes one pore, and each nearest-neighbor pair becomes one
    throat. The resulting graph is a regular square or cubic lattice. For the
    current synthetic geometry model, the throat core length is

    ``L_core = spacing - 2 * pore_radius``

    and the throat volume is approximated as

    ``V_throat = A_throat * L_core``.

    This makes the example useful for solver verification and scaling studies,
    while remaining intentionally simpler than an image-derived pore network.
    """

    dims = _normalize_shape(shape)
    ndim = len(dims)
    if spacing <= 0:
        raise ValueError("spacing must be positive")

    pore_radius = 0.2 * spacing if pore_radius is None else float(pore_radius)
    throat_radius = 0.1 * spacing if throat_radius is None else float(throat_radius)
    if pore_radius <= 0 or throat_radius <= 0:
        raise ValueError("pore_radius and throat_radius must be positive")
    if pore_radius >= 0.5 * spacing:
        raise ValueError("pore_radius must be smaller than half the pore spacing")
    if throat_radius >= 0.5 * spacing:
        raise ValueError("throat_radius must be smaller than half the pore spacing")

    if ndim == 2:
        nz = 1
        depth = float(spacing if thickness is None else thickness)
        if depth <= 0:
            raise ValueError("thickness must be positive for 2D meshes")
        shape3 = (dims[0], dims[1], nz)
        x = (np.arange(dims[0], dtype=float) + 0.5) * spacing
        y = (np.arange(dims[1], dtype=float) + 0.5) * spacing
        z = np.array([0.5 * depth], dtype=float)
        pore_volume_scalar = np.pi * pore_radius**2 * depth
        cross_sections = {
            "x": dims[1] * spacing * depth,
            "y": dims[0] * spacing * depth,
        }
        lengths = {
            "x": dims[0] * spacing,
            "y": dims[1] * spacing,
        }
        bulk_volume = dims[0] * dims[1] * spacing**2 * depth
    else:
        shape3 = (dims[0], dims[1], dims[2])
        x = (np.arange(dims[0], dtype=float) + 0.5) * spacing
        y = (np.arange(dims[1], dtype=float) + 0.5) * spacing
        z = (np.arange(dims[2], dtype=float) + 0.5) * spacing
        pore_volume_scalar = (4.0 / 3.0) * np.pi * pore_radius**3
        cross_sections = {
            "x": dims[1] * dims[2] * spacing**2,
            "y": dims[0] * dims[2] * spacing**2,
            "z": dims[0] * dims[1] * spacing**2,
        }
        lengths = {
            "x": dims[0] * spacing,
            "y": dims[1] * spacing,
            "z": dims[2] * spacing,
        }
        bulk_volume = dims[0] * dims[1] * dims[2] * spacing**3

    X, Y, Z = np.meshgrid(x, y, z, indexing="ij")
    pore_coords = np.column_stack([X.ravel(), Y.ravel(), Z.ravel()])
    throat_conns = _build_cartesian_connectivity(shape3, ndim=ndim)
    pore_labels = _build_boundary_labels(shape3, ndim=ndim)

    throat_area_scalar = np.pi * throat_radius**2
    throat_perimeter_scalar = 2.0 * np.pi * throat_radius
    throat_core_length_scalar = spacing - 2.0 * pore_radius
    if throat_core_length_scalar <= 0:  # pragma: no cover - guarded by pore_radius < spacing / 2
        raise ValueError(
            "pore_radius is too large relative to spacing; throat core length must stay positive"
        )

    pore_area_scalar = np.pi * pore_radius**2
    pore_perimeter_scalar = 2.0 * np.pi * pore_radius
    n_pores = pore_coords.shape[0]
    n_throats = throat_conns.shape[0]

    pore = {
        "volume": np.full(n_pores, pore_volume_scalar, dtype=float),
        "area": np.full(n_pores, pore_area_scalar, dtype=float),
        "perimeter": np.full(n_pores, pore_perimeter_scalar, dtype=float),
        "shape_factor": np.full(n_pores, _CIRCULAR_SHAPE_FACTOR, dtype=float),
        "radius_inscribed": np.full(n_pores, pore_radius, dtype=float),
        "diameter_inscribed": np.full(n_pores, 2.0 * pore_radius, dtype=float),
    }
    throat = {
        "volume": np.full(n_throats, throat_area_scalar * throat_core_length_scalar, dtype=float),
        "area": np.full(n_throats, throat_area_scalar, dtype=float),
        "perimeter": np.full(n_throats, throat_perimeter_scalar, dtype=float),
        "shape_factor": np.full(n_throats, _CIRCULAR_SHAPE_FACTOR, dtype=float),
        "radius_inscribed": np.full(n_throats, throat_radius, dtype=float),
        "diameter_inscribed": np.full(n_throats, 2.0 * throat_radius, dtype=float),
        "length": np.full(n_throats, spacing, dtype=float),
        "direct_length": np.full(n_throats, spacing, dtype=float),
        "pore1_length": np.full(n_throats, pore_radius, dtype=float),
        "core_length": np.full(n_throats, throat_core_length_scalar, dtype=float),
        "pore2_length": np.full(n_throats, pore_radius, dtype=float),
    }

    sample = SampleGeometry(
        bulk_volume=float(bulk_volume),
        lengths={k: float(v) for k, v in lengths.items()},
        cross_sections={k: float(v) for k, v in cross_sections.items()},
        units=units or {"length": "m", "pressure": "Pa"},
    )
    provenance = Provenance(
        source_kind="synthetic_mesh",
        extraction_method="cartesian_lattice",
        voxel_size_original=float(spacing),
        user_notes={"shape": list(dims)},
    )

    return Network(
        throat_conns=throat_conns,
        pore_coords=pore_coords,
        sample=sample,
        provenance=provenance,
        pore=pore,
        throat=throat,
        pore_labels=pore_labels,
        extra={
            "mesh_shape": tuple(dims),
            "mesh_spacing": float(spacing),
            "mesh_ndim": ndim,
        },
    )

---

FEM Manufactured Solutions

The mathematical definitions, convergence procedure, and interpretation caveats are documented in FEM Manufactured-Solution Verification.

voids.examples.mms

Manufactured Brinkman solutions and finite-element convergence studies.

BrinkmanMMSCase dataclass

Exact solution data for a manufactured Brinkman problem.

Parameters:

Name	Type	Description	Default
name	str	Stable case identifier.	required
dimension	Literal[2, 3]	Spatial dimension, either 2 or 3.	required
viscosity	float	Constant Brinkman diffusion coefficient :math:\nu.	required
reaction	float	Constant Darcy reaction coefficient :math:\gamma.	required
exact_solution_factory	MMSExactSolutionFactory	Callable receiving the imported ufl module and a DOLFINx mesh. It must return (velocity, pressure) as UFL expressions. The forcing is manufactured automatically as :math:-\nu\Delta u+\gamma u+\nabla p.	required
point_evaluator	MMSPointEvaluator | None	Optional NumPy evaluator for plotting. It receives coordinates with shape (dimension, npoints) and returns velocity with shape (dimension, npoints) and pressure with shape (npoints,).	None
description	str	Concise provenance suitable for notebook and metadata reporting.	''
reference	str	Concise provenance suitable for notebook and metadata reporting.	''

Notes

The convergence runner imposes the exact velocity on the complete boundary and compares pressure modulo an additive constant.

Source code in src/voids/examples/mms/_core.py

@dataclass(frozen=True, slots=True)
class BrinkmanMMSCase:
    r"""Exact solution data for a manufactured Brinkman problem.

    Parameters
    ----------
    name :
        Stable case identifier.
    dimension :
        Spatial dimension, either 2 or 3.
    viscosity :
        Constant Brinkman diffusion coefficient :math:`\nu`.
    reaction :
        Constant Darcy reaction coefficient :math:`\gamma`.
    exact_solution_factory :
        Callable receiving the imported ``ufl`` module and a DOLFINx mesh. It
        must return ``(velocity, pressure)`` as UFL expressions. The forcing is
        manufactured automatically as
        :math:`-\nu\Delta u+\gamma u+\nabla p`.
    point_evaluator :
        Optional NumPy evaluator for plotting. It receives coordinates with
        shape ``(dimension, npoints)`` and returns velocity with shape
        ``(dimension, npoints)`` and pressure with shape ``(npoints,)``.
    description, reference :
        Concise provenance suitable for notebook and metadata reporting.

    Notes
    -----
    The convergence runner imposes the exact velocity on the complete boundary
    and compares pressure modulo an additive constant.
    """

    name: str
    dimension: Literal[2, 3]
    viscosity: float
    reaction: float
    exact_solution_factory: MMSExactSolutionFactory = field(repr=False)
    point_evaluator: MMSPointEvaluator | None = field(default=None, repr=False)
    description: str = ""
    reference: str = ""

    def __post_init__(self) -> None:
        if not self.name.strip():
            raise ValueError("name must not be empty")
        if self.dimension not in {2, 3}:
            raise ValueError("dimension must be either 2 or 3")
        if self.viscosity <= 0.0 or not np.isfinite(self.viscosity):
            raise ValueError("viscosity must be positive and finite")
        if self.reaction < 0.0 or not np.isfinite(self.reaction):
            raise ValueError("reaction must be nonnegative and finite")
        if not callable(self.exact_solution_factory):
            raise TypeError("exact_solution_factory must be callable")
        if self.point_evaluator is not None and not callable(self.point_evaluator):
            raise TypeError("point_evaluator must be callable")

    def ufl_solution(self, ufl: Any, domain: Any) -> tuple[Any, Any]:
        """Return exact velocity and pressure UFL expressions on ``domain``."""

        velocity, pressure = self.exact_solution_factory(ufl, domain)
        return velocity, pressure

    def evaluate(self, coordinates: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
        """Evaluate the exact fields at NumPy coordinates for plotting."""

        if self.point_evaluator is None:
            raise NotImplementedError(f"Case {self.name!r} does not define a NumPy point evaluator")
        points = np.asarray(coordinates, dtype=float)
        if points.ndim == 1:
            points = points.reshape(self.dimension, 1)
        if points.ndim != 2 or points.shape[0] < self.dimension:
            raise ValueError(
                f"coordinates must have shape (dimension, npoints); received {points.shape}"
            )
        velocity, pressure = self.point_evaluator(points[: self.dimension])
        velocity_array = np.asarray(velocity, dtype=float)
        pressure_array = np.asarray(pressure, dtype=float)
        expected_velocity_shape = (self.dimension, points.shape[1])
        if velocity_array.shape != expected_velocity_shape:
            raise ValueError(
                "point_evaluator returned velocity with shape "
                f"{velocity_array.shape}; expected {expected_velocity_shape}"
            )
        if pressure_array.shape != (points.shape[1],):
            raise ValueError(
                "point_evaluator returned pressure with shape "
                f"{pressure_array.shape}; expected {(points.shape[1],)}"
            )
        return velocity_array, pressure_array

ufl_solution

ufl_solution(ufl, domain)

Return exact velocity and pressure UFL expressions on domain.

Source code in src/voids/examples/mms/_core.py

def ufl_solution(self, ufl: Any, domain: Any) -> tuple[Any, Any]:
    """Return exact velocity and pressure UFL expressions on ``domain``."""

    velocity, pressure = self.exact_solution_factory(ufl, domain)
    return velocity, pressure

evaluate

evaluate(coordinates)

Evaluate the exact fields at NumPy coordinates for plotting.

Source code in src/voids/examples/mms/_core.py

def evaluate(self, coordinates: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Evaluate the exact fields at NumPy coordinates for plotting."""

    if self.point_evaluator is None:
        raise NotImplementedError(f"Case {self.name!r} does not define a NumPy point evaluator")
    points = np.asarray(coordinates, dtype=float)
    if points.ndim == 1:
        points = points.reshape(self.dimension, 1)
    if points.ndim != 2 or points.shape[0] < self.dimension:
        raise ValueError(
            f"coordinates must have shape (dimension, npoints); received {points.shape}"
        )
    velocity, pressure = self.point_evaluator(points[: self.dimension])
    velocity_array = np.asarray(velocity, dtype=float)
    pressure_array = np.asarray(pressure, dtype=float)
    expected_velocity_shape = (self.dimension, points.shape[1])
    if velocity_array.shape != expected_velocity_shape:
        raise ValueError(
            "point_evaluator returned velocity with shape "
            f"{velocity_array.shape}; expected {expected_velocity_shape}"
        )
    if pressure_array.shape != (points.shape[1],):
        raise ValueError(
            "point_evaluator returned pressure with shape "
            f"{pressure_array.shape}; expected {(points.shape[1],)}"
        )
    return velocity_array, pressure_array

ConvergenceExpectation dataclass

Nominal smooth-solution orders used by an MMS method check.

Source code in src/voids/examples/mms/_core.py

@dataclass(frozen=True, slots=True)
class ConvergenceExpectation:
    """Nominal smooth-solution orders used by an MMS method check."""

    velocity_l2: float
    velocity_h1: float
    pressure_l2: float

    def as_dict(self) -> dict[str, float]:
        """Return expected orders keyed by error metric."""

        return {
            "velocity_l2": self.velocity_l2,
            "velocity_h1": self.velocity_h1,
            "pressure_l2": self.pressure_l2,
        }

as_dict

as_dict()

Return expected orders keyed by error metric.

Source code in src/voids/examples/mms/_core.py

def as_dict(self) -> dict[str, float]:
    """Return expected orders keyed by error metric."""

    return {
        "velocity_l2": self.velocity_l2,
        "velocity_h1": self.velocity_h1,
        "pressure_l2": self.pressure_l2,
    }

MMSConvergenceLevel dataclass

Errors and observed pairwise rates for one structured mesh level.

Source code in src/voids/examples/mms/_core.py

@dataclass(frozen=True, slots=True)
class MMSConvergenceLevel:
    """Errors and observed pairwise rates for one structured mesh level."""

    resolution: int
    h: float
    num_cells: int
    num_dofs: int
    solve_seconds: float
    velocity_l2_error: float
    velocity_h1_error: float
    pressure_l2_error: float
    divergence_l2: float
    rates: dict[str, float] = field(default_factory=dict)

    def errors(self) -> dict[str, float]:
        """Return the four reported errors as a metric mapping."""

        return {
            "velocity_l2": self.velocity_l2_error,
            "velocity_h1": self.velocity_h1_error,
            "pressure_l2": self.pressure_l2_error,
            "divergence_l2": self.divergence_l2,
        }

    def as_dict(self) -> dict[str, int | float | None]:
        """Return a row suitable for a table or CSV writer."""

        row: dict[str, int | float | None] = {
            "resolution": self.resolution,
            "h": self.h,
            "num_cells": self.num_cells,
            "num_dofs": self.num_dofs,
            "solve_seconds": self.solve_seconds,
            "velocity_l2_error": self.velocity_l2_error,
            "velocity_h1_error": self.velocity_h1_error,
            "pressure_l2_error": self.pressure_l2_error,
            "divergence_l2": self.divergence_l2,
        }
        for name in _ERROR_NAMES:
            row[f"{name}_rate"] = self.rates.get(name)
        return row

errors

errors()

Return the four reported errors as a metric mapping.

Source code in src/voids/examples/mms/_core.py

def errors(self) -> dict[str, float]:
    """Return the four reported errors as a metric mapping."""

    return {
        "velocity_l2": self.velocity_l2_error,
        "velocity_h1": self.velocity_h1_error,
        "pressure_l2": self.pressure_l2_error,
        "divergence_l2": self.divergence_l2,
    }

as_dict

as_dict()

Return a row suitable for a table or CSV writer.

Source code in src/voids/examples/mms/_core.py

def as_dict(self) -> dict[str, int | float | None]:
    """Return a row suitable for a table or CSV writer."""

    row: dict[str, int | float | None] = {
        "resolution": self.resolution,
        "h": self.h,
        "num_cells": self.num_cells,
        "num_dofs": self.num_dofs,
        "solve_seconds": self.solve_seconds,
        "velocity_l2_error": self.velocity_l2_error,
        "velocity_h1_error": self.velocity_h1_error,
        "pressure_l2_error": self.pressure_l2_error,
        "divergence_l2": self.divergence_l2,
    }
    for name in _ERROR_NAMES:
        row[f"{name}_rate"] = self.rates.get(name)
    return row

MMSConvergenceResult dataclass

Complete manufactured-solution refinement study.

Source code in src/voids/examples/mms/_core.py

@dataclass(slots=True)
class MMSConvergenceResult:
    """Complete manufactured-solution refinement study."""

    case: BrinkmanMMSCase
    method: MMSMethod
    levels: tuple[MMSConvergenceLevel, ...]
    expected_rates: ConvergenceExpectation
    finest_solution: MMSDiscreteSolution | None = None
    metadata: dict[str, Any] = field(default_factory=dict)

    @property
    def last_rates(self) -> dict[str, float]:
        """Return rates between the two finest mesh levels."""

        return dict(self.levels[-1].rates)

    def as_dicts(self) -> list[dict[str, int | float | None]]:
        """Return refinement rows without requiring pandas."""

        return [level.as_dict() for level in self.levels]

    def assert_expected_rates(self, *, absolute_tolerance: float = 0.35) -> None:
        """Raise if a finest-pair rate is below its nominal smooth rate.

        This is an asymptotic diagnostic, not a proof of convergence. Coarse
        meshes, boundary layers, algebraic solver error, or insufficient
        quadrature can all make a correct discretization fail this check.
        """

        if absolute_tolerance < 0.0 or not np.isfinite(absolute_tolerance):
            raise ValueError("absolute_tolerance must be nonnegative and finite")
        failures: list[str] = []
        for name, expected in self.expected_rates.as_dict().items():
            observed = self.last_rates.get(name, float("nan"))
            threshold = expected - absolute_tolerance
            if not np.isfinite(observed) or observed < threshold:
                failures.append(
                    f"{name}: observed {observed:.3f}, expected at least {threshold:.3f}"
                )
        if failures:
            details = "; ".join(failures)
            raise AssertionError(
                f"{self.case.name}/{self.method} did not meet the expected finest-pair "
                f"rates: {details}"
            )

last_rates property

last_rates

Return rates between the two finest mesh levels.

as_dicts

as_dicts()

Return refinement rows without requiring pandas.

Source code in src/voids/examples/mms/_core.py

def as_dicts(self) -> list[dict[str, int | float | None]]:
    """Return refinement rows without requiring pandas."""

    return [level.as_dict() for level in self.levels]

assert_expected_rates

assert_expected_rates(*, absolute_tolerance=0.35)

Raise if a finest-pair rate is below its nominal smooth rate.

This is an asymptotic diagnostic, not a proof of convergence. Coarse meshes, boundary layers, algebraic solver error, or insufficient quadrature can all make a correct discretization fail this check.

Source code in src/voids/examples/mms/_core.py

def assert_expected_rates(self, *, absolute_tolerance: float = 0.35) -> None:
    """Raise if a finest-pair rate is below its nominal smooth rate.

    This is an asymptotic diagnostic, not a proof of convergence. Coarse
    meshes, boundary layers, algebraic solver error, or insufficient
    quadrature can all make a correct discretization fail this check.
    """

    if absolute_tolerance < 0.0 or not np.isfinite(absolute_tolerance):
        raise ValueError("absolute_tolerance must be nonnegative and finite")
    failures: list[str] = []
    for name, expected in self.expected_rates.as_dict().items():
        observed = self.last_rates.get(name, float("nan"))
        threshold = expected - absolute_tolerance
        if not np.isfinite(observed) or observed < threshold:
            failures.append(
                f"{name}: observed {observed:.3f}, expected at least {threshold:.3f}"
            )
    if failures:
        details = "; ".join(failures)
        raise AssertionError(
            f"{self.case.name}/{self.method} did not meet the expected finest-pair "
            f"rates: {details}"
        )

MMSDiscreteSolution dataclass

Finest-mesh DOLFINx fields retained by a convergence run.

Source code in src/voids/examples/mms/_core.py

@dataclass(slots=True)
class MMSDiscreteSolution:
    """Finest-mesh DOLFINx fields retained by a convergence run."""

    mesh: Any
    velocity: Any
    pressure: Any

MMSPresentationReference dataclass

Exact configuration and reported values for one MMS presentation row.

Source code in src/voids/examples/mms/replication.py

@dataclass(frozen=True, slots=True)
class MMSPresentationReference:
    """Exact configuration and reported values for one MMS presentation row."""

    name: str
    case: BrinkmanMMSCase
    method: MMSMethod
    resolutions: tuple[int, ...]
    facet_law: MMSFacetLaw
    facet_size_mode: MMSFacetSizeMode
    quantities: tuple[ReferenceQuantity, ...]
    source: str
    face_refinement: int = 24

    def __post_init__(self) -> None:
        if not self.name.strip():
            raise ValueError("name must not be empty")
        if len(self.resolutions) < 2:
            raise ValueError("resolutions must contain at least two mesh levels")
        if any(value < 1 for value in self.resolutions):
            raise ValueError("all resolutions must be positive")
        if any(
            current <= previous for previous, current in zip(self.resolutions, self.resolutions[1:])
        ):
            raise ValueError("resolutions must be strictly increasing")
        if not self.quantities:
            raise ValueError("quantities must not be empty")
        if not self.source.strip():
            raise ValueError("source must not be empty")
        if self.face_refinement < 2:
            raise ValueError("face_refinement must be at least 2")
        if self.facet_size_mode not in {
            "cell_diameter",
            "facet_diameter",
            "representative",
        }:
            raise ValueError("facet_size_mode is not supported")

MMSPresentationRun dataclass

Live MMS refinement result paired with its baseline comparison.

Source code in src/voids/examples/mms/replication.py

@dataclass(slots=True)
class MMSPresentationRun:
    """Live MMS refinement result paired with its baseline comparison."""

    reference: MMSPresentationReference
    result: MMSConvergenceResult
    comparison: PresentationComparison

    def assert_matches(self) -> None:
        """Raise unless the live result reproduces every stored target."""

        self.comparison.assert_matches()

assert_matches

assert_matches()

Raise unless the live result reproduces every stored target.

Source code in src/voids/examples/mms/replication.py

def assert_matches(self) -> None:
    """Raise unless the live result reproduces every stored target."""

    self.comparison.assert_matches()

PresentationComparison dataclass

Comparison of one live solve with a supplied presentation baseline.

Source code in src/voids/examples/mms/replication.py

@dataclass(frozen=True, slots=True)
class PresentationComparison:
    """Comparison of one live solve with a supplied presentation baseline."""

    reference_name: str
    quantities: tuple[ReferenceComparison, ...]

    @property
    def passed(self) -> bool:
        """Return whether every reported quantity was reproduced."""

        return all(quantity.passed for quantity in self.quantities)

    def as_dicts(self) -> list[dict[str, str | float | bool]]:
        """Return comparison rows without requiring pandas."""

        return [quantity.as_dict() for quantity in self.quantities]

    def assert_matches(self) -> None:
        """Raise when any reported scalar falls outside its tolerance."""

        failures = [
            (
                f"{quantity.metric}: observed {quantity.observed:.6g}, "
                f"expected {quantity.expected:.6g}, "
                f"absolute error {quantity.absolute_error:.3g}"
            )
            for quantity in self.quantities
            if not quantity.passed
        ]
        if failures:
            raise AssertionError(
                f"{self.reference_name} did not reproduce the supplied baseline: "
                + "; ".join(failures)
            )

passed property

passed

Return whether every reported quantity was reproduced.

as_dicts

as_dicts()

Return comparison rows without requiring pandas.

Source code in src/voids/examples/mms/replication.py

def as_dicts(self) -> list[dict[str, str | float | bool]]:
    """Return comparison rows without requiring pandas."""

    return [quantity.as_dict() for quantity in self.quantities]

assert_matches

assert_matches()

Raise when any reported scalar falls outside its tolerance.

Source code in src/voids/examples/mms/replication.py

def assert_matches(self) -> None:
    """Raise when any reported scalar falls outside its tolerance."""

    failures = [
        (
            f"{quantity.metric}: observed {quantity.observed:.6g}, "
            f"expected {quantity.expected:.6g}, "
            f"absolute error {quantity.absolute_error:.3g}"
        )
        for quantity in self.quantities
        if not quantity.passed
    ]
    if failures:
        raise AssertionError(
            f"{self.reference_name} did not reproduce the supplied baseline: "
            + "; ".join(failures)
        )

ReferenceComparison dataclass

Observed-versus-reported comparison for one scalar quantity.

Source code in src/voids/examples/mms/replication.py

@dataclass(frozen=True, slots=True)
class ReferenceComparison:
    """Observed-versus-reported comparison for one scalar quantity."""

    metric: str
    observed: float
    expected: float
    absolute_tolerance: float
    relative_tolerance: float

    @property
    def absolute_error(self) -> float:
        """Return ``abs(observed - expected)``."""

        return abs(self.observed - self.expected)

    @property
    def relative_error(self) -> float:
        """Return the error relative to the reported magnitude."""

        if self.expected == 0.0:
            return 0.0 if self.observed == 0.0 else float("inf")
        return self.absolute_error / abs(self.expected)

    @property
    def passed(self) -> bool:
        """Return whether the observation lies within the stored tolerance."""

        return bool(
            np.isclose(
                self.observed,
                self.expected,
                atol=self.absolute_tolerance,
                rtol=self.relative_tolerance,
            )
        )

    def as_dict(self) -> dict[str, str | float | bool]:
        """Return a table-ready representation."""

        return {
            "metric": self.metric,
            "observed": self.observed,
            "expected": self.expected,
            "absolute_error": self.absolute_error,
            "relative_error": self.relative_error,
            "passed": self.passed,
        }

absolute_error property

absolute_error

Return abs(observed - expected).

relative_error property

relative_error

Return the error relative to the reported magnitude.

passed property

passed

Return whether the observation lies within the stored tolerance.

as_dict

as_dict()

Return a table-ready representation.

Source code in src/voids/examples/mms/replication.py

def as_dict(self) -> dict[str, str | float | bool]:
    """Return a table-ready representation."""

    return {
        "metric": self.metric,
        "observed": self.observed,
        "expected": self.expected,
        "absolute_error": self.absolute_error,
        "relative_error": self.relative_error,
        "passed": self.passed,
    }

ReferenceQuantity dataclass

One reported scalar value and its replication tolerance.

Source code in src/voids/examples/mms/replication.py

@dataclass(frozen=True, slots=True)
class ReferenceQuantity:
    """One reported scalar value and its replication tolerance."""

    metric: str
    expected: float
    absolute_tolerance: float = 0.0
    relative_tolerance: float = 0.0

    def __post_init__(self) -> None:
        if not self.metric.strip():
            raise ValueError("metric must not be empty")
        if not np.isfinite(self.expected):
            raise ValueError("expected must be finite")
        for name in ("absolute_tolerance", "relative_tolerance"):
            value = float(getattr(self, name))
            if value < 0.0 or not np.isfinite(value):
                raise ValueError(f"{name} must be nonnegative and finite")
        if self.absolute_tolerance == 0.0 and self.relative_tolerance == 0.0:
            raise ValueError("at least one comparison tolerance must be positive")

VugPresentationReference dataclass

Configuration and reported flux for a centered-vug presentation row.

Source code in src/voids/examples/mms/replication.py

@dataclass(frozen=True, slots=True)
class VugPresentationReference:
    """Configuration and reported flux for a centered-vug presentation row."""

    name: str
    benchmark: CenteredVugBenchmark
    method: MMSMethod
    facet_law: USFEMFacetLaw | None
    facet_size_mode: USFEMFacetSizeMode | None
    quantities: tuple[ReferenceQuantity, ...]
    source: str

    def __post_init__(self) -> None:
        if not self.name.strip():
            raise ValueError("name must not be empty")
        if self.benchmark.mesh_representation != "body_fitted":
            raise ValueError("presentation vug references require a body-fitted mesh")
        if not self.quantities:
            raise ValueError("quantities must not be empty")
        if self.method != "taylor_hood" and self.facet_size_mode is None:
            raise ValueError("USFEM vug references require a facet_size_mode")
        if not self.source.strip():
            raise ValueError("source must not be empty")

VugPresentationRun dataclass

Live centered-vug result paired with its baseline comparison.

Source code in src/voids/examples/mms/replication.py

@dataclass(slots=True)
class VugPresentationRun:
    """Live centered-vug result paired with its baseline comparison."""

    reference: VugPresentationReference
    result: FEMSinglePhaseResult
    comparison: PresentationComparison

    def assert_matches(self) -> None:
        """Raise unless the live result reproduces every stored target."""

        self.comparison.assert_matches()

assert_matches

assert_matches()

Raise unless the live result reproduces every stored target.

Source code in src/voids/examples/mms/replication.py

def assert_matches(self) -> None:
    """Raise unless the live result reproduces every stored target."""

    self.comparison.assert_matches()

BodyFittedVugMesh dataclass

DOLFINx mesh and physical tags generated for a centered vug.

Source code in src/voids/examples/mms/vug.py

@dataclass(slots=True)
class BodyFittedVugMesh:
    """DOLFINx mesh and physical tags generated for a centered vug."""

    mesh: Any
    cell_tags: Any
    facet_tags: Any
    physical_groups: dict[str, Any]

CenteredVugBenchmark dataclass

Centered circular/spherical vug benchmark on the unit domain.

The defaults reproduce the documented physical configuration: radius 0.25, matrix drag 1e7, vug drag 1, viscosity 1e-2, and pressure values 1 and -1. mesh_representation="body_fitted" uses Gmsh physical cell tags. The "structured" option classifies coefficient-map cells by their centers and is useful as a portable representation-sensitivity comparison.

resolution is the number of nominal elements per coordinate direction. For body-fitted meshes, the Gmsh target size is sqrt(dimension) / resolution, matching the mesh-size convention used in the reference vug studies.

Source code in src/voids/examples/mms/vug.py

@dataclass(frozen=True, slots=True)
class CenteredVugBenchmark:
    """Centered circular/spherical vug benchmark on the unit domain.

    The defaults reproduce the documented physical configuration: radius
    ``0.25``, matrix drag ``1e7``, vug drag ``1``, viscosity ``1e-2``, and
    pressure values ``1`` and ``-1``. ``mesh_representation="body_fitted"``
    uses Gmsh physical cell tags. The ``"structured"`` option classifies
    coefficient-map cells by their centers and is useful as a portable
    representation-sensitivity comparison.

    ``resolution`` is the number of nominal elements per coordinate direction.
    For body-fitted meshes, the Gmsh target size is
    ``sqrt(dimension) / resolution``, matching the mesh-size convention used in
    the reference vug studies.
    """

    dimension: Literal[2, 3] = 2
    resolution: int = 32
    radius: float = 0.25
    viscosity: float = 1.0e-2
    matrix_drag: float = 1.0e7
    vug_drag: float = 1.0
    pressure_inlet: float = 1.0
    pressure_outlet: float = -1.0
    mesh_representation: VugMeshRepresentation = "body_fitted"
    matrix_effective_viscosity: float | None = None
    vug_effective_viscosity: float | None = None

    def __post_init__(self) -> None:
        if self.dimension not in {2, 3}:
            raise ValueError("dimension must be either 2 or 3")
        if self.resolution < 2:
            raise ValueError("resolution must be at least 2")
        if not 0.0 <= self.radius < 0.5:
            raise ValueError("radius must lie between 0 (inclusive) and 0.5")
        for name in ("viscosity", "matrix_drag"):
            value = float(getattr(self, name))
            if value <= 0.0 or not np.isfinite(value):
                raise ValueError(f"{name} must be positive and finite")
        if self.vug_drag < 0.0 or not np.isfinite(self.vug_drag):
            raise ValueError("vug_drag must be non-negative and finite")
        for name in ("matrix_effective_viscosity", "vug_effective_viscosity"):
            value = getattr(self, name)
            if value is not None and (value <= 0.0 or not np.isfinite(value)):
                raise ValueError(f"{name} must be positive and finite when provided")
        if (
            not np.isfinite(self.pressure_inlet)
            or not np.isfinite(self.pressure_outlet)
            or self.pressure_inlet <= self.pressure_outlet
        ):
            raise ValueError("pressure_inlet must be finite and greater than pressure_outlet")
        if self.mesh_representation not in {"body_fitted", "structured"}:
            raise ValueError("mesh_representation must be either 'body_fitted' or 'structured'")

    @property
    def target_mesh_size(self) -> float:
        """Return the Gmsh target element diameter."""

        return float(np.sqrt(self.dimension) / self.resolution)

    @property
    def matrix_nu(self) -> float:
        """Return the matrix coefficient multiplying the Brinkman gradient term."""

        if self.matrix_effective_viscosity is None:
            return float(self.viscosity)
        return float(self.matrix_effective_viscosity)

    @property
    def vug_nu(self) -> float:
        """Return the vug coefficient multiplying the Brinkman gradient term."""

        if self.vug_effective_viscosity is None:
            return float(self.viscosity)
        return float(self.vug_effective_viscosity)

    def vug_mask(self) -> np.ndarray:
        """Return the structured cell-center classification of the vug."""

        coordinates = (np.arange(self.resolution, dtype=float) + 0.5) / self.resolution
        grids = np.meshgrid(
            *((coordinates,) * self.dimension),
            indexing="ij",
        )
        radius_squared = sum((grid - 0.5) ** 2 for grid in grids)
        return np.asarray(radius_squared <= self.radius**2, dtype=bool)

    @property
    def represented_fraction(self) -> float:
        """Return the structured cell-volume fraction classified as vug."""

        return float(np.mean(self.vug_mask()))

    @property
    def analytic_fraction(self) -> float:
        """Return the exact circular area or spherical volume fraction."""

        if self.dimension == 2:
            return float(np.pi * self.radius**2)
        return float((4.0 / 3.0) * np.pi * self.radius**3)

    def make_problem(self) -> FEMMapProblem:
        """Build the structured constant-porosity coefficient-map problem."""

        mask = self.vug_mask()
        drag = np.where(mask, self.vug_drag, self.matrix_drag)
        permeability = np.divide(
            self.viscosity,
            drag,
            out=np.full_like(drag, np.inf, dtype=float),
            where=drag > 0.0,
        )
        if not np.all(np.isfinite(permeability)):
            raise ValueError(
                "structured vug benchmarks require positive vug_drag; "
                "zero drag is supported only by the body-fitted formulation"
            )
        effective_viscosity = np.where(mask, self.vug_nu, self.matrix_nu)
        porosity = self.viscosity / effective_viscosity
        cell_size = (1.0 / self.resolution,) * self.dimension
        return FEMMapProblem(
            permeability_map=PermeabilityMap(
                permeability,
                cell_size=cell_size,
            ),
            porosity_map=PorosityMap(
                porosity,
                cell_size=cell_size,
            ),
            viscosity=self.viscosity,
        )

target_mesh_size property

target_mesh_size

Return the Gmsh target element diameter.

matrix_nu property

matrix_nu

Return the matrix coefficient multiplying the Brinkman gradient term.

vug_nu property

vug_nu

Return the vug coefficient multiplying the Brinkman gradient term.

represented_fraction property

represented_fraction

Return the structured cell-volume fraction classified as vug.

analytic_fraction property

analytic_fraction

Return the exact circular area or spherical volume fraction.

vug_mask

vug_mask()

Return the structured cell-center classification of the vug.

Source code in src/voids/examples/mms/vug.py

def vug_mask(self) -> np.ndarray:
    """Return the structured cell-center classification of the vug."""

    coordinates = (np.arange(self.resolution, dtype=float) + 0.5) / self.resolution
    grids = np.meshgrid(
        *((coordinates,) * self.dimension),
        indexing="ij",
    )
    radius_squared = sum((grid - 0.5) ** 2 for grid in grids)
    return np.asarray(radius_squared <= self.radius**2, dtype=bool)

make_problem

make_problem()

Build the structured constant-porosity coefficient-map problem.

Source code in src/voids/examples/mms/vug.py

def make_problem(self) -> FEMMapProblem:
    """Build the structured constant-porosity coefficient-map problem."""

    mask = self.vug_mask()
    drag = np.where(mask, self.vug_drag, self.matrix_drag)
    permeability = np.divide(
        self.viscosity,
        drag,
        out=np.full_like(drag, np.inf, dtype=float),
        where=drag > 0.0,
    )
    if not np.all(np.isfinite(permeability)):
        raise ValueError(
            "structured vug benchmarks require positive vug_drag; "
            "zero drag is supported only by the body-fitted formulation"
        )
    effective_viscosity = np.where(mask, self.vug_nu, self.matrix_nu)
    porosity = self.viscosity / effective_viscosity
    cell_size = (1.0 / self.resolution,) * self.dimension
    return FEMMapProblem(
        permeability_map=PermeabilityMap(
            permeability,
            cell_size=cell_size,
        ),
        porosity_map=PorosityMap(
            porosity,
            cell_size=cell_size,
        ),
        viscosity=self.viscosity,
    )

CenteredVugFlowCase2D dataclass

Physical 2D centered-vug case solved on a body-fitted unit-square mesh.

The continuum domain size is derived from image_shape and voxel_size_m. Gmsh coordinates are normalized to the unit square for numerical conditioning; the solver nondimensionalizes the equations and converts velocity, pressure, flux, and permeability back to SI units.

The Darcy--Brinkman branch follows the layered-domain model: its reaction coefficient is mu / K_matrix in the matrix and exactly zero in the vug. vug_permeability_m2 is used only by the Darcy--Darcy branch as a finite high-permeability closure. It is not an intrinsic permeability measurement of an open cavity.

Source code in src/voids/examples/mms/vug_flow.py

@dataclass(frozen=True, slots=True)
class CenteredVugFlowCase2D:
    """Physical 2D centered-vug case solved on a body-fitted unit-square mesh.

    The continuum domain size is derived from ``image_shape`` and
    ``voxel_size_m``. Gmsh coordinates are normalized to the unit square for
    numerical conditioning; the solver nondimensionalizes the equations and
    converts velocity, pressure, flux, and permeability back to SI units.

    The Darcy--Brinkman branch follows the layered-domain model: its reaction
    coefficient is ``mu / K_matrix`` in the matrix and exactly zero in the
    vug. ``vug_permeability_m2`` is used only by the Darcy--Darcy branch as a
    finite high-permeability closure. It is not an intrinsic permeability
    measurement of an open cavity.
    """

    area_fraction: float
    image_shape: tuple[int, int] = (500, 500)
    voxel_size_m: float = 15.0e-6
    matrix_porosity: float = 0.2
    matrix_permeability_md: float = 200.0
    vug_permeability_m2: float = 1.0e-8
    dynamic_viscosity_pa_s: float = 1.0e-3
    pressure_inlet_pa: float = 1.0
    pressure_outlet_pa: float = 0.0
    mesh_resolution: int = 100

    def __post_init__(self) -> None:
        if len(self.image_shape) != 2 or self.image_shape[0] != self.image_shape[1]:
            raise ValueError("image_shape must describe a square 2D image")
        if self.image_shape[0] < 2:
            raise ValueError("image_shape values must be at least 2")
        if self.voxel_size_m <= 0.0 or not np.isfinite(self.voxel_size_m):
            raise ValueError("voxel_size_m must be positive and finite")
        maximum_fraction = float(np.pi / 4.0)
        if (
            self.area_fraction < 0.0
            or self.area_fraction >= maximum_fraction
            or not np.isfinite(self.area_fraction)
        ):
            raise ValueError("area_fraction must lie between 0 (inclusive) and pi/4 (exclusive)")
        for name in (
            "matrix_porosity",
            "matrix_permeability_md",
            "vug_permeability_m2",
            "dynamic_viscosity_pa_s",
        ):
            value = float(getattr(self, name))
            if value <= 0.0 or not np.isfinite(value):
                raise ValueError(f"{name} must be positive and finite")
        if self.matrix_porosity > 1.0:
            raise ValueError("matrix_porosity must not exceed 1")
        if (
            not np.isfinite(self.pressure_inlet_pa)
            or not np.isfinite(self.pressure_outlet_pa)
            or self.pressure_inlet_pa <= self.pressure_outlet_pa
        ):
            raise ValueError("pressure_inlet_pa must be finite and greater than pressure_outlet_pa")
        if self.mesh_resolution < 8:
            raise ValueError("mesh_resolution must be at least 8")

    @property
    def side_length_m(self) -> float:
        """Return the square side length derived from the image metadata."""

        return float(self.image_shape[0] * self.voxel_size_m)

    @property
    def radius_fraction(self) -> float:
        """Return the circle radius divided by the square side length."""

        return float(np.sqrt(self.area_fraction / np.pi))

    @property
    def radius_m(self) -> float:
        """Return the physical centered-vug radius."""

        return float(self.radius_fraction * self.side_length_m)

    @property
    def matrix_permeability_m2(self) -> float:
        """Return the matrix permeability converted from mD to square metres."""

        return float(self.matrix_permeability_md * M2_PER_MILLIDARCY)

    @property
    def matrix_screening_length_m(self) -> float:
        """Return the matrix Brinkman screening length ``sqrt(K / phi)``."""

        return float(np.sqrt(self.matrix_permeability_m2 / self.matrix_porosity))

    @property
    def permeability_contrast(self) -> float:
        """Return the Darcy--Darcy closure contrast ``K_vug / K_matrix``."""

        return float(self.vug_permeability_m2 / self.matrix_permeability_m2)

    @property
    def pressure_drop_pa(self) -> float:
        """Return the applied pressure drop."""

        return float(self.pressure_inlet_pa - self.pressure_outlet_pa)

    @property
    def base_target_mesh_size_fraction(self) -> float:
        """Return the far-field target diameter divided by the side length."""

        return float(np.sqrt(2.0) / self.mesh_resolution)

    def make_benchmark(
        self,
        *,
        model: CenteredVugFlowModel = "darcy_brinkman",
    ) -> CenteredVugBenchmark:
        """Return the model-specific nondimensional mesh/coefficient definition."""

        if model not in {"darcy_brinkman", "darcy_darcy"}:
            raise ValueError("model must be either 'darcy_brinkman' or 'darcy_darcy'")
        darcy_number = self.matrix_permeability_m2 / self.side_length_m**2
        vug_drag = (
            0.0
            if model == "darcy_brinkman"
            else self.matrix_permeability_m2 / self.vug_permeability_m2
        )
        return CenteredVugBenchmark(
            dimension=2,
            resolution=self.mesh_resolution,
            radius=self.radius_fraction,
            viscosity=darcy_number,
            matrix_drag=1.0,
            vug_drag=vug_drag,
            pressure_inlet=1.0,
            pressure_outlet=0.0,
            mesh_representation="body_fitted",
            matrix_effective_viscosity=darcy_number / self.matrix_porosity,
            vug_effective_viscosity=darcy_number,
        )

side_length_m property

side_length_m

Return the square side length derived from the image metadata.

radius_fraction property

radius_fraction

Return the circle radius divided by the square side length.

radius_m property

radius_m

Return the physical centered-vug radius.

matrix_permeability_m2 property

matrix_permeability_m2

Return the matrix permeability converted from mD to square metres.

matrix_screening_length_m property

matrix_screening_length_m

Return the matrix Brinkman screening length sqrt(K / phi).

permeability_contrast property

permeability_contrast

Return the Darcy--Darcy closure contrast K_vug / K_matrix.

pressure_drop_pa property

pressure_drop_pa

Return the applied pressure drop.

base_target_mesh_size_fraction property

base_target_mesh_size_fraction

Return the far-field target diameter divided by the side length.

make_benchmark

make_benchmark(*, model='darcy_brinkman')

Return the model-specific nondimensional mesh/coefficient definition.

Source code in src/voids/examples/mms/vug_flow.py

def make_benchmark(
    self,
    *,
    model: CenteredVugFlowModel = "darcy_brinkman",
) -> CenteredVugBenchmark:
    """Return the model-specific nondimensional mesh/coefficient definition."""

    if model not in {"darcy_brinkman", "darcy_darcy"}:
        raise ValueError("model must be either 'darcy_brinkman' or 'darcy_darcy'")
    darcy_number = self.matrix_permeability_m2 / self.side_length_m**2
    vug_drag = (
        0.0
        if model == "darcy_brinkman"
        else self.matrix_permeability_m2 / self.vug_permeability_m2
    )
    return CenteredVugBenchmark(
        dimension=2,
        resolution=self.mesh_resolution,
        radius=self.radius_fraction,
        viscosity=darcy_number,
        matrix_drag=1.0,
        vug_drag=vug_drag,
        pressure_inlet=1.0,
        pressure_outlet=0.0,
        mesh_representation="body_fitted",
        matrix_effective_viscosity=darcy_number / self.matrix_porosity,
        vug_effective_viscosity=darcy_number,
    )

available_mms_methods

available_mms_methods()

Return the finite-element formulations supported by the MMS runner.

Source code in src/voids/examples/mms/_runner.py

def available_mms_methods() -> tuple[MMSMethod, ...]:
    """Return the finite-element formulations supported by the MMS runner."""

    return tuple(_METHOD_SPECS)

face3d_pressure_jump_coefficient

face3d_pressure_jump_coefficient(
    *, viscosity, reaction, resolution, face_refinement=24
)

Compute the scalar triangular-face subscale coefficient.

The reference problem is solved with continuous piecewise-linear finite elements on a uniformly refined right triangle. resolution is the unit-cube subdivision count, giving representative physical face diameter sqrt(2) / resolution.

Source code in src/voids/examples/mms/_runner.py

def face3d_pressure_jump_coefficient(
    *,
    viscosity: float,
    reaction: float,
    resolution: int,
    face_refinement: int = 24,
) -> float:
    """Compute the scalar triangular-face subscale coefficient.

    The reference problem is solved with continuous piecewise-linear finite
    elements on a uniformly refined right triangle. ``resolution`` is the
    unit-cube subdivision count, giving representative physical face diameter
    ``sqrt(2) / resolution``.
    """

    if viscosity <= 0.0 or not np.isfinite(viscosity):
        raise ValueError("viscosity must be positive and finite")
    if reaction < 0.0 or not np.isfinite(reaction):
        raise ValueError("reaction must be nonnegative and finite")
    if resolution < 1:
        raise ValueError("resolution must be positive")
    if face_refinement < 2:
        raise ValueError("face_refinement must be at least 2")
    scale = 1.0 / resolution
    face_diameter = np.sqrt(2.0) * scale
    alpha_squared = reaction * scale * scale / viscosity
    average = _reference_face_average(alpha_squared, face_refinement)
    return float(scale * scale / (viscosity * face_diameter) * average)

observed_rate

observed_rate(
    previous_h, previous_error, current_h, current_error
)

Return the two-level observed order log(e0/e1) / log(h0/h1).

Source code in src/voids/examples/mms/_runner.py

def observed_rate(
    previous_h: float,
    previous_error: float,
    current_h: float,
    current_error: float,
) -> float:
    """Return the two-level observed order ``log(e0/e1) / log(h0/h1)``."""

    values = (previous_h, previous_error, current_h, current_error)
    if any(value <= 0.0 or not np.isfinite(value) for value in values):
        return float("nan")
    if previous_h <= current_h:
        raise ValueError("previous_h must be greater than current_h")
    return float(np.log(previous_error / current_error) / np.log(previous_h / current_h))

run_mms_convergence

run_mms_convergence(
    case,
    *,
    method="taylor_hood",
    resolutions=(4, 8, 16),
    options=None,
    tau_factor=1.0,
    tau_gamma_cap=None,
    m_t=1.0 / 3.0,
    alpha_edge=1.0,
    facet_law="auto",
    facet_size_mode="cell_diameter",
    face_refinement=24,
    keep_solution=True,
    callback=None,
)

Run a structured-mesh Brinkman MMS refinement study.

Parameters:

Name	Type	Description	Default
case	BrinkmanMMSCase	Exact Brinkman solution. The body force is derived automatically.	required
method	MMSMethod	"taylor_hood" for CG2 x CG1, "usfem_p1dg0" for CG1 x DG0, or "usfem_p1dg1" for CG1 x DG1.	'taylor_hood'
resolutions	Sequence[int]	Strictly increasing numbers of subdivisions per coordinate direction. The reported refinement parameter is h = 1 / resolution.	(4, 8, 16)
options	FEniCSSolverOptions | None	Linear solver controls shared by all levels.	None
tau_factor	float	USFEM stabilization controls. facet_size_mode="cell_diameter" preserves the generic solver convention, "facet_diameter" uses physical edge length in 2D, and "representative" uses 1 / n in 2D or sqrt(2) / n in 3D. Set tau_factor=0 to disable the cell term. tau_gamma_cap optionally bounds gamma * tau_K. They are ignored by Taylor-Hood.	1.0
tau_gamma_cap	float	USFEM stabilization controls. facet_size_mode="cell_diameter" preserves the generic solver convention, "facet_diameter" uses physical edge length in 2D, and "representative" uses 1 / n in 2D or sqrt(2) / n in 3D. Set tau_factor=0 to disable the cell term. tau_gamma_cap optionally bounds gamma * tau_K. They are ignored by Taylor-Hood.	1.0
m_t	float	USFEM stabilization controls. facet_size_mode="cell_diameter" preserves the generic solver convention, "facet_diameter" uses physical edge length in 2D, and "representative" uses 1 / n in 2D or sqrt(2) / n in 3D. Set tau_factor=0 to disable the cell term. tau_gamma_cap optionally bounds gamma * tau_K. They are ignored by Taylor-Hood.	1.0
alpha_edge	float	USFEM stabilization controls. facet_size_mode="cell_diameter" preserves the generic solver convention, "facet_diameter" uses physical edge length in 2D, and "representative" uses 1 / n in 2D or sqrt(2) / n in 3D. Set tau_factor=0 to disable the cell term. tau_gamma_cap optionally bounds gamma * tau_K. They are ignored by Taylor-Hood.	1.0
facet_law	float	USFEM stabilization controls. facet_size_mode="cell_diameter" preserves the generic solver convention, "facet_diameter" uses physical edge length in 2D, and "representative" uses 1 / n in 2D or sqrt(2) / n in 3D. Set tau_factor=0 to disable the cell term. tau_gamma_cap optionally bounds gamma * tau_K. They are ignored by Taylor-Hood.	1.0
facet_size_mode	float	USFEM stabilization controls. facet_size_mode="cell_diameter" preserves the generic solver convention, "facet_diameter" uses physical edge length in 2D, and "representative" uses 1 / n in 2D or sqrt(2) / n in 3D. Set tau_factor=0 to disable the cell term. tau_gamma_cap optionally bounds gamma * tau_K. They are ignored by Taylor-Hood.	1.0
face_refinement	float	USFEM stabilization controls. facet_size_mode="cell_diameter" preserves the generic solver convention, "facet_diameter" uses physical edge length in 2D, and "representative" uses 1 / n in 2D or sqrt(2) / n in 3D. Set tau_factor=0 to disable the cell term. tau_gamma_cap optionally bounds gamma * tau_K. They are ignored by Taylor-Hood.	1.0
keep_solution	bool	Retain the finest DOLFINx velocity and pressure fields for plotting.	True
callback	Callable[[MMSConvergenceLevel], None] | None	Optional callable invoked after each completed level.	None

Returns:

Type	Description
MMSConvergenceResult	Errors, pairwise observed rates, nominal expected rates, solver metadata, and optionally the finest fields.

Source code in src/voids/examples/mms/_runner.py

def run_mms_convergence(
    case: BrinkmanMMSCase,
    *,
    method: MMSMethod = "taylor_hood",
    resolutions: Sequence[int] = (4, 8, 16),
    options: FEniCSSolverOptions | None = None,
    tau_factor: float = 1.0,
    tau_gamma_cap: float | None = None,
    m_t: float = 1.0 / 3.0,
    alpha_edge: float = 1.0,
    facet_law: MMSFacetLaw = "auto",
    facet_size_mode: MMSFacetSizeMode = "cell_diameter",
    face_refinement: int = 24,
    keep_solution: bool = True,
    callback: Callable[[MMSConvergenceLevel], None] | None = None,
) -> MMSConvergenceResult:
    """Run a structured-mesh Brinkman MMS refinement study.

    Parameters
    ----------
    case :
        Exact Brinkman solution. The body force is derived automatically.
    method :
        ``"taylor_hood"`` for CG2 x CG1, ``"usfem_p1dg0"`` for CG1 x
        DG0, or ``"usfem_p1dg1"`` for CG1 x DG1.
    resolutions :
        Strictly increasing numbers of subdivisions per coordinate direction.
        The reported refinement parameter is ``h = 1 / resolution``.
    options :
        Linear solver controls shared by all levels.
    tau_factor, tau_gamma_cap, m_t, alpha_edge, facet_law, facet_size_mode, face_refinement :
        USFEM stabilization controls. ``facet_size_mode="cell_diameter"``
        preserves the generic solver convention, ``"facet_diameter"`` uses
        physical edge length in 2D, and ``"representative"`` uses ``1 / n`` in
        2D or ``sqrt(2) / n`` in 3D. Set ``tau_factor=0`` to disable the cell
        term. ``tau_gamma_cap`` optionally bounds ``gamma * tau_K``. They are
        ignored by Taylor-Hood.
    keep_solution :
        Retain the finest DOLFINx velocity and pressure fields for plotting.
    callback :
        Optional callable invoked after each completed level.

    Returns
    -------
    MMSConvergenceResult
        Errors, pairwise observed rates, nominal expected rates, solver
        metadata, and optionally the finest fields.
    """

    if method not in _METHOD_SPECS:
        supported = ", ".join(_METHOD_SPECS)
        raise ValueError(f"method must be one of {supported}")
    normalized_resolutions = _validate_resolutions(resolutions)
    if tau_factor < 0.0 or not np.isfinite(tau_factor):
        raise ValueError("tau_factor must be nonnegative and finite")
    if tau_gamma_cap is not None and (
        tau_gamma_cap <= 0.0 or tau_gamma_cap >= 1.0 or not np.isfinite(tau_gamma_cap)
    ):
        raise ValueError("tau_gamma_cap must satisfy 0 < tau_gamma_cap < 1")
    if m_t <= 0.0 or not np.isfinite(m_t):
        raise ValueError("m_t must be positive and finite")
    if alpha_edge <= 0.0 or not np.isfinite(alpha_edge):
        raise ValueError("alpha_edge must be positive and finite")
    if facet_law not in {
        "auto",
        "classic",
        "reaction_diffusion",
        "shifted",
        "face3d",
    }:
        raise ValueError(
            "facet_law must be one of 'auto', 'classic', 'reaction_diffusion', "
            "'shifted', or 'face3d'"
        )
    if facet_size_mode not in {
        "cell_diameter",
        "facet_diameter",
        "representative",
    }:
        raise ValueError(
            "facet_size_mode must be one of 'cell_diameter', 'facet_diameter', or 'representative'"
        )
    if face_refinement < 2:
        raise ValueError("face_refinement must be at least 2")
    resolved_facet_law: MMSFacetLaw = (
        "face3d" if facet_law == "auto" and case.dimension == 3 else facet_law
    )
    if resolved_facet_law == "auto":
        resolved_facet_law = "reaction_diffusion"
    if (
        method != "taylor_hood"
        and resolved_facet_law in {"reaction_diffusion", "face3d"}
        and not np.isclose(alpha_edge, 1.0)
    ):
        warnings.warn(
            f"alpha_edge is ignored by the parameter-free {resolved_facet_law} facet law",
            RuntimeWarning,
            stacklevel=2,
        )

    levels: list[MMSConvergenceLevel] = []
    finest_solution: MMSDiscreteSolution | None = None
    metadata: dict[str, Any] = {}
    for resolution in normalized_resolutions:
        level, discrete_solution, metadata = _solve_level(
            case,
            method,
            resolution,
            options=options,
            tau_factor=tau_factor,
            tau_gamma_cap=tau_gamma_cap,
            m_t=m_t,
            alpha_edge=alpha_edge,
            facet_law=resolved_facet_law,
            facet_size_mode=facet_size_mode,
            face_refinement=face_refinement,
        )
        if levels:
            level = replace(level, rates=_rates(levels[-1], level))
        levels.append(level)
        finest_solution = discrete_solution
        if callback is not None:
            callback(level)

    return MMSConvergenceResult(
        case=case,
        method=method,
        levels=tuple(levels),
        expected_rates=_METHOD_SPECS[method].expected_rates,
        finest_solution=finest_solution if keep_solution else None,
        metadata=metadata,
    )

boundary_layer_case_2d

boundary_layer_case_2d(*, viscosity=0.01, reaction=1.0)

Return the two-dimensional Brinkman boundary-layer MMS case.

The unit-square exact solution is

.. math::

u_1 &= y - \frac{\exp((y-1)/\nu)-\exp(-1/\nu)} {1-\exp(-1/\nu)},\\ u_2 &= x - \frac{\exp((x-1)/\nu)-\exp(-1/\nu)} {1-\exp(-1/\nu)},\\ p &= x-y.

It is exactly divergence-free and develops boundary layers at the top and right boundaries as viscosity decreases.

Source code in src/voids/examples/mms/cases_2d.py

def boundary_layer_case_2d(
    *,
    viscosity: float = 1.0e-2,
    reaction: float = 1.0,
) -> BrinkmanMMSCase:
    r"""Return the two-dimensional Brinkman boundary-layer MMS case.

    The unit-square exact solution is

    .. math::

       u_1 &= y -
       \\frac{\\exp((y-1)/\\nu)-\\exp(-1/\\nu)}
            {1-\\exp(-1/\\nu)},\\\\
       u_2 &= x -
       \\frac{\\exp((x-1)/\\nu)-\\exp(-1/\\nu)}
            {1-\\exp(-1/\\nu)},\\\\
       p &= x-y.

    It is exactly divergence-free and develops boundary layers at the top and
    right boundaries as ``viscosity`` decreases.
    """

    nu = float(viscosity)
    gamma = float(reaction)

    def exact_solution_factory(ufl, domain):
        x = ufl.SpatialCoordinate(domain)
        exponential_floor = ufl.exp(-1.0 / nu)
        denominator = 1.0 - exponential_floor
        u_1 = x[1] - (ufl.exp((x[1] - 1.0) / nu) - exponential_floor) / denominator
        u_2 = x[0] - (ufl.exp((x[0] - 1.0) / nu) - exponential_floor) / denominator
        return ufl.as_vector((u_1, u_2)), x[0] - x[1]

    def point_evaluator(points: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
        x = points[0]
        y = points[1]
        exponential_floor = np.exp(-1.0 / nu)
        denominator = -np.expm1(-1.0 / nu)
        u_1 = y - (np.exp((y - 1.0) / nu) - exponential_floor) / denominator
        u_2 = x - (np.exp((x - 1.0) / nu) - exponential_floor) / denominator
        return np.vstack((u_1, u_2)), x - y

    return BrinkmanMMSCase(
        name="brinkman_boundary_layer_2d",
        dimension=2,
        viscosity=nu,
        reaction=gamma,
        exact_solution_factory=exact_solution_factory,
        point_evaluator=point_evaluator,
        description=(
            "Divergence-free unit-square Brinkman solution with exponential layers at y=1 and x=1."
        ),
        reference=(
            "Barrenechea and Valentin (2002), Numerische Mathematik 92, "
            "653-677, doi:10.1007/s002110100371."
        ),
    )

bubble_case_3d

bubble_case_3d(*, viscosity=0.01, reaction=1.0)

Return a smooth divergence-free three-dimensional Brinkman MMS case.

The velocity is constructed from mixed first derivatives of a polynomial boundary bubble. It therefore vanishes on the unit-cube boundary and is divergence-free by cancellation of mixed derivatives. The pressure is sin(2*pi*x) * sin(pi*y) * sin(pi*z).

Source code in src/voids/examples/mms/cases_3d.py

def bubble_case_3d(
    *,
    viscosity: float = 1.0e-2,
    reaction: float = 1.0,
) -> BrinkmanMMSCase:
    """Return a smooth divergence-free three-dimensional Brinkman MMS case.

    The velocity is constructed from mixed first derivatives of a polynomial
    boundary bubble. It therefore vanishes on the unit-cube boundary and is
    divergence-free by cancellation of mixed derivatives. The pressure is
    ``sin(2*pi*x) * sin(pi*y) * sin(pi*z)``.
    """

    nu = float(viscosity)
    gamma = float(reaction)

    def exact_solution_factory(ufl, domain):
        x = ufl.SpatialCoordinate(domain)
        bubble = (
            x[0] ** 2
            * (1.0 - x[0]) ** 2
            * x[1] ** 2
            * (1.0 - x[1]) ** 2
            * x[2] ** 2
            * (1.0 - x[2]) ** 2
        )
        phi = 32.0 * bubble * (1.0 + x[0] + 2.0 * x[1] + 3.0 * x[2])
        phi_gradient = ufl.grad(phi)
        phi_x = phi_gradient[0]
        phi_y = phi_gradient[1]
        phi_z = phi_gradient[2]
        velocity = ufl.as_vector(
            (
                3.0 * phi_y - 2.0 * phi_z,
                phi_z - 3.0 * phi_x,
                2.0 * phi_x - phi_y,
            )
        )
        pressure = ufl.sin(2.0 * np.pi * x[0]) * ufl.sin(np.pi * x[1]) * ufl.sin(np.pi * x[2])
        return velocity, pressure

    def point_evaluator(points: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
        x, y, z = points

        def bubble_1d(value: np.ndarray) -> np.ndarray:
            return value**2 * (1.0 - value) ** 2

        def bubble_1d_derivative(value: np.ndarray) -> np.ndarray:
            return 2.0 * value * (1.0 - value) * (1.0 - 2.0 * value)

        b_x, b_y, b_z = bubble_1d(x), bubble_1d(y), bubble_1d(z)
        db_x = bubble_1d_derivative(x)
        db_y = bubble_1d_derivative(y)
        db_z = bubble_1d_derivative(z)
        bubble = b_x * b_y * b_z
        linear = 1.0 + x + 2.0 * y + 3.0 * z
        phi_x = 32.0 * (db_x * b_y * b_z * linear + bubble)
        phi_y = 32.0 * (b_x * db_y * b_z * linear + 2.0 * bubble)
        phi_z = 32.0 * (b_x * b_y * db_z * linear + 3.0 * bubble)
        velocity = np.vstack(
            (
                3.0 * phi_y - 2.0 * phi_z,
                phi_z - 3.0 * phi_x,
                2.0 * phi_x - phi_y,
            )
        )
        pressure = np.sin(2.0 * np.pi * x) * np.sin(np.pi * y) * np.sin(np.pi * z)
        return velocity, pressure

    return BrinkmanMMSCase(
        name="brinkman_polynomial_bubble_3d",
        dimension=3,
        viscosity=nu,
        reaction=gamma,
        exact_solution_factory=exact_solution_factory,
        point_evaluator=point_evaluator,
        description=(
            "Smooth divergence-free unit-cube Brinkman solution generated from "
            "a polynomial boundary bubble."
        ),
        reference=(
            "Manufactured case used in the voids 3D USFEM verification study; "
            "forcing is generated from the documented strong residual."
        ),
    )

compare_mms_with_presentation

compare_mms_with_presentation(result, reference)

Compare a live MMS result with the exact supplied configuration and values.

Source code in src/voids/examples/mms/replication.py

def compare_mms_with_presentation(
    result: MMSConvergenceResult,
    reference: str | MMSPresentationReference,
) -> PresentationComparison:
    """Compare a live MMS result with the exact supplied configuration and values."""

    resolved = _resolve_mms_reference(reference)
    if result.method != resolved.method:
        raise ValueError(
            f"method mismatch: result uses {result.method}, reference uses {resolved.method}"
        )
    if result.case.name != resolved.case.name:
        raise ValueError(
            f"case mismatch: result uses {result.case.name}, reference uses {resolved.case.name}"
        )
    if not np.isclose(result.case.viscosity, resolved.case.viscosity) or not np.isclose(
        result.case.reaction,
        resolved.case.reaction,
    ):
        raise ValueError("case coefficient mismatch with presentation reference")
    result_resolutions = tuple(level.resolution for level in result.levels)
    if result_resolutions[-2:] != resolved.resolutions[-2:]:
        raise ValueError(
            "the result must end with the presentation's finest mesh pair "
            f"{resolved.resolutions[-2:]}; received {result_resolutions[-2:]}"
        )
    if result.method != "taylor_hood" and result.metadata.get("facet_law") != resolved.facet_law:
        raise ValueError(
            "facet-law mismatch: result uses "
            f"{result.metadata.get('facet_law')}, reference uses {resolved.facet_law}"
        )
    if (
        result.method != "taylor_hood"
        and result.metadata.get("facet_size_mode") != resolved.facet_size_mode
    ):
        raise ValueError(
            "facet-size mismatch: result uses "
            f"{result.metadata.get('facet_size_mode')}, "
            f"reference uses {resolved.facet_size_mode}"
        )

    finest = result.levels[-1]
    observed = {
        "velocity_l2_error": finest.velocity_l2_error,
        "velocity_h1_error": finest.velocity_h1_error,
        "pressure_l2_error": finest.pressure_l2_error,
        "divergence_l2": finest.divergence_l2,
        **{f"{name}_rate": value for name, value in result.last_rates.items()},
    }
    return _compare_quantities(resolved.name, resolved.quantities, observed)

compare_vug_with_presentation

compare_vug_with_presentation(result, reference)

Compare a live centered-vug result with the supplied report-scale values.

Source code in src/voids/examples/mms/replication.py

def compare_vug_with_presentation(
    result: FEMSinglePhaseResult,
    reference: str | VugPresentationReference,
) -> PresentationComparison:
    """Compare a live centered-vug result with the supplied report-scale values."""

    resolved = _resolve_vug_reference(reference)
    metadata = result.metadata
    for name, expected in (
        ("dimension", resolved.benchmark.dimension),
        ("resolution", resolved.benchmark.resolution),
        ("radius", resolved.benchmark.radius),
        ("matrix_drag", resolved.benchmark.matrix_drag),
        ("vug_drag", resolved.benchmark.vug_drag),
    ):
        actual = metadata.get(name)
        if actual is None or not np.isclose(
            float(actual),
            float(expected),
            rtol=1.0e-12,
            atol=0.0,
        ):
            raise ValueError(f"vug configuration mismatch for {name}: {actual!r} != {expected!r}")
    if resolved.method == "taylor_hood":
        if result.formulation != "brinkman_taylor_hood_p2p1":
            raise ValueError("vug formulation does not match the Taylor-Hood reference")
    else:
        if result.formulation != "brinkman_usfem_p1dg1":
            raise ValueError("vug formulation does not match the P1/DG1 reference")
        if metadata.get("facet_law") != resolved.facet_law:
            raise ValueError("vug facet law does not match the presentation reference")
        if metadata.get("facet_size_mode") != resolved.facet_size_mode:
            raise ValueError("vug facet size does not match the presentation reference")

    observed = {
        "flow_rate": result.flow_rate,
        "represented_vug_fraction": float(metadata["represented_vug_fraction"]),
    }
    return _compare_quantities(resolved.name, resolved.quantities, observed)

presentation_mms_references

presentation_mms_references()

Return the shipped MMS presentation-replication profiles.

Source code in src/voids/examples/mms/replication.py

def presentation_mms_references() -> tuple[MMSPresentationReference, ...]:
    """Return the shipped MMS presentation-replication profiles."""

    return _MMS_REFERENCES

presentation_vug_references

presentation_vug_references()

Return the shipped centered-vug presentation-replication profiles.

Source code in src/voids/examples/mms/replication.py

def presentation_vug_references() -> tuple[VugPresentationReference, ...]:
    """Return the shipped centered-vug presentation-replication profiles."""

    return _VUG_REFERENCES

run_presentation_mms

run_presentation_mms(
    reference,
    *,
    options=None,
    keep_solution=False,
    callback=None,
)

Run a full supplied MMS mesh sequence and compare its reported values.

Source code in src/voids/examples/mms/replication.py

def run_presentation_mms(
    reference: str | MMSPresentationReference,
    *,
    options: FEniCSSolverOptions | None = None,
    keep_solution: bool = False,
    callback: Any | None = None,
) -> MMSPresentationRun:
    """Run a full supplied MMS mesh sequence and compare its reported values."""

    resolved = _resolve_mms_reference(reference)
    result = run_mms_convergence(
        resolved.case,
        method=resolved.method,
        resolutions=resolved.resolutions,
        options=options,
        facet_law=resolved.facet_law,
        facet_size_mode=resolved.facet_size_mode,
        face_refinement=resolved.face_refinement,
        keep_solution=keep_solution,
        callback=callback,
    )
    comparison = compare_mms_with_presentation(result, resolved)
    return MMSPresentationRun(resolved, result, comparison)

run_presentation_vug

run_presentation_vug(reference, *, options=None)

Run a report-scale body-fitted vug case and compare its reported values.

Source code in src/voids/examples/mms/replication.py

def run_presentation_vug(
    reference: str | VugPresentationReference,
    *,
    options: FEniCSSolverOptions | None = None,
) -> VugPresentationRun:
    """Run a report-scale body-fitted vug case and compare its reported values."""

    resolved = _resolve_vug_reference(reference)
    result = run_centered_vug_benchmark(
        resolved.benchmark,
        method=resolved.method,
        options=options,
        facet_law=resolved.facet_law,
        facet_size_mode=(
            "facet_measure" if resolved.facet_size_mode is None else resolved.facet_size_mode
        ),
    )
    comparison = compare_vug_with_presentation(result, resolved)
    return VugPresentationRun(resolved, result, comparison)

make_body_fitted_centered_vug_mesh

make_body_fitted_centered_vug_mesh(benchmark)

Generate a tagged body-fitted centered-vug mesh with Gmsh.

Source code in src/voids/examples/mms/vug.py

def make_body_fitted_centered_vug_mesh(
    benchmark: CenteredVugBenchmark,
) -> BodyFittedVugMesh:
    """Generate a tagged body-fitted centered-vug mesh with Gmsh."""

    if benchmark.mesh_representation != "body_fitted":
        raise ValueError(
            "make_body_fitted_centered_vug_mesh requires mesh_representation='body_fitted'"
        )

    def build_model(gmsh: Any) -> None:
        if benchmark.dimension == 2:
            outer = gmsh.model.occ.addRectangle(0.0, 0.0, 0.0, 1.0, 1.0)
        else:
            outer = gmsh.model.occ.addBox(0.0, 0.0, 0.0, 1.0, 1.0, 1.0)
        if benchmark.radius > 0.0:
            if benchmark.dimension == 2:
                vug = gmsh.model.occ.addDisk(
                    0.5,
                    0.5,
                    0.0,
                    benchmark.radius,
                    benchmark.radius,
                )
            else:
                vug = gmsh.model.occ.addSphere(
                    0.5,
                    0.5,
                    0.5,
                    benchmark.radius,
                )
            gmsh.model.occ.fragment(
                [(benchmark.dimension, outer)],
                [(benchmark.dimension, vug)],
            )
        gmsh.model.occ.synchronize()
        if benchmark.radius > 0.0:
            matrix_entity, vug_entity = _classify_volume_entities(gmsh, benchmark)
        else:
            volume_entities = gmsh.model.getEntities(benchmark.dimension)
            if len(volume_entities) != 1:
                raise RuntimeError(
                    "matrix-only Gmsh model should contain exactly one volume entity"
                )
            matrix_entity = int(volume_entities[0][1])
            vug_entity = None
        add_physical_group(
            gmsh,
            benchmark.dimension,
            [matrix_entity],
            tag=_MATRIX_TAG,
            name="matrix",
        )
        if vug_entity is not None:
            add_physical_group(
                gmsh,
                benchmark.dimension,
                [vug_entity],
                tag=_VUG_TAG,
                name="vug",
            )
        outer_groups, interface = axis_aligned_boundary_entities(
            gmsh,
            benchmark.dimension,
            tolerance=1.0e-6,
        )
        missing = [name for name, entities in outer_groups.items() if not entities]
        if missing:
            raise RuntimeError(f"Gmsh did not create the expected outer facets: {missing}")
        if benchmark.radius > 0.0 and not interface:
            raise RuntimeError("Gmsh did not create a vug interface")
        if benchmark.radius == 0.0 and interface:
            raise RuntimeError("matrix-only Gmsh model unexpectedly contains interior facets")
        fdim = benchmark.dimension - 1
        for name, entities in outer_groups.items():
            marker = _FACET_MARKERS[name]
            add_physical_group(
                gmsh,
                fdim,
                entities,
                tag=marker,
                name=name,
            )
        if interface:
            add_physical_group(
                gmsh,
                fdim,
                interface,
                tag=_VUG_INTERFACE_TAG,
                name="vug_interface",
            )
        configure_uniform_mesh_size(gmsh, benchmark.target_mesh_size)
        gmsh.model.mesh.generate(benchmark.dimension)

    mesh_data = generate_dolfinx_gmsh_mesh(
        build_model,
        name=f"voids_centered_vug_{benchmark.dimension}d",
        geometric_dimension=benchmark.dimension,
    )

    if mesh_data.cell_tags is None or mesh_data.facet_tags is None:
        raise RuntimeError("Gmsh import did not preserve the required physical tags")
    return BodyFittedVugMesh(
        mesh=mesh_data.mesh,
        cell_tags=mesh_data.cell_tags,
        facet_tags=mesh_data.facet_tags,
        physical_groups=dict(mesh_data.physical_groups),
    )

run_centered_vug_benchmark

run_centered_vug_benchmark(
    benchmark,
    *,
    method="taylor_hood",
    options=None,
    tau_factor=1.0,
    tau_gamma_cap=None,
    m_t=1.0 / 3.0,
    alpha_edge=1.0,
    facet_law=None,
    facet_size_mode="facet_measure",
)

Run the pressure-driven centered-vug FEM benchmark.

Body-fitted USFEM defaults to the reaction-diffusion facet law in 2D and the shifted law in 3D, matching the documented benchmark choices. Pass facet_law explicitly whenever method-to-method comparisons need one common stabilization. The default facet size is the physical edge length in 2D and the square root of facet area in 3D, matching the body-fitted benchmark implementation. The latter is only a measure-based length, not an exact triangular-face diameter.

Source code in src/voids/examples/mms/vug.py

def run_centered_vug_benchmark(
    benchmark: CenteredVugBenchmark,
    *,
    method: MMSMethod = "taylor_hood",
    options: FEniCSSolverOptions | None = None,
    tau_factor: float = 1.0,
    tau_gamma_cap: float | None = None,
    m_t: float = 1.0 / 3.0,
    alpha_edge: float = 1.0,
    facet_law: USFEMFacetLaw | None = None,
    facet_size_mode: USFEMFacetSizeMode = "facet_measure",
) -> FEMSinglePhaseResult:
    """Run the pressure-driven centered-vug FEM benchmark.

    Body-fitted USFEM defaults to the reaction-diffusion facet law in 2D and
    the shifted law in 3D, matching the documented benchmark choices. Pass
    ``facet_law`` explicitly whenever method-to-method comparisons need one
    common stabilization. The default facet size is the physical edge length
    in 2D and the square root of facet area in 3D, matching the body-fitted
    benchmark implementation. The latter is only a measure-based length, not
    an exact triangular-face diameter.
    """

    if method not in {"taylor_hood", "usfem_p1dg0", "usfem_p1dg1"}:
        raise ValueError("method must be one of 'taylor_hood', 'usfem_p1dg0', or 'usfem_p1dg1'")
    resolved_facet_law: USFEMFacetLaw = facet_law or (
        "reaction_diffusion" if benchmark.dimension == 2 else "shifted"
    )
    _validate_usfem_controls(
        tau_factor=tau_factor,
        m_t=m_t,
        alpha_edge=alpha_edge,
        facet_law=resolved_facet_law,
        facet_size_mode=facet_size_mode,
        tau_gamma_cap=tau_gamma_cap,
    )
    uncapped_max_tau_gamma = (
        _centered_vug_p1dg0_uncapped_tau_gamma(
            benchmark,
            tau_factor=tau_factor,
            m_t=m_t,
        )
        if method == "usfem_p1dg0"
        else None
    )
    if (
        benchmark.mesh_representation == "body_fitted"
        and uncapped_max_tau_gamma is not None
        and tau_factor > 0.0
        and tau_gamma_cap is None
        and uncapped_max_tau_gamma >= 0.9
    ):
        warnings.warn(
            "The body-fitted CG1 x DG0 vug case has estimated "
            f"max(gamma * tau_K)={uncapped_max_tau_gamma:.3g}; the cell "
            "term can nearly cancel matrix drag. Set tau_gamma_cap below 1 "
            "or use tau_factor=0 as an explicit sensitivity branch.",
            RuntimeWarning,
            stacklevel=2,
        )
    if benchmark.mesh_representation == "body_fitted":
        result = _run_body_fitted(
            benchmark,
            method=method,
            options=options,
            tau_factor=tau_factor,
            tau_gamma_cap=tau_gamma_cap,
            m_t=m_t,
            alpha_edge=alpha_edge,
            facet_law=resolved_facet_law,
            facet_size_mode=facet_size_mode,
        )
        represented_fraction = float(result.metadata["represented_vug_fraction"])
    else:
        problem = benchmark.make_problem()
        if method == "taylor_hood":
            result = solve_brinkman_taylor_hood(
                problem,
                flow_axis="x",
                pressure_inlet=benchmark.pressure_inlet,
                pressure_outlet=benchmark.pressure_outlet,
                options=options,
            )
        elif method in {"usfem_p1dg0", "usfem_p1dg1"}:
            pressure_degree: Literal[0, 1] = 0 if method == "usfem_p1dg0" else 1
            result = solve_brinkman_usfem(
                problem,
                pressure_degree=pressure_degree,
                tau_factor=tau_factor,
                tau_gamma_cap=tau_gamma_cap,
                m_t=m_t,
                alpha_edge=alpha_edge,
                facet_law=resolved_facet_law,
                facet_size_mode=facet_size_mode,
                flow_axis="x",
                pressure_inlet=benchmark.pressure_inlet,
                pressure_outlet=benchmark.pressure_outlet,
                options=options,
            )
        represented_fraction = benchmark.represented_fraction
    result.metadata.update(
        {
            "benchmark": "centered_vug",
            "geometry_representation": benchmark.mesh_representation,
            "dimension": benchmark.dimension,
            "resolution": benchmark.resolution,
            "target_mesh_size": benchmark.target_mesh_size,
            "radius": benchmark.radius,
            "analytic_vug_fraction": benchmark.analytic_fraction,
            "represented_vug_fraction": represented_fraction,
            "matrix_drag": benchmark.matrix_drag,
            "vug_drag": benchmark.vug_drag,
            "matrix_effective_viscosity": benchmark.matrix_nu,
            "vug_effective_viscosity": benchmark.vug_nu,
            "p1dg0_uncapped_max_tau_gamma": uncapped_max_tau_gamma,
        }
    )
    return result

run_centered_vug_flow_case

run_centered_vug_flow_case(
    case, *, model, options=None, vms_constant=1.0
)

Solve one physical centered-vug case with Taylor--Hood P2/P1 fields.

model="darcy_brinkman" uses piecewise Brinkman viscosity mu/phi_matrix in the matrix and mu in the vug, with reaction mu/K_matrix in the matrix and zero in the vug. model="darcy_darcy" omits viscous diffusion, uses the configured finite vug permeability, and applies residual-based VMS stabilization with continuous P2 velocity and P1 pressure.

Source code in src/voids/examples/mms/vug_flow.py

def run_centered_vug_flow_case(
    case: CenteredVugFlowCase2D,
    *,
    model: CenteredVugFlowModel,
    options: FEniCSSolverOptions | None = None,
    vms_constant: float = 1.0,
) -> FEMSinglePhaseResult:
    """Solve one physical centered-vug case with Taylor--Hood P2/P1 fields.

    ``model="darcy_brinkman"`` uses piecewise Brinkman viscosity
    ``mu/phi_matrix`` in the matrix and ``mu`` in the vug, with reaction
    ``mu/K_matrix`` in the matrix and zero in the vug.
    ``model="darcy_darcy"`` omits viscous diffusion, uses the configured finite
    vug permeability, and applies residual-based VMS stabilization with
    continuous P2 velocity and P1 pressure.
    """

    if model not in {"darcy_brinkman", "darcy_darcy"}:
        raise ValueError("model must be either 'darcy_brinkman' or 'darcy_darcy'")
    if vms_constant <= 0.0 or not np.isfinite(vms_constant):
        raise ValueError("vms_constant must be positive and finite")

    solver_options = options or FEniCSSolverOptions()
    requested_dtype = _resolve_fem_linear_system_dtype(solver_options.linear_system_dtype)
    api = _require_dolfinx_core()
    selected_backend = _resolve_linear_backend(solver_options.linear_backend, api)
    if selected_backend == "petsc" and requested_dtype != "float64":
        raise ValueError("linear_system_dtype='float32' is not supported by the PETSc vug path")
    if selected_backend == "petsc":
        api = _require_dolfinx_petsc(api)

    benchmark = case.make_benchmark(model=model)
    context, represented_fraction, num_cells = _body_fitted_context(benchmark, api=api)
    mixed_space = _mixed_space(
        context.api,
        context.mesh,
        velocity_degree=2,
        pressure_family="Lagrange",
        pressure_degree=1,
    )
    u, p = context.api.ufl.TrialFunctions(mixed_space)
    v, q = context.api.ufl.TestFunctions(mixed_space)
    ufl = context.api.ufl
    gamma = context.coefficients["gamma"]
    if model == "darcy_brinkman":
        nu = context.coefficients["nu_eff"]
        form = (
            nu * ufl.inner(ufl.grad(u), ufl.grad(v)) * context.dx
            + ufl.inner(gamma * u, v) * context.dx
            - p * ufl.div(v) * context.dx
            + q * ufl.div(u) * context.dx
        )
        method = "Darcy-Brinkman Taylor-Hood CG2 x CG1"
        formulation = "centered_vug_darcy_brinkman_taylor_hood_p2p1"
    else:
        cell_diameter = ufl.CellDiameter(context.mesh)
        darcy_number = case.matrix_permeability_m2 / case.side_length_m**2
        tau_vms = ufl.min_value(
            vms_constant * cell_diameter**2 / darcy_number,
            vms_constant / gamma,
        )
        form = (
            ufl.inner(gamma * u, v) * context.dx
            - p * ufl.div(v) * context.dx
            + q * ufl.div(u) * context.dx
            + 0.5
            * ufl.inner(
                -gamma * v + ufl.grad(q),
                tau_vms * (gamma * u + ufl.grad(p)),
            )
            * context.dx
        )
        method = "Darcy-Darcy VMS Taylor-Hood CG2 x CG1"
        formulation = "centered_vug_darcy_darcy_vms_taylor_hood_p2p1"

    rhs = _pressure_boundary_load(
        context,
        v,
        flow_axis="x",
        pressure_inlet=1.0,
        pressure_outlet=0.0,
    )
    bcs = _side_wall_bcs(context, mixed_space, flow_axis="x")
    prefix_suffix = (
        f"centered_vug_flow_{model}_f{case.area_fraction:.3f}_n{case.mesh_resolution}"
    ).replace(".", "p")
    if selected_backend == "petsc":
        solution, solve_seconds, solver_metadata = _solve_mixed_problem(
            context,
            form=form,
            rhs=rhs,
            bcs=bcs,
            options=solver_options,
            prefix_suffix=prefix_suffix,
        )
    else:
        solution, solve_seconds, solver_metadata = _solve_mixed_problem_serial_direct(
            context,
            mixed_space=mixed_space,
            form=form,
            rhs=rhs,
            bcs=bcs,
            linear_backend=cast(Any, selected_backend),
            linear_system_dtype=requested_dtype,
            superlu_controls=solver_options.superlu_controls,
            umfpack_controls=solver_options.umfpack_controls,
            nvmath_cudss_controls=solver_options.nvmath_cudss_controls,
        )

    velocity_scale = (
        case.pressure_drop_pa
        * case.matrix_permeability_m2
        / (case.dynamic_viscosity_pa_s * case.side_length_m)
    )
    context.domain_length = case.side_length_m
    context.cross_section_area = 1.0
    result = _result_from_solution(
        context,
        solution,
        method=method,
        formulation=formulation,
        flow_axis="x",
        pressure_inlet=case.pressure_inlet_pa,
        pressure_outlet=case.pressure_outlet_pa,
        viscosity=case.dynamic_viscosity_pa_s,
        solve_seconds=solve_seconds,
        velocity_scale=velocity_scale,
        pressure_scale=case.pressure_drop_pa,
        metadata={
            "benchmark": "physical_centered_vug_flow_2d",
            "model": model,
            "linear_backend": selected_backend,
            "linear_system_dtype": requested_dtype,
            "velocity_degree": 2,
            "pressure_family": "Lagrange",
            "pressure_degree": 1,
            "pressure_constraint": "natural_traction",
            "returned_pressure_normalization": "zero_mean",
            "num_cells": num_cells,
            "image_shape": case.image_shape,
            "voxel_size_m": case.voxel_size_m,
            "side_length_m": case.side_length_m,
            "requested_vug_area_fraction": case.area_fraction,
            "represented_vug_area_fraction": represented_fraction,
            "vug_radius_m": case.radius_m,
            "matrix_porosity": case.matrix_porosity,
            "matrix_permeability_m2": case.matrix_permeability_m2,
            "matrix_permeability_md": case.matrix_permeability_md,
            "configured_darcy_darcy_vug_permeability_m2": (case.vug_permeability_m2),
            "vug_permeability_m2": (case.vug_permeability_m2 if model == "darcy_darcy" else None),
            "permeability_contrast": (
                case.permeability_contrast if model == "darcy_darcy" else None
            ),
            "matrix_drag_coefficient_dimensionless": benchmark.matrix_drag,
            "vug_drag_coefficient_dimensionless": benchmark.vug_drag,
            "vug_drag_pa_s_per_m2": (
                0.0
                if model == "darcy_brinkman"
                else case.dynamic_viscosity_pa_s / case.vug_permeability_m2
            ),
            "matrix_screening_length_m": case.matrix_screening_length_m,
            "base_target_mesh_size_m": (case.base_target_mesh_size_fraction * case.side_length_m),
            "mesh_size_policy": "nearly_uniform_body_fitted",
            "nondimensionalization": "matrix_darcy_velocity",
            "velocity_scale_m_per_s": velocity_scale,
            "pressure_scale_pa": case.pressure_drop_pa,
            "vms_constant": vms_constant if model == "darcy_darcy" else None,
            **solver_metadata,
        },
    )
    # The solve integrates velocity over a normalized outlet of length one.
    # Convert that mean velocity to physical 2D discharge per unit out-of-plane
    # depth while retaining the already computed physical permeability.
    result.flow_rate *= case.side_length_m
    result.cross_section_area = case.side_length_m
    return result

Two-dimensional cases

voids.examples.mms.cases_2d

boundary_layer_case_2d

boundary_layer_case_2d(*, viscosity=0.01, reaction=1.0)

Return the two-dimensional Brinkman boundary-layer MMS case.

The unit-square exact solution is

.. math::

u_1 &= y - \frac{\exp((y-1)/\nu)-\exp(-1/\nu)} {1-\exp(-1/\nu)},\\ u_2 &= x - \frac{\exp((x-1)/\nu)-\exp(-1/\nu)} {1-\exp(-1/\nu)},\\ p &= x-y.

It is exactly divergence-free and develops boundary layers at the top and right boundaries as viscosity decreases.

Source code in src/voids/examples/mms/cases_2d.py

def boundary_layer_case_2d(
    *,
    viscosity: float = 1.0e-2,
    reaction: float = 1.0,
) -> BrinkmanMMSCase:
    r"""Return the two-dimensional Brinkman boundary-layer MMS case.

    The unit-square exact solution is

    .. math::

       u_1 &= y -
       \\frac{\\exp((y-1)/\\nu)-\\exp(-1/\\nu)}
            {1-\\exp(-1/\\nu)},\\\\
       u_2 &= x -
       \\frac{\\exp((x-1)/\\nu)-\\exp(-1/\\nu)}
            {1-\\exp(-1/\\nu)},\\\\
       p &= x-y.

    It is exactly divergence-free and develops boundary layers at the top and
    right boundaries as ``viscosity`` decreases.
    """

    nu = float(viscosity)
    gamma = float(reaction)

    def exact_solution_factory(ufl, domain):
        x = ufl.SpatialCoordinate(domain)
        exponential_floor = ufl.exp(-1.0 / nu)
        denominator = 1.0 - exponential_floor
        u_1 = x[1] - (ufl.exp((x[1] - 1.0) / nu) - exponential_floor) / denominator
        u_2 = x[0] - (ufl.exp((x[0] - 1.0) / nu) - exponential_floor) / denominator
        return ufl.as_vector((u_1, u_2)), x[0] - x[1]

    def point_evaluator(points: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
        x = points[0]
        y = points[1]
        exponential_floor = np.exp(-1.0 / nu)
        denominator = -np.expm1(-1.0 / nu)
        u_1 = y - (np.exp((y - 1.0) / nu) - exponential_floor) / denominator
        u_2 = x - (np.exp((x - 1.0) / nu) - exponential_floor) / denominator
        return np.vstack((u_1, u_2)), x - y

    return BrinkmanMMSCase(
        name="brinkman_boundary_layer_2d",
        dimension=2,
        viscosity=nu,
        reaction=gamma,
        exact_solution_factory=exact_solution_factory,
        point_evaluator=point_evaluator,
        description=(
            "Divergence-free unit-square Brinkman solution with exponential layers at y=1 and x=1."
        ),
        reference=(
            "Barrenechea and Valentin (2002), Numerische Mathematik 92, "
            "653-677, doi:10.1007/s002110100371."
        ),
    )

Three-dimensional cases

voids.examples.mms.cases_3d

bubble_case_3d

bubble_case_3d(*, viscosity=0.01, reaction=1.0)

Return a smooth divergence-free three-dimensional Brinkman MMS case.

The velocity is constructed from mixed first derivatives of a polynomial boundary bubble. It therefore vanishes on the unit-cube boundary and is divergence-free by cancellation of mixed derivatives. The pressure is sin(2*pi*x) * sin(pi*y) * sin(pi*z).

Source code in src/voids/examples/mms/cases_3d.py

def bubble_case_3d(
    *,
    viscosity: float = 1.0e-2,
    reaction: float = 1.0,
) -> BrinkmanMMSCase:
    """Return a smooth divergence-free three-dimensional Brinkman MMS case.

    The velocity is constructed from mixed first derivatives of a polynomial
    boundary bubble. It therefore vanishes on the unit-cube boundary and is
    divergence-free by cancellation of mixed derivatives. The pressure is
    ``sin(2*pi*x) * sin(pi*y) * sin(pi*z)``.
    """

    nu = float(viscosity)
    gamma = float(reaction)

    def exact_solution_factory(ufl, domain):
        x = ufl.SpatialCoordinate(domain)
        bubble = (
            x[0] ** 2
            * (1.0 - x[0]) ** 2
            * x[1] ** 2
            * (1.0 - x[1]) ** 2
            * x[2] ** 2
            * (1.0 - x[2]) ** 2
        )
        phi = 32.0 * bubble * (1.0 + x[0] + 2.0 * x[1] + 3.0 * x[2])
        phi_gradient = ufl.grad(phi)
        phi_x = phi_gradient[0]
        phi_y = phi_gradient[1]
        phi_z = phi_gradient[2]
        velocity = ufl.as_vector(
            (
                3.0 * phi_y - 2.0 * phi_z,
                phi_z - 3.0 * phi_x,
                2.0 * phi_x - phi_y,
            )
        )
        pressure = ufl.sin(2.0 * np.pi * x[0]) * ufl.sin(np.pi * x[1]) * ufl.sin(np.pi * x[2])
        return velocity, pressure

    def point_evaluator(points: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
        x, y, z = points

        def bubble_1d(value: np.ndarray) -> np.ndarray:
            return value**2 * (1.0 - value) ** 2

        def bubble_1d_derivative(value: np.ndarray) -> np.ndarray:
            return 2.0 * value * (1.0 - value) * (1.0 - 2.0 * value)

        b_x, b_y, b_z = bubble_1d(x), bubble_1d(y), bubble_1d(z)
        db_x = bubble_1d_derivative(x)
        db_y = bubble_1d_derivative(y)
        db_z = bubble_1d_derivative(z)
        bubble = b_x * b_y * b_z
        linear = 1.0 + x + 2.0 * y + 3.0 * z
        phi_x = 32.0 * (db_x * b_y * b_z * linear + bubble)
        phi_y = 32.0 * (b_x * db_y * b_z * linear + 2.0 * bubble)
        phi_z = 32.0 * (b_x * b_y * db_z * linear + 3.0 * bubble)
        velocity = np.vstack(
            (
                3.0 * phi_y - 2.0 * phi_z,
                phi_z - 3.0 * phi_x,
                2.0 * phi_x - phi_y,
            )
        )
        pressure = np.sin(2.0 * np.pi * x) * np.sin(np.pi * y) * np.sin(np.pi * z)
        return velocity, pressure

    return BrinkmanMMSCase(
        name="brinkman_polynomial_bubble_3d",
        dimension=3,
        viscosity=nu,
        reaction=gamma,
        exact_solution_factory=exact_solution_factory,
        point_evaluator=point_evaluator,
        description=(
            "Smooth divergence-free unit-cube Brinkman solution generated from "
            "a polynomial boundary bubble."
        ),
        reference=(
            "Manufactured case used in the voids 3D USFEM verification study; "
            "forcing is generated from the documented strong residual."
        ),
    )

Centered-vug benchmark

voids.examples.mms.vug

CenteredVugBenchmark dataclass

Centered circular/spherical vug benchmark on the unit domain.

The defaults reproduce the documented physical configuration: radius 0.25, matrix drag 1e7, vug drag 1, viscosity 1e-2, and pressure values 1 and -1. mesh_representation="body_fitted" uses Gmsh physical cell tags. The "structured" option classifies coefficient-map cells by their centers and is useful as a portable representation-sensitivity comparison.

resolution is the number of nominal elements per coordinate direction. For body-fitted meshes, the Gmsh target size is sqrt(dimension) / resolution, matching the mesh-size convention used in the reference vug studies.

Source code in src/voids/examples/mms/vug.py

@dataclass(frozen=True, slots=True)
class CenteredVugBenchmark:
    """Centered circular/spherical vug benchmark on the unit domain.

    The defaults reproduce the documented physical configuration: radius
    ``0.25``, matrix drag ``1e7``, vug drag ``1``, viscosity ``1e-2``, and
    pressure values ``1`` and ``-1``. ``mesh_representation="body_fitted"``
    uses Gmsh physical cell tags. The ``"structured"`` option classifies
    coefficient-map cells by their centers and is useful as a portable
    representation-sensitivity comparison.

    ``resolution`` is the number of nominal elements per coordinate direction.
    For body-fitted meshes, the Gmsh target size is
    ``sqrt(dimension) / resolution``, matching the mesh-size convention used in
    the reference vug studies.
    """

    dimension: Literal[2, 3] = 2
    resolution: int = 32
    radius: float = 0.25
    viscosity: float = 1.0e-2
    matrix_drag: float = 1.0e7
    vug_drag: float = 1.0
    pressure_inlet: float = 1.0
    pressure_outlet: float = -1.0
    mesh_representation: VugMeshRepresentation = "body_fitted"
    matrix_effective_viscosity: float | None = None
    vug_effective_viscosity: float | None = None

    def __post_init__(self) -> None:
        if self.dimension not in {2, 3}:
            raise ValueError("dimension must be either 2 or 3")
        if self.resolution < 2:
            raise ValueError("resolution must be at least 2")
        if not 0.0 <= self.radius < 0.5:
            raise ValueError("radius must lie between 0 (inclusive) and 0.5")
        for name in ("viscosity", "matrix_drag"):
            value = float(getattr(self, name))
            if value <= 0.0 or not np.isfinite(value):
                raise ValueError(f"{name} must be positive and finite")
        if self.vug_drag < 0.0 or not np.isfinite(self.vug_drag):
            raise ValueError("vug_drag must be non-negative and finite")
        for name in ("matrix_effective_viscosity", "vug_effective_viscosity"):
            value = getattr(self, name)
            if value is not None and (value <= 0.0 or not np.isfinite(value)):
                raise ValueError(f"{name} must be positive and finite when provided")
        if (
            not np.isfinite(self.pressure_inlet)
            or not np.isfinite(self.pressure_outlet)
            or self.pressure_inlet <= self.pressure_outlet
        ):
            raise ValueError("pressure_inlet must be finite and greater than pressure_outlet")
        if self.mesh_representation not in {"body_fitted", "structured"}:
            raise ValueError("mesh_representation must be either 'body_fitted' or 'structured'")

    @property
    def target_mesh_size(self) -> float:
        """Return the Gmsh target element diameter."""

        return float(np.sqrt(self.dimension) / self.resolution)

    @property
    def matrix_nu(self) -> float:
        """Return the matrix coefficient multiplying the Brinkman gradient term."""

        if self.matrix_effective_viscosity is None:
            return float(self.viscosity)
        return float(self.matrix_effective_viscosity)

    @property
    def vug_nu(self) -> float:
        """Return the vug coefficient multiplying the Brinkman gradient term."""

        if self.vug_effective_viscosity is None:
            return float(self.viscosity)
        return float(self.vug_effective_viscosity)

    def vug_mask(self) -> np.ndarray:
        """Return the structured cell-center classification of the vug."""

        coordinates = (np.arange(self.resolution, dtype=float) + 0.5) / self.resolution
        grids = np.meshgrid(
            *((coordinates,) * self.dimension),
            indexing="ij",
        )
        radius_squared = sum((grid - 0.5) ** 2 for grid in grids)
        return np.asarray(radius_squared <= self.radius**2, dtype=bool)

    @property
    def represented_fraction(self) -> float:
        """Return the structured cell-volume fraction classified as vug."""

        return float(np.mean(self.vug_mask()))

    @property
    def analytic_fraction(self) -> float:
        """Return the exact circular area or spherical volume fraction."""

        if self.dimension == 2:
            return float(np.pi * self.radius**2)
        return float((4.0 / 3.0) * np.pi * self.radius**3)

    def make_problem(self) -> FEMMapProblem:
        """Build the structured constant-porosity coefficient-map problem."""

        mask = self.vug_mask()
        drag = np.where(mask, self.vug_drag, self.matrix_drag)
        permeability = np.divide(
            self.viscosity,
            drag,
            out=np.full_like(drag, np.inf, dtype=float),
            where=drag > 0.0,
        )
        if not np.all(np.isfinite(permeability)):
            raise ValueError(
                "structured vug benchmarks require positive vug_drag; "
                "zero drag is supported only by the body-fitted formulation"
            )
        effective_viscosity = np.where(mask, self.vug_nu, self.matrix_nu)
        porosity = self.viscosity / effective_viscosity
        cell_size = (1.0 / self.resolution,) * self.dimension
        return FEMMapProblem(
            permeability_map=PermeabilityMap(
                permeability,
                cell_size=cell_size,
            ),
            porosity_map=PorosityMap(
                porosity,
                cell_size=cell_size,
            ),
            viscosity=self.viscosity,
        )

target_mesh_size property

target_mesh_size

Return the Gmsh target element diameter.

matrix_nu property

matrix_nu

Return the matrix coefficient multiplying the Brinkman gradient term.

vug_nu property

vug_nu

Return the vug coefficient multiplying the Brinkman gradient term.

represented_fraction property

represented_fraction

Return the structured cell-volume fraction classified as vug.

analytic_fraction property

analytic_fraction

Return the exact circular area or spherical volume fraction.

vug_mask

vug_mask()

Return the structured cell-center classification of the vug.

Source code in src/voids/examples/mms/vug.py

def vug_mask(self) -> np.ndarray:
    """Return the structured cell-center classification of the vug."""

    coordinates = (np.arange(self.resolution, dtype=float) + 0.5) / self.resolution
    grids = np.meshgrid(
        *((coordinates,) * self.dimension),
        indexing="ij",
    )
    radius_squared = sum((grid - 0.5) ** 2 for grid in grids)
    return np.asarray(radius_squared <= self.radius**2, dtype=bool)

make_problem

make_problem()

Build the structured constant-porosity coefficient-map problem.

Source code in src/voids/examples/mms/vug.py

def make_problem(self) -> FEMMapProblem:
    """Build the structured constant-porosity coefficient-map problem."""

    mask = self.vug_mask()
    drag = np.where(mask, self.vug_drag, self.matrix_drag)
    permeability = np.divide(
        self.viscosity,
        drag,
        out=np.full_like(drag, np.inf, dtype=float),
        where=drag > 0.0,
    )
    if not np.all(np.isfinite(permeability)):
        raise ValueError(
            "structured vug benchmarks require positive vug_drag; "
            "zero drag is supported only by the body-fitted formulation"
        )
    effective_viscosity = np.where(mask, self.vug_nu, self.matrix_nu)
    porosity = self.viscosity / effective_viscosity
    cell_size = (1.0 / self.resolution,) * self.dimension
    return FEMMapProblem(
        permeability_map=PermeabilityMap(
            permeability,
            cell_size=cell_size,
        ),
        porosity_map=PorosityMap(
            porosity,
            cell_size=cell_size,
        ),
        viscosity=self.viscosity,
    )

BodyFittedVugMesh dataclass

DOLFINx mesh and physical tags generated for a centered vug.

Source code in src/voids/examples/mms/vug.py

@dataclass(slots=True)
class BodyFittedVugMesh:
    """DOLFINx mesh and physical tags generated for a centered vug."""

    mesh: Any
    cell_tags: Any
    facet_tags: Any
    physical_groups: dict[str, Any]

make_body_fitted_centered_vug_mesh

make_body_fitted_centered_vug_mesh(benchmark)

Generate a tagged body-fitted centered-vug mesh with Gmsh.

Source code in src/voids/examples/mms/vug.py

def make_body_fitted_centered_vug_mesh(
    benchmark: CenteredVugBenchmark,
) -> BodyFittedVugMesh:
    """Generate a tagged body-fitted centered-vug mesh with Gmsh."""

    if benchmark.mesh_representation != "body_fitted":
        raise ValueError(
            "make_body_fitted_centered_vug_mesh requires mesh_representation='body_fitted'"
        )

    def build_model(gmsh: Any) -> None:
        if benchmark.dimension == 2:
            outer = gmsh.model.occ.addRectangle(0.0, 0.0, 0.0, 1.0, 1.0)
        else:
            outer = gmsh.model.occ.addBox(0.0, 0.0, 0.0, 1.0, 1.0, 1.0)
        if benchmark.radius > 0.0:
            if benchmark.dimension == 2:
                vug = gmsh.model.occ.addDisk(
                    0.5,
                    0.5,
                    0.0,
                    benchmark.radius,
                    benchmark.radius,
                )
            else:
                vug = gmsh.model.occ.addSphere(
                    0.5,
                    0.5,
                    0.5,
                    benchmark.radius,
                )
            gmsh.model.occ.fragment(
                [(benchmark.dimension, outer)],
                [(benchmark.dimension, vug)],
            )
        gmsh.model.occ.synchronize()
        if benchmark.radius > 0.0:
            matrix_entity, vug_entity = _classify_volume_entities(gmsh, benchmark)
        else:
            volume_entities = gmsh.model.getEntities(benchmark.dimension)
            if len(volume_entities) != 1:
                raise RuntimeError(
                    "matrix-only Gmsh model should contain exactly one volume entity"
                )
            matrix_entity = int(volume_entities[0][1])
            vug_entity = None
        add_physical_group(
            gmsh,
            benchmark.dimension,
            [matrix_entity],
            tag=_MATRIX_TAG,
            name="matrix",
        )
        if vug_entity is not None:
            add_physical_group(
                gmsh,
                benchmark.dimension,
                [vug_entity],
                tag=_VUG_TAG,
                name="vug",
            )
        outer_groups, interface = axis_aligned_boundary_entities(
            gmsh,
            benchmark.dimension,
            tolerance=1.0e-6,
        )
        missing = [name for name, entities in outer_groups.items() if not entities]
        if missing:
            raise RuntimeError(f"Gmsh did not create the expected outer facets: {missing}")
        if benchmark.radius > 0.0 and not interface:
            raise RuntimeError("Gmsh did not create a vug interface")
        if benchmark.radius == 0.0 and interface:
            raise RuntimeError("matrix-only Gmsh model unexpectedly contains interior facets")
        fdim = benchmark.dimension - 1
        for name, entities in outer_groups.items():
            marker = _FACET_MARKERS[name]
            add_physical_group(
                gmsh,
                fdim,
                entities,
                tag=marker,
                name=name,
            )
        if interface:
            add_physical_group(
                gmsh,
                fdim,
                interface,
                tag=_VUG_INTERFACE_TAG,
                name="vug_interface",
            )
        configure_uniform_mesh_size(gmsh, benchmark.target_mesh_size)
        gmsh.model.mesh.generate(benchmark.dimension)

    mesh_data = generate_dolfinx_gmsh_mesh(
        build_model,
        name=f"voids_centered_vug_{benchmark.dimension}d",
        geometric_dimension=benchmark.dimension,
    )

    if mesh_data.cell_tags is None or mesh_data.facet_tags is None:
        raise RuntimeError("Gmsh import did not preserve the required physical tags")
    return BodyFittedVugMesh(
        mesh=mesh_data.mesh,
        cell_tags=mesh_data.cell_tags,
        facet_tags=mesh_data.facet_tags,
        physical_groups=dict(mesh_data.physical_groups),
    )

run_centered_vug_benchmark

run_centered_vug_benchmark(
    benchmark,
    *,
    method="taylor_hood",
    options=None,
    tau_factor=1.0,
    tau_gamma_cap=None,
    m_t=1.0 / 3.0,
    alpha_edge=1.0,
    facet_law=None,
    facet_size_mode="facet_measure",
)

Run the pressure-driven centered-vug FEM benchmark.

Body-fitted USFEM defaults to the reaction-diffusion facet law in 2D and the shifted law in 3D, matching the documented benchmark choices. Pass facet_law explicitly whenever method-to-method comparisons need one common stabilization. The default facet size is the physical edge length in 2D and the square root of facet area in 3D, matching the body-fitted benchmark implementation. The latter is only a measure-based length, not an exact triangular-face diameter.

Source code in src/voids/examples/mms/vug.py

def run_centered_vug_benchmark(
    benchmark: CenteredVugBenchmark,
    *,
    method: MMSMethod = "taylor_hood",
    options: FEniCSSolverOptions | None = None,
    tau_factor: float = 1.0,
    tau_gamma_cap: float | None = None,
    m_t: float = 1.0 / 3.0,
    alpha_edge: float = 1.0,
    facet_law: USFEMFacetLaw | None = None,
    facet_size_mode: USFEMFacetSizeMode = "facet_measure",
) -> FEMSinglePhaseResult:
    """Run the pressure-driven centered-vug FEM benchmark.

    Body-fitted USFEM defaults to the reaction-diffusion facet law in 2D and
    the shifted law in 3D, matching the documented benchmark choices. Pass
    ``facet_law`` explicitly whenever method-to-method comparisons need one
    common stabilization. The default facet size is the physical edge length
    in 2D and the square root of facet area in 3D, matching the body-fitted
    benchmark implementation. The latter is only a measure-based length, not
    an exact triangular-face diameter.
    """

    if method not in {"taylor_hood", "usfem_p1dg0", "usfem_p1dg1"}:
        raise ValueError("method must be one of 'taylor_hood', 'usfem_p1dg0', or 'usfem_p1dg1'")
    resolved_facet_law: USFEMFacetLaw = facet_law or (
        "reaction_diffusion" if benchmark.dimension == 2 else "shifted"
    )
    _validate_usfem_controls(
        tau_factor=tau_factor,
        m_t=m_t,
        alpha_edge=alpha_edge,
        facet_law=resolved_facet_law,
        facet_size_mode=facet_size_mode,
        tau_gamma_cap=tau_gamma_cap,
    )
    uncapped_max_tau_gamma = (
        _centered_vug_p1dg0_uncapped_tau_gamma(
            benchmark,
            tau_factor=tau_factor,
            m_t=m_t,
        )
        if method == "usfem_p1dg0"
        else None
    )
    if (
        benchmark.mesh_representation == "body_fitted"
        and uncapped_max_tau_gamma is not None
        and tau_factor > 0.0
        and tau_gamma_cap is None
        and uncapped_max_tau_gamma >= 0.9
    ):
        warnings.warn(
            "The body-fitted CG1 x DG0 vug case has estimated "
            f"max(gamma * tau_K)={uncapped_max_tau_gamma:.3g}; the cell "
            "term can nearly cancel matrix drag. Set tau_gamma_cap below 1 "
            "or use tau_factor=0 as an explicit sensitivity branch.",
            RuntimeWarning,
            stacklevel=2,
        )
    if benchmark.mesh_representation == "body_fitted":
        result = _run_body_fitted(
            benchmark,
            method=method,
            options=options,
            tau_factor=tau_factor,
            tau_gamma_cap=tau_gamma_cap,
            m_t=m_t,
            alpha_edge=alpha_edge,
            facet_law=resolved_facet_law,
            facet_size_mode=facet_size_mode,
        )
        represented_fraction = float(result.metadata["represented_vug_fraction"])
    else:
        problem = benchmark.make_problem()
        if method == "taylor_hood":
            result = solve_brinkman_taylor_hood(
                problem,
                flow_axis="x",
                pressure_inlet=benchmark.pressure_inlet,
                pressure_outlet=benchmark.pressure_outlet,
                options=options,
            )
        elif method in {"usfem_p1dg0", "usfem_p1dg1"}:
            pressure_degree: Literal[0, 1] = 0 if method == "usfem_p1dg0" else 1
            result = solve_brinkman_usfem(
                problem,
                pressure_degree=pressure_degree,
                tau_factor=tau_factor,
                tau_gamma_cap=tau_gamma_cap,
                m_t=m_t,
                alpha_edge=alpha_edge,
                facet_law=resolved_facet_law,
                facet_size_mode=facet_size_mode,
                flow_axis="x",
                pressure_inlet=benchmark.pressure_inlet,
                pressure_outlet=benchmark.pressure_outlet,
                options=options,
            )
        represented_fraction = benchmark.represented_fraction
    result.metadata.update(
        {
            "benchmark": "centered_vug",
            "geometry_representation": benchmark.mesh_representation,
            "dimension": benchmark.dimension,
            "resolution": benchmark.resolution,
            "target_mesh_size": benchmark.target_mesh_size,
            "radius": benchmark.radius,
            "analytic_vug_fraction": benchmark.analytic_fraction,
            "represented_vug_fraction": represented_fraction,
            "matrix_drag": benchmark.matrix_drag,
            "vug_drag": benchmark.vug_drag,
            "matrix_effective_viscosity": benchmark.matrix_nu,
            "vug_effective_viscosity": benchmark.vug_nu,
            "p1dg0_uncapped_max_tau_gamma": uncapped_max_tau_gamma,
        }
    )
    return result

Physical centered-vug flow family

voids.examples.mms.vug_flow

CenteredVugFlowCase2D dataclass

Physical 2D centered-vug case solved on a body-fitted unit-square mesh.

The continuum domain size is derived from image_shape and voxel_size_m. Gmsh coordinates are normalized to the unit square for numerical conditioning; the solver nondimensionalizes the equations and converts velocity, pressure, flux, and permeability back to SI units.

The Darcy--Brinkman branch follows the layered-domain model: its reaction coefficient is mu / K_matrix in the matrix and exactly zero in the vug. vug_permeability_m2 is used only by the Darcy--Darcy branch as a finite high-permeability closure. It is not an intrinsic permeability measurement of an open cavity.

Source code in src/voids/examples/mms/vug_flow.py

@dataclass(frozen=True, slots=True)
class CenteredVugFlowCase2D:
    """Physical 2D centered-vug case solved on a body-fitted unit-square mesh.

    The continuum domain size is derived from ``image_shape`` and
    ``voxel_size_m``. Gmsh coordinates are normalized to the unit square for
    numerical conditioning; the solver nondimensionalizes the equations and
    converts velocity, pressure, flux, and permeability back to SI units.

    The Darcy--Brinkman branch follows the layered-domain model: its reaction
    coefficient is ``mu / K_matrix`` in the matrix and exactly zero in the
    vug. ``vug_permeability_m2`` is used only by the Darcy--Darcy branch as a
    finite high-permeability closure. It is not an intrinsic permeability
    measurement of an open cavity.
    """

    area_fraction: float
    image_shape: tuple[int, int] = (500, 500)
    voxel_size_m: float = 15.0e-6
    matrix_porosity: float = 0.2
    matrix_permeability_md: float = 200.0
    vug_permeability_m2: float = 1.0e-8
    dynamic_viscosity_pa_s: float = 1.0e-3
    pressure_inlet_pa: float = 1.0
    pressure_outlet_pa: float = 0.0
    mesh_resolution: int = 100

    def __post_init__(self) -> None:
        if len(self.image_shape) != 2 or self.image_shape[0] != self.image_shape[1]:
            raise ValueError("image_shape must describe a square 2D image")
        if self.image_shape[0] < 2:
            raise ValueError("image_shape values must be at least 2")
        if self.voxel_size_m <= 0.0 or not np.isfinite(self.voxel_size_m):
            raise ValueError("voxel_size_m must be positive and finite")
        maximum_fraction = float(np.pi / 4.0)
        if (
            self.area_fraction < 0.0
            or self.area_fraction >= maximum_fraction
            or not np.isfinite(self.area_fraction)
        ):
            raise ValueError("area_fraction must lie between 0 (inclusive) and pi/4 (exclusive)")
        for name in (
            "matrix_porosity",
            "matrix_permeability_md",
            "vug_permeability_m2",
            "dynamic_viscosity_pa_s",
        ):
            value = float(getattr(self, name))
            if value <= 0.0 or not np.isfinite(value):
                raise ValueError(f"{name} must be positive and finite")
        if self.matrix_porosity > 1.0:
            raise ValueError("matrix_porosity must not exceed 1")
        if (
            not np.isfinite(self.pressure_inlet_pa)
            or not np.isfinite(self.pressure_outlet_pa)
            or self.pressure_inlet_pa <= self.pressure_outlet_pa
        ):
            raise ValueError("pressure_inlet_pa must be finite and greater than pressure_outlet_pa")
        if self.mesh_resolution < 8:
            raise ValueError("mesh_resolution must be at least 8")

    @property
    def side_length_m(self) -> float:
        """Return the square side length derived from the image metadata."""

        return float(self.image_shape[0] * self.voxel_size_m)

    @property
    def radius_fraction(self) -> float:
        """Return the circle radius divided by the square side length."""

        return float(np.sqrt(self.area_fraction / np.pi))

    @property
    def radius_m(self) -> float:
        """Return the physical centered-vug radius."""

        return float(self.radius_fraction * self.side_length_m)

    @property
    def matrix_permeability_m2(self) -> float:
        """Return the matrix permeability converted from mD to square metres."""

        return float(self.matrix_permeability_md * M2_PER_MILLIDARCY)

    @property
    def matrix_screening_length_m(self) -> float:
        """Return the matrix Brinkman screening length ``sqrt(K / phi)``."""

        return float(np.sqrt(self.matrix_permeability_m2 / self.matrix_porosity))

    @property
    def permeability_contrast(self) -> float:
        """Return the Darcy--Darcy closure contrast ``K_vug / K_matrix``."""

        return float(self.vug_permeability_m2 / self.matrix_permeability_m2)

    @property
    def pressure_drop_pa(self) -> float:
        """Return the applied pressure drop."""

        return float(self.pressure_inlet_pa - self.pressure_outlet_pa)

    @property
    def base_target_mesh_size_fraction(self) -> float:
        """Return the far-field target diameter divided by the side length."""

        return float(np.sqrt(2.0) / self.mesh_resolution)

    def make_benchmark(
        self,
        *,
        model: CenteredVugFlowModel = "darcy_brinkman",
    ) -> CenteredVugBenchmark:
        """Return the model-specific nondimensional mesh/coefficient definition."""

        if model not in {"darcy_brinkman", "darcy_darcy"}:
            raise ValueError("model must be either 'darcy_brinkman' or 'darcy_darcy'")
        darcy_number = self.matrix_permeability_m2 / self.side_length_m**2
        vug_drag = (
            0.0
            if model == "darcy_brinkman"
            else self.matrix_permeability_m2 / self.vug_permeability_m2
        )
        return CenteredVugBenchmark(
            dimension=2,
            resolution=self.mesh_resolution,
            radius=self.radius_fraction,
            viscosity=darcy_number,
            matrix_drag=1.0,
            vug_drag=vug_drag,
            pressure_inlet=1.0,
            pressure_outlet=0.0,
            mesh_representation="body_fitted",
            matrix_effective_viscosity=darcy_number / self.matrix_porosity,
            vug_effective_viscosity=darcy_number,
        )

side_length_m property

side_length_m

Return the square side length derived from the image metadata.

radius_fraction property

radius_fraction

Return the circle radius divided by the square side length.

radius_m property

radius_m

Return the physical centered-vug radius.

matrix_permeability_m2 property

matrix_permeability_m2

Return the matrix permeability converted from mD to square metres.

matrix_screening_length_m property

matrix_screening_length_m

Return the matrix Brinkman screening length sqrt(K / phi).

permeability_contrast property

permeability_contrast

Return the Darcy--Darcy closure contrast K_vug / K_matrix.

pressure_drop_pa property

pressure_drop_pa

Return the applied pressure drop.

base_target_mesh_size_fraction property

base_target_mesh_size_fraction

Return the far-field target diameter divided by the side length.

make_benchmark

make_benchmark(*, model='darcy_brinkman')

Return the model-specific nondimensional mesh/coefficient definition.

Source code in src/voids/examples/mms/vug_flow.py

def make_benchmark(
    self,
    *,
    model: CenteredVugFlowModel = "darcy_brinkman",
) -> CenteredVugBenchmark:
    """Return the model-specific nondimensional mesh/coefficient definition."""

    if model not in {"darcy_brinkman", "darcy_darcy"}:
        raise ValueError("model must be either 'darcy_brinkman' or 'darcy_darcy'")
    darcy_number = self.matrix_permeability_m2 / self.side_length_m**2
    vug_drag = (
        0.0
        if model == "darcy_brinkman"
        else self.matrix_permeability_m2 / self.vug_permeability_m2
    )
    return CenteredVugBenchmark(
        dimension=2,
        resolution=self.mesh_resolution,
        radius=self.radius_fraction,
        viscosity=darcy_number,
        matrix_drag=1.0,
        vug_drag=vug_drag,
        pressure_inlet=1.0,
        pressure_outlet=0.0,
        mesh_representation="body_fitted",
        matrix_effective_viscosity=darcy_number / self.matrix_porosity,
        vug_effective_viscosity=darcy_number,
    )

run_centered_vug_flow_case

run_centered_vug_flow_case(
    case, *, model, options=None, vms_constant=1.0
)

Solve one physical centered-vug case with Taylor--Hood P2/P1 fields.

model="darcy_brinkman" uses piecewise Brinkman viscosity mu/phi_matrix in the matrix and mu in the vug, with reaction mu/K_matrix in the matrix and zero in the vug. model="darcy_darcy" omits viscous diffusion, uses the configured finite vug permeability, and applies residual-based VMS stabilization with continuous P2 velocity and P1 pressure.

Source code in src/voids/examples/mms/vug_flow.py

def run_centered_vug_flow_case(
    case: CenteredVugFlowCase2D,
    *,
    model: CenteredVugFlowModel,
    options: FEniCSSolverOptions | None = None,
    vms_constant: float = 1.0,
) -> FEMSinglePhaseResult:
    """Solve one physical centered-vug case with Taylor--Hood P2/P1 fields.

    ``model="darcy_brinkman"`` uses piecewise Brinkman viscosity
    ``mu/phi_matrix`` in the matrix and ``mu`` in the vug, with reaction
    ``mu/K_matrix`` in the matrix and zero in the vug.
    ``model="darcy_darcy"`` omits viscous diffusion, uses the configured finite
    vug permeability, and applies residual-based VMS stabilization with
    continuous P2 velocity and P1 pressure.
    """

    if model not in {"darcy_brinkman", "darcy_darcy"}:
        raise ValueError("model must be either 'darcy_brinkman' or 'darcy_darcy'")
    if vms_constant <= 0.0 or not np.isfinite(vms_constant):
        raise ValueError("vms_constant must be positive and finite")

    solver_options = options or FEniCSSolverOptions()
    requested_dtype = _resolve_fem_linear_system_dtype(solver_options.linear_system_dtype)
    api = _require_dolfinx_core()
    selected_backend = _resolve_linear_backend(solver_options.linear_backend, api)
    if selected_backend == "petsc" and requested_dtype != "float64":
        raise ValueError("linear_system_dtype='float32' is not supported by the PETSc vug path")
    if selected_backend == "petsc":
        api = _require_dolfinx_petsc(api)

    benchmark = case.make_benchmark(model=model)
    context, represented_fraction, num_cells = _body_fitted_context(benchmark, api=api)
    mixed_space = _mixed_space(
        context.api,
        context.mesh,
        velocity_degree=2,
        pressure_family="Lagrange",
        pressure_degree=1,
    )
    u, p = context.api.ufl.TrialFunctions(mixed_space)
    v, q = context.api.ufl.TestFunctions(mixed_space)
    ufl = context.api.ufl
    gamma = context.coefficients["gamma"]
    if model == "darcy_brinkman":
        nu = context.coefficients["nu_eff"]
        form = (
            nu * ufl.inner(ufl.grad(u), ufl.grad(v)) * context.dx
            + ufl.inner(gamma * u, v) * context.dx
            - p * ufl.div(v) * context.dx
            + q * ufl.div(u) * context.dx
        )
        method = "Darcy-Brinkman Taylor-Hood CG2 x CG1"
        formulation = "centered_vug_darcy_brinkman_taylor_hood_p2p1"
    else:
        cell_diameter = ufl.CellDiameter(context.mesh)
        darcy_number = case.matrix_permeability_m2 / case.side_length_m**2
        tau_vms = ufl.min_value(
            vms_constant * cell_diameter**2 / darcy_number,
            vms_constant / gamma,
        )
        form = (
            ufl.inner(gamma * u, v) * context.dx
            - p * ufl.div(v) * context.dx
            + q * ufl.div(u) * context.dx
            + 0.5
            * ufl.inner(
                -gamma * v + ufl.grad(q),
                tau_vms * (gamma * u + ufl.grad(p)),
            )
            * context.dx
        )
        method = "Darcy-Darcy VMS Taylor-Hood CG2 x CG1"
        formulation = "centered_vug_darcy_darcy_vms_taylor_hood_p2p1"

    rhs = _pressure_boundary_load(
        context,
        v,
        flow_axis="x",
        pressure_inlet=1.0,
        pressure_outlet=0.0,
    )
    bcs = _side_wall_bcs(context, mixed_space, flow_axis="x")
    prefix_suffix = (
        f"centered_vug_flow_{model}_f{case.area_fraction:.3f}_n{case.mesh_resolution}"
    ).replace(".", "p")
    if selected_backend == "petsc":
        solution, solve_seconds, solver_metadata = _solve_mixed_problem(
            context,
            form=form,
            rhs=rhs,
            bcs=bcs,
            options=solver_options,
            prefix_suffix=prefix_suffix,
        )
    else:
        solution, solve_seconds, solver_metadata = _solve_mixed_problem_serial_direct(
            context,
            mixed_space=mixed_space,
            form=form,
            rhs=rhs,
            bcs=bcs,
            linear_backend=cast(Any, selected_backend),
            linear_system_dtype=requested_dtype,
            superlu_controls=solver_options.superlu_controls,
            umfpack_controls=solver_options.umfpack_controls,
            nvmath_cudss_controls=solver_options.nvmath_cudss_controls,
        )

    velocity_scale = (
        case.pressure_drop_pa
        * case.matrix_permeability_m2
        / (case.dynamic_viscosity_pa_s * case.side_length_m)
    )
    context.domain_length = case.side_length_m
    context.cross_section_area = 1.0
    result = _result_from_solution(
        context,
        solution,
        method=method,
        formulation=formulation,
        flow_axis="x",
        pressure_inlet=case.pressure_inlet_pa,
        pressure_outlet=case.pressure_outlet_pa,
        viscosity=case.dynamic_viscosity_pa_s,
        solve_seconds=solve_seconds,
        velocity_scale=velocity_scale,
        pressure_scale=case.pressure_drop_pa,
        metadata={
            "benchmark": "physical_centered_vug_flow_2d",
            "model": model,
            "linear_backend": selected_backend,
            "linear_system_dtype": requested_dtype,
            "velocity_degree": 2,
            "pressure_family": "Lagrange",
            "pressure_degree": 1,
            "pressure_constraint": "natural_traction",
            "returned_pressure_normalization": "zero_mean",
            "num_cells": num_cells,
            "image_shape": case.image_shape,
            "voxel_size_m": case.voxel_size_m,
            "side_length_m": case.side_length_m,
            "requested_vug_area_fraction": case.area_fraction,
            "represented_vug_area_fraction": represented_fraction,
            "vug_radius_m": case.radius_m,
            "matrix_porosity": case.matrix_porosity,
            "matrix_permeability_m2": case.matrix_permeability_m2,
            "matrix_permeability_md": case.matrix_permeability_md,
            "configured_darcy_darcy_vug_permeability_m2": (case.vug_permeability_m2),
            "vug_permeability_m2": (case.vug_permeability_m2 if model == "darcy_darcy" else None),
            "permeability_contrast": (
                case.permeability_contrast if model == "darcy_darcy" else None
            ),
            "matrix_drag_coefficient_dimensionless": benchmark.matrix_drag,
            "vug_drag_coefficient_dimensionless": benchmark.vug_drag,
            "vug_drag_pa_s_per_m2": (
                0.0
                if model == "darcy_brinkman"
                else case.dynamic_viscosity_pa_s / case.vug_permeability_m2
            ),
            "matrix_screening_length_m": case.matrix_screening_length_m,
            "base_target_mesh_size_m": (case.base_target_mesh_size_fraction * case.side_length_m),
            "mesh_size_policy": "nearly_uniform_body_fitted",
            "nondimensionalization": "matrix_darcy_velocity",
            "velocity_scale_m_per_s": velocity_scale,
            "pressure_scale_pa": case.pressure_drop_pa,
            "vms_constant": vms_constant if model == "darcy_darcy" else None,
            **solver_metadata,
        },
    )
    # The solve integrates velocity over a normalized outlet of length one.
    # Convert that mean velocity to physical 2D discharge per unit out-of-plane
    # depth while retaining the already computed physical permeability.
    result.flow_rate *= case.side_length_m
    result.cross_section_area = case.side_length_m
    return result

Presentation-replication profiles

voids.examples.mms.replication

ReferenceQuantity dataclass

One reported scalar value and its replication tolerance.

Source code in src/voids/examples/mms/replication.py

@dataclass(frozen=True, slots=True)
class ReferenceQuantity:
    """One reported scalar value and its replication tolerance."""

    metric: str
    expected: float
    absolute_tolerance: float = 0.0
    relative_tolerance: float = 0.0

    def __post_init__(self) -> None:
        if not self.metric.strip():
            raise ValueError("metric must not be empty")
        if not np.isfinite(self.expected):
            raise ValueError("expected must be finite")
        for name in ("absolute_tolerance", "relative_tolerance"):
            value = float(getattr(self, name))
            if value < 0.0 or not np.isfinite(value):
                raise ValueError(f"{name} must be nonnegative and finite")
        if self.absolute_tolerance == 0.0 and self.relative_tolerance == 0.0:
            raise ValueError("at least one comparison tolerance must be positive")

ReferenceComparison dataclass

Observed-versus-reported comparison for one scalar quantity.

Source code in src/voids/examples/mms/replication.py

@dataclass(frozen=True, slots=True)
class ReferenceComparison:
    """Observed-versus-reported comparison for one scalar quantity."""

    metric: str
    observed: float
    expected: float
    absolute_tolerance: float
    relative_tolerance: float

    @property
    def absolute_error(self) -> float:
        """Return ``abs(observed - expected)``."""

        return abs(self.observed - self.expected)

    @property
    def relative_error(self) -> float:
        """Return the error relative to the reported magnitude."""

        if self.expected == 0.0:
            return 0.0 if self.observed == 0.0 else float("inf")
        return self.absolute_error / abs(self.expected)

    @property
    def passed(self) -> bool:
        """Return whether the observation lies within the stored tolerance."""

        return bool(
            np.isclose(
                self.observed,
                self.expected,
                atol=self.absolute_tolerance,
                rtol=self.relative_tolerance,
            )
        )

    def as_dict(self) -> dict[str, str | float | bool]:
        """Return a table-ready representation."""

        return {
            "metric": self.metric,
            "observed": self.observed,
            "expected": self.expected,
            "absolute_error": self.absolute_error,
            "relative_error": self.relative_error,
            "passed": self.passed,
        }

absolute_error property

absolute_error

Return abs(observed - expected).

relative_error property

relative_error

Return the error relative to the reported magnitude.

passed property

passed

Return whether the observation lies within the stored tolerance.

as_dict

as_dict()

Return a table-ready representation.

Source code in src/voids/examples/mms/replication.py

def as_dict(self) -> dict[str, str | float | bool]:
    """Return a table-ready representation."""

    return {
        "metric": self.metric,
        "observed": self.observed,
        "expected": self.expected,
        "absolute_error": self.absolute_error,
        "relative_error": self.relative_error,
        "passed": self.passed,
    }

PresentationComparison dataclass

Comparison of one live solve with a supplied presentation baseline.

Source code in src/voids/examples/mms/replication.py

@dataclass(frozen=True, slots=True)
class PresentationComparison:
    """Comparison of one live solve with a supplied presentation baseline."""

    reference_name: str
    quantities: tuple[ReferenceComparison, ...]

    @property
    def passed(self) -> bool:
        """Return whether every reported quantity was reproduced."""

        return all(quantity.passed for quantity in self.quantities)

    def as_dicts(self) -> list[dict[str, str | float | bool]]:
        """Return comparison rows without requiring pandas."""

        return [quantity.as_dict() for quantity in self.quantities]

    def assert_matches(self) -> None:
        """Raise when any reported scalar falls outside its tolerance."""

        failures = [
            (
                f"{quantity.metric}: observed {quantity.observed:.6g}, "
                f"expected {quantity.expected:.6g}, "
                f"absolute error {quantity.absolute_error:.3g}"
            )
            for quantity in self.quantities
            if not quantity.passed
        ]
        if failures:
            raise AssertionError(
                f"{self.reference_name} did not reproduce the supplied baseline: "
                + "; ".join(failures)
            )

passed property

passed

Return whether every reported quantity was reproduced.

as_dicts

as_dicts()

Return comparison rows without requiring pandas.

Source code in src/voids/examples/mms/replication.py

def as_dicts(self) -> list[dict[str, str | float | bool]]:
    """Return comparison rows without requiring pandas."""

    return [quantity.as_dict() for quantity in self.quantities]

assert_matches

assert_matches()

Raise when any reported scalar falls outside its tolerance.

Source code in src/voids/examples/mms/replication.py

def assert_matches(self) -> None:
    """Raise when any reported scalar falls outside its tolerance."""

    failures = [
        (
            f"{quantity.metric}: observed {quantity.observed:.6g}, "
            f"expected {quantity.expected:.6g}, "
            f"absolute error {quantity.absolute_error:.3g}"
        )
        for quantity in self.quantities
        if not quantity.passed
    ]
    if failures:
        raise AssertionError(
            f"{self.reference_name} did not reproduce the supplied baseline: "
            + "; ".join(failures)
        )

MMSPresentationReference dataclass

Exact configuration and reported values for one MMS presentation row.

Source code in src/voids/examples/mms/replication.py

@dataclass(frozen=True, slots=True)
class MMSPresentationReference:
    """Exact configuration and reported values for one MMS presentation row."""

    name: str
    case: BrinkmanMMSCase
    method: MMSMethod
    resolutions: tuple[int, ...]
    facet_law: MMSFacetLaw
    facet_size_mode: MMSFacetSizeMode
    quantities: tuple[ReferenceQuantity, ...]
    source: str
    face_refinement: int = 24

    def __post_init__(self) -> None:
        if not self.name.strip():
            raise ValueError("name must not be empty")
        if len(self.resolutions) < 2:
            raise ValueError("resolutions must contain at least two mesh levels")
        if any(value < 1 for value in self.resolutions):
            raise ValueError("all resolutions must be positive")
        if any(
            current <= previous for previous, current in zip(self.resolutions, self.resolutions[1:])
        ):
            raise ValueError("resolutions must be strictly increasing")
        if not self.quantities:
            raise ValueError("quantities must not be empty")
        if not self.source.strip():
            raise ValueError("source must not be empty")
        if self.face_refinement < 2:
            raise ValueError("face_refinement must be at least 2")
        if self.facet_size_mode not in {
            "cell_diameter",
            "facet_diameter",
            "representative",
        }:
            raise ValueError("facet_size_mode is not supported")

VugPresentationReference dataclass

Configuration and reported flux for a centered-vug presentation row.

Source code in src/voids/examples/mms/replication.py

@dataclass(frozen=True, slots=True)
class VugPresentationReference:
    """Configuration and reported flux for a centered-vug presentation row."""

    name: str
    benchmark: CenteredVugBenchmark
    method: MMSMethod
    facet_law: USFEMFacetLaw | None
    facet_size_mode: USFEMFacetSizeMode | None
    quantities: tuple[ReferenceQuantity, ...]
    source: str

    def __post_init__(self) -> None:
        if not self.name.strip():
            raise ValueError("name must not be empty")
        if self.benchmark.mesh_representation != "body_fitted":
            raise ValueError("presentation vug references require a body-fitted mesh")
        if not self.quantities:
            raise ValueError("quantities must not be empty")
        if self.method != "taylor_hood" and self.facet_size_mode is None:
            raise ValueError("USFEM vug references require a facet_size_mode")
        if not self.source.strip():
            raise ValueError("source must not be empty")

MMSPresentationRun dataclass

Live MMS refinement result paired with its baseline comparison.

Source code in src/voids/examples/mms/replication.py

@dataclass(slots=True)
class MMSPresentationRun:
    """Live MMS refinement result paired with its baseline comparison."""

    reference: MMSPresentationReference
    result: MMSConvergenceResult
    comparison: PresentationComparison

    def assert_matches(self) -> None:
        """Raise unless the live result reproduces every stored target."""

        self.comparison.assert_matches()

assert_matches

assert_matches()

Raise unless the live result reproduces every stored target.

Source code in src/voids/examples/mms/replication.py

def assert_matches(self) -> None:
    """Raise unless the live result reproduces every stored target."""

    self.comparison.assert_matches()

VugPresentationRun dataclass

Live centered-vug result paired with its baseline comparison.

Source code in src/voids/examples/mms/replication.py

@dataclass(slots=True)
class VugPresentationRun:
    """Live centered-vug result paired with its baseline comparison."""

    reference: VugPresentationReference
    result: FEMSinglePhaseResult
    comparison: PresentationComparison

    def assert_matches(self) -> None:
        """Raise unless the live result reproduces every stored target."""

        self.comparison.assert_matches()

assert_matches

assert_matches()

Raise unless the live result reproduces every stored target.

Source code in src/voids/examples/mms/replication.py

def assert_matches(self) -> None:
    """Raise unless the live result reproduces every stored target."""

    self.comparison.assert_matches()

presentation_mms_references

presentation_mms_references()

Return the shipped MMS presentation-replication profiles.

Source code in src/voids/examples/mms/replication.py

def presentation_mms_references() -> tuple[MMSPresentationReference, ...]:
    """Return the shipped MMS presentation-replication profiles."""

    return _MMS_REFERENCES

presentation_vug_references

presentation_vug_references()

Return the shipped centered-vug presentation-replication profiles.

Source code in src/voids/examples/mms/replication.py

def presentation_vug_references() -> tuple[VugPresentationReference, ...]:
    """Return the shipped centered-vug presentation-replication profiles."""

    return _VUG_REFERENCES

compare_mms_with_presentation

compare_mms_with_presentation(result, reference)

Compare a live MMS result with the exact supplied configuration and values.

Source code in src/voids/examples/mms/replication.py

def compare_mms_with_presentation(
    result: MMSConvergenceResult,
    reference: str | MMSPresentationReference,
) -> PresentationComparison:
    """Compare a live MMS result with the exact supplied configuration and values."""

    resolved = _resolve_mms_reference(reference)
    if result.method != resolved.method:
        raise ValueError(
            f"method mismatch: result uses {result.method}, reference uses {resolved.method}"
        )
    if result.case.name != resolved.case.name:
        raise ValueError(
            f"case mismatch: result uses {result.case.name}, reference uses {resolved.case.name}"
        )
    if not np.isclose(result.case.viscosity, resolved.case.viscosity) or not np.isclose(
        result.case.reaction,
        resolved.case.reaction,
    ):
        raise ValueError("case coefficient mismatch with presentation reference")
    result_resolutions = tuple(level.resolution for level in result.levels)
    if result_resolutions[-2:] != resolved.resolutions[-2:]:
        raise ValueError(
            "the result must end with the presentation's finest mesh pair "
            f"{resolved.resolutions[-2:]}; received {result_resolutions[-2:]}"
        )
    if result.method != "taylor_hood" and result.metadata.get("facet_law") != resolved.facet_law:
        raise ValueError(
            "facet-law mismatch: result uses "
            f"{result.metadata.get('facet_law')}, reference uses {resolved.facet_law}"
        )
    if (
        result.method != "taylor_hood"
        and result.metadata.get("facet_size_mode") != resolved.facet_size_mode
    ):
        raise ValueError(
            "facet-size mismatch: result uses "
            f"{result.metadata.get('facet_size_mode')}, "
            f"reference uses {resolved.facet_size_mode}"
        )

    finest = result.levels[-1]
    observed = {
        "velocity_l2_error": finest.velocity_l2_error,
        "velocity_h1_error": finest.velocity_h1_error,
        "pressure_l2_error": finest.pressure_l2_error,
        "divergence_l2": finest.divergence_l2,
        **{f"{name}_rate": value for name, value in result.last_rates.items()},
    }
    return _compare_quantities(resolved.name, resolved.quantities, observed)

compare_vug_with_presentation

compare_vug_with_presentation(result, reference)

Compare a live centered-vug result with the supplied report-scale values.

Source code in src/voids/examples/mms/replication.py

def compare_vug_with_presentation(
    result: FEMSinglePhaseResult,
    reference: str | VugPresentationReference,
) -> PresentationComparison:
    """Compare a live centered-vug result with the supplied report-scale values."""

    resolved = _resolve_vug_reference(reference)
    metadata = result.metadata
    for name, expected in (
        ("dimension", resolved.benchmark.dimension),
        ("resolution", resolved.benchmark.resolution),
        ("radius", resolved.benchmark.radius),
        ("matrix_drag", resolved.benchmark.matrix_drag),
        ("vug_drag", resolved.benchmark.vug_drag),
    ):
        actual = metadata.get(name)
        if actual is None or not np.isclose(
            float(actual),
            float(expected),
            rtol=1.0e-12,
            atol=0.0,
        ):
            raise ValueError(f"vug configuration mismatch for {name}: {actual!r} != {expected!r}")
    if resolved.method == "taylor_hood":
        if result.formulation != "brinkman_taylor_hood_p2p1":
            raise ValueError("vug formulation does not match the Taylor-Hood reference")
    else:
        if result.formulation != "brinkman_usfem_p1dg1":
            raise ValueError("vug formulation does not match the P1/DG1 reference")
        if metadata.get("facet_law") != resolved.facet_law:
            raise ValueError("vug facet law does not match the presentation reference")
        if metadata.get("facet_size_mode") != resolved.facet_size_mode:
            raise ValueError("vug facet size does not match the presentation reference")

    observed = {
        "flow_rate": result.flow_rate,
        "represented_vug_fraction": float(metadata["represented_vug_fraction"]),
    }
    return _compare_quantities(resolved.name, resolved.quantities, observed)

run_presentation_mms

run_presentation_mms(
    reference,
    *,
    options=None,
    keep_solution=False,
    callback=None,
)

Run a full supplied MMS mesh sequence and compare its reported values.

Source code in src/voids/examples/mms/replication.py

def run_presentation_mms(
    reference: str | MMSPresentationReference,
    *,
    options: FEniCSSolverOptions | None = None,
    keep_solution: bool = False,
    callback: Any | None = None,
) -> MMSPresentationRun:
    """Run a full supplied MMS mesh sequence and compare its reported values."""

    resolved = _resolve_mms_reference(reference)
    result = run_mms_convergence(
        resolved.case,
        method=resolved.method,
        resolutions=resolved.resolutions,
        options=options,
        facet_law=resolved.facet_law,
        facet_size_mode=resolved.facet_size_mode,
        face_refinement=resolved.face_refinement,
        keep_solution=keep_solution,
        callback=callback,
    )
    comparison = compare_mms_with_presentation(result, resolved)
    return MMSPresentationRun(resolved, result, comparison)

run_presentation_vug

run_presentation_vug(reference, *, options=None)

Run a report-scale body-fitted vug case and compare its reported values.

Source code in src/voids/examples/mms/replication.py

def run_presentation_vug(
    reference: str | VugPresentationReference,
    *,
    options: FEniCSSolverOptions | None = None,
) -> VugPresentationRun:
    """Run a report-scale body-fitted vug case and compare its reported values."""

    resolved = _resolve_vug_reference(reference)
    result = run_centered_vug_benchmark(
        resolved.benchmark,
        method=resolved.method,
        options=options,
        facet_law=resolved.facet_law,
        facet_size_mode=(
            "facet_measure" if resolved.facet_size_mode is None else resolved.facet_size_mode
        ),
    )
    comparison = compare_vug_with_presentation(result, resolved)
    return VugPresentationRun(resolved, result, comparison)