diff --git a/benchmark/solvers/cases.jl b/benchmark/solvers/cases.jl index a6f452c..7f493fb 100644 --- a/benchmark/solvers/cases.jl +++ b/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,7 +197,7 @@ 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⁺) @@ -200,7 +205,11 @@ function hem_plate_c(n::Int, k::Int, φtop::Real=0.2) 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 diff --git a/benchmark/solvers/krylov.jl b/benchmark/solvers/krylov.jl index cb7840b..29813b2 100644 --- a/benchmark/solvers/krylov.jl +++ b/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,7 +179,7 @@ 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 @@ -182,15 +187,24 @@ function main(args) 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)) diff --git a/benchmark/solvers/run_s4b.sh b/benchmark/solvers/run_s4b.sh index 4278050..bd9c331 100755 --- a/benchmark/solvers/run_s4b.sh +++ b/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 diff --git a/src/pmultigrid.jl b/src/pmultigrid.jl index a5121ad..0492a19 100644 --- a/src/pmultigrid.jl +++ b/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) diff --git a/src/smoothers.jl b/src/smoothers.jl index 38f27f8..a72428b 100644 --- a/src/smoothers.jl +++ b/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 diff --git a/src/solvers.jl b/src/solvers.jl index 29b27d7..033da93 100644 --- a/src/solvers.jl +++ b/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." diff --git a/test/solvers/test_pmultigrid.jl b/test/solvers/test_pmultigrid.jl index 717f1c7..b52f623 100644 --- a/test/solvers/test_pmultigrid.jl +++ b/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