Lattice Boltzmann

The voids.lbm sub-package provides package-facing LBM namespaces for direct-image single-phase flow. The current Stokes-limit implementation wraps the optional XLB adapter owned by voids.lbm.singlephase.xlb.

For the binary image convention, lattice pressure relation, convergence metric, permeability conversion, and guidance on tuning LBM runs for a new sample, see Map-Based Single-Phase Solvers.

---

XLB Backend

voids.lbm.singlephase.xlb

Direct-image single-phase XLB adapters for binary segmented volumes.

XLBConvergenceWarning

Bases: RuntimeWarning

Warn that an XLB solve reached max_steps before steady convergence.

Source code in src/voids/lbm/singlephase/xlb.py

class XLBConvergenceWarning(RuntimeWarning):
    """Warn that an XLB solve reached ``max_steps`` before steady convergence."""

XLBOptions dataclass

Numerical controls for the direct-image XLB solver.

Attributes:

Name	Type	Description
formulation	str	Either "incompressible_navier_stokes" or "steady_stokes_limit".
backend	str	XLB compute backend. The current voids adapter supports only "jax".
precision_policy	str	XLB precision policy name, for example "FP32FP32".
collision_model	str	XLB collision operator label passed to the stepper.
streaming_scheme	str	XLB streaming scheme label passed to the stepper.
lattice_viscosity	float	Kinematic viscosity in lattice units.
pressure_inlet_lattice, pressure_outlet_lattice		Optional inlet and outlet lattice pressures. If both are provided they define the pressure BC directly.
pressure_drop_lattice	float | None	Optional lattice pressure drop. When set without explicit inlet/outlet pressures, it is applied relative to reference_density_lattice.
reference_density_lattice	float	Reference lattice density used to construct a baseline outlet pressure when only pressure_drop_lattice is provided.
rho_inlet, rho_outlet		Legacy density-based BC inputs retained for backward compatibility. They are converted internally to lattice pressure using p_lu = c_s^2 rho.
inlet_outlet_buffer_cells	int	Number of fluid reservoir layers inserted ahead of and behind the sample.
max_steps, min_steps, check_interval, steady_rtol		Iteration and convergence controls for the steady-state solve.

Notes

The current voids adapter uses XLB's JAX backend only. This keeps the dependency path compatible with CPU-only macOS and Linux environments.

The currently exposed XLB operator is the incompressible Navier-Stokes lattice-Boltzmann stepper. Setting formulation="steady_stokes_limit" does not switch to a different PDE solver; it selects conservative forcing and convergence defaults so the converged solution can be interpreted in the steady creeping-flow limit.

In this isothermal LBM setting, lattice pressure satisfies p_lu = c_s^2 rho. The preferred public inputs are therefore the lattice pressure fields pressure_inlet_lattice / pressure_outlet_lattice or the pressure drop pressure_drop_lattice. The legacy fields rho_inlet and rho_outlet remain supported for backward compatibility and are converted internally to pressure.

Source code in src/voids/lbm/singlephase/xlb.py

@dataclass(slots=True)
class XLBOptions:
    """Numerical controls for the direct-image XLB solver.

    Attributes
    ----------
    formulation :
        Either ``"incompressible_navier_stokes"`` or
        ``"steady_stokes_limit"``.
    backend :
        XLB compute backend. The current `voids` adapter supports only
        ``"jax"``.
    precision_policy :
        XLB precision policy name, for example ``"FP32FP32"``.
    collision_model :
        XLB collision operator label passed to the stepper.
    streaming_scheme :
        XLB streaming scheme label passed to the stepper.
    lattice_viscosity :
        Kinematic viscosity in lattice units.
    pressure_inlet_lattice, pressure_outlet_lattice :
        Optional inlet and outlet lattice pressures. If both are provided they
        define the pressure BC directly.
    pressure_drop_lattice :
        Optional lattice pressure drop. When set without explicit inlet/outlet
        pressures, it is applied relative to ``reference_density_lattice``.
    reference_density_lattice :
        Reference lattice density used to construct a baseline outlet pressure
        when only ``pressure_drop_lattice`` is provided.
    rho_inlet, rho_outlet :
        Legacy density-based BC inputs retained for backward compatibility.
        They are converted internally to lattice pressure using
        ``p_lu = c_s^2 rho``.
    inlet_outlet_buffer_cells :
        Number of fluid reservoir layers inserted ahead of and behind the
        sample.
    max_steps, min_steps, check_interval, steady_rtol :
        Iteration and convergence controls for the steady-state solve.

    Notes
    -----
    The current `voids` adapter uses XLB's JAX backend only. This keeps the
    dependency path compatible with CPU-only macOS and Linux environments.

    The currently exposed XLB operator is the incompressible Navier-Stokes
    lattice-Boltzmann stepper. Setting ``formulation="steady_stokes_limit"``
    does not switch to a different PDE solver; it selects conservative forcing
    and convergence defaults so the converged solution can be interpreted in the
    steady creeping-flow limit.

    In this isothermal LBM setting, lattice pressure satisfies
    ``p_lu = c_s^2 rho``. The preferred public inputs are therefore the lattice
    pressure fields ``pressure_inlet_lattice`` / ``pressure_outlet_lattice`` or
    the pressure drop ``pressure_drop_lattice``. The legacy fields ``rho_inlet``
    and ``rho_outlet`` remain supported for backward compatibility and are
    converted internally to pressure.
    """

    formulation: str = "incompressible_navier_stokes"
    backend: str = "jax"
    precision_policy: str = "FP32FP32"
    collision_model: str = "BGK"
    streaming_scheme: str = "pull"
    lattice_viscosity: float = 0.10
    pressure_inlet_lattice: float | None = None
    pressure_outlet_lattice: float | None = None
    pressure_drop_lattice: float | None = DEFAULT_PRESSURE_DROP_LATTICE
    reference_density_lattice: float = DEFAULT_REFERENCE_DENSITY_LATTICE
    rho_inlet: float | None = None
    rho_outlet: float | None = None
    inlet_outlet_buffer_cells: int = 6
    max_steps: int = 2000
    min_steps: int = 200
    check_interval: int = 100
    steady_rtol: float = 1.0e-3

    @classmethod
    def steady_stokes_defaults(cls, **overrides: float | int | str) -> "XLBOptions":
        """Return a conservative preset for the steady creeping-flow limit.

        Notes
        -----
        XLB does not currently expose a separate Stokes-only stepper in the
        installed package used by `voids`. This preset therefore still uses the
        incompressible Navier-Stokes LBM operator, but with a smaller lattice
        pressure drop and tighter steady-state controls so the converged solution is
        interpreted in the low-Reynolds, low-Mach limit.

        The buffer and convergence controls are intentionally stricter than the
        generic :class:`XLBOptions` defaults. They were selected from same-ROI
        DRP-317 sensitivity runs as a conservative direct-image permeability
        preset, not as a fit to experimental permeability.
        """

        values: dict[str, Any] = {
            "formulation": "steady_stokes_limit",
            "lattice_viscosity": 0.10,
            "pressure_drop_lattice": DEFAULT_STOKES_PRESSURE_DROP_LATTICE,
            "inlet_outlet_buffer_cells": 12,
            "max_steps": 8000,
            "min_steps": 1200,
            "check_interval": 100,
            "steady_rtol": 1.0e-4,
        }
        values.update(overrides)
        return cls(**values)

