Finite Elements

The voids.fem sub-package provides optional FEniCSx-backed finite-element single-phase solvers for porosity/permeability maps. These APIs require a compatible DOLFINx installation, such as the Pixi fem feature in this repository. The PyPI package does not install FEniCSx automatically.

The current single-phase FEM backends report effective permeability from the computed outlet flux and Darcy's law. They are numerical upscaling tools, not experimental validation claims by themselves.

The default FEniCSSolverOptions use PETSc LU with MUMPS. For high-contrast mixed Darcy-Brinkman maps, avoid treating PETSc's built-in LU backend as a scientific fallback unless it has been checked against a robust factorization package on the same problem class; it can return nonphysical permeabilities for these saddle-point systems.

The USFEM backend is especially sensitive to the external factorization package and workspace settings on larger 3-D maps. Record the PETSc options from the result metadata with any reported USFEM permeability.

Thread environment for FEM solves

For robust FEniCSx/PETSc direct solves, pin BLAS/OpenMP thread variables before Python imports NumPy, SciPy, PETSc, or DOLFINx. voids.fem applies conservative defaults when those variables are unset, but it does not override user-provided values. See the detailed warning in Map-Based Single-Phase Solvers.

For the governing equations, boundary conditions, spaces, stabilization terms, and permeability reporting convention, see Map-Based Single-Phase Solvers. For executable exact-solution refinement studies, expected convergence rates, and the centered-vug benchmark, see FEM Manufactured-Solution Verification.

---

Common Types

voids.fem.singlephase

Single-phase finite-element Darcy and Brinkman backends.

BrinkmanNondimensionalization dataclass

Velocity scaling choice for Brinkman nondimensional forms.

The internal unknowns are u* = u / U and p* = p / DeltaP. Use velocity_scale="viscous" for U = DeltaP L / mu. Use velocity_scale="unit_darcy" for U = DeltaP K / (mu L) on constant permeability maps. Results are converted back to physical velocity and pressure before permeability is reported.

Source code in src/voids/fem/singlephase/_common.py

@dataclass(frozen=True, slots=True)
class BrinkmanNondimensionalization:
    """Velocity scaling choice for Brinkman nondimensional forms.

    The internal unknowns are ``u* = u / U`` and ``p* = p / DeltaP``. Use
    ``velocity_scale="viscous"`` for ``U = DeltaP L / mu``. Use
    ``velocity_scale="unit_darcy"`` for ``U = DeltaP K / (mu L)`` on constant
    permeability maps. Results are converted back to physical velocity and
    pressure before permeability is reported.
    """

    velocity_scale: BrinkmanVelocityScale = "viscous"

FEMMapProblem dataclass

Porosity/permeability coefficient maps for FEM single-phase solves.

Parameters:

Name	Type	Description	Default
permeability_map	PermeabilityMap	Scalar cell-wise permeability map.	required
porosity_map	PorosityMap | None	Optional porosity map on the same grid. Brinkman solves use this field in nu_eff = mu / max(phi, porosity_floor). Darcy-only comparison solves do not use it.	None
viscosity	float	Dynamic viscosity mu.	1.0
porosity_floor	float	Lower bound used only in the Brinkman effective-viscosity coefficient.	1e-06
permeability_floor	float	Lower bound used in gamma = mu / max(K, permeability_floor).	1e-30

Source code in src/voids/fem/singlephase/_common.py

@dataclass(slots=True)
class FEMMapProblem:
    """Porosity/permeability coefficient maps for FEM single-phase solves.

    Parameters
    ----------
    permeability_map :
        Scalar cell-wise permeability map.
    porosity_map :
        Optional porosity map on the same grid. Brinkman solves use this field
        in ``nu_eff = mu / max(phi, porosity_floor)``. Darcy-only comparison
        solves do not use it.
    viscosity :
        Dynamic viscosity ``mu``.
    porosity_floor :
        Lower bound used only in the Brinkman effective-viscosity coefficient.
    permeability_floor :
        Lower bound used in ``gamma = mu / max(K, permeability_floor)``.
    """

    permeability_map: PermeabilityMap
    porosity_map: PorosityMap | None = None
    viscosity: float = 1.0
    porosity_floor: float = 1.0e-6
    permeability_floor: float = 1.0e-30

    def __post_init__(self) -> None:
        if self.viscosity <= 0.0 or not np.isfinite(self.viscosity):
            raise ValueError("viscosity must be positive and finite")
        if self.porosity_floor <= 0.0 or not np.isfinite(self.porosity_floor):
            raise ValueError("porosity_floor must be positive and finite")
        if self.permeability_floor <= 0.0 or not np.isfinite(self.permeability_floor):
            raise ValueError("permeability_floor must be positive and finite")
        if self.permeability_map.ndim not in {2, 3}:
            raise ValueError("permeability_map must be 2D or 3D")
        if self.porosity_map is not None:
            if self.porosity_map.shape != self.permeability_map.shape:
                raise ValueError("porosity_map and permeability_map must have the same shape")
            porosity_cell_size = tuple(
                float(v) for v in cast(tuple[float, ...], self.porosity_map.cell_size)
            )
            if porosity_cell_size != _cell_size_tuple(self.permeability_map):
                raise ValueError("porosity_map and permeability_map must have the same cell_size")

FEMSinglePhaseResult dataclass

Finite-element single-phase flow result.

Source code in src/voids/fem/singlephase/_common.py

@dataclass(slots=True)
class FEMSinglePhaseResult:
    """Finite-element single-phase flow result."""

    method: str
    formulation: str
    flow_axis: str
    permeability: float
    flow_rate: float
    pressure_inlet: float
    pressure_outlet: float
    pressure_drop: float
    viscosity: float
    domain_length: float
    cross_section_area: float
    solve_seconds: float
    velocity: Any
    pressure: Any
    metadata: dict[str, Any] = field(default_factory=dict)

FEniCSSolverOptions dataclass

Linear solver controls for FEniCSx linear problems.

The default linear_backend="auto" preserves the PETSc/MUMPS path on platforms with a full DOLFINx/PETSc stack. On native Windows, where that PETSc stack is not available in the conda-forge FEniCSx packages used by voids, auto uses DOLFINx assembly plus SciPy/SuperLU.

Use linear_backend="superlu", "scipy", "umfpack", "pardiso", or "nvmath_cudss" to request the serial DOLFINx-assembly/direct-sparse path explicitly. These paths use the same weak form and boundary conditions as PETSc; only the linear algebra backend changes. "scipy" is kept as a backward-compatible alias for the SciPy/SuperLU path. "umfpack" requires the optional scikits.umfpack package, "pardiso" requires the optional pypardiso package, and "nvmath_cudss" requires a CUDA-capable PyTorch/nvmath cuDSS runtime. linear_system_dtype selects the assembled value dtype for serial sparse backends. PETSc backends use the scalar type of the installed PETSc/DOLFINx stack and are therefore double precision in the supported Pixi FEM environment.

Source code in src/voids/fem/singlephase/_common.py

