JuMP front end: support spectral geometries; bump to 1.0.1

The spectral families (spectral1d/spectral2d) build their multigrid hierarchy
by basis truncation, so their amg does not accept dirichlet_nodes and their
Dirichlet subspace is the fixed zero-trace space. optimize! now detects this
(the amg method lacks the dirichlet_nodes kwarg) and maps whole-boundary
equality constraints onto the hierarchy's :dirichlet subspace, empty ones onto
:full, and rejects a partial-boundary Dirichlet set with an explanatory error
(basis truncation cannot honor per-node conditions). The prolongator attribute
is likewise rejected for these geometries. FEM paths are unchanged.

Tests cross-validate spectral2d and spectral1d JuMP models against the
classical assemble/mgb_solve to solver tolerance, and check the
partial-boundary and prolongator rejections. Docs and the MGBModel/jump.md
notes updated; the front end no longer says spectral geometries are unwired.

This is the first release to ship the JuMP and Gmsh modeling front ends, which
landed after the 1.0.0 tag.

Author: Sébastien Loisel

Project.toml

@@ -1,7 +1,7 @@
 name = "MultiGridBarrier"
 uuid = "9e2c1f1d-9131-4ad4-b32f-bd2a0b0ecd1e"
 authors = ["Sébastien Loisel"]
-version = "1.0.0"
+version = "1.0.1"
 
 [deps]
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"

docs/src/jump.md

@@ -208,8 +208,11 @@ notion of continuity is the geometry's connectivity `geom.t`, so slit domains
 built with explicit connectivity keep their slits. The model must be in
 conic form (epigraph slacks are yours to declare, as with any conic solver);
 pointwise equality requires `On`; variable bounds and products of variable
-expressions are rejected with explanatory errors. `dual` and spectral
-geometries are not wired up yet.
+expressions are rejected with explanatory errors. Spectral geometries work
+too, with one restriction inherited from their hierarchy: the spectral
+Dirichlet subspace is built by basis truncation, so a Dirichlet condition
+there must cover exactly the whole boundary (`find_boundary(geom)`). `dual`
+is not wired up yet.
 
 ## API reference
 

ext/MultiGridBarrierJuMPExt/MultiGridBarrierJuMPExt.jl

@@ -22,7 +22,7 @@ import MultiGridBarrier
 import MultiGridBarrier: Geometry, MultiGrid, Convex, amg, assemble, mgb_solve,
     convex_linear, convex_Euclidian_power, convex_piecewise, MGBConvergenceFailure,
     On, Broken, Continuous, Uniform, deriv, integral, set_start, mgb_solution,
-    solver_log
+    solver_log, find_boundary
 const MOI = JuMP.MOI
 
 # ---------------------------------------------------------------------------
@@ -804,10 +804,38 @@ function JuMP.optimize!(m::MGBModel{T}) where {T}
     n = m.nnodes
     V = size(geom.x, 1)
 
-    amg_kw = Dict{Symbol,Any}(:dirichlet_nodes => low.dirichlet_nodes)
-    haskey(m.attrs, "prolongator") && (amg_kw[:prolongator] = m.attrs["prolongator"])
-    isempty(low.dirichlet_nodes) && pop!(amg_kw, :dirichlet_nodes)
-    mg = amg(geom; amg_kw...)
+    sv = low.state_variables
+    if hasmethod(amg, Tuple{typeof(geom)}, (:dirichlet_nodes,))
+        amg_kw = Dict{Symbol,Any}(:dirichlet_nodes => low.dirichlet_nodes)
+        haskey(m.attrs, "prolongator") && (amg_kw[:prolongator] = m.attrs["prolongator"])
+        isempty(low.dirichlet_nodes) && pop!(amg_kw, :dirichlet_nodes)
+        mg = amg(geom; amg_kw...)
+    else
+        # Basis-truncation hierarchies (the spectral families): the Dirichlet
+        # subspace is the fixed zero-trace space, so per-variable node sets
+        # cannot be honored. Accept exactly whole-boundary conditions and map
+        # them onto the hierarchy's :dirichlet / :full subspaces.
+        haskey(m.attrs, "prolongator") &&
+            _argerror("the \"prolongator\" attribute is not supported by this geometry's amg")
+        bdset = Set(find_boundary(geom))
+        sv = copy(low.state_variables)
+        for i in 1:size(sv, 1)
+            sub = sv[i, 2]
+            startswith(String(sub), "dirichlet_") || continue
+            pairs = low.dirichlet_nodes[sub]
+            if isempty(pairs)
+                sv[i, 2] = :full
+            elseif Set(pairs) == bdset
+                sv[i, 2] = :dirichlet
+            else
+                _argerror("this geometry builds its Dirichlet subspace by basis " *
+                    "truncation, so a Dirichlet condition must cover exactly the " *
+                    "whole boundary (find_boundary(geom)); variable $(sv[i, 1]) is " *
+                    "constrained on $(length(pairs)) of $(length(bdset)) boundary nodes")
+            end
+        end
+        mg = amg(geom)
+    end
 
     pieces = Convex{T}[_piece(m, mg, low, c) for c in low.cones]
     if length(pieces) == 1 && low.cones[1].pairs === nothing
