Merge pull request #71 from michakraus/master

v0.4.1 bugfix and feature release

Project.toml

@@ -1,7 +1,7 @@
 name = "GeometricIntegrators"
 uuid = "dcce2d33-59f6-5b8d-9047-0defad88ae06"
 authors = ["Michael Kraus <[email protected]>"]
-version = "0.4.0"
+version = "0.4.1"
 
 [deps]
 DecFP = "55939f99-70c6-5e9b-8bb0-5071ed7d61fd"

docs/src/releasenotes.md

@@ -1,6 +1,19 @@
 
 # Release Notes
 
+## 0.4.1
+
+### New Features
+
+* Atomic solutions can now store a NamedTuple of internal variables of the integrator, including nonlinear solver output
+* Output of internal variables has been added to VPRK integrators
+* Add Gauss-Legendre tableaus for implicit partitioned Runge-Kutta methods
+
+### Fixes
+
+* Revision of integrator type hierarchy
+
+
 ## 0.4.0
 
 ### New Integrators

src/CommonFunctions.jl

@@ -33,4 +33,7 @@ module CommonFunctions
     eachsample() = nothing
     eachtimestep() = nothing
 
+    export nconstraints
+    nconstraints() = nothing
+
 end

src/Integrators.jl

@@ -33,8 +33,15 @@ module Integrators
     include("integrators/abstract_tableau.jl")
 
 
-    export Integrator, DeterministicIntegrator, StochasticIntegrator, IntegratorCache,
-           IntegratorConstructor
+    export Integrator, DeterministicIntegrator, StochasticIntegrator
+    export ODEIntegrator, DAEIntegrator, SDEIntegrator,
+           PODEIntegrator, PDAEIntegrator, PSDEIntegrator,
+           IODEIntegrator, IDAEIntegrator,
+           HODEIntegrator, HDAEIntegrator,
+           VODEIntegrator, VDAEIntegrator,
+           SPSDEIntegrator
+
+    export IntegratorCache, IntegratorConstructor
     export integrate, integrate!, integrate_step!, equation, timestep
     export NonlinearFunctionParameters, function_stages!
 

src/Solvers.jl

@@ -19,10 +19,11 @@ module Solvers
     export computeJacobian, check_jacobian, print_jacobian
 
     export NonlinearSolver, AbstractNewtonSolver,
-           NLsolveNewton,
-           NewtonSolver, QuasiNewtonSolver,
+           NLsolveNewton, NewtonSolver, QuasiNewtonSolver,
+           params, status,
            residual_initial!, residual_absolute!, residual_relative!,
            print_solver_status, check_solver_converged, check_solver_status,
+           get_solver_status, get_solver_status!,
            solve!
 
     include("solvers/nonlinear/nonlinear_solvers.jl")

src/Tableaus.jl

@@ -30,7 +30,8 @@ module Tableaus
            getCoefficientsLobIIIG2, getCoefficientsLobIIIG3, getCoefficientsLobIIIG4
 
     export getCoefficientsLobIII,  getCoefficientsLobIIIA, getCoefficientsLobIIIB,
-           getCoefficientsLobIIIC, getCoefficientsLobIIID, getCoefficientsLobIIIE
+           getCoefficientsLobIIIC, getCoefficientsLobIIID, getCoefficientsLobIIIE,
+           getCoefficientsLobIIIF, getCoefficientsLobIIIG
 
     include("tableaus/coefficients_lob.jl")
 
@@ -100,6 +101,10 @@ module Tableaus
     export getTableauSPARKGLRK, getTableauSPARKLobIIIAIIIB,
            getTableauSPARKGLRKLobIIIAIIIB
 
+    include("tableaus/tableaus_iprk.jl")
+
+    export getTableauIPGLRK           
+
     include("tableaus/tableaus_spark.jl")
 
     export getTableauVPGLRK,
@@ -124,15 +129,25 @@ module Tableaus
 
     include("tableaus/tableaus_vpark.jl")
 
