Linear Algebra

The voids.linalg sub-package handles matrix assembly, linear-system solving, boundary-condition imposition, and solver diagnostics.

---

Assembly

voids.linalg.assemble

assemble_pressure_system

assemble_pressure_system(net, throat_conductance)

Assemble the pore-pressure matrix for steady single-phase flow.

Parameters:

Name	Type	Description	Default
net	Network	Network defining the pore-throat topology.	required
throat_conductance	ndarray	Conductance array with shape (Nt,).	required

Returns:

Type	Description
csr_matrix	Symmetric matrix A with shape (Np, Np).

Raises:

Type	Description
ValueError	If the conductance array has the wrong shape or contains negative entries.

Notes

The assembled matrix is the conductance-weighted graph Laplacian. For a throat with conductance g_t connecting pores i and j, the local contribution is

A[i, i] += g_t A[j, j] += g_t A[i, j] -= g_t A[j, i] -= g_t

Source code in src/voids/linalg/assemble.py

def assemble_pressure_system(net: Network, throat_conductance: np.ndarray) -> sparse.csr_matrix:
    """Assemble the pore-pressure matrix for steady single-phase flow.

    Parameters
    ----------
    net :
        Network defining the pore-throat topology.
    throat_conductance :
        Conductance array with shape ``(Nt,)``.

    Returns
    -------
    scipy.sparse.csr_matrix
        Symmetric matrix ``A`` with shape ``(Np, Np)``.

    Raises
    ------
    ValueError
        If the conductance array has the wrong shape or contains negative
        entries.

    Notes
    -----
    The assembled matrix is the conductance-weighted graph Laplacian. For a
    throat with conductance ``g_t`` connecting pores ``i`` and ``j``, the local
    contribution is

    ``A[i, i] += g_t``
    ``A[j, j] += g_t``
    ``A[i, j] -= g_t``
    ``A[j, i] -= g_t``
    """

    g = np.asarray(throat_conductance, dtype=float)
    if g.shape != (net.Nt,):
        raise ValueError("throat_conductance must have shape (Nt,)")
    if (g < 0).any():
        raise ValueError("throat_conductance must be nonnegative")
    i = net.throat_conns[:, 0]
    j = net.throat_conns[:, 1]
    rows = np.concatenate([i, j])
    cols = np.concatenate([j, i])
    data = np.concatenate([-g, -g])
    diag = np.zeros(net.Np, dtype=float)
    np.add.at(diag, i, g)
    np.add.at(diag, j, g)
    rows = np.concatenate([rows, np.arange(net.Np)])
    cols = np.concatenate([cols, np.arange(net.Np)])
    data = np.concatenate([data, diag])
    return sparse.coo_matrix((data, (rows, cols)), shape=(net.Np, net.Np)).tocsr()

---

Solvers

voids.linalg.solve

solve_linear_system

solve_linear_system(
    A, b, *, method="direct", solver_parameters=None
)

Solve a sparse linear system with one of the supported SciPy backends.

Parameters:

Name	Type	Description	Default
A	csr_matrix	Sparse system matrix.	required
b	ndarray	Right-hand-side vector.	required
method	str	Solver backend. Supported values are "direct", "pardiso", "cg", and "gmres".	'direct'
solver_parameters	SolverParameters | None	Optional backend-specific solver options. For SciPy Krylov methods this maps directly to supported keyword arguments such as rtol, atol, restart, and maxiter. Setting {"preconditioner": "pyamg"} attaches a PyAMG preconditioner to cg or gmres.	None

Returns:

Type	Description
ndarray	Solution vector.
dict[str, Any]	Solver metadata containing the method name and the iterative solver status code info.

Raises:

Type	Description
ValueError	If method is not recognized.

Notes

The "pardiso" method uses the Intel MKL PARDISO solver via the pypardiso package. This is typically only available on Linux systems and may provide better performance for large systems compared to the default SciPy direct solver.

