Skip to content

Geometry Helpers

The voids.geom sub-package provides functions for normalizing characteristic pore and throat sizes and computing hydraulic conductances.


Characteristic Size

voids.geom.characteristic

area_equivalent_diameter

area_equivalent_diameter(area)

Return the circular-equivalent diameter associated with an area array.

Parameters:

Name Type Description Default
area ndarray

Cross-sectional area array.

required

Returns:

Type Description
ndarray

Diameter defined by d = 2 * sqrt(area / pi).

Source code in src/voids/geom/characteristic.py
def area_equivalent_diameter(area: np.ndarray) -> np.ndarray:
    """Return the circular-equivalent diameter associated with an area array.

    Parameters
    ----------
    area :
        Cross-sectional area array.

    Returns
    -------
    numpy.ndarray
        Diameter defined by ``d = 2 * sqrt(area / pi)``.
    """

    area = np.asarray(area, dtype=float)
    return np.asarray(2.0 * np.sqrt(area / np.pi))

normalize_characteristic_size

normalize_characteristic_size(values, *, field_name)

Normalize a size-like field to a characteristic diameter surrogate.

Parameters:

Name Type Description Default
values ndarray

Raw size-like field values.

required
field_name str | None

Source field name, used to convert radii and areas to diameters.

required

Returns:

Type Description
ndarray

Characteristic-size array interpreted as a diameter-like quantity.

Source code in src/voids/geom/characteristic.py
def normalize_characteristic_size(
    values: np.ndarray,
    *,
    field_name: str | None,
) -> np.ndarray:
    """Normalize a size-like field to a characteristic diameter surrogate.

    Parameters
    ----------
    values :
        Raw size-like field values.
    field_name :
        Source field name, used to convert radii and areas to diameters.

    Returns
    -------
    numpy.ndarray
        Characteristic-size array interpreted as a diameter-like quantity.
    """

    arr = np.asarray(values, dtype=float)
    if field_name == "radius_inscribed":
        return 2.0 * arr
    if field_name == "area":
        return area_equivalent_diameter(arr)
    return arr

characteristic_size

characteristic_size(
    store,
    *,
    expected_shape=None,
    fields=_CHARACTERISTIC_SIZE_FIELDS,
)

Return a preferred characteristic size array from a pore/throat store.

Parameters:

Name Type Description Default
store Mapping[str, object]

Mapping such as net.pore or net.throat.

required
expected_shape tuple[int, ...] | None

Optional expected array shape.

None
fields tuple[str, ...]

Ordered field priority. The default is diameter_equivalent -> diameter_inscribed -> radius_inscribed -> area.

_CHARACTERISTIC_SIZE_FIELDS

Returns:

Type Description
tuple

Pair (values, label) where values is a characteristic-size array and label is the originating field name.

Raises:

Type Description
KeyError

If none of the requested fields exists.

ValueError

If a selected field does not match expected_shape.

Source code in src/voids/geom/characteristic.py
def characteristic_size(
    store: Mapping[str, object],
    *,
    expected_shape: tuple[int, ...] | None = None,
    fields: tuple[str, ...] = _CHARACTERISTIC_SIZE_FIELDS,
) -> tuple[np.ndarray, str]:
    """Return a preferred characteristic size array from a pore/throat store.

    Parameters
    ----------
    store :
        Mapping such as ``net.pore`` or ``net.throat``.
    expected_shape :
        Optional expected array shape.
    fields :
        Ordered field priority. The default is
        ``diameter_equivalent -> diameter_inscribed -> radius_inscribed -> area``.

    Returns
    -------
    tuple
        Pair ``(values, label)`` where ``values`` is a characteristic-size array
        and ``label`` is the originating field name.

    Raises
    ------
    KeyError
        If none of the requested fields exists.
    ValueError
        If a selected field does not match ``expected_shape``.
    """

    for field_name in fields:
        if field_name not in store:
            continue
        arr = normalize_characteristic_size(
            np.asarray(store[field_name], dtype=float), field_name=field_name
        )
        if expected_shape is not None and arr.shape != expected_shape:
            raise ValueError(f"field '{field_name}' must have shape {expected_shape}")
        return arr, field_name
    raise KeyError(f"Need one of the characteristic size fields: {fields}")

