benchmark: S4-B hierarchy cases and h-multigrid configurations

benchmark/solvers/cases.jl

@@ -96,8 +96,8 @@ end
 
 # Mechanical block alone (the clean S1 measurement, D4 refined 2026-10-30): the EM geometries with the
 # neo-Hookean solid only, clamped at "fixedup" and sheared by a prescribed displacement of "topsuf".
-function mech_case(model, order, utop)
-    physmodel = NeoHookean3D(λ=10.0, μ=1.0)
+function mech_case(model, order, utop; λ=10.0)
+    physmodel = NeoHookean3D(λ=λ, μ=1.0)                  # λ = 10: ν = λ/(2(λ+μ)) = 0.45; λ = 1.5: ν = 0.3
     k = Kinematics(Mechano, Solid)
     Ω = Triangulation(model); dΩ = Measure(Ω, 2order)
     Du = BoundaryConditions(BoundaryCondition("fixedup", [0.0, 0.0, 0.0]), BoundaryCondition("topsuf", utop, Λ -> Λ))
@@ -136,6 +136,88 @@ function mech_cube(n::Int, p::Int=2)
     mech_case(m, p, [0.0, 0.0, 0.1])
 end
 
+# --- S4-B: the mechanical cube as a refined hierarchy (design/solvers-s4b-plan.md §4) ---------------------------------
+
+"Coarse Gmsh tetrahedral mesh of the unit cube with mesh size `h`; closed physical groups fixedup (x = 0) and topsuf (x = 1)."
+function cube_tet_mesh(h::Real)
+    file = joinpath(MESHDIR, "cube_tet_h$(h).msh")
+    isfile(file) && return file
+    gmsh = GridapGmsh.gmsh
+    gmsh.initialize(); gmsh.option.setNumber("General.Terminal", 0)
+    gmsh.model.add("cube")
+    gmsh.model.occ.addBox(0, 0, 0, 1, 1, 1); gmsh.model.occ.synchronize()
+    ents(d) = [e[2] for e in gmsh.model.getEntities(d)]
+    com(t) = gmsh.model.occ.getCenterOfMass(2, t)
+    for (name, x0) in (("fixedup", 0.0), ("topsuf", 1.0))
+        faces = [t for t in ents(2) if abs(com(t)[1] - x0) < 1e-9]
+        gmsh.model.addPhysicalGroup(2, faces, -1, name)
+        for d in (1, 0)                                   # the closure: edges and vertices under the same name
+            cl = unique([abs(e[2]) for e in gmsh.model.getBoundary([(2, t) for t in faces], false, false, d == 0) if e[1] == d])
+            gmsh.model.addPhysicalGroup(d, cl, -1, name)
+        end
+    end
+    gmsh.model.addPhysicalGroup(3, ents(3), -1, "solid")
+    gmsh.option.setNumber("Mesh.MeshSizeMax", h); gmsh.option.setNumber("Mesh.MeshSizeMin", h)
+    gmsh.model.mesh.generate(3)
+    gmsh.write(file); gmsh.finalize()
+    file
+end
+
+"""
+Mechanical cube on a hierarchy of `k` uniform refinements of a coarse mesh (`mech_cube`'s physics and boundary conditions,
+order 1): `hcube_tet(h, k)` from the Gmsh tetrahedra of size `h`, `hcube_hex(n, k)` from n³ Cartesian hexahedra. The case
+carries the hierarchy (`hier`) and the test-space builder (`space`) the h-multigrid solver needs; the Newton sees the fine
+model only.
+"""
+function hcube(coarse, k::Int; utop=[0.0, 0.0, 0.1], λ=10.0)
+    hier = MeshHierarchy(coarse, k + 1)
+    case = mech_case(fine_model(hier), 1, utop; λ)
+    Du = case.D
+    space(m) = TestFESpace(m, ReferenceFE(lagrangian, VectorValue{3,Float64}, 1), Du, conformity=:H1)
+    (; case..., hier, space)
+end
+hcube_tet(h::Real, k::Int, λ::Real=10.0) = hcube(GmshDiscreteModel(cube_tet_mesh(h)), k; λ)
+function hcube_hex(n::Int, k::Int, λ::Real=10.0)
+    m = CartesianDiscreteModel((0.0, 1.0, 0.0, 1.0, 0.0, 1.0), (n, n, n))
+    labels = get_face_labeling(m)
+    add_tag_from_tags!(labels, "fixedup", CartesianTags.face0YZ⁺)
+    add_tag_from_tags!(labels, "topsuf", CartesianTags.face1YZ⁺)
+    hcube(m, k; λ)
+end
+
+"""
+EM plate with constrained electrodes on a hierarchy (S4-B step 5): the 8 × 0.8 × 40 plate as (7n, 2n, 19n) Cartesian
+hexahedra refined `k` times, clamped at z = 0, electrodes on both faces y = 0 (φ = 0) and y = 0.8 (φ = φtop) with u_y = 0
+on them (closed Cartesian tags inherited by every level; `em_plate_c`'s mid-plane electrode is a vertex tag that the
+refinement does not inherit). φtop = 0.2 over the whole thickness: the field of `em_plate_c` (0.1 over half of it).
+The case carries the hierarchy and the test-space builders of the two fields (`spaces`).
+"""
+function hem_plate_c(n::Int, k::Int, φtop::Real=0.2)
+    m = CartesianDiscreteModel((0.0, 8.0, 0.0, 0.8, 0.0, 40.0), (7n, 2n, 19n))
+    labels = get_face_labeling(m)
+    add_tag_from_tags!(labels, "fixedup", CartesianTags.faceXY0⁺)
+    add_tag_from_tags!(labels, "botsuf", CartesianTags.faceX0Z⁺)
+    add_tag_from_tags!(labels, "topsuf", CartesianTags.faceX1Z⁺)
+    hier = MeshHierarchy(m, k + 1)
+    model = fine_model(hier)
+    physmodel = ElectroMechModel(IdealDielectric(ε=1.0), NeoHookean3D(λ=10.0, μ=1.0))
+    kin = (Kinematics(Mechano, Solid), Kinematics(Electro, Solid))
+    Ω = Triangulation(model); dΩ = Measure(Ω, 2)
+    Du = BoundaryConditions(BoundaryCondition("fixedup", [0.0, 0.0, 0.0]),
+                            BoundaryCondition("topsuf", [0.0, 0.0, 0.0]; mask=(false, true, false)),
+                            BoundaryCondition("botsuf", [0.0, 0.0, 0.0]; mask=(false, true, false)))
+    Dφ = BoundaryConditions(BoundaryCondition("botsuf", 0.0, Λ -> Λ), BoundaryCondition("topsuf", φtop, Λ -> Λ))
+    D = MultiFieldBC([Du, Dφ])
+    space_u(mm) = TestFESpace(mm, ReferenceFE(lagrangian, VectorValue{3,Float64}, 1), Du, conformity=:H1)
+    space_φ(mm) = TestFESpace(mm, ReferenceFE(lagrangian, Float64, 1), Dφ, conformity=:H1)
+    Vu, Vφ = space_u(model), space_φ(model)
+    V = MultiFieldFESpace([Vu, Vφ])
+    U = MultiFieldFESpace([TrialFESpace(Vu, Du, 1.0), TrialFESpace(Vφ, Dφ, 1.0)])
+    res(Λ) = ((u, φ), (v, vφ)) -> HyperFEM.residual(physmodel, kin, (u, φ), (v, vφ), dΩ)
+    jac(Λ) = ((u, φ), (du, dφ), (v, vφ)) -> HyperFEM.jacobian(physmodel, kin, (u, φ), (du, dφ), (v, vφ), dΩ)
+    (; res, jac, U, V, D, nls_kw=(maxiter=20, atol=1e-10, rtol=1e-8), nsteps=5, extra=(;), hier, spaces=(space_u, space_φ))
+end
+
 # Periodic unit cell of a magneto-active composite (the "periodic cell" study of 2026-09, magnetoactive_unit_cell_2D.jl,
 # adapted to the BoundaryConditions API and to both backends; 3D by the same construction): soft magnetically inert
 # matrix (Mooney–Rivlin + IdealMagnetic, χ = 0), rigid inclusion (×100) with remanent magnetisation Mr along x
