Revise integrator interface (again) (#185)

* Implement the new integrator interface for Euler, ERK, IRK, and VPRK methods.

* Implement the new integrator interface for splitting and composition integrators.

* Revise integrator interface in SPARK methods.

* Remove obsolete integrate_step! methods.

* Revise integrator interface in HPI methods.

* Revise integrator interface in DVI methods.

* Revise integrator interface in DIRK methods.

* Revise integrator interface in partitioned Runge-Kutta methods.

* Revise integrator interface in CGVI methods.

* Refactor Runge-Kutta methods for implicit ODEs.

* Revise integrator interface in partitioned Runge-Kutta methods for implicit ODEs.

* Refactor integrator interface and projection methods, and remove storage of previous solution in cache.

* Start to refactor extrapolation interface towards unification with integrator interface.

* Minor fix.

* Add minor comment in docs.

* Refactor integrator, remove solstep, and completely revise initial guess handling and extrapolation methods.

Author: michakraus

Project.toml

@@ -29,8 +29,8 @@ CompactBasisFunctions = "0.2"
 Documenter = "0.25, 0.26, 0.27, 1"
 ForwardDiff = "0.10"
 GenericLinearAlgebra = "0.2, 0.3"
-GeometricBase = "0.10"
-GeometricEquations = "0.16, 0.17"
+GeometricBase = "0.10.11"
+GeometricEquations = "0.18"
 GeometricSolutions = "0.3"
 OffsetArrays = "1"
 Parameters = "0.12"

docs/src/developer/projections.md

@@ -0,0 +1,8 @@
+# Projection Methods
+
+Projection methods are implemented in such a way that they can be combined with almost any one-step method in GeometricIntegrators.
+
+
+
+
+

src/Extrapolators.jl

@@ -7,13 +7,16 @@ module Extrapolators
     using GeometricEquations
 
     import Base: Callable
+    import GeometricBase: solutionstep!
 
     export Extrapolation,
            EulerExtrapolation,
            MidpointExtrapolation,
            HermiteExtrapolation
-    export extrapolate!
+    export NoInitialGuess
+    export extrapolate!, solutionstep!
 
+    include("extrapolation/vectorfields.jl")
     include("extrapolation/extrapolation.jl")
     include("extrapolation/aitken_neville.jl")
     include("extrapolation/euler.jl")

src/Integrators.jl

@@ -29,8 +29,8 @@ module Integrators
     import CompactBasisFunctions: nbasis
 
     import GeometricBase: description, reference, tableau, order
-    import GeometricBase: equations, nconstraints, parameters, timestep
-    import GeometricBase: integrate, integrate!
+    import GeometricBase: equations, initialguess, nconstraints, parameters, timestep
+    import GeometricBase: integrate, integrate!, solutionstep!
     import GeometricBase: reset!
     import GeometricBase.Utils: @big, @define, compensated_summation
 
@@ -97,14 +97,6 @@ module Integrators
     include("solutions/solution_step_constructors.jl")
 
 
-    export InitialGuess, NoInitialGuess
-    export initialguess!
-
-    include("initial_guess/initial_guess.jl")
-    include("initial_guess/hermite.jl")
-    include("initial_guess/midpoint.jl")
-
-
     export IntegratorCache, IntegratorConstructor
     export equation, timestep
 
@@ -117,7 +109,7 @@ module Integrators
 
 
     export GeometricIntegrator
-    export integrate, integrate!, integrate_step!
+    export integrate, integrate!, solutionstep!
 
     include("integrators/integrator.jl")
 

src/SPARK.jl

@@ -22,15 +22,18 @@ module SPARK
     import ..Integrators: GeometricIntegrator, Newton
     import ..Integrators: HDAEMethod, IDAEMethod, LDAEMethod, PDAEMethod
     import ..Integrators: SolutionStepPDAE
-    import ..Integrators: InitialGuess, Extrapolation, HermiteExtrapolation
-    import ..Integrators: initialguess!, initial_guess!, integrate_step!, residual!
+    import ..Integrators: StageVector
+    import ..Integrators: Extrapolation, HermiteExtrapolation
+    import ..Integrators: initial_guess!, integrate_step!, residual!, solutionstep!
     import ..Integrators: CacheDict, Cache, CacheType, IDAEIntegratorCache
     import ..Integrators: AbstractCoefficients,
                           @CoefficientsRK, @HeaderCoefficientsRK
     import ..Integrators: create_internal_stage_vector,
                           update!, update_vector_fields!, update_multiplier!,
-                          initialize!, initsolver, internal, nlsolution
-    import ..Integrators: equation, equations, timestep, eachstage, nstages
+                          initialize!, initsolver, internal, nlsolution,
+                          internal_variables
+    import ..Integrators: current, equation, equations, timestep, eachstage, nstages
+    import ..Integrators: cache, caches, iguess, method, problem, solver
 
 
     export CoefficientsARK, CoefficientsPRK, CoefficientsMRK, CoefficientsIRK,

src/extrapolation/euler.jl

@@ -36,9 +36,10 @@ end
 
 function extrapolate!(t₀::TT, x₀::AbstractVector{DT},
                       t₁::TT, x₁::AbstractVector{DT},
-                      v::Callable, params::OptionalParameters,
+                      problem::AbstractProblemODE,
                       extrap::EulerExtrapolation) where {DT,TT}
-    @assert size(x₀) == size(x₁)
+
+    @assert axes(x₀) == axes(x₁)
 
     local F   = collect(1:(extrap.s+1))
     local σ   = (t₁ - t₀) ./ F
@@ -53,7 +54,7 @@ function extrapolate!(t₀::TT, x₀::AbstractVector{DT},
             for k in axes(pts,1)
                 xᵢ[k] = pts[k,i]
             end
-            v(vᵢ, tᵢ, xᵢ, params)
+            initialguess(problem).v(vᵢ, tᵢ, xᵢ, parameters(problem))
             for k in axes(pts,1)
                 pts[k,i] += σ[i] * vᵢ[k]
             end
@@ -65,9 +66,8 @@ function extrapolate!(t₀::TT, x₀::AbstractVector{DT},
     return x₁
 end
 
-function extrapolate!(t₀, x₀::AbstractVector,
-                      t₁, x₁::AbstractVector,
-                      problem::AbstractProblemODE,
-                      extrap::EulerExtrapolation)
-    extrapolate!(t₀, x₀, t₁, x₁, functions(problem).v, parameters(problem), extrap)
+function solutionstep!(sol, history, problem::Union{AbstractProblemODE, SODEProblem}, extrap::EulerExtrapolation)
+    extrapolate!(history.t[1], history.q[1], sol.t, sol.q, problem, extrap)
+    update_vectorfields!(sol, problem)
+    return sol
 end

src/extrapolation/extrapolation.jl

@@ -2,5 +2,29 @@
 const default_extrapolation_stages = 5
 
 abstract type Extrapolation end
+# abstract type Extrapolation <: DeterministicMethod end
 
-function extrapolate! end
+
+struct NoInitialGuess <: Extrapolation end
+
+
+# """
+
+
+# """
+# function extrapolate! end
+
+
+
+# function extrapolate!(t₀, q₀, q̇₀, t₁, q₁, q̇₁, t₂, q₂, q̇₂, problem::AbstractProblemODE, extrap::Extrapolation; kwargs...)
+#     solution = (
+#         t = t₀, q = q₀, v = q̇₀,
+#     )
+
+#     history = (
+#         (t = t₁, q = q₁, v = q̇₁),
+#         (t = t₂, q = q₂, v = q̇₂),
+#     )
+
+#     solutionstep!(solution, history, problem, extrap; kwargs...)
+# end

src/extrapolation/hermite.jl

@@ -188,45 +188,40 @@ function extrapolate!(t₀::TT, x₀::AbstractArray{DT}, ẋ₀::AbstractArray{D
     return (xᵢ, ẋᵢ)
 end
 
+function solutionstep!(sol, history, problem::Union{AbstractProblemODE, SODEProblem}, extrap::HermiteExtrapolation; nowarn = false)
+    t₀, q₀, q̇₀ = history.t[2], history.q[2], history.v[2]
+    t₁, q₁, q̇₁ = history.t[1], history.q[1], history.v[1]
+
+    if q₀ == q₁
+        nowarn || @warn "Hermite Extrapolation: q's history[1] and history[2] are identical!"
+        sol.q .= q₁
+        sol.v .= q̇₁
+    else
+        extrapolate!(t₀, q₀, q̇₀, t₁, q₁, q̇₁, sol.t, sol.q, sol.v, extrap)
+    end
 
-function _vectorfield(t::TT, x::AbstractArray{DT}, v::Callable, params) where {DT,TT}
-    ẋ = zero(x)
-    v(ẋ, t, x, params)
-    return ẋ
+    return sol
 end
 
-function extrapolate!(t₀::TT, x₀::AbstractArray{DT},
-                      t₁::TT, x₁::AbstractArray{DT},
-                      tᵢ::TT, xᵢ::AbstractArray{DT},
-                      v::Function, params::OptionalParameters,
-                      extrap::HermiteExtrapolation) where {DT,TT}
-    ẋ₀ = _vectorfield(t₀, x₀, v, params)
-    ẋ₁ = _vectorfield(t₁, x₁, v, params)
-    extrapolate!(t₀, x₀, ẋ₀, t₁, x₁, ẋ₁, tᵢ, xᵢ, extrap)
-end
+function solutionstep!(sol, history, problem::Union{AbstractProblemPODE, AbstractProblemIODE}, extrap::HermiteExtrapolation; nowarn = false)
+    t₀, q₀, v₀, p₀, f₀ = history.t[2], history.q[2], history.v[2], history.p[2], history.f[2]
+    t₁, q₁, v₁, p₁, f₁ = history.t[1], history.q[1], history.v[1], history.p[1], history.f[1]
 
-function extrapolate!(t₀::TT, x₀::AbstractArray{DT},
-                      t₁::TT, x₁::AbstractArray{DT},
-                      tᵢ::TT, xᵢ::AbstractArray{DT}, ẋᵢ::AbstractArray{DT}, 
-                      v::Function, params::OptionalParameters,
-                      extrap::HermiteExtrapolation) where {DT,TT}
-    ẋ₀ = _vectorfield(t₀, x₀, v, params)
-    ẋ₁ = _vectorfield(t₁, x₁, v, params)
-    extrapolate!(t₀, x₀, ẋ₀, t₁, x₁, ẋ₁, tᵢ, xᵢ, ẋᵢ, extrap)
-end
+    if q₀ == q₁
+        nowarn || @warn "Hermite Extrapolation: q's history[1] and history[2] are identical!"
+        sol.q .= q₁
+        sol.v .= v₁
+    else
+        extrapolate!(t₀, q₀, v₀, t₁, q₁, v₁, sol.t, sol.q, sol.v, extrap)
+    end
 
-function extrapolate!(t₀::TT, x₀::AbstractArray{DT},
-                      t₁::TT, x₁::AbstractArray{DT},
-                      tᵢ::TT, xᵢ::AbstractArray{DT},
-                      problem::AbstractProblemODE,
-                      extrap::HermiteExtrapolation) where {DT,TT}
-    extrapolate!(t₀, x₀, t₁, x₁, tᵢ, xᵢ, functions(problem).v, parameters(problem), extrap)
-end
+    if p₀ == p₁
+        nowarn || @warn "Hermite Extrapolation: p's history[1] and history[2] are identical!"
+        sol.p .= p₁
+        sol.f .= f₁
+    else
+        extrapolate!(t₀, p₀, f₀, t₁, p₁, f₁, sol.t, sol.p, sol.f, extrap)
+    end
 
-function extrapolate!(t₀::TT, x₀::AbstractArray{DT},
-                      t₁::TT, x₁::AbstractArray{DT},
-                      tᵢ::TT, xᵢ::AbstractArray{DT}, ẋᵢ::AbstractArray{DT}, 
-                      problem::AbstractProblemODE,
-                      extrap::HermiteExtrapolation) where {DT,TT}
-    extrapolate!(t₀, x₀, t₁, x₁, tᵢ, xᵢ, ẋᵢ, functions(problem).v, parameters(problem), extrap)
+    return sol
 end