From d82d9bb295a5f59018761ba0224d80556b9e79dc Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Tue, 3 Dec 2024 16:01:16 +0100 Subject: [PATCH] Add nonlinear solver options as keyword argument for integrators. --- Project.toml | 2 +- src/integrators/integrator.jl | 9 +++++---- src/integrators/rk/integrators_dirk.jl | 4 ++-- src/integrators/solvers.jl | 4 ++-- src/integrators/splitting/composition_integrator.jl | 7 ++++--- src/projections/projection.jl | 3 ++- src/projections/standard_projection.jl | 4 ++-- src/spark/abstract.jl | 4 ++-- src/spark/integrators_slrk.jl | 4 ++-- test/integrators/euler_tests.jl | 9 +++++++++ 10 files changed, 31 insertions(+), 19 deletions(-) diff --git a/Project.toml b/Project.toml index 98d68ffd8..a6ab1f9ae 100644 --- a/Project.toml +++ b/Project.toml @@ -31,7 +31,7 @@ ForwardDiff = "0.10" GenericLinearAlgebra = "0.2, 0.3" GeometricBase = "0.10.11" GeometricEquations = "0.18" -GeometricSolutions = "0.3" +GeometricSolutions = "0.3, 0.4" OffsetArrays = "1" Parameters = "0.12" PrettyTables = "2" diff --git a/src/integrators/integrator.jl b/src/integrators/integrator.jl index 03b3ab531..0df1849a5 100644 --- a/src/integrators/integrator.jl +++ b/src/integrators/integrator.jl @@ -48,9 +48,10 @@ function GeometricIntegrator( integratormethod::GeometricMethod, solvermethod::SolverMethod, iguess::Extrapolation; + options = default_options(), method = initmethod(integratormethod, problem), caches = CacheDict(problem, method), - solver = initsolver(solvermethod, method, caches) + solver = initsolver(solvermethod, options, method, caches) ) GeometricIntegrator(problem, method, caches, solver, iguess) end @@ -190,8 +191,8 @@ function integrate!(sol::GeometricSolution, int::AbstractIntegrator) end -function integrate(integrator::AbstractIntegrator; kwargs...) - solution = Solution(problem(integrator); kwargs...) +function integrate(integrator::AbstractIntegrator) + solution = Solution(problem(integrator)) integrate!(solution, integrator) end @@ -201,7 +202,7 @@ function integrate(problem::AbstractProblem, method::GeometricMethod; kwargs...) end function integrate(problems::EnsembleProblem, method::GeometricMethod; kwargs...) - solutions = Solution(problems; kwargs...) + solutions = Solution(problems) for (problem, solution) in zip(problems, solutions) integrator = GeometricIntegrator(problem, method; kwargs...) diff --git a/src/integrators/rk/integrators_dirk.jl b/src/integrators/rk/integrators_dirk.jl index b682542d9..d95f750db 100644 --- a/src/integrators/rk/integrators_dirk.jl +++ b/src/integrators/rk/integrators_dirk.jl @@ -47,8 +47,8 @@ end Base.getindex(s::SingleStageSolvers, args...) = getindex(s.solvers, args...) -function initsolver(::NewtonMethod, method::DIRK, caches::CacheDict; kwargs...) - SingleStageSolvers([NewtonSolver(zero(cache(caches).x[i]), zero(cache(caches).x[i]); linesearch = Backtracking(), config = Options(min_iterations = 1)) for i in eachstage(method)]...) +function initsolver(::NewtonMethod, config::Options, method::DIRK, caches::CacheDict; kwargs...) + SingleStageSolvers([NewtonSolver(zero(cache(caches).x[i]), zero(cache(caches).x[i]); linesearch = Backtracking(), config = config) for i in eachstage(method)]...) end diff --git a/src/integrators/solvers.jl b/src/integrators/solvers.jl index bbdc38d7f..b6f3c7b32 100644 --- a/src/integrators/solvers.jl +++ b/src/integrators/solvers.jl @@ -7,10 +7,10 @@ default_options() = Options( f_abstol = 8eps(), ) -initsolver(::SolverMethod, ::GeometricMethod, ::CacheDict) = NoSolver() +initsolver(::SolverMethod, ::Options, ::GeometricMethod, ::CacheDict) = NoSolver() # create nonlinear solver -function initsolver(::NewtonMethod, ::GeometricMethod, caches::CacheDict; config = default_options(), kwargs...) +function initsolver(::NewtonMethod, config::Options, ::GeometricMethod, caches::CacheDict; kwargs...) x = zero(nlsolution(caches)) y = zero(nlsolution(caches)) NewtonSolver(x, y; linesearch = Backtracking(), config = config, kwargs...) diff --git a/src/integrators/splitting/composition_integrator.jl b/src/integrators/splitting/composition_integrator.jl index 15a883a17..fcda8881f 100644 --- a/src/integrators/splitting/composition_integrator.jl +++ b/src/integrators/splitting/composition_integrator.jl @@ -41,7 +41,7 @@ methods(c::Composition{<: GeometricMethod}, neqs) = Tuple(method(c) for _ in 1:n splitting(c::Composition) = c.splitting - +_options(methods::Tuple) = Tuple([default_options() for _ in methods]) _solvers(methods::Tuple) = Tuple([default_solver(m) for m in methods]) _iguesses(methods::Tuple) = Tuple([default_iguess(m) for m in methods]) @@ -66,16 +66,17 @@ struct CompositionIntegrator{ problem::SODEProblem, splitting::AbstractSplittingMethod, methods::Tuple; + options = _options(methods), solvers = _solvers(methods), initialguesses = _iguesses(methods)) - @assert length(methods) == length(solvers) == length(initialguesses) == _neqs(problem) + @assert length(methods) == length(options) == length(solvers) == length(initialguesses) == _neqs(problem) # get splitting indices and coefficients f, c = coefficients(problem, splitting) # construct composition integrators - subints = Tuple(GeometricIntegrator(SubstepProblem(problem, c[i], f[i]), methods[f[i]]) for i in eachindex(f,c)) + subints = Tuple(GeometricIntegrator(SubstepProblem(problem, c[i], f[i]), methods[f[i]]; options = options[f[i]], solver = solvers[f[i]]) for i in eachindex(f,c)) new{typeof(splitting), typeof(problem), typeof(subints)}(problem, splitting, subints) end diff --git a/src/projections/projection.jl b/src/projections/projection.jl index c4002a8e0..727db83b5 100644 --- a/src/projections/projection.jl +++ b/src/projections/projection.jl @@ -56,9 +56,10 @@ function ProjectionIntegrator( solvermethod::SolverMethod, iguess::Extrapolation, subint::AbstractIntegrator; + options = default_options(), method = initmethod(projectionmethod, problem), caches = CacheDict(problem, method), - solver = initsolver(solvermethod, method, caches) + solver = initsolver(solvermethod, options, method, caches) ) ProjectionIntegrator(problem, method, caches, solver, iguess, subint) end diff --git a/src/projections/standard_projection.jl b/src/projections/standard_projection.jl index 2dd6a329f..b6cb2d364 100644 --- a/src/projections/standard_projection.jl +++ b/src/projections/standard_projection.jl @@ -43,9 +43,9 @@ function split_nlsolution(x::AbstractVector, int::StandardProjectionIntegrator) end -function initsolver(::NewtonMethod, ::ProjectedMethod{<:StandardProjection}, caches::CacheDict; kwargs...) +function initsolver(::NewtonMethod, config::Options, ::ProjectedMethod{<:StandardProjection}, caches::CacheDict; kwargs...) x̄, x̃ = split_nlsolution(cache(caches)) - NewtonSolver(zero(x̃), zero(x̃); kwargs...) + NewtonSolver(zero(x̃), zero(x̃); config = config, kwargs...) end diff --git a/src/spark/abstract.jl b/src/spark/abstract.jl index 7ad492016..7b010fa0b 100644 --- a/src/spark/abstract.jl +++ b/src/spark/abstract.jl @@ -22,6 +22,6 @@ hasnullvector(method::AbstractSPARKMethod) = hasnullvector(tableau(method)) # create nonlinear solver -function initsolver(::Newton, method::AbstractSPARKMethod, caches::CacheDict) - NewtonSolver(zero(nlsolution(caches)), zero(nlsolution(caches)); linesearch = Backtracking(), config = Options(min_iterations = 1, x_abstol = 8eps(), f_abstol = 8eps())) +function initsolver(::Newton, config::Options, method::AbstractSPARKMethod, caches::CacheDict) + NewtonSolver(zero(nlsolution(caches)), zero(nlsolution(caches)); linesearch = Backtracking(), config = config) end diff --git a/src/spark/integrators_slrk.jl b/src/spark/integrators_slrk.jl index 6cf08ff42..374ce8d38 100644 --- a/src/spark/integrators_slrk.jl +++ b/src/spark/integrators_slrk.jl @@ -95,12 +95,12 @@ F^1_{n,i} + F^2_{n,i} &= \frac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , & i const IntegratorSLRK{DT,TT} = GeometricIntegrator{<:LDAEProblem{DT,TT}, <:SLRK} -# function Integrators.initsolver(::Newton, solstep::SolutionStepPDAE{DT}, problem::LDAEProblem, method::SLRK, caches::CacheDict) where {DT} +# function Integrators.initsolver(::Newton, config::Options, solstep::SolutionStepPDAE{DT}, problem::LDAEProblem, method::SLRK, caches::CacheDict) where {DT} # # create wrapper function f!(b,x) # f! = (b,x) -> residual!(b, x, solstep, problem, method, caches) # # create nonlinear solver -# NewtonSolver(zero(caches[DT].x), zero(caches[DT].x), f!; linesearch = Backtracking(), config = Options(min_iterations = 1, x_abstol = 8eps(), f_abstol = 8eps())) +# NewtonSolver(zero(caches[DT].x), zero(caches[DT].x), f!; linesearch = Backtracking(), config = config) # end diff --git a/test/integrators/euler_tests.jl b/test/integrators/euler_tests.jl index 14f3eb79f..66230df29 100644 --- a/test/integrators/euler_tests.jl +++ b/test/integrators/euler_tests.jl @@ -1,5 +1,6 @@ using GeometricIntegrators using GeometricProblems.HarmonicOscillator +using SimpleSolvers: Options using Test @@ -16,4 +17,12 @@ ref = exact_solution(ode) err = relative_maximum_error(sol, ref) @test err.q < 5E-2 + sol = integrate(ode, ImplicitEuler(); options = Options( + min_iterations = 1, + x_abstol = 2eps(), + f_abstol = 2eps(), + )) + err = relative_maximum_error(sol, ref) + @test err.q < 5E-2 + end