@@ -826,7 +854,7 @@ function JuMP.optimize!(m::MGBModel{T}) where {T}
     end
 
     prob = assemble(mg;
-        state_variables = low.state_variables,
+        state_variables = sv,
         D = low.D,
         f_grid = low.f_grid,
         g_grid = low.g_grid,

src/jump_frontend.jl

@@ -15,8 +15,10 @@
     MGBModel(geom::Geometry)
 
 Construct a JuMP model over a fixed MultiGridBarrier discretization. Requires
-`using JuMP` (which loads the modeling extension). Build it from any FEM
-`Geometry` (e.g. `fem2d_P2()`, `subdivide(fem2d_P1(), 4)`), then use the standard
+`using JuMP` (which loads the modeling extension). Build it from any FEM or
+spectral `Geometry` (e.g. `fem2d_P2()`, `subdivide(fem2d_P1(), 4)`,
+`spectral2d(n = 16)`; spectral Dirichlet conditions must cover the whole
+boundary), then use the standard
 JuMP macros (`@variable`, `@constraint`, `@objective`); `optimize!` lowers the
 model to `amg` → `assemble` → `mgb_solve`, constructing the AMG hierarchy
 automatically from the geometry and the Dirichlet constraints. All spatial data

test/test_jump.jl

@@ -386,6 +386,53 @@ _maxdiff(ref, cols...) =
         @test cv[1] < maximum(value(uu)) + 1e-2        # ... and binds
     end
 
+    @testset "spectral geometries (whole-boundary Dirichlet)" begin
+        # Spectral hierarchies build their Dirichlet subspace by basis
+        # truncation, so the JuMP path maps whole-boundary equality constraints
+        # onto :dirichlet and rejects partial-boundary sets. The 2D model
+        # mirrors the classical assemble defaults, so it must match to solver
+        # tolerance.
+        gs = spectral2d(n = 8)
+        bs = find_boundary(gs)
+        g2 = x -> x[1]^2 + x[2]^2
+        ms = MGBModel(gs); _quiet!(ms)
+        @variable(ms, us); @variable(ms, ss, Broken())
+        set_start(us, g2); set_start(ss, 100.0)
+        @constraint(ms, us == Coef(ms, g2), On(bs))
+        @constraint(ms, [deriv(us, :dx); deriv(us, :dy); ss] in EpiPower(1.5))
+        @objective(ms, Min, integral(Coef(ms, 0.5) * us + ss))
+        optimize!(ms)
+        ref = mgb_solve(assemble(amg(gs); p = 1.5); verbose = false)
+        @test maximum(abs.(value(us) .- ref.z[:, 1])) < _JUMP_MATCH_TOL
+
+        # 1D, with the same data given explicitly on both sides
+        g1d = spectral1d(n = 8)
+        m1 = MGBModel(g1d); _quiet!(m1)
+        @variable(m1, w); @variable(m1, r, Broken())
+        set_start(w, x -> x[1]^2); set_start(r, 100.0)
+        @constraint(m1, w == Coef(m1, x -> x[1]^2), On(find_boundary(g1d)))
+        @constraint(m1, [deriv(w, :dx); r] in EpiPower(1.5))
+        @objective(m1, Min, integral(Coef(m1, 0.5) * w + r))
+        optimize!(m1)
+        ref1 = mgb_solve(assemble(amg(g1d); p = 1.5,
+            f = x -> (0.5, 0.0, 1.0), g = x -> (x[1]^2, 100.0)); verbose = false)
+        @test maximum(abs.(value(w) .- ref1.z[:, 1])) < _JUMP_MATCH_TOL
+
+        # partial-boundary Dirichlet is rejected with the truncation
+        # explanation, and so is the prolongator attribute
+        mpb = MGBModel(gs); _quiet!(mpb)
+        @variable(mpb, up); @variable(mpb, sp, Broken())
+        set_start(sp, 100.0)
+        @constraint(mpb, up == Coef(mpb, 0.0), On(bs[1:5]))
+        @constraint(mpb, [deriv(up, :dx); deriv(up, :dy); sp] in EpiPower(2.0))
+        @objective(mpb, Min, integral(1.0 * sp))
+        @test_throws ArgumentError optimize!(mpb)
+        err = try optimize!(mpb); "" catch e; sprint(showerror, e) end
+        @test occursin("whole boundary", err)
+        set_attribute(mpb, "prolongator", :unused)
+        @test_throws ArgumentError optimize!(mpb)
+    end
+
     @testset "convergence failure surfaces as OTHER_ERROR" begin
         # maxit too small for the barrier to reach its target t: optimize! must
         # catch MGBConvergenceFailure, report it, and refuse to hand out values.