diff --git a/Project.toml b/Project.toml index 058c7ae4a..865fc43ff 100644 --- a/Project.toml +++ b/Project.toml @@ -83,7 +83,7 @@ NonlinearProblemLibrary = "0.1.2" OrdinaryDiffEq = "6.63" Pkg = "1" PrecompileTools = "1.2" -Printf = "1.9" +Printf = "1.10" Random = "1.91" RecursiveArrayTools = "3.2" Reexport = "1.2" @@ -92,7 +92,7 @@ SafeTestsets = "0.1" SciMLBase = "2.11" SciMLOperators = "0.3.7" SimpleNonlinearSolve = "1.0.2" -SparseArrays = "1.9" +SparseArrays = "1.10" SparseDiffTools = "2.14" SpeedMapping = "0.3" StableRNGs = "1" diff --git a/docs/src/solvers/FixedPointSolvers.md b/docs/src/solvers/FixedPointSolvers.md index ed507098c..704561905 100644 --- a/docs/src/solvers/FixedPointSolvers.md +++ b/docs/src/solvers/FixedPointSolvers.md @@ -47,3 +47,7 @@ In our tests, we have found the anderson method implemented here to NOT be the m robust. - `NLsolveJL(; method = :anderson)`: Anderson acceleration for fixed point problems. + +### SIAMFANLEquations.jl + + - `SIAMFANLEquationsJL(; method = :anderson)`: Anderson acceleration for fixed point problems. \ No newline at end of file diff --git a/ext/NonlinearSolveSIAMFANLEquationsExt.jl b/ext/NonlinearSolveSIAMFANLEquationsExt.jl index a5bbd5e1b..6ebe82268 100644 --- a/ext/NonlinearSolveSIAMFANLEquationsExt.jl +++ b/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,11 +52,16 @@ 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 @@ -104,6 +107,11 @@ 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 + f!, u = NonlinearSolve.__construct_f(prob; alias_u0, + can_handle_arbitrary_dims = Val(true), make_fixed_point = Val(true)) + sol = aasol(f!, u, m, zeros(T, N, 2*m+4), atol = atol, rtol = rtol, + maxit = maxiters, beta = beta) end else AJ!(J, u, x) = prob.f.jac(J, x, prob.p) diff --git a/src/extension_algs.jl b/src/extension_algs.jl index 195cced8d..8d7397ea0 100644 --- a/src/extension_algs.jl +++ b/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 diff --git a/test/wrappers/fixedpoint.jl b/test/wrappers/fixedpoint.jl index 56b646378..a66a2eaba 100644 --- a/test/wrappers/fixedpoint.jl +++ b/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))