module MultiGridBarrier

MultiGridBarrier solves nonlinear convex optimization problems in function spaces using a barrier (interior-point) method accelerated by a multigrid hierarchy constructed from your chosen discretization (FEM or spectral). The package provides simple, high-level entry points as well as a general solver that accept a "geometry" and optional keywords.

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:

\[\nabla \cdot \big(\|\nabla u\|_2^{p-2}\,\nabla u\big) = \tfrac{1}{p}\,f \quad \text{in } \Omega,\]

with the specified boundary conditions. This connection is obtained by integration by parts applied to the first variation of J(u).

Constrained linear reformulation with a slack variable

Introduce a slack $s(x) \geq \|\nabla u(x)\|_2^p$ and rewrite the objective using s:

\[\min_{u,\,s} \; \int_{\Omega} s + f\,u \, dx \quad \text{subject to}\quad s \ge \|\nabla u\|_2^p.\]

This is a convex optimization problem with a linear objective and convex constraints. In discrete form, we bundle the state into z, and apply a block "differential" operator D so that

\[D z = \begin{bmatrix} u \\ \nabla u \\ s \end{bmatrix}, \qquad c^\top = \begin{bmatrix} f & 0 & 1 \end{bmatrix}.\]

The problem becomes

\[\min_{z} \int_{\Omega} c(x)^\top \, (D z)(x) \, dx \quad \text{subject to}\quad (u,q,s) \in \mathcal{Q} := \{ s \ge \|q\|_2^p \}\ \text{pointwise},\]

which MultiGridBarrier solves by a barrier method. An illustrative (simple) barrier for $\mathcal{Q}$ is

\[\mathcal{F}(q,s) = -\log\!\big(s^{2/p} - \|q\|_2^2\big) - 2\log s,\]

and the method minimizes the barrier-augmented functional

\[\int_{\Omega} t\, c(x)^\top (D z)(x) + \mathcal{F}\!\big((D z)(x)\big) \, dx\]

for increasing barrier parameter t. Internally, the solve proceeds on a hierarchy of grids with damped Newton steps and line search, but these details are abstracted away.

How to use it (discretizations and solvers)

• Solve with a convenience wrapper (recommended to start):
  • sol = fem1d_solve(; kwargs...)
  • sol = fem2d_solve(; kwargs...)
  • sol = spectral1d_solve(; kwargs...)
  • sol = spectral2d_solve(; kwargs...)
• Or call the general solver directly:
  • sol = amgb(geometry; kwargs...) → AMGBSOL
• The solution can be plotted by calling plot(sol). If using amgb() directly, you must construct a suitable geometry object:
  • geometry = fem1d(; L=4) → 1D FEM on [-1, 1] with 2^L elements
  • geometry = fem2d(; L=2, K=...) → 2D FEM (quadratic + bubble triangles)
  • geometry = spectral1d(; n=16) → 1D spectral (Chebyshev/Clenshaw–Curtis)
  • geometry = spectral2d(; n=4) → 2D spectral (tensor Chebyshev)

Quick examples

# 1D FEM p-Laplace
z = fem1d_solve(L=5, p=1.0).z

# 2D spectral p-Laplace
z = spectral2d_solve(n=8, p=2.0).z

# 2D FEM with custom boundary data
g_custom(x) = [sin(π*x[1])*sin(π*x[2]), 10.0]
z = fem2d_solve(L=3; p=1.0, g=g_custom).z

# Time-dependent (implicit Euler)
sol = parabolic_solve(fem2d(L=3); h=0.1)
# plot(sol) animates the first component

Inputs and defaults (high level)

• p::Real = 1.0: exponent in the p-Laplace term
• g, f: boundary/initial data and forcing; either as functions g(x), f(x) or as grids g_grid, f_grid
• D and state_variables: symbolic specifications of which operators act on which variables (sane defaults provided based on the geometry’s dimension)
• Q: convex set (by default, a p-Laplace-compatible set via convex_Euclidian_power)
• verbose, logfile: visualization and logging
• Advanced control: tol, t, t_feasibility, line_search, stopping_criterion, finalize

What you get back

• Static solvers (amgb, *_solve) return an AMGBSOL with fields:
  • z::X (typically Matrix{T}): solution on the finest grid (nodes × components)
  • SOL_main, SOL_feasibility: per-phase diagnostics
  • log::String: textual log for debugging
  • geometry: the Geometry used to construct the multilevel operators
  The solution object supports plot(sol) to visualize the first component.
