Solvers and execution backends
Finite-element adapters, sparse solvers, parallel execution and separable operators.
Native spaces and interface forms
Owned, lazy native bindings for user-declared finite-element spaces.
Native coefficients remain in DOLFINx order. Conversion to PyMHM's existing nodal order uses cell incidence and Basix's reference interpolation lattice, never a global nearest-neighbour search. Moment-based elements, arbitrary custom bases and mixed fields retain their native representation unless an explicitly supported nodal component is selected.
NativeSpace
dataclass
Associate a portable mesh with one owned native function-space binding.
mapping converts the existing equispaced nodal ordering to flattened
native coefficients. Mixed spaces require an explicit component, and
H(div)/H(curl) spaces use native coefficients without pretending their
moment degrees of freedom are point values. close releases this
binding's references; separately held forms and spaces remain valid.
mesh: Any
property
Return the native mesh while this binding is open.
space: Any
property
Return the native space while this binding is open.
size: int
property
Count flattened, owned native coefficients on the serial mesh.
mapping: IntArray
property
Return the supported nodal coefficient map to the native space.
close() -> None
Release native references without destroying a borrowed communicator.
to_portable(coefficients: Any, *, component: int | None = None) -> FloatArray
Delegate conversion to the shared native-to-portable coefficient owner.
to_native(coefficients: Any, *, component: int | None = None) -> FloatArray
Delegate scattering to the shared portable-to-native coefficient owner.
transport_coupling(coupling: Any, *, component: int | None = None) -> FloatArray
Delegate test-row conversion without altering mathematical trace signs.
interpolate(expression: Any) -> FloatArray
Delegate interpolation in the executed native coefficient basis.
evaluate(coefficients: Any, points: Any, *, cells: Any = None, component: int | None = None) -> FloatArray
Delegate field evaluation, retaining explicitly selected one-sided cells.
descriptor(*, component: int | None = None) -> NativeFieldDescriptor
Describe an executed nodal field without retaining native resources.
NativeFieldDescriptor
dataclass
Persist a native nodal field's executed basis and coefficient layout.
Arrays are copied and immutable. The stored coefficient map transports an executed native vector into portable nodal order before rebinding, so replay does not depend on a later DOLFINx numbering. Only equispaced nodal Lagrange components have this automatic descriptor; other elements need an explicit evaluator and declared basis contract.
__post_init__() -> None
Freeze and digest the executed basis and its parent coefficient map.
__reduce__() -> tuple[Any, tuple[Any, ...]]
Restore immutable arrays and verify basis identity when unpickling.
create_native_mesh(mesh: Any) -> Any
Create a COMM_SELF DOLFINx mesh from a supported portable mesh.
The existing portable owner supplies geometric validation and incidence. Only the documented vertex permutation to Basix's reference cell changes. The shared COMM_SELF communicator is borrowed and is never freed here.
bind_space(mesh: Any, element: Any) -> NativeSpace
Bind a user-supplied Basix/UFL element or existing native FunctionSpace.
Existing native spaces must have been created from the same input topology and coordinates. Unrelated numbering is rejected instead of guessed. Any native element may be assembled; portable coefficient conversion is separately restricted to supported equispaced Lagrange components.
coefficient_map(binding: NativeSpace, *, component: int | None = None) -> IntArray
Map portable nodal coefficients into a selected native component.
Reference nodes are matched by exact integer lattice identities from
Basix. Native cell incidence supplies global numbering, including edge
permutations; no physical coordinate search or proximity merging occurs.
Vectors use node-major interleaved components. DG fields use cell-major
local nodes. Mixed fields are selected explicitly through component.
to_portable(binding: NativeSpace, coefficients: Any, *, component: int | None = None) -> FloatArray
Copy a native coefficient vector or basis columns into portable nodal order.
to_native(binding: NativeSpace, coefficients: Any, *, component: int | None = None) -> FloatArray
Scatter portable nodal coefficients or basis columns into a zero native array.
transport_coupling(binding: NativeSpace, coupling: Any, *, component: int | None = None) -> FloatArray
Scatter portable local test rows into native rows without changing trace signs.
interpolate(binding: NativeSpace, expression: Any) -> FloatArray
Interpolate a user callback, expression or function in native coefficients.
evaluate(binding: NativeSpace, coefficients: Any, points: Any, *, cells: Any = None, component: int | None = None) -> FloatArray
Evaluate native coefficients at physical points on declared owning cells.
cells contains portable fine-cell indices, one per point. Supplying it
preserves independent one-sided traces on discontinuous interfaces. When
omitted, DOLFINx finds a containing cell; points outside the mesh fail.
Scalar values return (points,) and vectors return (points, value_size).
describe_field(binding: NativeSpace, *, component: int | None = None) -> NativeFieldDescriptor
Describe a supported nodal component using its executed Basix basis.
evaluate_descriptor(descriptor: NativeFieldDescriptor, coefficients: Any, points: Any, *, cells: Any = None) -> FloatArray
Replay coefficients using the actual archived basis polynomial matrix.
Native geometry supplies cell location and physical pullback only. Basis values use the archived matrix directly, so regenerated DOF numbering or equivalent rotations of a newly computed basis cannot alter replay.
evaluate_descriptor_data(descriptor: NativeFieldDescriptor, coefficients: Any, points: Any, *, cells: Any = None) -> tuple[FloatArray, FloatArray]
Evaluate archived values and physical gradients in one geometry traversal.
Scalar gradients have shape (points, dimension); vector gradients have
shape (points, flattened_components, dimension). Each row differentiates
that component in physical Cartesian coordinates. The native mesh's actual
Basix coordinate element supplies its Jacobian. Explicit portable cells
preserve one-sided values and derivatives on discontinuous interfaces.
Geometry and basis evaluation are shared; no process-global cache retains
native resources. Gradient construction never imposes H(div) conformity.
evaluate_descriptor_gradient(descriptor: NativeFieldDescriptor, coefficients: Any, points: Any, *, cells: Any = None) -> FloatArray
Evaluate physical component gradients using the archived executed basis.
UFL pairings in executed edge/triangle bases with explicit trace transport.
Existing FaceSpace tabulation defines every polynomial. This adapter represents those same polynomials as UFL expressions on aligned native exterior facets. It never interpolates them into the local volume unknown's function space.
trace_pairings(context: Any, native: NativeSpace, expression: Callable[[Any, Any], Any], *, axis: Literal['columns', 'rows'] = 'columns', interface: Any = None) -> LinearForms
Apply a user UFL boundary form to each existing unsigned local trace mode.
expression(phi, ds) returns a local linear form. Scalar traces provide
scalar phi; vector traces provide a component-interleaved vector.
Trial/test pairings are compiled independently. Context equation binding
subsequently applies the mesh incidence sign for normal traces exactly once.
Piecewise edge/triangle polynomials require native facets aligned to their
partitions. Shared pressure nodes, facewise C0 modes and tangential frames
retain the basis declared by the interface owner. Arbitrary custom traces
may supply their own capability explicitly.
UFL determines quadrature from the supplied form; callers may set measure
metadata when coefficients require more accurate integration.
interface optionally selects another bound trace basis on the same
macro mesh, such as a private local conormal space alongside a global
pressure trace. Its forms remain unsigned; the caller chooses their
mathematical signs and applies any declared coefficient transport.
interface_equation(context: Any, builder: Callable[[Any, Any, Any], Any]) -> Equation
Compile a user-defined global UFL form in declared edge/triangle coordinates.
builder(lambda_, mu, dx) returns Equation(a, L) or (a, L) on
an embedded interval or triangular mesh. Its discontinuous space contains
the executed face basis, including independent degrees, facewise C0 modes,
shared pressure nodes and Cartesian/tangential components. An embedding
transports both forms back to
the declared global layout. Retained coordinates receive zero added terms.
Physical normal/traction signs remain those written by the user. This
adapter supplies no cross-mesh volume/interface UFL argument coupling.
Direct form compilation
Generic optional DOLFINx assembly of user-declared UFL blocks and pairings.
Rows belong to the test arguments and columns to the trial arguments. Each block is assembled independently; symmetry, boundary conditions, skeleton orientation and trial/test compatibility are explicit caller responsibilities. The module imports no optional FEM, UFL, MPI or PETSc package at import time.
assemble_form(form: Any, *, shape: tuple[int, ...] | None = None, form_compiler_options: Mapping[str, Any] | None = None, jit_options: Mapping[str, Any] | None = None, entity_maps: EntityMaps | None = None) -> FloatArray | sparse.csr_matrix
Assemble one linear or bilinear UFL form into owned real coefficients.
A linear form returns a binary64 vector. A bilinear form returns a binary64
SciPy CSR matrix with shape (test_dofs, trial_dofs). Distinct spaces and
rectangular matrices are accepted; a later condensation or solve decides
whether its chosen pivot and full equations are algebraically admissible.
Linear forms using either argument number zero or one are supported.
shape checks the declared coefficient dimensions. Supply it for zero
forms whose arguments were simplified away. Argument-preserving UFL zero
forms can infer dimensions from their native function spaces. A zero form
returns zeros without importing DOLFINx. Rank-zero nonzero functionals and
higher-rank forms are outside this vector/matrix contract.
All form domains must use a single-rank communicator, normally COMM_SELF.
Compiler and JIT options are copied before delegation; entity_maps
explicitly describes cross-mesh integration to DOLFINx. Supply the mapping
of meshes to index arrays accepted by DOLFINx 0.9, or the sequence of native
EntityMap objects accepted by DOLFINx 0.10. Only the container is copied;
its declared indices or native maps are forwarded without conversion.
This function
infers no trace map, normal orientation, quadrature, boundary elimination,
gauge or transposed coupling. Native matrices, vectors and compiled-form
references are released before returning their independent array copies.
assemble_pairing(forms: Iterable[Any], *, axis: Literal['columns', 'rows'] = 'columns', size: int | None = None, form_compiler_options: Mapping[str, Any] | None = None, jit_options: Mapping[str, Any] | None = None, entity_maps: EntityMaps | None = None) -> FloatArray
Assemble ordered linear basis pairings as independent columns or rows.
Every argument-preserving form must use the same native function space.
axis='columns' returns (size, number_of_forms) for a local test-side
coupling. axis='rows' returns (number_of_forms, size) for a global
equation paired with the local trial field. Argument number one is accepted
natively, so a trial-side linear form needs no symbolic replacement.
size provides the local dimension for an empty sequence or an initial
zero form whose arguments were erased by UFL. Otherwise the first assembled
vector determines it. Caller-supplied signs and form order are preserved;
row and column pairings are never inferred from one another. Compilation
options and cross-mesh maps have the same meaning as in assemble_form.
Reusable local form workspaces
Compile a bundle of UFL forms once, then bind it to compatible native meshes,
spaces and coefficient data. A worker can retain its own bindings and assembly
buffers while updating geometry, Functions and Constants between local problems.
Each thread owns the mesh and coefficient data it mutates; sharing compiled code
uses explicit bindings to those separate native objects.
All bound argument and coefficient spaces belong to the workspace's integration
mesh. Explicit entity_maps use DOLFINx 0.9 mappings or 0.10 native EntityMap
sequences; this version convention preserves the same workspace mesh contract.
Each assembly returns independent arrays; material matrices, loads and solver
factors remain specific to each local problem. Native integration controls cover
single-mesh scalar, vector and mixed forms in 2D, tensor diffusion in 3D, and
constrained Neumann reconstruction in both dimensions.
Reusable, data-independent UFL kernels and bounded native assembly workspaces.
Only kernels, native form bindings and zeroed assembly buffers are reused. Every assembly returns an independent binary64 array or CSR matrix. Materials, loads, condensed responses and solver factors remain specific to each problem. Optional DOLFINx/UFL imports occur only when the native API is invoked.
Closable
Bases: Protocol
Worker-owned resource with an idempotent explicit lifetime boundary.
close() -> None
Release retained native resources.
CompiledFormBundle
dataclass
Named data-independent DOLFINx kernels, with explicit rank/shape contracts.
Keep this object alive while any workspace uses its compiled code. Original UFL symbols identify the coefficient/constant/space maps when binding to a different mesh. Coefficient values and mesh coordinates are not compiled.
FormWorkspace
dataclass
Mutable native form bindings and buffers owned by one worker thread.
Instantiate lazily inside a spawn worker, or separately in each thread. Geometry updates preserve topology, finite-element spaces and integration entities. The caller must preserve a valid mesh, including its orientation. Native Functions/Constants may also be updated directly: each assembly freshly packs their current values. No assembled result aliases buffers.
assemble(name: str) -> FloatArray | sparse.csr_matrix
Assemble a named form, preserving independent owned coefficients.
close() -> None
Release buffers, native form bindings and mesh references idempotently.
__enter__() -> FormWorkspace
Enter an explicitly scoped native workspace lifetime.
__exit__(*exc: Any) -> None
Release native resources when the scope ends, including exceptions.
__getstate__() -> Any
Prevent native handles from crossing a process serialization boundary.
WorkspaceCache
Bases: Generic[Resource]
Bounded LRU of worker-owned closable workspaces or application records.
A record may contain a FormWorkspace plus proven geometry-only data. Keys identify structural compatibility; factories run only on misses. Eviction and explicit close release the whole record. Create separate caches inside each process/thread, rather than serializing or sharing native resources.
__init__(capacity: int = 4) -> None
Create a positive-capacity cache without importing a native backend.
get(key: Hashable, factory: Callable[[], Resource]) -> Resource
Return the most recently used resource, constructing only on a miss.
close() -> None
Close every retained resource once, then reject further cache use.
__enter__() -> WorkspaceCache[Resource]
Scope the lifetime of all worker-owned cached native records.
__exit__(*exc: Any) -> None
Release all records on normal completion or an exceptional exit.
__getstate__() -> Any
Require cache creation inside the process that owns native resources.
compile_form_bundle(forms: Mapping[str, Any], comm: Any, *, form_compiler_options: Mapping[str, Any] | None = None, jit_options: Mapping[str, Any] | None = None, shapes: Mapping[str, tuple[int, ...]] | None = None) -> CompiledFormBundle
Compile named linear/bilinear UFL forms once, independently of native data.
DOLFINx 0.9 or later is required. Forms retain their quadrature metadata;
compiler/JIT options are copied. The scalar type is real binary64. Nonzero
forms each use one integration domain with binary64 native geometry. Literal
zero forms skip native compilation and require shapes if UFL erased
their arguments. A single-rank communicator is required for local assembly.
create_workspace(bundle: CompiledFormBundle, mesh: Any, *, space_map: Mapping[Any, Any] | None = None, coefficient_map: Mapping[Any, Any] | None = None, constant_map: Mapping[Any, Any] | None = None, subdomains: Mapping[Any, Any] | None = None, entity_maps: EntityMaps | None = None) -> FormWorkspace
Bind compiled kernels to a mesh, native spaces, coefficients and constants.
Coefficient/constant maps use original UFL symbols as keys. Space maps use
argument numbers (test zero, trial one), or hashable original UFL spaces.
Omitting a map binds an original native prototype to itself. For a new mesh
supply spaces and native data on that mesh explicitly. subdomains contains
DOLFINx integration-domain
pairs; if omitted, original tagged UFL integrals are packed on mesh.
Tags and integration entities must describe the same topology.
The bound mesh preserves the compiled coordinate element. Native mesh,
Function and Constant storage uses binary64 precision. All native argument
and coefficient spaces belong to this bound integration mesh.
entity_maps preserves the explicit mapping accepted by DOLFINx 0.9 or
native EntityMap sequence accepted by DOLFINx 0.10, with a copied
container. It does not relax this workspace's same-mesh space contract.
Rows follow test arguments, columns trial arguments; no gauge, sign or boundary rule
is inferred. A workspace holds no solver or material-dependent response.
assemble_workspace_form(workspace: FormWorkspace, name: str) -> FloatArray | sparse.csr_matrix
Assemble current data into zeroed reusable buffers and return an owned copy.
update_workspace(workspace: FormWorkspace, *, geometry: Any = None, coefficients: Mapping[Any, Any] | None = None, constants: Mapping[Any, Any] | None = None) -> None
Update physical coordinates and declared native data without recompilation.
Array updates preserve the native coefficient/constant shape and require real finite values. A coefficient may instead be an interpolation callable; it is evaluated after the coordinate update, on the current physical mesh. Analytical UFL expressions involving SpatialCoordinate keep their original evaluation/quadrature and read updated coordinates directly. Function ghost entries are synchronized after updates. This operation never modifies FE spaces, topology, quadrature, trace orientation or integration entities.
workspace_key(mesh: Any, forms: Mapping[str, Any], *, space_map: Mapping[Any, Any] | None = None, subdomains: Mapping[Any, Any] | None = None, entity_maps: EntityMaps | None = None, form_compiler_options: Mapping[str, Any] | None = None) -> tuple[Any, ...]
Describe kernel, FE, geometry-layout and topology compatibility, excluding data.
Form signatures include approximation elements and quadrature metadata.
Geometry coordinates are intentionally absent; their layout and cell
connectivity remain fixed. Values of coefficients/constants are absent.
Explicit entity-map keys include their declared entity correspondence;
native EntityMap keys also include dimension and both topology identities.
Custom cache keys must enforce these same structural distinctions and any
additional application-specific reusable data contracts.
Fixed hybrid local-form adapter
Optional DOLFINx/UFL assembly of user-defined local variational problems.
The adapter separates finite-element assembly from the global skeleton. It does not construct a trace space or infer boundary conditions, kernels, stable mixed pairs, or macroface orientation. Local meshes must live on one MPI rank.
from_ufl(a: Any, load: Any, trace_forms: Iterable[Any], trace_dofs: Any, *, kernel: Any = None, constraint_forms: Iterable[Any] | None = None, coarse_basis: Any = None, form_compiler_options: Mapping[str, Any] | None = None, jit_options: Mapping[str, Any] | None = None, entity_maps: EntityMaps | None = None) -> LocalProblem
Assemble UFL forms as A u + B lambda = f on a serial local mesh.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
a
|
Any
|
Bilinear UFL form with trial and test arguments in the same space. Scalar, vector, mixed, and H(div) spaces are supported by the assembler. |
required |
load
|
Any
|
Linear UFL right-hand side in that test space. |
required |
trace_forms
|
Iterable[Any]
|
One signed linear form per global trace coefficient. Column |
required |
trace_dofs
|
Any
|
Distinct global skeleton indices, one per trace form. |
required |
kernel
|
Any
|
Explicit real coefficient array shaped |
None
|
constraint_forms
|
Iterable[Any] | None
|
One linear form per retained basis vector defining physical moments C.
Omitting these uses coefficient-space orthogonality C=Z. For a retained
scalar constant use |
None
|
coarse_basis
|
Any
|
Explicit real coefficient array shaped |
None
|
form_compiler_options
|
Mapping[str, Any] | None
|
Compiler and JIT settings forwarded to each form through
:func: |
None
|
jit_options
|
Mapping[str, Any] | None
|
Compiler and JIT settings forwarded to each form through
:func: |
None
|
entity_maps
|
EntityMaps | None
|
Explicit cross-mesh integration maps forwarded to DOLFINx: mesh/index
mappings for DOLFINx 0.9 or native |
None
|
Notes
Assembly uses DOLFINx's native CSR/vector interfaces and real binary64
kernels, without a PETSc requirement. Solver selection belongs to the
existing local condensation and global algebra APIs, after assembly.
MPI.COMM_SELF is the recommended mesh communicator; any communicator
containing exactly one rank is accepted. MPI-distributed local meshes and
complex forms are rejected. Local problems may subsequently be sent to CPU
workers because their independent array copies no longer hold DOLFINx or
MPI objects. Native buffers and compiled-form references are released after
each form is assembled. Argument-preserving UFL zero forms are supported.
No essential boundary condition is applied. In an H(div) formulation, a prescribed normal velocity is an essential condition and must be expressed with an appropriate augmented form or eliminated before local condensation. Merely attaching boundary pressure forms changes the hybrid formulation.
assemble_local_forms(forms: LocalForm) -> LocalProblem
Compile a LocalForm through the native DOLFINx/UFL adapter.
Signed trace forms and literal retained coefficients are passed unchanged
to from_ufl; physical moment_forms become its constraint forms.
The existing real-valued, same trial/test space, single-rank restrictions
apply. No boundary elimination, trace space, gauge or stabilization is
inferred. Native modules are imported only when this function executes.
primal_darcy_forms(space: Any, permeability: Any, source: Any) -> tuple[Any, Any]
Return primal Darcy forms (K grad p,grad v)=(f,v).
Pressure is the primal unknown, Darcy velocity is -K grad p and the
continuity equation is div(velocity)=source. Supply a positive scalar
or positive-definite tensor permeability. Boundary flux terms are separate.
mixed_darcy_forms(space: Any, inverse_permeability: Any, source: Any) -> tuple[Any, Any]
Return symmetric mixed Darcy forms for H(div) velocity and L2 pressure.
The equations are K^-1 u + grad p=0 and div u=source. The scalar
test equation is multiplied by minus one: the bilinear form is
(K^-1 u,v)-(p,div v)-(q,div u) and the load is -(source,q).
Choose an inf-sup stable pair such as Raviart-Thomas/DG. Pressure boundary
data contribute -<p_D,v.n> to the load; normal velocity is essential.
This helper only defines volume forms and does not impose MHM trace data.
brinkman_forms(space: Any, viscosity: Any, resistance: Any, force: Any, divergence_source: Any = 0.0) -> tuple[Any, Any]
Return symmetric Stokes-Brinkman velocity/pressure volume forms.
The strong momentum operator is -div(2 mu sym(grad u))+R u+grad p and
div u=g. The pressure test equation is negated to obtain a symmetric
saddle operator. An inf-sup stable H1 velocity/pressure pair, for example
Taylor-Hood P2/P1, is required. No stabilization or boundary conditions are
inferred. resistance=0 gives Stokes; positive resistance gives Brinkman.
elasticity_forms(space: Any, lame_lambda: Any, shear_modulus: Any, force: Any) -> tuple[Any, Any]
Return small-strain isotropic elasticity displacement volume forms.
Stress is 2 mu sym(grad u)+lambda div(u) I. In two dimensions these
Lamé parameters describe plane strain; a plane-stress reduction must be
supplied explicitly by the caller. Tractions and rigid-motion constraints
are assembled separately and are not inferred from the displacement space.
usfem_brinkman_forms(space: Any, viscosity: Any, resistance: Any, force: Any, *, smallest_resistance: Any, divergence_source: Any = 0.0, inverse_estimate: float = 1 / 3) -> tuple[Any, Any]
Return the symmetric-test version of the 2017 MHM-Brinkman USFEM forms.
This formulation uses the vector Laplacian and pseudostress, with volume
form nu grad(u):grad(v)+Theta u.v-p div(v)-q div(u). Its stabilization
is negative -kappa R(u,p).R(v,q) with the matching source term
-kappa f.R(v,q), where R(u,p)=-nu div(grad(u))+Theta u+grad(p).
The pressure test has been negated relative to the paper to obtain the
symmetric saddle convention used by this package.
smallest_resistance is the nonnegative smallest eigenvalue of Theta;
it must be supplied explicitly for heterogeneous or anisotropic media.
inverse_estimate is m=min(1/3,C_inverse)>0. The default m=1/3 is valid
for the piecewise-linear pair; higher degrees require a justified inverse
estimate. With cell diameter h, kappa=h²/(max(gamma*h²,4*nu/m)+4*nu/m).
For Stokes (gamma=0, nu>0), this gives kappa=m*h²/(8*nu).
Viscosity must be constant per cell for this strong residual to represent the diffusion operator. Continuous equal-order velocity/pressure spaces are intended. Boundary pseudotractions and a global pressure gauge must be handled consistently by the caller; this helper does not impose them.
Solvers and execution
Checked sparse solvers with explicit, lazily loaded optional backends.
Every accepted solution satisfies the requested columnwise Euclidean residual criterion. Optional backends never silently fall back to a different solver. Factorizations use double precision (real or complex). Difficult residuals are accumulated with extended precision or compensated sums before acceptance.
LinearSolveError
Bases: RuntimeError
A factorization, iteration, or independently checked residual failed.
SolverUnavailableError
Bases: ImportError
An explicitly requested optional solver cannot be loaded or executed.
LinearFactorization
dataclass
Reusable fixed-matrix factorization with checked solves and explicit cleanup.
Create with :func:factorize or :func:prepare_amgx and use as a context
manager, or call close explicitly. The latter prepares an iterative AMG
hierarchy rather than direct factors. Prepared backend state must not be
shared concurrently between threads. The matrix is copied when constructed.
matches(matrix: Any) -> bool
Check exact operator identity before reusing factors for an assembled system.
solve(rhs: Any, *, refinement_precision: Literal['double', 'extended'] = 'double', refinement_steps: int = 2) -> Array
Solve and refine up to refinement_steps times using the same factorization.
Corrections use an extended-precision residual where available and
compensated row sums elsewhere. Acceptance always uses the original
columnwise criterion; invalid output and failed refinement are rejected.
refinement_precision="extended" retains correction digits in a
wider solution array, while the factorization and correction solves use
their original precision. Original wider RHS digits are retained for
residuals; only backend operands are rounded to double precision.
This explicit mode requires a wider NumPy
long-double type; no backend or tolerance is substituted.
The nonnegative correction limit defaults to two; zero checks only the
initial backend solve. Increasing it permits additional defect corrections
but does not change acceptance or guarantee convergence.
close() -> None
Release backend resources once; further solves are rejected.
__enter__() -> LinearFactorization
Enter an open factorization context.
__exit__(exc_type: type[BaseException] | None, exc_value: BaseException | None, traceback: TracebackType | None) -> None
Release resources even when a solve raises an exception.
accurate_residual(matrix: sparse.csr_matrix, rhs: Array, solution: Array) -> Array
Evaluate rhs - matrix @ solution with accurate row accumulation.
The finite CSR operator may be rectangular. rhs and solution have
matching column counts, with shapes (rows,)/(columns,) or
(rows, nrhs)/(columns, nrhs). Real and complex coefficients are
supported. Extended precision is used when the platform supplies it;
otherwise each real and imaginary row sum is compensated independently.
This operation neither rescales the equations nor applies a tolerance.
check_linear_solution(matrix: Any, rhs: Any, solution: Any, *, rtol: float = 1e-10, atol: float = 0.0) -> Array
Check an external solver result against each original right-hand side.
Require ||A x-b|| <= max(atol, rtol*||b||) separately for every column,
using the same residual evaluation as the native solver adapters. This
function does not correct a result or select another solver. The operator
is represented in real or complex binary64; explicitly wider right-hand
sides and solution digits remain available to residual evaluation.
Passing this algebraic check does not establish discretization stability
or accuracy of physical fields.
preload_solver_backend(solver: str) -> None
Load only the selected solver's native libraries before thread limits.
This idempotent preparation constructs no explicit factors, AMG hierarchy,
communicator or GPU session. Call before a new threadpool_limits context so newly
loaded BLAS/OpenMP libraries receive the same policy during assembly and
solution. Portable SciPy solvers require no additional imports. Optional
dependencies remain lazy until explicitly selected; unavailable backends
raise the same diagnostic as the numerical adapters. Native parameter
profiles, precision and residual tolerances are unchanged. PETSc preparation
imports only its Python package; PETSc/MPI initialization remains owned by
the existing numerical adapter.
validate_invertible(matrix: Any) -> None
Check finite square input and numerical rank using equilibrated sparse LU.
This diagnostic rejects zero rows/columns, failed LU, or a pivot no larger
than eps * n * max(abs(U)) after row/column max-norm equilibration. It is
a floating-point rank diagnostic, not proof of exact singularity or a bound
on the condition number. In MHM it detects unsupported trace enrichment or
missing constraints even when a compatible RHS has a small residual.
The factorization is discarded; ordinary SciPy solves perform the same check
while reusing their LU. Use this separately for other global backends when
a backend-independent stability check is needed.
symmetric_equilibration(matrix: sparse.csr_matrix) -> tuple[sparse.csr_matrix, Array]
Return five infinity-row Ruiz congruences and their positive diagonal.
The relation is balanced = D @ matrix @ D. Symmetry is preserved;
neither rows nor equations are discarded, and the diagonal is not a
regularization. Input is a finite square CSR operator. A zero row is
incompatible with invertible direct solves. The returned vector contains
the diagonal of D; transform a right-hand side by D @ rhs and
recover a physical solution by D @ balanced_solution. This scaling
does not establish positive definiteness or remove a physical kernel.
factorize(matrix: Any, *, solver: str = 'scipy', rtol: float = 1e-10, atol: float = 0.0, equilibration: Literal['none', 'symmetric'] = 'none', rhs_columns: int = 1) -> LinearFactorization
Factor a fixed matrix once for repeated, independently checked solves.
scipy uses SuperLU; pypardiso uses MKL PARDISO (real matrices only).
pypardiso-symmetric explicitly selects its symmetric-indefinite factor
(mtype=-2), with an upper-triangle representation and stored zero diagonals.
pypardiso-symmetric-matching additionally selects the documented explicit
weighted-matching/pivot profile; neither profile substitutes for the other.
Symmetry is checked; residuals and corrections always use the full original
matrix, including any assembly roundoff in its lower triangle.
petsc requires MUMPS and uses its pivoted LU on PETSc COMM_SELF;
petsc-symmetric uses its real symmetric-indefinite pivoted LDLt through
PETSc's Cholesky interface with the SPD flag explicitly false.
cudss uses NVIDIA cuDSS through
nvmath and CuPy. Native resources are released by close or a context
manager. CuPy's QR and Krylov methods are only available in solve_linear.
No backend is substituted when the requested one is unavailable. SciPy LU
equilibrates rows/columns and checks numerical pivots before solving.
equilibration="symmetric" applies five infinity-row Ruiz congruences
before the selected factorization: A' = D A D, b' = D b, x = D y.
It requires a Hermitian operator, preserves the original residual criterion,
and is independent of the backend's own scaling. The default is "none".
rhs_columns is the positive cuDSS solve width. It permits one native
solve/transfer for all local response columns, instead of one per basis.
Other backends accept the hint without changing their multi-RHS behavior.
Subsequent cuDSS solves of another width reuse the same factors in chunks,
padding a final short chunk with zero loads. CUDA factors and cleanup stay
on their construction device and restore the caller's current device.
cuDSS host API entries are serialized within the process because its
analysis and factorization phases do not guarantee host thread safety.
Device transfers and asynchronous GPU work do not hold that entry lock.
PARDISO construction, factorization, solve and cleanup entries are also
serialized per process; each factor retains its independent lifetime.
require_hermitian(matrix: sparse.csr_matrix, operation: str = 'operation') -> None
Require a square CSR operator to be Hermitian within relative roundoff.
The maximum entry of A - A.conjugate().T must not exceed 1e-13
times the maximum absolute entry of A. Real symmetry uses the same
convention. operation identifies the caller in a rejection message.
Passing this check establishes neither invertibility nor positivity.
prepare_amgx(matrix: Any, *, rtol: float = 1e-10, atol: float = 0.0, maxiter: int | None = None) -> LinearFactorization
Prepare one GPU FGMRES/AMG hierarchy for repeated checked elliptic solves.
The real float64 operator is copied. It must be symmetric with positive diagonal entries; these checks do not establish positive definiteness. Each nonzero relative-tolerance RHS is divided by its largest absolute entry, and the returned solution is multiplied by the same value. This linear change of scale keeps tiny loads above AmgX's native absolute stopping floor. Original RHS columns and tolerances define acceptance.
The returned :class:LinearFactorization reuses this hierarchy, native
buffers and transfers for all RHS columns and true-residual corrections.
It owns process-global AmgX initialize/finalize and holds the package's
lifecycle lock until close; use a context manager and never overlap an
external pyamgx session in the same process. Separate spawned processes
can select separate GPUs before importing CUDA libraries. No numerical
matrix, response or hierarchy is reused for a different operator.
solve_linear(matrix: Any, rhs: Any, *, solver: str = 'scipy', rtol: float = 1e-10, atol: float = 0.0, maxiter: int | None = None, near_nullspace: Any = None, refinement_precision: Literal['double', 'extended'] = 'double', refinement_steps: int = 2, equilibration: Literal['none', 'symmetric'] = 'none') -> Array
Solve matrix @ x = rhs and require ||r|| <= max(atol, rtol*||b||).
The criterion is checked separately for every column. Inputs must be finite
and the square matrix must be invertible for direct solvers. Available
methods are scipy, cg, minres, gmres, pypardiso,
pypardiso-symmetric and pypardiso-symmetric-matching (real symmetric
indefinite), petsc (MUMPS LU), petsc-symmetric (real MUMPS LDLt),
cupy (GPU sparse QR), cudss, pyamg (SA-preconditioned CPU CG),
and amgx (GPU FGMRES with aggregation AMG). CG requires a Hermitian positive
definite matrix; MINRES requires real symmetry, permits indefiniteness, and
may solve compatible singular systems. GMRES accepts nonsymmetric systems.
Symmetry is checked for CG/MINRES; positive definiteness is the caller's
responsibility. Optional backends require a working native runtime.
The two AMG presets target positive-definite elliptic operators, not the
complete indefinite MHM saddle system. Symmetry and positive diagonals are
checked; these checks do not establish positive definiteness. PyAMG accepts
optional near_nullspace candidate columns (n,k), such as rigid modes for
constrained elasticity. Other backends reject this argument. The AmgX
adapter owns initialize/finalize and must not overlap an external pyamgx
session in the same process. Its native resource lifecycles are serialized.
GPU transfers are included in this call. PETSc requires MUMPS and uses
COMM_SELF, not a distributed matrix. Use :func:factorize to amortize
repeated direct solves. Direct backends accept equilibration="symmetric"
with the congruence and original-residual contract of :func:factorize.
Krylov methods may correct up to refinement_steps prematurely accepted iterates using
independently accumulated true residuals, reusing the preconditioner.
Reusable direct factors use the same nonnegative correction limit, defaulting
to two. Zero checks only the initial solve; extra corrections never relax the
requested criterion. AmgX reuses one prepared hierarchy for all columns and
corrections. Its relative-tolerance operands are normalized per column to
avoid an absolute native stopping floor on tiny loads. CuPy QR does not
expose this correction loop and rejects nondefault values.
refinement_precision="extended" retains the corrected solution in
extended precision while matrix storage and correction solves remain in
double precision. This explicit mixed-precision mode supports Krylov
methods and reusable direct factors (SciPy, PARDISO, PETSc and cuDSS); it
requires a wider NumPy long-double type. It is useful
when the requested residual lies below the double-precision solution's
rounding floor. No tolerance is changed, and unsupported platforms reject it.
Time updates of second-order coefficient equations without physical dispatch.
newmark_step(mass: Any, stiffness: Any, step_factor: LinearFactorization, mass_factor: LinearFactorization, duration: float, u: Array, v: Array, old: Array, new: Array) -> tuple[Array, Array]
Advance M u''+K u=f with beta=1/4 and gamma=1/2.
duration is a positive local step. u, v and the endpoint
forces old/new are matching vectors or matrices of independent
response columns, in one declared coefficient basis. The caller provides
factors of mass+duration**2/4*stiffness and mass and retains their
lifetime; this operation creates or closes no native resource. Real and
complex operators follow the same algebra. No material, trace or boundary
convention is inferred from these arrays.
The two checked factor solves preserve the coefficient operation order used for both local histories and slabwise constant trace response columns.
Exact real coordinates for user-defined complex coefficient equations.
The scalar ordering is (Re z0, Im z0, Re z1, Im z1, ...). No physical
operator, conjugation convention, damping or boundary condition is inferred.
realify_vector(value: Any) -> FloatArray
Interleave a finite complex vector's real and imaginary coordinates.
Real input has zero imaginary coordinates. The input must be rank one; the returned real array owns its data and preserves floating precision. This map represents coefficients rather than changing the definition of a complex variational form.
complexify_vector(value: Any) -> Any
Decode an even-length finite real vector in interleaved Re/Im order.
realify_operator(value: Any) -> Any
Represent complex multiplication by blocks [[Re,-Im],[Im,Re]].
Rectangular dense or sparse operators are accepted. The returned CSC
matrix has twice the row/column count and acts on :func:realify_vector
coordinates. For a sesquilinear form the caller first assembles its complex
matrix with the desired conjugated test convention; this function embeds
that matrix exactly and does not transpose or conjugate it.
Sparse scatter of user-integrated element forms in declared coefficient maps.
assemble_element_blocks(blocks: Any, rows: Any, columns: Any, shape: tuple[int, int]) -> Any
Scatter element matrices to a CSC operator without selecting a PDE.
blocks[e,i,j] pairs test coefficient rows[e,i] with trial
coefficient columns[e,j]. Repeated global coordinates sum, including
repeated coordinates within an element. Square matrices declare identical
trial/test maps and equal dimensions explicitly. Real/complex entries
keep their numeric dtype. The caller owns integration, orientation signs
and local basis transformations; none are inferred during this scatter.
Discrete point-source functionals with geometric sharing at mesh interfaces.
split_point_sources(mesh: TriangleMesh | CartesianMacroMesh | PolygonMesh, sources: Any) -> tuple[FloatArray, ...]
Share source strengths among cells in proportion to their incident angles.
Triangular, Cartesian and simple polygonal meshes are supported. Interior points belong to one cell; interface points are split by the incident angular sectors. At an interior edge this gives half to either cell, and at a vertex it uses the actual triangle angles. This is the limit of a radially symmetric regularization of the load, restricted to the domain. The total specified strength is preserved, including at domain corners, without duplicating wells shared by several macrocells.
point_load_vector(mesh: TriangleMesh | CartesianMacroMesh, degree: int, sources: Any) -> FloatArray
Assemble sum Q_j v(x_j) for continuous triangular Pk or Cartesian Qk tests.
A source on an internal fine edge/vertex is evaluated once; continuity of
the nodal finite-element test space makes the containing-cell choice
immaterial. Use :func:split_point_sources first to allocate sources across
a broken macro space. Two-dimensional Dirac sources are not H⁻¹ loads;
this discrete operation does not imply a finite continuum energy norm.
assemble_load(tests: Any, samples: Any, weights: Any, rows: Any, size: int) -> Any
Scatter sum_q weights[t,q]*tests[t,q,i]*samples[t,q] into test rows.
Tests have axes (cell, point, local test). Physical Jacobians and any conjugation are explicit parts of the supplied tests or weights; none is inferred. Samples/weights broadcast to (cell, point). Real and complex coefficients are supported. Continuous and discontinuous spaces differ only in the integer row map, whose repeated coordinates are summed.
assemble_mass(tests: Any, trials: Any, weights: Any, rows: Any, columns: Any, shape: tuple[int, int]) -> Any
Scatter the literal test/trial pairing with physical weighted measure.
Tabs have axes (cell, point, local basis), with possibly different trial
and test counts. weights broadcasts to (cell, point) and may contain
density, capacity or signed/Petrov coefficients. This operation makes no
positivity, Hermitian or time-scheme assumption. A sesquilinear form
supplies conjugated test values explicitly.
Ordered execution of independent local problems on portable CPU workers.
ExecutionConfig
dataclass
Runtime policy for ordered independent local computations.
This configuration changes scheduling and supported BLAS/OpenMP thread
counts only; it does not select numerical dtypes, operators or tolerances.
None workers use the executor's platform-appropriate CPU default.
None batch size uses that worker count in parallel and one in serial.
An explicit batch size bounds the number of items submitted together;
serial execution always consumes and solves one item at a time.
pipeline=True replaces parallel batch barriers with a bounded ordered
window; the default preserves atomic batch delivery and native limits.
native_threads=None preserves the existing native thread settings.
effective_workers: int
property
Return the actual pool size, including the Windows process limit.
effective_batch_size: int
property
Return the input bound; serial execution always has bound one.
__post_init__() -> None
Reject invalid scheduling limits before any input is consumed.
iter_local(function: Callable[[Input], Output], items: Iterable[Input], *, backend: Literal['serial', 'thread', 'process'] = 'serial', workers: int | None = None, native_threads: int | None = 1, batch_size: int | None = None, pipeline: bool = False) -> Generator[Output, None, None]
Yield independent local results in input order with bounded consumption.
Validation is immediate; computation starts on the first iteration.
Serial execution solves and yields one item before consuming the next.
Parallel execution consumes at most one batch, completes it, and yields its
results before consuming another batch. batch_size=None uses the
effective worker count. Parallel failures propagate before that batch is
yielded; already yielded earlier batches remain available to the caller.
With pipeline=True, parallel execution instead maintains at most
batch_size consumed inputs whose outputs have not yet been yielded.
It awaits only the next input's result and refills one slot after the
consumer resumes. Other workers can continue while the consumer reduces
each result. Successful earlier items may be yielded before a later input
or worker failure is observed; result order is always input order.
Serial execution is unchanged by pipeline.
Call close() or use contextlib.closing when ending iteration early
to join workers and release executor resources. The input iterable remains
caller-owned. Running tasks finish before closure returns; pending tasks
are cancelled after a failure.
Threads avoid serialization and suit native kernels that release the GIL.
A callable may implement an idempotent prepare_runtime() to load its
selected libraries before native thread limits and numerical work. It runs
once per serial/thread execution or spawned worker and creates no factors
or native sessions. Calls without this hook retain ordinary lazy execution.
A callable may implement close() to release its resident resources.
Serial/thread execution invokes it after all running jobs finish, including
failure or early iterator closure. Process execution invokes it on each
spawned copy at normal worker shutdown; the caller's copy is not closed.
Mutable native workspaces must be independent for each concurrent thread.
Processes use spawn on every platform: functions and inputs must be
picklable and callers must protect executable scripts with
if __name__ == '__main__':. No live PETSc, CUDA, or factorization object
should be passed between processes; construct these inside the worker.
The factory is serialized once per worker; individual tasks transfer only
their item. Large read-only mesh metadata in a factory is therefore not
repeatedly sent with every macrocell.
native_threads limits supported BLAS/OpenMP pools (default one) to avoid
worker/native-thread oversubscription. None preserves native settings.
In thread mode that limit is process-wide while a batch executes; do not
overlap calls that request different limits. The limit is restored before
yielding results in the default batch mode. Thread pipeline mode retains
the process-wide limit across consumer code until exhaustion, closure or
failure, and restores it only after workers join. Process pipeline mode
preserves the coordinator's native limits. Serial execution limits each
individual call.
map_local(function: Callable[[Input], Output], items: Iterable[Input], *, backend: Literal['serial', 'thread', 'process'] = 'serial', workers: int | None = None, native_threads: int | None = 1, batch_size: int | None = None, pipeline: bool = False) -> list[Output]
Collect :func:iter_local in order, preserving its execution contract.
The returned list holds all results; only input consumption and submitted
tasks are bounded by batch_size, with ordered rolling scheduling when
pipeline=True. Use :func:iter_local to release
results incrementally. Numerical values and dtypes are passed unchanged.
Distributed MHM assembly and direct solution on an MPI communicator.
Each rank owns its input cells, local matrices, responses and reconstructed fields. PETSc distributes the global matrix and vector by rows; off-rank contributions are assembled collectively. No rank gathers the global matrix or all local responses. MUMPS performs distributed pivoted LU of the saddle system.
DistributedHybridSolution
dataclass
Owned global coefficients and fields of the cells assigned to this rank.
ownership is PETSc's half-open row interval, including retained modes and
any gauge multipliers. local_trace contains exactly the skeleton values
needed by local cells. The residual uses the collective original physical
equations after excluding prescribed rows; gauges cannot repair incompatible
physical data beyond the componentwise evaluation bound for a declared
kernel. General retained modes receive no such allowance. linear_residual
also checks the augmented boundary system.
local_bases stores each owned cell's executed retained basis in the same
order as coarse and fields; persist it with coefficient vectors.
solve_distributed(factory: Callable[[Any], LocalProblem | LocalAssembly | LocalEquations | CompiledLocalEquations], local_items: Iterable[Any], *, trace_size: int, comm: Any, boundary_load: tuple[Any, Any] | None = None, fixed: dict[int, float] | None = None, moments: Sequence[tuple[Sequence[FloatArray], float]] | None = None, compiler: FormCompiler = compile_form, global_equation: Equation | None = None, local_solver: str = 'scipy', rtol: float = 1e-10, atol: float = 0.0) -> DistributedHybridSolution
Assemble and solve MHM collectively from nonreplicated rank-owned cells.
comm is an mpi4py communicator. Invoke this function on every rank,
including ranks with no cells. trace_size and fixed are replicated
small metadata and must agree between ranks. Skeleton numbering must cover
range(trace_size). Retained modes are numbered after traces, first by
rank, then by local input order. A distributed mesh partitioner is not
inferred: the application assigns each physical cell to exactly one rank.
LocalEquations use the common compiler and contribution owner, including
explicit D/g and distinct trial/test trace maps. global_equation supplies
this rank's additive forms in the physical trace/retained coordinates, before
gauge rows; replicated nonzero forms would be counted once per rank. Sparse
forms stay sparse through insertion. Nested MultiscaleProblem equations are
currently unsupported. Native local UFL meshes normally use COMM_SELF.
boundary_load=(indices, values) is an additive owned boundary-moment
contribution, not a replicated full load. Contributions are subtracted as in
HybridSystem. Every physical moment supplies this rank's weight vectors
and the same global target on all ranks. Fixed values are trace coefficients;
they are eliminated through PETSc row/column elimination with RHS lifting.
PETSc/MUMPS is mandatory; there is no sequential or iterative fallback. The numerical-library factorization detects singular pivots and the original equations are independently checked. This does not prove an inf-sup bound. Coefficients must be real and exactly representable in binary64; this backend does not retain wider real correction digits or complex coefficients.
Resident, batched double-precision local algebra for CUDA devices.
Dense batching is intended for many small local systems of equal size. Sparse finite-element assembly remains a separate operation; this module accepts both host arrays and matrices assembled directly as CuPy arrays. Factorizations, right-hand sides, residual checks and returned solutions stay on the device unless the caller explicitly requests host arrays.
BatchedFactorization
Reusable resident LU for real arrays of shape (batch, n, n).
The matrices are copied and equilibrated by rows and columns before LU. Residual acceptance is checked in original units, separately for every matrix and RHS column. Refinement reuses the same factors. This is dense algebra; quadratic storage and cubic local factorization limit useful sizes. A single object belongs to the CUDA device active at construction. Operations enter that device and restore the caller's context. Host inputs are copied to the owning device; resident arrays from another device are rejected. The object must not be shared concurrently between threads or processes.
__init__(matrices: Any, *, rtol: float = 1e-10, atol: float = 0.0) -> None
Upload or copy one equal-size matrix batch and factor it once.
solve(rhs: Any, *, host: bool = False) -> Any
Solve (batch,n,nrhs) loads; return resident arrays unless host=True.
close() -> None
Synchronize the owner and release arrays without clearing its application pool.
__enter__() -> BatchedFactorization
Enter an open device-factorization context.
__exit__(exc_type: type[BaseException] | None, exc_value: BaseException | None, traceback: TracebackType | None) -> None
Release this object's device arrays on normal and exceptional exit.
assemble_p1_batch(points: Any, cells: Any, *, diffusion: Any = 1.0, source: Any = 0.0) -> tuple[Any, Any, Any]
Assemble dense P1 simplex operators entirely on the current CUDA device.
points has shape (batch,npoints,d), d=2 or 3; all local meshes share
the integer connectivity cells[ncells,d+1]. Diffusion is a positive scalar
or SPD tensors broadcastable to (batch,ncells,d,d). Source values are
constant per fine cell, broadcastable to (batch,ncells). The returned
stiffness, consistent mass and source vectors are resident CuPy arrays.
Geometry, Jacobians, physical gradients, material contractions and assembly
scatter-add all run on the GPU. No coefficient callback or curved geometry
is inferred. Dense global-in-each-local-mesh storage limits useful sizes.
condense_batched(problems: Sequence[LocalProblem]) -> tuple[LocalResponse, ...]
Condense homogeneous-size groups on the GPU, preserving original order.
Matrices originate from each problem's shared condensation contract, including
nonsymmetric test/trial spaces. Grouping uses augmented size and RHS count.
Host conversion occurs once per group after all lifts have been computed;
the returned responses interoperate with the ordinary global assembler.
For repeated resident algebra use BatchedFactorization directly.
condense_multi_gpu(problems: Iterable[LocalProblem], *, devices: Sequence[int], batch_size: int = 8, solver: Literal['auto', 'batched', 'cudss'] = 'auto', dense_size_limit: int = 512) -> tuple[LocalResponse, ...]
Condense generic host local problems on explicitly assigned CUDA devices.
One dedicated thread owns each unique device. At most one batch per device is submitted; input is consumed incrementally and responses retain original order. Returned responses remain host data suitable for the ordinary global assembler. Local operators, moments, test/trial spaces and residual tolerances come from the shared condensation contract; no PDE or trace is inferred.
auto uses dense batched LU up to dense_size_limit augmented rows and
sparse cuDSS above that limit. Explicit batched and cudss select one
path for all nonempty systems. Dense groups reuse :func:condense_batched;
sparse systems reuse the shared checked factorization owner. Host/device
transfers, native setup, stream synchronization and executor shutdown finish
before return. Factors are released within the owning device context on
normal and exceptional exit. Resident cross-device arrays are not accepted
by the batched factorization; this interface receives host LocalProblem
objects. The bound controls submitted work, not total stored responses.
AMG block preconditioning of real symmetric mixed saddle systems.
SaddleBlockSolver
Reusable block-diagonal preconditioned GMRES for [A B; B.T -C].
The original matrix is never shifted. AMG acts on the positive augmented
preconditioner block P=A+alpha*B@B.T; the second block is
S=C+B.T@diag(P)^-1@B and uses a sparse direct factorization. Both blocks
must be positive definite. This supports a semidefinite A when the coupling
controls its kernel. It is not an arbitrary indefinite-matrix AMG adapter.
backend='pyamg' uses a retained CPU smoothed-aggregation V-cycle;
backend='amgx' retains a GPU aggregation hierarchy, transferring each
block vector. The outer Krylov matrix and coarse Schur block remain on CPU.
An AmgX context exclusively owns its process-global native lifecycle.
__init__(matrix: Any, split: int, *, backend: str = 'pyamg', augmentation: float | None = None, rtol: float = 1e-10, atol: float = 0.0, maxiter: int = 200, restart: int = 50) -> None
Prepare positive block hierarchies for an explicit contiguous primal/coarse split.
solve(rhs: Any) -> np.ndarray
Solve every RHS against the original saddle operator with strict residual checks.
as_factorization() -> LinearFactorization
Expose the prepared solve through HybridSystem's exact-matrix reuse contract.
close() -> None
Release retained coarse/GPU resources exactly once.
__enter__() -> SaddleBlockSolver
Enter an open block solver context.
__exit__(exc_type: type[BaseException] | None, exc_value: BaseException | None, traceback: TracebackType | None) -> None
Release all owned resources on normal or exceptional exit.
Matrix-free Cartesian Qk diffusion with a low-order refined AMG preconditioner.
The solved operator retains the explicitly separated coefficient and the Qk Gauss rule. A Q1 problem on the same equidistant nodal grid supplies only the preconditioner: it never replaces the high-order bilinear form or right-hand side.
TensorDiffusionOperator
dataclass
Store one-dimensional Galerkin factors and apply their exact Kronecker sum.
Terms are (y_factor, x_factor) and vectors use x-fast ordering. Restricting
both one-dimensional factors to their interior nodes imposes homogeneous
strong Dirichlet conditions without constructing a two-dimensional matrix.
apply(values: FloatArray) -> FloatArray
Evaluate the matrix action, preserving an explicitly wider input dtype.
interior() -> TensorDiffusionOperator
Remove the exterior nodal lines while retaining the physical load ordering.
assemble() -> sparse.csr_matrix
Materialize the same operator when a sparse matrix is explicitly required.
SeparableKrylovSolution
dataclass
Bases: ConformingQuadrilateralSolution
Conforming solution with the measured iteration count and true residual criterion.
tensor_diffusion_operator(counts: tuple[int, int], bounds: tuple[float, float, float, float], degree: int, *, permeability: SeparableField, source: SeparableField, order: int) -> TensorDiffusionOperator
Construct Qk factors without allocating mesh connectivity or a global matrix.
counts are the two cell counts and bounds are (x0,x1,y0,y1).
Positive scalar permeability is checked at the declared tensor Gauss points.
solve_separable_krylov(mesh: CartesianMacroMesh, *, degree: int, permeability: SeparableField, source: SeparableField, dirichlet: Any = 0.0, quadrature_order: int = 6, preconditioner_order: int = 2, rtol: float = 1e-10, maxiter: int = 1000, refinement_precision: Literal['double', 'extended'] = 'double') -> SeparableKrylovSolution
Solve the original Qk operator with CG and a Q1 refined-grid AMG V-cycle.
Both Galerkin matrices are symmetric positive definite after full exterior
Dirichlet elimination, for positive sampled permeability. Symmetric AMG
smoothing defines the CG preconditioner. This does not assert an arbitrary-
contrast or arbitrary-degree condition-number bound. The independently
accumulated true residual must satisfy ||b-Au|| <= rtol*||b||. Up to two
correction solves reuse the hierarchy; the explicit extended mode retains
correction digits in a wider long-double type and never changes this test.