diff --git a/Project.toml b/Project.toml index beae4b452..ba584681a 100644 --- a/Project.toml +++ b/Project.toml @@ -24,6 +24,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" RungeKutta = "fb486d5c-30a0-4a8a-8415-a8b4ace5a6f7" SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383" SimpleSolvers = "36b790f5-6434-4fbe-b711-1f64a1e2f6a2" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] @@ -45,6 +46,7 @@ QuadratureRules = "0.1" Reexport = "1.0" RungeKutta = "0.5" SimpleSolvers = "0.3" +StaticArrays = "1" julia = "1.6" [extras] diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index 4bf5d5b0b..d4729caf0 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -480,7 +480,7 @@ All implicit Runge-Kutta and partitioned Runge-Kutta methods can also be applied ### Integrators for Lagrangian ODEs Regular (non-degenerate) Lagragian ODEs can be integrated with Variational Partitioned Runge-Kutta ([`VPRK`](@ref)) -methods. +methods or Continuous Galerkin Variational Integrators ([`CGVI`](@ref)). | Function | Method | |:--------------------------------|:-------------------------------------------------------------------| diff --git a/src/Integrators.jl b/src/Integrators.jl index e1a4c3df9..dab79b2f9 100644 --- a/src/Integrators.jl +++ b/src/Integrators.jl @@ -17,6 +17,7 @@ module Integrators using RungeKutta.Tableaus using RungeKutta.PartitionedTableaus using SimpleSolvers + using StaticArrays using ..Discontinuities using ..Extrapolators @@ -121,26 +122,10 @@ module Integrators include("integrators/integrator.jl") - export ExplicitEuler, ImplicitEuler - include("integrators/euler/explicit_euler.jl") include("integrators/euler/implicit_euler.jl") - export AbstractIntegratorRK, AbstractIntegratorIRK, AbstractIntegratorPRK, IntegratorRK - - export ERKIntegrator - export IntegratorIRK - export IntegratorIRKimplicit - export IntegratorDIRK - #export IntegratorMidpointImplicit, IntegratorSRKimplicit - - export IntegratorEPRK - export IPRKIntegrator - export IntegratorIPRKimplicit - # export IntegratorFLRK - # export IntegratorPGLRK, CoefficientsPGLRK - export get_symplectic_conjugate_coefficients, symplecticize, check_symplecticity, symplecticity_conditions, check_symmetry, compute_symplecticity_error @@ -168,9 +153,6 @@ module Integrators export ExactSolution, - IntegratorExactODE, - IntegratorSplitting, - IntegratorComposition, AbstractTableauSplitting, Composition, Splitting, @@ -187,8 +169,6 @@ module Integrators include("integrators/splitting/composition_integrator.jl") - export IntegratorVPRK - include("integrators/vi/vi_methods.jl") include("integrators/vi/position_momentum_common.jl") include("integrators/vi/position_momentum_cache.jl") @@ -199,10 +179,7 @@ module Integrators include("integrators/vi/vprk_integrator.jl") - # export IntegratorCGVI, IntegratorDGVI, IntegratorDGVIEXP, - # IntegratorDGVIPI, IntegratorDGVIP0, IntegratorDGVIP1 - - # include("integrators/cgvi/integrators_cgvi.jl") + include("integrators/cgvi/integrators_cgvi.jl") # include("integrators/dgvi/integrators_dgvi.jl") # include("integrators/dgvi/integrators_dgvi_experimental.jl") # include("integrators/dgvi/integrators_dgvi_path_integral.jl") @@ -210,10 +187,6 @@ module Integrators # include("integrators/dgvi/integrators_dgvi_projection_final.jl") - export IntegratorDVIA, IntegratorDVIB, - IntegratorCMDVI, IntegratorCTDVI, - IntegratorDVRK - include("integrators/dvi/dvi_common.jl") include("integrators/dvi/dvi_cache.jl") include("integrators/dvi/dvi_euler.jl") @@ -222,8 +195,6 @@ module Integrators include("integrators/dvi/dvrk.jl") - export HPImidpointIntegrator, HPItrapezoidalIntegrator - include("integrators/hpi/hpi_common.jl") include("integrators/hpi/hpi_cache.jl") include("integrators/hpi/hpi_midpoint.jl") diff --git a/src/integrators/cgvi/integrators_cgvi.jl b/src/integrators/cgvi/integrators_cgvi.jl index 4cb0b477d..3114a741f 100644 --- a/src/integrators/cgvi/integrators_cgvi.jl +++ b/src/integrators/cgvi/integrators_cgvi.jl @@ -1,11 +1,6 @@ @doc raw""" -`ParametersCGVI`: Parameters for right-hand side function of continuous Galerkin variational Integrator. +Continuous Galerkin Variational Integrator. -### Parameters - -* `Θ`: function of the noncanonical one-form (∂L/∂v) -* `f`: function of the force (∂L/∂q) -* `Δt`: time step * `b`: weights of the quadrature rule * `c`: nodes of the quadrature rule * `x`: nodes of the basis @@ -13,76 +8,89 @@ * `a`: derivative matrix * `r₀`: reconstruction coefficients at the beginning of the interval * `r₁`: reconstruction coefficients at the end of the interval -* `t`: initial time -* `q`: solution of q at time t -* `p`: solution of p at time t + """ -mutable struct ParametersCGVI{DT, TT, D, S, R, ET <: NamedTuple} <: Parameters{DT,TT} - equs::ET - Δt::TT +struct CGVI{T, NBASIS, NNODES, NDOF, basisType <: Basis{T}} <: LODEMethod + basis::basisType + quadrature::QuadratureRule{T,NNODES} - b::Vector{TT} - c::Vector{TT} + b::SVector{NNODES,T} + c::SVector{NNODES,T} - x::Vector{TT} - m::Matrix{TT} - a::Matrix{TT} - r₀::Vector{TT} - r₁::Vector{TT} + x::SVector{NBASIS,T} - t::TT - q::Vector{DT} - p::Vector{DT} + m::SMatrix{NNODES, NBASIS, T, NDOF} + a::SMatrix{NNODES, NBASIS, T, NDOF} - function ParametersCGVI{DT,D}(equs::ET, Δt::TT, b, c, x, m, a, r₀, r₁) where {DT, TT, D, ET <: NamedTuple} - new{DT,TT,D,length(x),length(c),ET}(equs, Δt, b, c, x, m, a, r₀, r₁, zero(TT), zeros(DT,D), zeros(DT,D)) - end + r₀::SVector{NBASIS,T} + r₁::SVector{NBASIS,T} + + function CGVI(basis::Basis{T}, quadrature::QuadratureRule{T}) where {T} + # get number of quadrature nodes and number of basis functions + NNODES = QuadratureRules.nnodes(quadrature) + NBASIS = CompactBasisFunctions.nbasis(basis) + + # get quadrature nodes and weights + quad_weights = QuadratureRules.weights(quadrature) + quad_nodes = QuadratureRules.nodes(quadrature) - function ParametersCGVI{DT,D}(equs::NamedTuple, Δt::TT, basis::Basis{TT}, quadrature::QuadratureRule{TT}) where {DT,TT,D} # compute coefficients - r₀ = zeros(TT, nbasis(basis)) - r₁ = zeros(TT, nbasis(basis)) - m = zeros(TT, nnodes(quadrature), nbasis(basis)) - a = zeros(TT, nnodes(quadrature), nbasis(basis)) + r₀ = zeros(T, NBASIS) + r₁ = zeros(T, NBASIS) + m = zeros(T, NNODES, NBASIS) + a = zeros(T, NNODES, NBASIS) for i in eachindex(basis) - r₀[i] = basis[zero(TT), i] - r₁[i] = basis[one(TT), i] - for j in eachindex(quadrature) - m[j,i] = basis[nodes(quadrature)[j], i] - a[j,i] = basis'[nodes(quadrature)[j], i] + r₀[i] = basis[zero(T), i] + r₁[i] = basis[one(T), i] + for j in eachindex(quad_nodes) + m[j,i] = basis[quad_nodes[j], i] + a[j,i] = basis'[quad_nodes[j], i] end end - if get_config(:verbosity) > 1 - println() - println(" Continuous Galerkin Variational Integrator") - println(" ==========================================") - println() - println(" c = ", nodes(quadrature)) - println(" b = ", weights(quadrature)) - println(" x = ", nodes(basis)) - println(" m = ", m) - println(" a = ", a) - println(" r₀ = ", r₀) - println(" r₁ = ", r₁) - println() - end - - ParametersCGVI{DT,D}(equs, Δt, weights(quadrature), nodes(quadrature), grid(basis), m, a, r₀, r₁) + new{T, NBASIS, NNODES, NBASIS * NNODES, typeof(basis)}(basis, quadrature, quad_weights, quad_nodes, CompactBasisFunctions.grid(basis), m, a, r₀, r₁) end end - -function update_params!(params::ParametersCGVI, sol::SolutionStepPODE) - # set time for nonlinear solver and copy previous solution - params.t = sol.t - params.q .= sol.q - params.p .= sol.p +basis(method::CGVI) = method.basis +quadrature(method::CGVI) = method.quadrature + +nbasis(::CGVI{T,NB,NN}) where {T,NB,NN} = NB +nnodes(::CGVI{T,NB,NN}) where {T,NB,NN} = NN + +isexplicit(::Union{CGVI, Type{<:CGVI}}) = false +isimplicit(::Union{CGVI, Type{<:CGVI}}) = true +issymmetric(::Union{CGVI, Type{<:CGVI}}) = missing +issymplectic(::Union{CGVI, Type{<:CGVI}}) = true + +isiodemethod(::Union{CGVI, Type{<:CGVI}}) = true + +default_solver(::CGVI) = Newton() +default_iguess(::CGVI) = HermiteExtrapolation() + +function Base.show(io::IO, method::CGVI) + print(io, "\n") + print(io, " Continuous Galerkin Variational Integrator", "\n") + print(io, " ==========================================", "\n") + print(io, "\n") + print(io, " c = ", method.c, "\n") + print(io, " b = ", method.b, "\n") + print(io, " x = ", method.x, "\n") + print(io, " m = ", method.m, "\n") + print(io, " a = ", method.a, "\n") + print(io, " r₀ = ", method.r₀, "\n") + print(io, " r₁ = ", method.r₁, "\n") + print(io, "\n") end -struct IntegratorCacheCGVI{ST,D,S,R} <: IODEIntegratorCache{ST,D} +struct CGVICache{ST,D,S,R} <: IODEIntegratorCache{ST,D} + x::Vector{ST} + + q̄::Vector{ST} + p̄::Vector{ST} + q̃::Vector{ST} p̃::Vector{ST} ṽ::Vector{ST} @@ -96,7 +104,12 @@ struct IntegratorCacheCGVI{ST,D,S,R} <: IODEIntegratorCache{ST,D} F::Vector{Vector{ST}} - function IntegratorCacheCGVI{ST,D,S,R}() where {ST,D,S,R} + function CGVICache{ST,D,S,R}() where {ST,D,S,R} + x = zeros(ST, D*(S+1)) + + q̄ = zeros(ST,D) + p̄ = zeros(ST,D) + # create temporary vectors q̃ = zeros(ST,D) p̃ = zeros(ST,D) @@ -111,133 +124,58 @@ struct IntegratorCacheCGVI{ST,D,S,R} <: IODEIntegratorCache{ST,D} V = create_internal_stage_vector(ST,D,R) F = create_internal_stage_vector(ST,D,R) - new(q̃, p̃, ṽ, f̃, s̃, X, Q, P, V, F) + new(x, q̄, p̄, q̃, p̃, ṽ, f̃, s̃, X, Q, P, V, F) end end -function IntegratorCache{ST}(params::ParametersCGVI{DT,TT,D,S,R}; kwargs...) where {ST,DT,TT,D,S,R} - IntegratorCacheCGVI{ST,D,S,R}(; kwargs...) +function reset!(cache::CGVICache, t, q, p) + copyto!(cache.q̄, q) + copyto!(cache.p̄, p) end -@inline CacheType(ST, params::ParametersCGVI{DT,TT,D,S,R}) where {DT,TT,D,S,R} = IntegratorCacheCGVI{ST,D,S,R} - - -"Continuous Galerkin Variational Integrator." -struct IntegratorCGVI{DT, TT, D, S, R, - BT <: Basis, - PT <: ParametersCGVI{DT,TT,D,S,R}, - ST <: NonlinearSolver, - IT <: InitialGuessIODE{TT}} <: IODEIntegrator{DT,TT} - basis::BT - quadrature::QuadratureRule{TT,R} - - params::PT - solver::ST - iguess::IT - caches::CacheDict{PT} - - function IntegratorCGVI(basis::BT, quadrature::QuadratureRule{TT,R}, params::ParametersCGVI{DT,TT,D,S}, - solver::ST, iguess::IT, caches) where {DT,TT,D,S,R,BT,ST,IT} - new{DT, TT, D, S, R, BT, typeof(params), ST, IT}(basis, quadrature, params, solver, iguess, caches) - end - - function IntegratorCGVI{DT,D}(equations::NamedTuple, basis::Basis{TT}, quadrature::QuadratureRule{TT,R}, Δt::TT; - interpolation=HermiteExtrapolation{DT}) where {DT,TT,D,R} - - # get number of stages - S = nbasis(basis) +nlsolution(cache::CGVICache) = cache.x - # create params - params = ParametersCGVI{DT,D}(equations, Δt, basis, quadrature) - - # create cache dict - caches = CacheDict(params) - - # create nonlinear solver - solver = create_nonlinear_solver(DT, D*(S+1), params, caches) - - # create initial guess - iguess = InitialGuessIODE(get_config(:ig_extrapolation), equations[:v̄], equations[:f̄], Δt) - - # create integrator - IntegratorCGVI(basis, quadrature, params, solver, iguess, caches) - end - - function IntegratorCGVI(problem::Union{IODEProblem{DT}, LODEProblem{DT}}, basis::Basis, quadrature::QuadratureRule; kwargs...) where {DT} - IntegratorCGVI{DT, ndims(problem)}(functions(problem), basis, quadrature, timestep(problem); kwargs...) - end +function Cache{ST}(problem::AbstractProblemIODE, method::CGVI; kwargs...) where {ST} + CGVICache{ST, ndims(problem), nbasis(method), nnodes(method)}(; kwargs...) end -@inline GeometricBase.equation(integrator::IntegratorCGVI, i::Symbol) = integrator.params.equs[i] -@inline GeometricBase.equations(integrator::IntegratorCGVI) = integrator.params.equs -@inline GeometricBase.timestep(integrator::IntegratorCGVI) = integrator.params.Δt -@inline Base.ndims(::IntegratorCGVI{DT,TT,D}) where {DT,TT,D} = D +@inline CacheType(ST, problem::AbstractProblemIODE, method::CGVI) = CGVICache{ST, ndims(problem), nbasis(method), nnodes(method)} -function IntegratorCache(int::IntegratorCGVI{DT,TT}) where {DT,TT} - IntegratorCacheCGVI{DT, TT, ndims(int), nbasis(int.basis), nnodes(int.quadrature)}() -end - +function initial_guess!(int::GeometricIntegrator{<:CGVI}) + # set some local variables for convenience + local D = ndims(int) + local S = nbasis(method(int)) + local x = nlsolution(int) -function initialize!(int::IntegratorCGVI, sol::SolutionStepPODE) - sol.t̄ = sol.t - timestep(int) - - equation(int, :v̄)(sol.v, sol.t, sol.q) - equation(int, :f̄)(sol.f, sol.t, sol.q, sol.v) - - initialize!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f, - sol.t̄, sol.q̄, sol.p̄, sol.v̄, sol.f̄) -end - - -function initial_guess!(int::IntegratorCGVI{DT,TT,D,S,R}, sol::SolutionStepPODE{DT,TT}, - cache::IntegratorCacheCGVI{DT}=int.caches[DT]) where {DT,TT,D,S,R} - int.solver.x .= 0 - - for i in eachindex(int.basis) - evaluate!(int.iguess, sol.q̄, sol.p̄, sol.v̄, sol.f̄, - sol.q, sol.p, sol.v, sol.f, - cache.q̃, cache.p̃, - cache.ṽ, cache.f̃, - grid(int.basis)[i], grid(int.basis)[i]) + for i in eachindex(basis(method(int))) + initialguess!(method(int).x[i], cache(int).q̃, cache(int).p̃, solstep(int), problem(int), iguess(int)) for k in 1:D - int.solver.x[D*(i-1)+k] = cache.q̃[k] + x[D*(i-1)+k] = cache(int).q̃[k] end end - evaluate!(int.iguess, sol.q̄, sol.p̄, sol.v̄, sol.f̄, - sol.q, sol.p, sol.v, sol.f, - cache.q̃, cache.p̃, - one(TT), one(TT)) + initialguess!(one(timestep(int)), cache(int).q̃, cache(int).p̃, solstep(int), problem(int), iguess(int)) for k in 1:D - int.solver.x[D*S+k] = cache.p̃[k] + x[D*S+k] = cache(int).p̃[k] end end -"Compute stages of variational partitioned Runge-Kutta methods." -function residual!(x::Vector{ST}, b::Vector{ST}, params::ParametersCGVI{DT,TT,D,S,R}, - caches::CacheDict) where {ST,DT,TT,D,S,R} - @assert length(x) == length(b) +function components!(x::AbstractVector{ST}, int::GeometricIntegrator{<:CGVI}) where {ST} + # set some local variables for convenience and clarity + local D = ndims(int) + local S = nbasis(method(int)) + local q = cache(int, ST).q̃ + local p = cache(int, ST).p̃ + local Q = cache(int, ST).Q + local V = cache(int, ST).V + local P = cache(int, ST).P + local F = cache(int, ST).F + local X = cache(int, ST).X - # get cache for internal stages - cache = caches[ST] - - # compute stages from nonlinear solver solution x - components!(x, cache, params) - - # compute rhs b of nonlinear solver - compute_rhs!(b, cache.X, cache.Q, cache.P, cache.F, cache.p̃, params) -end - - -function components!(x, cache::IntegratorCacheCGVI, params::ParametersCGVI) - components!(x, cache.X, cache.Q, cache.V, cache.P, cache.F, cache.q̃, cache.p̃, params) -end - -function components!(x, X, Q, V, P, F, q, p, params::ParametersCGVI{DT,TT,D,S,R}) where {DT,TT,D,S,R} # copy x to X for i in eachindex(X) @@ -251,31 +189,12 @@ function components!(x, X, Q, V, P, F, q, p, params::ParametersCGVI{DT,TT,D,S,R} p[k] = x[D*S+k] end - # compute Q - compute_stages_q!(X, Q, q, params) - - # compute V - compute_stages_v!(X, V, params) - - # compute P and F - compute_stages_p!(Q, V, P, F, params) -end - - -function compute_stages_q!(X::Vector{Vector{ST}}, Q::Vector{Vector{ST}}, q::Vector{ST}, params::ParametersCGVI{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R} - @assert R == length(Q) - @assert S == length(X) - - local y::ST - # compute Q for i in eachindex(Q) - @assert D == length(Q[i]) for k in eachindex(Q[i]) - y = 0 + y = zero(ST) for j in eachindex(X) - @assert D == length(X[j]) - y += params.m[i,j] * X[j][k] + y += method(int).m[i,j] * X[j][k] end Q[i][k] = y end @@ -283,112 +202,108 @@ function compute_stages_q!(X::Vector{Vector{ST}}, Q::Vector{Vector{ST}}, q::Vect # compute q for k in eachindex(q) - y = 0 + y = zero(ST) for i in eachindex(X) - y += params.r₁[i] * X[i][k] + y += method(int).r₁[i] * X[i][k] end q[k] = y end -end - - -function compute_stages_v!(X::Vector{Vector{ST}}, V::Vector{Vector{ST}}, params::ParametersCGVI{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R} - @assert R == length(V) - @assert S == length(X) - - local y::ST # compute V for i in eachindex(V) - @assert D == length(V[i]) for k in eachindex(V[i]) - y = 0 + y = zero(ST) for j in eachindex(X) - @assert D == length(X[j]) - y += params.a[i,j] * X[j][k] + y += method(int).a[i,j] * X[j][k] end - V[i][k] = y / params.Δt + V[i][k] = y / timestep(int) end end -end - - -function compute_stages_p!(Q::Vector{Vector{ST}}, V::Vector{Vector{ST}}, - P::Vector{Vector{ST}}, F::Vector{Vector{ST}}, - params::ParametersCGVI{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R} - - @assert R == length(Q) == length(V) == length(P) == length(F) - - local tᵢ::TT # compute P=ϑ(Q,V) and F=f(Q,V) for i in eachindex(Q,V,P,F) - @assert D == length(Q[i]) == length(V[i]) == length(P[i]) == length(F[i]) - tᵢ = params.t + params.Δt * params.c[i] - params.equs[:ϑ](P[i], tᵢ, Q[i], V[i]) - params.equs[:f](F[i], tᵢ, Q[i], V[i]) + tᵢ = solstep(int).t + timestep(int) * method(int).c[i] + equations(int).ϑ(P[i], tᵢ, Q[i], V[i], parameters(solstep(int))) + equations(int).f(F[i], tᵢ, Q[i], V[i], parameters(solstep(int))) end end -function compute_rhs!(b::Vector{ST}, X::Vector{Vector{ST}}, Q::Vector{Vector{ST}}, - P::Vector{Vector{ST}}, F::Vector{Vector{ST}}, - p::Vector{ST}, params::ParametersCGVI{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R} - local y::ST - local z::ST +function residual!(b::Vector{ST}, int::GeometricIntegrator{<:CGVI}) where {ST} + # set some local variables for convenience and clarity + local D = ndims(int) + local S = nbasis(method(int)) + local q̄ = cache(int, ST).q̄ + local p̄ = cache(int, ST).p̄ + local p̃ = cache(int, ST).p̃ + local P = cache(int, ST).P + local F = cache(int, ST).F + local X = cache(int, ST).X # compute b = - [(P-AF)] - for i in 1:S - for k in 1:D - z = 0 + for i in eachindex(method(int).r₀, method(int).r₁) + for k in eachindex(p̃, p̄) + z = zero(ST) for j in eachindex(P,F) - z += params.b[j] * params.m[j,i] * F[j][k] - z += params.b[j] * params.a[j,i] * P[j][k] / params.Δt + z += method(int).b[j] * method(int).m[j,i] * F[j][k] * timestep(int) + z += method(int).b[j] * method(int).a[j,i] * P[j][k] end - b[D*(i-1)+k] = z - (params.r₁[i] * p[k] - params.r₀[i] * params.p[k]) / params.Δt + b[D*(i-1)+k] = (method(int).r₁[i] * p̃[k] - method(int).r₀[i] * p̄[k]) - z end end # compute b = - [(q-r₀Q)] - for k in eachindex(params.q) - y = 0 + for k in eachindex(q̄) + y = zero(ST) for j in eachindex(X) - y += params.r₀[j] * X[j][k] + y += method(int).r₀[j] * X[j][k] end - b[D*S+k] = y - params.q[k] + b[D*S+k] = q̄[k] - y end end -function update_solution!(sol::SolutionStepPODE, cache::IntegratorCacheCGVI) - sol.q .= cache.q̃ - sol.p .= cache.p̃ +# Compute stages of Variational Partitioned Runge-Kutta methods. +function residual!(b::AbstractVector{ST}, x::AbstractVector{ST}, int::GeometricIntegrator{<:CGVI}) where {ST} + @assert axes(x) == axes(b) + + # copy previous solution from solstep to cache + reset!(cache(int, ST), current(solstep(int))...) + + # compute stages from nonlinear solver solution x + components!(x, int) + + # compute residual vector + residual!(b, int) end -function integrate_step!(int::IntegratorCGVI{DT,TT}, sol::SolutionStepPODE{DT,TT}, - cache::IntegratorCacheCGVI{DT}=int.caches[DT]) where {DT,TT} - # update nonlinear solver parameters from cache - update_params!(int.params, sol) +function update!(x::AbstractVector{DT}, int::GeometricIntegrator{<:CGVI}) where {DT} + # copy previous solution from solstep to cache + reset!(cache(int, DT), current(solstep(int))...) + + # compute vector field at internal stages + components!(x, int) - # compute initial guess - initial_guess!(int, sol, cache) + # compute final update + solstep(int).q .= cache(int, DT).q̃ + solstep(int).p .= cache(int, DT).p̃ +end - # reset cache - reset!(sol) + +function integrate_step!(int::GeometricIntegrator{<:CGVI, <:AbstractProblemIODE}) + # copy previous solution from solstep to cache + reset!(cache(int), current(solstep(int))...) # call nonlinear solver - solve!(int.solver) + solve!(nlsolution(int), (b,x) -> residual!(b, x, int), solver(int)) # print solver status - print_solver_status(int.solver.status, int.solver.params) + # print_solver_status(int.solver.status, int.solver.params) # check if solution contains NaNs or error bounds are violated - check_solver_status(int.solver.status, int.solver.params) - - # compute vector fields at internal stages - components!(int.solver.x, cache, int.params) + # check_solver_status(int.solver.status, int.solver.params) # compute final update - update_solution!(sol, cache) + update!(nlsolution(int), int) end diff --git a/src/integrators/method_list.jl b/src/integrators/method_list.jl index 6c8fcd423..80fd92fdc 100644 --- a/src/integrators/method_list.jl +++ b/src/integrators/method_list.jl @@ -8,11 +8,17 @@ meta_methods = ( # FLRK, VPRK, DVRK, + CGVI, ProjectedMethod, HPImidpoint, HPItrapezoidal, ) +euler_methods = ( + ExplicitEuler, + ImplicitEuler, +) + explicit_rungekutta_methods = ( ForwardEuler, ExplicitEulerRK, @@ -126,8 +132,10 @@ variational_partitioned_runge_kutta_families = ( ) variational_integrators = ( + HPImidpoint, + HPItrapezoidal, PMVImidpoint, - PMVItrapezoidal, + PMVItrapezoidal, ) degenerate_variational_integrators = ( @@ -152,6 +160,7 @@ splitting_methods = ( method_groups = ( meta_methods, + euler_methods, runge_kutta_methods, runge_kutta_families, partitioned_runge_kutta_methods, diff --git a/test/integrators/galerkin_integrators_tests.jl b/test/integrators/galerkin_integrators_tests.jl index 255a60220..a7490b0c8 100644 --- a/test/integrators/galerkin_integrators_tests.jl +++ b/test/integrators/galerkin_integrators_tests.jl @@ -1,13 +1,12 @@ using GeometricIntegrators using GeometricProblems.HarmonicOscillator -using Test - using CompactBasisFunctions using QuadratureRules +using Test -using GeometricProblems.HarmonicOscillator: reference_solution iode = iodeproblem() +pref = exact_solution(podeproblem()) QGau4 = GaussLegendreQuadrature(4) BGau4 = Lagrange(QuadratureRules.nodes(QGau4)) @@ -15,29 +14,23 @@ BGau4 = Lagrange(QuadratureRules.nodes(QGau4)) ### CGVI Integrators ### -cgint = IntegratorCGVI(iode, BGau4, QGau4) -cgsol = integrate(iode, cgint) -@test relative_maximum_error(cgsol.q, reference_solution) < 1E-7 +cgsol = integrate(iode, CGVI(BGau4, QGau4)) +@test relative_maximum_error(cgsol.q, pref.q) < 1E-7 ### DGVI Integrators ### -dgint = IntegratorDGVI(iode, BGau4, QGau4) -dgsol = integrate(iode, dgint) -@test relative_maximum_error(dgsol.q, reference_solution) < 1E-7 +# dgsol = integrate(iode, DGVI(BGau4, QGau4)) +# @test relative_maximum_error(dgsol.q, pref.q) < 1E-7 -dgint = IntegratorDGVIP0(iode, BGau4, QGau4) -dgsol = integrate(iode, dgint) -@test relative_maximum_error(dgsol.q, reference_solution) < 1E-7 +# dgsol = integrate(iode, DGVIP0(BGau4, QGau4)) +# @test relative_maximum_error(dgsol.q, pref.q) < 1E-7 -dgint = IntegratorDGVIP1(iode, BGau4, QGau4) -dgsol = integrate(iode, dgint) -@test relative_maximum_error(dgsol.q, reference_solution) < 1E-7 +# dgsol = integrate(iode, DGVIP1(BGau4, QGau4)) +# @test relative_maximum_error(dgsol.q, pref.q) < 1E-7 -dgint = IntegratorDGVIEXP(iode, BGau4, QGau4) -dgsol = integrate(iode, dgint) -@test relative_maximum_error(dgsol.q, reference_solution) < 1E-7 +# dgsol = integrate(iode, DGVIEXP(BGau4, QGau4)) +# @test relative_maximum_error(dgsol.q, pref.q) < 1E-7 -dgint = IntegratorDGVIPI(iode, BGau4, QGau4, Discontinuity(PathIntegralLinear(), LobattoLegendreQuadrature(2))) -dgsol = integrate(iode, dgint) -@test relative_maximum_error(dgsol.q, reference_solution) < 1E-7 +# dgsol = integrate(iode, DGVIPI(BGau4, QGau4, Discontinuity(PathIntegralLinear(), LobattoLegendreQuadrature(2)))) +# @test relative_maximum_error(dgsol.q, pref.q) < 1E-7 diff --git a/test/runtests.jl b/test/runtests.jl index e5b756d27..237ba048b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,15 +3,15 @@ using SafeTestsets @safetestset "Extrapolation Methods " begin include("extrapolation_tests.jl") end @safetestset "Solution Tests " begin include("solutions/solutions_tests.jl") end @safetestset "Initial Guesses " begin include("integrators/initial_guess_tests.jl") end + @safetestset "Euler Integrators " begin include("integrators/euler_tests.jl") end @safetestset "Runge-Kutta Integrators " begin include("integrators/rk_integrators_tests.jl") end @safetestset "Runge-Kutta Integrators for Implicit Equations " begin include("integrators/rk_implicit_integrators_tests.jl") end @safetestset "Splitting Integrators " begin include("integrators/splitting_integrators_tests.jl") end + @safetestset "Variational Integrators " begin include("integrators/variational_integrators_tests.jl") end @safetestset "Degenerate Variational Integrators " begin include("integrators/dvi_integrators_tests.jl") end - -# @safetestset "Galerkin Variational Integrators " begin include("integrators/galerkin_integrators_tests.jl") end - +@safetestset "Galerkin Variational Integrators " begin include("integrators/galerkin_integrators_tests.jl") end @safetestset "Hamilton-Pontryagin Integrators " begin include("integrators/hamilton_pontryagin_integrators_tests.jl") end @safetestset "Projection Methods " begin include("projections/projections_tests.jl") end @@ -21,8 +21,7 @@ using SafeTestsets @safetestset "Ensemble Integrator Tests " begin include("integrators/ensemble_integrators_tests.jl") end @safetestset "Method Tests " begin include("methods/methods_tests.jl") end +@safetestset "Simulation Tests " begin include("simulations/simulations_tests.jl") end @safetestset "SPARK Integrators " begin include("spark/spark_integrators_tests.jl") end @safetestset "SPARK Tableau Tests " begin include("spark/spark_tableaus_tests.jl") end - -@safetestset "Simulation Tests " begin include("simulations/simulations_tests.jl") end