Hydraulic Geometry

voids.geom.hydraulic

generic_poiseuille_conductance

generic_poiseuille_conductance(
    net, viscosity, *, throat_viscosity=None
)

Compute throat conductance using a circular Poiseuille approximation.

Parameters:

Name Type Description Default
net Network

Network containing throat geometry or precomputed hydraulic conductance.

required
viscosity float | ndarray | None

Dynamic viscosity.

required

Returns:

Type Description
ndarray

Conductance array for all throats.

Raises:

Type Description
ValueError

If viscosity is non-positive or if precomputed conductance is negative.

KeyError

If required geometry is missing.

Notes

When no precomputed conductance is supplied, the model uses

g = pi * r**4 / (8 * mu * L)

with radius inferred from either throat.diameter_inscribed or throat.area.

Source code in src/voids/geom/hydraulic.py
def generic_poiseuille_conductance(
    net: Network,
    viscosity: float | np.ndarray | None,
    *,
    throat_viscosity: float | np.ndarray | None = None,
) -> np.ndarray:
    """Compute throat conductance using a circular Poiseuille approximation.

    Parameters
    ----------
    net :
        Network containing throat geometry or precomputed hydraulic conductance.
    viscosity :
        Dynamic viscosity.

    Returns
    -------
    numpy.ndarray
        Conductance array for all throats.

    Raises
    ------
    ValueError
        If viscosity is non-positive or if precomputed conductance is negative.
    KeyError
        If required geometry is missing.

    Notes
    -----
    When no precomputed conductance is supplied, the model uses

    ``g = pi * r**4 / (8 * mu * L)``

    with radius inferred from either ``throat.diameter_inscribed`` or
    ``throat.area``.
    """

    selected_viscosity = throat_viscosity if throat_viscosity is not None else viscosity
    if selected_viscosity is None:
        raise ValueError("Need either viscosity or throat_viscosity")
    mu_t = _broadcast_viscosity(selected_viscosity, (net.Nt,))
    if "hydraulic_conductance" in net.throat:
        g = np.asarray(net.throat["hydraulic_conductance"], dtype=float)
        if (g < 0).any():
            raise ValueError("throat.hydraulic_conductance contains negative values")
        return g.copy()

    _require(net, "throat", ("length",))
    L = np.asarray(net.throat["length"], dtype=float)
    if "diameter_inscribed" in net.throat:
        d = np.asarray(net.throat["diameter_inscribed"], dtype=float)
    elif "area" in net.throat:
        d = _diameter_from_area(np.asarray(net.throat["area"], dtype=float))
    else:
        raise KeyError(
            "Need throat.diameter_inscribed or throat.area (or precomputed hydraulic_conductance)"
        )
    r = 0.5 * d
    return np.asarray((np.pi * r**4) / (8.0 * mu_t * L), dtype=float)

valvatne_blunt_throat_conductance

valvatne_blunt_throat_conductance(
    net, viscosity, *, throat_viscosity=None
)

Compute shape-aware throat conductance using throat geometry only.

Parameters:

Name Type Description Default
net Network

Network containing throat length and cross-sectional geometry.

required
viscosity float | ndarray | None

Dynamic viscosity.

required

Returns:

Type Description
ndarray

Shape-aware throat conductance array.

Raises:

Type Description
ValueError

If viscosity is not positive.

KeyError

If required throat geometry is unavailable.