steady_stokes_defaults classmethod

steady_stokes_defaults(**overrides)

Return a conservative preset for the steady creeping-flow limit.

Notes

XLB does not currently expose a separate Stokes-only stepper in the installed package used by voids. This preset therefore still uses the incompressible Navier-Stokes LBM operator, but with a smaller lattice pressure drop and tighter steady-state controls so the converged solution is interpreted in the low-Reynolds, low-Mach limit.

The buffer and convergence controls are intentionally stricter than the generic :class:XLBOptions defaults. They were selected from same-ROI DRP-317 sensitivity runs as a conservative direct-image permeability preset, not as a fit to experimental permeability.

Source code in src/voids/lbm/singlephase/xlb.py

@classmethod
def steady_stokes_defaults(cls, **overrides: float | int | str) -> "XLBOptions":
    """Return a conservative preset for the steady creeping-flow limit.

    Notes
    -----
    XLB does not currently expose a separate Stokes-only stepper in the
    installed package used by `voids`. This preset therefore still uses the
    incompressible Navier-Stokes LBM operator, but with a smaller lattice
    pressure drop and tighter steady-state controls so the converged solution is
    interpreted in the low-Reynolds, low-Mach limit.

    The buffer and convergence controls are intentionally stricter than the
    generic :class:`XLBOptions` defaults. They were selected from same-ROI
    DRP-317 sensitivity runs as a conservative direct-image permeability
    preset, not as a fit to experimental permeability.
    """

    values: dict[str, Any] = {
        "formulation": "steady_stokes_limit",
        "lattice_viscosity": 0.10,
        "pressure_drop_lattice": DEFAULT_STOKES_PRESSURE_DROP_LATTICE,
        "inlet_outlet_buffer_cells": 12,
        "max_steps": 8000,
        "min_steps": 1200,
        "check_interval": 100,
        "steady_rtol": 1.0e-4,
    }
    values.update(overrides)
    return cls(**values)

XLBDirectSimulationResult dataclass

Store direct-image LBM outputs from an XLB run.

Attributes:

Name	Type	Description
lattice_pressure_inlet, lattice_pressure_outlet, lattice_pressure_drop		Resolved inlet, outlet, and differential pressure in lattice units.
lattice_density_inlet, lattice_density_outlet		Equivalent lattice densities associated with the pressure BCs through p_lu = c_s^2 rho.
permeability	float	Apparent permeability mapped back to physical units.
max_mach_lattice, reynolds_voxel_max		Low-inertia diagnostics useful when interpreting a run as a creeping-flow reference.

Source code in src/voids/lbm/singlephase/xlb.py

