API Reference

This page contains the API reference for the poromics package.

Metrics (poromics)

Convenience pipelines for computing transport properties. These are the main entry points for most users.

poromics

Modules:

cli –
julia_helpers –
simulation –
utils –

Classes:

PermeabilityResult –

Results from an LBM flow / permeability simulation.
SimulationResult –

Base class for simulation results.
TortuosityResult –

Results from a tortuosity / effective diffusivity simulation.

Functions:

permeability_lbm –

Compute absolute permeability using LBM (D3Q19 MRT).
tortuosity_fd –

Compute tortuosity via Julia FD solver.
tortuosity_lbm –

Compute tortuosity and effective diffusivity using LBM (D3Q7 BGK).

PermeabilityResult

PermeabilityResult(
    im,
    axis,
    porosity,
    k,
    u_darcy,
    u_pore,
    velocity,
    kinematic_pressure,
    pressure=None,
    *,
    converged=True,
    n_iterations=None,
    _velocity_lu=None,
    _rho_lu=None,
    _rho_out=None,
    _k_lu=None,
    _u_darcy_lu=None,
    _u_pore_lu=None,
)

Bases: SimulationResult

Results from an LBM flow / permeability simulation.

Use rescale to reinterpret the simulation for a different voxel size, fluid viscosity, or fluid density without rerunning the solver.

Parameters:

im
(ndarray) –

The boolean image used in the simulation.
axis
(int) –

The axis along which flow was computed.
porosity
(float) –

Pore volume fraction.
k
(float) –

Permeability in m².
u_darcy
(float) –

Darcy (superficial) velocity in m/s.
u_pore
(float) –

Mean pore-space velocity in m/s.
velocity
((ndarray, shape(nx, ny, nz, 3))) –

Steady-state velocity field in m/s.
kinematic_pressure
((ndarray, shape(nx, ny, nz))) –

Gauge kinematic pressure (p / rho) in m²/s². Density-free; always populated.
pressure
((ndarray or None, shape(nx, ny, nz)), default: None ) –

Gauge pressure in Pa. Populated only when a fluid density was provided (rho in permeability_lbm or rescale).
converged
(bool, default: True ) –

Whether the solver reached the requested tolerance. False means the reported k / velocity are from a pre-steady-state field and should not be trusted without further iteration.
n_iterations
(int or None, default: None ) –

Iterations the solver took. None for non-iterative backends.

Methods:

plot_pressure –

Plot gauge pressure field on a 2D slice.
plot_velocity –

Plot velocity magnitude and streamlines on a 2D slice.
rescale –

Reinterpret the simulation for different physical parameters.

Source code in src/poromics/_metrics.py

def __init__(self, im, axis, porosity, k, u_darcy, u_pore,
             velocity, kinematic_pressure, pressure=None, *,
             converged=True, n_iterations=None,
             _velocity_lu=None, _rho_lu=None, _rho_out=None,
             _k_lu=None, _u_darcy_lu=None,
             _u_pore_lu=None):  # fmt: skip
    super().__init__(im, axis, porosity)
    self.k = k
    self.u_darcy = u_darcy
    self.u_pore = u_pore
    self.velocity = velocity
    self.kinematic_pressure = kinematic_pressure
    self.pressure = pressure
    self.converged = converged
    self.n_iterations = n_iterations
    self._velocity_lu = _velocity_lu
    self._rho_lu = _rho_lu
    self._rho_out = _rho_out
    self._k_lu = _k_lu
    self._u_darcy_lu = _u_darcy_lu
    self._u_pore_lu = _u_pore_lu

plot_pressure

plot_pressure(z=None, ax=None, kinematic=False)

Plot gauge pressure field on a 2D slice.

Parameters:

z
(int or None, default: None ) –

Slice index along the third axis. Defaults to nz // 2.
ax
(Axes or None, default: None ) –

If provided, plot on this axes. If None, create a new figure.
kinematic
(bool, default: False ) –

Plot kinematic_pressure (m²/s²) instead of pressure (Pa). Useful when no fluid density was provided.

Returns:

fig ( Figure ) –
ax ( Axes ) –

Source code in src/poromics/_metrics.py

def plot_pressure(self, z=None, ax=None, kinematic=False):
    """Plot gauge pressure field on a 2D slice.

    Parameters
    ----------
    z : int or None
        Slice index along the third axis. Defaults to nz // 2.
    ax : matplotlib.axes.Axes or None
        If provided, plot on this axes. If None, create a new
        figure.
    kinematic : bool
        Plot ``kinematic_pressure`` (m²/s²) instead of
        ``pressure`` (Pa). Useful when no fluid density was
        provided.

    Returns
    -------
    fig : matplotlib.figure.Figure
    ax : matplotlib.axes.Axes
    """
    import matplotlib.pyplot as plt

    if kinematic:
        field, title, label = (
            self.kinematic_pressure,
            "Kinematic pressure (m²/s²)",
            "p/ρ (m²/s²)",
        )
    else:
        if self.pressure is None:
            raise RuntimeError(
                "Pressure in Pa is unavailable because no fluid "
                "density was provided. Pass rho=... to "
                "permeability_lbm / rescale, or call with "
                "kinematic=True."
            )
        field, title, label = self.pressure, "Pressure field (Pa)", "P (Pa)"

    im = np.atleast_3d(self.im)
    P = np.atleast_3d(field)
    if z is None:
        z = P.shape[2] // 2
    p_slice = P[:, :, z].astype(float)
    p_slice[~im[:, :, z]] = np.nan

    if ax is None:
        fig, ax = plt.subplots(figsize=(4, 4))
    else:
        fig = ax.figure
    import matplotlib.ticker as ticker

    mappable = ax.imshow(p_slice, cmap="coolwarm", interpolation="nearest")
    ax.set_title(title)
    cb = fig.colorbar(mappable, ax=ax, fraction=0.046, pad=0.04, label=label)
    cb.ax.yaxis.set_major_formatter(ticker.ScalarFormatter(useMathText=True))
    cb.ax.ticklabel_format(style="scientific", scilimits=(-1, 1))
    plt.tight_layout()
    return fig, ax

plot_velocity

plot_velocity(z=None, ax=None)

Plot velocity magnitude and streamlines on a 2D slice.

Parameters:

z
(int or None, default: None ) –

Slice index along the third axis. Defaults to nz // 2.
ax
(Axes or None, default: None ) –

If provided, plot velocity magnitude only on this axes. If None, create a side-by-side figure (magnitude + streamlines).

Returns:

fig ( Figure ) –
axes ( ndarray of Axes (if ax is None) or Axes ) –

Source code in src/poromics/_metrics.py