-    export getTableauVSPARKMidpointProjection,
+    export getTableauVSPARKLobIIIAIIIBProjection,
+           getTableauVSPARKMidpointProjection,
            getTableauVSPARKSymmetricProjection,
+           getTableauVSPARKSymmetricLobProjection,
            getTableauVSPARKSymplecticProjection,
+           getTableauVSPARKGLRKpLobIIIAIIIB,
            getTableauVSPARKGLRKpMidpoint,
            getTableauVSPARKGLRKpSymmetric,
+           getTableauVSPARKGLRKpSymmetricLob,
            getTableauVSPARKGLRKpSymplectic,
+           getTableauVSPARKLobIIIAIIIB2pLobIIIAIIIB,
+           getTableauVSPARKLobIIIAIIIB3pLobIIIAIIIB,
+           getTableauVSPARKLobIIIAIIIB4pLobIIIAIIIB,
            getTableauVSPARKLobIIIAIIIB2pSymmetric,
            getTableauVSPARKLobIIIAIIIB3pSymmetric,
-           getTableauVSPARKLobIIIAIIIB4pSymmetric
+           getTableauVSPARKLobIIIAIIIB4pSymmetric,
+           getTableauVSPARKLobIIIAIIIB2pSymmetricLob,
+           getTableauVSPARKLobIIIAIIIB3pSymmetricLob,
+           getTableauVSPARKLobIIIAIIIB4pSymmetricLob
 
     include("tableaus/tableaus_vspark_primary.jl")
 

src/equations/dae.jl

