[ANN] MultiGridBarrier.jl 1.0: quasi-optimal solvers for nonsmooth convex variational problems

Package Announcements
1

A nonsmooth p-Laplace solution computed with MultiGridBarrier.jl

I am happy to announce MultiGridBarrier.jl has reached 1.0, registered in General.

MultiGridBarrier solves convex variational problems and the corresponding nonlinear PDEs and boundary-value problems, including the nonsmooth ones that defeat most solvers: the p-Laplacian for every p ∈ [1, ∞], total-variation (ROF) denoising, obstacle problems, minimal surfaces, and friends. The core algorithm is the multigrid barrier method: an interior-point (barrier) method driven by a multigrid hierarchy. When the usual regularity conditions are satisfied, the solvers are provably quasi-optimal, i.e. the cost is nearly linear in the number of unknowns. The animation above is a p = 1 solution; that case is genuinely nonsmooth, and Newton-type methods generally fail on it.

Quickstart

using Pkg; Pkg.add("MultiGridBarrier")
using MultiGridBarrier
geom = fem2d_P2()                               # a 2D P2 triangular mesh
sol  = mgb_solve(assemble(amg(geom); p = 1.0))  # solve a nonsmooth p = 1 problem
plot(sol)

Every problem follows the same four-step pattern: build a mesh, attach a multigrid hierarchy with amg, assemble, mgb_solve. The mesh constructors cover finite elements in 1D/2D/3D (simplicial P1/P2 and tensor-product Q_k at any order k, all isoparametric) plus Chebyshev spectral discretizations. A Zoo module ships six classic problems as one-liners, e.g. mgb_solve(Zoo.rof(amg(geom))) for total-variation denoising.

Model in JuMP syntax

Loading JuMP activates a modeling front end (a package extension), so you can state variational problems with the standard macros. The p-Laplacian, in conic form with an epigraph slack:

using MultiGridBarrier, JuMP
geom = subdivide(fem2d_P2(), 3)
m = MGBModel(geom)
@variable(m, u)                       # a conforming FEM function
@variable(m, s, Broken())             # epigraph slack, one dof per node
set_start(u, x -> x[1]^2 + x[2]^2)    # initial iterate and Dirichlet lift
set_start(s, 100.0)
@constraint(m, u == Coef(m, x -> x[1]^2 + x[2]^2), On(find_boundary(geom)))
@constraint(m, [deriv(u, :dx); deriv(u, :dy); s] in EpiPower(1.5))  # s ≥ ‖∇u‖^1.5
@objective(m, Min, integral(0.5 * u + s))
optimize!(m)
plot(mgb_solution(m))

No MOI model is built; optimize! lowers the model directly to the multigrid barrier pipeline. Constraints can be restricted to named regions (On), spatial data can be given as functions, constants, or raw per-node vectors, and value(u) returns the solution in the same nodal format. FEM and spectral geometries both work.

Other features

  • Mesh import from Gmsh (package extension): gmsh_import("part.msh") reads triangles, quads, and hexes at any order, and converts physical groups into named regions ready for boundary conditions.

  • GPU: using CUDA, CUDSS_jll, then mgb_solve(prob; device = CUDADevice) runs the solve on the GPU.

  • Topological meshes: slit domains, branch cuts, and glued manifolds via explicit connectivity.

  • Time-dependent problems with parabolic_solve.

Links

The underlying theory is developed in Numerische Mathematik 146:369–400 (2020) (doi) and 158:281–302 (2026) (doi), and in the DD29 proceedings (PDF); citations welcome if you use the package in research.

Requires Julia 1.10+. Feedback, issues, and PRs are very welcome.

3 Likes