def plot_velocity(self, z=None, ax=None):
    """Plot velocity magnitude and streamlines on a 2D slice.

    Parameters
    ----------
    z : int or None
        Slice index along the third axis. Defaults to nz // 2.
    ax : matplotlib.axes.Axes or None
        If provided, plot velocity magnitude only on this axes.
        If None, create a side-by-side figure (magnitude +
        streamlines).

    Returns
    -------
    fig : matplotlib.figure.Figure
    axes : ndarray of Axes (if ax is None) or Axes
    """
    import matplotlib.pyplot as plt

    im = np.atleast_3d(self.im)
    vel = np.atleast_3d(self.velocity)
    if z is None:
        z = im.shape[2] // 2
    v = vel[:, :, z, :]
    speed = np.linalg.norm(v, axis=-1)
    solid_mask = ~im[:, :, z]
    speed_masked = speed.copy()
    speed_masked[solid_mask] = np.nan

    import matplotlib.ticker as ticker

    if ax is not None:
        mappable = ax.imshow(speed_masked, cmap="viridis", interpolation="nearest")
        ax.set_title("Velocity magnitude")
        cb = ax.figure.colorbar(
            mappable, ax=ax, fraction=0.046, pad=0.04, label="|u| (m/s)"
        )
        cb.ax.yaxis.set_major_formatter(ticker.ScalarFormatter(useMathText=True))
        cb.ax.ticklabel_format(style="scientific", scilimits=(-1, 1))
        return ax.figure, ax

    from mpl_toolkits.axes_grid1 import make_axes_locatable

    fig, axes = plt.subplots(1, 2, figsize=(7.5, 4))
    mappable0 = axes[0].imshow(speed_masked, cmap="viridis", interpolation="nearest")
    axes[0].set_title("Velocity magnitude")
    div0 = make_axes_locatable(axes[0])
    cax0 = div0.append_axes("right", size="5%", pad=0.05)
    cb = fig.colorbar(mappable0, cax=cax0, label="|u| (m/s)")
    cb.ax.yaxis.set_major_formatter(ticker.ScalarFormatter(useMathText=True))
    cb.ax.ticklabel_format(style="scientific", scilimits=(-1, 1))

    nx, ny = im.shape[0], im.shape[1]
    X, Y = np.meshgrid(np.arange(ny), np.arange(nx))
    lw = 3 * speed / speed.max() + 0.1 if speed.max() > 0 else np.ones_like(speed)
    axes[1].imshow(im[:, :, z], cmap="gray", interpolation="nearest", alpha=0.3)
    axes[1].streamplot(
        X,
        Y,
        v[:, :, 1],
        v[:, :, 0],
        color=speed,
        cmap="viridis",
        density=2,
        linewidth=lw,
        arrowstyle="fancy",
    )
    axes[1].set_xlim(-0.5, ny - 0.5)
    axes[1].set_ylim(nx - 0.5, -0.5)
    axes[1].set_title("Velocity streamlines")
    div1 = make_axes_locatable(axes[1])
    phantom = div1.append_axes("right", size="5%", pad=0.05)
    phantom.set_visible(False)
    fig.tight_layout()
    return fig, axes

rescale

rescale(voxel_size, nu, rho=None)

Reinterpret the simulation for different physical parameters.

Returns a new PermeabilityResult with recomputed physical quantities. The underlying lattice-unit fields are unchanged, so no solver re-run is needed. This is valid because permeability depends only on geometry, and the Stokes equations are linear.

Parameters:

voxel_size
(float) –

Physical voxel edge length in metres.
nu
(float) –

Kinematic viscosity in m²/s.
rho
(float or None, default: None ) –

Fluid density in kg/m³. Required for pressure conversion to Pa. If None, the pressure field is omitted (set to None); kinematic_pressure is always populated.

Returns:

result ( PermeabilityResult ) –

Source code in src/poromics/_metrics.py

def rescale(self, voxel_size, nu, rho=None):
    """Reinterpret the simulation for different physical parameters.

    Returns a new ``PermeabilityResult`` with recomputed physical
    quantities. The underlying lattice-unit fields are unchanged,
    so no solver re-run is needed. This is valid because
    permeability depends only on geometry, and the Stokes equations
    are linear.

    Parameters
    ----------
    voxel_size : float
        Physical voxel edge length in metres.
    nu : float
        Kinematic viscosity in m²/s.
    rho : float or None
        Fluid density in kg/m³. Required for pressure conversion
        to Pa. If None, the ``pressure`` field is omitted (set to
        None); ``kinematic_pressure`` is always populated.

    Returns
    -------
    result : PermeabilityResult
    """
    if self._velocity_lu is None:
        raise RuntimeError("Lattice-unit data not available for rescaling.")
    dt = _d3q19.nu * voxel_size**2 / nu
    lu_to_phys = voxel_size / dt
    k = self._k_lu * voxel_size**2
    u_darcy = self._u_darcy_lu * lu_to_phys
    u_pore = self._u_pore_lu * lu_to_phys
    velocity = self._velocity_lu * lu_to_phys
    gauge_rho = self._rho_lu - self._rho_out
    kinematic_pressure = _d3q19.cs2 * lu_to_phys**2 * gauge_rho
    pressure = rho * kinematic_pressure if rho is not None else None
    return PermeabilityResult(
        im=self.im,
        axis=self.axis,
        porosity=self.porosity,
        k=k,
        u_darcy=u_darcy,
        u_pore=u_pore,
        velocity=velocity,
        kinematic_pressure=kinematic_pressure,
        pressure=pressure,
        converged=self.converged,
        n_iterations=self.n_iterations,
        _velocity_lu=self._velocity_lu,
        _rho_lu=self._rho_lu,
        _rho_out=self._rho_out,
        _k_lu=self._k_lu,
        _u_darcy_lu=self._u_darcy_lu,
        _u_pore_lu=self._u_pore_lu,
    )

SimulationResult

SimulationResult(im, axis, porosity)

Base class for simulation results.

Parameters:

im
(ndarray) –

The boolean image used in the simulation.
axis
(int) –

The axis along which the simulation was run.
porosity
(float) –

Pore volume fraction.

Source code in src/poromics/_metrics.py

def __init__(self, im, axis, porosity):
    self.im = im
    self.axis = axis
    self.porosity = porosity

TortuosityResult

TortuosityResult(
    im,
    axis,
    porosity,
    tau,
    D_eff,
    c,
    formation_factor=None,
    D=None,
    *,
    converged=True,
    n_iterations=None,
)

Bases: SimulationResult

Results from a tortuosity / effective diffusivity simulation.

Tortuosity, effective diffusivity, and the concentration field are all dimensionless (or normalized), so no rescale method is provided. See PermeabilityResult.rescale for the flow analogue.

Parameters:

im
(ndarray) –

The boolean image used in the simulation.
axis
(int) –

The axis along which diffusion was computed.
porosity
(float) –

Pore volume fraction.
tau
(float) –

Tortuosity factor (>= 1).
D_eff
(float) –

Normalized effective diffusivity (D_eff / D_0).
c
(ndarray) –

Steady-state concentration field.
formation_factor
(float, default: None ) –

Formation factor F = 1 / D_eff_norm.
D
(float or ndarray, default: None ) –

Bulk diffusivity (float for uniform, ndarray for spatially variable).
converged
(bool, default: True ) –

Whether the solver reached the requested tolerance. False means the reported tau / D_eff are from a pre-steady-state field and should not be trusted without further iteration.
n_iterations
(int or None, default: None ) –

Iterations the solver took. None for non-iterative backends.

Methods:

plot_concentration –

Plot concentration field on a 2D slice.

Source code in src/poromics/_metrics.py

def __init__(self, im, axis, porosity, tau, D_eff, c,
             formation_factor=None, D=None, *,
             converged=True, n_iterations=None):  # fmt: skip
    super().__init__(im, axis, porosity)
    self.tau = tau
    self.D_eff = D_eff
    self.c = c
    self.formation_factor = formation_factor
    self.D = D
    self.converged = converged
    self.n_iterations = n_iterations

plot_concentration

plot_concentration(z=None, ax=None)

Plot concentration field on a 2D slice.

Parameters:

z
(int or None, default: None ) –

Slice index along the third axis. Defaults to nz // 2.
ax
(Axes or None, default: None ) –