@dataclass(slots=True)
class XLBDirectSimulationResult:
    """Store direct-image LBM outputs from an XLB run.

    Attributes
    ----------
    lattice_pressure_inlet, lattice_pressure_outlet, lattice_pressure_drop :
        Resolved inlet, outlet, and differential pressure in lattice units.
    lattice_density_inlet, lattice_density_outlet :
        Equivalent lattice densities associated with the pressure BCs through
        ``p_lu = c_s^2 rho``.
    permeability :
        Apparent permeability mapped back to physical units.
    max_mach_lattice, reynolds_voxel_max :
        Low-inertia diagnostics useful when interpreting a run as a creeping-flow
        reference.
    """

    flow_axis: str
    voxel_size: float
    image_porosity: float
    sample_lengths: dict[str, float]
    sample_cross_sections: dict[str, float]
    lattice_viscosity: float
    lattice_pressure_inlet: float
    lattice_pressure_outlet: float
    lattice_density_inlet: float
    lattice_density_outlet: float
    lattice_pressure_drop: float
    inlet_outlet_buffer_cells: int
    omega: float
    superficial_velocity_lattice: float
    superficial_velocity_profile_lattice: np.ndarray
    velocity_lattice: np.ndarray
    axial_velocity_lattice: np.ndarray
    converged: bool
    n_steps: int
    convergence_metric: float
    permeability: float
    backend: str
    backend_version: str | None
    formulation: str
    velocity_set: str
    collision_model: str
    streaming_scheme: str
    max_speed_lattice: float
    max_mach_lattice: float
    reynolds_voxel_max: float

solve_binary_volume_with_xlb

solve_binary_volume_with_xlb(
    phases, *, voxel_size, flow_axis=None, options=None
)

Solve a binary segmented volume directly with XLB and estimate permeability.

Parameters:

Name	Type	Description	Default
phases	ndarray	Binary segmented volume with void=1 and solid=0.	required
voxel_size	float	Physical voxel size used when mapping lattice permeability back to physical units.	required
flow_axis	str | None	Requested flow axis. When omitted, the canonical sample axis inferred by :func:voids.image.network_extraction.infer_sample_axes is used.	None
options	XLBOptions | None	XLB numerical controls. This low-level interface expects lattice-unit pressure inputs through pressure_inlet_lattice, pressure_outlet_lattice, or pressure_drop_lattice.	None

Returns:

Type	Description
XLBDirectSimulationResult	Direct-image XLB solution summary and permeability diagnostics.

Raises:

Type	Description
ValueError	If the input volume, flow axis, or XLB numerical controls are invalid.
RuntimeError	If the run completes with a non-physical permeability estimate.

Warns:

Type	Description
XLBConvergenceWarning	If max_steps is reached before the steady-state velocity criterion is satisfied. The result is still returned with converged=False and the final n_steps / convergence_metric diagnostics.

Notes

The current adapter uses XLB's incompressible Navier-Stokes lattice Boltzmann stepper with BGK-style collision and pressure / bounce-back boundary conditions. If options.formulation == "steady_stokes_limit", the same stepper is still used; the solution is simply driven more gently and interpreted in the low-Reynolds, low-Mach steady limit.

The inlet and outlet conditions are pressure boundary conditions. The public API accepts lattice pressure values or a lattice pressure drop, which are converted internally to the density values expected by the isothermal LBM boundary operator through p_lu = c_s^2 rho.

The permeability conversion is based on lattice units:

K_phys = nu_lu * U_lu * L_lu * dx_phys**2 / dp_lu

where U_lu is the superficial sample velocity, nu_lu is the lattice kinematic viscosity, L_lu is the voxel count along the flow axis, and dp_lu = p_in,lu - p_out,lu.

This keeps the solve focused on permeability, which is the transport quantity comparable to the PNM solve without requiring a full physical pressure-unit calibration of the lattice simulation.

Source code in src/voids/lbm/singlephase/xlb.py

