MatrixFreeSolver: linear_indefinite; S5 benchmark cases

benchmark/solvers/cases.jl

@@ -63,12 +63,16 @@ function em_box(nx::Int, ny::Int, nz::Int)
 end
 
 "hMRE beam (Short_magnet geometry refined by r on the beam, rv elsewhere), mechanics with the magnetic load."
-function hmre(r::Real, rv::Real=2)
+function hmre(r::Real, rv::Real=2; hybrid::Bool=false, Bfactor::Real=1.0)
     model = GmshDiscreteModel(joinpath(MESHDIR, "hmre_short_r$(Float64(r))_rv$(Float64(rv)).msh"))
     params = [10.0, 12480.64495232286, 1.0, 2.0, 5195.545287237134, 0.2602717127043121, 2.0]
     modelmech = NonlinearMooneyRivlin2D(λ=(params[1] + params[2]) * 1e2, μ1=params[1] * 1e3, μ2=params[2], α1=params[3], α2=params[4]) +
                 TransverseIsotropy2D(μ=params[5], α1=params[6], α2=params[7])
-    μ0 = 4π * 1e-7; Banorm = 80.0E-9; αr = 10e-3 / μ0
+    # hybrid = true (S5): the mechanical model in HessianRegularization with a TangentSwitch (the Newton's HybridTangent);
+    # Bfactor scales the applied field (the instability sweep)
+    switch = hybrid ? TangentSwitch() : nothing
+    hybrid && (modelmech = HessianRegularization(modelmech; δ=1e-6, switch))
+    μ0 = 4π * 1e-7; Banorm = 80.0E-9 * Bfactor; αr = 10e-3 / μ0
     reffeu = ReferenceFE(lagrangian, VectorValue{2,Float64}, 1)
     Ωs = Triangulation(model, tags=["Beam"]); dΩs = Measure(Ωs, 2)
     Du = BoundaryConditions(BoundaryCondition("Fixed_Solid", [0.0, 0.0]), BoundaryCondition("Fixed_Solid_x", [0.0, 0.0]; mask=(true, false)))
@@ -91,7 +95,8 @@ function hmre(r::Real, rv::Real=2)
     ∂Ψmag∂F(Fn, N, Ba) = (-αr / μ0) * (Ba ⊗ ((invU(Fn)) * N))
     res(Λ) = (u, v) -> ∫((∇(v)' ⊙ (∂Ψs∂F ∘ (F ∘ (∇(u)'), Nh))))dΩs + Λ * ∫((∂Ψmag∂F ∘ (F ∘ (∇(uh⁻)'), Nh, Bah)) ⊙ (∇(v)'))dΩs
     jac(Λ) = (u, du, v) -> ∫(∇(v)' ⊙ ((∂Ψs∂FF ∘ (F ∘ (∇(u)'), Nh)) ⊙ (∇(du)')))dΩs
-    (; res, jac, U, V, D=Du, nls_kw=(maxiter=20, atol=1e-8, rtol=1e-2), nsteps=20, extra=(xh=uh, xh⁻=uh⁻))
+    extra = hybrid ? (xh=uh, xh⁻=uh⁻, tangent_switch=switch) : (xh=uh, xh⁻=uh⁻)
+    (; res, jac, U, V, D=Du, nls_kw=(maxiter=20, atol=1e-8, rtol=1e-2), nsteps=20, extra, switch)
 end
 
 # Mechanical block alone (the clean S1 measurement, D4 refined 2026-10-30): the EM geometries with the
@@ -192,15 +197,19 @@ on them (closed Cartesian tags inherited by every level; `em_plate_c`'s mid-plan
 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)
+function hem_plate_c(n::Int, k::Int, φtop::Real=0.2; hybrid::Bool=false)
     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))
+    # hybrid = true (S5): the coupled model wrapped in HessianRegularization with a TangentSwitch — the Newton's
+    # HybridTangent projects the total mechanical tangent ∂²Ψ/∂F∂F(F, E), elastic and Maxwell parts, in E
+    switch = hybrid ? TangentSwitch() : nothing
+    em = ElectroMechModel(IdealDielectric(ε=1.0), NeoHookean3D(λ=10.0, μ=1.0))
+    physmodel = hybrid ? HessianRegularization(em; δ=1e-6, switch) : em
     kin = (Kinematics(Mechano, Solid), Kinematics(Electro, Solid))
     Ω = Triangulation(model); dΩ = Measure(Ω, 2)
     Du = BoundaryConditions(BoundaryCondition("fixedup", [0.0, 0.0, 0.0]),
@@ -215,7 +224,8 @@ function hem_plate_c(n::Int, k::Int, φtop::Real=0.2)
     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_φ))