If provided, plot concentration only on this axes. If None, create a side-by-side figure (pore structure + concentration).

Returns:

fig ( Figure ) –
axes ( ndarray of Axes (if ax is None) or Axes ) –

Source code in src/poromics/_metrics.py

def plot_concentration(self, z=None, ax=None):
    """Plot concentration field on a 2D slice.

    Parameters
    ----------
    z : int or None
        Slice index along the third axis. Defaults to nz // 2.
    ax : matplotlib.axes.Axes or None
        If provided, plot concentration only on this axes.
        If None, create a side-by-side figure (pore structure +
        concentration).

    Returns
    -------
    fig : matplotlib.figure.Figure
    axes : ndarray of Axes (if ax is None) or Axes
    """
    import matplotlib.pyplot as plt

    im = np.atleast_3d(self.im)
    c = np.atleast_3d(self.c)
    if z is None:
        z = im.shape[2] // 2

    c_slice = c[:, :, z].astype(float)
    c_slice[~im[:, :, z]] = np.nan

    import matplotlib.ticker as ticker

    if ax is not None:
        mappable = ax.imshow(c_slice, cmap="viridis", interpolation="nearest")
        ax.set_title("Concentration field")
        cb = ax.figure.colorbar(mappable, ax=ax, fraction=0.046, pad=0.04, label="c")
        cb.ax.yaxis.set_major_formatter(ticker.ScalarFormatter(useMathText=True))
        cb.ax.ticklabel_format(style="scientific", scilimits=(-1, 1))
        return ax.figure, ax

    from mpl_toolkits.axes_grid1 import make_axes_locatable

    fig, axes = plt.subplots(1, 2, figsize=(7.5, 4))
    axes[0].imshow(im[:, :, z], cmap="gray", interpolation="nearest")
    axes[0].set_title("Pore structure")
    div0 = make_axes_locatable(axes[0])
    phantom = div0.append_axes("right", size="5%", pad=0.05)
    phantom.set_visible(False)

    mappable = axes[1].imshow(c_slice, cmap="viridis", interpolation="nearest")
    axes[1].set_title("Concentration field")
    div1 = make_axes_locatable(axes[1])
    cax1 = div1.append_axes("right", size="5%", pad=0.05)
    cb = fig.colorbar(mappable, cax=cax1, label="c")
    cb.ax.yaxis.set_major_formatter(ticker.ScalarFormatter(useMathText=True))
    cb.ax.ticklabel_format(style="scientific", scilimits=(-1, 1))
    fig.tight_layout()
    return fig, axes

permeability_lbm

permeability_lbm(
    im,
    *,
    axis: int,
    nu: float = 1e-06,
    rho: float | None = None,
    voxel_size: float,
    tol: float = 0.001,
    n_steps: int = 100000,
    sparse: bool = False,
    verbose: bool = True,
) -> PermeabilityResult

Compute absolute permeability using LBM (D3Q19 MRT).

Solves creeping (Stokes) flow on the pore space of a 3D binary image using the Lattice Boltzmann Method. Permeability is extracted via Darcy's law.

Parameters:

im
((ndarray, shape(nx, ny, nz))) –

Binary image. True (or 1) = pore, False (or 0) = solid.
axis
(int) –

Axis along which to apply the pressure gradient (0=x, 1=y, 2=z).
nu
(float, default: 1e-06 ) –

Kinematic viscosity in m²/s. Default 1e-6 (water at ~20 °C).
rho
(float or None, default: None ) –

Fluid density in kg/m³. Required to return pressure in pascals; if None, pressure is set to None and only kinematic_pressure (m²/s²) is populated. Permeability and velocity are independent of density (Stokes regime).
voxel_size
(float) –

Physical voxel edge length in metres.
tol
(float, default: 0.001 ) –

Convergence tolerance on relative velocity change.
n_steps
(int, default: 100000 ) –

Maximum number of LBM iterations.
sparse
(bool, default: False ) –

Use Taichi sparse storage.
verbose
(bool, default: True ) –

Show a rich progress bar. Default False.

Returns:

result ( PermeabilityResult ) –

Source code in src/poromics/_metrics.py

def permeability_lbm(
    im,
    *,
    axis: int,
    nu: float = 1e-6,
    rho: float | None = None,
    voxel_size: float,
    tol: float = 1e-3,
    n_steps: int = 100_000,
    sparse: bool = False,
    verbose: bool = True,
) -> PermeabilityResult:
    """Compute absolute permeability using LBM (D3Q19 MRT).

    Solves creeping (Stokes) flow on the pore space of a 3D binary image
    using the Lattice Boltzmann Method. Permeability is extracted via
    Darcy's law.

    Parameters
    ----------
    im : ndarray, shape (nx, ny, nz)
        Binary image. True (or 1) = pore, False (or 0) = solid.
    axis : int
        Axis along which to apply the pressure gradient
        (0=x, 1=y, 2=z).
    nu : float
        Kinematic viscosity in m²/s. Default 1e-6 (water at ~20 °C).
    rho : float or None
        Fluid density in kg/m³. Required to return ``pressure`` in
        pascals; if None, ``pressure`` is set to None and only
        ``kinematic_pressure`` (m²/s²) is populated. Permeability and
        velocity are independent of density (Stokes regime).
    voxel_size : float
        Physical voxel edge length in metres.
    tol : float
        Convergence tolerance on relative velocity change.
    n_steps : int
        Maximum number of LBM iterations.
    sparse : bool
        Use Taichi sparse storage.
    verbose : bool
        Show a rich progress bar. Default False.

    Returns
    -------
    result : PermeabilityResult
    """
    im = np.atleast_3d(np.asarray(im, dtype=bool))
    n_pore_before = int(im.sum())
    im = _trim_nonpercolating(im, axis)
    if im.sum() == 0:
        raise RuntimeError("No percolating paths along the given axis found in the image.")
    n_removed = n_pore_before - int(im.sum())
    if n_removed > 0:
        logger.warning(f"Trimmed {n_removed} non-percolating pore voxels from the image.")
    solver = TransientFlow(im, axis=axis, nu=nu, rho=rho,
                           voxel_size=voxel_size, sparse=sparse)  # fmt: skip
    solver.run(n_steps=n_steps, tol=tol, verbose=verbose)
    if not solver.converged:
        logger.warning(
            f"LBM flow solver did not converge after "
            f"{solver.n_iterations} steps (tol={tol:.2e}). "
            f"Reported k / velocity are from a pre-steady-state field; "
            f"increase n_steps or loosen tol."
        )

    # Work in lattice units for Darcy's law, then convert
    v_lu = solver._solver.get_velocity()
    rho_lu = solver._solver.get_density()
    pore_mask = im
    porosity = float(pore_mask.sum()) / pore_mask.size
    L = im.shape[axis]

    v_flow_lu = v_lu[..., axis]
    u_darcy_lu = float(np.mean(v_flow_lu))
    u_pore_lu = float(np.mean(v_flow_lu[pore_mask]))
    grad_P_lu = (solver._rho_in - solver._rho_out) * _d3q19.cs2 / L
    k_lu = u_darcy_lu * _d3q19.nu / grad_P_lu

    # Convert to physical units
    dx = voxel_size
    k = k_lu * dx**2
    lu_to_phys = dx / solver.dt
    u_darcy = u_darcy_lu * lu_to_phys
    u_pore = u_pore_lu * lu_to_phys

    kinematic_pressure = solver.kinematic_pressure
    pressure = solver.pressure if rho is not None else None

    return PermeabilityResult(
        im=im,
        axis=axis,
        porosity=porosity,
        k=k,
        u_darcy=u_darcy,
        u_pore=u_pore,
        velocity=solver.velocity,
        kinematic_pressure=kinematic_pressure,
        pressure=pressure,
        converged=solver.converged,
        n_iterations=solver.n_iterations,
        _velocity_lu=v_lu,
        _rho_lu=rho_lu,
        _rho_out=solver._rho_out,
        _k_lu=k_lu,
        _u_darcy_lu=u_darcy_lu,
        _u_pore_lu=u_pore_lu,
    )