@dataclass(slots=True)
class FEniCSSolverOptions:
    """Linear solver controls for FEniCSx linear problems.

    The default ``linear_backend="auto"`` preserves the PETSc/MUMPS path on
    platforms with a full DOLFINx/PETSc stack. On native Windows, where that
    PETSc stack is not available in the conda-forge FEniCSx packages used by
    ``voids``, ``auto`` uses DOLFINx assembly plus SciPy/SuperLU.

    Use ``linear_backend="superlu"``, ``"scipy"``, ``"umfpack"``,
    ``"pardiso"``, or ``"nvmath_cudss"`` to request the serial
    DOLFINx-assembly/direct-sparse path explicitly. These paths use the same
    weak form and boundary conditions as PETSc; only the linear algebra backend
    changes. ``"scipy"`` is kept as a backward-compatible alias for the
    SciPy/SuperLU path. ``"umfpack"`` requires the optional
    ``scikits.umfpack`` package, ``"pardiso"`` requires the optional
    ``pypardiso`` package, and ``"nvmath_cudss"`` requires a CUDA-capable
    PyTorch/nvmath cuDSS runtime. ``linear_system_dtype`` selects the assembled
    value dtype for serial sparse backends. PETSc backends use the scalar type
    of the installed PETSc/DOLFINx stack and are therefore double precision in
    the supported Pixi FEM environment.
    """

    linear_backend: LinearSolverBackend = "auto"
    solver_preset: FEMSolverPreset = "direct_reference"
    linear_system_dtype: LinearSystemDType = "float64"
    superlu_controls: dict[str, Any] = field(default_factory=dict)
    umfpack_controls: dict[str, Any] = field(default_factory=dict)
    nvmath_cudss_controls: dict[str, Any] = field(default_factory=dict)
    iterative_solver_controls: dict[str, Any] = field(default_factory=dict)
    petsc_options: dict[str, Any] = field(
        default_factory=lambda: {
            "ksp_type": "preonly",
            "pc_type": "lu",
            "pc_factor_mat_solver_type": "mumps",
            "pc_factor_shift_type": "nonzero",
            "ksp_error_if_not_converged": True,
        }
    )
    petsc_options_prefix: str = "voids_fem_"

    @classmethod
    def direct_lu(
        cls,
        backend: str = "mumps",
        *,
        linear_backend: LinearSolverBackend = "petsc",
        petsc_options_prefix: str = "voids_fem_",
        shift_amount: float | None = 1.0e-12,
        mumps_memory_relaxation_percent: int | None = None,
        mumps_workspace_mb: int | None = None,
    ) -> FEniCSSolverOptions:
        """Create PETSc options for a direct sparse LU solve.

        Parameters
        ----------
        backend :
            PETSc factorization package, for example ``"mumps"`` or
            ``"superlu_dist"``.
        linear_backend :
            Linear algebra backend. This builder configures PETSc options, so
            the default is ``"petsc"``.
        petsc_options_prefix :
            Prefix used by DOLFINx for PETSc runtime options.
        shift_amount :
            Nonzero diagonal shift used during factorization. Pass ``None`` to
            omit the shift options.
        mumps_memory_relaxation_percent, mumps_workspace_mb :
            Optional MUMPS memory controls. They are added only when the backend
            is ``"mumps"``.
        """

        options: dict[str, Any] = {
            "ksp_type": "preonly",
            "pc_type": "lu",
            "pc_factor_mat_solver_type": backend,
            "ksp_error_if_not_converged": True,
        }
        if shift_amount is not None:
            options["pc_factor_shift_type"] = "nonzero"
            options["pc_factor_shift_amount"] = float(shift_amount)
        if backend == "mumps":
            if mumps_memory_relaxation_percent is not None:
                options["mat_mumps_icntl_14"] = int(mumps_memory_relaxation_percent)
            if mumps_workspace_mb is not None:
                options["mat_mumps_icntl_23"] = int(mumps_workspace_mb)
        return cls(
            linear_backend=linear_backend,
            solver_preset="direct_reference",
            petsc_options=options,
            petsc_options_prefix=petsc_options_prefix,
        )

    @classmethod
    def direct_reference(
        cls,
        backend: str = "mumps",
        *,
        petsc_options_prefix: str = "voids_fem_",
        shift_amount: float | None = 1.0e-12,
        mumps_memory_relaxation_percent: int | None = None,
        mumps_workspace_mb: int | None = None,
    ) -> FEniCSSolverOptions:
        """Create the conservative PETSc direct-solve preset.

        This preset is the stable baseline for permeability comparisons: one
        monolithic sparse LU factorization through PETSc, normally MUMPS. It is
        appropriate as the direct reference before testing faster iterative or
        distributed configurations.
        """

        return cls.direct_lu(
            backend,
            petsc_options_prefix=petsc_options_prefix,
            shift_amount=shift_amount,
            mumps_memory_relaxation_percent=mumps_memory_relaxation_percent,
            mumps_workspace_mb=mumps_workspace_mb,
        )

    @classmethod
    def direct_parallel(
        cls,
        backend: str = "mumps",
        *,
        petsc_options_prefix: str = "voids_fem_",
        shift_amount: float | None = 1.0e-12,
        mumps_memory_relaxation_percent: int | None = 500,
        mumps_workspace_mb: int | None = None,
    ) -> FEniCSSolverOptions:
        """Create an MPI-oriented PETSc direct-solve preset.

        This still uses a direct factorization, but records a distinct preset
        name and applies a larger default MUMPS memory-relaxation factor. It is
        intended for runs launched under MPI, where MUMPS or SuperLU_DIST can
        distribute the factorization.
        """

        options = cls.direct_lu(
            backend,
            petsc_options_prefix=petsc_options_prefix,
            shift_amount=shift_amount,
            mumps_memory_relaxation_percent=(
                mumps_memory_relaxation_percent if backend == "mumps" else None
            ),
            mumps_workspace_mb=mumps_workspace_mb if backend == "mumps" else None,
        )
        options.solver_preset = "direct_parallel"
        return options

    @classmethod
    def iterative_fieldsplit_experimental(
        cls,
        *,
        petsc_options_prefix: str = "voids_fem_",
        ksp_type: str = "fgmres",
        rtol: float = 1.0e-8,
        max_it: int = 500,
    ) -> FEniCSSolverOptions:
        """Create an experimental PETSc field-split preset for mixed systems.

        The preset is deliberately labelled experimental because scalable
        saddle-point preconditioning depends on the formulation, coefficient
        contrast, mesh, and PETSc build. Permeability results from this preset
        should be compared against a direct reference before being used for
        scientific conclusions.
        """

        return cls(
            linear_backend="petsc",
            solver_preset="iterative_fieldsplit_experimental",
            petsc_options={
                "ksp_type": ksp_type,
                "ksp_rtol": float(rtol),
                "ksp_max_it": int(max_it),
                "ksp_error_if_not_converged": True,
                "pc_type": "fieldsplit",
                "pc_fieldsplit_type": "schur",
                "pc_fieldsplit_schur_fact_type": "upper",
                "pc_fieldsplit_detect_saddle_point": True,
                "fieldsplit_0_ksp_type": "preonly",
                "fieldsplit_0_pc_type": "hypre",
                "fieldsplit_1_ksp_type": "preonly",
                "fieldsplit_1_pc_type": "jacobi",
            },
            petsc_options_prefix=petsc_options_prefix,
        )

    @classmethod
    def iterative_block_lgmres_experimental(
        cls,
        *,
        petsc_options_prefix: str = "voids_fem_",
        rtol: float = 1.0e-8,
        atol: float = 1.0e-10,
        max_it: int = 3000,
        block_lu_backend: str = "superlu_dist",
    ) -> FEniCSSolverOptions:
        """Create an experimental block-LGMRES preset for USFEM block solves.

        This preset is intended for ``solve_brinkman_usfem_block`` with
        ``matrix_kind="nest"`` and ``preconditioner="none"``. It uses a
        multiplicative velocity/pressure field split as the outer
        preconditioner, with direct LU subsolves on the two diagonal operator
        blocks. It is a correctness-oriented iterative baseline, not a scalable
        multigrid preconditioner.
        """

        return cls(
            linear_backend="petsc",
            solver_preset="iterative_block_lgmres_experimental",
            petsc_options={
                "ksp_type": "lgmres",
                "ksp_rtol": float(rtol),
                "ksp_atol": float(atol),
                "ksp_max_it": int(max_it),
                "ksp_norm_type": "unpreconditioned",
                "ksp_error_if_not_converged": True,
                "pc_type": "fieldsplit",
                "pc_fieldsplit_type": "multiplicative",
                "fieldsplit_u_0_ksp_type": "preonly",
                "fieldsplit_u_0_pc_type": "lu",
                "fieldsplit_u_0_pc_factor_mat_solver_type": block_lu_backend,
                "fieldsplit_p_1_ksp_type": "preonly",
                "fieldsplit_p_1_pc_type": "lu",
                "fieldsplit_p_1_pc_factor_mat_solver_type": block_lu_backend,
            },
            petsc_options_prefix=petsc_options_prefix,
        )

    @classmethod
    def usfem_schurdiag_cudss_experimental(
        cls,
        *,
        dtype: Literal["float32", "float64"] = "float64",
        device_ids: int | Sequence[int] | Literal["all"] | None = None,
        ir_steps: int = 5,
        use_matching: bool = True,
        rtol: float = 1.0e-8,
        atol: float = 0.0,
        max_it: int = 1000,
        restart: int = 200,
        velocity_solver: Literal["amg", "exact"] = "amg",
        schur_drop_rel: float = 0.0,
        error_if_not_converged: bool = True,
        reordering_alg: str | int | None = None,
        matching_alg: str | int | None = None,
        factorization_alg: str | int | None = None,
        solve_alg: str | int | None = None,
        pivot_type: Literal["col", "row", "none"] | None = None,
        pivot_threshold: float | None = None,
        pivot_epsilon: float | None = None,
        pivot_epsilon_alg: str | int | None = None,
        nd_nlevels: int | None = None,
        host_nthreads: int | None = None,
        threading_lib: str | None = None,
        hybrid_mode: bool | None = None,
        hybrid_device_memory_limit: int | None = None,
        hybrid_execute_mode: bool | None = None,
        use_cuda_register_memory: bool | None = None,
        use_superpanels: bool | None = None,
        deterministic_mode: bool | None = None,
        residual_rtol: float | None = None,
        controls: Mapping[str, Any] | None = None,
    ) -> FEniCSSolverOptions:
        """Create the experimental USFEM Schurdiag/cuDSS iterative preset.

        This preset is consumed by :func:`solve_brinkman_usfem_block`. It uses
        PETSc/DOLFINx only to assemble the USFEM velocity/pressure block forms,
        converts the blocks to SciPy CSR matrices in a single Python process,
        and solves the coupled system with SciPy GMRES. The preconditioner is a
        lower-Schur block form: the velocity block is approximately inverted by
        PyAMG by default, while the diagonal-velocity Schur approximation
        ``S_hat = A_pp - A_pu diag(A_uu)^-1 A_up`` is factored once with
        nvmath/cuDSS and reused for every pressure correction.

        The preset is experimental and should be compared against a same-map
        direct reference before scientific use. Local 3D USFEM Brinkman probes
        found ``dtype="float64"`` reliable for permeability and field recovery;
        ``dtype="float32"`` reduced memory but did not reliably preserve the
        pressure field on high-contrast maps.
        """

        if velocity_solver not in {"amg", "exact"}:
            raise ValueError("velocity_solver must be either 'amg' or 'exact'")
        if rtol <= 0.0 or not np.isfinite(rtol):
            raise ValueError("rtol must be positive and finite")
        if atol < 0.0 or not np.isfinite(atol):
            raise ValueError("atol must be non-negative and finite")
        if max_it <= 0:
            raise ValueError("max_it must be positive")
        if restart <= 0:
            raise ValueError("restart must be positive")
        if schur_drop_rel < 0.0 or not np.isfinite(schur_drop_rel):
            raise ValueError("schur_drop_rel must be non-negative and finite")

        resolved_controls = _nvmath_cudss_controls_from_arguments(
            dtype=dtype,
            device_ids=device_ids,
            ir_steps=ir_steps,
            use_matching=use_matching,
            reordering_alg=reordering_alg,
            matching_alg=matching_alg,
            factorization_alg=factorization_alg,
            solve_alg=solve_alg,
            pivot_type=pivot_type,
            pivot_threshold=pivot_threshold,
            pivot_epsilon=pivot_epsilon,
            pivot_epsilon_alg=pivot_epsilon_alg,
            nd_nlevels=nd_nlevels,
            host_nthreads=host_nthreads,
            threading_lib=threading_lib,
            hybrid_mode=hybrid_mode,
            hybrid_device_memory_limit=hybrid_device_memory_limit,
            hybrid_execute_mode=hybrid_execute_mode,
            use_cuda_register_memory=use_cuda_register_memory,
            use_superpanels=use_superpanels,
            deterministic_mode=deterministic_mode,
            check_residual=False,
            residual_rtol=residual_rtol,
            controls=controls,
        )
        return cls(
            linear_backend="petsc",
            solver_preset="iterative_schurdiag_cudss_experimental",
            linear_system_dtype=cast(LinearSystemDType, resolved_controls["dtype"]),
            nvmath_cudss_controls=resolved_controls,
            iterative_solver_controls={
                "gmres_rtol": float(rtol),
                "gmres_atol": float(atol),
                "gmres_maxiter": int(max_it),
                "gmres_restart": int(restart),
                "velocity_solver": velocity_solver,
                "schurdiag_drop_rel": float(schur_drop_rel),
                "error_if_not_converged": bool(error_if_not_converged),
            },
            petsc_options={},
        )

    @classmethod
    def scipy_direct(
        cls,
        *,
        linear_system_dtype: LinearSystemDType = "float64",
        permc_spec: str | None = None,
        diag_pivot_thresh: float | None = None,
        relax: int | None = None,
        panel_size: int | None = None,
        equil: bool | None = None,
        controls: Mapping[str, Any] | None = None,
    ) -> FEniCSSolverOptions:
        """Create options for the serial DOLFINx-assembly/SciPy SuperLU backend."""

        return cls(
            linear_backend="scipy",
            solver_preset="direct_reference",
            linear_system_dtype=linear_system_dtype,
            superlu_controls=_superlu_controls_from_arguments(
                permc_spec=permc_spec,
                diag_pivot_thresh=diag_pivot_thresh,
                relax=relax,
                panel_size=panel_size,
                equil=equil,
                controls=controls,
            ),
        )

    @classmethod
    def superlu_direct(
        cls,
        *,
        linear_system_dtype: LinearSystemDType = "float64",
        permc_spec: str | None = None,
        diag_pivot_thresh: float | None = None,
        relax: int | None = None,
        panel_size: int | None = None,
        equil: bool | None = None,
        controls: Mapping[str, Any] | None = None,
    ) -> FEniCSSolverOptions:
        """Create options for the serial DOLFINx-assembly/SuperLU backend."""

        return cls(
            linear_backend="superlu",
            solver_preset="direct_reference",
            linear_system_dtype=linear_system_dtype,
            superlu_controls=_superlu_controls_from_arguments(
                permc_spec=permc_spec,
                diag_pivot_thresh=diag_pivot_thresh,
                relax=relax,
                panel_size=panel_size,
                equil=equil,
                controls=controls,
            ),
        )

    @classmethod
    def umfpack_direct(
        cls,
        *,
        linear_system_dtype: LinearSystemDType = "float64",
        ordering: str | int | float | None = None,
        strategy: str | int | float | None = None,
        pivot_tolerance: float | None = None,
        sym_pivot_tolerance: float | None = None,
        scale: str | int | float | None = None,
        block_size: int | None = None,
        controls: Mapping[str, Any] | None = None,
    ) -> FEniCSSolverOptions:
        """Create options for the serial DOLFINx-assembly/UMFPACK backend."""

        umfpack_controls: dict[str, Any] = dict(controls or {})
        if ordering is not None:
            umfpack_controls["ordering"] = ordering
        if strategy is not None:
            umfpack_controls["strategy"] = strategy
        if pivot_tolerance is not None:
            umfpack_controls["pivot_tolerance"] = float(pivot_tolerance)
        if sym_pivot_tolerance is not None:
            umfpack_controls["sym_pivot_tolerance"] = float(sym_pivot_tolerance)
        if scale is not None:
            umfpack_controls["scale"] = scale
        if block_size is not None:
            umfpack_controls["block_size"] = int(block_size)
        return cls(
            linear_backend="umfpack",
            solver_preset="direct_reference",
            linear_system_dtype=linear_system_dtype,
            umfpack_controls=umfpack_controls,
        )

    @classmethod
    def nvmath_cudss_direct(
        cls,
        *,
        dtype: Literal["float32", "float64"] = "float64",
        device_ids: int | Sequence[int] | Literal["all"] | None = None,
        ir_steps: int = 5,
        use_matching: bool = True,
        reordering_alg: str | int | None = None,
        matching_alg: str | int | None = None,
        factorization_alg: str | int | None = None,
        solve_alg: str | int | None = None,
        pivot_type: Literal["col", "row", "none"] | None = None,
        pivot_threshold: float | None = None,
        pivot_epsilon: float | None = None,
        pivot_epsilon_alg: str | int | None = None,
        nd_nlevels: int | None = None,
        host_nthreads: int | None = None,
        threading_lib: str | None = None,
        hybrid_mode: bool | None = None,
        hybrid_device_memory_limit: int | None = None,
        hybrid_execute_mode: bool | None = None,
        use_cuda_register_memory: bool | None = None,
        use_superpanels: bool | None = None,
        deterministic_mode: bool | None = None,
        check_residual: bool = True,
        residual_rtol: float | None = None,
        controls: Mapping[str, Any] | None = None,
    ) -> FEniCSSolverOptions:
        """Create options for the optional CUDA cuDSS direct-solver backend.

        The backend uses serial DOLFINx assembly, converts the assembled system
        to SciPy CSR format, copies the sparse matrix and right-hand side to
        CUDA tensors, and solves the system through ``nvmath.bindings.cudss``.
        It requires PyTorch with CUDA support and the nvmath/cuDSS runtime at
        solve time. These controls affect only the numerical linear solve; they
        do not change the Brinkman/Darcy weak form, porosity, permeability,
        viscosity, pressure drop, or dimensional scaling.

        Parameters
        ----------
        dtype :
            Floating-point value precision used for the cuDSS matrix, right-hand
            side, and solution. ``"float64"`` is the conservative default.
            ``"float32"`` can reduce memory use and wall time, but it must be
            accepted only after the residual check and a same-map reference
            comparison because high permeability contrasts can be ill-conditioned.
        device_ids :
            CUDA device selection. ``None`` uses PyTorch's current CUDA device,
            an integer selects one device, a sequence such as ``(0, 1)`` requests
            a single-node multi-GPU cuDSS handle, and ``"all"`` uses all CUDA
            devices visible to PyTorch. Some cuDSS controls are not supported by
            all single-GPU and multi-GPU execution paths.
        ir_steps :
            cuDSS iterative-refinement step count
            (``ConfigParam.IR_N_STEPS``). ``voids`` sets this to ``5`` by
            default, although a fresh cuDSS config in the tested stack reports
            ``IR_N_STEPS = 0``. Increasing this value can improve single-precision
            residuals, but the effect is problem-dependent and not guaranteed to
            be monotonic.
        use_matching :
            Whether to request cuDSS matching/scaling
            (``ConfigParam.USE_MATCHING``). Matching is usually helpful for
            general sparse matrices because it can reduce pivot perturbations;
            it is enabled by default.
        reordering_alg, matching_alg, factorization_alg, solve_alg, pivot_epsilon_alg :
            Optional cuDSS algorithm selectors for
            ``REORDERING_ALG``, ``MATCHING_ALG``, ``FACTORIZATION_ALG``,
            ``SOLVE_ALG``, and ``PIVOT_EPSILON_ALG``. Values may be integers
            ``0`` through ``5`` or strings such as ``"default"``, ``"alg_1"``,
            or ``"3"``. Their exact meaning and support are defined by the
            installed cuDSS version; unsupported combinations raise a cuDSS
            error during analysis, factorization, or solve.
        pivot_type :
            Optional cuDSS pivoting mode (``ConfigParam.PIVOT_TYPE``): ``"col"``
            for column pivoting, ``"row"`` for row pivoting, or ``"none"`` to
            request no pivoting. Leave as ``None`` to use the cuDSS default.
        pivot_threshold :
            Optional non-negative pivoting threshold
            (``ConfigParam.PIVOT_THRESHOLD``). Larger values can make pivoting
            more conservative, but may increase cost or be ignored depending on
            the selected cuDSS algorithm.
        pivot_epsilon :
            Optional non-negative pivot perturbation/floor
            (``ConfigParam.PIVOT_EPSILON``). This can stabilize factorizations
            with very small pivots, especially in single precision. Treat this
            as a solver-stabilization experiment: high-contrast FEM and TPFA
            maps can pass a loose residual check while still failing same-map
            field comparisons in single precision.
        nd_nlevels :
            Optional non-negative nested-dissection level control
            (``ConfigParam.ND_NLEVELS``) used by cuDSS reordering algorithms that
            support it.
        host_nthreads :
            Optional positive host-thread count for cuDSS host-side work
            (``ConfigParam.HOST_NTHREADS``). This affects execution only when a
            cuDSS threading-layer library is loaded.
        threading_lib :
            Optional path to a cuDSS threading-layer library. Use ``"auto"`` to
            request the packaged ``libcudss_mtlayer_gomp`` library when
            available. If host threading is requested and this is omitted,
            ``voids`` uses ``CUDSS_THREADING_LIB`` when set and otherwise
            auto-loads the packaged threading layer when ``host_nthreads`` or
            ``hybrid_execute_mode=True`` is requested.
        hybrid_mode :
            Optional request to enable cuDSS hybrid memory mode
            (``ConfigParam.HYBRID_MODE``), which can spill factorization data to
            host memory when device memory is limiting.
        hybrid_device_memory_limit :
            Optional positive device-memory limit in bytes for cuDSS hybrid
            memory mode (``ConfigParam.HYBRID_DEVICE_MEMORY_LIMIT``). ``voids``
            applies this after analysis and before factorization. Multi-GPU
            hybrid-memory limit handling is runtime dependent in the low-level
            nvmath binding; compare against a small same-configuration probe
            before relying on it.
        hybrid_execute_mode :
            Optional request to enable cuDSS hybrid execute mode
            (``ConfigParam.HYBRID_EXECUTE_MODE``). Support is runtime/backend
            dependent.
        use_cuda_register_memory :
            Optional request to register host memory with CUDA
            (``ConfigParam.USE_CUDA_REGISTER_MEMORY``) for hybrid-memory runs.
        use_superpanels :
            Optional flag for cuDSS superpanel optimization
            (``ConfigParam.USE_SUPERPANELS``). Leave as ``None`` to use the cuDSS
            default.
        deterministic_mode :
            Optional request for deterministic cuDSS execution
            (``ConfigParam.DETERMINISTIC_MODE``). This is version/backend
            dependent and may raise ``NOT_SUPPORTED`` in some cuDSS paths.
        check_residual :
            If true, ``voids`` computes the assembled-system relative residual
            after cuDSS returns and raises an error when it exceeds
            ``residual_rtol``. Keep this enabled for lower-precision studies.
        residual_rtol :
            Relative residual acceptance tolerance used when
            ``check_residual=True``. If omitted, ``voids`` uses ``1.0e-8`` for
            ``float64`` and ``1.0e-4`` for ``float32``.
        controls :
            Optional low-level control mapping. Keys use the same names as the
            keyword arguments, with hyphens normalized to underscores. Explicit
            keyword arguments override the same keys in this mapping.

        Returns
        -------
        FEniCSSolverOptions
            Solver options with ``linear_backend="nvmath_cudss"`` and resolved
            cuDSS controls stored in ``nvmath_cudss_controls``.

        Notes
        -----
        This backend is experimental and hardware/runtime dependent. It should
        be compared against a same-map direct reference before using the result
        for scientific claims.
        """

        resolved_controls = _nvmath_cudss_controls_from_arguments(
            dtype=dtype,
            device_ids=device_ids,
            ir_steps=ir_steps,
            use_matching=use_matching,
            reordering_alg=reordering_alg,
            matching_alg=matching_alg,
            factorization_alg=factorization_alg,
            solve_alg=solve_alg,
            pivot_type=pivot_type,
            pivot_threshold=pivot_threshold,
            pivot_epsilon=pivot_epsilon,
            pivot_epsilon_alg=pivot_epsilon_alg,
            nd_nlevels=nd_nlevels,
            host_nthreads=host_nthreads,
            threading_lib=threading_lib,
            hybrid_mode=hybrid_mode,
            hybrid_device_memory_limit=hybrid_device_memory_limit,
            hybrid_execute_mode=hybrid_execute_mode,
            use_cuda_register_memory=use_cuda_register_memory,
            use_superpanels=use_superpanels,
            deterministic_mode=deterministic_mode,
            check_residual=check_residual,
            residual_rtol=residual_rtol,
            controls=controls,
        )
        return cls(
            linear_backend="nvmath_cudss",
            solver_preset="direct_reference",
            linear_system_dtype=cast(LinearSystemDType, resolved_controls["dtype"]),
            nvmath_cudss_controls=resolved_controls,
        )

    @classmethod
    def pardiso_direct(
        cls,
        *,
        linear_system_dtype: LinearSystemDType = "float64",
    ) -> FEniCSSolverOptions:
        """Create options for the serial DOLFINx-assembly/PARDISO backend."""

        return cls(
            linear_backend="pardiso",
            solver_preset="direct_reference",
            linear_system_dtype=linear_system_dtype,
        )