+    extra = hybrid ? (; tangent_switch=switch) : (;)
+    (; res, jac, U, V, D, nls_kw=(maxiter=20, atol=1e-10, rtol=1e-8), nsteps=5, extra, hier, spaces=(space_u, space_φ), switch)
 end
 
 # Periodic unit cell of a magneto-active composite (the "periodic cell" study of 2026-09, magnetoactive_unit_cell_2D.jl,
@@ -265,7 +275,7 @@ 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_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."
+"Case from command-line words: em_plate n | em_plate_c n | hem_plate_c n k [φtop] | hem_plate_ch n k [φtop] (S5: projected mechanical tangent) | 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])"
@@ -277,8 +287,10 @@ function build_case(words)
     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 == "hem_plate_c" && return hem_plate_c(parse(Int, words[2]), parse(Int, words[3]), parse(Float64, get(words, 4, "0.2"))), "hem_plate_c_n$(words[2])_k$(words[3])" * (length(words) ≥ 4 ? "_phi$(words[4])" : "")
+    kind == "hem_plate_ch" && return hem_plate_c(parse(Int, words[2]), parse(Int, words[3]), parse(Float64, get(words, 4, "0.2")); hybrid=true), "hem_plate_ch_n$(words[2])_k$(words[3])" * (length(words) ≥ 4 ? "_phi$(words[4])" : "")
     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])"
+    kind == "hmre" && return hmre(parse(Float64, words[2]), parse(Float64, get(words, 3, "2")); Bfactor=parse(Float64, get(words, 4, "1"))), "hmre_r$(words[2])" * (length(words) ≥ 4 ? "_B$(words[4])" : "")
+    kind == "hmre_h" && return hmre(parse(Float64, words[2]), parse(Float64, get(words, 3, "2")); hybrid=true, Bfactor=parse(Float64, get(words, 4, "1"))), "hmre_h_r$(words[2])" * (length(words) ≥ 4 ? "_B$(words[4])" : "")
     error("unknown case $kind")
 end

benchmark/solvers/krylov.jl

@@ -46,7 +46,12 @@ function parse_args(args)
     ewmax = haskey(kw, "ewmax") ? parse(Float64, kw["ewmax"]) : nothing   # ηmax of the forcing term (_ew); library default if absent
     probe = get(kw, "probe", "0") == "1"          # mf_*: time the matrix-free product vs the SpMV and the preconditioner application
     restart = parse(Int, get(kw, "restart", "50"))   # mf_*: GMRES restart (PETSc's default is 30)
-    words, backend, config, rtol, nsteps, bs, opts, ewmax, probe, restart
+    grow = parse(Float64, get(kw, "grow", "1"))      # LoadStepping: growth of ΔΛ after an easy step (1: never)
+    easy = parse(Int, get(kw, "easy", "10"))         # LoadStepping: a step with ≤ easy Newton iterations is easy
+    atolp = parse(Float64, get(kw, "atolp", "1e-6"))  # HybridTangent: absolute tolerance of the projected mode (S5 cases)
+    nlverbose = parse(Int, get(kw, "verbose", "1"))   # Newton verbosity (2: one row per iteration)
+    maxit = parse(Int, get(kw, "maxit", "0"))         # Newton maxiter (0: the case's)
+    words, backend, config, rtol, nsteps, bs, opts, ewmax, probe, restart, grow, easy, atolp, nlverbose, maxit
 end
 
 "Matrix-free product vs SpMV of the assembled tangent, and one preconditioner application, at the final state (10 repetitions each)."