tortuosity_fd

tortuosity_fd(
    im,
    *,
    axis: int,
    D: ndarray | None = None,
    rtol: float = 1e-05,
    gpu: bool = True,
    verbose: bool = True,
) -> TortuosityResult

Compute tortuosity via Julia FD solver.

By default (POROMICS_JULIA_SUBPROCESS=1), runs Julia in a persistent subprocess to avoid LLVM symbol collisions with Taichi. The worker stays alive so Julia's JIT cache is reused across calls.

Set POROMICS_JULIA_SUBPROCESS=0 to run Julia in-process (faster if Taichi is not used). An error is raised if Taichi was already initialized in the same process.

Parameters:

im
(ndarray) –

The input image.
axis
(int) –

The axis along which to compute tortuosity (0=x, 1=y, 2=z).
D
(ndarray, default: None ) –

Diffusivity field. If None, uniform diffusivity of 1.0 is assumed.
rtol
(float, default: 1e-05 ) –

Relative tolerance for the solver.
gpu
(bool, default: True ) –

If True (default), use GPU for computation. Poromics picks a backend based on the platform (Darwin/arm64 → Metal, Linux/x86_64 → CUDA, Windows/AMD64 → CUDA) and installs it lazily on first call. Set POROMICS_GPU_BACKEND={metal,cuda,amdgpu} to override. When no supported backend is available or the device isn't functional, falls back to CPU with a warning. Set to False to force CPU.
verbose
(bool, default: True ) –

If True, print additional information during the solution process.

Returns:

result ( TortuosityResult ) –

Raises:

RuntimeError –

If no percolating paths are found along the specified axis.

Source code in src/poromics/_metrics.py

def tortuosity_fd(
    im,
    *,
    axis: int,
    D: np.ndarray | None = None,
    rtol: float = 1e-5,
    gpu: bool = True,
    verbose: bool = True,
) -> "TortuosityResult":
    """Compute tortuosity via Julia FD solver.

    By default (``POROMICS_JULIA_SUBPROCESS=1``), runs Julia in a
    persistent subprocess to avoid LLVM symbol collisions with Taichi.
    The worker stays alive so Julia's JIT cache is reused across calls.

    Set ``POROMICS_JULIA_SUBPROCESS=0`` to run Julia in-process (faster
    if Taichi is not used). An error is raised if Taichi was already
    initialized in the same process.

    Parameters
    ----------
    im : ndarray
        The input image.
    axis : int
        The axis along which to compute tortuosity (0=x, 1=y, 2=z).
    D : ndarray, optional
        Diffusivity field. If None, uniform diffusivity of 1.0 is assumed.
    rtol : float
        Relative tolerance for the solver.
    gpu : bool
        If True (default), use GPU for computation. Poromics picks a backend
        based on the platform (Darwin/arm64 → Metal, Linux/x86_64 → CUDA,
        Windows/AMD64 → CUDA) and installs it lazily on first call. Set
        ``POROMICS_GPU_BACKEND={metal,cuda,amdgpu}`` to override. When no
        supported backend is available or the device isn't functional, falls
        back to CPU with a warning. Set to False to force CPU.
    verbose : bool
        If True, print additional information during the solution process.

    Returns
    -------
    result : TortuosityResult

    Raises
    ------
    RuntimeError
        If no percolating paths are found along the specified axis.
    """
    im = np.atleast_3d(np.asarray(im, dtype=bool))
    n_pore_before = int(im.sum())
    im = _trim_nonpercolating(im, axis)
    if im.sum() == 0:
        raise RuntimeError("No percolating paths along the given axis found in the image.")
    n_removed = n_pore_before - int(im.sum())
    if n_removed > 0:
        logger.warning(f"Trimmed {n_removed} non-percolating pore voxels from the image.")
        if D is not None:
            D = np.array(D, copy=True)
            D[~im] = 0.0
    if not _JULIA_SUBPROCESS:
        return _tortuosity_fd_inprocess(im, axis, D, rtol, gpu, verbose)
    payload = {
        "im": im,
        "axis": axis,
        "D": D,
        "rtol": rtol,
        "gpu": gpu,
        "verbose": verbose,
    }
    return TortuosityResult(**_julia_call(payload))

tortuosity_lbm

tortuosity_lbm(
    im,
    *,
    axis: int,
    D: float = 1e-09,
    voxel_size: float,
    tol: float = 0.001,
    n_steps: int = 100000,
    n_flux_slices: int = 5,
    sparse: bool = False,
    verbose: bool = True,
) -> TortuosityResult

Compute tortuosity and effective diffusivity using LBM (D3Q7 BGK).

Solves the steady-state diffusion equation on the pore space of a 3D binary image using the Lattice Boltzmann Method.

Parameters:

im
((ndarray, shape(nx, ny, nz))) –

Binary image. True (or 1) = pore, False (or 0) = solid.
axis
(int) –

Axis along which to apply the concentration gradient (0=x, 1=y, 2=z).
D
(float, default: 1e-09 ) –

Bulk diffusivity in m²/s. Default 1e-9.
voxel_size
(float) –

Physical voxel edge length in metres.
tol
(float, default: 0.001 ) –

Convergence tolerance on relative concentration change.
n_steps
(int, default: 100000 ) –

Maximum number of LBM iterations.
n_flux_slices
(int, default: 5 ) –

Number of interior slices to average when computing the flux used to derive D_eff and tau. Default 5. Use 1 for a single-midplane estimate.
sparse
(bool, default: False ) –

Use Taichi sparse storage.
verbose
(bool, default: True ) –

Show a rich progress bar. Default False.

Returns:

result ( TortuosityResult ) –

Source code in src/poromics/_metrics.py

