# \[ANN\] MultiGridBarrier.jl: quasi-optimal solvers for nonsmooth convex variational problems
**URL:** https://discourse.julialang.org/t/ann-multigridbarrier-jl-quasi-optimal-solvers-for-nonsmooth-convex-variational-problems/138305
**Category:** Package Announcements
**Created:** [July 17, 2026, 8:20am UTC](https://discourse.julialang.org/t/ann-multigridbarrier-jl-quasi-optimal-solvers-for-nonsmooth-convex-variational-problems/138305 "2026-07-17T08:20:04Z")
**Posts on this page:** 2
**Page:** 1
### Author:  [@loisel](https://discourse.julialang.org/u/loisel)
#### Post date: [July 17, 2026, 8:20am UTC](https://discourse.julialang.org/t/ann-multigridbarrier-jl-quasi-optimal-solvers-for-nonsmooth-convex-variational-problems/138305/1 "2026-07-17T08:20:04Z")
[MultiGridBarrier.jl](https://github.com/sloisel/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](https://jump.dev) modeling layer, activated by loading JuMP (a package extension). The p-Laplacian, in conic form with an epigraph slack:
```julia
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](https://gmsh.info) (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:
```julia
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; a `device` keyword overrides).
- **Plotting** : loading PyPlot extends `plot` to 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](https://sloisel.github.io/MultiGridBarrier.jl/stable/)
- Software paper (PDF): [https://sloisel.github.io/MultiGridBarrier.jl/paper.pdf](https://sloisel.github.io/MultiGridBarrier.jl/paper.pdf)
- Repository: [GitHub - sloisel/MultiGridBarrier.jl: MultiGrid Barrier method · GitHub](https://github.com/sloisel/MultiGridBarrier.jl)
- Archived at Zenodo: [doi:10.5281/zenodo.20610950](https://doi.org/10.5281/zenodo.20610950)
The underlying theory is developed in Numerische Mathematik 146:369–400 (2020) ([doi](https://doi.org/10.1007/s00211-020-01141-z)) and 158:281–302 (2026) ([doi](https://doi.org/10.1007/s00211-025-01508-0)), and in the DD29 proceedings ([PDF](https://www.ddm.org/DD29/proceedings/ID74_pages.pdf)); citations welcome if you use the package in research.
Requires Julia 1.10+. Feedback, issues, and PRs are very welcome.
---
### Author:  [@loisel](https://discourse.julialang.org/u/loisel)
#### Post date: [July 17, 2026, 8:20am UTC](https://discourse.julialang.org/t/ann-multigridbarrier-jl-quasi-optimal-solvers-for-nonsmooth-convex-variational-problems/138305/2 "2026-07-17T08:20:19Z")
**v1.1.0 is out.** Perhaps the biggest change is a repositioning of the package’s front doors: the JuMP modeling layer is promoted to the main, high-level front end, Gmsh to the main mesh generator, and the native API to the low-level layer underneath them (still fully supported — its guided tour is now the [API Guide](https://sloisel.github.io/MultiGridBarrier.jl/stable/api_guide/)).
Highlights:
- **Constraint duals** : `dual(cr)` for every constraint type — contact pressures for inequalities and cones, per-node boundary reactions for Dirichlet conditions.
- **Standard JuMP accessors** throughout (`set_silent`, `solution_summary`, `is_solved_and_feasible`, `result` keyword, `MOI.OPTIMAL` on success), plus spatial data directly in the constraint algebra (`u == g`, `u >= phi_vals`, `f * u`) and two-sided constraints `lo <= u <= hi`.
- **Pure Lagrange P2** elements with a `barrier_nodes` collocation keyword restoring optimal rates.
- **Feasibility phase redesigned** , with infeasibility certification (`MOI.INFEASIBLE`).
- Gmsh: topology-preserving `subdivide`, exact node-tag connectivity (slit domains survive import).
- Console silence: nothing on stdout/stderr but the opt-in progress bar; diagnostics go to the solve log.
Full release notes (covering everything since v1.0.0): [Release v1.1.0 · sloisel/MultiGridBarrier.jl · GitHub](https://github.com/sloisel/MultiGridBarrier.jl/releases/tag/v1.1.0)