direct_lu classmethod

direct_lu(
    backend="mumps",
    *,
    linear_backend="petsc",
    petsc_options_prefix="voids_fem_",
    shift_amount=1e-12,
    mumps_memory_relaxation_percent=None,
    mumps_workspace_mb=None,
)

Create PETSc options for a direct sparse LU solve.

Parameters:

Name	Type	Description	Default
backend	str	PETSc factorization package, for example "mumps" or "superlu_dist".	'mumps'
linear_backend	LinearSolverBackend	Linear algebra backend. This builder configures PETSc options, so the default is "petsc".	'petsc'
petsc_options_prefix	str	Prefix used by DOLFINx for PETSc runtime options.	'voids_fem_'
shift_amount	float | None	Nonzero diagonal shift used during factorization. Pass None to omit the shift options.	1e-12
mumps_memory_relaxation_percent	int | None	Optional MUMPS memory controls. They are added only when the backend is "mumps".	None
mumps_workspace_mb	int | None	Optional MUMPS memory controls. They are added only when the backend is "mumps".	None

Source code in src/voids/fem/singlephase/_common.py

@classmethod
def direct_lu(
    cls,
    backend: str = "mumps",
    *,
    linear_backend: LinearSolverBackend = "petsc",
    petsc_options_prefix: str = "voids_fem_",
    shift_amount: float | None = 1.0e-12,
    mumps_memory_relaxation_percent: int | None = None,
    mumps_workspace_mb: int | None = None,
) -> FEniCSSolverOptions:
    """Create PETSc options for a direct sparse LU solve.

    Parameters
    ----------
    backend :
        PETSc factorization package, for example ``"mumps"`` or
        ``"superlu_dist"``.
    linear_backend :
        Linear algebra backend. This builder configures PETSc options, so
        the default is ``"petsc"``.
    petsc_options_prefix :
        Prefix used by DOLFINx for PETSc runtime options.
    shift_amount :
        Nonzero diagonal shift used during factorization. Pass ``None`` to
        omit the shift options.
    mumps_memory_relaxation_percent, mumps_workspace_mb :
        Optional MUMPS memory controls. They are added only when the backend
        is ``"mumps"``.
    """

    options: dict[str, Any] = {
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": backend,
        "ksp_error_if_not_converged": True,
    }
    if shift_amount is not None:
        options["pc_factor_shift_type"] = "nonzero"
        options["pc_factor_shift_amount"] = float(shift_amount)
    if backend == "mumps":
        if mumps_memory_relaxation_percent is not None:
            options["mat_mumps_icntl_14"] = int(mumps_memory_relaxation_percent)
        if mumps_workspace_mb is not None:
            options["mat_mumps_icntl_23"] = int(mumps_workspace_mb)
    return cls(
        linear_backend=linear_backend,
        solver_preset="direct_reference",
        petsc_options=options,
        petsc_options_prefix=petsc_options_prefix,
    )

direct_reference classmethod

direct_reference(
    backend="mumps",
    *,
    petsc_options_prefix="voids_fem_",
    shift_amount=1e-12,
    mumps_memory_relaxation_percent=None,
    mumps_workspace_mb=None,
)

Create the conservative PETSc direct-solve preset.

This preset is the stable baseline for permeability comparisons: one monolithic sparse LU factorization through PETSc, normally MUMPS. It is appropriate as the direct reference before testing faster iterative or distributed configurations.

Source code in src/voids/fem/singlephase/_common.py

@classmethod
def direct_reference(
    cls,
    backend: str = "mumps",
    *,
    petsc_options_prefix: str = "voids_fem_",
    shift_amount: float | None = 1.0e-12,
    mumps_memory_relaxation_percent: int | None = None,
    mumps_workspace_mb: int | None = None,
) -> FEniCSSolverOptions:
    """Create the conservative PETSc direct-solve preset.

    This preset is the stable baseline for permeability comparisons: one
    monolithic sparse LU factorization through PETSc, normally MUMPS. It is
    appropriate as the direct reference before testing faster iterative or
    distributed configurations.
    """

    return cls.direct_lu(
        backend,
        petsc_options_prefix=petsc_options_prefix,
        shift_amount=shift_amount,
        mumps_memory_relaxation_percent=mumps_memory_relaxation_percent,
        mumps_workspace_mb=mumps_workspace_mb,
    )

direct_parallel classmethod

direct_parallel(
    backend="mumps",
    *,
    petsc_options_prefix="voids_fem_",
    shift_amount=1e-12,
    mumps_memory_relaxation_percent=500,
    mumps_workspace_mb=None,
)

Create an MPI-oriented PETSc direct-solve preset.

This still uses a direct factorization, but records a distinct preset name and applies a larger default MUMPS memory-relaxation factor. It is intended for runs launched under MPI, where MUMPS or SuperLU_DIST can distribute the factorization.

Source code in src/voids/fem/singlephase/_common.py

@classmethod
def direct_parallel(
    cls,
    backend: str = "mumps",
    *,
    petsc_options_prefix: str = "voids_fem_",
    shift_amount: float | None = 1.0e-12,
    mumps_memory_relaxation_percent: int | None = 500,
    mumps_workspace_mb: int | None = None,
) -> FEniCSSolverOptions:
    """Create an MPI-oriented PETSc direct-solve preset.

    This still uses a direct factorization, but records a distinct preset
    name and applies a larger default MUMPS memory-relaxation factor. It is
    intended for runs launched under MPI, where MUMPS or SuperLU_DIST can
    distribute the factorization.
    """

    options = cls.direct_lu(
        backend,
        petsc_options_prefix=petsc_options_prefix,
        shift_amount=shift_amount,
        mumps_memory_relaxation_percent=(
            mumps_memory_relaxation_percent if backend == "mumps" else None
        ),
        mumps_workspace_mb=mumps_workspace_mb if backend == "mumps" else None,
    )
    options.solver_preset = "direct_parallel"
    return options

iterative_fieldsplit_experimental classmethod

iterative_fieldsplit_experimental(
    *,
    petsc_options_prefix="voids_fem_",
    ksp_type="fgmres",
    rtol=1e-08,
    max_it=500,
)

Create an experimental PETSc field-split preset for mixed systems.

The preset is deliberately labelled experimental because scalable saddle-point preconditioning depends on the formulation, coefficient contrast, mesh, and PETSc build. Permeability results from this preset should be compared against a direct reference before being used for scientific conclusions.

Source code in src/voids/fem/singlephase/_common.py

@classmethod
def iterative_fieldsplit_experimental(
    cls,
    *,
    petsc_options_prefix: str = "voids_fem_",
    ksp_type: str = "fgmres",
    rtol: float = 1.0e-8,
    max_it: int = 500,
) -> FEniCSSolverOptions:
    """Create an experimental PETSc field-split preset for mixed systems.

    The preset is deliberately labelled experimental because scalable
    saddle-point preconditioning depends on the formulation, coefficient
    contrast, mesh, and PETSc build. Permeability results from this preset
    should be compared against a direct reference before being used for
    scientific conclusions.
    """

    return cls(
        linear_backend="petsc",
        solver_preset="iterative_fieldsplit_experimental",
        petsc_options={
            "ksp_type": ksp_type,
            "ksp_rtol": float(rtol),
            "ksp_max_it": int(max_it),
            "ksp_error_if_not_converged": True,
            "pc_type": "fieldsplit",
            "pc_fieldsplit_type": "schur",
            "pc_fieldsplit_schur_fact_type": "upper",
            "pc_fieldsplit_detect_saddle_point": True,
            "fieldsplit_0_ksp_type": "preonly",
            "fieldsplit_0_pc_type": "hypre",
            "fieldsplit_1_ksp_type": "preonly",
            "fieldsplit_1_pc_type": "jacobi",
        },
        petsc_options_prefix=petsc_options_prefix,
    )

iterative_block_lgmres_experimental classmethod

iterative_block_lgmres_experimental(
    *,
    petsc_options_prefix="voids_fem_",
    rtol=1e-08,
    atol=1e-10,
    max_it=3000,
    block_lu_backend="superlu_dist",
)

Create an experimental block-LGMRES preset for USFEM block solves.

This preset is intended for solve_brinkman_usfem_block with matrix_kind="nest" and preconditioner="none". It uses a multiplicative velocity/pressure field split as the outer preconditioner, with direct LU subsolves on the two diagonal operator blocks. It is a correctness-oriented iterative baseline, not a scalable multigrid preconditioner.

Source code in src/voids/fem/singlephase/_common.py

@classmethod
def iterative_block_lgmres_experimental(
    cls,
    *,
    petsc_options_prefix: str = "voids_fem_",
    rtol: float = 1.0e-8,
    atol: float = 1.0e-10,
    max_it: int = 3000,
    block_lu_backend: str = "superlu_dist",
) -> FEniCSSolverOptions:
    """Create an experimental block-LGMRES preset for USFEM block solves.

    This preset is intended for ``solve_brinkman_usfem_block`` with
    ``matrix_kind="nest"`` and ``preconditioner="none"``. It uses a
    multiplicative velocity/pressure field split as the outer
    preconditioner, with direct LU subsolves on the two diagonal operator
    blocks. It is a correctness-oriented iterative baseline, not a scalable
    multigrid preconditioner.
    """

    return cls(
        linear_backend="petsc",
        solver_preset="iterative_block_lgmres_experimental",
        petsc_options={
            "ksp_type": "lgmres",
            "ksp_rtol": float(rtol),
            "ksp_atol": float(atol),
            "ksp_max_it": int(max_it),
            "ksp_norm_type": "unpreconditioned",
            "ksp_error_if_not_converged": True,
            "pc_type": "fieldsplit",
            "pc_fieldsplit_type": "multiplicative",
            "fieldsplit_u_0_ksp_type": "preonly",
            "fieldsplit_u_0_pc_type": "lu",
            "fieldsplit_u_0_pc_factor_mat_solver_type": block_lu_backend,
            "fieldsplit_p_1_ksp_type": "preonly",
            "fieldsplit_p_1_pc_type": "lu",
            "fieldsplit_p_1_pc_factor_mat_solver_type": block_lu_backend,
        },
        petsc_options_prefix=petsc_options_prefix,
    )

usfem_schurdiag_cudss_experimental classmethod

usfem_schurdiag_cudss_experimental(
    *,
    dtype="float64",
    device_ids=None,
    ir_steps=5,
    use_matching=True,
    rtol=1e-08,
    atol=0.0,
    max_it=1000,
    restart=200,
    velocity_solver="amg",
    schur_drop_rel=0.0,
    error_if_not_converged=True,
    reordering_alg=None,
    matching_alg=None,
    factorization_alg=None,
    solve_alg=None,
    pivot_type=None,
    pivot_threshold=None,
    pivot_epsilon=None,
    pivot_epsilon_alg=None,
    nd_nlevels=None,
    host_nthreads=None,
    threading_lib=None,
    hybrid_mode=None,
    hybrid_device_memory_limit=None,
    hybrid_execute_mode=None,
    use_cuda_register_memory=None,
    use_superpanels=None,
    deterministic_mode=None,
    residual_rtol=None,
    controls=None,
)

Create the experimental USFEM Schurdiag/cuDSS iterative preset.

This preset is consumed by :func:solve_brinkman_usfem_block. It uses PETSc/DOLFINx only to assemble the USFEM velocity/pressure block forms, converts the blocks to SciPy CSR matrices in a single Python process, and solves the coupled system with SciPy GMRES. The preconditioner is a lower-Schur block form: the velocity block is approximately inverted by PyAMG by default, while the diagonal-velocity Schur approximation S_hat = A_pp - A_pu diag(A_uu)^-1 A_up is factored once with nvmath/cuDSS and reused for every pressure correction.

The preset is experimental and should be compared against a same-map direct reference before scientific use. Local 3D USFEM Brinkman probes found dtype="float64" reliable for permeability and field recovery; dtype="float32" reduced memory but did not reliably preserve the pressure field on high-contrast maps.

Source code in src/voids/fem/singlephase/_common.py

@classmethod
def usfem_schurdiag_cudss_experimental(
    cls,
    *,
    dtype: Literal["float32", "float64"] = "float64",
    device_ids: int | Sequence[int] | Literal["all"] | None = None,
    ir_steps: int = 5,
    use_matching: bool = True,
    rtol: float = 1.0e-8,
    atol: float = 0.0,
    max_it: int = 1000,
    restart: int = 200,
    velocity_solver: Literal["amg", "exact"] = "amg",
    schur_drop_rel: float = 0.0,
    error_if_not_converged: bool = True,
    reordering_alg: str | int | None = None,
    matching_alg: str | int | None = None,
    factorization_alg: str | int | None = None,
    solve_alg: str | int | None = None,
    pivot_type: Literal["col", "row", "none"] | None = None,
    pivot_threshold: float | None = None,
    pivot_epsilon: float | None = None,
    pivot_epsilon_alg: str | int | None = None,
    nd_nlevels: int | None = None,
    host_nthreads: int | None = None,
    threading_lib: str | None = None,
    hybrid_mode: bool | None = None,
    hybrid_device_memory_limit: int | None = None,
    hybrid_execute_mode: bool | None = None,
    use_cuda_register_memory: bool | None = None,
    use_superpanels: bool | None = None,
    deterministic_mode: bool | None = None,
    residual_rtol: float | None = None,
    controls: Mapping[str, Any] | None = None,
) -> FEniCSSolverOptions:
    """Create the experimental USFEM Schurdiag/cuDSS iterative preset.

    This preset is consumed by :func:`solve_brinkman_usfem_block`. It uses
    PETSc/DOLFINx only to assemble the USFEM velocity/pressure block forms,
    converts the blocks to SciPy CSR matrices in a single Python process,
    and solves the coupled system with SciPy GMRES. The preconditioner is a
    lower-Schur block form: the velocity block is approximately inverted by
    PyAMG by default, while the diagonal-velocity Schur approximation
    ``S_hat = A_pp - A_pu diag(A_uu)^-1 A_up`` is factored once with
    nvmath/cuDSS and reused for every pressure correction.

    The preset is experimental and should be compared against a same-map
    direct reference before scientific use. Local 3D USFEM Brinkman probes
    found ``dtype="float64"`` reliable for permeability and field recovery;
    ``dtype="float32"`` reduced memory but did not reliably preserve the
    pressure field on high-contrast maps.
    """

    if velocity_solver not in {"amg", "exact"}:
        raise ValueError("velocity_solver must be either 'amg' or 'exact'")
    if rtol <= 0.0 or not np.isfinite(rtol):
        raise ValueError("rtol must be positive and finite")
    if atol < 0.0 or not np.isfinite(atol):
        raise ValueError("atol must be non-negative and finite")
    if max_it <= 0:
        raise ValueError("max_it must be positive")
    if restart <= 0:
        raise ValueError("restart must be positive")
    if schur_drop_rel < 0.0 or not np.isfinite(schur_drop_rel):
        raise ValueError("schur_drop_rel must be non-negative and finite")

    resolved_controls = _nvmath_cudss_controls_from_arguments(
        dtype=dtype,
        device_ids=device_ids,
        ir_steps=ir_steps,
        use_matching=use_matching,
        reordering_alg=reordering_alg,
        matching_alg=matching_alg,
        factorization_alg=factorization_alg,
        solve_alg=solve_alg,
        pivot_type=pivot_type,
        pivot_threshold=pivot_threshold,
        pivot_epsilon=pivot_epsilon,
        pivot_epsilon_alg=pivot_epsilon_alg,
        nd_nlevels=nd_nlevels,
        host_nthreads=host_nthreads,
        threading_lib=threading_lib,
        hybrid_mode=hybrid_mode,
        hybrid_device_memory_limit=hybrid_device_memory_limit,
        hybrid_execute_mode=hybrid_execute_mode,
        use_cuda_register_memory=use_cuda_register_memory,
        use_superpanels=use_superpanels,
        deterministic_mode=deterministic_mode,
        check_residual=False,
        residual_rtol=residual_rtol,
        controls=controls,
    )
    return cls(
        linear_backend="petsc",
        solver_preset="iterative_schurdiag_cudss_experimental",
        linear_system_dtype=cast(LinearSystemDType, resolved_controls["dtype"]),
        nvmath_cudss_controls=resolved_controls,
        iterative_solver_controls={
            "gmres_rtol": float(rtol),
            "gmres_atol": float(atol),
            "gmres_maxiter": int(max_it),
            "gmres_restart": int(restart),
            "velocity_solver": velocity_solver,
            "schurdiag_drop_rel": float(schur_drop_rel),
            "error_if_not_converged": bool(error_if_not_converged),
        },
        petsc_options={},
    )

