API reference

module MultiGridBarrier

MultiGridBarrier solves nonlinear convex optimization problems in function spaces using a barrier (interior-point) method accelerated by a multigrid hierarchy. The package exposes:

• Single-level mesh constructors: fem1d, fem2d, fem2d_P1, fem2d_P2, fem3d, spectral1d, spectral2d. Each returns a Geometry. fem1d/fem2d also build embedded manifolds — a curve or surface in a higher ambient dimension — via ambient=Val(e) (intrinsic gradient operators, surface/arc-length measure); see the manual's "Embedded manifolds".
• Topological connectivity for slit domains / glued manifolds: tensor_dofmap builds full-node connectivity from corner connectivity (no coordinates); pass the result as the t= keyword of fem1d/fem2d/fem3d so geometrically-coincident nodes stay distinct.
• The hierarchy builder: amg(geom) wraps a Geometry and returns a MultiGrid (algebraic-multigrid hierarchy on the fine mesh).
• A mesh-refinement utility: subdivide(geom, L) returns a refined Geometry. Compose with AMG via amg(subdivide(geom, L)).
• A legacy geometric-subdivision hierarchy: geometric_mg(geom, L). Still available for callers that specifically want geometric transfers; new code should prefer amg.
• Problem assembly: assemble(mg::MultiGrid; kwargs...) -> MGBProblem.
• The main solver: mgb_solve(prob::MGBProblem; kwargs...) -> MGBSOL.
• A time-dependent solver: parabolic_solve(mg::MultiGrid; kwargs...).

A gentle introduction via the p-Laplacian

For a domain Ω ⊂ ℝᵈ and p ≥ 1, consider the variational problem

\[\min_{u} \; J(u) = \int_{\Omega} \|\nabla u\|_2^p + f\,u \, dx\]

subject to appropriate boundary conditions (e.g., homogeneous Dirichlet). The Euler–Lagrange equation gives the p-Laplace PDE.

Typical workflow

using MultiGridBarrier, PyPlot          # PyPlot enables the plotting extension
geom = fem2d_P2()                       # single-level mesh
mg   = amg(geom)                        # build AMG hierarchy
prob = assemble(mg; p = 1.5)            # assemble the (native) MGBProblem
sol  = mgb_solve(prob)
plot(sol)

CUDA GPU acceleration (optional extension)

CUDA support is provided as a Julia package extension (MultiGridBarrierCUDAExt). Load CUDA and CUDSS_jll to enable it, then select the GPU with device = CUDADevice:

using CUDA, CUDSS_jll
using MultiGridBarrier
sol = mgb_solve(assemble(amg(fem2d_P2()); p = 1.5); device = CUDADevice)

mgb_solve moves the assembled problem to the GPU, solves there, and moves the solution back, so the returned MGBSOL is always in native CPU types. When a functional GPU is present the default device becomes CUDADevice; pass device = CPUDevice to force the CPU. The lower-level native_to_cuda / cuda_to_native converters remain available.

See also

• Mesh constructors: fem1d, fem2d, fem2d_P1, fem2d_P2, fem3d, spectral1d, spectral2d; tensor_dofmap (+ the t= keyword) for slit-domain / glued-manifold connectivity.
• Hierarchy: amg, subdivide (and the legacy geometric_mg).
• Solvers: mgb_solve, parabolic_solve.
• Convex: convex_Euclidian_power, convex_linear, convex_piecewise, intersect.
• Visualization & sampling: plot (via the PyPlot extension: using MultiGridBarrier, PyPlot), interpolate. The solver core has no Python dependency; plotting and the pyamg prolongators are opt-in extensions.
• Backend selection: device kwarg of mgb_solve, CPUDevice, CUDADevice, native_to_device, device_to_native.
• Problem assembly: assemble, MGBProblem.
• CUDA (lower level): native_to_cuda, cuda_to_native, clear_cudss_cache!.

FEM2D_P1{T}

Discretization tag for 2D P1 triangular FEM. Stores the geometry's own per-element corner tensor K of shape (3, N, 2) — K[v, e, d] is coordinate d of vertex v of triangle e, equal to the enclosing Geometry's mesh tensor x, so after subdivide it is the fine mesh, not the original input. Informational: the hierarchy builders work from the Geometry's x and connectivity t.

FEM2D_P2{T,B}

2D simplicial P2 FEM discretization descriptor. The Bool parameter B selects the space: B = true is P2 enriched with a per-element cubic bubble (7 nodes per triangle; the historical default, whose nodal quadrature weights are strictly positive), B = false is plain Lagrange P2 (6 nodes per triangle; the nodal quadrature is the edge-midpoint rule, whose corner weights are zero). FEM2D_P2{T} (either variant) is the family; concrete descriptors are FEM2D_P2{T,true} / FEM2D_P2{T,false}.

Stores the geometry's own full node mesh Kfull of shape (V, N, 2), V = B ? 7 : 6 — equal to the enclosing Geometry's mesh tensor x, so after subdivide it is the fine mesh, not the original input — and the corner triangulation K of shape (3, N, 2) extracted from it (local slots 1, 3, 5). Informational: the hierarchy builders work from the Geometry's x and connectivity t. The per-triangle node layout is corner1, edge(1,2), corner2, edge(2,3), corner3, edge(3,1)[, centroid].

For the pure-P2 variant the zero corner weights mean a variational problem's slack must not live in the fully broken :full space: a corner slack value would carry no cost and the central path escapes along it. amg and geometric_mg therefore provide the :broken_P1 subspace (per-element linear functions, 3 DOFs per triangle, parametrized by the edge-midpoint values), which assemble uses as the default :s space on pure-P2 geometries: every :broken_P1 direction moves the positively weighted midpoint values, so the slack is determined.

Geometry{T,X<:AbstractArray{T,3},W,M_op,Discretization}

