Update for GeometricBase v0.10 and GeometricEquations v0.16.

Project.toml

@@ -33,8 +33,8 @@ DecFP = "1"
 Documenter = "0.25, 0.26, 0.27, 1"
 ForwardDiff = "0.10"
 GenericLinearAlgebra = "0.2, 0.3"
-GeometricBase = "0.9, 0.10"
-GeometricEquations = "0.14, 0.15"
+GeometricBase = "0.10"
+GeometricEquations = "0.16"
 GeometricSolutions = "0.3"
 HDF5 = "0.15, 0.16, 0.17"
 NLsolve = "4"

src/SPARK.jl

@@ -21,7 +21,7 @@ module SPARK
 
     import ..Integrators: GeometricIntegrator, Newton
     import ..Integrators: HDAEMethod, IDAEMethod, LDAEMethod, PDAEMethod
-    import ..Integrators: SolutionStepPDAE, SolutionVector
+    import ..Integrators: SolutionStepPDAE
     import ..Integrators: InitialGuess, Extrapolation, HermiteExtrapolation
     import ..Integrators: initialguess!, initial_guess!, integrate_step!, residual!
     import ..Integrators: CacheDict, Cache, CacheType, IDAEIntegratorCache

src/initial_guess/initial_guess.jl

@@ -20,10 +20,10 @@ function initialguess!(t, q, p, q̇, ṗ, solstep::SolutionStepPODE, ::Union{Abs
     (t = t, q = q, p = p, v = q̇, f = ṗ)
 end
 
-function initialguess!(t₀, q₀, p₀, λ₀, q̇₀, ṗ₀, u₀, g₀, t₁, q₁, p₁, λ₁, q̇₁, ṗ₁, u₁, g₁, t, q, p, iguess::Union{Extrapolation,InitialGuess}; kwargs...)
+function initialguess!(t₀, q₀, p₀, λ₀, μ₀, q̇₀, ṗ₀, u₀, g₀, t₁, q₁, p₁, λ₁, μ₁, q̇₁, ṗ₁, u₁, g₁, t, q, p, iguess::Union{Extrapolation,InitialGuess}; kwargs...)
     initialguess!(t₀, q₀, p₀, q̇₀, ṗ₀, t₁, q₁, p₁, q̇₁, ṗ₁, t, q, p, iguess; kwargs...)
 end
 
-function initialguess!(t₀, q₀, p₀, λ₀, q̇₀, ṗ₀, u₀, g₀, t₁, q₁, p₁, λ₁, q̇₁, ṗ₁, u₁, g₁, t, q, p, q̇, ṗ, iguess::Union{Extrapolation,InitialGuess}; kwargs...)
+function initialguess!(t₀, q₀, p₀, λ₀, μ₀, q̇₀, ṗ₀, u₀, g₀, t₁, q₁, p₁, λ₁, μ₁, q̇₁, ṗ₁, u₁, g₁, t, q, p, q̇, ṗ, iguess::Union{Extrapolation,InitialGuess}; kwargs...)
     initialguess!(t₀, q₀, p₀, q̇₀, ṗ₀, t₁, q₁, p₁, q̇₁, ṗ₁, t, q, p, q̇, ṗ, iguess; kwargs...)
 end

src/integrators/VPRK.jl

@@ -12,7 +12,7 @@ module VPRK
 
     using ..Methods
 
-    import ..Solutions: SolutionStepPODE, SolutionStepPDAE, SolutionPDAE, SolutionVector
+    import ..Solutions: SolutionStepPODE, SolutionStepPDAE, SolutionPDAE
 
     import ..Integrators
 

src/integrators/rk/integrators_irk.jl

@@ -114,7 +114,7 @@ end
 
 nlsolution(cache::IRKCache) = cache.x
 
-reset!(cache::IRKCache, t, q, λ = missing) = copyto!(cache.q̄, q)
+reset!(cache::IRKCache, t, q, λ = missing, μ = missing) = copyto!(cache.q̄, q)
 
 function Cache{ST}(problem::AbstractProblem, method::IRKMethod; kwargs...) where {ST}
     S = nstages(tableau(method))

src/integrators/rk/updates.jl

@@ -64,7 +64,7 @@ function update_vector!(Δq::AbstractVector, Δp::AbstractVector, V::StageVector
 end
 
 
-function update_multiplier!(λ::SolutionVector{T}, Λ::StageVector{T}, b::AbstractVector) where {T}
+function update_multiplier!(λ::AbstractVector{T}, Λ::StageVector{T}, b::AbstractVector) where {T}
     for Λᵢ in Λ
         @assert length(λ) == length(Λᵢ)
     end
@@ -104,7 +104,7 @@ end
 #     end
 # end
 
-# function update_solution!(x::SolutionVector{T}, ẋ::Matrix{T}, b::AbstractVector, Δt) where {T}
+# function update_solution!(x::AbstractVector{T}, ẋ::Matrix{T}, b::AbstractVector, Δt) where {T}
 #     @assert length(x) == size(ẋ, 1)
 #     @assert length(b) == size(ẋ, 2)
 
@@ -120,7 +120,7 @@ end
 # end
 
 
-# function update_solution!(x::SolutionVector{T}, ẋ::Vector{Vector{T}}, b::AbstractVector, Δt) where {T}
+# function update_solution!(x::AbstractVector{T}, ẋ::Vector{Vector{T}}, b::AbstractVector, Δt) where {T}
 #     @assert length(b) == length(ẋ)
 #     @assert length(x) == length(ẋ[1])
 
@@ -140,12 +140,12 @@ end
 #     update_solution!(x, xₑᵣᵣ, ẋ, b̂, Δt)
 # end
 
-# function update_solution!(x::SolutionVector{T}, ẋ::Union{Matrix{T},Vector{Vector{T}}}, b::AbstractVector, b̂::AbstractVector, Δt) where {T}
+# function update_solution!(x::AbstractVector{T}, ẋ::Union{Matrix{T},Vector{Vector{T}}}, b::AbstractVector, b̂::AbstractVector, Δt) where {T}
 #     update_solution!(x, ẋ, b, Δt)
 #     update_solution!(x, ẋ, b̂, Δt)
 # end
 
-# function update_multiplier!(λ::SolutionVector{T}, Λ::Vector{Vector{T}}, b::AbstractVector) where {T}
+# function update_multiplier!(λ::AbstractVector{T}, Λ::Vector{Vector{T}}, b::AbstractVector) where {T}
 #     for Λᵢ in Λ
 #         @assert length(λ) == length(Λᵢ)
 #     end

src/integrators/vprk/integrators_vprk_pinternal.jl

@@ -97,7 +97,7 @@ end
 
 
 function compute_projection_vprk!(x::Vector{ST},
-                q::SolutionVector{ST}, p::SolutionVector{ST}, v::SolutionVector{ST}, λ::SolutionVector{ST},
+                q::AbstractVector{ST}, p::AbstractVector{ST}, v::AbstractVector{ST}, λ::AbstractVector{ST},
                 Q::Vector{Vector{ST}}, V::Vector{Vector{ST}}, U::Vector{Vector{ST}}, G::Vector{Vector{ST}},
                 params::ParametersVPRKpInternal{DT,TT,D,S}) where {ST,DT,TT,D,S}
 

src/integrators/vprk/integrators_vprk_pvariational.jl

@@ -118,7 +118,7 @@ end
 
 
 function compute_projection!(
-                x::Vector{ST}, q::SolutionVector{ST}, p::SolutionVector{ST}, v::SolutionVector{ST}, λ::SolutionVector{ST},
+                x::Vector{ST}, q::AbstractVector{ST}, p::AbstractVector{ST}, v::AbstractVector{ST}, λ::AbstractVector{ST},
                 U::Vector{Vector{ST}}, G::Vector{Vector{ST}},
                 params::ParametersVPRKpVariational{DT,TT,D,S}) where {ST,DT,TT,D,S}
 

src/solutions/solution_step_constructors.jl

@@ -38,8 +38,7 @@ end
 
 # Create SolutionStep for a PDAEProblem.
 function SolutionStep(problem::Union{PDAEProblem, HDAEProblem, IDAEProblem, LDAEProblem}, extrap::Extrapolation = default_extrapolation(); kwargs...)
-    ics = initial_conditions(problem)
-    solstep = SolutionStepPDAE(ics.t, ics.q, ics.p, ics.λ, parameters(problem); kwargs...)
+    solstep = SolutionStepPDAE(initial_conditions(problem)..., parameters(problem); kwargs...)
     initialize!(solstep, problem, extrap)
     return solstep
 end

src/solutions/solution_step_dae.jl

@@ -14,8 +14,10 @@ Solution step for a [`DAEProblem`](@ref).
 * `t̄`: time of previous time steps
 * `q`: current solution of q
 * `λ`: current solution of λ
+* `μ`: current solution of μ
 * `q̄`: previous solutions of q
 * `λ̄`: previous solutions of λ
+* `μ̄`: previous solution of μ
 * `v`: vector field of q
 * `v̄`: vector field of q̄
 * `u`: projective vector field of q
@@ -43,11 +45,13 @@ struct SolutionStepDAE{
 
     q::AT
     λ::ΛT
+    μ::ΛT
     v::VT
     u::VT
 
     q̄::AT
     λ̄::ΛT
+    μ̄::ΛT
     v̄::VT
     ū::VT
 
@@ -58,31 +62,34 @@ struct SolutionStepDAE{
 
     parameters::paramsType
 
-    function SolutionStepDAE(t::TT, q::AT, λ::ΛT, parameters; nhistory = 2, internal::IT = NamedTuple()) where {DT, TT, AT <: AbstractArray{DT}, ΛT <: AbstractArray{DT}, IT}
+    function SolutionStepDAE(t::TT, q::AT, λ::ΛT, μ::ΛT, parameters; nhistory = 2, internal::IT = NamedTuple()) where {DT, TT, AT <: AbstractArray{DT}, ΛT <: AbstractArray{DT}, IT}
         # TODO: nhistory should default to 1 and set to higher values by integrator / initial guess method
         @assert nhistory ≥ 1
 
         history = (
             t = OffsetVector([zero(t) for _ in 0:nhistory], 0:nhistory),
             q = OffsetVector([zero(q) for _ in 0:nhistory], 0:nhistory),
             λ = OffsetVector([zero(λ) for _ in 0:nhistory], 0:nhistory),
+            μ = OffsetVector([zero(μ) for _ in 0:nhistory], 0:nhistory),
             v = OffsetVector([vectorfield(q) for _ in 0:nhistory], 0:nhistory),
             u = OffsetVector([vectorfield(q) for _ in 0:nhistory], 0:nhistory),
         )
 
         q = history.q[0]
         λ = history.λ[0]
+        μ = history.μ[0]
         v = history.v[0]
         u = history.u[0]
 
         q̄ = history.q[1]
         λ̄ = history.λ[1]
+        μ̄ = history.μ[1]
         v̄ = history.v[1]
         ū = history.u[1]
 
         q̃ = zero(q)
 
-        new{DT, TT, AT, ΛT, typeof(v), typeof(history), IT, typeof(parameters), nhistory}(q, λ, v, u, q̄, λ̄, v̄, ū, q̃, history, internal, parameters)
+        new{DT, TT, AT, ΛT, typeof(v), typeof(history), IT, typeof(parameters), nhistory}(q, λ, μ, v, u, q̄, λ̄, μ̄, v̄, ū, q̃, history, internal, parameters)
     end
 end
 
@@ -112,13 +119,14 @@ end
 
 nhistory(::SolutionStepDAE{DT,TT,AT,ΛT,VT,HT,IT,PT,NT}) where {DT,TT,AT,ΛT,VT,HT,IT,PT,NT} = NT
 
-current(solstep::SolutionStepDAE) = (t = solstep.t, q = solstep.q, λ = solstep.λ)
-previous(solstep::SolutionStepDAE) = (t = solstep.t̄, q = solstep.q̄, λ = solstep.λ̄)
+current(solstep::SolutionStepDAE) = (t = solstep.t, q = solstep.q, λ = solstep.λ, μ = solstep.μ)
+previous(solstep::SolutionStepDAE) = (t = solstep.t̄, q = solstep.q̄, λ = solstep.λ̄, μ = solstep.μ̄)
 history(solstep::SolutionStepDAE) = solstep.history
 history(solstep::SolutionStepDAE, i::Int) = (
     t = history(solstep).t[i],
     q = history(solstep).q[i],
     λ = history(solstep).λ[i],
+    μ = history(solstep).μ[i],
     v = history(solstep).v[i],
     u = history(solstep).u[i])
 internal(solstep::SolutionStepDAE) = solstep.internal
@@ -134,6 +142,7 @@ function initialize!(solstep::SolutionStepDAE, problem::DAEProblem, extrap::Extr
     solstep.t  = initial_conditions(problem).t
     solstep.q .= initial_conditions(problem).q
     solstep.λ .= initial_conditions(problem).λ
+    solstep.μ .= initial_conditions(problem).μ
     solstep.q̃ .= 0
 
     update_vector_fields!(solstep, problem)
@@ -147,37 +156,50 @@ function initialize!(solstep::SolutionStepDAE, problem::DAEProblem, extrap::Extr
     return solstep
 end
 
-function update!(solstep::SolutionStepDAE{DT,TT,AT,ΛT}, Δq::AT) where {DT,TT,AT,ΛT}
+function update!(solstep::SolutionStepDAE{DT,TT,AT,ΛT}, Δq::Union{AT,AbstractArray}) where {DT,TT,AT,ΛT}
     for k in eachindex(Δq)
         solstep.q[k], solstep.q̃[k] = compensated_summation(Δq[k], solstep.q[k], solstep.q̃[k])
     end
     return solstep
 end
 
-function update!(solstep::SolutionStepDAE{DT,TT,AT,ΛT}, q̇::AT, Δt::TT) where {DT,TT,AT,ΛT}
+function update!(solstep::SolutionStepDAE{DT,TT,AT,ΛT}, q̇::Union{AT,AbstractArray}, Δt::TT) where {DT,TT,AT,ΛT}
     for k in eachindex(q̇)
         solstep.q[k], solstep.q̃[k] = compensated_summation(Δt * q̇[k], solstep.q[k], solstep.q̃[k])
     end
     return solstep
 end
 
-function update!(solstep::SolutionStepDAE{DT,TT,AT,ΛT}, Δq::AT, λ::ΛT) where {DT,TT,AT,ΛT}
+function update!(solstep::SolutionStepDAE{DT,TT,AT,ΛT}, Δq::Union{AT,AbstractArray}, λ::Union{ΛT,AbstractArray}) where {DT,TT,AT,ΛT}
     update!(solstep, Δq)
     copyto!(solstep.λ, λ)
     return solstep
 end
 
-function update!(solstep::SolutionStepDAE{DT,TT,AT,ΛT}, q̇::AT, λ::ΛT, Δt::TT) where {DT,TT,AT,ΛT}
+function update!(solstep::SolutionStepDAE{DT,TT,AT,ΛT}, Δq::Union{AT,AbstractArray}, λ::Union{ΛT,AbstractArray}, μ::Union{ΛT,AbstractArray}) where {DT,TT,AT,ΛT}
+    update!(solstep, Δq, λ)
+    copyto!(solstep.μ, μ)
+    return solstep
+end
+
+function update!(solstep::SolutionStepDAE{DT,TT,AT,ΛT}, q̇::Union{AT,AbstractArray}, λ::Union{ΛT,AbstractArray}, Δt::TT) where {DT,TT,AT,ΛT}
     update!(solstep, q̇, Δt)
     copyto!(solstep.λ, λ)
     return solstep
 end
 
+function update!(solstep::SolutionStepDAE{DT,TT,AT,ΛT}, q̇::Union{AT,AbstractArray}, λ::Union{ΛT,AbstractArray}, μ::Union{ΛT,AbstractArray}, Δt::TT) where {DT,TT,AT,ΛT}
+    update!(solstep, q̇, λ, Δt)
+    copyto!(solstep.μ, μ)
+    return solstep
+end
+
 function GeometricBase.reset!(solstep::SolutionStepDAE, Δt)
     for i in backwardhistory(solstep)
         history(solstep).t[i]  = history(solstep).t[i-1]
         history(solstep).q[i] .= history(solstep).q[i-1]
         history(solstep).λ[i] .= history(solstep).λ[i-1]
+        history(solstep).μ[i] .= history(solstep).μ[i-1]
         history(solstep).v[i] .= history(solstep).v[i-1]
         history(solstep).u[i] .= history(solstep).u[i-1]
     end
@@ -192,13 +214,15 @@ function Base.copy!(solstep::SolutionStepDAE, sol::NamedTuple)
     solstep.t  = sol.t
     solstep.q .= sol.q
     solstep.λ .= sol.λ
+    solstep.μ .= sol.μ
     solstep.q̃ .= 0
     solstep.v .= 0
 
     # for i in backwardhistory(solstep)
     #     solstep.t̄[i]  = 0
     #     solstep.q̄[i] .= 0
     #     solstep.λ̄[i] .= 0
+    #     solstep.μ̄[i] .= 0
     #     solstep.v̄[i] .= 0
     #     solstep.ū[i] .= 0
     # end