scipy_direct classmethod

scipy_direct(
    *,
    linear_system_dtype="float64",
    permc_spec=None,
    diag_pivot_thresh=None,
    relax=None,
    panel_size=None,
    equil=None,
    controls=None,
)

Create options for the serial DOLFINx-assembly/SciPy SuperLU backend.

Source code in src/voids/fem/singlephase/_common.py

@classmethod
def scipy_direct(
    cls,
    *,
    linear_system_dtype: LinearSystemDType = "float64",
    permc_spec: str | None = None,
    diag_pivot_thresh: float | None = None,
    relax: int | None = None,
    panel_size: int | None = None,
    equil: bool | None = None,
    controls: Mapping[str, Any] | None = None,
) -> FEniCSSolverOptions:
    """Create options for the serial DOLFINx-assembly/SciPy SuperLU backend."""

    return cls(
        linear_backend="scipy",
        solver_preset="direct_reference",
        linear_system_dtype=linear_system_dtype,
        superlu_controls=_superlu_controls_from_arguments(
            permc_spec=permc_spec,
            diag_pivot_thresh=diag_pivot_thresh,
            relax=relax,
            panel_size=panel_size,
            equil=equil,
            controls=controls,
        ),
    )

superlu_direct classmethod

superlu_direct(
    *,
    linear_system_dtype="float64",
    permc_spec=None,
    diag_pivot_thresh=None,
    relax=None,
    panel_size=None,
    equil=None,
    controls=None,
)

Create options for the serial DOLFINx-assembly/SuperLU backend.

Source code in src/voids/fem/singlephase/_common.py

@classmethod
def superlu_direct(
    cls,
    *,
    linear_system_dtype: LinearSystemDType = "float64",
    permc_spec: str | None = None,
    diag_pivot_thresh: float | None = None,
    relax: int | None = None,
    panel_size: int | None = None,
    equil: bool | None = None,
    controls: Mapping[str, Any] | None = None,
) -> FEniCSSolverOptions:
    """Create options for the serial DOLFINx-assembly/SuperLU backend."""

    return cls(
        linear_backend="superlu",
        solver_preset="direct_reference",
        linear_system_dtype=linear_system_dtype,
        superlu_controls=_superlu_controls_from_arguments(
            permc_spec=permc_spec,
            diag_pivot_thresh=diag_pivot_thresh,
            relax=relax,
            panel_size=panel_size,
            equil=equil,
            controls=controls,
        ),
    )

umfpack_direct classmethod

umfpack_direct(
    *,
    linear_system_dtype="float64",
    ordering=None,
    strategy=None,
    pivot_tolerance=None,
    sym_pivot_tolerance=None,
    scale=None,
    block_size=None,
    controls=None,
)

Create options for the serial DOLFINx-assembly/UMFPACK backend.

Source code in src/voids/fem/singlephase/_common.py

@classmethod
def umfpack_direct(
    cls,
    *,
    linear_system_dtype: LinearSystemDType = "float64",
    ordering: str | int | float | None = None,
    strategy: str | int | float | None = None,
    pivot_tolerance: float | None = None,
    sym_pivot_tolerance: float | None = None,
    scale: str | int | float | None = None,
    block_size: int | None = None,
    controls: Mapping[str, Any] | None = None,
) -> FEniCSSolverOptions:
    """Create options for the serial DOLFINx-assembly/UMFPACK backend."""

    umfpack_controls: dict[str, Any] = dict(controls or {})
    if ordering is not None:
        umfpack_controls["ordering"] = ordering
    if strategy is not None:
        umfpack_controls["strategy"] = strategy
    if pivot_tolerance is not None:
        umfpack_controls["pivot_tolerance"] = float(pivot_tolerance)
    if sym_pivot_tolerance is not None:
        umfpack_controls["sym_pivot_tolerance"] = float(sym_pivot_tolerance)
    if scale is not None:
        umfpack_controls["scale"] = scale
    if block_size is not None:
        umfpack_controls["block_size"] = int(block_size)
    return cls(
        linear_backend="umfpack",
        solver_preset="direct_reference",
        linear_system_dtype=linear_system_dtype,
        umfpack_controls=umfpack_controls,
    )

nvmath_cudss_direct classmethod

nvmath_cudss_direct(
    *,
    dtype="float64",
    device_ids=None,
    ir_steps=5,
    use_matching=True,
    reordering_alg=None,
    matching_alg=None,
    factorization_alg=None,
    solve_alg=None,
    pivot_type=None,
    pivot_threshold=None,
    pivot_epsilon=None,
    pivot_epsilon_alg=None,
    nd_nlevels=None,
    host_nthreads=None,
    threading_lib=None,
    hybrid_mode=None,
    hybrid_device_memory_limit=None,
    hybrid_execute_mode=None,
    use_cuda_register_memory=None,
    use_superpanels=None,
    deterministic_mode=None,
    check_residual=True,
    residual_rtol=None,
    controls=None,
)

Create options for the optional CUDA cuDSS direct-solver backend.

The backend uses serial DOLFINx assembly, converts the assembled system to SciPy CSR format, copies the sparse matrix and right-hand side to CUDA tensors, and solves the system through nvmath.bindings.cudss. It requires PyTorch with CUDA support and the nvmath/cuDSS runtime at solve time. These controls affect only the numerical linear solve; they do not change the Brinkman/Darcy weak form, porosity, permeability, viscosity, pressure drop, or dimensional scaling.

Parameters:

Name	Type	Description	Default
dtype	Literal['float32', 'float64']	Floating-point value precision used for the cuDSS matrix, right-hand side, and solution. "float64" is the conservative default. "float32" can reduce memory use and wall time, but it must be accepted only after the residual check and a same-map reference comparison because high permeability contrasts can be ill-conditioned.	'float64'
device_ids	int | Sequence[int] | Literal['all'] | None	CUDA device selection. None uses PyTorch's current CUDA device, an integer selects one device, a sequence such as (0, 1) requests a single-node multi-GPU cuDSS handle, and "all" uses all CUDA devices visible to PyTorch. Some cuDSS controls are not supported by all single-GPU and multi-GPU execution paths.	None
ir_steps	int	cuDSS iterative-refinement step count (ConfigParam.IR_N_STEPS). voids sets this to 5 by default, although a fresh cuDSS config in the tested stack reports IR_N_STEPS = 0. Increasing this value can improve single-precision residuals, but the effect is problem-dependent and not guaranteed to be monotonic.	5
use_matching	bool	Whether to request cuDSS matching/scaling (ConfigParam.USE_MATCHING). Matching is usually helpful for general sparse matrices because it can reduce pivot perturbations; it is enabled by default.	True
reordering_alg	str | int | None	Optional cuDSS algorithm selectors for REORDERING_ALG, MATCHING_ALG, FACTORIZATION_ALG, SOLVE_ALG, and PIVOT_EPSILON_ALG. Values may be integers 0 through 5 or strings such as "default", "alg_1", or "3". Their exact meaning and support are defined by the installed cuDSS version; unsupported combinations raise a cuDSS error during analysis, factorization, or solve.	None
matching_alg	str | int | None	Optional cuDSS algorithm selectors for REORDERING_ALG, MATCHING_ALG, FACTORIZATION_ALG, SOLVE_ALG, and PIVOT_EPSILON_ALG. Values may be integers 0 through 5 or strings such as "default", "alg_1", or "3". Their exact meaning and support are defined by the installed cuDSS version; unsupported combinations raise a cuDSS error during analysis, factorization, or solve.	None
factorization_alg	str | int | None	Optional cuDSS algorithm selectors for REORDERING_ALG, MATCHING_ALG, FACTORIZATION_ALG, SOLVE_ALG, and PIVOT_EPSILON_ALG. Values may be integers 0 through 5 or strings such as "default", "alg_1", or "3". Their exact meaning and support are defined by the installed cuDSS version; unsupported combinations raise a cuDSS error during analysis, factorization, or solve.	None
solve_alg	str | int | None	Optional cuDSS algorithm selectors for REORDERING_ALG, MATCHING_ALG, FACTORIZATION_ALG, SOLVE_ALG, and PIVOT_EPSILON_ALG. Values may be integers 0 through 5 or strings such as "default", "alg_1", or "3". Their exact meaning and support are defined by the installed cuDSS version; unsupported combinations raise a cuDSS error during analysis, factorization, or solve.	None
pivot_epsilon_alg	str | int | None	Optional cuDSS algorithm selectors for REORDERING_ALG, MATCHING_ALG, FACTORIZATION_ALG, SOLVE_ALG, and PIVOT_EPSILON_ALG. Values may be integers 0 through 5 or strings such as "default", "alg_1", or "3". Their exact meaning and support are defined by the installed cuDSS version; unsupported combinations raise a cuDSS error during analysis, factorization, or solve.	None
pivot_type	Literal['col', 'row', 'none'] | None	Optional cuDSS pivoting mode (ConfigParam.PIVOT_TYPE): "col" for column pivoting, "row" for row pivoting, or "none" to request no pivoting. Leave as None to use the cuDSS default.	None
pivot_threshold	float | None	Optional non-negative pivoting threshold (ConfigParam.PIVOT_THRESHOLD). Larger values can make pivoting more conservative, but may increase cost or be ignored depending on the selected cuDSS algorithm.	None
pivot_epsilon	float | None	Optional non-negative pivot perturbation/floor (ConfigParam.PIVOT_EPSILON). This can stabilize factorizations with very small pivots, especially in single precision. Treat this as a solver-stabilization experiment: high-contrast FEM and TPFA maps can pass a loose residual check while still failing same-map field comparisons in single precision.	None
nd_nlevels	int | None	Optional non-negative nested-dissection level control (ConfigParam.ND_NLEVELS) used by cuDSS reordering algorithms that support it.	None
host_nthreads	int | None	Optional positive host-thread count for cuDSS host-side work (ConfigParam.HOST_NTHREADS). This affects execution only when a cuDSS threading-layer library is loaded.	None
threading_lib	str | None	Optional path to a cuDSS threading-layer library. Use "auto" to request the packaged libcudss_mtlayer_gomp library when available. If host threading is requested and this is omitted, voids uses CUDSS_THREADING_LIB when set and otherwise auto-loads the packaged threading layer when host_nthreads or hybrid_execute_mode=True is requested.	None
hybrid_mode	bool | None	Optional request to enable cuDSS hybrid memory mode (ConfigParam.HYBRID_MODE), which can spill factorization data to host memory when device memory is limiting.	None
hybrid_device_memory_limit	int | None	Optional positive device-memory limit in bytes for cuDSS hybrid memory mode (ConfigParam.HYBRID_DEVICE_MEMORY_LIMIT). voids applies this after analysis and before factorization. Multi-GPU hybrid-memory limit handling is runtime dependent in the low-level nvmath binding; compare against a small same-configuration probe before relying on it.	None
hybrid_execute_mode	bool | None	Optional request to enable cuDSS hybrid execute mode (ConfigParam.HYBRID_EXECUTE_MODE). Support is runtime/backend dependent.	None
use_cuda_register_memory	bool | None	Optional request to register host memory with CUDA (ConfigParam.USE_CUDA_REGISTER_MEMORY) for hybrid-memory runs.	None
use_superpanels	bool | None	Optional flag for cuDSS superpanel optimization (ConfigParam.USE_SUPERPANELS). Leave as None to use the cuDSS default.	None
deterministic_mode	bool | None	Optional request for deterministic cuDSS execution (ConfigParam.DETERMINISTIC_MODE). This is version/backend dependent and may raise NOT_SUPPORTED in some cuDSS paths.	None
check_residual	bool	If true, voids computes the assembled-system relative residual after cuDSS returns and raises an error when it exceeds residual_rtol. Keep this enabled for lower-precision studies.	True
residual_rtol	float | None	Relative residual acceptance tolerance used when check_residual=True. If omitted, voids uses 1.0e-8 for float64 and 1.0e-4 for float32.	None
controls	Mapping[str, Any] | None	Optional low-level control mapping. Keys use the same names as the keyword arguments, with hyphens normalized to underscores. Explicit keyword arguments override the same keys in this mapping.	None

Returns:

Type	Description
FEniCSSolverOptions	Solver options with linear_backend="nvmath_cudss" and resolved cuDSS controls stored in nvmath_cudss_controls.

Notes

This backend is experimental and hardware/runtime dependent. It should be compared against a same-map direct reference before using the result for scientific claims.

Source code in src/voids/fem/singlephase/_common.py