@@ -183,16 +265,19 @@ function hmre_cell(dim::Int, h::Real, order::Int=1, R::Real=(dim == 2 ? 0.175 :
     (; res, jac, U, V, D, nls_kw=(maxiter=20, atol=1e-10, rtol=1e-8), nsteps=10, extra=(;))
 end
 
-"Case from command-line words: em_plate n | em_box nx ny nz | hmre r [rv] | mech_plate n | mech_box nx ny nz | hmre_cell dim h [order] | hmre_cell_c dim h order contrast."
+"Case from command-line words: em_plate n | em_plate_c n | hem_plate_c n k | em_box nx ny nz | hmre r [rv] | mech_plate n | mech_box nx ny nz | hmre_cell dim h [order] | hmre_cell_c dim h order contrast."
 function build_case(words)
     kind = words[1]
     kind == "hmre_cell_c" && return hmre_cell(parse(Int, words[2]), parse(Float64, words[3]), parse(Int, words[4]); contrast=parse(Float64, words[5])), "hmre_cell$(words[2])d_h$(words[3])_p$(words[4])_c$(words[5])"
     kind == "hmre_cell" && return hmre_cell(parse(Int, words[2]), parse(Float64, words[3]), parse(Int, get(words, 4, "1"))), "hmre_cell$(words[2])d_h$(words[3])_p$(get(words, 4, "1"))"
     kind == "mech_plate" && return mech_plate(parse(Int, words[2])), "mech_plate_n$(words[2])"
     kind == "mech_box" && return mech_box(parse.(Int, words[2:4])..., parse(Int, get(words, 5, "2"))), "mech_box_$(join(words[2:4], 'x'))" * (get(words, 5, "2") == "2" ? "" : "_p$(words[5])")
+    kind == "hcube_tet" && return hcube_tet(parse(Float64, words[2]), parse(Int, words[3]), parse(Float64, get(words, 4, "10"))), "hcube_tet_h$(words[2])_k$(words[3])" * (length(words) ≥ 4 ? "_lam$(words[4])" : "")
+    kind == "hcube_hex" && return hcube_hex(parse(Int, words[2]), parse(Int, words[3]), parse(Float64, get(words, 4, "10"))), "hcube_hex_$(words[2])_k$(words[3])" * (length(words) ≥ 4 ? "_lam$(words[4])" : "")
     kind == "mech_cube" && return mech_cube(parse(Int, words[2]), parse(Int, get(words, 3, "2"))), "mech_cube_$(words[2])_p$(get(words, 3, "2"))"
     kind == "em_plate" && return em_plate(parse(Int, words[2])), "em_plate_n$(words[2])"
     kind == "em_plate_c" && return em_plate_c(parse(Int, words[2])), "em_plate_c_n$(words[2])"
+    kind == "hem_plate_c" && return hem_plate_c(parse(Int, words[2]), parse(Int, words[3])), "hem_plate_c_n$(words[2])_k$(words[3])"
     kind == "em_box" && return em_box(parse.(Int, words[2:4])...), "em_box_$(join(words[2:4], 'x'))"
     kind == "hmre" && return hmre(parse(Float64, words[2]), parse(Float64, get(words, 3, "2"))), "hmre_r$(words[2])"
     error("unknown case $kind")

benchmark/solvers/krylov.jl

@@ -121,7 +121,33 @@ function make_solver(config::String, N, rtol, U, bs, opts, restart=50; case=noth
         pw == "fs" && !(U isa MultiFieldFESpace) && error("mf_*_fs needs a multi-field case")
         return MatrixFreeSolver(; krylov, rtol, maxiter=300, restart, preconditioner=prec)
     end
-    if parts[1] == "fs"                        # S2: fs_<split>[_<fact>_<prec>][_<ksp>][_lagN][_lublocks]
+    if parts[1] == "hmg"                       # S4-B: hmg_<krylov>[_lagN][_degK][_rngA][_lu|_mumps]  (coarse GAMG V-cycle by default)
+        hasproperty(case, :hier) || error("hmg_* needs a hierarchy case (hcube_tet, hcube_hex)")
+        krylov = Symbol(parts[2])
+        lag = findfirst(startswith("lag"), parts); reuse = lag === nothing ? 1 : parse(Int, parts[lag][4:end])
+        deg = findfirst(startswith("deg"), parts); degree = deg === nothing ? 2 : parse(Int, parts[deg][4:end])
+        rng = findfirst(startswith("rng"), parts); range = rng === nothing ? 20.0 : parse(Float64, parts[rng][4:end])
+        Uc = TrialFESpace(case.space(level_model(case.hier, 1)))
+        coarse = "lu" in parts ? LUSolver() : "mumps" in parts ? PETScSolver(:preonly, :mumps) :
+                 PETScSolver(:preonly, :gamg; nullspace=near_nullspace(Uc), blocksize=bs, options=opts)
+        return HMultigridSolver(case.hier, case.space; krylov, rtol, maxiter=300, restart, reuse, degree, range, coarse)
+    end
+    if parts[1] == "fshmg"                     # S4-B step 5: fshmg_<split>[_lagN][_lu]  (one h-multigrid V-cycle per block, coarse GAMG by default)
+        hasproperty(case, :hier) && hasproperty(case, :spaces) || error("fshmg_* needs a multi-field hierarchy case (hem_plate_c)")
+        split_ = Symbol(parts[2])
+        lag = findfirst(startswith("lag"), parts); reuse = lag === nothing ? 1 : parse(Int, parts[lag][4:end])
+        blocks = map(case.spaces) do space
+            Vc = space(level_model(case.hier, 1))
+            Uc = TrialFESpace(Vc)
+            bl = HyperFEM.Solvers.nodal_blocks(Vc)
+            ncomp = allequal(length, bl) ? length(bl[1]) : 1         # 3 for an unmasked displacement, 1 with masked components or for the potential
+            coarse = "lu" in parts ? LUSolver() :
+                     PETScSolver(:preonly, :gamg; nullspace=near_nullspace(Uc), blocksize=ncomp, options=opts)
+            HMultigridSolver(case.hier, space; krylov=:cg, reuse, coarse)
+        end
+        return FieldSplitSolver(U, blocks; split=split_, rtol, maxiter=300, restart)
+    end
+    if parts[1] == "fs"                        # S2: fs_<split>[_<fact>_<prec>][_<ksp>][_lublocks]
         U isa MultiFieldFESpace || error("fs_* needs a multi-field case")
         split_ = Symbol(parts[2])
         fact = length(parts) ≥ 4 && parts[3] in ("upper", "lower", "full", "diag") ? Symbol(parts[3]) : :upper

benchmark/solvers/run_s4b.sh

@@ -0,0 +1,21 @@
+#!/bin/bash
+# S4-B (design/solvers-s4b-plan.md, step 3): h-multigrid against CG + GAMG on the same fine mesh, through krylov.jl.
+#   bash benchmark/solvers/run_s4b.sh tet H K CFG...    # hcube_tet H K (Gmsh tetrahedra of size H refined K times)
+#   bash benchmark/solvers/run_s4b.sh hex N K CFG...    # hcube_hex N K (n³ Cartesian hexahedra refined K times)
+# Stop between two runs with:  touch benchmark/results/STOP
+cd "$(dirname "$0")/../.." || exit 1
+MODE=$1; shift
+T=${T:-6}
+case "$MODE" in
+  tet) H=$1; K=$2; shift 2; CASE="hcube_tet $H $K" ;;
+  hex) N=$1; K=$2; shift 2; CASE="hcube_hex $N $K" ;;
+  tetl) H=$1; K=$2; L=$3; shift 3; CASE="hcube_tet $H $K $L" ;;     # λ of the neo-Hookean law (10: ν = 0.45; 1.5: ν = 0.3)
+  hexl) N=$1; K=$2; L=$3; shift 3; CASE="hcube_hex $N $K $L" ;;
+  em) N=$1; K=$2; shift 2; CASE="hem_plate_c $N $K" ;;                # step 5: EM plate (u, φ), fs_additive vs fshmg_additive
+  *) echo "unknown mode $MODE"; exit 1 ;;
+esac
+for cfg in "$@"; do
+  [ -f benchmark/results/STOP ] && rm -f benchmark/results/STOP && echo "stopped" && exit 0
+  echo "=== $CASE compiled $cfg threads=$T  $(date +%H:%M:%S)"
+  julia -t $T --project=benchmark/solvers benchmark/solvers/krylov.jl $CASE compiled $cfg bs=3 2>&1 | grep --line-buffered -v "^Info\|Warning\|└\|┌"
+done