• The time-dependent solver parabolic_solve returns a ParabolicSOL with fields:
  • geometry, ts::Vector, u::Array(nodes × components × timesteps)
  Call plot(parabolic_sol) to animate using ts (see plot docs for timing options).

Utilities

• interpolate(geometry, z, points): evaluate the discrete solution at arbitrary points
• plot(sol) or plot(geometry, z): plot 1D curves or 2D surfaces
• plot(geometry, ts, U; frame_time=..., embed_limit=..., printer=...): animate a time sequence at absolute times ts (seconds), e.g., from parabolic_solve
• Convex set helpers: convex_Euclidian_power, convex_linear, convex_piecewise, intersect

Errors and diagnostics

• Throws AMGBConvergenceFailure if the feasibility subproblem or the main solve cannot converge
• Set verbose=true for a progress bar; inspect SOL_main/feasibility and log for details

See also

• Discretizations: fem1d, fem2d, spectral1d, spectral2d
• Solvers: amgb, fem1d_solve, fem2d_solve, spectral1d_solve, spectral2d_solve, parabolic_solve
• Convex: convex_Euclidian_power, convex_linear, convex_piecewise, intersect
• Visualization & sampling: plot, interpolate

AMGBConvergenceFailure <: Exception

Thrown when the AMGB solver fails to converge (feasibility or main phase). Includes a descriptive message about the failure.

Convex{T}

Container for a convex constraint set used by AMGB.

Fields:

• barrier(x, y): barrier of the set
• cobarrier(x, yhat): barrier with slack for feasibility
• slack(x, y): initial slack value

Construct via helpers like convex_linear, convex_Euclidian_power, convex_piecewise, or intersect.

FEM1D{T}

1D FEM geometry descriptor. Field: L::Int (levels). Use with amgb.

FEM2D{T}

2D FEM geometry descriptor for quadratic+bubble triangles. Fields: K::Matrix{T} (3n×2 mesh), L::Int (levels). Use with amgb.

Geometry{T,X,W,M,Discretization}

Container for discretization geometry and the multigrid transfer machinery used by AMGB.

Constructed by high-level front-ends like fem1d, fem2d, spectral1d, and spectral2d. It collects the physical/sample points, quadrature weights, per-level subspace embeddings, discrete operators (e.g. identity and derivatives), and intergrid transfer operators (refine/coarsen).

Type parameters

• T: scalar numeric type (e.g. Float64)
• X: type of the point storage for x (typically Matrix{T}; may be any AbstractArray with size (n_nodes, dim))
• W: type of the weight storage for w (typically Vector{T}; may be any AbstractVector)
• M: matrix type used for linear operators (e.g. SparseMatrixCSC{T,Int} or Matrix{T})
• Discretization: front-end descriptor (e.g. FEM1D{T}, FEM2D{T}, SPECTRAL1D{T}, SPECTRAL2D{T})

Fields

• discretization::Discretization: Discretization descriptor that encodes dimension and grid construction
• x::X: Sample/mesh points on the finest level; size is (n_nodes, dim)
• w::W: Quadrature weights matching x (length n_nodes)
• subspaces::Dict{Symbol,Vector{M}}: Per-level selection/embedding matrices for function spaces (keys commonly include :dirichlet, :full, :uniform). Each value is a vector of length L with one matrix per level.
• operators::Dict{Symbol,M}: Discrete operators defined on the finest level (e.g. :id, :dx, :dy). Operators at other levels are obtained via coarsen_fine * operator * refine_fine inside amg.
• refine::Vector{M}: Level-to-level refinement (prolongation) matrices for the primary state space
• coarsen::Vector{M}: Level-to-level coarsening (restriction) matrices for the primary state space

Notes

• Geometry is consumed by amg to build an AMG hierarchy and by utilities like interpolate and plot.
• The length of refine/coarsen equals the number of levels L; the last entry is typically the identity.

SPECTRAL1D{T}

1D spectral geometry descriptor (Chebyshev). Field: n::Int (nodes). Use with amgb.

SPECTRAL2D{T}

2D spectral geometry descriptor (tensor Chebyshev). Field: n::Int (nodes per dim). Use with amgb.