def solve_binary_volume_with_xlb(
    phases: np.ndarray,
    *,
    voxel_size: float,
    flow_axis: str | None = None,
    options: XLBOptions | None = None,
) -> XLBDirectSimulationResult:
    """Solve a binary segmented volume directly with XLB and estimate permeability.

    Parameters
    ----------
    phases :
        Binary segmented volume with ``void=1`` and ``solid=0``.
    voxel_size :
        Physical voxel size used when mapping lattice permeability back to
        physical units.
    flow_axis :
        Requested flow axis. When omitted, the canonical sample axis inferred by
        :func:`voids.image.network_extraction.infer_sample_axes` is used.
    options :
        XLB numerical controls. This low-level interface expects lattice-unit
        pressure inputs through ``pressure_inlet_lattice``,
        ``pressure_outlet_lattice``, or ``pressure_drop_lattice``.

    Returns
    -------
    XLBDirectSimulationResult
        Direct-image XLB solution summary and permeability diagnostics.

    Raises
    ------
    ValueError
        If the input volume, flow axis, or XLB numerical controls are invalid.
    RuntimeError
        If the run completes with a non-physical permeability estimate.

    Warns
    -----
    XLBConvergenceWarning
        If ``max_steps`` is reached before the steady-state velocity criterion
        is satisfied. The result is still returned with ``converged=False`` and
        the final ``n_steps`` / ``convergence_metric`` diagnostics.

    Notes
    -----
    The current adapter uses XLB's incompressible Navier-Stokes lattice
    Boltzmann stepper with BGK-style collision and pressure / bounce-back
    boundary conditions. If ``options.formulation == "steady_stokes_limit"``,
    the same stepper is still used; the solution is simply driven more gently
    and interpreted in the low-Reynolds, low-Mach steady limit.

    The inlet and outlet conditions are pressure boundary conditions. The public
    API accepts lattice pressure values or a lattice pressure drop, which are
    converted internally to the density values expected by the isothermal LBM
    boundary operator through ``p_lu = c_s^2 rho``.

    The permeability conversion is based on lattice units:

    ``K_phys = nu_lu * U_lu * L_lu * dx_phys**2 / dp_lu``

    where ``U_lu`` is the superficial sample velocity, ``nu_lu`` is the lattice
    kinematic viscosity, ``L_lu`` is the voxel count along the flow axis, and
    ``dp_lu = p_in,lu - p_out,lu``.

    This keeps the solve focused on permeability, which is the transport
    quantity comparable to the PNM solve without requiring a full physical
    pressure-unit calibration of the lattice simulation.
    """

    arr = _as_binary_volume(phases)
    xlb_options = options or XLBOptions()

    if xlb_options.backend.lower() != "jax":
        raise ValueError("The current `voids` XLB adapter supports only backend='jax'")
    if xlb_options.formulation not in {"incompressible_navier_stokes", "steady_stokes_limit"}:
        raise ValueError(
            "formulation must be 'incompressible_navier_stokes' or 'steady_stokes_limit'"
        )
    if xlb_options.max_steps <= 0:
        raise ValueError("max_steps must be positive")
    if xlb_options.min_steps < 0:
        raise ValueError("min_steps must be non-negative")
    if xlb_options.check_interval <= 0:
        raise ValueError("check_interval must be positive")
    if xlb_options.steady_rtol <= 0:
        raise ValueError("steady_rtol must be positive")
    if xlb_options.lattice_viscosity <= 0:
        raise ValueError("lattice_viscosity must be positive")
    if xlb_options.inlet_outlet_buffer_cells < 0:
        raise ValueError("inlet_outlet_buffer_cells must be non-negative")
    pressure_inlet_lattice, pressure_outlet_lattice = _resolve_lattice_pressure_bc(
        xlb_options,
        cs2=ISOTHERMAL_LATTICE_CS2,
    )
    density_inlet_lattice = pressure_inlet_lattice / ISOTHERMAL_LATTICE_CS2
    density_outlet_lattice = pressure_outlet_lattice / ISOTHERMAL_LATTICE_CS2
    density_drop_lattice = density_inlet_lattice - density_outlet_lattice
    if xlb_options.formulation == "steady_stokes_limit" and density_drop_lattice > 5.0e-4:
        warnings.warn(
            "steady_stokes_limit is being used with a relatively large lattice "
            "pressure drop. The run still uses the same incompressible LBM stepper, "
            "so keep the imposed driving small if the goal is a steady Stokes-limit "
            "reference.",
            RuntimeWarning,
            stacklevel=2,
        )

    _, axis_lengths, axis_areas, inferred_axis = infer_sample_axes(arr.shape, voxel_size=voxel_size)
    axis = inferred_axis if flow_axis is None else flow_axis
    axis_index = _axis_to_index(axis, arr.ndim)

    aligned_void_sample = np.moveaxis(np.asarray(arr, dtype=bool), axis_index, 0)
    if not np.any(aligned_void_sample[0]):
        raise ValueError("The inlet plane contains no void voxels for the requested flow axis")
    if not np.any(aligned_void_sample[-1]):
        raise ValueError("The outlet plane contains no void voxels for the requested flow axis")

    buffer_cells = int(xlb_options.inlet_outlet_buffer_cells)
    if buffer_cells > 0:
        aligned_void = np.pad(
            aligned_void_sample,
            pad_width=((buffer_cells, buffer_cells),) + ((0, 0),) * (aligned_void_sample.ndim - 1),
            mode="constant",
            constant_values=True,
        )
    else:
        aligned_void = aligned_void_sample

    xlb_api = _import_xlb()
    jax = xlb_api["jax"]
    xlb = xlb_api["xlb"]
    ComputeBackend = xlb_api["ComputeBackend"]
    PrecisionPolicy = xlb_api["PrecisionPolicy"]
    IncompressibleNavierStokesStepper = xlb_api["IncompressibleNavierStokesStepper"]
    HalfwayBounceBackBC = xlb_api["HalfwayBounceBackBC"]
    RegularizedBC = xlb_api["RegularizedBC"]
    velocity_set_cls = xlb_api["D2Q9"] if arr.ndim == 2 else xlb_api["D3Q19"]

    precision_policy = getattr(PrecisionPolicy, xlb_options.precision_policy, None)
    if precision_policy is None:
        raise ValueError(f"Unknown XLB precision policy {xlb_options.precision_policy!r}")

    compute_backend = ComputeBackend.JAX
    velocity_set = velocity_set_cls(precision_policy, compute_backend)
    xlb.init(velocity_set, compute_backend, precision_policy)
    grid = xlb.grid.grid_factory(aligned_void.shape, compute_backend=compute_backend)

    sealed_side_mask = np.zeros_like(aligned_void, dtype=bool)
    for side_axis in range(1, aligned_void.ndim):
        lower: list[slice | int] = [slice(None)] * aligned_void.ndim
        upper: list[slice | int] = [slice(None)] * aligned_void.ndim
        lower[side_axis] = 0
        upper[side_axis] = -1
        sealed_side_mask[tuple(lower)] = True
        sealed_side_mask[tuple(upper)] = True

    inlet_mask = np.zeros_like(aligned_void, dtype=bool)
    outlet_mask = np.zeros_like(aligned_void, dtype=bool)
    # Pressure BCs are imposed on planar reservoir faces restricted to void voxels,
    # with side-wall edges excluded because those cells belong to the sealed sample
    # jacket.  Intersecting with ``aligned_void`` prevents solid voxels from being
    # assigned a pressure BC, which would otherwise "open" them and corrupt the
    # bounce-back mask assignment below.
    inlet_mask[0, ...] = aligned_void[0, ...] & ~sealed_side_mask[0, ...]
    outlet_mask[-1, ...] = aligned_void[-1, ...] & ~sealed_side_mask[-1, ...]
    if not np.any(inlet_mask):
        raise ValueError(
            "The trimmed inlet plane has no interior void voxels for the requested flow axis"
        )
    if not np.any(outlet_mask):
        raise ValueError(
            "The trimmed outlet plane has no interior void voxels for the requested flow axis"
        )

    bounceback_mask = ~aligned_void
    bounceback_mask |= sealed_side_mask
    bounceback_mask &= ~(inlet_mask | outlet_mask)

    inlet_indices = _mask_to_indices(inlet_mask)
    outlet_indices = _mask_to_indices(outlet_mask)
    bounceback_indices = _mask_to_indices(bounceback_mask)

    boundary_conditions = [
        RegularizedBC(
            "pressure",
            prescribed_value=float(pressure_inlet_lattice / float(np.asarray(velocity_set.cs2))),
            indices=inlet_indices,
        ),
        RegularizedBC(
            "pressure",
            prescribed_value=float(pressure_outlet_lattice / float(np.asarray(velocity_set.cs2))),
            indices=outlet_indices,
        ),
    ]
    if bounceback_indices is not None:
        boundary_conditions.append(HalfwayBounceBackBC(indices=bounceback_indices))

    stepper = IncompressibleNavierStokesStepper(
        grid=grid,
        boundary_conditions=boundary_conditions,
        collision_type=xlb_options.collision_model,
        streaming_scheme=xlb_options.streaming_scheme,
    )
    f_0, f_1, bc_mask, missing_mask = stepper.prepare_fields()

    omega_float = 1.0 / (3.0 * float(xlb_options.lattice_viscosity) + 0.5)
    omega = np.asarray(omega_float, dtype=precision_policy.compute_precision.jax_dtype)

    velocity_aligned = np.zeros((arr.ndim, *aligned_void_sample.shape), dtype=float)
    axial_velocity_aligned = np.zeros_like(aligned_void, dtype=float)
    superficial_profile = np.zeros(aligned_void.shape[0], dtype=float)
    superficial_velocity = 0.0
    max_speed_lattice = 0.0
    convergence_metric = np.inf
    previous_superficial_velocity: float | None = None
    converged = False
    n_steps = 0

    # f_current is a JAX array; typed as `object` because `jaxlib` types are
    # not available in the default mypy environment (jax is an optional dep).
    def _measure_current_state(
        f_current: object,
    ) -> tuple[np.ndarray, np.ndarray, np.ndarray, float, float]:
        jax.block_until_ready(f_current)
        _, u = stepper.macroscopic(f_current)
        u_full = np.asarray(u, dtype=float)
        axial_full = np.asarray(u_full[0], dtype=float)
        speed_full = np.linalg.norm(u_full, axis=0)
        sample_slice = slice(buffer_cells, buffer_cells + aligned_void_sample.shape[0])
        velocity = np.asarray(u_full[:, sample_slice, ...], dtype=float)
        axial = axial_full[sample_slice, ...]
        speed = speed_full[sample_slice, ...]
        profile = _superficial_velocity_profile(axial, aligned_void_sample)
        interior_profile = profile[1:-1] if profile.size > 2 else profile
        mean_superficial_velocity = float(np.mean(interior_profile))
        max_speed = float(np.max(speed[aligned_void_sample]))
        return velocity, axial, profile, mean_superficial_velocity, max_speed

    for step in range(xlb_options.max_steps):
        f_0, f_1 = stepper(f_0, f_1, bc_mask, missing_mask, omega, step)
        f_0, f_1 = f_1, f_0
        n_steps = step + 1

        if n_steps % xlb_options.check_interval != 0 and n_steps != xlb_options.max_steps:
            continue

        (
            velocity_aligned,
            axial_velocity_aligned,
            superficial_profile,
            superficial_velocity,
            max_speed_lattice,
        ) = _measure_current_state(f_0)
        if previous_superficial_velocity is not None:
            convergence_metric = abs(superficial_velocity - previous_superficial_velocity) / max(
                abs(previous_superficial_velocity),
                1.0e-30,
            )
        previous_superficial_velocity = superficial_velocity

        if (
            n_steps >= xlb_options.min_steps
            and np.isfinite(convergence_metric)
            and convergence_metric < xlb_options.steady_rtol
        ):
            converged = True
            break

    if not np.any(axial_velocity_aligned):
        (
            velocity_aligned,
            axial_velocity_aligned,
            superficial_profile,
            superficial_velocity,
            max_speed_lattice,
        ) = _measure_current_state(f_0)
        if previous_superficial_velocity is not None:
            convergence_metric = abs(superficial_velocity - previous_superficial_velocity) / max(
                abs(previous_superficial_velocity),
                1.0e-30,
            )

    lattice_pressure_drop = float(pressure_inlet_lattice) - float(pressure_outlet_lattice)
    permeability = (
        float(xlb_options.lattice_viscosity)
        * float(superficial_velocity)
        * float(axis_lengths[axis])
        * float(voxel_size)
        / float(lattice_pressure_drop)
    )
    cs_lattice = float(np.sqrt(np.asarray(velocity_set.cs2)))
    max_mach_lattice = float(max_speed_lattice) / cs_lattice
    reynolds_voxel_max = float(max_speed_lattice) / float(xlb_options.lattice_viscosity)
    if not converged:
        warnings.warn(
            "XLB direct-image solve did not satisfy the steady-state tolerance "
            f"after {n_steps} steps; the reported permeability may be biased. "
            f"Last relative velocity change: {convergence_metric:.3e}.",
            XLBConvergenceWarning,
            stacklevel=2,
        )
    if not np.isfinite(permeability) or permeability <= 0.0:
        raise RuntimeError(
            "XLB produced a non-physical permeability estimate "
            f"({permeability:.6e} m^2-equivalent). "
            "This usually indicates incompatible boundary conditions, an "
            "insufficient inlet/outlet buffer, or a run that is too short to "
            "reach steady state."
        )

    velocity_spatial_original = np.moveaxis(velocity_aligned, 1, axis_index + 1)
    aligned_axes = [axis_index, *[idx for idx in range(arr.ndim) if idx != axis_index]]
    velocity_original = np.empty_like(velocity_spatial_original)
    for aligned_component, original_component in enumerate(aligned_axes):
        velocity_original[original_component] = velocity_spatial_original[aligned_component]

    axial_velocity_original = np.asarray(velocity_original[axis_index], dtype=float)

    return XLBDirectSimulationResult(
        flow_axis=axis,
        voxel_size=float(voxel_size),
        image_porosity=float(arr.mean()),
        sample_lengths=axis_lengths,
        sample_cross_sections=axis_areas,
        lattice_viscosity=float(xlb_options.lattice_viscosity),
        lattice_pressure_inlet=float(pressure_inlet_lattice),
        lattice_pressure_outlet=float(pressure_outlet_lattice),
        lattice_density_inlet=float(pressure_inlet_lattice / float(np.asarray(velocity_set.cs2))),
        lattice_density_outlet=float(pressure_outlet_lattice / float(np.asarray(velocity_set.cs2))),
        lattice_pressure_drop=float(lattice_pressure_drop),
        inlet_outlet_buffer_cells=buffer_cells,
        omega=omega_float,
        superficial_velocity_lattice=float(superficial_velocity),
        superficial_velocity_profile_lattice=np.asarray(superficial_profile, dtype=float),
        velocity_lattice=np.asarray(velocity_original, dtype=float),
        axial_velocity_lattice=np.asarray(axial_velocity_original, dtype=float),
        converged=bool(converged),
        n_steps=int(n_steps),
        convergence_metric=float(convergence_metric),
        permeability=float(permeability),
        backend="jax",
        backend_version=getattr(xlb, "__version__", None),
        formulation=str(xlb_options.formulation),
        velocity_set=str(velocity_set_cls.__name__),
        collision_model=str(xlb_options.collision_model),
        streaming_scheme=str(xlb_options.streaming_scheme),
        max_speed_lattice=float(max_speed_lattice),
        max_mach_lattice=float(max_mach_lattice),
        reynolds_voxel_max=float(reynolds_voxel_max),
    )