def tortuosity_lbm(
    im,
    *,
    axis: int,
    D: float = 1e-9,
    voxel_size: float,
    tol: float = 1e-3,
    n_steps: int = 100_000,
    n_flux_slices: int = 5,
    sparse: bool = False,
    verbose: bool = True,
) -> TortuosityResult:
    """Compute tortuosity and effective diffusivity using LBM (D3Q7 BGK).

    Solves the steady-state diffusion equation on the pore space of a
    3D binary image using the Lattice Boltzmann Method.

    Parameters
    ----------
    im : ndarray, shape (nx, ny, nz)
        Binary image. True (or 1) = pore, False (or 0) = solid.
    axis : int
        Axis along which to apply the concentration gradient
        (0=x, 1=y, 2=z).
    D : float
        Bulk diffusivity in m²/s. Default 1e-9.
    voxel_size : float
        Physical voxel edge length in metres.
    tol : float
        Convergence tolerance on relative concentration change.
    n_steps : int
        Maximum number of LBM iterations.
    n_flux_slices : int
        Number of interior slices to average when computing the flux
        used to derive ``D_eff`` and ``tau``. Default 5. Use 1 for a
        single-midplane estimate.
    sparse : bool
        Use Taichi sparse storage.
    verbose : bool
        Show a rich progress bar. Default False.

    Returns
    -------
    result : TortuosityResult
    """
    im = np.atleast_3d(np.asarray(im, dtype=bool))
    n_pore_before = int(im.sum())
    im = _trim_nonpercolating(im, axis)
    if im.sum() == 0:
        raise RuntimeError("No percolating paths along the given axis found in the image.")
    n_removed = n_pore_before - int(im.sum())
    if n_removed > 0:
        logger.warning(f"Trimmed {n_removed} non-percolating pore voxels from the image.")
    solver = TransientDiffusion(im, axis=axis, D=D, voxel_size=voxel_size, sparse=sparse)
    solver.run(n_steps=n_steps, tol=tol, verbose=verbose)
    if not solver.converged:
        logger.warning(
            f"LBM diffusion solver did not converge after "
            f"{solver.n_iterations} steps (tol={tol:.2e}). "
            f"Reported tau / D_eff are from a pre-steady-state field; "
            f"increase n_steps or loosen tol."
        )
    c = solver.concentration

    porosity = float(im.sum()) / im.size
    L = im.shape[axis]
    J_mean = solver.flux(axis, n_slices=n_flux_slices)
    D_eff_lu = J_mean * L  # delta_c = 1.0
    D_eff_norm = D_eff_lu / _d3q7.D
    if D_eff_norm > 0:
        formation_factor = 1.0 / D_eff_norm
    else:
        formation_factor = float("inf")
    tau = formation_factor * porosity

    return TortuosityResult(
        im=im,
        axis=axis,
        porosity=porosity,
        tau=tau,
        D_eff=D_eff_norm,
        formation_factor=formation_factor,
        c=c,
        D=D,
        converged=solver.converged,
        n_iterations=solver.n_iterations,
    )

Simulation solvers (poromics.simulation)

Lower-level solver classes that accept physical SI units and expose intermediate state for custom post-processing.

simulation

Classes:

TransientDiffusion –

Transient diffusion solver on 3D voxel images.
TransientFlow –

Transient single-phase flow solver on 3D voxel images.

TransientDiffusion

TransientDiffusion(
    im,
    axis,
    D,
    voxel_size,
    c_in=1.0,
    c_out=0.0,
    sparse=False,
)

Transient diffusion solver on 3D voxel images.

Solves the time-dependent diffusion equation on the pore space using the Lattice Boltzmann Method (D3Q7 BGK). Boundary conditions are Dirichlet (fixed concentration) on the inlet/outlet faces and periodic on the remaining faces.

Parameters:

im
((ndarray, shape(nx, ny, nz))) –

Binary image. True (or 1) = pore, False (or 0) = solid.
axis
(int) –

Axis along which to apply the concentration gradient (0=x, 1=y, 2=z).
D
(float) –

Bulk diffusivity in m²/s.
voxel_size
(float) –

Physical voxel edge length in metres.
c_in
(float, default: 1.0 ) –

Inlet concentration. Default 1.0.
c_out
(float, default: 0.0 ) –

Outlet concentration. Default 0.0.
sparse
(bool, default: False ) –

Use Taichi sparse storage. Default False.

Examples:

>>> solver = TransientDiffusion(im, axis=0, D=1e-9, voxel_size=1e-6)
>>> solver.run(tol=1e-2)
>>> c = solver.concentration

Methods:

flux –

Mean diffusive flux averaged across interior slices (lattice units).
run –

Run the solver to steady state.
step –

Advance by one time step.

Attributes:

concentration –

Concentration field, shape (nx, ny, nz).
converged –

Whether the solver has converged (set by run()).
dt –

Physical time step in seconds.
n_iterations –

Total number of time steps taken.
voxel_size –

Voxel edge length in metres.

Source code in src/poromics/simulation/_diffusion.py

def __init__(self, im, axis, D, voxel_size, c_in=1.0, c_out=0.0,sparse=False):  # fmt: skip
    if axis not in (0, 1, 2):
        raise ValueError(f"axis must be 0, 1, or 2, got {axis}")
    if D <= 0:
        raise ValueError(f"D must be > 0, got {D}")
    if voxel_size <= 0:
        raise ValueError(f"voxel_size must be > 0, got {voxel_size}")
    im_arr = np.atleast_3d(np.asarray(im))
    solid = (im_arr == 0).astype(np.int8)
    if solid.sum() == solid.size:
        raise RuntimeError("Image has no pore voxels.")

    self._axis = axis
    self._D = D
    self._voxel_size = float(voxel_size)
    self._dt = _d3q7.D * voxel_size**2 / D
    self._n_iterations = 0
    self._converged = False

    from .._lbm._taichi_helpers import ensure_taichi
    from .._lbm._diffusion_solver import _D3Q7Solver

    ti = ensure_taichi()
    self._solver = _D3Q7Solver(ti, solid, D=_d3q7.D, sparse=sparse)
    inlet, outlet = axis_to_face[axis]
    self._solver.set_bc(inlet, c_in)
    self._solver.set_bc(outlet, c_out)
    self._solver.init_fields()

concentration `property`

concentration

Concentration field, shape (nx, ny, nz).

Units match c_in and c_out (default: dimensionless).

converged `property`

converged

Whether the solver has converged (set by run()).

dt `property`

dt

Physical time step in seconds.

n_iterations `property`

n_iterations

Total number of time steps taken.

voxel_size `property`

voxel_size

Voxel edge length in metres.

flux

flux(axis=None, n_slices=5)

Mean diffusive flux averaged across interior slices (lattice units).

Uses Fick's law on the concentration gradient rather than distribution moments (which are unreliable at Dirichlet faces). Averages -D_lu * dc/dx_lu across n_slices slices evenly spaced in the central 60% of the domain, which smooths out local noise while avoiding boundary-layer artifacts. Equals the normalized effective diffusivity D_eff/D_0 for a unit concentration drop across the domain.

Parameters:

axis
(int or None, default: None ) –

Axis along which to compute flux. Defaults to the axis used in the constructor.
n_slices
(int, default: 5 ) –

Number of interior slices to average. Default 5.

Returns:

J_mean ( float ) –

Source code in src/poromics/simulation/_diffusion.py

def flux(self, axis=None, n_slices=5):
    """Mean diffusive flux averaged across interior slices (lattice units).

    Uses Fick's law on the concentration gradient rather than
    distribution moments (which are unreliable at Dirichlet faces).
    Averages ``-D_lu * dc/dx_lu`` across ``n_slices`` slices evenly
    spaced in the central 60% of the domain, which smooths out
    local noise while avoiding boundary-layer artifacts. Equals
    the normalized effective diffusivity D_eff/D_0 for a unit
    concentration drop across the domain.

    Parameters
    ----------
    axis : int or None
        Axis along which to compute flux. Defaults to the
        axis used in the constructor.
    n_slices : int
        Number of interior slices to average. Default 5.

    Returns
    -------
    J_mean : float
    """
    if axis is None:
        axis = self._axis
    return self._solver.compute_flux(axis, n_slices=n_slices)

run

run(
    n_steps=100000, tol=0.001, log_every=500, verbose=False
)

Run the solver to steady state.

Parameters:

n_steps
(int, default: 100000 ) –

Maximum number of iterations.
tol
(float or None, default: 0.001 ) –

Convergence tolerance on relative concentration change. None disables early stopping.
log_every
(int, default: 500 ) –

Log convergence every this many steps.
verbose
(bool, default: False ) –

