Transport, reaction and diffusion

Stationary and transient scalar problems, conservative transport and stabilization.

All API families

User-defined formulations compose these public numerical owners with LocalEquations and Equation; their editable providers are listed in the variational guide. Physical field records provide optional interpretation of the resulting coefficients. Named FieldDefinition views require no physical solver object.

Conservative RAD uses the Robin multiplier \(\lambda_K=(-K\nabla u+\beta u/2)\cdot n_K\), while its physical flux is \(j=-K\nabla u+\beta u\). The low-level p1_operators function uses the conventional advective volume form; triangle_transport_operators and tetra_transport_operators supply the conservative skew form. Scalar UNUSUAL subtracts the complete reaction–diffusion residual pairing and its source term; SUPG instead tests the conservative residual along the velocity. The original local/global equations and analytical assumptions are given by Harder, Paredes and Valentin (2015) and Araya et al. (2024). The API below documents the implemented operators and boundary conventions.

Three-dimensional tetrahedral and polyhedral examples share the editable examples/formulations/scalar_3d.py four-block provider. The importable define_transport_3d composition selects only geometry and trace integration; scalar_constraints_3d builds the complete physical integral when the declared boundary/operator admits a constant gauge. The tetrahedral notebook and polyhedral notebook show assembly and recovery without a prepared method solver. Selective local counts come from the actual coefficient samples and boundary tangency, without coordinator operator assembly or numerical nullspace thresholds.

Small, auditable triangular finite-element assembly and quadrature kernels.

triangle_quadrature(order: int = 4) -> tuple[FloatArray, FloatArray]

Return barycentric Duffy-product Gauss points and unit-sum weights.

p1_geometry(mesh: TriangleMesh) -> tuple[FloatArray, FloatArray]

Return affine P1 gradients per cell and positive integration areas.

p1_operators(mesh: TriangleMesh, diffusion: Any = 1.0, source: Any = 0.0, reaction: Any = 0.0, advection: Any = (0.0, 0.0), order: int = 4) -> tuple[Any, Any, FloatArray]

Assemble diffusion, mass and load with conventional advective transport.

The bilinear form is (K grad u,grad v)+(beta.grad u,v)+(c u,v). This operator is not the skew conservative Robin formulation of RAD-MHM.

face_integration(mesh: TriangleMesh, cell: int, fine: TriangleMesh, skeleton: SkeletonSpace) -> tuple[FloatArray, FloatArray]

Integrate skeleton basis against P1 traces and fine boundary flux DOFs.

Quadrature is split at both fine-mesh vertices and skeleton breaks. The second result maps trace coefficients to integrated outward RT0 boundary fluxes, ordered by fine.boundary_faces.

boundary_data(skeleton: SkeletonSpace, dirichlet: Any, neumann: dict[int, Any] | None = None, *, neumann_components: dict[int, dict[int, Any]] | None = None, order: int = 5) -> tuple[FloatArray, dict[int, float]]

Integrate primary fields and L2-project full or componentwise boundary fluxes.

neumann_components[face][component] prescribes a Cartesian component of the multiplier as a scalar field. Other components retain Dirichlet data. Full and componentwise Neumann declarations cannot share a face. Components refer to the skeleton's interleaved Cartesian ordering; no normal/tangential rotation or outward-sign conversion is performed here.

rt0_operators(mesh: TriangleMesh, permeability: Any = 1.0, source: Any = 0.0, order: int = 4) -> tuple[Any, Any, FloatArray]

Assemble RT0 flux mass, integrated divergence and P0 source moments.

A flux DOF is the integral across a globally oriented face. In a cell the basis associated with side (i,j) is sign*(x-opposite)/(2*area).

rt0_evaluate(mesh: TriangleMesh, flux: FloatArray, barycentric: FloatArray) -> FloatArray

Evaluate RT0 at common (q,3) or cellwise (t,q,3) barycentric points.

Physical positive residual scales for scalar finite-element forms.

streamline_scale(fine: TriangleMesh, tensor: Any, beta: Any, strong_reaction: Any) -> Any