intersect(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.

amgb(geometry::Geometry{T,X,W,Mat,Discretization}=fem1d(); kwargs...) where {T, X, W, Mat, Discretization}

Algebraic MultiGrid Barrier (AMGB) solver for nonlinear convex optimization problems in function spaces using multigrid barrier methods.

This is the main high-level entry point for solving p-Laplace and related problems using the barrier method with multigrid acceleration. The solver operates in two phases:

1. Feasibility phase: Finds an interior point for the constraint set (if needed)
2. Main optimization phase: Solves the barrier-augmented optimization problem

Arguments

• geometry: Discretization geometry (default: fem1d()). Options:
  • fem1d(L=n): 1D finite elements with 2^L elements
  • fem2d(L=n, K=mesh): 2D finite elements
  • spectral1d(n=m): 1D spectral with m nodes
  • spectral2d(n=m): 2D spectral with m×m nodes

Keyword Arguments

Problem Specification

• dim::Integer = amg_dim(geometry.discretization): Problem dimension (1 or 2), auto-detected from geometry
• state_variables::Matrix{Symbol} = [:u :dirichlet; :s :full]: Solution components and their function spaces
• D::Matrix{Symbol} = default_D[dim]: Differential operators to apply to state variables
• x = geometry.x: Mesh/sample points where f and g are evaluated when they are functions

Problem Data

• p::T = 1.0: Exponent for p-Laplace operator (p ≥ 1)
• g::Function = default_g(T)[dim]: Boundary conditions/initial guess (function of spatial coordinates)
• g_grid::X (same container type as x): Alternative to g, directly provide values on grid (default: g evaluated at x)
• f::Function = default_f(T)[dim]: Forcing term/cost functional (function of spatial coordinates)
• f_grid::X (same container type as x): Alternative to f, directly provide values on grid (default: f evaluated at x)
• Q::Convex{T} = convex_Euclidian_power(T, idx=2:dim+2, p=x->p): Convex constraint set

Output Control

• verbose::Bool = true: Display progress bar during solving
• logfile = devnull: IO stream for logging (default: no file logging)

Solver Control (passthrough)

Additional keyword arguments are forwarded to the internal solvers:

• tol = sqrt(eps(T)): Barrier tolerance; the method stops once 1/t > 1/tol
• t = T(0.1): Initial barrier parameter for the main solve
• t_feasibility = t: Initial barrier parameter for the feasibility solve
• maxit = 10000: Maximum number of outer (barrier) iterations
• kappa = T(10.0): Multiplier controlling barrier growth (t ← min(kappa*t, ...))
• c0 = T(0): Base offset added to the objective (c0 + t*c)
• early_stop = z->false: Callback z -> Bool; if true, halt early (e.g. after feasibility is reached)
• max_newton = ceil((log2(-log2(eps(T))))+2): Max Newton iterations per inner solve
• stopping_criterion = stopping_inexact(sqrt(minimum(M[1].w))/2, T(0.5)): Newton stopping rule
  • use stopping_exact(theta) or stopping_inexact(lambda_tol, theta)
• line_search = linesearch_backtracking(T): Line search strategy (linesearch_backtracking or linesearch_illinois)
• finalize = stopping_exact(T(0.5)): Extra final Newton pass with stricter stopping
• progress = x->nothing: Progress callback in [0,1]; when verbose=true this is set internally
• printlog = (args...)->nothing: Logging callback; overridden by logfile when using amgb

Default Values

The defaults for f, g, and D depend on the problem dimension:

1D Problems

• f(x) = [0.5, 0.0, 1.0] - Forcing term
• g(x) = [x[1], 2] - Boundary conditions
• D = [:u :id; :u :dx; :s :id] - Identity, derivative, identity

2D Problems

• f(x) = [0.5, 0.0, 0.0, 1.0] - Forcing term
• g(x) = [x[1]²+x[2]², 100.0] - Boundary conditions
• D = [:u :id; :u :dx; :u :dy; :s :id] - Identity, x-derivative, y-derivative, identity

Returns

Solution object with fields:

• z: Solution matrix of size (n_nodes, n_components) containing the computed solution
• SOL_feasibility: Feasibility phase results (nothing if the initial point was already feasible), otherwise a solution object (see below)
• SOL_main: Main optimization phase results as a solution object (see below)
• log: String containing detailed iteration log for debugging
• geometry: The input geometry object

Each solution object (SOL_feasibility and SOL_main) is a NamedTuple containing:

• z: Solution vector (flattened; for feasibility phase includes auxiliary slack variable)
• z_unfinalized: Solution before final refinement step
• c: Cost functional used in this phase
• its: Iteration counts across levels and barrier steps (L×k matrix where L is number of levels, k is number of barrier iterations)
• ts: Sequence of barrier parameters t used (length k)
• kappas: Step size multipliers used at each iteration (length k)
• times: Wall-clock timestamps for each iteration (length k)
• t_begin, t_end, t_elapsed: Timing information for this phase
• passed: Boolean array indicating phase 1 success at each level
• c_dot_Dz: Values of ⟨c, D*z⟩ at each barrier iteration (length k)

Algorithm Overview

The AMGB method combines:

1. Interior point method: Uses logarithmic barriers to handle constraints
2. Multigrid acceleration: Solves on hierarchy of grids from coarse to fine
3. Damped Newton iteration: Inner solver with line search for robustness

The solver automatically handles:

• Construction of appropriate discretization and multigrid hierarchy
• Feasibility restoration when initial point is infeasible
• Adaptive barrier parameter updates with step size control
• Convergence monitoring across multiple grid levels
• Progress reporting (when verbose=true) and logging (to logfile if specified)

Errors

Throws AMGBConvergenceFailure if:

• The feasibility problem cannot be solved (problem may be infeasible)
• The main optimization fails to converge within maxit iterations
• Newton iteration fails at any grid level

Examples

# Solve 1D p-Laplace problem with p=1.5 using FEM
z = amgb(fem1d(L=4); p=1.5).z

# Solve 2D problem with spectral elements
z = amgb(spectral2d(n=8); p=2.0).z

# Custom boundary conditions
g_custom(x) = [sin(π*x[1])*sin(π*x[2]), 10.0]
z = amgb(fem2d(L=3); g=g_custom).z

# Get detailed solution information
sol = amgb(fem1d(L=3); verbose=true)
println("Iterations: ", sum(sol.SOL_main.its))
println("Final barrier parameter: ", sol.SOL_main.ts[end])

# Visualize the first component
plot(sol)

# Log iterations to a file
open("solver.log", "w") do io
    amgb(fem2d(L=2); logfile=io, verbose=false)
end

See Also

• fem1d_solve, fem2d_solve, spectral1d_solve, spectral2d_solve: Convenience wrappers for specific discretizations
• Convex: Constraint set specification type

apply_D(D,z) = hcat([D[k]*z for k in 1:length(D)]...)

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

Create a convex set defined by Euclidean norm power constraints.

Constructs a Convex{T} object 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

• 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

Returns

Convex{T} object with logarithmic barrier for the power cone

Mathematical Details

The barrier function is:

• For p = 2: -log(s² - ‖q‖²)
• For p ≠ 2: -log(s^(2/p) - ‖q‖²) - μ(p)*log(s) where μ(p) = 0 if p∈{1,2}, 1 if p<2, 2 if p>2

Examples

# Standard p-Laplace constraint: s ≥ ‖∇u‖^p
Q = convex_Euclidian_power(; idx=2:4, p=x->1.5)

# Spatially varying exponent
p_var(x) = 1.0 + 0.5 * x[1]  # variable p
Q = convex_Euclidian_power(; p=p_var)

# Second-order cone constraint: s ≥ ‖q‖₂
Q = convex_Euclidian_power(; p=x->1.0)

convex_linear(::Type{T}=Float64; idx=Colon(), A=(x)->I, b=(x)->T(0))

Create a convex set defined by linear inequality constraints.

Defines F(x, y) = A(x) * y[idx] + b(x). The interior of the set is given by F(x, y) > 0 (a logarithmic barrier is applied to each component of F). The boundary F(x, y) = 0 corresponds to constraint activation.

Arguments

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

Keyword Arguments

• 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

Returns

Convex{T} object with appropriate barrier functions

Examples

# Box constraints in 2D: -1 ≤ y ≤ 1
A_box(x) = [1.0 0.0; 0.0 1.0; -1.0 0.0; 0.0 -1.0]
b_box(x) = [1.0, 1.0, 1.0, 1.0]
Q = convex_linear(; A=A_box, b=b_box, idx=1:2)

# Single linear constraint: y[1] + 2*y[2] ≤ 3
# Choose F = 3 - (y1 + 2*y2) > 0
A_single(x) = [-1.0 -2.0]
b_single(x) = [3.0]
Q = convex_linear(; A=A_single, b=b_single, idx=1:2)

convex_piecewise(::Type{T}=Float64; Q::Vector{Convex{T}}, select::Function=(tr=fill(true,length(Q));x->tr)) where {T}

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

Arguments

• Q::Vector{Convex{T}}: a vector of convex pieces to be combined.
• select::Function: a function x -> Vector{Bool} indicating which pieces are active at x. Default: all pieces active everywhere (equivalent to intersection).

Semantics

For sel = select(x), the resulting convex domain has:

• barrier(x, y) = ∑(Q[k].barrier(x, y) for k where sel[k])
• cobarrier(x, yhat) = ∑(Q[k].cobarrier(x, yhat) for k where sel[k])
• slack(x, y) = max(Q[k].slack(x, y) for k where sel[k])

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

Use cases

1. Intersections (default): All pieces active everywhere creates Q₁ ∩ Q₂ ∩ ...
2. Spatial switching: Different constraints in different regions
3. Conditional constraints: Activate constraints based on solution state

Examples

# Intersection (using default select)
U = convex_Euclidian_power(Float64; idx=[1, 3], p = x->2)
V = convex_linear(Float64; A = x->A_matrix, b = x->b_vector)
Qint = convex_piecewise(Float64; Q = [U, V])  # U ∩ V everywhere

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

# Conditional activation
Q_base = convex_linear(Float64; A = x->I, b = x->-ones(2))
Q_extra = convex_Euclidian_power(Float64; p = x->3)
select(x) = [true, norm(x) > 0.5]  # extra constraint outside radius 0.5
Qcond = convex_piecewise(Float64; Q = [Q_base, Q_extra], select = select)

See also: intersect, convex_linear, convex_Euclidian_power.

fem1d(::Type{T}=Float64; L=4, kwargs...)

Construct 1D FEM geometry (piecewise linear) on [-1, 1]. Returns a Geometry suitable for use with amgb. Keyword L sets 2^L elements.

fem1d_solve(::Type{T}=Float64;rest...) where {T} = amgb(fem1d(T;rest...);rest...)

fem2d(::Type{T}=Float64; L=2, K=T[-1 -1;1 -1;-1 1;1 -1;1 1;-1 1], kwargs...)

Construct 2D FEM geometry (quadratic + bubble) on a triangular mesh. Returns a Geometry suitable for use with amgb. Keywords: L levels, K 3n×2 vertices.

fem2d_solve(::Type{T}=Float64;rest...) where {T} = amgb(fem2d(T;rest...);rest...)

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): piecewise-linear interpolation
• 1D spectral (SPECTRAL1D): spectral polynomial interpolation
• 2D spectral (SPECTRAL2D): tensor-product spectral interpolation