@classmethod
def nvmath_cudss_direct(
    cls,
    *,
    dtype: Literal["float32", "float64"] = "float64",
    device_ids: int | Sequence[int] | Literal["all"] | None = None,
    ir_steps: int = 5,
    use_matching: bool = True,
    reordering_alg: str | int | None = None,
    matching_alg: str | int | None = None,
    factorization_alg: str | int | None = None,
    solve_alg: str | int | None = None,
    pivot_type: Literal["col", "row", "none"] | None = None,
    pivot_threshold: float | None = None,
    pivot_epsilon: float | None = None,
    pivot_epsilon_alg: str | int | None = None,
    nd_nlevels: int | None = None,
    host_nthreads: int | None = None,
    threading_lib: str | None = None,
    hybrid_mode: bool | None = None,
    hybrid_device_memory_limit: int | None = None,
    hybrid_execute_mode: bool | None = None,
    use_cuda_register_memory: bool | None = None,
    use_superpanels: bool | None = None,
    deterministic_mode: bool | None = None,
    check_residual: bool = True,
    residual_rtol: float | None = None,
    controls: Mapping[str, Any] | None = None,
) -> FEniCSSolverOptions:
    """Create options for the optional CUDA cuDSS direct-solver backend.

    The backend uses serial DOLFINx assembly, converts the assembled system
    to SciPy CSR format, copies the sparse matrix and right-hand side to
    CUDA tensors, and solves the system through ``nvmath.bindings.cudss``.
    It requires PyTorch with CUDA support and the nvmath/cuDSS runtime at
    solve time. These controls affect only the numerical linear solve; they
    do not change the Brinkman/Darcy weak form, porosity, permeability,
    viscosity, pressure drop, or dimensional scaling.

    Parameters
    ----------
    dtype :
        Floating-point value precision used for the cuDSS matrix, right-hand
        side, and solution. ``"float64"`` is the conservative default.
        ``"float32"`` can reduce memory use and wall time, but it must be
        accepted only after the residual check and a same-map reference
        comparison because high permeability contrasts can be ill-conditioned.
    device_ids :
        CUDA device selection. ``None`` uses PyTorch's current CUDA device,
        an integer selects one device, a sequence such as ``(0, 1)`` requests
        a single-node multi-GPU cuDSS handle, and ``"all"`` uses all CUDA
        devices visible to PyTorch. Some cuDSS controls are not supported by
        all single-GPU and multi-GPU execution paths.
    ir_steps :
        cuDSS iterative-refinement step count
        (``ConfigParam.IR_N_STEPS``). ``voids`` sets this to ``5`` by
        default, although a fresh cuDSS config in the tested stack reports
        ``IR_N_STEPS = 0``. Increasing this value can improve single-precision
        residuals, but the effect is problem-dependent and not guaranteed to
        be monotonic.
    use_matching :
        Whether to request cuDSS matching/scaling
        (``ConfigParam.USE_MATCHING``). Matching is usually helpful for
        general sparse matrices because it can reduce pivot perturbations;
        it is enabled by default.
    reordering_alg, matching_alg, factorization_alg, solve_alg, pivot_epsilon_alg :
        Optional cuDSS algorithm selectors for
        ``REORDERING_ALG``, ``MATCHING_ALG``, ``FACTORIZATION_ALG``,
        ``SOLVE_ALG``, and ``PIVOT_EPSILON_ALG``. Values may be integers
        ``0`` through ``5`` or strings such as ``"default"``, ``"alg_1"``,
        or ``"3"``. Their exact meaning and support are defined by the
        installed cuDSS version; unsupported combinations raise a cuDSS
        error during analysis, factorization, or solve.
    pivot_type :
        Optional cuDSS pivoting mode (``ConfigParam.PIVOT_TYPE``): ``"col"``
        for column pivoting, ``"row"`` for row pivoting, or ``"none"`` to
        request no pivoting. Leave as ``None`` to use the cuDSS default.
    pivot_threshold :
        Optional non-negative pivoting threshold
        (``ConfigParam.PIVOT_THRESHOLD``). Larger values can make pivoting
        more conservative, but may increase cost or be ignored depending on
        the selected cuDSS algorithm.
    pivot_epsilon :
        Optional non-negative pivot perturbation/floor
        (``ConfigParam.PIVOT_EPSILON``). This can stabilize factorizations
        with very small pivots, especially in single precision. Treat this
        as a solver-stabilization experiment: high-contrast FEM and TPFA
        maps can pass a loose residual check while still failing same-map
        field comparisons in single precision.
    nd_nlevels :
        Optional non-negative nested-dissection level control
        (``ConfigParam.ND_NLEVELS``) used by cuDSS reordering algorithms that
        support it.
    host_nthreads :
        Optional positive host-thread count for cuDSS host-side work
        (``ConfigParam.HOST_NTHREADS``). This affects execution only when a
        cuDSS threading-layer library is loaded.
    threading_lib :
        Optional path to a cuDSS threading-layer library. Use ``"auto"`` to
        request the packaged ``libcudss_mtlayer_gomp`` library when
        available. If host threading is requested and this is omitted,
        ``voids`` uses ``CUDSS_THREADING_LIB`` when set and otherwise
        auto-loads the packaged threading layer when ``host_nthreads`` or
        ``hybrid_execute_mode=True`` is requested.
    hybrid_mode :
        Optional request to enable cuDSS hybrid memory mode
        (``ConfigParam.HYBRID_MODE``), which can spill factorization data to
        host memory when device memory is limiting.
    hybrid_device_memory_limit :
        Optional positive device-memory limit in bytes for cuDSS hybrid
        memory mode (``ConfigParam.HYBRID_DEVICE_MEMORY_LIMIT``). ``voids``
        applies this after analysis and before factorization. Multi-GPU
        hybrid-memory limit handling is runtime dependent in the low-level
        nvmath binding; compare against a small same-configuration probe
        before relying on it.
    hybrid_execute_mode :
        Optional request to enable cuDSS hybrid execute mode
        (``ConfigParam.HYBRID_EXECUTE_MODE``). Support is runtime/backend
        dependent.
    use_cuda_register_memory :
        Optional request to register host memory with CUDA
        (``ConfigParam.USE_CUDA_REGISTER_MEMORY``) for hybrid-memory runs.
    use_superpanels :
        Optional flag for cuDSS superpanel optimization
        (``ConfigParam.USE_SUPERPANELS``). Leave as ``None`` to use the cuDSS
        default.
    deterministic_mode :
        Optional request for deterministic cuDSS execution
        (``ConfigParam.DETERMINISTIC_MODE``). This is version/backend
        dependent and may raise ``NOT_SUPPORTED`` in some cuDSS paths.
    check_residual :
        If true, ``voids`` computes the assembled-system relative residual
        after cuDSS returns and raises an error when it exceeds
        ``residual_rtol``. Keep this enabled for lower-precision studies.
    residual_rtol :
        Relative residual acceptance tolerance used when
        ``check_residual=True``. If omitted, ``voids`` uses ``1.0e-8`` for
        ``float64`` and ``1.0e-4`` for ``float32``.
    controls :
        Optional low-level control mapping. Keys use the same names as the
        keyword arguments, with hyphens normalized to underscores. Explicit
        keyword arguments override the same keys in this mapping.

    Returns
    -------
    FEniCSSolverOptions
        Solver options with ``linear_backend="nvmath_cudss"`` and resolved
        cuDSS controls stored in ``nvmath_cudss_controls``.

    Notes
    -----
    This backend is experimental and hardware/runtime dependent. It should
    be compared against a same-map direct reference before using the result
    for scientific claims.
    """

    resolved_controls = _nvmath_cudss_controls_from_arguments(
        dtype=dtype,
        device_ids=device_ids,
        ir_steps=ir_steps,
        use_matching=use_matching,
        reordering_alg=reordering_alg,
        matching_alg=matching_alg,
        factorization_alg=factorization_alg,
        solve_alg=solve_alg,
        pivot_type=pivot_type,
        pivot_threshold=pivot_threshold,
        pivot_epsilon=pivot_epsilon,
        pivot_epsilon_alg=pivot_epsilon_alg,
        nd_nlevels=nd_nlevels,
        host_nthreads=host_nthreads,
        threading_lib=threading_lib,
        hybrid_mode=hybrid_mode,
        hybrid_device_memory_limit=hybrid_device_memory_limit,
        hybrid_execute_mode=hybrid_execute_mode,
        use_cuda_register_memory=use_cuda_register_memory,
        use_superpanels=use_superpanels,
        deterministic_mode=deterministic_mode,
        check_residual=check_residual,
        residual_rtol=residual_rtol,
        controls=controls,
    )
    return cls(
        linear_backend="nvmath_cudss",
        solver_preset="direct_reference",
        linear_system_dtype=cast(LinearSystemDType, resolved_controls["dtype"]),
        nvmath_cudss_controls=resolved_controls,
    )

pardiso_direct classmethod

pardiso_direct(*, linear_system_dtype='float64')

Create options for the serial DOLFINx-assembly/PARDISO backend.

Source code in src/voids/fem/singlephase/_common.py

@classmethod
def pardiso_direct(
    cls,
    *,
    linear_system_dtype: LinearSystemDType = "float64",
) -> FEniCSSolverOptions:
    """Create options for the serial DOLFINx-assembly/PARDISO backend."""

    return cls(
        linear_backend="pardiso",
        solver_preset="direct_reference",
        linear_system_dtype=linear_system_dtype,
    )

FEMUpscalingResult dataclass

Principal-direction FEM micro-continuum permeability result.

Source code in src/voids/fem/singlephase/upscaling.py

@dataclass(slots=True)
class FEMUpscalingResult:
    """Principal-direction FEM micro-continuum permeability result."""

    results: dict[str, FEMSinglePhaseResult]
    backend: str

    @property
    def permeability(self) -> dict[str, float]:
        """Return effective permeability by principal axis."""

        return {axis: result.permeability for axis, result in self.results.items()}

    @property
    def solve_seconds(self) -> dict[str, float]:
        """Return wall-clock solve time by principal axis."""

        return {axis: result.solve_seconds for axis, result in self.results.items()}

permeability property

permeability

Return effective permeability by principal axis.

solve_seconds property

solve_seconds

Return wall-clock solve time by principal axis.

solve_brinkman_taylor_hood

solve_brinkman_taylor_hood(
    problem,
    *,
    flow_axis="x",
    pressure_inlet=1.0,
    pressure_outlet=0.0,
    options=None,
    nondimensional=False,
)

Solve the Darcy-Brinkman micro-continuum model with Taylor-Hood elements.

The weak form uses CG2 velocity and CG1 pressure:

(mu / phi grad u, grad v) + (mu / K u, v) - (p, div v) + (q, div u) = boundary pressure work.

K and phi are piecewise-constant maps supplied through :class:~voids.fem.singlephase.FEMMapProblem.

Source code in src/voids/fem/singlephase/taylorhood.py

def solve_brinkman_taylor_hood(
    problem: FEMMapProblem,
    *,
    flow_axis: str = "x",
    pressure_inlet: float = 1.0,
    pressure_outlet: float = 0.0,
    options: FEniCSSolverOptions | None = None,
    nondimensional: bool | BrinkmanNondimensionalization = False,
) -> FEMSinglePhaseResult:
    """Solve the Darcy-Brinkman micro-continuum model with Taylor-Hood elements.

    The weak form uses CG2 velocity and CG1 pressure:

    ``(mu / phi grad u, grad v) + (mu / K u, v)
    - (p, div v) + (q, div u) = boundary pressure work``.

    ``K`` and ``phi`` are piecewise-constant maps supplied through
    :class:`~voids.fem.singlephase.FEMMapProblem`.
    """

    nondimensional_options = _resolve_brinkman_nondimensionalization(nondimensional)
    scales = None
    if nondimensional_options is not None:
        _validate_pressure_drop(pressure_inlet, pressure_outlet)
        context_for_scales = _build_context(
            problem,
            flow_axis=flow_axis,
            api=_require_dolfinx_core(),
        )
        scales = _brinkman_nondimensional_scales(
            context_for_scales,
            problem,
            pressure_inlet=pressure_inlet,
            pressure_outlet=pressure_outlet,
            velocity_scale=nondimensional_options.velocity_scale,
        )

    def form_builder(context, u, p, v, q):
        ufl = context.api.ufl
        if scales is None:
            gamma = context.coefficients["gamma"]
            nu_eff = context.coefficients["nu_eff"]
        else:
            gamma, nu_eff = _brinkman_nondimensional_coefficients(context, problem, scales)
        return (
            nu_eff * ufl.inner(ufl.grad(u), ufl.grad(v)) * context.dx
            + ufl.inner(gamma * u, v) * context.dx
            - p * ufl.div(v) * context.dx
            + q * ufl.div(u) * context.dx
        )

    return _solve_with_form_builder(
        problem,
        flow_axis=flow_axis,
        pressure_inlet=pressure_inlet,
        pressure_outlet=pressure_outlet,
        options=options,
        velocity_degree=2,
        pressure_family="Lagrange",
        method="Darcy-Brinkman Taylor-Hood CG2 x CG1",
        formulation="brinkman_taylor_hood_p2p1",
        prefix_suffix=f"brinkman_taylor_hood_{flow_axis}",
        form_builder=form_builder,
        boundary_pressure_inlet=None if scales is None else 1.0,
        boundary_pressure_outlet=None if scales is None else 0.0,
        velocity_scale=1.0 if scales is None else scales.velocity_scale,
        pressure_scale=1.0 if scales is None else scales.pressure_scale,
        extra_metadata=_nondimensional_metadata(scales),
    )

solve_darcy_taylor_hood

solve_darcy_taylor_hood(
    problem,
    *,
    flow_axis="x",
    pressure_inlet=1.0,
    pressure_outlet=0.0,
    options=None,
)

Solve the mixed Darcy-Darcy comparison model with Taylor-Hood elements.

The weak form uses a CG2 velocity approximation and a CG1 pressure approximation:

(mu / K u, v) - (p, div v) + (q, div u) = boundary pressure work.

Side boundaries impose only zero normal velocity on faces transverse to the flow axis. Inlet and outlet pressures are applied as natural traction terms.

Source code in src/voids/fem/singlephase/taylorhood.py

def solve_darcy_taylor_hood(
    problem: FEMMapProblem,
    *,
    flow_axis: str = "x",
    pressure_inlet: float = 1.0,
    pressure_outlet: float = 0.0,
    options: FEniCSSolverOptions | None = None,
) -> FEMSinglePhaseResult:
    """Solve the mixed Darcy-Darcy comparison model with Taylor-Hood elements.

    The weak form uses a CG2 velocity approximation and a CG1 pressure
    approximation:

    ``(mu / K u, v) - (p, div v) + (q, div u) = boundary pressure work``.

    Side boundaries impose only zero normal velocity on faces transverse to the
    flow axis. Inlet and outlet pressures are applied as natural traction terms.
    """

    def form_builder(context, u, p, v, q):
        ufl = context.api.ufl
        gamma = context.coefficients["gamma"]
        return (
            ufl.inner(gamma * u, v) * context.dx
            - p * ufl.div(v) * context.dx
            + q * ufl.div(u) * context.dx
        )

    return _solve_with_form_builder(
        problem,
        flow_axis=flow_axis,
        pressure_inlet=pressure_inlet,
        pressure_outlet=pressure_outlet,
        options=options,
        velocity_degree=2,
        pressure_family="Lagrange",
        method="Darcy-Darcy Taylor-Hood CG2 x CG1",
        formulation="darcy_taylor_hood_p2p1",
        prefix_suffix=f"darcy_taylor_hood_{flow_axis}",
        form_builder=form_builder,
    )

upscale_permeability_fem

upscale_permeability_fem(
    problem,
    *,
    backend="taylor_hood_brinkman",
    axes=None,
    pressure_inlet=1.0,
    pressure_outlet=0.0,
    options=None,
    backend_kwargs=None,
)

Compute principal-direction FEM permeability estimates.

Parameters:

Name	Type	Description	Default
problem	FEMMapProblem	Porosity/permeability map problem.	required
backend	FEMBackend	Either a backend name or a compatible solver callable. Supported names are "taylor_hood_brinkman", "taylor_hood_darcy", and "usfem_brinkman".	'taylor_hood_brinkman'
axes	tuple[str, ...] | None	Principal axes to solve. By default all axes supported by the map dimensionality are solved.	None
pressure_inlet	float	Natural pressure values imposed on opposite faces of each flow axis.	1.0
pressure_outlet	float	Natural pressure values imposed on opposite faces of each flow axis.	1.0
options	FEniCSSolverOptions | None	Linear solver options passed to the FEM backend. The default preserves PETSc where available and uses a serial SciPy direct solve on native Windows when PETSc is unavailable.	None
backend_kwargs	dict[str, object] | None	Additional backend-specific keyword arguments.	None

Source code in src/voids/fem/singlephase/upscaling.py

def upscale_permeability_fem(
    problem: FEMMapProblem,
    *,
    backend: FEMBackend = "taylor_hood_brinkman",
    axes: tuple[str, ...] | None = None,
    pressure_inlet: float = 1.0,
    pressure_outlet: float = 0.0,
    options: FEniCSSolverOptions | None = None,
    backend_kwargs: dict[str, object] | None = None,
) -> FEMUpscalingResult:
    """Compute principal-direction FEM permeability estimates.

    Parameters
    ----------
    problem :
        Porosity/permeability map problem.
    backend :
        Either a backend name or a compatible solver callable. Supported names
        are ``"taylor_hood_brinkman"``, ``"taylor_hood_darcy"``, and
        ``"usfem_brinkman"``.
    axes :
        Principal axes to solve. By default all axes supported by the map
        dimensionality are solved.
    pressure_inlet, pressure_outlet :
        Natural pressure values imposed on opposite faces of each flow axis.
    options :
        Linear solver options passed to the FEM backend. The default preserves
        PETSc where available and uses a serial SciPy direct solve on native
        Windows when PETSc is unavailable.
    backend_kwargs :
        Additional backend-specific keyword arguments.
    """

    solve_axes = axes or _default_axes(problem.permeability_map.ndim)
    solver = _backend_from_name(backend) if isinstance(backend, str) else backend
    kwargs = dict(backend_kwargs or {})
    results = {
        axis: solver(
            problem,
            flow_axis=axis,
            pressure_inlet=pressure_inlet,
            pressure_outlet=pressure_outlet,
            options=options,
            **kwargs,
        )
        for axis in solve_axes
    }
    backend_name = backend if isinstance(backend, str) else getattr(backend, "__name__", "callable")
    return FEMUpscalingResult(results=results, backend=backend_name)

