MultiGridBarrier.jl 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.
(This thread replaces the earlier versioned announcement; new releases will be announced as replies below.)
Model in JuMP syntax
The main front end is a JuMP modeling layer, activated by loading JuMP (a package extension). The p-Laplacian, in conic form with an epigraph slack:
using Pkg; Pkg.add(["MultiGridBarrier", "JuMP", "PyPlot"])
using MultiGridBarrier, JuMP, PyPlot
g = x -> x[1]^2 + x[2]^2 # boundary data
geom = subdivide(fem2d_P2(), 3)
m = MGBModel(geom)
set_silent(m)
@variable(m, u) # a conforming FEM function
@variable(m, s, Broken(), start = 100.0) # epigraph slack, one dof per node
set_start(u, g) # initial iterate and Dirichlet lift
@constraint(m, u == g, 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, with the multigrid hierarchy derived automatically from the geometry and the Dirichlet constraints. Spatial data — functions, constants, or raw per-node vectors — appears directly in the constraint algebra (u == g above), constraints restrict to regions with On, and value(u) returns the solution as a nodal vector. Constraint duals are available after the solve: dual(cr) returns the contact pressure of an obstacle constraint or the boundary reaction of a Dirichlet condition. The usual JuMP workflow works unchanged (set_silent, solution_summary, is_solved_and_feasible, termination_status == OPTIMAL, …), and FEM and spectral geometries both work.
Meshes from Gmsh
Real-world geometry comes from Gmsh (also a package extension): gmsh_import() converts the current Gmsh model or a .msh/.geo file into a mesh, and Gmsh physical groups become named node sets that plug directly into boundary conditions and On regions — @constraint(m, u == 0.0, On(gm.regions["boundary"])). Triangles import as P1/P2; quadrilaterals and hexahedra at any order, straight or curved.
The low-level API
Under the JuMP layer sits a native pipeline — each step a plain function returning a plain object:
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)
This is where the package’s full breadth lives: finite elements in 1D/2D/3D (simplicial P1/P2 and tensor-product Q_k at any order, all isoparametric), Chebyshev spectral discretizations, embedded manifolds (curves in R²/R³, surfaces in R³), slit domains and branch cuts via explicit connectivity, time-dependent problems with parabolic_solve, and custom convex constraints. A Zoo module ships classic problems as one-liners, e.g. mgb_solve(Zoo.rof(amg(geom))) for total-variation denoising.
Other features
- GPU:
using CUDA, CUDSS_jll— solves run on the GPU automatically (no code changes; adevicekeyword overrides). - Plotting: loading PyPlot extends
plotto every solution and geometry type — matplotlib in 1D/2D, PyVista in 3D. Python is entirely optional: the core package has no Python dependency.
Links
- Documentation: Home · MultiGridBarrier.jl 1.1.0
- Software paper (PDF): https://sloisel.github.io/MultiGridBarrier.jl/paper.pdf
- Repository: GitHub - sloisel/MultiGridBarrier.jl: MultiGrid Barrier method · GitHub
- Archived at Zenodo: doi:10.5281/zenodo.20610950
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.