@@ -174,23 +179,32 @@ function make_solver(config::String, N, rtol, U, bs, opts, restart=50; case=noth
 end
 
 function main(args)
-    words, backend, config, rtol, nsteps_arg, bs, opts, ewmax, probe, restart = parse_args(args)
+    words, backend, config, rtol, nsteps_arg, bs, opts, ewmax, probe, restart, grow, easy, atolp, nlverbose, maxit = parse_args(args)
     t0 = time()
     t_case = @elapsed ((case, tag) = build_case(words))
     nsteps = nsteps_arg > 0 ? nsteps_arg : case.nsteps
     t_null = @elapsed (N = near_nullspace(case.U))
     ew = endswith(config, "_ew")
     ls = make_solver(ew ? config[1:end-3] : config, N, rtol, case.U, bs, opts, restart; case)
     startswith(config, "mf_") && (backend = CompiledBackend(jacobian=:matrixfree))   # the Newton receives the operator
-    nls = NewtonRaphson(ls; case.nls_kw..., verbose=1, forcing=!ew ? nothing : ewmax === nothing ? EisenstatWalker() : EisenstatWalker(ηmax=ewmax, η0=min(1e-2, ewmax)))
+    # S5: a case with a TangentSwitch (hem_plate_ch) runs the hybrid Newton — the projected tangent after a failed line
+    # search or an indefinite tangent reported by the linear solver; atolp= is its absolute tolerance in projected mode
+    tangent = hasproperty(case, :switch) && case.switch !== nothing ? HybridTangent(; atol_projected=atolp) : ConsistentTangent()
+    nlkw = maxit > 0 ? merge(case.nls_kw, (; maxiter=maxit)) : case.nls_kw
+    nls = NewtonRaphson(ls; nlkw..., verbose=nlverbose, tangent,
+                        forcing=!ew ? nothing : ewmax === nothing ? EisenstatWalker() : EisenstatWalker(ηmax=ewmax, η0=min(1e-2, ewmax)))
     t_build = @elapsed (m = StaticNonlinearModel(case.res, case.jac, case.U, case.V, case.D; nls, backend, case.extra...))
     n = length(m.caches.x)
-    @printf("%s  backend=%s  config=%s  rtol=%.0e  bs=%d  opts=\"%s\"  threads=%d  DOFs=%d  nullspace %d modes (%.2f s)  build %.1f s  rss %.2f GB\n",
-            tag, backend isa Symbol ? backend : "compiled/$(backend.jacobian)", config, rtol, bs, opts, Threads.nthreads(), n, size(N, 2), t_null, t_build, rss())
+    @printf("%s  backend=%s  config=%s  rtol=%.0e  bs=%d  opts=\"%s\"  threads=%d  DOFs=%d  nullspace %d modes (%.2f s)  build %.1f s  rss %.2f GB%s%s\n",
+            tag, backend isa Symbol ? backend : "compiled/$(backend.jacobian)", config, rtol, bs, opts, Threads.nthreads(), n, size(N, 2), t_null, t_build, rss(),
+            grow > 1 ? "  grow=$grow easy=$easy" : "", tangent isa HybridTangent ? "  hybrid atolp=$atolp" : "")
     flush(stdout)
-    st = @timed solve!(m; stepping=(nsteps=nsteps, maxbisec=2))
+    st = @timed solve!(m; stepping=LoadStepping(; nsteps, maxbisec=2, grow, easy))
     x, flag = st.value
     log = nls.log
+    # S5: the converged state is saved with the row's tag (comparison of the hybrid Newton with the direct solver)
+    savex = get(ENV, "HFK_SAVE_X", "")
+    isempty(savex) || (mkpath(dirname(savex)); open(f -> write(f, x), savex, "w"); @printf("  state saved: %s (%d dofs, ‖x‖ = %.10e)\n", savex, length(x), norm(x)))
     nit = log.niters; nsolves = log.nsolves
     klin = log.nlinear; kmean = klin / max(nit, 1)
     klast = collect(linear_history(nls))

benchmark/solvers/run_s4b.sh

@@ -12,6 +12,8 @@ case "$MODE" in
   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
+  emphi) N=$1; K=$2; P=$3; shift 3; CASE="hem_plate_c $N $K $P" ;;     # S5: EM plate with the electrode voltage P (pull-in sweep)
+  emh) N=$1; K=$2; P=$3; shift 3; CASE="hem_plate_ch $N $K $P" ;;     # S5: the same with the hybrid Newton (projected mechanical tangent)
   *) echo "unknown mode $MODE"; exit 1 ;;
 esac
 for cfg in "$@"; do

src/pmultigrid.jl

@@ -339,3 +339,5 @@ end
 
 preconditioner_builds(st::PMultigridState) = preconditioner_builds(st.lv)
 preconditioner_ok(st::PMultigridState) = st.ok