upscale_principal_permeabilities_fem

upscale_principal_permeabilities_fem(
    problem,
    *,
    backend="taylor_hood_brinkman",
    axes=None,
    pressure_inlet=1.0,
    pressure_outlet=0.0,
    options=None,
    backend_kwargs=None,
)

Return only the principal FEM permeability values.

Source code in src/voids/fem/singlephase/upscaling.py

def upscale_principal_permeabilities_fem(
    problem: FEMMapProblem,
    *,
    backend: FEMBackend = "taylor_hood_brinkman",
    axes: tuple[str, ...] | None = None,
    pressure_inlet: float = 1.0,
    pressure_outlet: float = 0.0,
    options: FEniCSSolverOptions | None = None,
    backend_kwargs: dict[str, object] | None = None,
) -> dict[str, float]:
    """Return only the principal FEM permeability values."""

    return upscale_permeability_fem(
        problem,
        backend=backend,
        axes=axes,
        pressure_inlet=pressure_inlet,
        pressure_outlet=pressure_outlet,
        options=options,
        backend_kwargs=backend_kwargs,
    ).permeability

solve_brinkman_usfem

solve_brinkman_usfem(
    problem,
    *,
    flow_axis="x",
    pressure_inlet=1.0,
    pressure_outlet=0.0,
    tau_factor=1.0,
    m_t=1.0 / 3.0,
    alpha_edge=1.0,
    facet_law="reaction_diffusion",
    facet_size_mode="cell_diameter",
    tau_gamma_cap=None,
    pressure_degree=1,
    options=None,
    nondimensional=False,
)

Solve a stabilized Darcy-Brinkman micro-continuum model.

The formulation uses CG1 velocity and discontinuous pressure fields. DG1 is the backward-compatible default; pass pressure_degree=0 for the low-order CG1 x DG0 pair used in manufactured-solution and vug studies. It augments the Brinkman weak form with a residual-based cell stabilization term and an interior pressure-jump penalty. The coefficients are intended for porosity/permeability maps obtained from a segmented image.

Set tau_factor=0 to disable the cell residual while retaining the pressure-jump term. In high-contrast CG1 x DG0 studies, tau_gamma_cap can enforce gamma * tau_K <= tau_gamma_cap and avoid near-cancellation of the physical drag. facet_size_mode="facet_measure" uses edge length in 2D and the square root of facet area in 3D; the latter is a measure-based length, not an exact triangular-facet diameter.

Source code in src/voids/fem/singlephase/usfem.py

def solve_brinkman_usfem(
    problem: FEMMapProblem,
    *,
    flow_axis: str = "x",
    pressure_inlet: float = 1.0,
    pressure_outlet: float = 0.0,
    tau_factor: float = 1.0,
    m_t: float = 1.0 / 3.0,
    alpha_edge: float = 1.0,
    facet_law: USFEMFacetLaw = "reaction_diffusion",
    facet_size_mode: USFEMFacetSizeMode = "cell_diameter",
    tau_gamma_cap: float | None = None,
    pressure_degree: Literal[0, 1] = 1,
    options: FEniCSSolverOptions | None = None,
    nondimensional: bool | BrinkmanNondimensionalization = False,
) -> FEMSinglePhaseResult:
    """Solve a stabilized Darcy-Brinkman micro-continuum model.

    The formulation uses CG1 velocity and discontinuous pressure fields. DG1 is
    the backward-compatible default; pass ``pressure_degree=0`` for the
    low-order CG1 x DG0 pair used in manufactured-solution and vug studies. It
    augments the Brinkman weak form with a residual-based cell stabilization term and an
    interior pressure-jump penalty. The coefficients are intended for
    porosity/permeability maps obtained from a segmented image.

    Set ``tau_factor=0`` to disable the cell residual while retaining the
    pressure-jump term. In high-contrast CG1 x DG0 studies,
    ``tau_gamma_cap`` can enforce ``gamma * tau_K <= tau_gamma_cap`` and avoid
    near-cancellation of the physical drag. ``facet_size_mode="facet_measure"``
    uses edge length in 2D and the square root of facet area in 3D; the latter
    is a measure-based length, not an exact triangular-facet diameter.
    """

    _validate_usfem_controls(
        tau_factor=tau_factor,
        m_t=m_t,
        alpha_edge=alpha_edge,
        facet_law=facet_law,
        facet_size_mode=facet_size_mode,
        tau_gamma_cap=tau_gamma_cap,
    )
    if pressure_degree not in {0, 1}:
        raise ValueError("pressure_degree must be either 0 or 1")
    uncapped_max_tau_gamma = (
        _p1dg0_max_uncapped_tau_gamma(
            problem,
            tau_factor=tau_factor,
            m_t=m_t,
        )
        if pressure_degree == 0
        else None
    )
    if (
        uncapped_max_tau_gamma is not None
        and tau_factor > 0.0
        and tau_gamma_cap is None
        and uncapped_max_tau_gamma >= 0.9
    ):
        warnings.warn(
            "The uncapped CG1 x DG0 cell term has estimated "
            f"max(gamma * tau_K)={uncapped_max_tau_gamma:.3g}; it can nearly "
            "cancel physical drag. Set tau_gamma_cap below 1 or use "
            "tau_factor=0 as an explicit sensitivity branch.",
            RuntimeWarning,
            stacklevel=2,
        )
    nondimensional_options = _resolve_brinkman_nondimensionalization(nondimensional)
    scales = None
    if nondimensional_options is not None:
        _validate_pressure_drop(pressure_inlet, pressure_outlet)
        context_for_scales = _build_context(
            problem,
            flow_axis=flow_axis,
            api=_require_dolfinx_core(),
        )
        scales = _brinkman_nondimensional_scales(
            context_for_scales,
            problem,
            pressure_inlet=pressure_inlet,
            pressure_outlet=pressure_outlet,
            velocity_scale=nondimensional_options.velocity_scale,
        )

    def form_builder(
        context: _FEMContext,
        u: _UFLExpression,
        p: _UFLExpression,
        v: _UFLExpression,
        q: _UFLExpression,
    ) -> _UFLExpression:
        ufl = cast(_UFLAlgebra, context.api.ufl)
        dx = cast(_UFLExpression, context.dx)
        dS = cast(_UFLExpression, context.dS)
        gamma_override: _UFLExpression | None = None
        nu_eff_override: _UFLExpression | None = None
        if scales is not None:
            gamma_raw, nu_eff_raw = _brinkman_nondimensional_coefficients(
                context,
                problem,
                scales,
            )
            gamma_override = cast(_UFLExpression, gamma_raw)
            nu_eff_override = cast(_UFLExpression, nu_eff_raw)
        gamma, nu_eff, tau, tau_f = _usfem_stabilization_terms(
            context,
            tau_factor=tau_factor,
            m_t=m_t,
            alpha_edge=alpha_edge,
            facet_law=facet_law,
            gamma=gamma_override,
            nu_eff=nu_eff_override,
            facet_size=_facet_size_expression(context, facet_size_mode),
            tau_gamma_cap=tau_gamma_cap,
        )
        residual_u = gamma * u + ufl.grad(p) - nu_eff * ufl.div(ufl.grad(u))
        residual_vq = gamma * v - ufl.grad(q) - nu_eff * ufl.div(ufl.grad(v))
        return (
            nu_eff * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
            + ufl.inner(gamma * u, v) * dx
            - p * ufl.div(v) * dx
            + q * ufl.div(u) * dx
            + tau_f * ufl.jump(p) * ufl.jump(q) * dS
            - tau * ufl.inner(residual_u, residual_vq) * dx
        )

    result = _solve_with_form_builder(
        problem,
        flow_axis=flow_axis,
        pressure_inlet=pressure_inlet,
        pressure_outlet=pressure_outlet,
        options=options,
        velocity_degree=1,
        pressure_family="DG",
        pressure_degree=pressure_degree,
        method=f"Darcy-Brinkman USFEM CG1 x DG{pressure_degree}",
        formulation=f"brinkman_usfem_p1dg{pressure_degree}",
        prefix_suffix=f"brinkman_usfem_{flow_axis}",
        form_builder=form_builder,
        boundary_pressure_inlet=None if scales is None else 1.0,
        boundary_pressure_outlet=None if scales is None else 0.0,
        velocity_scale=1.0 if scales is None else scales.velocity_scale,
        pressure_scale=1.0 if scales is None else scales.pressure_scale,
        extra_metadata=_nondimensional_metadata(scales),
    )
    result.metadata.update(
        {
            "tau_factor": float(tau_factor),
            "m_t": float(m_t),
            "alpha_edge": float(alpha_edge),
            "alpha_edge_active": facet_law in {"classic", "shifted"},
            "facet_law": facet_law,
            "facet_size_mode": facet_size_mode,
            "tau_gamma_cap": (None if tau_gamma_cap is None else float(tau_gamma_cap)),
            "p1dg0_uncapped_max_tau_gamma": uncapped_max_tau_gamma,
        }
    )
    return result

solve_brinkman_usfem_block

solve_brinkman_usfem_block(
    problem,
    *,
    flow_axis="x",
    pressure_inlet=1.0,
    pressure_outlet=0.0,
    tau_factor=1.0,
    m_t=1.0 / 3.0,
    alpha_edge=1.0,
    facet_law="reaction_diffusion",
    facet_size_mode="cell_diameter",
    tau_gamma_cap=None,
    options=None,
    matrix_kind="mpi",
    preconditioner="none",
    nondimensional=False,
)

Solve USFEM through explicit velocity/pressure block forms.

This experimental path expands the same CG1 velocity/DG1 pressure USFEM bilinear form into separate (u, p) blocks. matrix_kind="mpi" assembles the block forms into an ordinary monolithic PETSc matrix and is useful for parity checks against :func:solve_brinkman_usfem. matrix_kind="nest" preserves the block structure for PETSc field-split experiments. Pass preconditioner="diagonal" only for iterative field-split experiments, not for monolithic direct solves. Results from the nested mode should be treated as experimental until compared against a direct reference on the same coefficient map.

Source code in src/voids/fem/singlephase/usfem.py

def solve_brinkman_usfem_block(
    problem: FEMMapProblem,
    *,
    flow_axis: str = "x",
    pressure_inlet: float = 1.0,
    pressure_outlet: float = 0.0,
    tau_factor: float = 1.0,
    m_t: float = 1.0 / 3.0,
    alpha_edge: float = 1.0,
    facet_law: USFEMFacetLaw = "reaction_diffusion",
    facet_size_mode: USFEMFacetSizeMode = "cell_diameter",
    tau_gamma_cap: float | None = None,
    options: FEniCSSolverOptions | None = None,
    matrix_kind: Literal["mpi", "nest"] = "mpi",
    preconditioner: Literal["none", "diagonal"] = "none",
    nondimensional: bool | BrinkmanNondimensionalization = False,
) -> FEMSinglePhaseResult:
    """Solve USFEM through explicit velocity/pressure block forms.

    This experimental path expands the same CG1 velocity/DG1 pressure USFEM
    bilinear form into separate ``(u, p)`` blocks. ``matrix_kind="mpi"``
    assembles the block forms into an ordinary monolithic PETSc matrix and is
    useful for parity checks against :func:`solve_brinkman_usfem`.
    ``matrix_kind="nest"`` preserves the block structure for PETSc field-split
    experiments. Pass ``preconditioner="diagonal"`` only for iterative
    field-split experiments, not for monolithic direct solves. Results from the
    nested mode should be treated as experimental until compared against a
    direct reference on the same coefficient map.
    """

    if matrix_kind not in {"mpi", "nest"}:
        raise ValueError("matrix_kind must be either 'mpi' or 'nest'")
    if preconditioner not in {"none", "diagonal"}:
        raise ValueError("preconditioner must be either 'none' or 'diagonal'")
    _validate_pressure_drop(pressure_inlet, pressure_outlet)
    _validate_usfem_controls(
        tau_factor=tau_factor,
        m_t=m_t,
        alpha_edge=alpha_edge,
        facet_law=facet_law,
        facet_size_mode=facet_size_mode,
        tau_gamma_cap=tau_gamma_cap,
    )
    nondimensional_options = _resolve_brinkman_nondimensionalization(nondimensional)

    solver_options = options or FEniCSSolverOptions()
    api = _require_dolfinx_core()
    selected_linear_backend = _resolve_linear_backend(solver_options.linear_backend, api)
    if selected_linear_backend != "petsc":
        raise NotImplementedError(
            "solve_brinkman_usfem_block currently supports only the PETSc backend; "
            "use solve_brinkman_usfem for serial direct solves."
        )
    uses_schurdiag_cudss = solver_options.solver_preset == "iterative_schurdiag_cudss_experimental"
    if uses_schurdiag_cudss and preconditioner != "none":
        raise ValueError(
            "preconditioner must be 'none' for the "
            "iterative_schurdiag_cudss_experimental preset; the preset defines "
            "its own lower-Schur preconditioner."
        )
    api = _require_dolfinx_petsc(api)
    context = _build_context(problem, flow_axis=flow_axis, api=api)
    scales = None
    if nondimensional_options is not None:
        scales = _brinkman_nondimensional_scales(
            context,
            problem,
            pressure_inlet=pressure_inlet,
            pressure_outlet=pressure_outlet,
            velocity_scale=nondimensional_options.velocity_scale,
        )

    ufl = cast(_UFLAlgebra, context.api.ufl)
    fem = context.api.fem
    dx = cast(_UFLExpression, context.dx)
    dS = cast(_UFLExpression, context.dS)
    velocity_element = context.api.basix_ufl.element(
        "Lagrange",
        context.mesh.basix_cell(),
        1,
        shape=(context.mesh.geometry.dim,),
    )
    pressure_element = context.api.basix_ufl.element("DG", context.mesh.basix_cell(), 1)
    velocity_space = fem.functionspace(context.mesh, velocity_element)
    pressure_space = fem.functionspace(context.mesh, pressure_element)
    velocity = fem.Function(velocity_space)
    velocity.name = "u"
    pressure = fem.Function(pressure_space)
    pressure.name = "p"

    u = ufl.TrialFunction(velocity_space)
    p = ufl.TrialFunction(pressure_space)
    v = ufl.TestFunction(velocity_space)
    q = ufl.TestFunction(pressure_space)
    gamma_override: _UFLExpression | None = None
    nu_eff_override: _UFLExpression | None = None
    if scales is not None:
        gamma_scaled, nu_eff_scaled = _brinkman_nondimensional_coefficients(
            context,
            problem,
            scales,
        )
        gamma_override = cast(_UFLExpression, gamma_scaled)
        nu_eff_override = cast(_UFLExpression, nu_eff_scaled)
    gamma, nu_eff, tau, tau_f = _usfem_stabilization_terms(
        context,
        tau_factor=tau_factor,
        m_t=m_t,
        alpha_edge=alpha_edge,
        facet_law=facet_law,
        gamma=gamma_override,
        nu_eff=nu_eff_override,
        facet_size=_facet_size_expression(context, facet_size_mode),
        tau_gamma_cap=tau_gamma_cap,
    )

    def residual_velocity_part(w: _UFLExpression) -> _UFLExpression:
        return gamma * w - nu_eff * ufl.div(ufl.grad(w))

    a00 = (
        nu_eff * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
        + ufl.inner(gamma * u, v) * dx
        - tau
        * ufl.inner(
            residual_velocity_part(u),
            residual_velocity_part(v),
        )
        * dx
    )
    a01 = -p * ufl.div(v) * dx - tau * ufl.inner(ufl.grad(p), residual_velocity_part(v)) * dx
    a10 = q * ufl.div(u) * dx + tau * ufl.inner(residual_velocity_part(u), ufl.grad(q)) * dx
    a11 = tau_f * ufl.jump(p) * ufl.jump(q) * dS + tau * ufl.inner(ufl.grad(p), ufl.grad(q)) * dx
    rhs_velocity = _pressure_boundary_load(
        context,
        v,
        flow_axis=flow_axis,
        pressure_inlet=pressure_inlet if scales is None else 1.0,
        pressure_outlet=pressure_outlet if scales is None else 0.0,
    )
    rhs_pressure = _ufl_constant(context, 0.0) * q * dx

    bcs = _velocity_side_wall_bcs(context, velocity_space, flow_axis=flow_axis)
    preconditioner_forms = None
    if preconditioner == "diagonal":
        preconditioner_forms = [[a00, None], [None, a11]]

    if uses_schurdiag_cudss:
        solution, solve_seconds, solver_metadata = _solve_usfem_schurdiag_cudss(
            context,
            forms=[[a00, a01], [a10, a11]],
            rhs=[rhs_velocity, rhs_pressure],
            bcs=bcs,
            solution_functions=[velocity, pressure],
            options=solver_options,
        )
    else:
        solution, solve_seconds, solver_metadata = _solve_block_problem_petsc(
            context,
            forms=[[a00, a01], [a10, a11]],
            rhs=[rhs_velocity, rhs_pressure],
            bcs=bcs,
            solution_functions=[velocity, pressure],
            options=solver_options,
            prefix_suffix=f"brinkman_usfem_block_{flow_axis}",
            matrix_kind=matrix_kind,
            preconditioner_forms=preconditioner_forms,
        )
    if scales is not None:
        solution[0].x.array[:] *= scales.velocity_scale
        solution[0].x.scatter_forward()
        solution[1].x.array[:] *= scales.pressure_scale
        solution[1].x.scatter_forward()
    result = _result_from_velocity_pressure(
        context,
        solution[0],
        solution[1],
        method="Darcy-Brinkman USFEM block CG1 x DG1",
        formulation="brinkman_usfem_p1dg1_block",
        flow_axis=flow_axis,
        pressure_inlet=pressure_inlet,
        pressure_outlet=pressure_outlet,
        viscosity=problem.viscosity,
        solve_seconds=solve_seconds,
        metadata={
            "linear_backend": selected_linear_backend,
            "solver_preset": solver_options.solver_preset,
            "velocity_degree": 1,
            "pressure_family": "DG",
            "pressure_constraint": "natural_traction",
            "returned_pressure_normalization": "zero_mean",
            "porosity_floor": problem.porosity_floor,
            "permeability_floor": problem.permeability_floor,
            "petsc_options": dict(solver_options.petsc_options),
            "petsc_options_prefix": solver_options.petsc_options_prefix,
            "nvmath_cudss_controls": dict(solver_options.nvmath_cudss_controls),
            "iterative_solver_controls": dict(solver_options.iterative_solver_controls),
            "block_matrix_kind": matrix_kind,
            "block_preconditioner": preconditioner,
            "thread_environment": _thread_environment_metadata(),
            **_nondimensional_metadata(scales),
            **_mpi_metadata(context),
            **solver_metadata,
        },
    )
    result.metadata.update(
        {
            "tau_factor": float(tau_factor),
            "m_t": float(m_t),
            "alpha_edge": float(alpha_edge),
            "alpha_edge_active": facet_law in {"classic", "shifted"},
            "facet_law": facet_law,
            "facet_size_mode": facet_size_mode,
            "tau_gamma_cap": (None if tau_gamma_cap is None else float(tau_gamma_cap)),
        }
    )
    return result