---

Stokes-Limit XLB Backend

voids.lbm.singlephase.stokes

XLBConvergenceWarning

Bases: RuntimeWarning

Warn that an XLB solve reached max_steps before steady convergence.

Source code in src/voids/lbm/singlephase/xlb.py

class XLBConvergenceWarning(RuntimeWarning):
    """Warn that an XLB solve reached ``max_steps`` before steady convergence."""

XLBDirectSimulationResult dataclass

Store direct-image LBM outputs from an XLB run.

Attributes:

Name	Type	Description
lattice_pressure_inlet, lattice_pressure_outlet, lattice_pressure_drop		Resolved inlet, outlet, and differential pressure in lattice units.
lattice_density_inlet, lattice_density_outlet		Equivalent lattice densities associated with the pressure BCs through p_lu = c_s^2 rho.
permeability	float	Apparent permeability mapped back to physical units.
max_mach_lattice, reynolds_voxel_max		Low-inertia diagnostics useful when interpreting a run as a creeping-flow reference.

Source code in src/voids/lbm/singlephase/xlb.py

@dataclass(slots=True)
class XLBDirectSimulationResult:
    """Store direct-image LBM outputs from an XLB run.

    Attributes
    ----------
    lattice_pressure_inlet, lattice_pressure_outlet, lattice_pressure_drop :
        Resolved inlet, outlet, and differential pressure in lattice units.
    lattice_density_inlet, lattice_density_outlet :
        Equivalent lattice densities associated with the pressure BCs through
        ``p_lu = c_s^2 rho``.
    permeability :
        Apparent permeability mapped back to physical units.
    max_mach_lattice, reynolds_voxel_max :
        Low-inertia diagnostics useful when interpreting a run as a creeping-flow
        reference.
    """

    flow_axis: str
    voxel_size: float
    image_porosity: float
    sample_lengths: dict[str, float]
    sample_cross_sections: dict[str, float]
    lattice_viscosity: float
    lattice_pressure_inlet: float
    lattice_pressure_outlet: float
    lattice_density_inlet: float
    lattice_density_outlet: float
    lattice_pressure_drop: float
    inlet_outlet_buffer_cells: int
    omega: float
    superficial_velocity_lattice: float
    superficial_velocity_profile_lattice: np.ndarray
    velocity_lattice: np.ndarray
    axial_velocity_lattice: np.ndarray
    converged: bool
    n_steps: int
    convergence_metric: float
    permeability: float
    backend: str
    backend_version: str | None
    formulation: str
    velocity_set: str
    collision_model: str
    streaming_scheme: str
    max_speed_lattice: float
    max_mach_lattice: float
    reynolds_voxel_max: float

