Skip to content

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.

Common Types

voids.fem.singlephase

Single-phase finite-element Darcy and Brinkman backends.

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 
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
@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 
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
@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's direct sparse solver.

Use linear_backend="scipy" or "umfpack" to request the serial DOLFINx-assembly/direct-sparse path explicitly on any platform. These paths use the same weak form and boundary conditions as PETSc; only the linear algebra backend changes. "umfpack" requires the optional scikits.umfpack package.

Source code in src/voids/fem/singlephase/_common.py 
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
@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's direct sparse solver.

    Use ``linear_backend="scipy"`` or ``"umfpack"`` to request the serial
    DOLFINx-assembly/direct-sparse path explicitly on any platform. These paths
    use the same weak form and boundary conditions as PETSc; only the linear
    algebra backend changes. ``"umfpack"`` requires the optional
    ``scikits.umfpack`` package.
    """

    linear_backend: LinearSolverBackend = "auto"
    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,
            petsc_options=options,
            petsc_options_prefix=petsc_options_prefix,
        )

    @classmethod
    def scipy_direct(cls) -> FEniCSSolverOptions:
        """Create options for the serial DOLFINx-assembly/SciPy direct backend."""

        return cls(linear_backend="scipy")

    @classmethod
    def umfpack_direct(cls) -> FEniCSSolverOptions:
        """Create options for the serial DOLFINx-assembly/UMFPACK backend."""

        return cls(linear_backend="umfpack")

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 
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
@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,
        petsc_options=options,
        petsc_options_prefix=petsc_options_prefix,
    )

scipy_direct classmethod 

scipy_direct()

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

Source code in src/voids/fem/singlephase/_common.py 
118
119
120
121
122
@classmethod
def scipy_direct(cls) -> FEniCSSolverOptions:
    """Create options for the serial DOLFINx-assembly/SciPy direct backend."""

    return cls(linear_backend="scipy")

umfpack_direct classmethod 

umfpack_direct()

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

Source code in src/voids/fem/singlephase/_common.py 
124
125
126
127
128
@classmethod
def umfpack_direct(cls) -> FEniCSSolverOptions:
    """Create options for the serial DOLFINx-assembly/UMFPACK backend."""

    return cls(linear_backend="umfpack")

FEMUpscalingResult dataclass

Principal-direction FEM micro-continuum permeability result.

Source code in src/voids/fem/singlephase/upscaling.py 
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
@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,
)

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 
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
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,
) -> 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`.
    """

    def form_builder(context, u, p, v, q):
        ufl = context.api.ufl
        gamma = context.coefficients["gamma"]
        nu_eff = context.coefficients["nu_eff"]
        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,
    )

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 
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
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 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
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 
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
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,
    options=None,
)

Solve a stabilized Darcy-Brinkman micro-continuum model.

The formulation uses CG1 velocity and DG1 pressure fields. 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.

Source code in src/voids/fem/singlephase/usfem.py 
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
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,
    options: FEniCSSolverOptions | None = None,
) -> FEMSinglePhaseResult:
    """Solve a stabilized Darcy-Brinkman micro-continuum model.

    The formulation uses CG1 velocity and DG1 pressure fields. 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.
    """

    if tau_factor <= 0.0:
        raise ValueError("tau_factor must be positive")
    if m_t <= 0.0:
        raise ValueError("m_t must be positive")
    if alpha_edge <= 0.0:
        raise ValueError("alpha_edge must be positive")

    def form_builder(context, u, p, v, q):
        ufl = context.api.ufl
        fem = context.api.fem
        gamma = context.coefficients["gamma"]
        nu_eff = context.coefficients["nu_eff"]
        h = ufl.CellDiameter(context.mesh)
        h_f = ufl.avg(h)
        tau = fem.Constant(context.mesh, float(tau_factor)) * _paper_tau(
            context,
            h,
            gamma,
            nu_eff,
            m_t=m_t,
        )
        tau_f = _interior_pressure_tau(
            context,
            h_f,
            gamma,
            nu_eff,
            alpha_edge=alpha_edge,
        )
        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)) * context.dx
            + ufl.inner(gamma * u, v) * context.dx
            - p * ufl.div(v) * context.dx
            + q * ufl.div(u) * context.dx
            + tau_f * ufl.jump(p) * ufl.jump(q) * context.dS
            - tau * ufl.inner(residual_u, residual_vq) * context.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",
        method="Darcy-Brinkman USFEM CG1 x DG1",
        formulation="brinkman_usfem_p1dg1",
        prefix_suffix=f"brinkman_usfem_{flow_axis}",
        form_builder=form_builder,
    )
    result.metadata.update(
        {
            "tau_factor": float(tau_factor),
            "m_t": float(m_t),
            "alpha_edge": float(alpha_edge),
        }
    )
    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 
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
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,
)

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 
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
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,
) -> 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`.
    """

    def form_builder(context, u, p, v, q):
        ufl = context.api.ufl
        gamma = context.coefficients["gamma"]
        nu_eff = context.coefficients["nu_eff"]
        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,
    )

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,
    options=None,
)

Solve a stabilized Darcy-Brinkman micro-continuum model.

The formulation uses CG1 velocity and DG1 pressure fields. 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.

Source code in src/voids/fem/singlephase/usfem.py 
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
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,
    options: FEniCSSolverOptions | None = None,
) -> FEMSinglePhaseResult:
    """Solve a stabilized Darcy-Brinkman micro-continuum model.

    The formulation uses CG1 velocity and DG1 pressure fields. 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.
    """

    if tau_factor <= 0.0:
        raise ValueError("tau_factor must be positive")
    if m_t <= 0.0:
        raise ValueError("m_t must be positive")
    if alpha_edge <= 0.0:
        raise ValueError("alpha_edge must be positive")

    def form_builder(context, u, p, v, q):
        ufl = context.api.ufl
        fem = context.api.fem
        gamma = context.coefficients["gamma"]
        nu_eff = context.coefficients["nu_eff"]
        h = ufl.CellDiameter(context.mesh)
        h_f = ufl.avg(h)
        tau = fem.Constant(context.mesh, float(tau_factor)) * _paper_tau(
            context,
            h,
            gamma,
            nu_eff,
            m_t=m_t,
        )
        tau_f = _interior_pressure_tau(
            context,
            h_f,
            gamma,
            nu_eff,
            alpha_edge=alpha_edge,
        )
        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)) * context.dx
            + ufl.inner(gamma * u, v) * context.dx
            - p * ufl.div(v) * context.dx
            + q * ufl.div(u) * context.dx
            + tau_f * ufl.jump(p) * ufl.jump(q) * context.dS
            - tau * ufl.inner(residual_u, residual_vq) * context.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",
        method="Darcy-Brinkman USFEM CG1 x DG1",
        formulation="brinkman_usfem_p1dg1",
        prefix_suffix=f"brinkman_usfem_{flow_axis}",
        form_builder=form_builder,
    )
    result.metadata.update(
        {
            "tau_factor": float(tau_factor),
            "m_t": float(m_t),
            "alpha_edge": float(alpha_edge),
        }
    )
    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 
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
@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 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
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 
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
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