Source code in src/voids/geom/hydraulic.py
def valvatne_blunt_throat_conductance(
    net: Network,
    viscosity: float | np.ndarray | None,
    *,
    throat_viscosity: float | np.ndarray | None = None,
) -> np.ndarray:
    """Compute shape-aware throat conductance using throat geometry only.

    Parameters
    ----------
    net :
        Network containing throat length and cross-sectional geometry.
    viscosity :
        Dynamic viscosity.

    Returns
    -------
    numpy.ndarray
        Shape-aware throat conductance array.

    Raises
    ------
    ValueError
        If viscosity is not positive.
    KeyError
        If required throat geometry is unavailable.
    """

    selected_viscosity = throat_viscosity if throat_viscosity is not None else viscosity
    if selected_viscosity is None:
        raise ValueError("Need either viscosity or throat_viscosity")
    _broadcast_viscosity(selected_viscosity, (net.Nt,))
    if "hydraulic_conductance" in net.throat:
        return generic_poiseuille_conductance(net, viscosity, throat_viscosity=throat_viscosity)
    return _throat_only_shape_factor_conductance(
        net,
        viscosity,
        throat_viscosity=throat_viscosity,
    )

valvatne_blunt_conductance

valvatne_blunt_conductance(
    net,
    viscosity,
    *,
    pore_viscosity=None,
    throat_viscosity=None,
)

Compute a shape-factor-aware single-phase conductance following Valvatne-Blunt.

Parameters:

Name Type Description Default
net Network

Network containing throat and, ideally, pore geometry.

required
viscosity float | ndarray | None

Dynamic viscosity.

required

Returns:

Type Description
ndarray

Throat conductance array.

Raises:

Type Description
ValueError

If viscosity is not positive.

Notes

This implements the single-phase geometric closure used in the Imperial College Valvatne-Blunt style network model:

  • segment conductance is evaluated as g = k * G * A**2 / (mu * L)
  • k = 3/5 for triangular ducts
  • k = 0.5623 for square ducts
  • k = 1/2 for circular ducts

The selection logic is:

  1. If throat.hydraulic_conductance is explicitly present, return it.
  2. Else, if conduit lengths and pore/throat shape data are available, compute a harmonic pore1-core-pore2 conductance.
  3. Else, if throat-only shape data are available, use a throat-only model.
  4. Else, warn and fall back to circular Poiseuille conductance.

This is still a single-phase closure; corner films and multiphase occupancy are intentionally out of scope here.

Source code in src/voids/geom/hydraulic.py
def valvatne_blunt_conductance(
    net: Network,
    viscosity: float | np.ndarray | None,
    *,
    pore_viscosity: float | np.ndarray | None = None,
    throat_viscosity: float | np.ndarray | None = None,
) -> np.ndarray:
    """Compute a shape-factor-aware single-phase conductance following Valvatne-Blunt.

    Parameters
    ----------
    net :
        Network containing throat and, ideally, pore geometry.
    viscosity :
        Dynamic viscosity.

    Returns
    -------
    numpy.ndarray
        Throat conductance array.

    Raises
    ------
    ValueError
        If viscosity is not positive.

    Notes
    -----
    This implements the single-phase geometric closure used in the Imperial
    College Valvatne-Blunt style network model:

    - segment conductance is evaluated as ``g = k * G * A**2 / (mu * L)``
    - ``k = 3/5`` for triangular ducts
    - ``k = 0.5623`` for square ducts
    - ``k = 1/2`` for circular ducts

    The selection logic is:

    1. If ``throat.hydraulic_conductance`` is explicitly present, return it.
    2. Else, if conduit lengths and pore/throat shape data are available, compute a
       harmonic pore1-core-pore2 conductance.
    3. Else, if throat-only shape data are available, use a throat-only model.
    4. Else, warn and fall back to circular Poiseuille conductance.

    This is still a single-phase closure; corner films and multiphase occupancy are
    intentionally out of scope here.
    """

    _resolve_pore_throat_viscosities(
        net,
        viscosity,
        pore_viscosity=pore_viscosity,
        throat_viscosity=throat_viscosity,
    )

    if "hydraulic_conductance" in net.throat:
        return generic_poiseuille_conductance(net, viscosity, throat_viscosity=throat_viscosity)

    try:
        return _valvatne_conduit_baseline(
            net,
            viscosity,
            pore_viscosity=pore_viscosity,
            throat_viscosity=throat_viscosity,
        )
    except KeyError:
        pass

    try:
        return _throat_only_shape_factor_conductance(
            net,
            viscosity,
            throat_viscosity=throat_viscosity,
        )
    except KeyError:
        warnings.warn(
            "Insufficient geometry for shape-factor model; falling back to generic_poiseuille",
            RuntimeWarning,
            stacklevel=2,
        )
        return generic_poiseuille_conductance(net, viscosity, throat_viscosity=throat_viscosity)