Show a progress bar. Default False.

Source code in src/poromics/simulation/_diffusion.py

def run(self, n_steps=100_000, tol=1e-3, log_every=500,verbose=False):  # fmt: skip
    """Run the solver to steady state.

    Parameters
    ----------
    n_steps : int
        Maximum number of iterations.
    tol : float or None
        Convergence tolerance on relative concentration change.
        None disables early stopping.
    log_every : int
        Log convergence every this many steps.
    verbose : bool
        Show a progress bar. Default False.
    """
    self._converged = False
    t_start = time.time()
    c_prev = self._solver.get_concentration()
    pbar = None
    if verbose:
        pbar = make_progress(n_steps, tol, "Diffusion")
    try:
        for step in range(n_steps + 1):
            self._solver.step()
            self._n_iterations += 1
            if step % log_every != 0:
                continue
            elapsed = time.time() - t_start
            c_now = self._solver.get_concentration()
            c_total = np.sum(np.abs(c_now))
            c_change = np.sum(np.abs(c_now - c_prev))
            ratio = c_change / c_total if c_total > 0 else 0.0
            logger.info(
                f"Step {step:>6d}/{n_steps}  "
                f"|c|={c_total:.4e}  "
                f"delta={ratio:.2e}  "
                f"elapsed={elapsed:.1f}s"
            )
            if pbar is not None:
                update_progress(pbar, step, ratio, tol, n_steps)
            if tol is not None and step > 0 and c_total > 0 and ratio < tol:
                logger.info(
                    f"Converged at step {step} (delta|c|/|c|={ratio:.2e} < tol={tol:.2e})"
                )
                if pbar is not None:
                    pbar.n = 100
                    pbar.set_postfix_str("converged")
                    pbar.refresh()
                self._converged = True
                return
            c_prev = c_now
    finally:
        if pbar is not None:
            pbar.close()

step

step()

Advance by one time step.

Source code in src/poromics/simulation/_diffusion.py

def step(self):
    """Advance by one time step."""
    self._solver.step()
    self._n_iterations += 1

TransientFlow

TransientFlow(
    im,
    axis,
    nu,
    voxel_size,
    rho=None,
    rho_in=1.0,
    rho_out=0.99,
    sparse=False,
)

Transient single-phase flow solver on 3D voxel images.

Solves the incompressible Navier-Stokes equations (Stokes regime) on the pore space using the Lattice Boltzmann Method (D3Q19 MRT). Boundary conditions are fixed-pressure on the inlet/outlet faces and periodic on the remaining faces.

Parameters:

im
((ndarray, shape(nx, ny, nz))) –

Binary image. True (or 1) = pore, False (or 0) = solid.
axis
(int) –

Axis along which to apply the pressure gradient (0=x, 1=y, 2=z).
nu
(float) –

Kinematic viscosity in m²/s.
voxel_size
(float) –

Physical voxel edge length in metres.
rho
(float or None, default: None ) –

Fluid density in kg/m³. Required to expose pressure in pascals; if None, only kinematic_pressure (m²/s²) is available. Permeability and velocity are independent of density (Stokes regime), so None is fine for those.
rho_in
(float, default: 1.0 ) –

Inlet lattice density (pressure BC). Default 1.0.
rho_out
(float, default: 0.99 ) –

Outlet lattice density (pressure BC). Default 0.99.
sparse
(bool, default: False ) –

Use Taichi sparse storage. Default False.

Examples:

>>> solver = TransientFlow(im, axis=0, nu=1e-6, voxel_size=1e-6,
...                        rho=1000.0)
>>> solver.run(tol=1e-3)
>>> v = solver.velocity            # m/s
>>> P = solver.pressure            # Pa (gauge, needs rho)
>>> Pk = solver.kinematic_pressure # m²/s² (density-free)

Methods:

run –

Run the solver to steady state.
step –

Advance by one time step.

Attributes:

converged –

Whether the solver has converged (set by run()).
dt –

Physical time step in seconds.
kinematic_pressure –

Kinematic gauge pressure field in m²/s², shape (nx, ny, nz).
n_iterations –

Total number of time steps taken.
pressure –

Gauge pressure field in Pa, shape (nx, ny, nz).
velocity –

Velocity field in m/s, shape (nx, ny, nz, 3).
voxel_size –

Voxel edge length in metres.

Source code in src/poromics/simulation/_flow.py

def __init__(self, im, axis, nu, voxel_size, rho=None, rho_in=1.0,rho_out=0.99, sparse=False):  # fmt: skip
    if axis not in (0, 1, 2):
        raise ValueError(f"axis must be 0, 1, or 2, got {axis}")
    if nu <= 0:
        raise ValueError(f"nu must be > 0, got {nu}")
    if voxel_size <= 0:
        raise ValueError(f"voxel_size must be > 0, got {voxel_size}")
    if rho is not None and rho <= 0:
        raise ValueError(f"rho must be > 0 or None, got {rho}")
    im_arr = np.atleast_3d(np.asarray(im))
    solid = (im_arr == 0).astype(np.int8)
    if solid.sum() == solid.size:
        raise RuntimeError("Image has no pore voxels.")

    self._axis = axis
    self._nu = nu
    self._voxel_size = float(voxel_size)
    self._rho = rho
    self._rho_in = rho_in
    self._rho_out = rho_out
    self._dt = _d3q19.nu * voxel_size**2 / nu
    self._n_iterations = 0
    self._converged = False

    from .._lbm._taichi_helpers import ensure_taichi
    from .._lbm._flow_solver import _D3Q19Solver

    ti = ensure_taichi()
    self._solver = _D3Q19Solver(ti, solid, nu=_d3q19.nu, sparse=sparse)
    inlet, outlet = axis_to_face[axis]
    self._solver.set_bc_rho(inlet, rho_in)
    self._solver.set_bc_rho(outlet, rho_out)
    self._solver.init_fields()

converged `property`

converged

Whether the solver has converged (set by run()).

dt `property`

dt

Physical time step in seconds.

kinematic_pressure `property`

kinematic_pressure

Kinematic gauge pressure field in m²/s², shape (nx, ny, nz).

Equals p / rho, i.e. the density-free part of the LBM equation of state p = rho * cs². Always available, regardless of whether a fluid density was provided.

n_iterations `property`

n_iterations

Total number of time steps taken.

pressure `property`

pressure

Gauge pressure field in Pa, shape (nx, ny, nz).

Relative to the outlet face. Requires the fluid density rho to have been provided in the constructor; otherwise use kinematic_pressure (m²/s²) instead.

velocity `property`

velocity

Velocity field in m/s, shape (nx, ny, nz, 3).

voxel_size `property`

voxel_size

Voxel edge length in metres.

run

run(
    n_steps=100000, tol=0.001, log_every=500, verbose=False
)

Run the solver to steady state.

Parameters:

n_steps
(int, default: 100000 ) –

Maximum number of iterations.
tol
(float or None, default: 0.001 ) –

Convergence tolerance on relative velocity change. None disables early stopping.
log_every
(int, default: 500 ) –

Log convergence every this many steps.
verbose
(bool, default: False ) –

Show a progress bar. Default False.

Source code in src/poromics/simulation/_flow.py

