From 322725a45950766c47bb206c9e57287aa7f513e4 Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Thu, 28 Dec 2023 15:32:58 +0800 Subject: [PATCH 1/7] Add SIAMFANLEquations fixed point solvers Signed-off-by: ErikQQY <2283984853@qq.com> --- ext/NonlinearSolveSIAMFANLEquationsExt.jl | 21 ++++++++++++++------- src/extension_algs.jl | 11 +++++++++-- test/wrappers/fixedpoint.jl | 7 ++++++- 3 files changed, 29 insertions(+), 10 deletions(-) diff --git a/ext/NonlinearSolveSIAMFANLEquationsExt.jl b/ext/NonlinearSolveSIAMFANLEquationsExt.jl index a5bbd5e1b..54cdb2e41 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,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) diff --git a/src/extension_algs.jl b/src/extension_algs.jl index 195cced8d..ad4c02e2f 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)) From 9f461f18b6afa7e876abd285080059186dbeebef Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Thu, 28 Dec 2023 15:44:18 +0800 Subject: [PATCH 2/7] Downgrade CI Signed-off-by: ErikQQY <2283984853@qq.com> --- Project.toml | 2 +- ext/NonlinearSolveSIAMFANLEquationsExt.jl | 3 +-- src/extension_algs.jl | 2 +- 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 058c7ae4a..13b9da8cf 100644 --- a/Project.toml +++ b/Project.toml @@ -92,7 +92,7 @@ SafeTestsets = "0.1" SciMLBase = "2.11" SciMLOperators = "0.3.7" SimpleNonlinearSolve = "1.0.2" -SparseArrays = "1.9" +SparseArrays = "1.11.0" SparseDiffTools = "2.14" SpeedMapping = "0.3" StableRNGs = "1" diff --git a/ext/NonlinearSolveSIAMFANLEquationsExt.jl b/ext/NonlinearSolveSIAMFANLEquationsExt.jl index 54cdb2e41..32ed3cbc4 100644 --- a/ext/NonlinearSolveSIAMFANLEquationsExt.jl +++ b/ext/NonlinearSolveSIAMFANLEquationsExt.jl @@ -108,8 +108,7 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, arg 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, + sol = aasol(f!, u, m, zeros(eltype(u), N, 2*m+4), atol, rtol, maxit = maxiters, beta = beta) end else diff --git a/src/extension_algs.jl b/src/extension_algs.jl index ad4c02e2f..8d7397ea0 100644 --- a/src/extension_algs.jl +++ b/src/extension_algs.jl @@ -377,7 +377,7 @@ 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 + - `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. From 9dfc6b163d04d6a93d114efb704da6bebd4c381c Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Thu, 28 Dec 2023 15:45:55 +0800 Subject: [PATCH 3/7] Downgrade CI Signed-off-by: ErikQQY <2283984853@qq.com> --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 13b9da8cf..57697e370 100644 --- a/Project.toml +++ b/Project.toml @@ -92,7 +92,7 @@ SafeTestsets = "0.1" SciMLBase = "2.11" SciMLOperators = "0.3.7" SimpleNonlinearSolve = "1.0.2" -SparseArrays = "1.11.0" +SparseArrays = "1.10" SparseDiffTools = "2.14" SpeedMapping = "0.3" StableRNGs = "1" From 61355db2f3fc354b69790c4495958527b26b5567 Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Thu, 28 Dec 2023 16:05:41 +0800 Subject: [PATCH 4/7] should specify atol and rtol Signed-off-by: ErikQQY <2283984853@qq.com> --- ext/NonlinearSolveSIAMFANLEquationsExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/NonlinearSolveSIAMFANLEquationsExt.jl b/ext/NonlinearSolveSIAMFANLEquationsExt.jl index 32ed3cbc4..7386d41a3 100644 --- a/ext/NonlinearSolveSIAMFANLEquationsExt.jl +++ b/ext/NonlinearSolveSIAMFANLEquationsExt.jl @@ -108,7 +108,7 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, arg sol = ptcsol(f!, u, FS, FPS; atol, rtol, maxit = maxiters, delta0 = delta, printerr = ShT) elseif method == :anderson - sol = aasol(f!, u, m, zeros(eltype(u), N, 2*m+4), atol, rtol, + sol = aasol(f!, u, m, zeros(T, N, 2*m+4), atol = atol, rtol = rtol, maxit = maxiters, beta = beta) end else From 20e6a8ba8a483a9c692db1398215e0bfcd49a395 Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Thu, 28 Dec 2023 16:21:59 +0800 Subject: [PATCH 5/7] Fix construct_f for fixed point case Signed-off-by: ErikQQY <2283984853@qq.com> --- ext/NonlinearSolveSIAMFANLEquationsExt.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/ext/NonlinearSolveSIAMFANLEquationsExt.jl b/ext/NonlinearSolveSIAMFANLEquationsExt.jl index 7386d41a3..dc6b9cade 100644 --- a/ext/NonlinearSolveSIAMFANLEquationsExt.jl +++ b/ext/NonlinearSolveSIAMFANLEquationsExt.jl @@ -67,7 +67,7 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, arg end f!, u = NonlinearSolve.__construct_f(prob; alias_u0, - can_handle_arbitrary_dims = Val(true), make_fixed_point = Val(true)) + can_handle_arbitrary_dims = Val(true)) # Allocate ahead for function N = length(u) @@ -108,6 +108,8 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, arg 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 From ed779680da3243a9aac38dd83667f5cbd65fd39d Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Thu, 28 Dec 2023 16:39:44 +0800 Subject: [PATCH 6/7] Fix stats for pseudotransient Signed-off-by: ErikQQY <2283984853@qq.com> --- Project.toml | 2 +- ext/NonlinearSolveSIAMFANLEquationsExt.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 57697e370..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" diff --git a/ext/NonlinearSolveSIAMFANLEquationsExt.jl b/ext/NonlinearSolveSIAMFANLEquationsExt.jl index dc6b9cade..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) || (method === :anderson) && 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 From a72da2537f4f6dedfea56f92d6014b92ba375772 Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Thu, 28 Dec 2023 17:24:15 +0800 Subject: [PATCH 7/7] Add docs for new solvers Signed-off-by: ErikQQY <2283984853@qq.com> --- docs/src/solvers/FixedPointSolvers.md | 4 ++++ 1 file changed, 4 insertions(+) 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