Add nonlinear solver options as keyword argument for integrators.

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"

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...)

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
 
 

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...)

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

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

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
 
 

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

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
 
 

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