Note: 2D 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 = fem1d(L=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)
z = exp.(-geom.x[:,1].^2 .- geom.x[:,2].^2)
points = [0.0 0.0; 0.5 0.5; -0.5 0.5]
vals = interpolate(geom, z, points)

parabolic_solve(geometry::Geometry{T,X,W,Mat,Discretization}=fem2d(); kwargs...) where {T, X, W, Mat, Discretization}

Solve time-dependent p-Laplace problems using implicit Euler timestepping.

Solves the parabolic PDE:

\[u_t - \nabla \cdot (\|\nabla u\|_2^{p-2}\nabla u) = -f_1\]

using implicit Euler discretization and barrier methods.

Arguments

• geometry: Discretization geometry (default: fem2d()).

Keyword Arguments

Discretization

• state_variables: State variables (default: [:u :dirichlet; :s1 :full; :s2 :full]).
• D: Differential operators (default depends on spatial dimension).
• dim::Int: Spatial dimension (auto-detected from geometry).

Time Integration

• t0::T=0: Initial time.
• t1::T=1: Final time.
• h::T=0.2: Time step size.
• ts::AbstractVector{T}=t0:h:t1: Time grid; override to provide a custom, nonuniform, nondecreasing sequence.

Problem Parameters

• p::T=1: Exponent for the p-Laplacian.
• f1: Source term function of signature (t, x) -> T (default: (t,x)->T(0.5)).
• g: Initial/boundary conditions function of signature (t, x) -> Vector{T} (default depends on dimension).
• Q: Convex constraints (default: appropriate for p-Laplace).