def run(self, n_steps=100_000, tol=1e-3, log_every=500,verbose=False):  # fmt: skip
    """Run the solver to steady state.

    Parameters
    ----------
    n_steps : int
        Maximum number of iterations.
    tol : float or None
        Convergence tolerance on relative velocity change.
        None disables early stopping.
    log_every : int
        Log convergence every this many steps.
    verbose : bool
        Show a progress bar. Default False.
    """
    self._converged = False
    t_start = time.time()
    v_prev = self._solver.get_velocity()
    pbar = None
    if verbose:
        pbar = make_progress(n_steps, tol, "Flow")
    try:
        for step in range(n_steps + 1):
            self._solver.step()
            self._n_iterations += 1
            if step % log_every != 0:
                continue
            elapsed = time.time() - t_start
            v_now = self._solver.get_velocity()
            v_total = np.sum(np.abs(v_now))
            v_change = np.sum(np.abs(v_now - v_prev))
            ratio = v_change / v_total if v_total > 0 else 0.0
            logger.info(
                f"Step {step:>6d}/{n_steps}  "
                f"|v|={v_total:.4e}  "
                f"delta={ratio:.2e}  "
                f"elapsed={elapsed:.1f}s"
            )
            if pbar is not None:
                update_progress(pbar, step, ratio, tol, n_steps)
            if tol is not None and step > 0 and v_total > 0 and ratio < tol:
                logger.info(
                    f"Converged at step {step} (delta|v|/|v|={ratio:.2e} < tol={tol:.2e})"
                )
                if pbar is not None:
                    pbar.n = 100
                    pbar.set_postfix_str("converged")
                    pbar.refresh()
                self._converged = True
                return
            v_prev = v_now
    finally:
        if pbar is not None:
            pbar.close()

step

step()

Advance by one time step.

Source code in src/poromics/simulation/_flow.py

def step(self):
    """Advance by one time step."""
    self._solver.step()
    self._n_iterations += 1

Julia helpers (poromics.julia_helpers)

Helper functions for interacting with Julia and the Tortuosity.jl package. These are used internally by tortuosity_fd. Feel free to explore them if you're interested in the implementation details.

julia_helpers

Functions:

ensure_gpu_backend –

Install and load a GPU backend in Julia; verify the device is functional.
ensure_julia_deps_ready –

Ensures Julia and Tortuosity.jl are installed.
import_backend –

Imports Tortuosity.jl package from Julia.
import_package –

Imports a package in Julia and returns the module.
init_julia –

Initializes Julia and returns the Main module.
install_backend –

Installs Julia dependencies for Poromics.
install_julia –

Installs Julia using juliapkg.
is_backend_installed –

Checks if Tortuosity.jl is installed.
is_julia_installed –

Checks that Julia is installed.
remove_julia_env –

Removes the active Julia environment directory.

ensure_gpu_backend

ensure_gpu_backend(Main: Any) -> str | None

Install and load a GPU backend in Julia; verify the device is functional.

Returns the backend package name (e.g. "Metal") when a usable GPU is available, or None when the caller should fall back to CPU. All fallback paths emit a warning explaining why.

Parameters:

Main
(juliacall Main module) –

Returns:

backend ( str or None ) –

"Metal", "CUDA", "AMDGPU", or None.

Source code in src/poromics/julia_helpers.py

def ensure_gpu_backend(Main: Any) -> str | None:
    """Install and load a GPU backend in Julia; verify the device is functional.

    Returns the backend package name (e.g. ``"Metal"``) when a usable GPU is
    available, or ``None`` when the caller should fall back to CPU. All
    fallback paths emit a warning explaining why.

    Parameters
    ----------
    Main : juliacall Main module

    Returns
    -------
    backend : str or None
        ``"Metal"``, ``"CUDA"``, ``"AMDGPU"``, or ``None``.
    """
    from juliacall import JuliaError

    backend = _detect_gpu_backend()
    if backend is None:
        logger.warning(
            f"No GPU backend mapping for {platform.system()}/{platform.machine()}. "
            "Set POROMICS_GPU_BACKEND={metal,cuda,amdgpu} to override. "
            "Falling back to CPU."
        )
        return None
    try:
        Main.seval(f"using {backend}")
    except JuliaError:
        try:
            Main.seval(f'import Pkg; Pkg.add("{backend}"); using {backend}')
        except JuliaError as e:
            logger.warning(f"Failed to install/load {backend}.jl: {e}. Falling back to CPU.")
            return None
    # Upstream Tortuosity{Backend}Ext.__init__ registers the backend only when
    # `.functional()` is true. When the package loads but the hardware is not
    # usable, Tortuosity errors later with a "no GPU backend is registered"
    # message that doesn't mention the real cause. See
    # https://github.com/ma-sadeghi/Tortuosity.jl/issues/79 — we pre-check here
    # so the user gets a meaningful warning before falling back to CPU.
    try:
        functional = bool(Main.seval(f"{backend}.functional()"))
    except JuliaError:
        functional = False
    if not functional:
        logger.warning(
            f"{backend}.jl loaded but {backend}.functional() returned false "
            "(no working GPU device detected). Falling back to CPU."
        )
        return None
    return backend

ensure_julia_deps_ready

ensure_julia_deps_ready(
    quiet: bool = False, retry: bool = True
) -> None

Ensures Julia and Tortuosity.jl are installed.

Args: quiet: If True, suppresses output during installation. Default is False. retry: If True, retries the installation if it fails. Default is True.

Raises: ImportError: If Julia or Tortuosity.jl cannot be installed.

Source code in src/poromics/julia_helpers.py

def ensure_julia_deps_ready(quiet: bool = False, retry: bool = True) -> None:
    """Ensures Julia and Tortuosity.jl are installed.

    Args:
        quiet: If True, suppresses output during installation. Default is False.
        retry: If True, retries the installation if it fails. Default is True.

    Raises:
        ImportError: If Julia or Tortuosity.jl cannot be installed.
    """

    def _ensure_julia_deps_ready(quiet):
        if not is_julia_installed(error=False):
            logger.warning("Julia not found, installing Julia...")
            install_julia(quiet=quiet)
        Main = init_julia(quiet=quiet)
        if not is_backend_installed(Main=Main, error=False):
            logger.warning("Julia dependencies not found, installing Tortuosity.jl...")
            install_backend(quiet=quiet)

    def _reset_julia_env(quiet):
        remove_julia_env()
        if quiet:
            with suppress_output():
                juliapkg.resolve(force=True)
        else:
            juliapkg.resolve(force=True)

    try:
        _ensure_julia_deps_ready(quiet)
    except Exception:
        if retry:
            _reset_julia_env(quiet)
            _ensure_julia_deps_ready(quiet)
            return
        raise

import_backend

import_backend(Main: Any = None) -> Any

Imports Tortuosity.jl package from Julia.

Args: Main: Julia Main module. Default is None. If None, the Main module will be initialized.

Returns: backend: Handle to the Tortuosity.jl package.

Raises: ImportError: If Julia is not installed or the package is not found.

Source code in src/poromics/julia_helpers.py

def import_backend(Main: Any = None) -> Any:
    """Imports Tortuosity.jl package from Julia.

    Args:
        Main: Julia Main module. Default is None. If None, the Main module will
            be initialized.

    Returns:
        backend: Handle to the Tortuosity.jl package.

    Raises:
        ImportError: If Julia is not installed or the package is not found.
    """
    Main = init_julia() if Main is None else Main
    is_backend_installed(Main=Main, error=True)
    return import_package("Tortuosity", Main)

import_package

import_package(
    package_name: str, Main: Any, error: bool = False
) -> Any

Imports a package in Julia and returns the module.

Args: package_name: Name of the Julia package to import. Main: Julia Main module. error: If True, raises an error if the package is not found. Default is False.