Return the common positive residual time scale for spatial and time-step forms.

Conservative scalar volume forms and explicit strong-residual stabilization.

ScalarTransportOperators dataclass

Integrated A/load and physical nodal moments, with exact sampled constant-mode flags.

triangle_transport_operators(fine: TriangleMesh, *, degree: int, diffusion: Any = 1.0, velocity: Any = (0.0, 0.0), velocity_divergence: Any = 0.0, diffusion_divergence: Any = (0.0, 0.0), reaction: Any = 0.0, source: Any = 0.0, stabilization: str = 'galerkin', order: int = 6, unusual_parameters: UnusualParameters | None = None) -> ScalarTransportOperators

Integrate conservative skew transport and Galerkin, SUPG or UNUSUAL residuals.

Trial strong residual is -div(K grad(u))+beta.grad(u)+(c+div(beta))*u. SUPG tests against tau beta.grad(v); UNUSUAL subtracts tau times the complete reaction-diffusion residual pair with zero advection. The weak skew operator has c+div(beta)/2 and requires that sampled energy coefficient nonnegative. Matrices contain only volume and an explicitly supplied numerical-normal advection callback. Boundary enforcement, trace maps and kernels are caller declarations. Material-cut quadrature does not resolve interfaces for strong residual methods; those interfaces must fit the local approximation mesh.

Coefficient and inverse-inequality contracts for scalar UNUSUAL stabilization.

UnusualParameters dataclass

Declare elementwise bounds for the reaction--diffusion residual method.

diffusion_lower bounds the smallest eigenvalue of A from below and reaction_upper bounds the nonnegative reaction from above throughout each fine cell. Scalars are global bounds; callables receive fine-cell centroids and must return bounds valid on the whole corresponding cell. Smooth coefficient callbacks require these declarations. Constant and resolved Cartesian materials admit bounds obtained directly from their data.

inverse_constant optionally prescribes m in Eq. (15) of Martins, Pacazuca and Valentin (CILAMCE 2025). It must satisfy 0<m<=1/3 and the discrete elementwise coercivity check in :func:unusual_scale. Without an explicit value, a conservative operator-inverse value up to 1/3 is computed. Thus constant-A P1 uses exactly the published m=1/3. Higher-degree and variable-A constants depend on the physical element and operator, not only the degree.

__post_init__() -> None

Reject inverse parameters outside the stated admissible interval.

unusual_scale(basis: FloatArray, gradient: FloatArray, diffusion_residual: FloatArray, weights: FloatArray, diameters: FloatArray, tensor: FloatArray, reaction: FloatArray, centers: FloatArray, parameters: UnusualParameters) -> FloatArray

Compute elementwise Eq. (15), including its zero-reaction limit.

Arrays have leading axes (cell, quadrature point); gradients additionally have (basis, dimension), and diffusion_residual is div(A grad(phi)). Bounds are checked at quadrature points. The automatic inverse parameter obeys mh²||div(A grad(v))||² <= A_min²*||grad(v)||² on the local polynomial space, with both norms defined by the supplied positive integration rule. Adequate quadrature and valid bounds between samples remain caller obligations. An explicit m is accepted only when the negative residual is strictly smaller than the physical diffusion--reaction energy for every polynomial on each cell (modulo constants when reaction vanishes). This generalized eigenvalue check admits stable explicit choices beyond the conservative automatic inverse bound without claiming those bounds are identical.

Eq. (15) is evaluated as mh²/[max(sigma_maxmh²,2A_min)+2*A_min], avoiding division by zero when sigma vanishes. This is the continuous diffusion limit of that equation, with no artificial reaction floor.

Conservative tetrahedral scalar transport integration with explicit derivatives.

These numerical forms expose diffusion, skew transport, reaction and full SUPG residuals. They do not choose local kernels, macroface spaces or a solver.

coefficient_derivatives_3d(points: FloatArray, diffusion: Any, velocity: Any, velocity_divergence: Any, diffusion_divergence: Any, stabilization: str) -> tuple[Any, Any]

Require declared variable derivatives and reject contradictory constant data.