Source code in src/voids/linalg/solve.py

def solve_linear_system(
    A: sparse.csr_matrix,
    b: np.ndarray,
    *,
    method: str = "direct",
    solver_parameters: SolverParameters | None = None,
) -> tuple[np.ndarray, dict[str, str | float | int]]:
    """Solve a sparse linear system with one of the supported SciPy backends.

    Parameters
    ----------
    A :
        Sparse system matrix.
    b :
        Right-hand-side vector.
    method :
        Solver backend. Supported values are ``"direct"``, ``"pardiso"``,
        ``"cg"``, and ``"gmres"``.
    solver_parameters :
        Optional backend-specific solver options. For SciPy Krylov methods this
        maps directly to supported keyword arguments such as ``rtol``,
        ``atol``, ``restart``, and ``maxiter``. Setting
        ``{"preconditioner": "pyamg"}`` attaches a PyAMG preconditioner to
        ``cg`` or ``gmres``.

    Returns
    -------
    numpy.ndarray
        Solution vector.
    dict[str, Any]
        Solver metadata containing the method name and the iterative solver
        status code ``info``.

    Raises
    ------
    ValueError
        If ``method`` is not recognized.

    Notes
    -----
    The ``"pardiso"`` method uses the Intel MKL PARDISO solver via the
    ``pypardiso`` package. This is typically only available on Linux systems
    and may provide better performance for large systems compared to the
    default SciPy direct solver.
    """

    if method == "direct":
        x = spsolve(A, b)
        return np.asarray(x, dtype=float), {"method": method, "info": 0}
    if method == "pardiso":
        pardiso_spsolve = _import_pypardiso()
        x = pardiso_spsolve(A, b)
        return np.asarray(x, dtype=float), {"method": method, "info": 0}
    if method == "cg":
        parameters = dict(solver_parameters or {})
        preconditioner, preconditioner_info = _build_preconditioner(A, solver_parameters=parameters)
        cg_kwargs = {
            key: parameters[key] for key in ("rtol", "atol", "maxiter", "M") if key in parameters
        }
        if preconditioner is not None and "M" not in cg_kwargs:
            cg_kwargs["M"] = preconditioner
        x, info = cg(A, b, **cg_kwargs)
        return np.asarray(x, dtype=float), {
            "method": method,
            "info": int(info),
            **preconditioner_info,
        }
    if method == "gmres":
        parameters = dict(solver_parameters or {})
        preconditioner, preconditioner_info = _build_preconditioner(A, solver_parameters=parameters)
        gmres_kwargs = {
            key: parameters[key]
            for key in ("rtol", "atol", "restart", "maxiter", "M")
            if key in parameters
        }
        if preconditioner is not None and "M" not in gmres_kwargs:
            gmres_kwargs["M"] = preconditioner
        x, info = gmres(A, b, **gmres_kwargs)
        return np.asarray(x, dtype=float), {
            "method": method,
            "info": int(info),
            **preconditioner_info,
        }
    raise ValueError(f"Unknown solver method '{method}'")

---

Boundary Conditions

voids.linalg.bc

apply_dirichlet_rowcol

apply_dirichlet_rowcol(A, b, values, mask)

Apply Dirichlet conditions by row and column elimination.

Parameters:

Name	Type	Description	Default
A	csr_matrix	System matrix with shape (N, N).	required
b	ndarray	Right-hand-side vector with shape (N,).	required
values	ndarray	Full-length vector of prescribed values. Only entries selected by mask are enforced.	required
mask	ndarray	Boolean array selecting the Dirichlet degrees of freedom.	required

Returns:

Type	Description
csr_matrix	Modified system matrix in CSR format.
ndarray	Modified right-hand-side vector.

Raises:

Type	Description
ValueError	If values, mask, and b do not have the same shape.

Notes

For each constrained degree of freedom k, the routine enforces

A[k, :] = 0 A[:, k] = 0 A[k, k] = 1 b[k] = values[k]