+"An indefinite tangent on some level (a non-positive Ritz value of a smoother with finite blocks), the S5 signal."
+preconditioner_indefinite(st::PMultigridState) = !st.ok && levels_ok(st.lv) && any(smoother_indefinite, st.smoothers)

src/smoothers.jl

@@ -39,6 +39,14 @@ function chebyshev_damping(c::Chebyshev)
     1 / cosh(c.degree * acosh(θ / δ))
 end
 
+"""
+Whether the smoother's last interval estimate found an **indefinite** operator: finite nodal blocks and a finite
+Ritz estimate, but a non-positive Ritz value (a limit point, pull-in, buckling). Distinct from a non-finite tangent
+(`ok == false` as well, `smoother_indefinite == false`), which is a plain failure; the Newton with a `HybridTangent`
+switches to the projected tangent on the indefinite signal (S5).
+"""
+smoother_indefinite(st) = !st.ok && st.blocks.ok && isfinite(st.λest) && isfinite(st.λmin) && (st.λmin ≤ 0 || st.λest ≤ 0)
+
 "Smoother state on one level: the operator, its nodal blocks, the interval and the work vectors."
 mutable struct ChebyshevState{MF}
     cfg::Chebyshev

src/solvers.jl

@@ -8,7 +8,7 @@
 # `preconditioners.jl`, applied with `ldiv!`.
 #
 # Hooks read by a Newton solver (HyperFEM's `NewtonRaphson` through its HyperFEMKernels extension):