XLBOptions dataclass

Numerical controls for the direct-image XLB solver.

Attributes:

Name	Type	Description
formulation	str	Either "incompressible_navier_stokes" or "steady_stokes_limit".
backend	str	XLB compute backend. The current voids adapter supports only "jax".
precision_policy	str	XLB precision policy name, for example "FP32FP32".
collision_model	str	XLB collision operator label passed to the stepper.
streaming_scheme	str	XLB streaming scheme label passed to the stepper.
lattice_viscosity	float	Kinematic viscosity in lattice units.
pressure_inlet_lattice, pressure_outlet_lattice		Optional inlet and outlet lattice pressures. If both are provided they define the pressure BC directly.
pressure_drop_lattice	float | None	Optional lattice pressure drop. When set without explicit inlet/outlet pressures, it is applied relative to reference_density_lattice.
reference_density_lattice	float	Reference lattice density used to construct a baseline outlet pressure when only pressure_drop_lattice is provided.
rho_inlet, rho_outlet		Legacy density-based BC inputs retained for backward compatibility. They are converted internally to lattice pressure using p_lu = c_s^2 rho.
inlet_outlet_buffer_cells	int	Number of fluid reservoir layers inserted ahead of and behind the sample.
max_steps, min_steps, check_interval, steady_rtol		Iteration and convergence controls for the steady-state solve.