tetra_transport_operators(mesh: TetraMesh, degree: int, *, diffusion: Any = 1.0, velocity: Any = (0.0, 0.0, 0.0), reaction: Any = 0.0, source: Any = 0.0, velocity_divergence: Any = None, diffusion_divergence: Any = None, stabilization: str = 'galerkin', order: int = 6) -> tuple[Any, FloatArray, FloatArray, bool]

Assemble skew conservative RAD and optional full-residual SUPG.

The associated boundary flux is (-K grad(u)+beta*u/2).n. The source equation is -div(K grad(u))+div(beta*u)+reaction*u=f. Coercivity requires reaction+div(beta)/2 >= 0; the implementation checks integration points. Returned physical moments integrate the nodal functions over the domain. Galerkin requests first derivatives only; SUPG also tabulates its strong-residual Hessian.

Explicit one-sided coefficients associated with a macro partition.

MacroCoefficient dataclass

Explicit one-sided coefficient callbacks, one for each macroelement.

Entries follow mesh.cells. This wrapper distinguishes a tuple of local fields from an ordinary constant vector or tensor. A field is evaluated exclusively on the corresponding macroelement, including its boundary.

for_cell(cell: int, count: int) -> Any

Select a field after checking the coefficient-to-macro association.

Hydrodynamic material tensors and exact one-sided vector-field evaluation.

RT0VectorField

Evaluate one macro-local RT0 field with an exact geometric membership test.

A nearest-centroid query supplies candidates only; barycentric tests select the containing triangle and an exhaustive fallback handles skew meshes. Points on an internal fine edge choose one adjacent fine cell. Macro sides remain explicit because each instance owns only one macroelement.

__init__(mesh: TriangleMesh, flux: Any) -> None

Validate integrated face fluxes and cache affine RT0 coefficients.

locate(points: Any) -> IntArray

Locate every point, rejecting exterior or nonfinite coordinates.

__call__(points: FloatArray) -> FloatArray

Evaluate physical Darcy velocity q=a+b*x inside each fine cell.

divergence(points: FloatArray) -> FloatArray

Return the exact piecewise-constant RT0 divergence.

HydrodynamicDispersion dataclass

D=(alpha_m+alpha_t|v|)I+(alpha_l-alpha_t)v tensor v/|v|.

The final term is defined as zero at zero velocity. Molecular diffusion is strictly positive; longitudinal and transverse dispersivities are nonnegative with alpha_l>=alpha_t. Coordinates, time and coefficient units must agree.

__post_init__() -> None

Enforce ellipticity and the physical ordering of dispersivities.

__call__(points: FloatArray) -> FloatArray

Evaluate the symmetric positive-definite dispersion tensor.

Use equation (5.5) of Harder, Paredes and Valentin (2015).

divergence(points: FloatArray) -> FloatArray

Differentiate D within each fine cell, without smoothing material jumps.

Since grad(v)=b I in two dimensions, div(D)=(2 alpha_l-alpha_t)bv/|v|. Its value is set to zero at an isolated velocity zero, where |v| is not differentiable; this is an almost-everywhere derivative for broken SUPG.

PrimalDispersion dataclass

Bases: HydrodynamicDispersion

Integrate dispersion on the same material intersections as the primal Darcy solve.

triangle_quadrature(mesh: TriangleMesh, order: int) -> tuple[FloatArray, FloatArray, FloatArray]

Retain exact cut material samples and unsmoothed fine pressure derivatives.

One-sided polynomial and primal Darcy fields for passive transport.

PolynomialVectorField

Bases: _TriangleLocator

An exact cellwise Pk representation of an arbitrary physical vector field.

values contains the physical vector at equispaced Pk nodes in each fine triangle, in :func:pymhm.fem.scalar.triangle.multiindices order. For RT or BDM fields these nodal coordinates are an equivalent representation of the executed vector polynomial, rather than a projection or a smoothing operation. Degree is a positive integer. H(div) continuity and equilibrium belong to the supplied coefficients and their originating discretization; this low-level constructor does not certify either property.

