From cd6b16400209eb03120348b7c8f7e20fbeffe9dd Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Fri, 1 Mar 2024 15:37:49 +0100 Subject: [PATCH] Update for GeometricBase v0.10 and GeometricEquations v0.16. --- Project.toml | 4 +- src/SPARK.jl | 2 +- src/initial_guess/initial_guess.jl | 4 +- src/integrators/VPRK.jl | 2 +- src/integrators/rk/integrators_irk.jl | 2 +- src/integrators/rk/updates.jl | 10 ++-- .../vprk/integrators_vprk_pinternal.jl | 2 +- .../vprk/integrators_vprk_pvariational.jl | 2 +- src/solutions/solution_step_constructors.jl | 3 +- src/solutions/solution_step_dae.jl | 40 ++++++++++--- src/solutions/solution_step_pdae.jl | 40 ++++++++++--- .../deterministic_solutions_tests.jl | 34 +++++------ test/solutions/solution_step_tests.jl | 56 ++++++++++++++----- 13 files changed, 137 insertions(+), 64 deletions(-) diff --git a/Project.toml b/Project.toml index 5548d174f..bb47309f1 100644 --- a/Project.toml +++ b/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" diff --git a/src/SPARK.jl b/src/SPARK.jl index 14986b854..cbaf6fe09 100644 --- a/src/SPARK.jl +++ b/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 diff --git a/src/initial_guess/initial_guess.jl b/src/initial_guess/initial_guess.jl index ea618271f..5e0c277e6 100644 --- a/src/initial_guess/initial_guess.jl +++ b/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 diff --git a/src/integrators/VPRK.jl b/src/integrators/VPRK.jl index 42838c50a..20c7b3571 100644 --- a/src/integrators/VPRK.jl +++ b/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 diff --git a/src/integrators/rk/integrators_irk.jl b/src/integrators/rk/integrators_irk.jl index d196642ad..0eb086c41 100644 --- a/src/integrators/rk/integrators_irk.jl +++ b/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)) diff --git a/src/integrators/rk/updates.jl b/src/integrators/rk/updates.jl index ec23c1d65..24ba5f6c4 100644 --- a/src/integrators/rk/updates.jl +++ b/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 diff --git a/src/integrators/vprk/integrators_vprk_pinternal.jl b/src/integrators/vprk/integrators_vprk_pinternal.jl index e401d45e6..1cd463f69 100644 --- a/src/integrators/vprk/integrators_vprk_pinternal.jl +++ b/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} diff --git a/src/integrators/vprk/integrators_vprk_pvariational.jl b/src/integrators/vprk/integrators_vprk_pvariational.jl index dc87be846..7b0a6e64b 100644 --- a/src/integrators/vprk/integrators_vprk_pvariational.jl +++ b/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} diff --git a/src/solutions/solution_step_constructors.jl b/src/solutions/solution_step_constructors.jl index f4250e9d9..9c8cdc918 100644 --- a/src/solutions/solution_step_constructors.jl +++ b/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 diff --git a/src/solutions/solution_step_dae.jl b/src/solutions/solution_step_dae.jl index 499ef3a61..5a55d7821 100644 --- a/src/solutions/solution_step_dae.jl +++ b/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,7 +62,7 @@ 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 @@ -66,23 +70,26 @@ struct SolutionStepDAE{ 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,6 +214,7 @@ 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 @@ -199,6 +222,7 @@ function Base.copy!(solstep::SolutionStepDAE, sol::NamedTuple) # solstep.t̄[i] = 0 # solstep.q̄[i] .= 0 # solstep.λ̄[i] .= 0 + # solstep.μ̄[i] .= 0 # solstep.v̄[i] .= 0 # solstep.ū[i] .= 0 # end diff --git a/src/solutions/solution_step_pdae.jl b/src/solutions/solution_step_pdae.jl index 60aa146a1..4fe302965 100644 --- a/src/solutions/solution_step_pdae.jl +++ b/src/solutions/solution_step_pdae.jl @@ -15,9 +15,11 @@ Solution step for a [`PDAEProblem`](@ref). * `q`: current solution of q * `p`: current solution of p * `λ`: current solution of λ +* `μ`: current solution of μ * `q̄`: previous solution of q * `p̄`: previous solution of p * `λ̄`: previous solution of λ +* `μ̄`: previous solution of μ * `v`: vector field of q * `f`: vector field of p * `v̄`: vector field of q̄ @@ -52,6 +54,7 @@ mutable struct SolutionStepPDAE{ q::AT p::AT λ::ΛT + μ::ΛT v::VT f::FT u::VT @@ -60,6 +63,7 @@ mutable struct SolutionStepPDAE{ q̄::AT p̄::AT λ̄::ΛT + μ̄::ΛT v̄::VT f̄::FT ū::VT @@ -73,7 +77,7 @@ mutable struct SolutionStepPDAE{ parameters::paramsType - function SolutionStepPDAE(t::TT, q::AT, p::AT, λ::ΛT, parameters; nhistory = 2, internal::IT = NamedTuple()) where {DT, TT, AT <: AbstractArray{DT}, ΛT <: AbstractArray{DT}, IT} + function SolutionStepPDAE(t::TT, q::AT, p::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 @@ -82,6 +86,7 @@ mutable struct SolutionStepPDAE{ q = OffsetVector([zero(q) for _ in 0:nhistory], 0:nhistory), p = OffsetVector([zero(p) 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), f = OffsetVector([vectorfield(p) for _ in 0:nhistory], 0:nhistory), u = OffsetVector([vectorfield(q) for _ in 0:nhistory], 0:nhistory), @@ -91,6 +96,7 @@ mutable struct SolutionStepPDAE{ q = history.q[0] p = history.p[0] λ = history.λ[0] + μ = history.μ[0] v = history.v[0] f = history.f[0] u = history.u[0] @@ -99,6 +105,7 @@ mutable struct SolutionStepPDAE{ q̄ = history.q[1] p̄ = history.p[1] λ̄ = history.λ[1] + μ̄ = history.μ[1] v̄ = history.v[1] f̄ = history.f[1] ū = history.u[1] @@ -107,7 +114,7 @@ mutable struct SolutionStepPDAE{ q̃ = zero(q) p̃ = zero(p) - new{DT, TT, AT, ΛT, typeof(v), typeof(f), typeof(history), IT, typeof(parameters), nhistory}(q, p, λ, v, f, u, g, q̄, p̄, λ̄, v̄, f̄, ū, ḡ, q̃, p̃, history, internal, parameters) + new{DT, TT, AT, ΛT, typeof(v), typeof(f), typeof(history), IT, typeof(parameters), nhistory}(q, p, λ, μ, v, f, u, g, q̄, p̄, λ̄, μ̄, v̄, f̄, ū, ḡ, q̃, p̃, history, internal, parameters) end end @@ -137,14 +144,15 @@ end nhistory(::SolutionStepPDAE{DT,TT,AT,ΛT,VT,FT,HT,IT,PT,NT}) where {DT,TT,AT,ΛT,VT,FT,HT,IT,PT,NT} = NT -current(solstep::SolutionStepPDAE) = (t = solstep.t, q = solstep.q, p = solstep.p, λ = solstep.λ) -previous(solstep::SolutionStepPDAE) = (t = solstep.t̄, q = solstep.q̄, p = solstep.p̄, λ = solstep.λ̄) +current(solstep::SolutionStepPDAE) = (t = solstep.t, q = solstep.q, p = solstep.p, λ = solstep.λ, μ = solstep.μ) +previous(solstep::SolutionStepPDAE) = (t = solstep.t̄, q = solstep.q̄, p = solstep.p̄, λ = solstep.λ̄, μ = solstep.μ̄) history(solstep::SolutionStepPDAE) = solstep.history history(solstep::SolutionStepPDAE, i::Int) = ( t = history(solstep).t[i], q = history(solstep).q[i], p = history(solstep).p[i], λ = history(solstep).λ[i], + μ = history(solstep).μ[i], v = history(solstep).v[i], f = history(solstep).f[i], u = history(solstep).u[i], @@ -178,6 +186,7 @@ function initialize!(solstep::SolutionStepPDAE, problem::Union{PDAEProblem, HDAE solstep.q .= initial_conditions(problem).q solstep.p .= initial_conditions(problem).p solstep.λ .= initial_conditions(problem).λ + solstep.μ .= initial_conditions(problem).μ solstep.q̃ .= 0 solstep.p̃ .= 0 @@ -195,7 +204,7 @@ function initialize!(solstep::SolutionStepPDAE, problem::Union{PDAEProblem, HDAE return solstep end -function update!(solstep::SolutionStepPDAE{DT,TT,AT,ΛT}, Δq::AT, Δp::AT) where {DT,TT,AT,ΛT} +function update!(solstep::SolutionStepPDAE{DT,TT,AT,ΛT}, Δq::Union{AT,AbstractArray}, Δp::Union{AT,AbstractArray}) where {DT,TT,AT,ΛT} for k in eachindex(Δq,Δp) solstep.q[k], solstep.q̃[k] = compensated_summation(Δq[k], solstep.q[k], solstep.q̃[k]) solstep.p[k], solstep.p̃[k] = compensated_summation(Δp[k], solstep.p[k], solstep.p̃[k]) @@ -203,7 +212,7 @@ function update!(solstep::SolutionStepPDAE{DT,TT,AT,ΛT}, Δq::AT, Δp::AT) wher return solstep end -function update!(solstep::SolutionStepPDAE{DT,TT,AT,ΛT}, q̇::AT, ṗ::AT, Δt::TT) where {DT,TT,AT,ΛT} +function update!(solstep::SolutionStepPDAE{DT,TT,AT,ΛT}, q̇::Union{AT,AbstractArray}, ṗ::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]) solstep.p[k], solstep.p̃[k] = compensated_summation(Δt * ṗ[k], solstep.p[k], solstep.p̃[k]) @@ -211,24 +220,37 @@ function update!(solstep::SolutionStepPDAE{DT,TT,AT,ΛT}, q̇::AT, ṗ::AT, Δt: return solstep end -function update!(solstep::SolutionStepPDAE{DT,TT,AT,ΛT}, Δq::AT, Δp::AT, λ::ΛT) where {DT,TT,AT,ΛT} +function update!(solstep::SolutionStepPDAE{DT,TT,AT,ΛT}, Δq::Union{AT,AbstractArray}, Δp::Union{AT,AbstractArray}, λ::Union{ΛT,AbstractArray}) where {DT,TT,AT,ΛT} update!(solstep, Δq, Δp) copyto!(solstep.λ, λ) return solstep end -function update!(solstep::SolutionStepPDAE{DT,TT,AT,ΛT}, q̇::AT, ṗ::AT, λ::ΛT, Δt::TT) where {DT,TT,AT,ΛT} +function update!(solstep::SolutionStepPDAE{DT,TT,AT,ΛT}, Δq::Union{AT,AbstractArray}, Δp::Union{AT,AbstractArray}, λ::Union{ΛT,AbstractArray}, μ::Union{ΛT,AbstractArray}) where {DT,TT,AT,ΛT} + update!(solstep, Δq, Δp, λ) + copyto!(solstep.μ, μ) + return solstep +end + +function update!(solstep::SolutionStepPDAE{DT,TT,AT,ΛT}, q̇::Union{AT,AbstractArray}, ṗ::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::SolutionStepPDAE{DT,TT,AT,ΛT}, q̇::Union{AT,AbstractArray}, ṗ::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::SolutionStepPDAE, Δ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).p[i] .= history(solstep).p[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).f[i] .= history(solstep).f[i-1] history(solstep).u[i] .= history(solstep).u[i-1] @@ -246,6 +268,7 @@ function Base.copy!(solstep::SolutionStepPDAE, sol::NamedTuple) solstep.q .= sol.q solstep.p .= sol.p solstep.λ .= sol.λ + solstep.μ .= sol.μ solstep.q̃ .= 0 solstep.p̃ .= 0 solstep.v .= 0 @@ -258,6 +281,7 @@ function Base.copy!(solstep::SolutionStepPDAE, sol::NamedTuple) # solstep.q̄[i] .= 0 # solstep.p̄[i] .= 0 # solstep.λ̄[i] .= 0 + # solstep.μ̄[i] .= 0 # solstep.v̄[i] .= 0 # solstep.f̄[i] .= 0 # solstep.ū[i] .= 0 diff --git a/test/solutions/deterministic_solutions_tests.jl b/test/solutions/deterministic_solutions_tests.jl index 36732f278..fb67f015b 100644 --- a/test/solutions/deterministic_solutions_tests.jl +++ b/test/solutions/deterministic_solutions_tests.jl @@ -15,7 +15,7 @@ q0 = rand(1) p0 = rand(1) z0 = rand(3) λ0 = rand(1) -μ0 = rand(2) +μ0 = rand(1) t1 = 1.0 x1 = [rand(2) for j = 1:ni] @@ -24,7 +24,7 @@ q1 = [rand(1) for j = 1:ni] p1 = [rand(1) for j = 1:ni] z1 = [rand(3) for j = 1:ni] λ1 = [rand(1) for j = 1:ni] -μ1 = [rand(2) for j = 1:ni] +μ1 = [rand(1) for j = 1:ni] t2 = t1 + (t1 - t0) x2 = [rand(2) for j = 1:ni] @@ -33,7 +33,7 @@ q2 = [rand(1) for j = 1:ni] p2 = [rand(1) for j = 1:ni] z2 = [rand(3) for j = 1:ni] λ2 = [rand(1) for j = 1:ni] -μ2 = [rand(2) for j = 1:ni] +μ2 = [rand(1) for j = 1:ni] tx = zero(x0) ty = zero(y0) @@ -49,7 +49,7 @@ qs = [rand(1) for i = 1:nt] ps = [rand(1) for i = 1:nt] zs = [rand(3) for i = 1:nt] λs = [rand(1) for i = 1:nt] -μs = [rand(2) for i = 1:nt] +μs = [rand(1) for i = 1:nt] Xs = [rand(2) for i = 1:nt, j = 1:ni] Ys = [rand(2) for i = 1:nt, j = 1:ni] @@ -57,7 +57,7 @@ Qs = [rand(1) for i = 1:nt, j = 1:ni] Ps = [rand(1) for i = 1:nt, j = 1:ni] Zs = [rand(3) for i = 1:nt, j = 1:ni] Λs = [rand(1) for i = 1:nt, j = 1:ni] -Ms = [rand(2) for i = 1:nt, j = 1:ni] +Ms = [rand(1) for i = 1:nt, j = 1:ni] ode = odeproblem() dae = daeproblem() @@ -74,7 +74,7 @@ pdae = pdaeproblem() test_interface(sol) # TODO reactivate - sol0 = Solution(similar(ode, ics=(q=x0,))) + sol0 = Solution(similar(ode, ics = (q = StateVariable(x0),))) @test typeof(sol0) <: SolutionODE @test sol != sol0 @@ -85,8 +85,8 @@ pdae = pdaeproblem() @test solstep.q == x0 # test set/get solution - sol1 = Solution(similar(ode, ics=(q=x0,))) - sol2 = Solution(similar(ode, ics=(q=x0,))) + sol1 = Solution(similar(ode, ics = (q = StateVariable(x0),))) + sol2 = Solution(similar(ode, ics = (q = StateVariable(x0),))) for i in eachindex(xs) solstep.q .= xs[i] sol1[i] = (q = copy(xs[i]),) @@ -116,7 +116,7 @@ end test_interface(sol) - sol0 = Solution(similar(pode, ics=(q=q0, p=p0))) + sol0 = Solution(similar(pode, ics = (q = StateVariable(q0), p = StateVariable(p0)))) @test typeof(sol0) <: SolutionPODE @test sol != sol0 @@ -128,8 +128,8 @@ end @test solstep.p == p0 # test set/get solution - sol1 = Solution(similar(pode, ics=(q=q0, p=p0))) - sol2 = Solution(similar(pode, ics=(q=q0, p=p0))) + sol1 = Solution(similar(pode, ics = (q = StateVariable(q0), p = StateVariable(p0)))) + sol2 = Solution(similar(pode, ics = (q = StateVariable(q0), p = StateVariable(p0)))) for i = 1:nt solstep.q .= qs[i] solstep.p .= ps[i] @@ -161,7 +161,7 @@ end test_interface(sol) - sol0 = Solution(similar(dae, ics=(q=x0, λ=λ0))) + sol0 = Solution(similar(dae, ics = (q = StateVariable(x0), λ = AlgebraicVariable(λ0), μ = AlgebraicVariable(μ0)))) @test typeof(sol0) <: SolutionDAE @test sol != sol0 @@ -173,8 +173,8 @@ end @test solstep.λ == λ0 # test set/get solution - sol1 = Solution(similar(dae, ics=(q=x0, λ=λ0))) - sol2 = Solution(similar(dae, ics=(q=x0, λ=λ0))) + sol1 = Solution(similar(dae, ics = (q = StateVariable(x0), λ = AlgebraicVariable(λ0), μ = AlgebraicVariable(μ0)))) + sol2 = Solution(similar(dae, ics = (q = StateVariable(x0), λ = AlgebraicVariable(λ0), μ = AlgebraicVariable(μ0)))) for i = 1:nt solstep.q .= xs[i] solstep.λ .= λs[i] @@ -206,7 +206,7 @@ end test_interface(sol) - sol0 = Solution(similar(pdae, ics=(q=q0, p=p0, λ=λ0))) + sol0 = Solution(similar(pdae, ics = (q = StateVariable(q0), p = StateVariable(p0), λ = AlgebraicVariable(λ0), μ = AlgebraicVariable(μ0)))) @test typeof(sol0) <: SolutionPDAE @test sol != sol0 @@ -219,8 +219,8 @@ end @test solstep.λ == λ0 # test set/get solution - sol1 = Solution(similar(pdae, ics=(q=q0, p=p0, λ=λ0))) - sol2 = Solution(similar(pdae, ics=(q=q0, p=p0, λ=λ0))) + sol1 = Solution(similar(pdae, ics = (q = StateVariable(q0), p = StateVariable(p0), λ = AlgebraicVariable(λ0), μ = AlgebraicVariable(μ0)))) + sol2 = Solution(similar(pdae, ics = (q = StateVariable(q0), p = StateVariable(p0), λ = AlgebraicVariable(λ0), μ = AlgebraicVariable(μ0)))) for i = 1:nt solstep.q .= qs[i] solstep.p .= ps[i] diff --git a/test/solutions/solution_step_tests.jl b/test/solutions/solution_step_tests.jl index 3c7c2add8..468d617b7 100644 --- a/test/solutions/solution_step_tests.jl +++ b/test/solutions/solution_step_tests.jl @@ -11,6 +11,8 @@ q0 = rand(1) p0 = q0.^2 λ0 = rand(1) λ1 = rand(1) +μ0 = rand(1) +μ1 = rand(1) Δx = rand(2) Δq = rand(1) Δp = rand(1) @@ -45,7 +47,7 @@ end @testset "$(rpad("ODE Solution Step",80))" begin - solstep = SolutionStepODE(t0, x0, parameters(ode)) + solstep = SolutionStepODE(t0, StateVariable(x0), parameters(ode)) @test solstep.t == history(solstep).t[0] == current(solstep).t @test solstep.q == history(solstep).q[0] == current(solstep).q @@ -96,7 +98,7 @@ end @testset "$(rpad("PODE Solution Step",80))" begin - solstep = SolutionStepPODE(t0, q0, p0, parameters(pode)) + solstep = SolutionStepPODE(t0, StateVariable(q0), StateVariable(p0), parameters(pode)) @test solstep.t == history(solstep).t[0] == current(solstep).t @test solstep.q == history(solstep).q[0] == current(solstep).q @@ -156,40 +158,52 @@ end @testset "$(rpad("DAE Solution Step",80))" begin - solstep = SolutionStepDAE(t0, x0, λ0, parameters(dae)) + solstep = SolutionStepDAE(t0, StateVariable(x0), AlgebraicVariable(λ0), AlgebraicVariable(μ0), parameters(dae)) @test solstep.t == history(solstep).t[0] == current(solstep).t @test solstep.q == history(solstep).q[0] == current(solstep).q @test solstep.λ == history(solstep).λ[0] == current(solstep).λ + @test solstep.μ == history(solstep).μ[0] == current(solstep).μ @test solstep.t̄ == history(solstep).t[1] == previous(solstep).t @test solstep.q̄ == history(solstep).q[1] == previous(solstep).q @test solstep.λ̄ == history(solstep).λ[1] == previous(solstep).λ + @test solstep.μ̄ == history(solstep).μ[1] == previous(solstep).μ solstep.t = t0 solstep.q .= x0 solstep.λ .= λ0 + solstep.μ .= μ0 solstep.t̄ = zero(t0) solstep.q̄ .= zero(x0) solstep.λ̄ .= zero(λ0) + solstep.μ̄ .= zero(μ0) - @test current(solstep) == (t = t0, q = x0, λ = λ0) - @test previous(solstep) == (t = zero(t0), q = zero(x0), λ = zero(λ0)) + @test current(solstep) == (t = t0, q = x0, λ = λ0, μ = μ0) + @test previous(solstep) == (t = zero(t0), q = zero(x0), λ = zero(λ0), μ = zero(μ0)) reset!(solstep, Δt) - @test current(solstep) == (t = t0 + Δt, q = x0, λ = λ0) - @test previous(solstep) == (t = t0, q = x0, λ = λ0) + @test current(solstep) == (t = t0 + Δt, q = x0, λ = λ0, μ = μ0) + @test previous(solstep) == (t = t0, q = x0, λ = λ0, μ = μ0) - update!(solstep, Δx, λ1) + update!(solstep, Δx, AlgebraicVariable(λ1)) @test solstep.t == t0 + Δt @test solstep.q == x0 .+ Δx @test solstep.λ == λ1 + @test solstep.μ == μ0 - copy!(solstep, (t = t1, q = [-2π,2π], λ = λ1)) + update!(solstep, Δx, λ1, μ1) + @test solstep.t == t0 + Δt + @test solstep.q == x0 .+ Δx .+ Δx + @test solstep.λ == λ1 + @test solstep.μ == μ1 + + copy!(solstep, (t = t1, q = [-2π,2π], λ = λ1, μ = μ1)) @test solstep.t == t1 @test solstep.q == [-2π,2π] @test solstep.λ == λ1 + @test solstep.μ == μ1 cut_periodic_solution!(solstep, (q = [2π, 0.],)) @test solstep.q == [0., 2π] @@ -213,42 +227,54 @@ end @testset "$(rpad("PDAE Solution Step",80))" begin - solstep = SolutionStepPDAE(t0, q0, p0, λ0, parameters(pdae)) + solstep = SolutionStepPDAE(t0, StateVariable(q0), StateVariable(p0), AlgebraicVariable(λ0), AlgebraicVariable(μ0), parameters(pdae)) @test solstep.t == history(solstep).t[0] == current(solstep).t @test solstep.q == history(solstep).q[0] == current(solstep).q @test solstep.p == history(solstep).p[0] == current(solstep).p @test solstep.λ == history(solstep).λ[0] == current(solstep).λ + @test solstep.μ == history(solstep).μ[0] == current(solstep).μ @test solstep.t̄ == history(solstep).t[1] == previous(solstep).t @test solstep.q̄ == history(solstep).q[1] == previous(solstep).q @test solstep.p̄ == history(solstep).p[1] == previous(solstep).p @test solstep.λ̄ == history(solstep).λ[1] == previous(solstep).λ + @test solstep.μ̄ == history(solstep).μ[1] == previous(solstep).μ solstep.t = t0 solstep.q .= q0 solstep.p .= p0 solstep.λ .= λ0 + solstep.μ .= μ0 solstep.t̄ = zero(t0) solstep.q̄ .= zero(q0) solstep.p̄ .= zero(p0) solstep.λ̄ .= zero(λ0) + solstep.μ̄ .= zero(μ0) - @test current(solstep) == (t = t0, q = q0, p = p0, λ = λ0) - @test previous(solstep) == (t = zero(t0), q = zero(q0), p = zero(p0), λ = zero(λ0)) + @test current(solstep) == (t = t0, q = q0, p = p0, λ = λ0, μ = μ0) + @test previous(solstep) == (t = zero(t0), q = zero(q0), p = zero(p0), λ = zero(λ0), μ = zero(μ0)) reset!(solstep, Δt) - @test current(solstep) == (t = t0 + Δt, q = q0, p = p0, λ = λ0) - @test previous(solstep) == (t = t0, q = q0, p = p0, λ = λ0) + @test current(solstep) == (t = t0 + Δt, q = q0, p = p0, λ = λ0, μ = μ0) + @test previous(solstep) == (t = t0, q = q0, p = p0, λ = λ0, μ = μ0) update!(solstep, Δq, Δp, λ1) @test solstep.t == t0 + Δt @test solstep.q == q0 .+ Δq @test solstep.p == p0 .+ Δp @test solstep.λ == λ1 + @test solstep.μ == μ0 + + update!(solstep, Δq, Δp, λ1, μ1) + @test solstep.t == t0 + Δt + @test solstep.q == q0 .+ Δq .+ Δq + @test solstep.p == p0 .+ Δp .+ Δp + @test solstep.λ == λ1 + @test solstep.μ == μ1 - copy!(solstep, (t = t1, q = [2π], p = [2π], λ = λ1)) + copy!(solstep, (t = t1, q = [2π], p = [2π], λ = λ1, μ = μ1)) @test solstep.t == t1 @test solstep.q == [2π] @test solstep.p == [2π]