valvatne_blunt_baseline_conductance

valvatne_blunt_baseline_conductance(
    net,
    viscosity,
    *,
    pore_viscosity=None,
    throat_viscosity=None,
)

Backward-compatible alias for :func:valvatne_blunt_conductance.

Source code in src/voids/geom/hydraulic.py
def valvatne_blunt_baseline_conductance(
    net: Network,
    viscosity: float | np.ndarray | None,
    *,
    pore_viscosity: float | np.ndarray | None = None,
    throat_viscosity: float | np.ndarray | None = None,
) -> np.ndarray:
    """Backward-compatible alias for :func:`valvatne_blunt_conductance`."""

    return valvatne_blunt_conductance(
        net,
        viscosity,
        pore_viscosity=pore_viscosity,
        throat_viscosity=throat_viscosity,
    )

available_conductance_models

available_conductance_models()

Return the names of built-in hydraulic conductance models.

Returns:

Type Description
tuple of str

Available model names.

Source code in src/voids/geom/hydraulic.py
def available_conductance_models() -> tuple[str, ...]:
    """Return the names of built-in hydraulic conductance models.

    Returns
    -------
    tuple of str
        Available model names.
    """

    return (
        "generic_poiseuille",
        "valvatne_blunt_throat",
        "valvatne_blunt",
        "valvatne_blunt_baseline",
    )

throat_conductance

throat_conductance(
    net,
    viscosity,
    model="generic_poiseuille",
    *,
    pore_viscosity=None,
    throat_viscosity=None,
)

Dispatch to a throat hydraulic conductance model.

Parameters:

Name Type Description Default
net Network

Network containing the required geometry.

required
viscosity float | ndarray | None

Dynamic viscosity.

required
model str

Conductance model name.

'generic_poiseuille'

Returns:

Type Description
ndarray

Throat conductance array.

Raises:

Type Description
ValueError

If model is unknown.

Source code in src/voids/geom/hydraulic.py
def throat_conductance(
    net: Network,
    viscosity: float | np.ndarray | None,
    model: str = "generic_poiseuille",
    *,
    pore_viscosity: float | np.ndarray | None = None,
    throat_viscosity: float | np.ndarray | None = None,
) -> np.ndarray:
    """Dispatch to a throat hydraulic conductance model.

    Parameters
    ----------
    net :
        Network containing the required geometry.
    viscosity :
        Dynamic viscosity.
    model :
        Conductance model name.

    Returns
    -------
    numpy.ndarray
        Throat conductance array.

    Raises
    ------
    ValueError
        If ``model`` is unknown.
    """

    if model == "generic_poiseuille":
        return generic_poiseuille_conductance(net, viscosity, throat_viscosity=throat_viscosity)
    if model == "valvatne_blunt_throat":
        return valvatne_blunt_throat_conductance(
            net,
            viscosity,
            throat_viscosity=throat_viscosity,
        )
    if model == "valvatne_blunt":
        return valvatne_blunt_conductance(
            net,
            viscosity,
            pore_viscosity=pore_viscosity,
            throat_viscosity=throat_viscosity,
        )
    if model == "valvatne_blunt_baseline":
        return valvatne_blunt_conductance(
            net,
            viscosity,
            pore_viscosity=pore_viscosity,
            throat_viscosity=throat_viscosity,
        )
    raise ValueError(f"Unknown conductance model '{model}'")

throat_conductance_with_sensitivities

throat_conductance_with_sensitivities(
    net,
    viscosity,
    model="generic_poiseuille",
    *,
    pore_viscosity=None,
    throat_viscosity=None,
    pore_dviscosity_dpressure=None,
    throat_dviscosity_dpressure=None,
)

Return throat conductance and endpoint pressure sensitivities.

Notes