Notes

The current voids adapter uses XLB's JAX backend only. This keeps the dependency path compatible with CPU-only macOS and Linux environments.

The currently exposed XLB operator is the incompressible Navier-Stokes lattice-Boltzmann stepper. Setting formulation="steady_stokes_limit" does not switch to a different PDE solver; it selects conservative forcing and convergence defaults so the converged solution can be interpreted in the steady creeping-flow limit.

In this isothermal LBM setting, lattice pressure satisfies p_lu = c_s^2 rho. The preferred public inputs are therefore the lattice pressure fields pressure_inlet_lattice / pressure_outlet_lattice or the pressure drop pressure_drop_lattice. The legacy fields rho_inlet and rho_outlet remain supported for backward compatibility and are converted internally to pressure.

Source code in src/voids/lbm/singlephase/xlb.py

@dataclass(slots=True)
class XLBOptions:
    """Numerical controls for the direct-image XLB solver.

    Attributes
    ----------
    formulation :
        Either ``"incompressible_navier_stokes"`` or
        ``"steady_stokes_limit"``.
    backend :
        XLB compute backend. The current `voids` adapter supports only
        ``"jax"``.
    precision_policy :
        XLB precision policy name, for example ``"FP32FP32"``.
    collision_model :
        XLB collision operator label passed to the stepper.
    streaming_scheme :
        XLB streaming scheme label passed to the stepper.
    lattice_viscosity :
        Kinematic viscosity in lattice units.
    pressure_inlet_lattice, pressure_outlet_lattice :
        Optional inlet and outlet lattice pressures. If both are provided they
        define the pressure BC directly.
    pressure_drop_lattice :
        Optional lattice pressure drop. When set without explicit inlet/outlet
        pressures, it is applied relative to ``reference_density_lattice``.
    reference_density_lattice :
        Reference lattice density used to construct a baseline outlet pressure
        when only ``pressure_drop_lattice`` is provided.
    rho_inlet, rho_outlet :
        Legacy density-based BC inputs retained for backward compatibility.
        They are converted internally to lattice pressure using
        ``p_lu = c_s^2 rho``.
    inlet_outlet_buffer_cells :
        Number of fluid reservoir layers inserted ahead of and behind the
        sample.
    max_steps, min_steps, check_interval, steady_rtol :
        Iteration and convergence controls for the steady-state solve.

    Notes
    -----
    The current `voids` adapter uses XLB's JAX backend only. This keeps the
    dependency path compatible with CPU-only macOS and Linux environments.

    The currently exposed XLB operator is the incompressible Navier-Stokes
    lattice-Boltzmann stepper. Setting ``formulation="steady_stokes_limit"``
    does not switch to a different PDE solver; it selects conservative forcing
    and convergence defaults so the converged solution can be interpreted in the
    steady creeping-flow limit.

    In this isothermal LBM setting, lattice pressure satisfies
    ``p_lu = c_s^2 rho``. The preferred public inputs are therefore the lattice
    pressure fields ``pressure_inlet_lattice`` / ``pressure_outlet_lattice`` or
    the pressure drop ``pressure_drop_lattice``. The legacy fields ``rho_inlet``
    and ``rho_outlet`` remain supported for backward compatibility and are
    converted internally to pressure.
    """

    formulation: str = "incompressible_navier_stokes"
    backend: str = "jax"
    precision_policy: str = "FP32FP32"
    collision_model: str = "BGK"
    streaming_scheme: str = "pull"
    lattice_viscosity: float = 0.10
    pressure_inlet_lattice: float | None = None
    pressure_outlet_lattice: float | None = None
    pressure_drop_lattice: float | None = DEFAULT_PRESSURE_DROP_LATTICE
    reference_density_lattice: float = DEFAULT_REFERENCE_DENSITY_LATTICE
    rho_inlet: float | None = None
    rho_outlet: float | None = None
    inlet_outlet_buffer_cells: int = 6
    max_steps: int = 2000
    min_steps: int = 200
    check_interval: int = 100
    steady_rtol: float = 1.0e-3

    @classmethod
    def steady_stokes_defaults(cls, **overrides: float | int | str) -> "XLBOptions":
        """Return a conservative preset for the steady creeping-flow limit.

        Notes
        -----
        XLB does not currently expose a separate Stokes-only stepper in the
        installed package used by `voids`. This preset therefore still uses the
        incompressible Navier-Stokes LBM operator, but with a smaller lattice
        pressure drop and tighter steady-state controls so the converged solution is
        interpreted in the low-Reynolds, low-Mach limit.

        The buffer and convergence controls are intentionally stricter than the
        generic :class:`XLBOptions` defaults. They were selected from same-ROI
        DRP-317 sensitivity runs as a conservative direct-image permeability
        preset, not as a fit to experimental permeability.
        """

        values: dict[str, Any] = {
            "formulation": "steady_stokes_limit",
            "lattice_viscosity": 0.10,
            "pressure_drop_lattice": DEFAULT_STOKES_PRESSURE_DROP_LATTICE,
            "inlet_outlet_buffer_cells": 12,
            "max_steps": 8000,
            "min_steps": 1200,
            "check_interval": 100,
            "steady_rtol": 1.0e-4,
        }
        values.update(overrides)
        return cls(**values)