after first subtracting the eliminated column contribution from the unconstrained rows of b.

Source code in src/voids/linalg/bc.py

def apply_dirichlet_rowcol(
    A: sparse.csr_matrix, b: np.ndarray, values: np.ndarray, mask: np.ndarray
) -> tuple[sparse.csr_matrix, np.ndarray]:
    """Apply Dirichlet conditions by row and column elimination.

    Parameters
    ----------
    A :
        System matrix with shape ``(N, N)``.
    b :
        Right-hand-side vector with shape ``(N,)``.
    values :
        Full-length vector of prescribed values. Only entries selected by
        ``mask`` are enforced.
    mask :
        Boolean array selecting the Dirichlet degrees of freedom.

    Returns
    -------
    scipy.sparse.csr_matrix
        Modified system matrix in CSR format.
    numpy.ndarray
        Modified right-hand-side vector.

    Raises
    ------
    ValueError
        If ``values``, ``mask``, and ``b`` do not have the same shape.

    Notes
    -----
    For each constrained degree of freedom ``k``, the routine enforces

    ``A[k, :] = 0``
    ``A[:, k] = 0``
    ``A[k, k] = 1``
    ``b[k] = values[k]``

    after first subtracting the eliminated column contribution from the
    unconstrained rows of ``b``.
    """

    A = A.tolil(copy=True)
    b2 = np.asarray(b, dtype=float).copy()
    values = np.asarray(values, dtype=float)
    mask = np.asarray(mask, dtype=bool)
    if values.shape != b2.shape or mask.shape != b2.shape:
        raise ValueError("values, mask and b must have the same shape")
    idx = np.flatnonzero(mask)
    if idx.size == 0:
        return A.tocsr(), b2

    A_csr = A.tocsr()
    b2 = b2 - A_csr[:, idx] @ values[idx]

    for k in idx:
        A[:, k] = 0.0
        A[k, :] = 0.0
        A[k, k] = 1.0
        b2[k] = values[k]
    return A.tocsr(), b2

---

Backends

voids.linalg.backends

SciPyBackend dataclass

Namespace collecting SciPy sparse constructors and solvers.

Attributes:

Name	Type	Description
coo_matrix, csr_matrix		Sparse matrix constructors.
spsolve, cg, gmres		Direct and iterative sparse linear solvers used by the package.

Source code in src/voids/linalg/backends.py

@dataclass(frozen=True, slots=True)
class SciPyBackend:
    """Namespace collecting SciPy sparse constructors and solvers.

    Attributes
    ----------
    coo_matrix, csr_matrix :
        Sparse matrix constructors.
    spsolve, cg, gmres :
        Direct and iterative sparse linear solvers used by the package.
    """

    coo_matrix = staticmethod(sparse.coo_matrix)
    csr_matrix = staticmethod(sparse.csr_matrix)
    spsolve = staticmethod(spsolve)
    cg = staticmethod(cg)
    gmres = staticmethod(gmres)

---

Diagnostics

voids.linalg.diagnostics

residual_norm

residual_norm(A, x, b)

Return the Euclidean norm of the linear-system residual.

Parameters:

Name	Type	Description	Default
A	csr_matrix	System matrix.	required
x	ndarray	Trial or converged solution vector.	required
b	ndarray	Right-hand-side vector.	required

Returns:

Type	Description
float	Value of ||A x - b||_2.

Source code in src/voids/linalg/diagnostics.py

def residual_norm(A: sparse.csr_matrix, x: np.ndarray, b: np.ndarray) -> float:
    """Return the Euclidean norm of the linear-system residual.

    Parameters
    ----------
    A :
        System matrix.
    x :
        Trial or converged solution vector.
    b :
        Right-hand-side vector.

    Returns
    -------
    float
        Value of ``||A x - b||_2``.
    """

    r = A @ x - b
    return float(np.linalg.norm(r))