@@ -132,8 +132,8 @@ function Base.similar(dae::DAE, t₀::TT, q₀::AbstractArray{DT}, λ₀::Abstra
     DAE(dae.v, dae.u, dae.ϕ, t₀, q₀, λ₀; h=h, parameters=parameters, periodicity=periodicity)
 end
 
-@inline Base.ndims(dae::DAE) = dae.d
-
+@inline Base.ndims(equation::DAE) = equation.d
+@inline CommonFunctions.nconstraints(equation::DAE) = equation.m
 @inline CommonFunctions.periodicity(equation::DAE) = equation.periodicity
 
 function get_function_tuple(equation::DAE{DT,TT,VT,UT,ϕT,HT,Nothing}) where {DT, TT, VT, UT, ϕT, HT}

src/equations/hdae.jl

@@ -123,8 +123,8 @@ function Base.similar(dae::HDAE, t₀::TT, q₀::AbstractArray{DT}, p₀::Abstra
          parameters=parameters, periodicity=periodicity)
 end
 
-@inline Base.ndims(dae::HDAE) = dae.d
-
+@inline Base.ndims(equation::HDAE) = equation.d
+@inline CommonFunctions.nconstraints(equation::HDAE) = equation.m
 @inline CommonFunctions.periodicity(equation::HDAE) = equation.periodicity
 
 function get_function_tuple(equation::HDAE{DT,TT,VT,FT,UT,GT,U̅T,G̅T,ϕT,ψT,HT,Nothing}) where {DT, TT, VT, FT, UT, GT, U̅T, G̅T, ϕT, ψT, HT}

src/equations/idae.jl

@@ -122,8 +122,8 @@ function Base.similar(dae::IDAE, t₀::TT, q₀::AbstractArray{DT}, p₀::Abstra
     IDAE(dae.ϑ, dae.f, dae.u, dae.g, dae.ϕ, t₀, q₀, p₀, λ₀; h=h, v=v, parameters=parameters, periodicity=periodicity)
 end
 
-@inline Base.ndims(dae::IDAE) = dae.d
-
+@inline Base.ndims(equation::IDAE) = equation.d
+@inline CommonFunctions.nconstraints(equation::IDAE) = equation.m
 @inline CommonFunctions.periodicity(equation::IDAE) = equation.periodicity
 
 function get_function_tuple(equation::IDAE{DT,TT,ϑT,FT,UT,GT,ϕT,HT,VT,Nothing}) where {DT, TT, ϑT, FT, UT, GT, ϕT, HT, VT}

src/equations/pdae.jl

@@ -117,8 +117,8 @@ function Base.similar(dae::PDAE, t₀::TT, q₀::AbstractArray{DT}, p₀::Abstra
     PDAE(dae.v, dae.f, dae.u, dae.g, dae.ϕ, t₀, q₀, p₀, λ₀; h=h, parameters=parameters, periodicity=periodicity)
 end
 
-@inline Base.ndims(dae::PDAE) = dae.d
-
+@inline Base.ndims(equation::PDAE) = equation.d
+@inline CommonFunctions.nconstraints(equation::PDAE) = equation.m
 @inline CommonFunctions.periodicity(equation::PDAE) = equation.periodicity
 
 function get_function_tuple(equation::PDAE{DT,TT,VT,FT,UT,GT,ϕT,HT,Nothing}) where {DT, TT, VT, FT, UT, GT, ϕT, HT}

src/equations/spdae.jl

@@ -105,8 +105,8 @@ function Base.similar(dae::SPDAE, t₀::TT, q₀::AbstractArray{DT}, p₀::Abstr
          parameters=parameters, periodicity=periodicity)
 end
 
-@inline Base.ndims(dae::SPDAE) = dae.d
-
+@inline Base.ndims(equation::SPDAE) = equation.d
+@inline CommonFunctions.nconstraints(equation::SPDAE) = equation.m
 @inline CommonFunctions.periodicity(equation::SPDAE) = equation.periodicity
 
 function get_function_tuple(equation::SPDAE{DT,TT,VT,FT,ϕT,ψT,Nothing}) where {DT, TT, VT, FT, ϕT, ψT}

src/equations/vdae.jl

@@ -175,8 +175,8 @@ function Base.similar(dae::VDAE, t₀::TT, q₀::AbstractArray{DT}, p₀::Abstra
          h=h, v=v, Ω=Ω, ∇H=∇H, parameters=parameters, periodicity=periodicity)
 end
 
-@inline Base.ndims(dae::VDAE) = dae.d
-
+@inline Base.ndims(equation::VDAE) = equation.d
+@inline CommonFunctions.nconstraints(equation::VDAE) = equation.m
 @inline CommonFunctions.periodicity(equation::VDAE) = equation.periodicity
 
 

src/integrators/SPARK.jl

@@ -13,9 +13,8 @@ module SPARK
 
     import ..Integrators
 
-    import ..Integrators: DeterministicIntegrator, IntegratorPRK, Parameters,
-                          IDAEIntegratorCache, InitialGuessIODE, InitialGuessPODE
-    import ..Integrators: IntegratorCache, CacheDict, CacheType
+    import ..Integrators: PDAEIntegrator, InitialGuessIODE, InitialGuessPODE, Parameters
+    import ..Integrators: IDAEIntegratorCache, IntegratorCache, CacheDict, CacheType
     import ..Integrators: AbstractTableau, AbstractTableauERK, AbstractTableauIRK,
                           AbstractCoefficients, CoefficientsRK,
                           @CoefficientsRK, @HeaderCoefficientsRK

src/integrators/Stochastic.jl

@@ -9,12 +9,13 @@ module Stochastic
     using ..Utils
 
     import ..Integrators
+    import ..Integrators: nstages, noisedims
 
     import ..Equations: SDE, PSDE, SPSDE, get_function_tuple
     import ..Solutions: AtomicSolutionSDE, AtomicSolutionPSDE, SolutionVector
     import ..Solutions: update!
 
-    import ..Integrators: StochasticIntegrator, Parameters
+    import ..Integrators: StochasticIntegrator, SDEIntegrator, PSDEIntegrator, Parameters
     import ..Integrators: SDEIntegratorCache, PSDEIntegratorCache,
                           IntegratorCache, CacheDict, CacheType
     import ..Integrators: CoefficientsRK, AbstractTableauERK, AbstractTableauIRK

src/integrators/VPRK.jl

@@ -14,9 +14,8 @@ module VPRK
 
     import ..Integrators
 
-    import ..Integrators: DeterministicIntegrator, IntegratorPRK, Parameters,
-                          IODEIntegratorCache, InitialGuessIODE
-    import ..Integrators: IntegratorCache, CacheDict, CacheType
+    import ..Integrators: IODEIntegrator, IODEIntegratorCache, InitialGuessIODE, IntegratorPRK
+    import ..Integrators: IntegratorCache, CacheDict, CacheType, Parameters
     import ..Integrators: AbstractTableauPRK, AbstractCoefficients,
                           CoefficientsRK, CoefficientsPGLRK,
                           @CoefficientsRK, @HeaderTableau, @HeaderCoefficientsRK,

src/integrators/abstract_integrator.jl

@@ -6,11 +6,11 @@ abstract type StochasticIntegrator{dType, tType} <: Integrator{dType, tType} end
 
 abstract type ODEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end
 abstract type DAEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end
-abstract type IODEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end
-abstract type IDAEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end
 abstract type PODEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end
 abstract type PDAEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end
 
+abstract type IODEIntegrator{dType, tType} <: PODEIntegrator{dType, tType} end
+abstract type IDAEIntegrator{dType, tType} <: PDAEIntegrator{dType, tType} end
 abstract type HODEIntegrator{dType, tType} <: PODEIntegrator{dType, tType} end
 abstract type HDAEIntegrator{dType, tType} <: PDAEIntegrator{dType, tType} end
 abstract type VODEIntegrator{dType, tType} <: IODEIntegrator{dType, tType} end
@@ -23,14 +23,41 @@ abstract type SPSDEIntegrator{dType, tType} <: StochasticIntegrator{dType, tType
 equation(integrator::Integrator) = error("equation() not implemented for ", typeof(integrator))
 timestep(integrator::Integrator) = error("timestep() not implemented for ", typeof(integrator))
 Base.ndims(integrator::Integrator) = error("ndims() not implemented for ", typeof(integrator))
+CommonFunctions.nconstraints(integrator::Integrator) = error("nconstraints() not implemented for ", typeof(integrator))
+noisedims(integrator::Integrator) = error("noisedims() not implemented for ", typeof(integrator))
+nstages(integrator::Integrator) = error("nstages() not implemented for ", typeof(integrator))
 
 eachdim(integrator::Integrator) = 1:ndims(integrator)
 
+get_internal_variables(::Integrator) = NamedTuple()
+
+
+Solutions.AtomicSolution(integrator::ODEIntegrator{DT,TT}) where {DT,TT} =
+    AtomicSolutionODE(DT, TT, ndims(integrator), get_internal_variables(integrator))
+
+Solutions.AtomicSolution(integrator::PODEIntegrator{DT,TT}) where {DT,TT} =
+    AtomicSolutionPODE(DT, TT, ndims(integrator), get_internal_variables(integrator))
+
+Solutions.AtomicSolution(integrator::DAEIntegrator{DT,TT}) where {DT,TT} =
+    AtomicSolutionDAE(DT, TT, ndims(integrator), nconstraints(integrator), get_internal_variables(integrator))
+
+Solutions.AtomicSolution(integrator::PDAEIntegrator{DT,TT}) where {DT,TT} =
+    AtomicSolutionPDAE(DT, TT, ndims(integrator), nconstraints(integrator), get_internal_variables(integrator))
+
+Solutions.AtomicSolution(integrator::SDEIntegrator{DT,TT}) where {DT,TT} =
+    AtomicSolutionSDE(DT, TT, ndims(integrator), noisedims(integrator), get_internal_variables(integrator))
+
+Solutions.AtomicSolution(integrator::PSDEIntegrator{DT,TT}) where {DT,TT} =
+    AtomicSolutionPSDE(DT, TT, ndims(integrator), noisedims(integrator), get_internal_variables(integrator))
+
+Solutions.AtomicSolution(integrator::SPSDEIntegrator{DT,TT}) where {DT,TT} =
+    AtomicSolutionPSDE(DT, TT, ndims(integrator), noisedims(integrator), get_internal_variables(integrator))
+
 
 abstract type Parameters{DT,TT} end
 
-function_stages!(x::Vector{DT}, b::Vector{DT}, params::PT) where {DT, TT, PT <: Parameters{DT,TT}} = error("function_stages!() not implemented for ", PT)
-solution_stages!(x::Vector{DT}, y::Vector{DT}, params::PT) where {DT, TT, PT <: Parameters{DT,TT}} = error("solution_stages!() not implemented for ", PT)
+function_stages!(::Vector{DT}, ::Vector{DT}, ::PT) where {DT, TT, PT <: Parameters{DT,TT}} = error("function_stages!() not implemented for ", PT)
+solution_stages!(::Vector{DT}, ::Vector{DT}, ::PT) where {DT, TT, PT <: Parameters{DT,TT}} = error("solution_stages!() not implemented for ", PT)
 
 initialize!(::Integrator, ::AtomicSolution) = nothing
 

src/integrators/cgvi/integrators_cgvi.jl

@@ -127,7 +127,7 @@ struct IntegratorCGVI{DT, TT, D, S, R,
                       BT <: Basis,
                       PT <: ParametersCGVI{DT,TT,D,S,R},
                       ST <: NonlinearSolver{DT},
-                      IT <: InitialGuessIODE{DT,TT}} <: DeterministicIntegrator{DT,TT}
+                      IT <: InitialGuessIODE{DT,TT}} <: IODEIntegrator{DT,TT}
     basis::BT
     quadrature::Quadrature{TT,R}
 
@@ -171,6 +171,8 @@ end
 @inline equation(integrator::IntegratorCGVI, i::Symbol) = integrator.params.equs[i]
 @inline equations(integrator::IntegratorCGVI) = integrator.params.equs
 @inline timestep(integrator::IntegratorCGVI) = integrator.params.Δt
+@inline Base.ndims(::IntegratorCGVI{DT,TT,D}) where {DT,TT,D} = D
+
 
 function IntegratorCache(int::IntegratorCGVI{DT,TT}) where {DT,TT}
     IntegratorCacheCGVI{DT, TT, ndims(int), nbasis(int.basis), nnodes(int.quadrature)}()

src/integrators/dgvi/integrators_dgvi.jl

@@ -380,7 +380,7 @@ struct IntegratorDGVI{DT, TT, D, S, R,
                 BT <: Basis,
                 PT <: ParametersDGVI{DT,TT,D,S},
                 ST <: NonlinearSolver{DT},
-                IT <: InitialGuessODE{DT,TT}} <: DeterministicIntegrator{DT,TT}
+                IT <: InitialGuessODE{DT,TT}} <: IODEIntegrator{DT,TT}
     basis::BT
     quadrature::Quadrature{TT,R}
 
@@ -424,6 +424,8 @@ end
 @inline equation(integrator::IntegratorDGVI, i::Symbol) = integrator.params.equs[i]
 @inline equations(integrator::IntegratorDGVI) = integrator.params.equs
 @inline timestep(integrator::IntegratorDGVI) = integrator.params.Δt
+@inline Base.ndims(::IntegratorDGVI{DT,TT,D}) where {DT,TT,D} = D
+
 
 function IntegratorCache(int::IntegratorDGVI{DT,TT}) where {DT,TT}
     IntegratorCacheDGVI{DT, TT, ndims(int), nbasis(int.basis), nnodes(int.quadrature)}()

src/integrators/dgvi/integrators_dgvi_experimental.jl

@@ -124,7 +124,7 @@ struct IntegratorDGVIEXP{DT, TT, D, S, R,
                 BT <: Basis,
                 PT <: ParametersDGVIEXP{DT,TT,D,S},
                 ST <: NonlinearSolver{DT},
-                IT <: InitialGuessODE{DT,TT}} <: DeterministicIntegrator{DT,TT}
+                IT <: InitialGuessODE{DT,TT}} <: IODEIntegrator{DT,TT}
     basis::BT
     quadrature::Quadrature{TT,R}
 
@@ -177,6 +177,7 @@ end
 @inline equation(integrator::IntegratorDGVIEXP, i::Symbol) = integrator.params.equs[i]
 @inline equations(integrator::IntegratorDGVIEXP) = integrator.params.equs
 @inline timestep(integrator::IntegratorDGVIEXP) = integrator.params.Δt
+@inline Base.ndims(::IntegratorDGVIEXP{DT,TT,D}) where {DT,TT,D} = D
 
 
 function update_params!(params::ParametersDGVIEXP, int::IntegratorDGVIEXP)

src/integrators/dgvi/integrators_dgvi_path_integral.jl

@@ -266,7 +266,7 @@ struct IntegratorDGVIPI{DT, TT, D, S, R,
                 JT <: Discontinuity,
                 PT <: ParametersDGVIPI{DT,TT,D,S},
                 ST <: NonlinearSolver{DT},
-                IT <: InitialGuessODE{DT,TT}} <: DeterministicIntegrator{DT,TT}
+                IT <: InitialGuessODE{DT,TT}} <: IODEIntegrator{DT,TT}
     basis::BT
     quadrature::Quadrature{TT,R}
     jump::JT