---

Taylor-Hood Backends

voids.fem.singlephase.taylorhood

solve_darcy_taylor_hood

solve_darcy_taylor_hood(
    problem,
    *,
    flow_axis="x",
    pressure_inlet=1.0,
    pressure_outlet=0.0,
    options=None,
)

Solve the mixed Darcy-Darcy comparison model with Taylor-Hood elements.

The weak form uses a CG2 velocity approximation and a CG1 pressure approximation:

(mu / K u, v) - (p, div v) + (q, div u) = boundary pressure work.

Side boundaries impose only zero normal velocity on faces transverse to the flow axis. Inlet and outlet pressures are applied as natural traction terms.

Source code in src/voids/fem/singlephase/taylorhood.py

def solve_darcy_taylor_hood(
    problem: FEMMapProblem,
    *,
    flow_axis: str = "x",
    pressure_inlet: float = 1.0,
    pressure_outlet: float = 0.0,
    options: FEniCSSolverOptions | None = None,
) -> FEMSinglePhaseResult:
    """Solve the mixed Darcy-Darcy comparison model with Taylor-Hood elements.

    The weak form uses a CG2 velocity approximation and a CG1 pressure
    approximation:

    ``(mu / K u, v) - (p, div v) + (q, div u) = boundary pressure work``.

    Side boundaries impose only zero normal velocity on faces transverse to the
    flow axis. Inlet and outlet pressures are applied as natural traction terms.
    """

    def form_builder(context, u, p, v, q):
        ufl = context.api.ufl
        gamma = context.coefficients["gamma"]
        return (
            ufl.inner(gamma * u, v) * context.dx
            - p * ufl.div(v) * context.dx
            + q * ufl.div(u) * context.dx
        )

    return _solve_with_form_builder(
        problem,
        flow_axis=flow_axis,
        pressure_inlet=pressure_inlet,
        pressure_outlet=pressure_outlet,
        options=options,
        velocity_degree=2,
        pressure_family="Lagrange",
        method="Darcy-Darcy Taylor-Hood CG2 x CG1",
        formulation="darcy_taylor_hood_p2p1",
        prefix_suffix=f"darcy_taylor_hood_{flow_axis}",
        form_builder=form_builder,
    )

solve_brinkman_taylor_hood

solve_brinkman_taylor_hood(
    problem,
    *,
    flow_axis="x",
    pressure_inlet=1.0,
    pressure_outlet=0.0,
    options=None,
    nondimensional=False,
)

Solve the Darcy-Brinkman micro-continuum model with Taylor-Hood elements.

The weak form uses CG2 velocity and CG1 pressure:

(mu / phi grad u, grad v) + (mu / K u, v) - (p, div v) + (q, div u) = boundary pressure work.

K and phi are piecewise-constant maps supplied through :class:~voids.fem.singlephase.FEMMapProblem.

Source code in src/voids/fem/singlephase/taylorhood.py

def solve_brinkman_taylor_hood(
    problem: FEMMapProblem,
    *,
    flow_axis: str = "x",
    pressure_inlet: float = 1.0,
    pressure_outlet: float = 0.0,
    options: FEniCSSolverOptions | None = None,
    nondimensional: bool | BrinkmanNondimensionalization = False,
) -> FEMSinglePhaseResult:
    """Solve the Darcy-Brinkman micro-continuum model with Taylor-Hood elements.

    The weak form uses CG2 velocity and CG1 pressure:

    ``(mu / phi grad u, grad v) + (mu / K u, v)
    - (p, div v) + (q, div u) = boundary pressure work``.

    ``K`` and ``phi`` are piecewise-constant maps supplied through
    :class:`~voids.fem.singlephase.FEMMapProblem`.
    """

    nondimensional_options = _resolve_brinkman_nondimensionalization(nondimensional)
    scales = None
    if nondimensional_options is not None:
        _validate_pressure_drop(pressure_inlet, pressure_outlet)
        context_for_scales = _build_context(
            problem,
            flow_axis=flow_axis,
            api=_require_dolfinx_core(),
        )
        scales = _brinkman_nondimensional_scales(
            context_for_scales,
            problem,
            pressure_inlet=pressure_inlet,
            pressure_outlet=pressure_outlet,
            velocity_scale=nondimensional_options.velocity_scale,
        )

    def form_builder(context, u, p, v, q):
        ufl = context.api.ufl
        if scales is None:
            gamma = context.coefficients["gamma"]
            nu_eff = context.coefficients["nu_eff"]
        else:
            gamma, nu_eff = _brinkman_nondimensional_coefficients(context, problem, scales)
        return (
            nu_eff * ufl.inner(ufl.grad(u), ufl.grad(v)) * context.dx
            + ufl.inner(gamma * u, v) * context.dx
            - p * ufl.div(v) * context.dx
            + q * ufl.div(u) * context.dx
        )

    return _solve_with_form_builder(
        problem,
        flow_axis=flow_axis,
        pressure_inlet=pressure_inlet,
        pressure_outlet=pressure_outlet,
        options=options,
        velocity_degree=2,
        pressure_family="Lagrange",
        method="Darcy-Brinkman Taylor-Hood CG2 x CG1",
        formulation="brinkman_taylor_hood_p2p1",
        prefix_suffix=f"brinkman_taylor_hood_{flow_axis}",
        form_builder=form_builder,
        boundary_pressure_inlet=None if scales is None else 1.0,
        boundary_pressure_outlet=None if scales is None else 0.0,
        velocity_scale=1.0 if scales is None else scales.velocity_scale,
        pressure_scale=1.0 if scales is None else scales.pressure_scale,
        extra_metadata=_nondimensional_metadata(scales),
    )

---

USFEM Backends

voids.fem.singlephase.usfem

solve_brinkman_usfem

solve_brinkman_usfem(
    problem,
    *,
    flow_axis="x",
    pressure_inlet=1.0,
    pressure_outlet=0.0,
    tau_factor=1.0,
    m_t=1.0 / 3.0,
    alpha_edge=1.0,
    facet_law="reaction_diffusion",
    facet_size_mode="cell_diameter",
    tau_gamma_cap=None,
    pressure_degree=1,
    options=None,
    nondimensional=False,
)

Solve a stabilized Darcy-Brinkman micro-continuum model.

The formulation uses CG1 velocity and discontinuous pressure fields. DG1 is the backward-compatible default; pass pressure_degree=0 for the low-order CG1 x DG0 pair used in manufactured-solution and vug studies. It augments the Brinkman weak form with a residual-based cell stabilization term and an interior pressure-jump penalty. The coefficients are intended for porosity/permeability maps obtained from a segmented image.

Set tau_factor=0 to disable the cell residual while retaining the pressure-jump term. In high-contrast CG1 x DG0 studies, tau_gamma_cap can enforce gamma * tau_K <= tau_gamma_cap and avoid near-cancellation of the physical drag. facet_size_mode="facet_measure" uses edge length in 2D and the square root of facet area in 3D; the latter is a measure-based length, not an exact triangular-facet diameter.

Source code in src/voids/fem/singlephase/usfem.py

def solve_brinkman_usfem(
    problem: FEMMapProblem,
    *,
    flow_axis: str = "x",
    pressure_inlet: float = 1.0,
    pressure_outlet: float = 0.0,
    tau_factor: float = 1.0,
    m_t: float = 1.0 / 3.0,
    alpha_edge: float = 1.0,
    facet_law: USFEMFacetLaw = "reaction_diffusion",
    facet_size_mode: USFEMFacetSizeMode = "cell_diameter",
    tau_gamma_cap: float | None = None,
    pressure_degree: Literal[0, 1] = 1,
    options: FEniCSSolverOptions | None = None,
    nondimensional: bool | BrinkmanNondimensionalization = False,
) -> FEMSinglePhaseResult:
    """Solve a stabilized Darcy-Brinkman micro-continuum model.

    The formulation uses CG1 velocity and discontinuous pressure fields. DG1 is
    the backward-compatible default; pass ``pressure_degree=0`` for the
    low-order CG1 x DG0 pair used in manufactured-solution and vug studies. It
    augments the Brinkman weak form with a residual-based cell stabilization term and an
    interior pressure-jump penalty. The coefficients are intended for
    porosity/permeability maps obtained from a segmented image.

    Set ``tau_factor=0`` to disable the cell residual while retaining the
    pressure-jump term. In high-contrast CG1 x DG0 studies,
    ``tau_gamma_cap`` can enforce ``gamma * tau_K <= tau_gamma_cap`` and avoid
    near-cancellation of the physical drag. ``facet_size_mode="facet_measure"``
    uses edge length in 2D and the square root of facet area in 3D; the latter
    is a measure-based length, not an exact triangular-facet diameter.
    """

    _validate_usfem_controls(
        tau_factor=tau_factor,
        m_t=m_t,
        alpha_edge=alpha_edge,
        facet_law=facet_law,
        facet_size_mode=facet_size_mode,
        tau_gamma_cap=tau_gamma_cap,
    )
    if pressure_degree not in {0, 1}:
        raise ValueError("pressure_degree must be either 0 or 1")
    uncapped_max_tau_gamma = (
        _p1dg0_max_uncapped_tau_gamma(
            problem,
            tau_factor=tau_factor,
            m_t=m_t,
        )
        if pressure_degree == 0
        else None
    )
    if (
        uncapped_max_tau_gamma is not None
        and tau_factor > 0.0
        and tau_gamma_cap is None
        and uncapped_max_tau_gamma >= 0.9
    ):
        warnings.warn(
            "The uncapped CG1 x DG0 cell term has estimated "
            f"max(gamma * tau_K)={uncapped_max_tau_gamma:.3g}; it can nearly "
            "cancel physical drag. Set tau_gamma_cap below 1 or use "
            "tau_factor=0 as an explicit sensitivity branch.",
            RuntimeWarning,
            stacklevel=2,
        )
    nondimensional_options = _resolve_brinkman_nondimensionalization(nondimensional)
    scales = None
    if nondimensional_options is not None:
        _validate_pressure_drop(pressure_inlet, pressure_outlet)
        context_for_scales = _build_context(
            problem,
            flow_axis=flow_axis,
            api=_require_dolfinx_core(),
        )
        scales = _brinkman_nondimensional_scales(
            context_for_scales,
            problem,
            pressure_inlet=pressure_inlet,
            pressure_outlet=pressure_outlet,
            velocity_scale=nondimensional_options.velocity_scale,
        )

    def form_builder(
        context: _FEMContext,
        u: _UFLExpression,
        p: _UFLExpression,
        v: _UFLExpression,
        q: _UFLExpression,
    ) -> _UFLExpression:
        ufl = cast(_UFLAlgebra, context.api.ufl)
        dx = cast(_UFLExpression, context.dx)
        dS = cast(_UFLExpression, context.dS)
        gamma_override: _UFLExpression | None = None
        nu_eff_override: _UFLExpression | None = None
        if scales is not None:
            gamma_raw, nu_eff_raw = _brinkman_nondimensional_coefficients(
                context,
                problem,
                scales,
            )
            gamma_override = cast(_UFLExpression, gamma_raw)
            nu_eff_override = cast(_UFLExpression, nu_eff_raw)
        gamma, nu_eff, tau, tau_f = _usfem_stabilization_terms(
            context,
            tau_factor=tau_factor,
            m_t=m_t,
            alpha_edge=alpha_edge,
            facet_law=facet_law,
            gamma=gamma_override,
            nu_eff=nu_eff_override,
            facet_size=_facet_size_expression(context, facet_size_mode),
            tau_gamma_cap=tau_gamma_cap,
        )
        residual_u = gamma * u + ufl.grad(p) - nu_eff * ufl.div(ufl.grad(u))
        residual_vq = gamma * v - ufl.grad(q) - nu_eff * ufl.div(ufl.grad(v))
        return (
            nu_eff * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
            + ufl.inner(gamma * u, v) * dx
            - p * ufl.div(v) * dx
            + q * ufl.div(u) * dx
            + tau_f * ufl.jump(p) * ufl.jump(q) * dS
            - tau * ufl.inner(residual_u, residual_vq) * dx
        )

    result = _solve_with_form_builder(
        problem,
        flow_axis=flow_axis,
        pressure_inlet=pressure_inlet,
        pressure_outlet=pressure_outlet,
        options=options,
        velocity_degree=1,
        pressure_family="DG",
        pressure_degree=pressure_degree,
        method=f"Darcy-Brinkman USFEM CG1 x DG{pressure_degree}",
        formulation=f"brinkman_usfem_p1dg{pressure_degree}",
        prefix_suffix=f"brinkman_usfem_{flow_axis}",
        form_builder=form_builder,
        boundary_pressure_inlet=None if scales is None else 1.0,
        boundary_pressure_outlet=None if scales is None else 0.0,
        velocity_scale=1.0 if scales is None else scales.velocity_scale,
        pressure_scale=1.0 if scales is None else scales.pressure_scale,
        extra_metadata=_nondimensional_metadata(scales),
    )
    result.metadata.update(
        {
            "tau_factor": float(tau_factor),
            "m_t": float(m_t),
            "alpha_edge": float(alpha_edge),
            "alpha_edge_active": facet_law in {"classic", "shifted"},
            "facet_law": facet_law,
            "facet_size_mode": facet_size_mode,
            "tau_gamma_cap": (None if tau_gamma_cap is None else float(tau_gamma_cap)),
            "p1dg0_uncapped_max_tau_gamma": uncapped_max_tau_gamma,
        }
    )
    return result

