|
| 1 | +""" |
| 2 | + EisenstatWalkerForcing2(; η₀ = 0.5, ηₘₐₓ = 0.9, γ = 0.9, α = 2, safeguard = true, safeguard_threshold = 0.1) |
| 3 | +
|
| 4 | +Algorithm 2 from the classical work by Eisenstat and Walker (1996) as described by formula (2.6): |
| 5 | + ηₖ = γ * (||rₖ|| / ||rₖ₋₁||)^α |
| 6 | +
|
| 7 | +Here the variables denote: |
| 8 | + rₖ residual at iteration k |
| 9 | + η₀ ∈ [0,1) initial value for η |
| 10 | + ηₘₐₓ ∈ [0,1) maximum value for η |
| 11 | + γ ∈ [0,1) correction factor |
| 12 | + α ∈ [1,2) correction exponent |
| 13 | +
|
| 14 | +Furthermore, the proposed safeguard is implemented: |
| 15 | + ηₖ = max(ηₖ, γ*ηₖ₋₁^α) if γ*ηₖ₋₁^α > safeguard_threshold |
| 16 | +to prevent ηₖ from shrinking too fast. |
| 17 | +""" |
| 18 | +@concrete struct EisenstatWalkerForcing2 |
| 19 | + η₀ |
| 20 | + ηₘₐₓ |
| 21 | + γ |
| 22 | + α |
| 23 | + safeguard |
| 24 | + safeguard_threshold |
| 25 | +end |
| 26 | + |
| 27 | +function EisenstatWalkerForcing2(; η₀ = 0.5, ηₘₐₓ = 0.9, γ = 0.9, α = 2, safeguard = true, safeguard_threshold = 0.1) |
| 28 | + EisenstatWalkerForcing2(η₀, ηₘₐₓ, γ, α, safeguard, safeguard_threshold) |
| 29 | +end |
| 30 | + |
| 31 | + |
| 32 | +@concrete mutable struct EisenstatWalkerForcing2Cache |
| 33 | + p::EisenstatWalkerForcing2 |
| 34 | + η |
| 35 | + rnorm |
| 36 | + rnorm_prev |
| 37 | + internalnorm |
| 38 | + verbosity |
| 39 | +end |
| 40 | + |
| 41 | + |
| 42 | + |
| 43 | +function pre_step_forcing!(cache::EisenstatWalkerForcing2Cache, descend_cache::NonlinearSolveBase.NewtonDescentCache, J, u, fu, iter) |
| 44 | + @SciMLMessage("Eisenstat-Walker forcing residual norm $(cache.rnorm) with rate estimate $(cache.rnorm / cache.rnorm_prev).", cache.verbosity, :forcing) |
| 45 | + |
| 46 | + # On the first iteration we initialize η with the default initial value and stop. |
| 47 | + if iter == 0 |
| 48 | + cache.η = cache.p.η₀ |
| 49 | + @SciMLMessage("Eisenstat-Walker initial iteration to η=$(cache.η).", cache.verbosity, :forcing) |
| 50 | + LinearSolve.update_tolerances!(descend_cache.lincache; reltol=cache.η) |
| 51 | + return nothing |
| 52 | + end |
| 53 | + |
| 54 | + # Store previous |
| 55 | + ηprev = cache.η |
| 56 | + |
| 57 | + # Formula (2.6) |
| 58 | + # ||r|| > 0 should be guaranteed by the convergence criterion |
| 59 | + (; rnorm, rnorm_prev) = cache |
| 60 | + (; α, γ) = cache.p |
| 61 | + cache.η = γ * (rnorm / rnorm_prev)^α |
| 62 | + |
| 63 | + # Safeguard 2 to prevent over-solving |
| 64 | + if cache.p.safeguard |
| 65 | + ηsg = γ*ηprev^α |
| 66 | + if ηsg > cache.p.safeguard_threshold && ηsg > cache.η |
| 67 | + cache.η = ηsg |
| 68 | + end |
| 69 | + end |
| 70 | + |
| 71 | + # Far away from the root we also need to respect η ∈ [0,1) |
| 72 | + cache.η = clamp(cache.η, 0.0, cache.p.ηₘₐₓ) |
| 73 | + |
| 74 | + @SciMLMessage("Eisenstat-Walker iter $iter update to η=$(cache.η).", cache.verbosity, :forcing) |
| 75 | + |
| 76 | + # Communicate new relative tolerance to linear solve |
| 77 | + LinearSolve.update_tolerances!(descend_cache.lincache; reltol=cache.η) |
| 78 | + |
| 79 | + return nothing |
| 80 | +end |
| 81 | + |
| 82 | + |
| 83 | + |
| 84 | +function post_step_forcing!(cache::EisenstatWalkerForcing2Cache, J, u, fu, δu, iter) |
| 85 | + # Cache previous residual norm |
| 86 | + cache.rnorm_prev = cache.rnorm |
| 87 | + cache.rnorm = cache.internalnorm(fu) |
| 88 | + |
| 89 | + # @SciMLMessage("Eisenstat-Walker sanity check: $(cache.internalnorm(fu + J*δu)) ≤ $(cache.η * cache.internalnorm(fu)).", cache.verbosity, :linear_verbosity) |
| 90 | +end |
| 91 | + |
| 92 | + |
| 93 | + |
| 94 | +function InternalAPI.init( |
| 95 | + prob::AbstractNonlinearProblem, alg::EisenstatWalkerForcing2, f, fu, u, p, |
| 96 | + args...; verbose, internalnorm::F = L2_NORM, kwargs... |
| 97 | +) where {F} |
| 98 | + fu_norm = internalnorm(fu) |
| 99 | + |
| 100 | + return EisenstatWalkerForcing2Cache( |
| 101 | + alg, alg.η₀, fu_norm, fu_norm, internalnorm, verbose |
| 102 | + ) |
| 103 | +end |
| 104 | + |
| 105 | + |
| 106 | + |
| 107 | +function InternalAPI.reinit!( |
| 108 | + cache::EisenstatWalkerForcing2Cache; p = cache.p, kwargs... |
| 109 | +) |
| 110 | + cache.p = p |
| 111 | +end |
0 commit comments