@@ -321,6 +321,7 @@ end
 @inline equation(integrator::IntegratorDGVIPI, i::Symbol) = integrator.params.equs[i]
 @inline equations(integrator::IntegratorDGVIPI) = integrator.params.equs
 @inline timestep(integrator::IntegratorDGVIPI) = integrator.params.Δt
+@inline Base.ndims(::IntegratorDGVIPI{DT,TT,D}) where {DT,TT,D} = D
 
 
 function update_params!(params::ParametersDGVIPI, int::IntegratorDGVIPI)

src/integrators/dgvi/integrators_dgvi_projection_final.jl

@@ -124,7 +124,7 @@ struct IntegratorDGVIP1{DT, TT, D, S, R,
                 BT <: Basis,
                 PT <: ParametersDGVIP1{DT,TT,D,S},
                 ST <: NonlinearSolver{DT},
-                IT <: InitialGuessODE{DT,TT}} <: DeterministicIntegrator{DT,TT}
+                IT <: InitialGuessODE{DT,TT}} <: IODEIntegrator{DT,TT}
     basis::BT
     quadrature::Quadrature{TT,R}
 
@@ -177,6 +177,7 @@ end
 @inline equation(integrator::IntegratorDGVIP1, i::Symbol) = integrator.params.equs[i]
 @inline equations(integrator::IntegratorDGVIP1) = integrator.params.equs
 @inline timestep(integrator::IntegratorDGVIP1) = integrator.params.Δt
