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 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", "umfpack", "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 "direct" method uses :func:scipy.sparse.linalg.spsolve. The "umfpack" method requests SuiteSparse/UMFPACK explicitly through scikit-umfpack. The "pardiso" method uses Intel MKL PARDISO through pypardiso; this is typically only available on Linux systems.

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 backends.

    Parameters
    ----------
    A :
        Sparse system matrix.
    b :
        Right-hand-side vector.
    method :
        Solver backend. Supported values are ``"direct"``, ``"umfpack"``,
        ``"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 ``"direct"`` method uses :func:`scipy.sparse.linalg.spsolve`. The
    ``"umfpack"`` method requests SuiteSparse/UMFPACK explicitly through
    ``scikit-umfpack``. The ``"pardiso"`` method uses Intel MKL PARDISO through
    ``pypardiso``; this is typically only available on Linux systems.
    """

    if method == "direct":
        x = spsolve(A, b)
        return np.asarray(x, dtype=float), {
            "method": method,
            "backend": "scipy.sparse.linalg.spsolve",
            "info": 0,
        }
    if method == "umfpack":
        umfpack_spsolve = _import_umfpack()
        x = umfpack_spsolve(
            sparse.csc_matrix(A, dtype=float),
            np.ascontiguousarray(np.asarray(b, dtype=float)),
        )
        return np.asarray(x, dtype=float), {
            "method": method,
            "backend": "scikits.umfpack.spsolve",
            "info": 0,
        }
    if method == "pardiso":
        pardiso_spsolve = _import_pypardiso()
        x = pardiso_spsolve(A, b)
        return np.asarray(x, dtype=float), {"method": method, "backend": "pypardiso", "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))