__init__(mesh: TriangleMesh, values: Any, degree: int) -> None

Validate and copy physical cell coefficients, retaining independent sides.

evaluate(points: Any, cells: Any = None) -> FloatArray

Evaluate each polynomial on its declared side, without nearest-node sampling.

__call__(points: FloatArray) -> FloatArray

Evaluate containing fine cells, choosing one side at an interior edge.

gradient(points: Any, cells: Any = None) -> FloatArray

Return grad(q)[component, derivative] by exact polynomial differentiation.

divergence(points: FloatArray, cells: Any = None) -> FloatArray

Return the physical polynomial divergence, without source projection.

PrimalDarcyVelocity

Bases: _TriangleLocator

The explicit raw field -K grad(p_h), with material derivatives provided.

The field can have normal jumps inside a macrocell and is not asserted to be H(div)-conforming or fine-cell conservative. permeability_gradient has axes (point, tensor row, tensor column, derivative). For a piecewise constant material this broken derivative is zero. Material traces use the specified incident triangle, including Cartesian/planar interfaces.

Transport uses the raw vector in volume integrals and the Darcy skeletal multiplier as the numerical normal velocity. The conservative Galerkin form is integrated directly; it does not replace distributional interface derivatives by the broken divergence. The boundary term contains half this numerical normal velocity, preserving the Robin multiplier convention of Harder, Paredes and Valentin (2015). This explicitly declared discrete flux pair is distinct from an H(div) vector field. Strong-residual stabilization is not provided for this pair.

__init__(solution: DarcySolution, cell: int, permeability_gradient: Any) -> None

Retain actual primal geometry and require an explicit broken tensor derivative.

evaluate(points: Any, cells: Any = None, *, material_values: Any = None) -> FloatArray

Evaluate -K grad(p_h), optionally retaining exact material-quadrature samples.

__call__(points: FloatArray) -> FloatArray

Evaluate one containing side of the unsmoothed raw Darcy field.

gradient(points: Any, cells: Any = None) -> FloatArray

Differentiate the raw velocity, retaining all tensor derivative components.

divergence(points: FloatArray, cells: Any = None) -> FloatArray

Return broken div(-K grad(p_h)); interface distributions are separate.

normal_trace(face: int, points: FloatArray) -> FloatArray

Return the Darcy numerical normal velocity in the global macroface orientation.

advection_boundary_matrix(degree: int, order: int, *, faces: tuple[int, ...] | None = None) -> Any

Integrate half the numerical normal velocity on every macro boundary side.

The volume form is -(v_raw*u, grad(test)). This term supplies the half-advection Robin boundary contribution of Harder, Paredes and Valentin (2015). Integration splits at both fine-edge and Darcy-trace breakpoints, independently of the concentration skeleton; no pressure-gradient normal is substituted.

polynomial_darcy_velocity(solution: Any, cell: int) -> PolynomialVectorField

Convert an executed triangular RT/BDM field to its equivalent vector polynomial.

Supported solutions are RTDarcySolution, BDMDarcySolution and MomentFluxSolution. Their own field evaluators preserve canonical moments, Piola transforms and orientation. Reconstruction conservation remains exactly the conservation imposed by the chosen reconstruction; no extra constraint is introduced by this conversion.

Physical scalar transport trajectories and stated discrete mass diagnostics.

TransientTransportResult dataclass

Time history, mass diagnostics and the number of distinct offline operators.

times includes the initial time and the retained output times; solutions excludes the initial state. integration_times records the complete time grid, including steps discarded by output_steps. The balance residuals sum the original discrete physical equations on each macrocell, including essential-boundary reaction forces. They measure the discrete weak balance, not independent error or fine-cell conservation. When requested, original_residual_norms and original_rhs_norms include every executed step, in integration_times[1:] order. They check the full original physical rows and free trace equations, not just their macro sums, using the shared original-equation convention.

total_mass() -> FloatArray

Integrate capacity*u at the initial time and every retained output time.

One-sided polynomial traces and local scalar boundary integration.

PreparedScalarTrace dataclass