Single-level container for discretization geometry. Holds only the fine-level mesh and operators — no multigrid hierarchy. Use amg(geom) to attach an algebraic-multigrid hierarchy and return a MultiGrid. (The legacy geometric_mg(geom, L) builds a geometric-subdivision hierarchy instead; new code should prefer amg.)

Type parameters

• T: scalar numeric type (e.g. Float64)
• X<:AbstractArray{T,3}: type of the mesh tensor x (typically Array{T,3}).
• W: type of the weight storage w (typically Vector{T}).
• M_op: matrix type for operators (e.g. SparseMatrixCSC{T,Int}, BlockDiag{T}).
• Discretization: front-end descriptor (e.g. FEM1D{e,T}, FEM2D_P2{T,B}, SPECTRAL1D{T}).

Fields

• discretization::Discretization: discretization descriptor encoding dimension and grid info.
• t::Matrix{Int}: cached element connectivity of shape (V, N) — t[v, e] is the global node id of local node v in element e, so broken nodes that coincide across elements share an id and maximum(t) is the number of distinct nodes. The 4-arg constructor derives t from x by geometric dedup; the inner 5-arg constructor (and fem1d/fem2d/fem3d with t=) take it verbatim, letting coincident-but-distinct nodes (slits, branch cuts, glued manifolds) stay topologically separate. The amg path (find_boundary, amg) reads t directly; coordinate dedup builds t only in the 4-arg constructor and the legacy geometric_mg path. Always a CPU Matrix{Int}, even for GPU geometries.
• x::X: mesh as a 3-tensor of shape (V, N, D) — V vertices per element, N elements, D spatial dimensions. Reshape-compatible with the legacy flat layout via reshape(x, V*N, D) (zero-copy). Spectral discretizations use N = 1 (a single notional "element" comprising every Chebyshev node).
• w::W: quadrature weights matching the flattened node order (length V*N).
• operators::Dict{Symbol,M_op}: fine-level discrete operators (e.g. :id, :dx, :dy).

MGBConvergenceFailure <: Exception

Thrown when the MGB solver fails to converge (feasibility or main phase). Carries a descriptive message and a machine-readable code:

• :infeasible: the feasibility subproblem converged to a minimizer with positive constraint violation strictly inside the bounding box, certifying that the problem is infeasible (see mgb_solve, section "Feasibility phase").
• :feasibility_Rmax: no strictly feasible point was found with nodal values bounded by feasibility_Rmax; the problem is infeasible, or its feasible points have values exceeding the cap.
• :stall: the t-ramp of the barrier method stopped making progress before reaching the target tolerance (the step refinement collapsed).
• :iteration_limit: the t-ramp hit the maxit iteration cap.
• :failure: any other convergence failure (the default of the one-argument constructor).

Front ends can dispatch on code; the JuMP extension maps it to the corresponding MOI.TerminationStatusCode.

MGBProblem{T,MT,FT,GT,QT,GeoT}

A fully assembled, array-valued, closure-free convex problem ready for the MGB solver. Produced by assemble; consumed by mgb_solve. Every closure that specifies the problem (f, g, and the constraint data A, b, p, select) has been lowered to per-vertex grids during assembly, so an MGBProblem is pure data: moving it between backends (native_to_device) is plain array movement, and a MGBProblem{T,...,CuArray,...} differs from its CPU sibling only by array type.

Fields

• M: the (main, feasibility) AMG hierarchy pair (_prepare_amg output).
• f: linear-term grid (f_grid).
• g: Dirichlet/initial-data grid (g_grid).
• Q::Convex: the convex constraint (its args are per-vertex grids).
• geometry: the fine-level Geometry (for the returned solution / plotting).

MultiGrid{T,M_R,G<:Geometry{T}}

