Add SIAMFANLEquations fixed point solvers

Signed-off-by: ErikQQY <2283984853@qq.com>

Author: ErikQQY

ext/NonlinearSolveSIAMFANLEquationsExt.jl

@@ -24,7 +24,7 @@ end
 # pseudo transient continuation has a fixed cost per iteration, iteration statistics are
 # not interesting here.
 @inline function __siam_fanl_equations_stats_mapping(method, sol)
-    method === :pseudotransient && return nothing
+    (method === :pseudotransient) || (method === :anderson) && return nothing
     return SciMLBase.NLStats(sum(sol.stats.ifun), sum(sol.stats.ijac), 0, 0,
         sum(sol.stats.iarm))
 end
@@ -36,16 +36,14 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, arg
     @assert (termination_condition ===
              nothing)||(termination_condition isa AbsNormTerminationMode) "SIAMFANLEquationsJL does not support termination conditions!"
 
-    (; method, delta, linsolve) = alg
-
-    iip = SciMLBase.isinplace(prob)
+    (; method, delta, linsolve, m, beta) = alg
 
     T = eltype(prob.u0)
     atol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, T)
     rtol = NonlinearSolve.DEFAULT_TOLERANCE(reltol, T)
 
     if prob.u0 isa Number
-        f = (u) -> prob.f(u, prob.p)
+        f = method == :anderson ? (du, u) -> (du = prob.f(u, prob.p)) : ((u) -> prob.f(u, prob.p))
 
         if method == :newton
             sol = nsolsc(f, prob.u0; maxit = maxiters, atol, rtol, printerr = ShT)
@@ -54,17 +52,22 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, arg
                 printerr = ShT)
         elseif method == :secant
             sol = secant(f, prob.u0; maxit = maxiters, atol, rtol, printerr = ShT)
+        elseif method == :anderson
+            f, u = NonlinearSolve.__construct_f(prob; alias_u0,
+                make_fixed_point = Val(true), can_handle_arbitrary_dims = Val(true))
+            sol = aasol(f, [prob.u0], m, __zeros_like(u, 1, 2*m+4); maxit = maxiters,
+                atol, rtol, beta = beta)
         end
 
         retcode = __siam_fanl_equations_retcode_mapping(sol)
         stats = __siam_fanl_equations_stats_mapping(method, sol)
-        resid = NonlinearSolve.evaluate_f(prob, sol.solution)
+        resid = NonlinearSolve.evaluate_f(prob, sol.solution[1])
         return SciMLBase.build_solution(prob, alg, sol.solution, resid; retcode,
             stats, original = sol)
     end
 
     f!, u = NonlinearSolve.__construct_f(prob; alias_u0,
-        can_handle_arbitrary_dims = Val(true))
+        can_handle_arbitrary_dims = Val(true), make_fixed_point = Val(true))
 
     # Allocate ahead for function
     N = length(u)
@@ -104,6 +107,10 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, arg
         elseif method == :pseudotransient
             sol = ptcsol(f!, u, FS, FPS; atol, rtol, maxit = maxiters,
                 delta0 = delta, printerr = ShT)
+        elseif method == :anderson
+            println(u[:])
+            sol = aasol(f!, u[:], m, zeros(eltype(u), N, 2*m+4), atol, rtol,
+                maxit = maxiters, beta = beta)
         end
     else
         AJ!(J, u, x) = prob.f.jac(J, x, prob.p)

src/extension_algs.jl

@@ -377,23 +377,30 @@ end
   - `method`: the choice of method for solving the nonlinear system.
   - `delta`: initial pseudo time step, default is 1e-3.
   - `linsolve` : JFNK linear solvers, choices are `gmres` and `bicgstab`.
+  - `m`: Depth for Anderson acceleration, default as 0 for Picard iteration
+  - `beta`: Anderson mixing parameter, change f(x) to (1-beta)x+beta*f(x),
+    equivalent to accelerating damped Picard iteration.
 
 ### Submethod Choice
 
   - `:newton`: Classical Newton method.
   - `:pseudotransient`: Pseudo transient method.
   - `:secant`: Secant method for scalar equations.
+  - `:anderson`: Anderson acceleration for fixed point iterations.
 """
 @concrete struct SIAMFANLEquationsJL{L <: Union{Symbol, Nothing}} <:
                  AbstractNonlinearSolveAlgorithm
     method::Symbol
     delta
     linsolve::L
+    m::Int
+    beta
 end
 
-function SIAMFANLEquationsJL(; method = :newton, delta = 1e-3, linsolve = nothing)
+function SIAMFANLEquationsJL(; method = :newton, delta = 1e-3, linsolve = nothing, m = 0,
+        beta = 1.0)
     if Base.get_extension(@__MODULE__, :NonlinearSolveSIAMFANLEquationsExt) === nothing
         error("SIAMFANLEquationsJL requires SIAMFANLEquations.jl to be loaded")
     end
-    return SIAMFANLEquationsJL(method, delta, linsolve)
+    return SIAMFANLEquationsJL(method, delta, linsolve, m, beta)
 end

test/wrappers/fixedpoint.jl

@@ -1,4 +1,4 @@
-using NonlinearSolve, FixedPointAcceleration, SpeedMapping, NLsolve, LinearAlgebra, Test
+using NonlinearSolve, FixedPointAcceleration, SpeedMapping, NLsolve, SIAMFANLEquations, LinearAlgebra, Test
 
 # Simple Scalar Problem
 @testset "Simple Scalar Problem" begin
@@ -14,6 +14,7 @@ using NonlinearSolve, FixedPointAcceleration, SpeedMapping, NLsolve, LinearAlgeb
     @test abs(solve(prob, SpeedMappingJL(; stabilize = true)).resid) ≤ 1e-10
 
     @test abs(solve(prob, NLsolveJL(; method = :anderson)).resid) ≤ 1e-10
+    @test abs(solve(prob, SIAMFANLEquationsJL(; method = :anderson)).resid) ≤ 1e-10
 end
 
 # Simple Vector Problem
@@ -28,6 +29,7 @@ end
     @test maximum(abs.(solve(prob, SpeedMappingJL()).resid)) ≤ 1e-10
     @test maximum(abs.(solve(prob, SpeedMappingJL(; orders = [3, 2])).resid)) ≤ 1e-10
     @test maximum(abs.(solve(prob, SpeedMappingJL(; stabilize = true)).resid)) ≤ 1e-10
+    @test maximum(abs.(solve(prob, SIAMFANLEquationsJL(; method = :anderson)).resid)) ≤ 1e-10
 
     @test_broken maximum(abs.(solve(prob, NLsolveJL(; method = :anderson)).resid)) ≤ 1e-10
 end
@@ -67,6 +69,9 @@ end
     sol = solve(prob, NLsolveJL(; method = :anderson))
     @test_broken sol.u' * A[:, 3] ≈ 32.916472867168096
 
+    sol = solve(prob, SIAMFANLEquationsJL(; method = :anderson))
+    @test sol.u' * A[:, 3] ≈ 32.916472867168096
+
     # Non vector inputs
     function power_method_nonvec!(du, u, A)
         mul!(vec(du), A, vec(u))