steady_stokes_defaults classmethod

steady_stokes_defaults(**overrides)

Return a conservative preset for the steady creeping-flow limit.

Notes

XLB does not currently expose a separate Stokes-only stepper in the installed package used by voids. This preset therefore still uses the incompressible Navier-Stokes LBM operator, but with a smaller lattice pressure drop and tighter steady-state controls so the converged solution is interpreted in the low-Reynolds, low-Mach limit.

The buffer and convergence controls are intentionally stricter than the generic :class:XLBOptions defaults. They were selected from same-ROI DRP-317 sensitivity runs as a conservative direct-image permeability preset, not as a fit to experimental permeability.

Source code in src/voids/lbm/singlephase/xlb.py

@classmethod
def steady_stokes_defaults(cls, **overrides: float | int | str) -> "XLBOptions":
    """Return a conservative preset for the steady creeping-flow limit.

    Notes
    -----
    XLB does not currently expose a separate Stokes-only stepper in the
    installed package used by `voids`. This preset therefore still uses the
    incompressible Navier-Stokes LBM operator, but with a smaller lattice
    pressure drop and tighter steady-state controls so the converged solution is
    interpreted in the low-Reynolds, low-Mach limit.

    The buffer and convergence controls are intentionally stricter than the
    generic :class:`XLBOptions` defaults. They were selected from same-ROI
    DRP-317 sensitivity runs as a conservative direct-image permeability
    preset, not as a fit to experimental permeability.
    """

    values: dict[str, Any] = {
        "formulation": "steady_stokes_limit",
        "lattice_viscosity": 0.10,
        "pressure_drop_lattice": DEFAULT_STOKES_PRESSURE_DROP_LATTICE,
        "inlet_outlet_buffer_cells": 12,
        "max_steps": 8000,
        "min_steps": 1200,
        "check_interval": 100,
        "steady_rtol": 1.0e-4,
    }
    values.update(overrides)
    return cls(**values)

steady_stokes_options

steady_stokes_options(**overrides)

Return XLB options configured for the steady Stokes-limit formulation.

Source code in src/voids/lbm/singlephase/stokes.py

def steady_stokes_options(**overrides: Any) -> XLBOptions:
    """Return XLB options configured for the steady Stokes-limit formulation."""

    return XLBOptions.steady_stokes_defaults(**overrides)

solve_binary_volume_stokes

solve_binary_volume_stokes(
    phases, *, voxel_size, flow_axis=None, options=None
)

Solve Stokes-limit flow on a binary image using the XLB backend.

This is the package-facing LBM namespace for direct-image single-phase solves. Benchmark utilities consume this lower-level adapter for verification workflows. The input follows the current XLB adapter convention: void=1 and solid=0.

Source code in src/voids/lbm/singlephase/stokes.py

def solve_binary_volume_stokes(
    phases: np.ndarray,
    *,
    voxel_size: float,
    flow_axis: str | None = None,
    options: XLBOptions | None = None,
) -> XLBDirectSimulationResult:
    """Solve Stokes-limit flow on a binary image using the XLB backend.

    This is the package-facing LBM namespace for direct-image single-phase
    solves. Benchmark utilities consume this lower-level adapter for
    verification workflows. The input follows the current XLB adapter
    convention: ``void=1`` and ``solid=0``.
    """

    return solve_binary_volume_with_xlb(
        phases,
        voxel_size=voxel_size,
        flow_axis=flow_axis,
        options=options or steady_stokes_options(),
    )