A Geometry plus a multigrid hierarchy. For each state-variable subspace symbol (:dirichlet, :full, :uniform, a user-named Dirichlet class, or a discretization-specific extra such as FEM2D_P2's :broken_P1) the hierarchy stores a single family of fine-level prolongations R[X][l]: the matrix that lifts level-l-X-subspace coefficients directly to the fine broken basis (the R_{ℓ,X} of the paper). R[X][end] is the fine-level subspace embedding itself (identity for :full, the constant column for :uniform, the continuous zero-trace embedding for Dirichlet classes, the per-element-linear embedding for :broken_P1).

The per-level (level → level+1) transfer/restriction operators are an implementation detail of amg/geometric_mg: they are composed into the R[X][l] at construction time and not retained, since the solver only ever evaluates the barrier at the fine level and restricts to a coarse search space via R[X][l].

Fields

• geometry::G: the fine-level Geometry.
• R::Dict{Symbol,Vector{M_R}}: R[X][l] is the level-l → fine prolongation in subspace X.

Constructors

• MultiGrid(geometry, subspaces, refine) — stretches the per-subspace hierarchies to a common depth and composes R[X][l] = (refine[X] chain l→L) · subspaces[X][l]. Accepts either per-subspace Dict transfers or plain Vectors (replicated across every subspace key).
• MultiGrid(geometry, R) — store precomposed prolongations directly (e.g. the Kronecker construction in spectral2d).

SPECTRAL1D{T}

1D spectral discretization descriptor (Chebyshev). Field: n::Int (nodes).

SPECTRAL2D{T}

2D spectral discretization descriptor (tensor Chebyshev). Field: n::Int (nodes per dim).

TensorFEM{d, e, T}

Discretization descriptor for a d-dimensional tensor-product Q_k Lagrange FEM embedded in ambient dimension e ≥ d (an embedded manifold when e > d). Here d is the intrinsic dimension (number of reference axes) and e the ambient / embedding dimension (= size(geom.x, 3)); both are integers, T the scalar type.

Fields

• k::Int: polynomial order (each element carries (k+1)^d Lagrange-Chebyshev nodes).
• K::Array{T,3}: the Q1 corner mesh tensor of shape (2^d, N, e) — K[v, j, c] is ambient coordinate c of corner v of element j. Corners are in tensor-product order over {-1,+1}^d (axis-1 fastest). Informational — it records the corner input; the hierarchy builders work from the stored full node tensor geom.x.

FEM1D = TensorFEM{1}, FEM2D = TensorFEM{2}, FEM3D = TensorFEM{3} are aliases fixing only the intrinsic dimension (leaving e, T free): e.g. TensorFEM{2,2} is a planar quad mesh, TensorFEM{2,3} a surface in ℝ³, TensorFEM{1,3} a curve in ℝ³.

intersect(mg::MultiGrid, U::Convex{T}, rest...) where {T}

Return the intersection of convex domains as a single Convex{T}. Equivalent to convex_piecewise with all pieces active at all vertices.

amg(geom::Geometry; prolongator=amg_ruge_stuben(max_coarse=2), dirichlet_nodes=Dict(:dirichlet => find_boundary(geom))) -> MultiGrid   # FEM
amg(geom::Geometry) -> MultiGrid                                                                          # spectral

Build an algebraic-multigrid hierarchy on top of geom, returning a MultiGrid. Dispatched per discretization; the hierarchy's fine level matches geom.

Keyword arguments (FEM discretizations: FEM1D, FEM2D_P1, FEM2D_P2, FEM3D)

• prolongator = amg_ruge_stuben(max_coarse=2): the algebraic-multigrid hierarchy builder used to coarsen the auxiliary P1 problem. A prolongator is a callable mapping a SparseMatrixCSC{Float64,Int} stiffness to the vector of level prolongation matrices (finest → coarsest); these set the depth of the hierarchy. Three factory functions construct one while capturing the underlying AMG parameters:

  • amg_ruge_stuben(; kwargs...) — classical Ruge–Stüben (the default), via AlgebraicMultigrid.ruge_stuben. The depth parameter max_coarse is set through the factory, e.g. amg(geom; prolongator=amg_ruge_stuben(max_coarse=4)).
  • amg_smoothed_aggregation(; kwargs...) — smoothed aggregation, via AlgebraicMultigrid.smoothed_aggregation.
  • amg_pyamg(; solver=:rootnode, kwargs...) — the Python pyamg package (:rootnode energy-minimization, :smoothed_aggregation, or :ruge_stuben); provided by the MultiGridBarrierPyAMGExt extension (load PyCall).

  All kwargs are forwarded to the underlying solver.

• dirichlet_nodes::Dict{Symbol,Vector{Tuple{Int,Int}}} = Dict(:dirichlet => find_boundary(geom)): one entry per zero-trace ("dirichlet-style") subspace. Each key is a subspace symbol you may assign to a state-variable component via state_variables, and its value is the set of mesh nodes constrained to zero trace in that subspace, given as (vertex, element) index pairs (the format find_boundary returns).

  The default builds a single :dirichlet subspace constraining the whole boundary ∂Ω. To impose different Dirichlet boundaries on different components, name them and supply a node set for each, e.g. for z = (u, s) with u clamped on the north edge and s on the south:

  state_variables = [:u :dirichlet_north
                     :s :dirichlet_south]
  mg = amg(geom; dirichlet_nodes = Dict(:dirichlet_north => north_pairs,
                                        :dirichlet_south => south_pairs))

  Per entry, pass a subset of ∂Ω for mixed Dirichlet/Neumann conditions, Tuple{Int,Int}[] for pure Neumann, or a single pinned node to break the constant nullspace. The reserved subspaces :full (all-corners Neumann), :uniform (global constants) and, on FEM2D_P2 geometries, :broken_P1 (per-element linears, the pure-P2 slack space) are always available and must not appear as keys. These select which nodes are constrained; the boundary values (the Dirichlet lift g) are supplied separately to mgb_solve.

• auxiliary_postprocess::Function = identity (opt-in; the tensor-product family FEM1D/FEM2D/FEM3D, and FEM2D_P1 — not FEM2D_P2): a unary function applied to the all-corners (Neumann) auxiliary stiffness before it is fed to prolongator. Use to swap the geometric Galerkin matrix for a graph-Laplacian-style operator that AMG coarsens on graph topology alone. Example: pass the combinatorial Laplacian of the same sparsity (off-diag = -1, diag = degree) when running aggregation-based prolongators (amg_smoothed_aggregation, amg_pyamg(:rootnode)) on highly anisotropic problems near p=1 — that regime can otherwise blow up Newton iteration counts on the central path. The default identity is robust across the bench (the geometric stiffness encodes useful per-row scaling that RS uses), so this is an opt-in escape hatch rather than a recommended default.

Spectral discretizations (SPECTRAL1D, SPECTRAL2D)

amg(geom) takes no keyword arguments. The zero-trace subspace is built by basis truncation rather than node masking, so dirichlet_nodes does not apply.

See also find_boundary, geometric_mg, subdivide.

amg_ruge_stuben(; kwargs...) -> prolongator

Build a prolongator that runs classical Ruge–Stüben AMG via AlgebraicMultigrid.ruge_stuben (the package default). Any kwargs (e.g. max_coarse=4) are forwarded. Pass the result to amg(geom; prolongator=...).

amg_smoothed_aggregation(; kwargs...) -> prolongator

Build a prolongator that runs smoothed-aggregation AMG via AlgebraicMultigrid.smoothed_aggregation. Any kwargs (e.g. max_coarse=4) are forwarded. Pass the result to amg(geom; prolongator=...).

assemble(mg::MultiGrid{T}; kwargs...) -> MGBProblem

Lower a problem specification to a closure-free, CPU MGBProblem. This is the single place where the problem closures are evaluated to grids: f/g via map_rows, and the constraint closures inside Q at its construction. The result is a backend-agnostic, native data structure; to solve on a GPU, hand it to mgb_solve(prob; device=CUDADevice), which moves it to the device and the solution back.

The five MGBProblem fields are built from mg and the keyword arguments as:

Keyword Arguments

• dim::Integer = amg_dim(mg.geometry.discretization): spatial dimension; auto-detected.
• state_variables, D = default_D(dim): solution components / function spaces and the differential operators applied to them; together they define M. The default is [:u :dirichlet; :s :full], except on pure-P2 geometries (fem2d_P2(bubble=false)), where the slack defaults to the :broken_P1 subspace (zero corner quadrature weights make a :full slack ill-posed there; see FEM2D_P2).
• x = _xflat(mg): sample points where f/g are evaluated, one row per node.
• p::T = T(1.0): p-Laplace exponent for the default Q.
• f/f_grid, g/g_grid: forcing and boundary/initial data, as closures (lowered via map_rows) or as pre-built grids.
• Q::Convex{T}: convex constraint; defaults to a p-Laplace power-cone barrier.
• M: supply an AMG hierarchy pair directly instead of building it from mg.

Any extra (solver-control) keywords are accepted and ignored.

clear_cudss_cache!()

Destroy all cached cuDSS factorizations and free associated GPU memory.

convex_Euclidian_power(::Type{T}=Float64; mg, idx=Colon(), A=(x)->I, b=(x)->T(0), p=x->T(2), ...)

Create a convex set defined by Euclidean norm power constraints, with GPU support.

Constructs a Convex{T} representing the power cone: {y : s ≥ ‖q‖₂^p} where [q; s] = A(x)*y[idx] + b(x)

This is the fundamental constraint for p-Laplace problems where we need s ≥ ‖∇u‖^p for some scalar field u.

Arguments

• T::Type=Float64: Numeric type for computations

Keyword Arguments

• mg::MultiGrid: Required. The multigrid hierarchy (provides the fine grid).
• idx=Colon(): Indices of y to which transformation applies
• A::Function: Matrix function x -> A(x) for linear transformation
• b::Function: Vector function x -> b(x) for affine shift
• p::Function: Exponent function x -> p(x) where p(x) ≥ 1
• A_grid, b_grid, p_grid: Pre-computed fine grids; they default to sampling A, b, p at the mesh nodes, so pass them to skip the closures entirely.

Returns

A Convex{T} whose args carry the pre-computed fine parameter grids (A, b, p, μ); the barrier functors receive (args_rows..., y) per vertex via broadcasting (see Convex).

Mathematical Details

The barrier function is -log(s^(2/p) - ‖q‖²) - μ(p)*log(s), where μ(p) = 0 if p∈{1,2}, 1 if p<2, 2 if p>2. In particular:

• For p = 1: -log(s² - ‖q‖²)
• For p = 2: -log(s - ‖q‖²)

Examples

# Standard p-Laplace constraint with GPU support. With the default 1D layout
# Dz = (u, ∂u/∂x, s), idx = SVector(2, 3) selects (q, s) = (∂u/∂x, s).
mg = amg(fem1d(Float32; nodes=collect(range(-1f0, 1f0, length=33))))
Q = convex_Euclidian_power(Float32; mg=mg, idx=SVector(2, 3), p=x->1.5f0)

# Q is a single Convex{Float32}

convex_linear(::Type{T}=Float64; mg, idx=Colon(), A=(x)->I, b=(x)->T(0), A_grid=<sampled from A>, b_grid=<sampled from b>)

Create a convex set defined by linear inequality constraints, with GPU support.

Constructs a Convex{T} representing linear constraints. Defines F(y) = A * y[idx] + b where A, b are pre-computed per vertex. The interior is F > 0 (logarithmic barrier applied to each component).

Arguments

• T::Type=Float64: Numeric type for computations

Keyword Arguments

• mg::MultiGrid: Required. The multigrid hierarchy (provides the fine grid).
• idx=Colon(): Indices of y to which constraints apply (default: all)
• A::Function: Matrix function x -> A(x) for constraint coefficients
• b::Function: Vector function x -> b(x) for constraint bounds. A scalar-valued b is expanded to the constraint count.
• A_grid, b_grid: Pre-computed fine grids; they default to sampling A and b at the mesh nodes, so pass them to skip the closures entirely.

Returns

A Convex{T} whose barriers capture the pre-computed fine grids.

Examples

mg = amg(fem1d(Float32; nodes=collect(range(-1f0, 1f0, length=33))))

# Box constraints: -1 ≤ y ≤ 1
A_box(x) = SMatrix{4,2,Float32}(1,0,-1,0, 0,1,0,-1)
b_box(x) = SVector{4,Float32}(1, 1, 1, 1)
Q = convex_linear(Float32; mg=mg, A=A_box, b=b_box, idx=SVector(1, 2))

convex_piecewise(::Type{T}=Float64; Q::Tuple{Vararg{Convex{T}}}, mg, select::Function=x->(true,...), select_grid=<sampled from select>) where {T}

Build a single Convex{T} that combines multiple convex domains with spatial selectivity.

Arguments

• Q::Tuple{Vararg{Convex{T}}}: tuple of convex pieces, one Convex{T} each.
• mg::MultiGrid: multigrid hierarchy (provides the fine grid).
• select::Function: a function x -> Tuple{Bool,...} indicating which pieces are active at spatial position x.
• select_grid: pre-computed fine selection grid; defaults to sampling select at the mesh nodes, so pass it to skip the closure entirely.

Semantics

At each vertex, over the pieces k that are active there (select(x)[k] true):

• barrier = ∑ₖ Q[k].barrier
• cobarrier = ∑ₖ Q[k].cobarrier
• slack = maxₖ Q[k].slack

The slack is the maximum over active pieces, ensuring a single slack value suffices for feasibility.

Examples

# Intersection of two convex domains
U = convex_Euclidian_power(Float64; mg=mg, idx=SVector(1, 3), p=x->2)
V = convex_linear(Float64; mg=mg, A=x->A_matrix, b=x->b_vector)
select_both(x) = (true, true)
Qint = convex_piecewise(Float64; Q=(U, V), mg=mg, select=select_both)

# Region-dependent constraints
Q_left = convex_Euclidian_power(Float64; mg=mg, p=x->1.5)
Q_right = convex_Euclidian_power(Float64; mg=mg, p=x->2.0)
select(x) = (x[1] < 0, x[1] >= 0)
Qreg = convex_piecewise(Float64; Q=(Q_left, Q_right), mg=mg, select=select)

See also: intersect, convex_linear, convex_Euclidian_power.

cuda_to_native(x)

Convert a CUDA Geometry or MGBSOL back to native CPU types.

fem1d(::Type{T}=Float64; nodes=T[-1,1], k=1, K=<from nodes>, ambient=Val(1), t=nothing) -> Geometry

Construct a single-level 1D Qk FEM Geometry. Each element carries k+1 Lagrange-Chebyshev nodes; k=1 reproduces the legacy P1 discretization exactly. nodes is the strictly increasing vector of element endpoints. It defaults to [-1, 1]; the resulting mesh has one straight Qk element on each [nodes[i], nodes[i+1]].

The map is isoparametric: pass K as the full (k+1, n_e, 1) Lagrange-node tensor to give each element a nontrivial 1D parametrization (interior nodes off the affine positions), or as the (2, n_e, 1) endpoint tensor for straight elements. K may be passed on its own — nodes is only needed when K is not given.

The curve may be embedded in a higher ambient dimension e ∈ {2,3}: pass ambient=Val(2) or ambient=Val(3) (default Val(1)) together with a mesh K of that many coordinate columns, and fem1d builds a 1-manifold (curve) in ℝ^e of type TensorFEM{1,e,T}. The operators :dx, :dy[, :dz] are then the e ambient components of the intrinsic (arc-length) gradient — each tangent to the curve — and the quadrature weights are the arc-length measure √det(JᵀJ). Closed curves (circles, loops) glue automatically by coordinate dedup; this reduces exactly to the e=1 case on the diagonal.

Pass t (a (k+1, n_e) Integer matrix; t[v,e] is the global id of local node v in element e) to supply full-node connectivity explicitly instead of recovering it by coordinate dedup — needed when geometrically-coincident nodes must stay topologically distinct (slit domains, branch cuts, glued manifolds). Build one from corner connectivity with tensor_dofmap.

Attach a hierarchy with amg(geom).

fem2d(::Type{T}=Float64; k=1, K=<unit square>, ambient=Val(2), t=nothing) -> Geometry

Construct a single-level 2D Q_k FEM Geometry on quadrilaterals (k=1 is bilinear Q1). The map is isoparametric: pass K as the full ((k+1)^2, N, 2) Lagrange-node tensor (tensor order, axis-1 fastest) for curved quads, or as the (4, N, 2) Q1-corner tensor — order (-1,-1), (+1,-1), (-1,+1), (+1,+1) — for straight quads.

The quads may be embedded in ℝ³: pass ambient=Val(3) (default Val(2)) together with a mesh K of 3 coordinate columns, and fem2d builds a 2-manifold (surface) of type TensorFEM{2,3,T}. The operators :dx, :dy, :dz are then the three ambient components of the intrinsic tangential gradient ∇_Γ — tangent to the surface, so n·∇_Γ = 0 holds by construction — and the weights are the surface-area measure √det(JᵀJ). Closed surfaces (spheres, tori) glue by coordinate dedup. This reduces exactly to the planar discretization when e=2.

Pass t (a ((k+1)^2, N) Integer matrix; t[v,e] is the global id of local node v in element e) to supply full-node connectivity explicitly instead of recovering it by coordinate dedup — needed for slit domains, branch cuts, and glued manifolds, where geometrically-coincident nodes must stay topologically distinct. Build one from corner connectivity with tensor_dofmap.

Attach a hierarchy with amg(geom).

fem2d_P1(::Type{T}=Float64; K=<default unit-square>, t=<from K>) -> Geometry

Construct a single-level 2D FEM P1 Geometry on the doubled-per-element fine triangulation K. Use amg(geom) to attach an algebraic-multigrid hierarchy. (The legacy geometric_mg(geom, L) builds geometric-subdivision transfers instead.)

Arguments

• K::Array{T,3} (3 × N × 2): per-triangle corner tensor; K[v, t, d] is coordinate d of vertex v of triangle t.
• t::AbstractMatrix{<:Integer} (3 × N): optional corner connectivity. By default it is recovered from K by coordinate deduplication. Pass it explicitly when coincident vertices must remain topologically distinct.

Boundary conditions are chosen later, at amg time, via its dirichlet_nodes keyword (the default constrains the whole boundary).

fem2d_P2(::Type{T}=Float64; bubble=true, K=<default unit square>, t=<from K>) -> Geometry

Construct a single-level 2D simplicial P2 FEM Geometry. With bubble = true (the default) the element is P2 enriched with a per-element cubic bubble — 7 nodes per triangle, laid out as corner1, edge(1,2), corner2, edge(2,3), corner3, edge(3,1), centroid — whose nodal quadrature weights are strictly positive. With bubble = false the element is plain Lagrange P2 (the same layout without the centroid) and the nodal quadrature is the edge-midpoint rule, whose corner weights are zero; the slack of a variational problem then lives in the :broken_P1 subspace (see FEM2D_P2), which assemble selects automatically. Use amg(geom) to attach an algebraic-multigrid hierarchy. (The legacy geometric_mg(geom, L) builds geometric-subdivision transfers instead.)

Arguments

• bubble::Bool: selects the variant. When K is supplied its row count (7 or 6) determines the variant directly; passing a contradicting explicit bubble is an error.
• K::Array{T,3} (V × N × 2, V = bubble ? 7 : 6): per-triangle node mesh.
• t::AbstractMatrix{<:Integer} (V × N): optional full-node connectivity. By default it is recovered from K by coordinate deduplication. Pass it explicitly when coincident nodes must remain topologically distinct.

The element geometry is isoparametric: the map from the reference triangle is built from all V node positions via the shape functions, so displacing the edge (and centroid) nodes off the straight midpoints (barycenter) genuinely curves the element (with a node-varying Jacobian). Triangles must be orientation-preserving and non-self-intersecting — construction errors if any element's det J ≤ 0 at a node.

fem3d(::Type{T}=Float64; k=3, K=<unit cube>, t=nothing) -> Geometry

Construct a single-level 3D Q_k FEM Geometry on hexahedra. The map is isoparametric: pass K as the full ((k+1)^3, N, 3) Lagrange-node tensor (tensor order, axis-1 fastest) so displacing edge/face/interior nodes curves the hex (node-varying Jacobian), or as the (8, N, 3) Q1-corner tensor for straight hexes. The default is a single straight unit cube.

Pass t (a ((k+1)^3, N) Integer matrix; t[v,e] is the global id of local node v in element e) to supply full-node connectivity explicitly instead of recovering it by coordinate dedup (slit domains / glued manifolds). tensor_dofmap builds one from corner connectivity, though for k≥3 in 3D it does not yet number shared face-interior grids (it throws) — there, supply t by hand or use the dedup default.

Attach a hierarchy with amg(geom).

find_boundary(geom::Geometry) -> Vector{Tuple{Int,Int}}

Return the (v, e) index pairs of the mesh nodes on ∂Ω, in the per-element layout of geom.x (the 3-tensor of shape (V, N, D)): v is the local vertex index within element e. Duplicates are present — a corner shared by k elements contributes its k pairs (one per element that owns it).

amg(geom; dirichlet_nodes=Dict(:sym => set, …)) consumes these (v, e) pairs: each value set is a Vector{Tuple{Int,Int}} of the nodes constrained to zero trace in subspace :sym. A geometric position is treated as Dirichlet for that subspace iff any pair at that position is in set, so you can pass either the full set returned by find_boundary (or a subset) or a sparser representative-only set.

Defined for each FEM discretization (the tensor-product family FEM1D/FEM2D/FEM3D, plus FEM2D_P1 and FEM2D_P2). For spectral discretizations the zero-trace subspace is built by basis truncation rather than by node masking; the spectral amg does not accept dirichlet_nodes and find_boundary returns the perimeter Chebyshev nodes (paired with the single notional element index 1) for reference only.

Empty or singleton node sets are allowed; the resulting AMG problem may still be well-posed if the variational form carries a mass term.

find_boundary(geom::Geometry{...,<:FEM2D_P2{T}}) -> Vector{Tuple{Int,Int}}

(v, t) index pairs into geom.x for every P2(+bubble) DOF (corner vertex or edge midpoint) on ∂Ω. Centroids (bubble variant) never appear (they are interior by construction). Duplicates are present (a corner shared by k triangles contributes its k pairs; a boundary-edge midpoint contributes its single pair).

find_boundary(geom::Geometry{...,FEM2D_P1{T}}) -> Vector{Tuple{Int,Int}}

(v, t) index pairs into geom.x for every vertex on ∂Ω. A corner shared by k triangles contributes its k pairs (one per triangle that owns it).

find_boundary(geom::Geometry{...,SPECTRAL1D{T}}) -> Vector{Tuple{Int,Int}}

The two (v, t) pairs [(1, 1), (n, 1)] for the endpoint Chebyshev nodes. Spectral geometries use a single notional element (size(geom.x, 2) == 1) since there is no real per-element structure; the flattened view has one row per unique node. Informational only: the spectral amg builds the zero-trace subspace by basis truncation, not node masking; it does not accept dirichlet_nodes.

find_boundary(geom::Geometry{...,SPECTRAL2D{T}}) -> Vector{Tuple{Int,Int}}

(v, 1) pairs (single notional element) for every Chebyshev node on the perimeter of the tensor-product spectral 2D mesh. Informational only: the spectral amg builds the zero-trace subspace by basis truncation, not node masking; it does not accept dirichlet_nodes.

find_boundary(geom::Geometry{...,TensorFEM{d,e,T}}) -> Vector{Tuple{Int,Int}}

(v, e) index pairs into geom.x for every Q_k DOF on ∂Ω. A (d-1)-face used by exactly one element is on the boundary; every DOF on such a face is returned (corner / edge / face-interior, including the k≥2 interior-of-face nodes).

geometric_mg(geom::Geometry, L::Int) -> MultiGrid

Build a geometric-subdivision multigrid hierarchy of L levels on top of geom. The returned MultiGrid's geometry is the finest mesh (after L-1 levels of subdivision).

The spectral discretizations have no geometric subdivision: their geometric_mg ignores L and returns the same hierarchy as amg (so subdivide is a no-op for them).

interpolate(M::Geometry, z::Vector, t)

Interpolate a solution vector at specified points.

Given a solution z on the mesh M, evaluates the solution at new points t using the appropriate interpolation method for the discretization.

Supported discretizations

• 1D FEM (FEM1D): exact per-element degree-k Lagrange (Q_k) interpolation (piecewise linear at k = 1); evaluation points outside the mesh are clamped to the boundary values
• 1D spectral (SPECTRAL1D): spectral polynomial interpolation
• 2D spectral (SPECTRAL2D): tensor-product spectral interpolation

Note: 2D/3D FEM interpolation is not currently provided.

Arguments

• M::Geometry: The geometry containing grid and basis information
• z::Vector: Solution vector on the finest grid (length must match number of DOFs)
• t: Evaluation points. Format depends on dimension:
  • 1D: scalar or Vector{T} of x-coordinates
  • 2D spectral: Matrix{T} where each row is [x, y]

Returns

Interpolated values at the specified points. Shape matches input t.

Examples

# 1D interpolation (FEM)
geom = subdivide(fem1d(; nodes=collect(range(-1.0, 1.0, length=3))), 3)
z = sin.(π .* vec(geom.x))
y = interpolate(geom, z, 0.5)
y_vec = interpolate(geom, z, [-0.5, 0.0, 0.5])

# 2D interpolation (spectral)
geom = spectral2d(n=4)
xf = reshape(geom.x, :, size(geom.x, 3))   # flat (n_nodes, 2) view
z = exp.(-xf[:,1].^2 .- xf[:,2].^2)
points = [0.0 0.0; 0.5 0.5; -0.5 0.5]
vals = interpolate(geom, z, points)

linesearch_backtracking(::Type{T}=Float64; beta=T(0.5)) where {T}

Create a backtracking line search function for Newton methods.

Arguments

• T : numeric type for computations (default: Float64).

Keyword arguments

• beta : backtracking parameter for step size reduction (default: 0.5).
• c1 : Armijo sufficient-decrease constant c₁ (default: 0.1).

Returns

A line search function ls(x, y, g, n, F0, F1; printlog) where:

• x : current point (vector of type T).
• y : current objective value F0(x).
• g : current gradient F1(x).
• n : search direction (typically Newton direction H\g).
• F0 : objective function.
• F1 : gradient function.
• printlog : logging function.

Returns (xnext, ynext, gnext) where xnext = x - s*n for some step size s.

Algorithm

Implements the Armijo backtracking line search with sufficient decrease condition: F(x - s*n) ≤ F(x) - c₁ * s * ⟨∇F(x), n⟩, where c₁ is the keyword c1 (default 0.1). The step size starts at s = 1 and is reduced by factor beta until the condition is satisfied or numerical limits are reached.

Notes

This is a robust and commonly used line search that guarantees sufficient decrease in the objective function, making it suitable for general nonlinear optimization.

linesearch_illinois(::Type{T}=Float64; beta=T(0.5)) where {T}

Create an Illinois-based line search function for Newton methods.

Arguments

• T : numeric type for computations (default: Float64).

Keyword arguments

• beta : backtracking parameter for step size reduction when Illinois fails (default: 0.5).

Returns

A line search function ls(x, y, g, n, F0, F1; printlog) where:

• x : current point (vector of type T).
• y : current objective value F0(x).
• g : current gradient F1(x).
• n : Newton direction (typically H\g where H is the Hessian).
• F0 : objective function.
• F1 : gradient function.
• printlog : logging function.

Returns (xnext, ynext, gnext) where xnext = x - s*n for some step size s.

Algorithm

The Illinois algorithm finds a root of φ(s) = ⟨∇F(x - s*n), n⟩, which corresponds to the exact line search condition. If the Illinois solver fails or encounters numerical issues, the step size is reduced by factor beta and the process repeats.

Notes

This line search strategy aims for the exact minimizer along the search direction, making it potentially more aggressive than backtracking but also more expensive per iteration.

mgb_solve(prob::MGBProblem; device=..., kwargs...) -> MGBSOL

MultiGrid Barrier (MGB) solver for nonlinear convex optimization problems on a multigrid hierarchy. Operates in a feasibility phase followed by a main optimization phase, with damped Newton inner solves and line search.

Solves an assembled MGBProblem — build one with assemble (where the problem-specification keywords p, f, g, Q, state_variables, D, … live), or take one from the Zoo. The problem is moved to device, solved there, and the solution moved back, so the returned MGBSOL is always in native CPU types regardless of device.

Keyword Arguments

Backend

• device::Type{<:Device}: compute backend, CPUDevice or CUDADevice. The default is CPUDevice, unless the CUDA extension is loaded (using CUDA, CUDSS_jll) and a functional GPU is present, in which case solves run on CUDADevice automatically. The problem is moved to the device, solved there, and the solution moved back; the returned MGBSOL is always native.

Output Control

• verbose::Bool = true: progress bar.
• logfile = devnull: optional log stream.

Solver Control (forwarded to mgb_driver)

• tol, t, t_feasibility, feasibility_Rmax, maxit, kappa, early_stop, max_newton, stopping_criterion, line_search, finalize, progress, printlog.
• barrier_nodes: the constraint-collocation set — the nodes over which the main-phase barrier is flat-averaged — as a Vector{Bool} nodal mask, a vector of node indices, or Colon() (all nodes, the historical average). Defaults to the mask of nodes with nonzero quadrature weight: all nodes on every discretization except pure P2 (fem2d_P2(bubble=false)), where the triangle corners (zero midpoint-rule weight) drop out and the constraints are collocated at the edge midpoints only. The feasibility phase always uses the all-nodes barrier (see mgb_driver), so unselected nodes are constrained during phase I but unconstrained in the main solve.

Feasibility phase

If the initial point is infeasible, a phase-I subproblem minimizes the integral of a slack variable inside a bounding box |nodal values| < R (the box keeps the phase-I domain bounded, so its barrier has a minimizer; without it, iterates can drift to infinity). R starts at max(10, 10·max|g|) and is grown tenfold whenever the phase-I minimizer presses against the box, up to feasibility_Rmax (default 1/√eps(T)). If the phase-I minimizer has positive violation but sits strictly inside the box, the problem is reported infeasible via MGBConvergenceFailure (for a convex problem, an interior phase-I minimizer is global, so a larger box cannot help).

The handoff to the main phase does not stop at the first strictly feasible iterate: feasibility is tested only between completed (centered) t-steps, and once it first holds at t = tfirst the ramp continues until `t ≥ 2·tfirst`. Along the phase-I central path the slack suboptimality is bounded by the duality gap ν/t, so this certifies a feasibility margin of at least half the maximum achievable; a barely-feasible handoff point would sit on the main barrier's wall, where the Hessian is numerically singular.

Returns

An MGBSOL whose z is the fine-level solution matrix and whose geometry is the fine-level Geometry (the MultiGrid itself is not stored).

Examples

sol = mgb_solve(assemble(amg(fem1d(; nodes = collect(range(-1.0, 1.0, length=33)))); p = 1.5))
sol = mgb_solve(assemble(amg(subdivide(fem2d_P2(), 3)); p = 1.5))
sol = mgb_solve(assemble(amg(spectral2d(n = 8)); p = 2.0); device = CPUDevice)
sol = mgb_solve(Zoo.p_harmonic(amg(fem2d_P2()); p = 1.5); tol = 1e-4)

native_to_cuda(prob::MGBProblem) -> MGBProblem
native_to_cuda(geom::Geometry) -> Geometry

Convert a native (CPU) assembled MGBProblem — or one of its pieces: a Geometry, an AMG hierarchy, a Convex, raw arrays — to CUDA GPU types, faithfully preserving matrix types: BlockDiag operators stay block (driving the structured batched-GEMM Hessian assembly), sparse matrices become CuSparseMatrixCSR. FEM geometries always carry BlockDiag operators, so the solve is structured; only the spectral discretizations use dense operators. Requires using CUDA, CUDSS_jll.

parabolic_solve(mg::MultiGrid; kwargs...) -> ParabolicSOL

Solve time-dependent p-Laplace problems using implicit Euler timestepping on the multigrid hierarchy mg.

Keyword Arguments

Discretization

• state_variables = [:u :dirichlet; :s1 :full; :s2 :full].
• D: differential operators (default depends on spatial dimension).
• dim::Int: spatial dimension (auto-detected from mg).

Time Integration

• t0=T(0), t1=T(1), h=T(0.2), ts=t0:h:t1.

Problem Parameters

• p::T=T(1): exponent for the p-Laplacian.
• f1: source term (t, x) -> T (default: (t,x)->T(0.5)).
• g: initial/boundary conditions (t, x) -> Vector{T}.
• Q: convex constraints.

Advanced (grid-level) overrides

• f1_grid: per-node samples of f1 at every ts[j] (one column per time).
• g_grid: a function j -> grid, the sampled g at time ts[j].
• f_grid: a function (z, j) -> grid, the linear-term grid for step j given the previous state z; together with f_default (the per-row layout of that linear functional) the default implements implicit Euler for the p-Laplace flow.

Output Control

• verbose::Bool=true: progress bar.

Additional Parameters

• rest...: forwarded to mgb_solve for each time step.

Examples

sol = parabolic_solve(amg(fem2d_P2()); p=2.0, h=0.1)
sol = parabolic_solve(amg(fem1d(; nodes=collect(range(-1.0, 1.0, length=33)))); p=1.5, h=0.05)

spectral1d(::Type{T}=Float64; n=16) -> Geometry

Construct a 1D spectral single-level Geometry with n Chebyshev nodes. Use amg(geom) to attach a multigrid hierarchy.

spectral2d(::Type{T}=Float64; n=4) -> Geometry

Construct a 2D tensor-Chebyshev spectral single-level Geometry. Use amg(geom) to attach the spectral multigrid hierarchy.

stopping_exact(theta::T) where {T}

Create an exact stopping criterion for Newton methods.

Arguments

• theta : tolerance parameter for gradient norm relative decrease (type T).

Returns

A stopping criterion function with signature: stop(ymin, ynext, gmin, gnext, n, ndecmin, ndec) -> Bool

where:

• ymin : minimum objective value seen so far.
• ynext : current objective value.
• gmin : minimum gradient norm seen so far.
• gnext : current gradient vector.
• n : current Newton direction.
• ndecmin : square root of minimum Newton decrement seen so far.
• ndec : square root of current Newton decrement.

Algorithm

Returns true (stop) if both conditions hold:

1. No objective improvement: ynext ≥ ymin
2. Gradient norm stagnation: ‖gnext‖ ≥ theta * gmin

Notes

This criterion is "exact" in the sense that it requires both objective and gradient stagnation before stopping, making it suitable for high-precision optimization. Typical values of theta are in the range [0.1, 0.9].

stopping_inexact(lambda_tol::T, theta::T) where {T}

Create an inexact stopping criterion for Newton methods that combines Newton decrement and exact stopping conditions.

Arguments

• lambda_tol : tolerance for the Newton decrement (type T).
• theta : tolerance parameter for the exact stopping criterion (type T).

Returns

A stopping criterion function with signature: stop(ymin, ynext, gmin, gnext, n, ndecmin, ndec) -> Bool

where:

• ymin : minimum objective value seen so far.
• ynext : current objective value.
• gmin : minimum gradient norm seen so far.
• gnext : current gradient vector.
• n : current Newton direction.
• ndecmin : square root of minimum Newton decrement seen so far.
• ndec : square root of current Newton decrement (√(gᵀH⁻¹g)).

Algorithm

Returns true (stop) if either condition holds:

1. Newton decrement condition: ndec < lambda_tol
2. Exact stopping condition: stopping_exact(theta) is satisfied

Notes

This criterion is "inexact" because it allows early termination based on the Newton decrement, which provides a quadratic convergence estimate. The Newton decrement λ = √(gᵀH⁻¹g) approximates the distance to the optimum in the Newton metric. Typical values: lambda_tol ∈ [1e-6, 1e-3], theta ∈ [0.1, 0.9].

subdivide(geom::Geometry, L::Int) -> Geometry

Refine geom's mesh by L-1 levels of geometric subdivision and return a new fine-mesh Geometry, discarding the transfer operators that the geometric MG construction would otherwise produce. For FEM discretizations the fine operators are BlockDiag, so a subsequent amg(subdivide(geom, L)) runs the batched-GEMM (structured) Hessian assembly via the D_fine path.