-# `linear_converged`, `linear_iterations`, `set_tolerance!`, `tolerance`, `preconditioner_builds`.
+# `linear_converged`, `linear_iterations`, `linear_indefinite`, `set_tolerance!`, `tolerance`, `preconditioner_builds`.
 # They are generic functions of this package, not exported.
 
 """
@@ -75,6 +75,7 @@ mutable struct MatrixFreeSolverNS{S<:MatrixFreeSolver,W,P,V} <: Algebra.Numerica
     rtol::Float64            # tolerance in force for the next solves (`set_tolerance!`)
     iters::Int
     converged::Bool
+    indefinite::Bool         # the last solve met an indefinite operator (a negative curvature in CG, an indefinite level)
     rnorm::Float64           # relative residual of the last solve when it had to be checked (NaN otherwise)
     nsetup::Int
 end
@@ -92,7 +93,7 @@ function Algebra.numerical_setup(ss::MatrixFreeSolverSS, mf::MatrixFreeJacobian)
     length(mf.x) == n || error("HyperFEMKernels: the MatrixFreeJacobian has not been prepared (call jacobian!(mf, op, x) before numerical_setup)")
     prec = s.preconditioner === nothing ? nothing : preconditioner_setup(s.preconditioner, mf)
     ws = _workspace(s, n)
-    MatrixFreeSolverNS(s, mf, ws, prec, similar(ws.x), s.rtol, 0, true, NaN, 1)
+    MatrixFreeSolverNS(s, mf, ws, prec, similar(ws.x), s.rtol, 0, true, false, NaN, 1)
 end
 
 function Algebra.numerical_setup!(ns::MatrixFreeSolverNS, mf::MatrixFreeJacobian)
@@ -104,13 +105,18 @@ end
 
 "Whether a preconditioner state can be applied (a p-multigrid with a non-finite level cannot); `true` by default."
 preconditioner_ok(::Any) = true
+"Whether a preconditioner state found an indefinite tangent (the S5 instability signal); `false` by default."
+preconditioner_indefinite(::Any) = false
 
 function Algebra.solve!(x::AbstractVector, ns::MatrixFreeSolverNS, b::AbstractVector)
     s, ws, P = ns.solver, ns.workspace, ns.prec
+    ns.indefinite = false
     if P !== nothing && !preconditioner_ok(P)
-        # e.g. a NaN tangent on a coarse level of the p-multigrid: no Krylov iteration, a failed solve the Newton handles
+        # e.g. a NaN tangent on a coarse level of the p-multigrid: no Krylov iteration, a failed solve the Newton handles;
+        # an indefinite level is reported as such (the Newton with a HybridTangent switches tangent, S5)
         fill!(x, 0.0)
         ns.iters, ns.converged, ns.rnorm = 0, false, NaN
+        ns.indefinite = preconditioner_indefinite(P)
         return x
     end
     ldiv = P !== nothing
@@ -124,7 +130,9 @@ function Algebra.solve!(x::AbstractVector, ns::MatrixFreeSolverNS, b::AbstractVe
         elseif s.krylov === :fgmres
             Krylov.fgmres!(ws, ns.mf, b; N, ldiv, restart=true, kw...)
         else
-            Krylov.cg!(ws, ns.mf, b; M=N, ldiv, kw...)
+            # linesearch = true: a non-positive curvature pᵀAp ≤ 0 stops CG (status "nonpositive curvature detected") with
+            # the last positive-curvature iterate instead of diverging — reported below as a failed, indefinite solve
+            Krylov.cg!(ws, ns.mf, b; M=N, ldiv, linesearch=true, kw...)
         end
     finally
         s.blas_threads > 0 && BLAS.set_num_threads(nblas)
@@ -133,6 +141,11 @@ function Algebra.solve!(x::AbstractVector, ns::MatrixFreeSolverNS, b::AbstractVe
     ns.iters = ws.stats.niter
     ns.converged = ws.stats.solved
     ns.rnorm = NaN
+    if s.krylov === :cg && occursin("curvature", ws.stats.status)
+        ns.converged = false
+        ns.indefinite = true
+        return x
+    end
     if !ns.converged
         # Krylov.jl stops without its residual test on an Arnoldi breakdown (an exact preconditioner
         # solves in one iteration) or at the iteration limit: decide on the true residual
@@ -151,6 +164,8 @@ end
 linear_converged(ns::MatrixFreeSolverNS) = ns.converged
 "Krylov iterations of the last `solve!`."
 linear_iterations(ns::MatrixFreeSolverNS) = ns.iters
+"Whether the last `solve!` met an indefinite operator (CG's non-positive curvature, or an indefinite multigrid level): S5."
+linear_indefinite(ns::MatrixFreeSolverNS) = ns.indefinite
 "Relative tolerance of the next solves (`nothing` restores the solver's own)."
 set_tolerance!(ns::MatrixFreeSolverNS, rtol) = (ns.rtol = rtol === nothing ? ns.solver.rtol : _check_rtol(rtol); ns)
 "Relative tolerance in force."

test/solvers/test_pmultigrid.jl

@@ -127,6 +127,7 @@ end
     jacobian!(mf, op, xb); numerical_setup!(ns, mf)
     solve!(dx, ns, b)
     @test !HK.linear_converged(ns) && HK.linear_iterations(ns) == 0
+    @test !HK.linear_indefinite(ns)                 # non-finite is not the instability signal
     jacobian!(mf, op, x); numerical_setup!(ns, mf)
     solve!(dx, ns, b); @test HK.linear_converged(ns)
 end
@@ -243,4 +244,17 @@ end
     @test ns.prec.smoothers[1].λmin < 0
     dx = zeros(size(A, 1)); solve!(dx, ns, randn(size(A, 1)))
     @test !HK.linear_converged(ns) && HK.linear_iterations(ns) == 0
+    @test HK.linear_indefinite(ns)                  # the S5 signal: HyperFEM's HybridTangent switches tangent on it
+    # without a preconditioner CG itself meets the negative curvature (Krylov.jl's linesearch mode): failed, indefinite,
+    # and the iterate is finite (the last positive-curvature one)
+    ls0 = MatrixFreeSolver(; krylov=:cg, rtol=1e-8, maxiter=500)
+    ns0 = numerical_setup(symbolic_setup(ls0, mf), mf)
+    solve!(dx, ns0, randn(size(A, 1)))
+    @test !HK.linear_converged(ns0) && HK.linear_indefinite(ns0) && all(isfinite, dx)
+    # a definite tangent through the same setups: converged, not indefinite
+    op2 = KernelOperator(c.lres, c.ljac, U, V; jacobian=:matrixfree)       # linear elasticity: SPD
+    x2 = zeros(num_free_dofs(U)); mf2 = allocate_jacobian(op2, x2); jacobian!(mf2, op2, x2)
+    ns2 = numerical_setup(symbolic_setup(ls0, mf2), mf2)
+    solve!(dx, ns2, randn(size(A, 1)))
+    @test HK.linear_converged(ns2) && !HK.linear_indefinite(ns2)
 end