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

2

Hey, thanks for your work on this, I was curious and wanted to explore it. I noticed that loading the package pulls in the whole PyCall/Conda/matplotlib stack, because the module does using PyPlot at top level. It caught me off guard to see Python involved in a solver package at all, and as far as I can tell it’s only there for the plot helpers.

Would you consider moving the visualization behind an extension like you already do for CUDA, Gmsh, and JuMP? That way the core loads with just its Julia deps and Python only comes into play if someone wants to use your plotting features?

3

@kahliburke Thanks for your suggestion. Yes, I was thinking about exactly that. Actually, Python is also optionally used in the AMG construction but maybe that too could be an extension module.

4

@kahliburke I’ve fixed it now!