Read-only face geometry for repeated one-sided nodal trace evaluations.

Intervals retain the original fine-edge order and endpoint tolerance. A sorted interval index visits only overlapping pieces; the original last assignment wins at shared endpoints, including tolerance overlaps. Arrays are snapshots of the geometry, while each evaluation reads current field coefficients. No process-global geometry cache or field averaging is used.

selections(parameter: FloatArray) -> Iterator[tuple[FloatArray, IntArray, IntArray]]

Yield overlapping sample indices in the original fine-boundary order.

evaluate(coefficients: FloatArray, parameter: FloatArray) -> FloatArray

Evaluate without rebuilding edge geometry or changing endpoint ownership.

edge_basis(degree: int, parameter: FloatArray) -> FloatArray

Evaluate edge nodes ordered as endpoints followed by interior nodes.

edge_pieces(macro: TriangleMesh, fine: TriangleMesh, face: int, degree: int) -> Iterator[tuple[FloatArray, IntArray]]

Yield oriented macro parameters and nodal DOFs of matching fine edges.

prepare_scalar_trace(macro: TriangleMesh, fine: TriangleMesh, face: int, degree: int) -> PreparedScalarTrace

Prepare immutable oriented edge pieces and their exact tolerance intervals.

scalar_trace(macro: TriangleMesh, fine: TriangleMesh, face: int, degree: int, coefficients: FloatArray, parameter: FloatArray) -> FloatArray

Evaluate a continuous local trace without averaging across macrofaces.

diffusive_boundary_matrix(macro: TriangleMesh, fine: TriangleMesh, faces: tuple[int, ...], degree: int, velocity: Any, order: int) -> Any

Assemble +half beta.n boundary mass for prescribed diffusive flux.

On these exterior faces the retained multiplier is -K grad(u).n, rather than the interior Robin multiplier (-K grad(u)+beta*u/2).n.

strong_boundary_dofs(macro: TriangleMesh, fine: TriangleMesh, faces: tuple[int, ...], degree: int, value: Any) -> tuple[IntArray, FloatArray]

Collect unique local nodal Dirichlet DOFs and evaluate the boundary data.

One-sided scalar trace evaluation and reduced jump integration on macrofaces.

The caller supplies executed local responses and the explicit coefficient alpha*lower/(2H). This operation does not solve a method or choose its global Equation. Adjacent fields preserve independent one-sided traces.

FaceJumpForm dataclass

Common face quadrature, signed jumps and projected Dirichlet moments.

boundary_value(parameter: FloatArray) -> FloatArray

Evaluate the piecewise Pk orthogonal boundary projection used in enrichment.

scalar_trace_matrix(macro: Any, fine: TriangleMesh, face: int, degree: int, parameter: FloatArray, prepared: PreparedScalarTrace | None = None) -> Any

Represent one-sided local nodal traces as a sparse evaluation operator.

face_jump_form(system: HybridSystem, skeleton: SkeletonSpace, face: int, degree: int, order: int, alpha: float, lower: float, dirichlet: Any) -> FaceJumpForm

Assemble one macroface jump on the common incident fine-edge partition.

Face-residual marking and fixed-macro-mesh RAD adaptivity.

The face indicator follows Section 4 of Harder, Paredes and Valentin (2015).

TransportBounds dataclass

Declared global bounds for the transport face indicator.

The indicator follows Harder, Paredes and Valentin (2015).

diffusion bounds the largest eigenvalue of K, velocity bounds the Euclidean magnitude of beta, and effective_reaction bounds abs(c+div(beta)/2). Callers supply analytical bounds, not sampled maxima.

__post_init__() -> None

Require finite positive diffusion and nonnegative remaining bounds.

scale(skeleton: SkeletonSpace) -> float

Compute max(sqrt(Kmax), diameterbeta_max, diameter²c_eff_max).

FaceIndicators dataclass

One nonnegative indicator per skeletal segment, including exterior zeros.

The indicator follows Harder, Paredes and Valentin (2015).

total: float property

Return the Euclidean norm of the segment indicators.