+@inline Base.ndims(::IntegratorDGVIP1{DT,TT,D}) where {DT,TT,D} = D
 
 
 function update_params!(params::ParametersDGVIP1, int::IntegratorDGVIP1)

src/integrators/dgvi/integrators_dgvi_projection_initial.jl

@@ -127,7 +127,7 @@ struct IntegratorDGVIP0{DT, TT, D, S, R,
                 BT <: Basis,
                 PT <: ParametersDGVIP0{DT,TT,D,S},
                 ST <: NonlinearSolver{DT},
-                IT <: InitialGuessODE{DT,TT}} <: DeterministicIntegrator{DT,TT}
+                IT <: InitialGuessODE{DT,TT}} <: IODEIntegrator{DT,TT}
     basis::BT
     quadrature::Quadrature{TT,R}
 
@@ -188,6 +188,7 @@ end
 @inline equation(integrator::IntegratorDGVIP0, i::Symbol) = integrator.params.equs[i]
 @inline equations(integrator::IntegratorDGVIP0) = integrator.params.equs
 @inline timestep(integrator::IntegratorDGVIP0) = integrator.params.Δt
+@inline Base.ndims(::IntegratorDGVIP0{DT,TT,D}) where {DT,TT,D} = D
 
 
 function update_params!(params::ParametersDGVIP0, int::IntegratorDGVIP0)