Returns: package: Handle to the imported package.

Raises: ImportError: If the package is not found and error is True.

Source code in src/poromics/julia_helpers.py

def import_package(package_name: str, Main: Any, error: bool = False) -> Any:
    """Imports a package in Julia and returns the module.

    Args:
        package_name: Name of the Julia package to import.
        Main: Julia Main module.
        error: If True, raises an error if the package is not found. Default is False.

    Returns:
        package: Handle to the imported package.

    Raises:
        ImportError: If the package is not found and error is True.
    """
    from juliacall import JuliaError

    try:
        Main.seval(f"using {package_name}")
        return eval(f"Main.{package_name}")
    except JuliaError as e:
        if error:
            raise e
    return None

init_julia

init_julia(quiet: bool = False) -> Any

Initializes Julia and returns the Main module.

Args: quiet: If True, suppresses the output of Julia initialization. Default is False.

Returns: Main: The Julia Main module.

Raises: ImportError: If Julia is not installed.

Source code in src/poromics/julia_helpers.py

def init_julia(quiet: bool = False) -> Any:
    """Initializes Julia and returns the Main module.

    Args:
        quiet: If True, suppresses the output of Julia initialization. Default is False.

    Returns:
        Main: The Julia Main module.

    Raises:
        ImportError: If Julia is not installed.
    """
    is_julia_installed(error=True)
    if not can_skip_resolve():
        logger.warning("Julia is installed, but needs to be resolved...")
    if quiet:
        with suppress_output():
            from juliacall import Main  # type: ignore
    else:
        from juliacall import Main  # type: ignore

    return Main

install_backend

install_backend(quiet: bool = False) -> None

Installs Julia dependencies for Poromics.

Args: quiet: If True, suppresses output during installation. Default is False.

Raises: ImportError: If Julia is not installed.

Source code in src/poromics/julia_helpers.py

def install_backend(quiet: bool = False) -> None:
    """Installs Julia dependencies for Poromics.

    Args:
        quiet: If True, suppresses output during installation. Default is False.

    Raises:
        ImportError: If Julia is not installed.
    """
    is_julia_installed(error=True)

    if quiet:
        with suppress_output():
            juliapkg.resolve()
    else:
        juliapkg.resolve()

install_julia

install_julia(quiet: bool = False) -> None

Installs Julia using juliapkg.

Args: quiet: If True, suppresses output during installation. Default is False.

Source code in src/poromics/julia_helpers.py

def install_julia(quiet: bool = False) -> None:
    """Installs Julia using juliapkg.

    Args:
        quiet: If True, suppresses output during installation. Default is False.
    """
    # Importing juliacall automatically installs Julia using juliapkg
    if quiet:
        with suppress_output():
            import juliacall  # noqa: F401
    else:
        import juliacall  # noqa: F401

is_backend_installed

is_backend_installed(
    Main: Any = None, error: bool = False
) -> bool

Checks if Tortuosity.jl is installed.

Args: Main: Julia Main module. Default is None. If None, it will be initialized. error: If True, raises an error if backend is not found. Default is False.

Returns: flag: True if the package is installed, False otherwise.

Raises: ImportError: If Julia is not installed or backend is not found and error is True.

Source code in src/poromics/julia_helpers.py

def is_backend_installed(Main: Any = None, error: bool = False) -> bool:
    """Checks if Tortuosity.jl is installed.

    Args:
        Main: Julia Main module. Default is None. If None, it will be initialized.
        error: If True, raises an error if backend is not found. Default is False.

    Returns:
        flag: True if the package is installed, False otherwise.

    Raises:
        ImportError: If Julia is not installed or backend is not found and error is True.
    """
    Main = init_julia() if Main is None else Main
    if import_package("Tortuosity", Main, error=False) is not None:
        return True
    msg = "Tortuosity.jl not found, run 'python -m poromics install'"
    if error:
        raise ImportError(msg)
    return False

is_julia_installed

is_julia_installed(error: bool = False) -> bool

Checks that Julia is installed.

Args: error: If True, raises an error if Julia is not found. Default is False.

Returns: flag: True if Julia is installed, False otherwise.

Raises: ImportError: If Julia is not installed and error is True.

Source code in src/poromics/julia_helpers.py

def is_julia_installed(error: bool = False) -> bool:
    """Checks that Julia is installed.

    Args:
        error: If True, raises an error if Julia is not found. Default is False.

    Returns:
        flag: True if Julia is installed, False otherwise.

    Raises:
        ImportError: If Julia is not installed and error is True.
    """
    # Look for system-wide Julia executable
    try:
        find_julia()
        return True
    except Exception:
        pass
    # Look for local Julia executable (e.g., installed by juliapkg)
    if can_skip_resolve():
        return True
    msg = "Julia not found. Visit https://github.com/JuliaLang/juliaup and install Julia."
    if error:
        raise ImportError(msg)
    return False

remove_julia_env

remove_julia_env() -> None

Removes the active Julia environment directory.

When Julia or its dependencies are corrupted, this is a possible fix.

Source code in src/poromics/julia_helpers.py

def remove_julia_env() -> None:
    """Removes the active Julia environment directory.

    When Julia or its dependencies are corrupted, this is a possible fix.
    """
    path_julia_env = Path(juliapkg.project())

    if path_julia_env.exists():
        logger.warning(f"Removing Julia environment directory: {path_julia_env}")
        shutil.rmtree(path_julia_env)
    else:
        logger.warning("Julia environment directory not found.")

API Reference

poromics

PermeabilityResult

im

axis

porosity

k

u_darcy

u_pore

velocity

kinematic_pressure

pressure

converged

n_iterations

plot_pressure

z

ax

kinematic

plot_velocity

z

ax

rescale

voxel_size

nu

rho

SimulationResult

im

axis

porosity

TortuosityResult

im

axis

porosity

tau

D_eff

c

formation_factor

D

converged

n_iterations

plot_concentration

z

ax

permeability_lbm

im

axis

nu

rho

voxel_size

tol

n_steps

sparse

verbose

tortuosity_fd

im

axis

D

rtol

gpu

verbose

tortuosity_lbm

im

axis

D

voxel_size

tol

n_steps

n_flux_slices

sparse

verbose

simulation

TransientDiffusion

im

axis

D

voxel_size

c_in

c_out

sparse

concentration property

`im`

`axis`

`porosity`

`k`

`u_darcy`

`u_pore`

`velocity`

`kinematic_pressure`

`pressure`

`converged`

`n_iterations`

`z`

`ax`

`kinematic`

`z`

`ax`

`voxel_size`

`nu`

`rho`

`im`

`axis`

`porosity`

`im`

`axis`

`porosity`

`tau`

`D_eff`

`c`

`formation_factor`

`D`

`converged`

`n_iterations`

`z`

`ax`

`im`

`axis`

`nu`

`rho`

`voxel_size`

`tol`

`n_steps`

`sparse`

`verbose`

`im`

`axis`

`D`

`rtol`

`gpu`

`verbose`

`im`

`axis`

`D`

`voxel_size`

`tol`

`n_steps`

`n_flux_slices`

`sparse`

`verbose`

`im`

`axis`

`D`

`voxel_size`

`c_in`

`c_out`

`sparse`

concentration `property`

converged `property`

dt `property`

n_iterations `property`

voxel_size `property`

`axis`

`n_slices`

`n_steps`

`tol`

`log_every`

`verbose`

`im`

`axis`

`nu`

`voxel_size`