The returned arrays are (g, dg_dpi, dg_dpj) where i and j are the pore indices in net.throat_conns[:, 0] and net.throat_conns[:, 1].

Source code in src/voids/geom/hydraulic.py
def throat_conductance_with_sensitivities(
    net: Network,
    viscosity: float | np.ndarray | None,
    model: str = "generic_poiseuille",
    *,
    pore_viscosity: float | np.ndarray | None = None,
    throat_viscosity: float | np.ndarray | None = None,
    pore_dviscosity_dpressure: float | np.ndarray | None = None,
    throat_dviscosity_dpressure: float | np.ndarray | None = None,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Return throat conductance and endpoint pressure sensitivities.

    Notes
    -----
    The returned arrays are ``(g, dg_dpi, dg_dpj)`` where ``i`` and ``j`` are
    the pore indices in ``net.throat_conns[:, 0]`` and ``net.throat_conns[:, 1]``.
    """

    conns = net.throat_conns
    i_idx = conns[:, 0]
    j_idx = conns[:, 1]

    if "hydraulic_conductance" in net.throat:
        g = generic_poiseuille_conductance(net, viscosity, throat_viscosity=throat_viscosity)
        zeros = np.zeros_like(g, dtype=float)
        return g, zeros, zeros

    if model == "generic_poiseuille":
        selected_viscosity = throat_viscosity if throat_viscosity is not None else viscosity
        if selected_viscosity is None:
            raise ValueError("Need either viscosity or throat_viscosity")
        mu_t = _broadcast_viscosity(selected_viscosity, (net.Nt,))
        dmu_t = (
            _broadcast_finite(
                throat_dviscosity_dpressure,
                (net.Nt,),
                name="throat_dviscosity_dpressure",
            )
            if throat_dviscosity_dpressure is not None
            else np.zeros(net.Nt, dtype=float)
        )
        _require(net, "throat", ("length",))
        L = np.asarray(net.throat["length"], dtype=float)
        if "diameter_inscribed" in net.throat:
            d = np.asarray(net.throat["diameter_inscribed"], dtype=float)
        elif "area" in net.throat:
            d = _diameter_from_area(np.asarray(net.throat["area"], dtype=float))
        else:
            raise KeyError(
                "Need throat.diameter_inscribed or throat.area (or precomputed hydraulic_conductance)"
            )
        r = 0.5 * d
        g = (np.pi * r**4) / (8.0 * mu_t * L)
        dg_dmu = -g / mu_t
        factor = 0.5 * dmu_t
        return g, dg_dmu * factor, dg_dmu * factor

    if model == "valvatne_blunt_throat":
        selected_viscosity = throat_viscosity if throat_viscosity is not None else viscosity
        if selected_viscosity is None:
            raise ValueError("Need either viscosity or throat_viscosity")
        mu_t = _broadcast_viscosity(selected_viscosity, (net.Nt,))
        dmu_t = (
            _broadcast_finite(
                throat_dviscosity_dpressure,
                (net.Nt,),
                name="throat_dviscosity_dpressure",
            )
            if throat_dviscosity_dpressure is not None
            else np.zeros(net.Nt, dtype=float)
        )
        _require(net, "throat", ("length",))
        A = _get_entity_area(net, "throat")
        G = _get_entity_shape_factor(net, "throat", area=A)
        L = np.asarray(net.throat["length"], dtype=float)
        g = _segment_conductance_valvatne_blunt(A, G, L, mu_t)
        dg_dmu = np.zeros_like(g, dtype=float)
        positive = np.isfinite(g) & (g > 0.0)
        dg_dmu[positive] = -g[positive] / mu_t[positive]
        factor = 0.5 * dmu_t
        return g, dg_dmu * factor, dg_dmu * factor

    if model in {"valvatne_blunt", "valvatne_blunt_baseline"}:
        mu_p, mu_t = _resolve_pore_throat_viscosities(
            net,
            viscosity,
            pore_viscosity=pore_viscosity,
            throat_viscosity=throat_viscosity,
        )
        dmu_p = (
            _broadcast_finite(
                pore_dviscosity_dpressure,
                (net.Np,),
                name="pore_dviscosity_dpressure",
            )
            if pore_dviscosity_dpressure is not None
            else np.zeros(net.Np, dtype=float)
        )
        dmu_t = (
            _broadcast_finite(
                throat_dviscosity_dpressure,
                (net.Nt,),
                name="throat_dviscosity_dpressure",
            )
            if throat_dviscosity_dpressure is not None
            else np.zeros(net.Nt, dtype=float)
        )
        try:
            if not _conduit_lengths_available(net):
                raise KeyError("Missing conduit lengths (pore1_length, core_length, pore2_length)")
            At = _get_entity_area(net, "throat")
            Gt = _get_entity_shape_factor(net, "throat", area=At)
            Lt = np.asarray(net.throat["core_length"], dtype=float)
            gt = _segment_conductance_valvatne_blunt(At, Gt, Lt, mu_t)
            dgt_dmu = np.zeros_like(gt, dtype=float)
            positive_t = np.isfinite(gt) & (gt > 0.0)
            dgt_dmu[positive_t] = -gt[positive_t] / mu_t[positive_t]
            dgt_dpi = dgt_dmu * (0.5 * dmu_t)
            dgt_dpj = dgt_dmu * (0.5 * dmu_t)

            Ap = _get_entity_area(net, "pore")
            Gp = _get_entity_shape_factor(net, "pore", area=Ap)
            L1 = np.asarray(net.throat["pore1_length"], dtype=float)
            L2 = np.asarray(net.throat["pore2_length"], dtype=float)
            g1 = _segment_conductance_valvatne_blunt(Ap[i_idx], Gp[i_idx], L1, mu_p[i_idx])
            g2 = _segment_conductance_valvatne_blunt(Ap[j_idx], Gp[j_idx], L2, mu_p[j_idx])
            dg1_dmu = np.zeros_like(g1, dtype=float)
            dg2_dmu = np.zeros_like(g2, dtype=float)
            positive_1 = np.isfinite(g1) & (g1 > 0.0)
            positive_2 = np.isfinite(g2) & (g2 > 0.0)
            dg1_dmu[positive_1] = -g1[positive_1] / mu_p[i_idx][positive_1]
            dg2_dmu[positive_2] = -g2[positive_2] / mu_p[j_idx][positive_2]
            dg1_dpi = dg1_dmu * dmu_p[i_idx]
            dg2_dpj = dg2_dmu * dmu_p[j_idx]

            g = _harmonic_combine_segments(g1, gt, g2)
            dg_dpi = _harmonic_sensitivity(g, (g1, gt, g2), (dg1_dpi, dgt_dpi, np.zeros_like(g2)))
            dg_dpj = _harmonic_sensitivity(g, (g1, gt, g2), (np.zeros_like(g1), dgt_dpj, dg2_dpj))
            return g, dg_dpi, dg_dpj
        except KeyError:
            try:
                A = _get_entity_area(net, "throat")
                G = _get_entity_shape_factor(net, "throat", area=A)
                L = np.asarray(net.throat["length"], dtype=float)
                g = _segment_conductance_valvatne_blunt(A, G, L, mu_t)
                dg_dmu = np.zeros_like(g, dtype=float)
                positive = np.isfinite(g) & (g > 0.0)
                dg_dmu[positive] = -g[positive] / mu_t[positive]
                factor = 0.5 * dmu_t
                return g, dg_dmu * factor, dg_dmu * factor
            except KeyError:
                warnings.warn(
                    "Insufficient geometry for shape-factor model; falling back to generic_poiseuille",
                    RuntimeWarning,
                    stacklevel=2,
                )
                return throat_conductance_with_sensitivities(
                    net,
                    viscosity,
                    model="generic_poiseuille",
                    pore_viscosity=pore_viscosity,
                    throat_viscosity=throat_viscosity,
                    pore_dviscosity_dpressure=pore_dviscosity_dpressure,
                    throat_dviscosity_dpressure=throat_dviscosity_dpressure,
                )

    raise ValueError(f"Unknown conductance model '{model}'")