mark(theta: float = 0.75) -> tuple[np.ndarray, ...]

Mark eta >= theta*max(eta), stopping for a zero indicator field.

AdaptiveTransportResult dataclass

Solved adaptive states and indicators on their actual skeletal partitions.

estimate_transport_faces(solution: ScalarSolution, bounds: TransportBounds, *, order: int | None = None) -> FaceIndicators

Integrate the skeletal jump estimator.

Use equations (4.3)--(4.4) of Harder, Paredes and Valentin (2015).

On an interior macroface F, R_F=-jump(u)/2 and each segment Fhat receives C*norm(R_F,Fhat)/sqrt(length(F)). The denominator is the original macroface length, not the refined segment length. Dirichlet conditions must be imposed strongly on the essential boundary. Prescribed natural faces are also supported: their multiplier variations vanish, so neither essential nor natural exterior faces contribute to this interior-jump indicator. solve_transport preserves the distinction between its Robin multiplier and a prescribed physical diffusive_flux; the latter has its own half-advection boundary mass. In Section 5.1 of Harder, Paredes and Valentin (2015) the velocity is tangent to the natural walls, and the two flux conventions coincide there.

This is a refinement indicator, not a computable constant-one upper bound or a maximum-principle certificate. Local-discretization errors and the approximation of nonrepresentable boundary data are not measured. A solution with weak Dirichlet imposition is rejected.

Quadrature splits at both one-sided fine-edge partitions. The returned values therefore preserve broken macro traces and exactly integrate their polynomial squared jump for the default order.

refine_skeleton_faces(skeleton: SkeletonSpace, marked: tuple[np.ndarray, ...]) -> SkeletonSpace

Bisect marked segments, preserving degrees, continuity and vector components.

solve_adaptive_transport(skeleton: SkeletonSpace, bounds: TransportBounds, *, iterations: int = 5, theta: float = 0.75, tolerance: float = 0.0, max_trace_dofs: int | None = None, solve_step: Callable[..., ScalarSolution] | None = None, **options: Any) -> AdaptiveTransportResult

Solve, estimate, mark and bisect faces while keeping macro/local meshes fixed.

iterations counts solves, including the initial space. options are arguments of solve_transport; strong Dirichlet imposition is mandatory on the essential faces. Mixed data supplied through neumann (the Robin multiplier) or diffusive_flux retain their physical meaning at every state. Exterior faces remain unmarked because their data are prescribed. The local resolution must remain adequate as the trace space grows; the shared rank diagnostics reject incompatible enrichments. Neither local refinement nor a solver tolerance is changed automatically.

solve_step(mesh, skeleton=current, **options) may consume user-declared public variational equations instead of the built-in transport solver. It must return a ScalarSolution in the supplied mesh and trace space, with the same strong essential boundary and physical Robin/diffusive-flux data. None retains the built-in solver. The indicator, marking, enrichment and stopping rules are unchanged by this assembly choice.

Broken Pk scalar field reconstructed from a hybrid system.

l2_error(exact: Any, order: int = 5) -> float

Integrate scalar error using quadrature independent of assembly.

Broken tetrahedral scalar field and full conservative physical flux.

evaluate returns u and -K grad(u)+beta*u. This physical flux differs from the skeletal Robin multiplier by beta*u/2.

evaluate(cell: int, bary: FloatArray) -> tuple[FloatArray, FloatArray]

Evaluate scalar and physical flux without averaging across macro interfaces.

l2_error(exact: Any, order: int = 7) -> float

Integrate the scalar error with a positive tetrahedral rule.

h1_seminorm_error(exact_gradient: Any, order: int = 7) -> float

Integrate the complete broken physical gradient error.

flux_l2_error(exact_flux: Any, order: int = 7) -> float

Integrate the complete conservative physical flux error.

References

• Christopher Harder, Diego Paredes, and Frédéric Valentin (2015). On a Multiscale Hybrid-Mixed Method for Advective-Reactive Dominated Problems with Heterogeneous Coefficients, Multiscale Modeling & Simulation 13(2), 491–518. DOI: 10.1137/130938499.