Output Control

• verbose::Bool=true: Show a progress bar during time stepping.

Additional Parameters

• rest...: Passed through to amgb for each time step.

Returns

A ParabolicSOL with fields:

• geometry: the Geometry used.
• ts::Vector{T}: time stamps (seconds).
• u::Array{T,3}: solution tensor of size (n_nodes, n_components, n_timesteps).

Animate with plot(sol) (or plot(sol, k) for component k). To save to a file, use the plotting printer, e.g. plot(sol; printer=anim->anim.save("out.mp4")).

Mathematical Formulation

The implicit Euler scheme $u_t \approx (u_{k+1}-u_k)/h$ gives:

\[u_{k+1} - h\nabla \cdot (\|\nabla u_{k+1}\|^{p-2}\nabla u_{k+1}) = u_k - h f_1.\]

We minimize the functional:

\[J(u) = \int_\Omega \tfrac{1}{2}u^2 + \tfrac{h}{p}\|\nabla u\|^p + (h f_1 - u_k)u \, dx.\]

With slack variables $s_1 \ge u^2$ and $s_2 \ge \|\nabla u\|^p$, this becomes:

\[\min \int_\Omega \tfrac{1}{2}s_1 + \tfrac{h}{p}s_2 + (h f_1 - u_k)u \, dx.\]