solve_brinkman_usfem_block

solve_brinkman_usfem_block(
    problem,
    *,
    flow_axis="x",
    pressure_inlet=1.0,
    pressure_outlet=0.0,
    tau_factor=1.0,
    m_t=1.0 / 3.0,
    alpha_edge=1.0,
    facet_law="reaction_diffusion",
    facet_size_mode="cell_diameter",
    tau_gamma_cap=None,
    options=None,
    matrix_kind="mpi",
    preconditioner="none",
    nondimensional=False,
)

Solve USFEM through explicit velocity/pressure block forms.

This experimental path expands the same CG1 velocity/DG1 pressure USFEM bilinear form into separate (u, p) blocks. matrix_kind="mpi" assembles the block forms into an ordinary monolithic PETSc matrix and is useful for parity checks against :func:solve_brinkman_usfem. matrix_kind="nest" preserves the block structure for PETSc field-split experiments. Pass preconditioner="diagonal" only for iterative field-split experiments, not for monolithic direct solves. Results from the nested mode should be treated as experimental until compared against a direct reference on the same coefficient map.

Source code in src/voids/fem/singlephase/usfem.py

def solve_brinkman_usfem_block(
    problem: FEMMapProblem,
    *,
    flow_axis: str = "x",
    pressure_inlet: float = 1.0,
    pressure_outlet: float = 0.0,
    tau_factor: float = 1.0,
    m_t: float = 1.0 / 3.0,
    alpha_edge: float = 1.0,
    facet_law: USFEMFacetLaw = "reaction_diffusion",
    facet_size_mode: USFEMFacetSizeMode = "cell_diameter",
    tau_gamma_cap: float | None = None,
    options: FEniCSSolverOptions | None = None,
    matrix_kind: Literal["mpi", "nest"] = "mpi",
    preconditioner: Literal["none", "diagonal"] = "none",
    nondimensional: bool | BrinkmanNondimensionalization = False,
) -> FEMSinglePhaseResult:
    """Solve USFEM through explicit velocity/pressure block forms.

    This experimental path expands the same CG1 velocity/DG1 pressure USFEM
    bilinear form into separate ``(u, p)`` blocks. ``matrix_kind="mpi"``
    assembles the block forms into an ordinary monolithic PETSc matrix and is
    useful for parity checks against :func:`solve_brinkman_usfem`.
    ``matrix_kind="nest"`` preserves the block structure for PETSc field-split
    experiments. Pass ``preconditioner="diagonal"`` only for iterative
    field-split experiments, not for monolithic direct solves. Results from the
    nested mode should be treated as experimental until compared against a
    direct reference on the same coefficient map.
    """

    if matrix_kind not in {"mpi", "nest"}:
        raise ValueError("matrix_kind must be either 'mpi' or 'nest'")
    if preconditioner not in {"none", "diagonal"}:
        raise ValueError("preconditioner must be either 'none' or 'diagonal'")
    _validate_pressure_drop(pressure_inlet, pressure_outlet)
    _validate_usfem_controls(
        tau_factor=tau_factor,
        m_t=m_t,
        alpha_edge=alpha_edge,
        facet_law=facet_law,
        facet_size_mode=facet_size_mode,
        tau_gamma_cap=tau_gamma_cap,
    )
    nondimensional_options = _resolve_brinkman_nondimensionalization(nondimensional)

    solver_options = options or FEniCSSolverOptions()
    api = _require_dolfinx_core()
    selected_linear_backend = _resolve_linear_backend(solver_options.linear_backend, api)
    if selected_linear_backend != "petsc":
        raise NotImplementedError(
            "solve_brinkman_usfem_block currently supports only the PETSc backend; "
            "use solve_brinkman_usfem for serial direct solves."
        )
    uses_schurdiag_cudss = solver_options.solver_preset == "iterative_schurdiag_cudss_experimental"
    if uses_schurdiag_cudss and preconditioner != "none":
        raise ValueError(
            "preconditioner must be 'none' for the "
            "iterative_schurdiag_cudss_experimental preset; the preset defines "
            "its own lower-Schur preconditioner."
        )
    api = _require_dolfinx_petsc(api)
    context = _build_context(problem, flow_axis=flow_axis, api=api)
    scales = None
    if nondimensional_options is not None:
        scales = _brinkman_nondimensional_scales(
            context,
            problem,
            pressure_inlet=pressure_inlet,
            pressure_outlet=pressure_outlet,
            velocity_scale=nondimensional_options.velocity_scale,
        )

    ufl = cast(_UFLAlgebra, context.api.ufl)
    fem = context.api.fem
    dx = cast(_UFLExpression, context.dx)
    dS = cast(_UFLExpression, context.dS)
    velocity_element = context.api.basix_ufl.element(
        "Lagrange",
        context.mesh.basix_cell(),
        1,
        shape=(context.mesh.geometry.dim,),
    )
    pressure_element = context.api.basix_ufl.element("DG", context.mesh.basix_cell(), 1)
    velocity_space = fem.functionspace(context.mesh, velocity_element)
    pressure_space = fem.functionspace(context.mesh, pressure_element)
    velocity = fem.Function(velocity_space)
    velocity.name = "u"
    pressure = fem.Function(pressure_space)
    pressure.name = "p"

    u = ufl.TrialFunction(velocity_space)
    p = ufl.TrialFunction(pressure_space)
    v = ufl.TestFunction(velocity_space)
    q = ufl.TestFunction(pressure_space)
    gamma_override: _UFLExpression | None = None
    nu_eff_override: _UFLExpression | None = None
    if scales is not None:
        gamma_scaled, nu_eff_scaled = _brinkman_nondimensional_coefficients(
            context,
            problem,
            scales,
        )
        gamma_override = cast(_UFLExpression, gamma_scaled)
        nu_eff_override = cast(_UFLExpression, nu_eff_scaled)
    gamma, nu_eff, tau, tau_f = _usfem_stabilization_terms(
        context,
        tau_factor=tau_factor,
        m_t=m_t,
        alpha_edge=alpha_edge,
        facet_law=facet_law,
        gamma=gamma_override,
        nu_eff=nu_eff_override,
        facet_size=_facet_size_expression(context, facet_size_mode),
        tau_gamma_cap=tau_gamma_cap,
    )

    def residual_velocity_part(w: _UFLExpression) -> _UFLExpression:
        return gamma * w - nu_eff * ufl.div(ufl.grad(w))

    a00 = (
        nu_eff * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
        + ufl.inner(gamma * u, v) * dx
        - tau
        * ufl.inner(
            residual_velocity_part(u),
            residual_velocity_part(v),
        )
        * dx
    )
    a01 = -p * ufl.div(v) * dx - tau * ufl.inner(ufl.grad(p), residual_velocity_part(v)) * dx
    a10 = q * ufl.div(u) * dx + tau * ufl.inner(residual_velocity_part(u), ufl.grad(q)) * dx
    a11 = tau_f * ufl.jump(p) * ufl.jump(q) * dS + tau * ufl.inner(ufl.grad(p), ufl.grad(q)) * dx
    rhs_velocity = _pressure_boundary_load(
        context,
        v,
        flow_axis=flow_axis,
        pressure_inlet=pressure_inlet if scales is None else 1.0,
        pressure_outlet=pressure_outlet if scales is None else 0.0,
    )
    rhs_pressure = _ufl_constant(context, 0.0) * q * dx

    bcs = _velocity_side_wall_bcs(context, velocity_space, flow_axis=flow_axis)
    preconditioner_forms = None
    if preconditioner == "diagonal":
        preconditioner_forms = [[a00, None], [None, a11]]

    if uses_schurdiag_cudss:
        solution, solve_seconds, solver_metadata = _solve_usfem_schurdiag_cudss(
            context,
            forms=[[a00, a01], [a10, a11]],
            rhs=[rhs_velocity, rhs_pressure],
            bcs=bcs,
            solution_functions=[velocity, pressure],
            options=solver_options,
        )
    else:
        solution, solve_seconds, solver_metadata = _solve_block_problem_petsc(
            context,
            forms=[[a00, a01], [a10, a11]],
            rhs=[rhs_velocity, rhs_pressure],
            bcs=bcs,
            solution_functions=[velocity, pressure],
            options=solver_options,
            prefix_suffix=f"brinkman_usfem_block_{flow_axis}",
            matrix_kind=matrix_kind,
            preconditioner_forms=preconditioner_forms,
        )
    if scales is not None:
        solution[0].x.array[:] *= scales.velocity_scale
        solution[0].x.scatter_forward()
        solution[1].x.array[:] *= scales.pressure_scale
        solution[1].x.scatter_forward()
    result = _result_from_velocity_pressure(
        context,
        solution[0],
        solution[1],
        method="Darcy-Brinkman USFEM block CG1 x DG1",
        formulation="brinkman_usfem_p1dg1_block",
        flow_axis=flow_axis,
        pressure_inlet=pressure_inlet,
        pressure_outlet=pressure_outlet,
        viscosity=problem.viscosity,
        solve_seconds=solve_seconds,
        metadata={
            "linear_backend": selected_linear_backend,
            "solver_preset": solver_options.solver_preset,
            "velocity_degree": 1,
            "pressure_family": "DG",
            "pressure_constraint": "natural_traction",
            "returned_pressure_normalization": "zero_mean",
            "porosity_floor": problem.porosity_floor,
            "permeability_floor": problem.permeability_floor,
            "petsc_options": dict(solver_options.petsc_options),
            "petsc_options_prefix": solver_options.petsc_options_prefix,
            "nvmath_cudss_controls": dict(solver_options.nvmath_cudss_controls),
            "iterative_solver_controls": dict(solver_options.iterative_solver_controls),
            "block_matrix_kind": matrix_kind,
            "block_preconditioner": preconditioner,
            "thread_environment": _thread_environment_metadata(),
            **_nondimensional_metadata(scales),
            **_mpi_metadata(context),
            **solver_metadata,
        },
    )
    result.metadata.update(
        {
            "tau_factor": float(tau_factor),
            "m_t": float(m_t),
            "alpha_edge": float(alpha_edge),
            "alpha_edge_active": facet_law in {"classic", "shifted"},
            "facet_law": facet_law,
            "facet_size_mode": facet_size_mode,
            "tau_gamma_cap": (None if tau_gamma_cap is None else float(tau_gamma_cap)),
        }
    )
    return result

---

Upscaling

voids.fem.singlephase.upscaling

FEMUpscalingResult dataclass

Principal-direction FEM micro-continuum permeability result.

Source code in src/voids/fem/singlephase/upscaling.py

@dataclass(slots=True)
class FEMUpscalingResult:
    """Principal-direction FEM micro-continuum permeability result."""

    results: dict[str, FEMSinglePhaseResult]
    backend: str

    @property
    def permeability(self) -> dict[str, float]:
        """Return effective permeability by principal axis."""

        return {axis: result.permeability for axis, result in self.results.items()}

    @property
    def solve_seconds(self) -> dict[str, float]:
        """Return wall-clock solve time by principal axis."""

        return {axis: result.solve_seconds for axis, result in self.results.items()}

permeability property

permeability

Return effective permeability by principal axis.

solve_seconds property

solve_seconds

Return wall-clock solve time by principal axis.

upscale_permeability_fem

upscale_permeability_fem(
    problem,
    *,
    backend="taylor_hood_brinkman",
    axes=None,
    pressure_inlet=1.0,
    pressure_outlet=0.0,
    options=None,
    backend_kwargs=None,
)

Compute principal-direction FEM permeability estimates.

Parameters:

Name	Type	Description	Default
problem	FEMMapProblem	Porosity/permeability map problem.	required
backend	FEMBackend	Either a backend name or a compatible solver callable. Supported names are "taylor_hood_brinkman", "taylor_hood_darcy", and "usfem_brinkman".	'taylor_hood_brinkman'
axes	tuple[str, ...] | None	Principal axes to solve. By default all axes supported by the map dimensionality are solved.	None
pressure_inlet	float	Natural pressure values imposed on opposite faces of each flow axis.	1.0
pressure_outlet	float	Natural pressure values imposed on opposite faces of each flow axis.	1.0
options	FEniCSSolverOptions | None	Linear solver options passed to the FEM backend. The default preserves PETSc where available and uses a serial SciPy direct solve on native Windows when PETSc is unavailable.	None
backend_kwargs	dict[str, object] | None	Additional backend-specific keyword arguments.	None

Source code in src/voids/fem/singlephase/upscaling.py

def upscale_permeability_fem(
    problem: FEMMapProblem,
    *,
    backend: FEMBackend = "taylor_hood_brinkman",
    axes: tuple[str, ...] | None = None,
    pressure_inlet: float = 1.0,
    pressure_outlet: float = 0.0,
    options: FEniCSSolverOptions | None = None,
    backend_kwargs: dict[str, object] | None = None,
) -> FEMUpscalingResult:
    """Compute principal-direction FEM permeability estimates.

    Parameters
    ----------
    problem :
        Porosity/permeability map problem.
    backend :
        Either a backend name or a compatible solver callable. Supported names
        are ``"taylor_hood_brinkman"``, ``"taylor_hood_darcy"``, and
        ``"usfem_brinkman"``.
    axes :
        Principal axes to solve. By default all axes supported by the map
        dimensionality are solved.
    pressure_inlet, pressure_outlet :
        Natural pressure values imposed on opposite faces of each flow axis.
    options :
        Linear solver options passed to the FEM backend. The default preserves
        PETSc where available and uses a serial SciPy direct solve on native
        Windows when PETSc is unavailable.
    backend_kwargs :
        Additional backend-specific keyword arguments.
    """

    solve_axes = axes or _default_axes(problem.permeability_map.ndim)
    solver = _backend_from_name(backend) if isinstance(backend, str) else backend
    kwargs = dict(backend_kwargs or {})
    results = {
        axis: solver(
            problem,
            flow_axis=axis,
            pressure_inlet=pressure_inlet,
            pressure_outlet=pressure_outlet,
            options=options,
            **kwargs,
        )
        for axis in solve_axes
    }
    backend_name = backend if isinstance(backend, str) else getattr(backend, "__name__", "callable")
    return FEMUpscalingResult(results=results, backend=backend_name)

upscale_principal_permeabilities_fem

upscale_principal_permeabilities_fem(
    problem,
    *,
    backend="taylor_hood_brinkman",
    axes=None,
    pressure_inlet=1.0,
    pressure_outlet=0.0,
    options=None,
    backend_kwargs=None,
)

Return only the principal FEM permeability values.

Source code in src/voids/fem/singlephase/upscaling.py

def upscale_principal_permeabilities_fem(
    problem: FEMMapProblem,
    *,
    backend: FEMBackend = "taylor_hood_brinkman",
    axes: tuple[str, ...] | None = None,
    pressure_inlet: float = 1.0,
    pressure_outlet: float = 0.0,
    options: FEniCSSolverOptions | None = None,
    backend_kwargs: dict[str, object] | None = None,
) -> dict[str, float]:
    """Return only the principal FEM permeability values."""

    return upscale_permeability_fem(
        problem,
        backend=backend,
        axes=axes,
        pressure_inlet=pressure_inlet,
        pressure_outlet=pressure_outlet,
        options=options,
        backend_kwargs=backend_kwargs,
    ).permeability