Examples

# Basic 2D heat equation (p=2)
sol = parabolic_solve(; p=2.0, h=0.1)

# 1D p-Laplace with custom parameters
sol = parabolic_solve(fem1d(L=5); p=1.5, h=0.05, t1=2.0)

# Spectral discretization
sol = parabolic_solve(spectral2d(n=8); verbose=true)

# Custom initial condition
g_init(t, x) = [exp(-10*(x[1]^2 + x[2]^2)), 0, 0]
sol = parabolic_solve(; g=g_init)

See Also

• amgb: Single time step solver
• plot: Animation and plotting function

spectral1d(::Type{T}=Float64; n=16, kwargs...)

Construct 1D spectral geometry with n Chebyshev nodes (degree n-1). Returns a Geometry suitable for use with amgb.

spectral1d_solve(::Type{T}=Float64;rest...) where {T} = amgb(spectral1d(T;rest...);rest...)

spectral2d(::Type{T}=Float64; n=4, kwargs...)

Construct 2D spectral geometry with n×n Chebyshev nodes on [-1,1]^2. Returns a Geometry suitable for use with amgb.

spectral2d_solve(::Type{T}=Float64;rest...) where {T} = amgb(spectral2d(T;rest...);rest...)

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].

plot(sol::AMGBSOL, k::Int=1; kwargs...)
plot(sol::ParabolicSOL, k::Int=1; kwargs...)
plot(M::Geometry, z::Vector; kwargs...)
plot(M::Geometry, ts::AbstractVector, U::Matrix; frame_time=..., embed_limit=..., printer=...)

Visualize solutions and time sequences on meshes.

• 1D problems: Line plot. For spectral methods, you can specify evaluation points with x=-1:0.01:1.
• 2D FEM: Triangulated surface plot using the mesh structure.
• 2D spectral: 3D surface plot. You can specify evaluation grids with x=-1:0.01:1, y=-1:0.01:1.

Time sequences (animation):

• Call plot(M, ts, U; frame_time=1/30, printer=anim->nothing) where U has columns as frames and ts are absolute times in seconds (non-uniform allowed).
• Or simply call plot(sol) where sol is a ParabolicSOL returned by parabolic_solve (uses sol.ts).
• Animation advances at a fixed frame rate given by frame_time (seconds per video frame). For irregular ts, each video frame shows the latest data frame with timestamp ≤ current video time.
• The printer callback receives the Matplotlib animation object; use it to display or save (e.g., anim.save("out.mp4")).
• embed_limit controls the maximum embedded HTML5 video size in megabytes.

When sol is a solution object returned by amgb or the *_solve helpers, plot(sol,k) plots the kth component sol.z[:, k] using sol.geometry. plot(sol) uses the default k=1.

All other keyword arguments are passed to the underlying PyPlot functions.