Update for GeometricIntegrators v0.8.0's revision of equation types (…

…and some minor fixes along the way).

Project.toml

@@ -11,6 +11,7 @@ Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 Documenter = "0.23, 0.24, 0.25, 0.26"
@@ -20,6 +21,7 @@ Parameters = "0.12"
 Reexport = "0.2, 1.0"
 Requires = "1"
 RuntimeGeneratedFunctions = "0.4, 0.5"
+Symbolics = "0.1"
 julia = "^1.5"
 
 [extras]

src/diagnostics.jl

@@ -86,7 +86,7 @@ module Diagnostics
         invds = SDataSeries(T, ntime(q))
         try
             for i in eachindex(invds)
-                invds[i] = invariant(t[i], q[:,i], p[:,i])
+                invds[i] = invariant(t[i], q[i], p[i])
             end
         catch ex
             if isa(ex, DomainError)

src/harmonic_oscillator.jl

@@ -22,10 +22,13 @@ module HarmonicOscillator
     const p = (k=k, ω=ω)
 
 
+    ϑ₁(t,q) = q[2]
+    ϑ₂(t,q) = zero(eltype(q))
+
     function ϑ(q)
         p = zero(q)
-        p[1] = q[2]
-        p[2] = 0
+        p[1] = ϑ₁(0,q)
+        p[2] = ϑ₂(0,q)
         return p
     end
 
@@ -39,6 +42,10 @@ module HarmonicOscillator
         p[1]^2 / 2 + k * q[1]^2 / 2
     end
 
+    function lagrangian(t, q, v, params)
+        ϑ₁(t,q) * v[1] + ϑ₂(t,q) * v[2] - hamiltonian(t, q, params)
+    end
+
 
     t₀=0.0
     q₀=[0.5, 0.0]
@@ -62,7 +69,7 @@ module HarmonicOscillator
 
     function harmonic_oscillator_ode(x₀=q₀, params=p)
         # @assert size(x₀,1) == 2
-        ODE(oscillator_ode_v, x₀; parameters=params, h=hamiltonian)
+        ODE(oscillator_ode_v, x₀; invariants=(h=hamiltonian,), parameters=params)
     end
 
 
@@ -81,7 +88,7 @@ module HarmonicOscillator
         # @assert length(q₀) == length(p₀)
         # @assert all([length(q) == length(p) == 1 for (q,p) in zip(q₀,p₀)])
         # @assert size(q₀,1) == size(p₀,1) == 1
-        PODE(oscillator_pode_v, oscillator_pode_f, q₀, p₀; parameters=params)
+        PODE(oscillator_pode_v, oscillator_pode_f, q₀, p₀; invariants=(h=hamiltonian,), parameters=params)
     end
 
 
@@ -160,7 +167,7 @@ module HarmonicOscillator
         # @assert size(q₀,1) == size(p₀,1) == 2
         IODE(oscillator_iode_ϑ, oscillator_iode_f,
              oscillator_iode_g, q₀, p₀;
-             parameters=params,
+             invariants=(h=hamiltonian,), parameters=params,
              v̄=oscillator_iode_v)
     end
 
@@ -188,8 +195,8 @@ module HarmonicOscillator
         # @assert size(λ₀,1) == 1
         # @assert all([length(z) == 3 for z in z₀])
         # @assert all([length(λ) == 1 for λ in λ₀])
-        DAE(oscillator_dae_v, oscillator_dae_u, oscillator_dae_ϕ, z₀, λ₀;
-            v̄=oscillator_ode_v, parameters=params)
+        DAE(oscillator_dae_v, oscillator_dae_u, oscillator_dae_ϕ,
+            z₀, λ₀; invariants=(h=hamiltonian,), parameters=params, v̄=oscillator_ode_v)
     end
 
 
@@ -232,9 +239,8 @@ module HarmonicOscillator
     function harmonic_oscillator_idae(q₀=q₀, p₀=ϑ(q₀), λ₀=zero(q₀), params=p)
         # @assert size(q₀,1) == size(p₀,1) == size(λ₀,1) == 2
         IDAE(oscillator_iode_ϑ, oscillator_iode_f,
-             oscillator_idae_u, oscillator_idae_g,
-             oscillator_idae_ϕ, q₀, p₀, λ₀;
-             parameters=params, v̄=oscillator_iode_v)
+             oscillator_idae_u, oscillator_idae_g, oscillator_idae_ϕ,
+             q₀, p₀, λ₀; invariants=(h=hamiltonian,), parameters=params, v̄=oscillator_iode_v)
     end
 
     function oscillator_pdae_v(t, q, p, v, params)
@@ -254,19 +260,16 @@ module HarmonicOscillator
     function harmonic_oscillator_pdae(q₀=q₀, p₀=ϑ(q₀), λ₀=zero(q₀), params=p)
         # @assert size(q₀,1) == size(p₀,1) == 2
         PDAE(oscillator_pdae_v, oscillator_pdae_f,
-             oscillator_idae_u, oscillator_idae_g,
-             oscillator_idae_ϕ, q₀, p₀, λ₀;
-             parameters=params)
+             oscillator_idae_u,  oscillator_idae_g, oscillator_idae_ϕ,
+             q₀, p₀, λ₀; invariants=(h=hamiltonian,), parameters=params)
     end
 
     function harmonic_oscillator_hdae(q₀=q₀, p₀=ϑ(q₀), λ₀=zero(q₀), params=p)
         # @assert size(q₀,1) == size(p₀,1) == 2
-        HDAE(oscillator_pdae_v, oscillator_pdae_f,
-             oscillator_idae_u, oscillator_idae_g,
-             oscillator_idae_ū, oscillator_idae_ḡ,
-             oscillator_idae_ϕ, oscillator_idae_ψ,
-             hamiltonian, q₀, p₀, λ₀;
-             parameters=params)
+        HDAE(oscillator_pdae_v, oscillator_pdae_f, 
+             oscillator_idae_u, oscillator_idae_g, oscillator_idae_ϕ,
+             oscillator_idae_ū, oscillator_idae_ḡ, oscillator_idae_ψ,
+             hamiltonian, q₀, p₀, λ₀; parameters=params)
     end
 
 

src/lotka_volterra_2d_common.jl

@@ -55,7 +55,14 @@ function hamiltonian(t, q, params)
     a₁*q[1] + a₂*q[2] + b₁*log(q[1]) + b₂*log(q[2])
 end
 
-hamiltonian(t, q, p, params) = hamiltonian(t, q, params)
+hamiltonian_iode(t, q, v, params) = hamiltonian(t, q, params)
+
+hamiltonian_pode(t, q, p, params) = hamiltonian(t, q, params)
+
+function lagrangian(t, q, v, params)
+    ϑ₁(t,q) * v[1] + ϑ₂(t,q) * v[2] - hamiltonian(t, q, params)
+end
+
 
 function dHd₁(t, q, params)
     @unpack a₁, b₁ = params
@@ -89,6 +96,7 @@ g₂(t, q, v) = dϑ₂dx₁(t,q) * v[1] + dϑ₂dx₂(t,q) * v[2]
 lotka_volterra_2d_ϑ(t, q, Θ, params) = ϑ(t, q, Θ)
 lotka_volterra_2d_ϑ(t, q, v, Θ, params) = ϑ(t, q, Θ)
 lotka_volterra_2d_ω(t, q, Ω, params) = ω(t, q, Ω)
+lotka_volterra_2d_ω(t, q, v, Ω, params) = ω(t, q, Ω)
 
 
 function lotka_volterra_2d_dH(t, q, dH, params)
@@ -200,3 +208,9 @@ function lotka_volterra_2d_ψ(t, q, p, v, f, ψ, params)
     ψ[2] = f[2] - g₂(t,q,v)
     nothing
 end
+
+function lotka_volterra_2d_ψ_lode(t, q, p, v, f, ψ, params)
+    ψ[1] = f₁(t,q,v) - g₁(t,q,v) - dHd₁(t, q, params)
+    ψ[2] = f₂(t,q,v) - g₂(t,q,v) - dHd₂(t, q, params)
+    nothing
+end

src/lotka_volterra_2d_equations.jl

@@ -61,101 +61,97 @@ compute_energy_error(t,q,params) = compute_invariant_error(t,q, (t,q) -> hamilto
 
 "Creates an ODE object for the Lotka-Volterra 2D model."
 function lotka_volterra_2d_ode(q₀=q₀; params=parameters)
-    ODE(lotka_volterra_2d_v, q₀; parameters=params, h=hamiltonian)
+    ODE(lotka_volterra_2d_v, q₀; parameters=params, invariants=(h=hamiltonian,))
 end
 
 "Creates a Hamiltonian ODE object for the Lotka-Volterra 2D model."
 function lotka_volterra_2d_hode(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀); params=parameters)
-    HODE(lotka_volterra_2d_v, lotka_volterra_2d_f, hamiltonian, q₀, p₀;
+    HODE(lotka_volterra_2d_v, lotka_volterra_2d_f, hamiltonian_pode, q₀, p₀;
          parameters=params)
 end
 
 "Creates an implicit ODE object for the Lotka-Volterra 2D model."
 function lotka_volterra_2d_iode(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀); params=parameters)
     IODE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f,
          lotka_volterra_2d_g, q₀, p₀;
-         parameters=params, h=hamiltonian, v̄=lotka_volterra_2d_v)
+         parameters=params, invariants=(h=hamiltonian_iode,), v̄=lotka_volterra_2d_v)
 end
 
 "Creates a partitioned ODE object for the Lotka-Volterra 2D model."
 function lotka_volterra_2d_pode(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀); params=parameters)
-    PODE(lotka_volterra_2d_v, lotka_volterra_2d_f,
-         q₀, p₀; parameters=params, h=hamiltonian)
+    PODE(lotka_volterra_2d_v, lotka_volterra_2d_f, q₀, p₀;
+         parameters=params, invariants=(h=hamiltonian_pode,))
 end
 
 "Creates a variational ODE object for the Lotka-Volterra 2D model."
 function lotka_volterra_2d_lode(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀); params=parameters)
     LODE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f,
-         lotka_volterra_2d_g, q₀, p₀;
-         parameters=params, h=hamiltonian, v̄=lotka_volterra_2d_v,
-         Ω=lotka_volterra_2d_ω, ∇H=lotka_volterra_2d_dH)
+         lotka_volterra_2d_g, lagrangian, lotka_volterra_2d_ω, q₀, p₀;
+         parameters=params, invariants=(h=hamiltonian_iode,), v̄=lotka_volterra_2d_v)
 end
 
 "Creates a DAE object for the Lotka-Volterra 2D model."
 function lotka_volterra_2d_dae(q₀=vcat(q₀,v₀), λ₀=zero(q₀); params=parameters)
     DAE(lotka_volterra_2d_v_dae, lotka_volterra_2d_u_dae, lotka_volterra_2d_ϕ_dae, q₀, λ₀;
-        parameters=params, h=hamiltonian)
+        parameters=params, invariants=(h=hamiltonian,))
 end
 
 "Creates a Hamiltonian DAE object for the Lotka-Volterra 2D model."
 function lotka_volterra_2d_hdae(q₀=q₀, p₀=ϑ(0, q₀), λ₀=zero(q₀); params=parameters)
-    HDAE(lotka_volterra_2d_v, lotka_volterra_2d_f,
-         lotka_volterra_2d_u, lotka_volterra_2d_g,
-         lotka_volterra_2d_u̅, lotka_volterra_2d_g̅,
-         lotka_volterra_2d_ϕ, lotka_volterra_2d_ψ,
-         hamiltonian, q₀, p₀, λ₀;
-         parameters=params)
+    HDAE(lotka_volterra_2d_v, lotka_volterra_2d_f, 
+         lotka_volterra_2d_u, lotka_volterra_2d_g, lotka_volterra_2d_ϕ,
+         lotka_volterra_2d_u̅, lotka_volterra_2d_g̅, lotka_volterra_2d_ψ,
+         hamiltonian_pode, q₀, p₀, λ₀; parameters=params)
 end
 
 "Creates an implicit DAE object for the Lotka-Volterra 2D model."
 function lotka_volterra_2d_idae(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀), λ₀=zero(q₀); params=parameters)
     IDAE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f,
-         lotka_volterra_2d_u, lotka_volterra_2d_g,
-         lotka_volterra_2d_ϕ, q₀, p₀, λ₀;
-         parameters=params, h=hamiltonian, v̄=lotka_volterra_2d_v)
+         lotka_volterra_2d_u, lotka_volterra_2d_g, lotka_volterra_2d_ϕ,
+         q₀, p₀, λ₀; parameters=params, invariants=(h=hamiltonian_iode,),
+         v̄=lotka_volterra_2d_v)
 end
 
 "Creates an implicit DAE object for the Lotka-Volterra 2D model."
 function lotka_volterra_2d_idae_spark(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀), λ₀=zero(q₀); params=parameters)
     IDAE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f_ham,
-         lotka_volterra_2d_u, lotka_volterra_2d_g,
-         lotka_volterra_2d_ϕ, q₀, p₀, λ₀;
-         parameters=params, h=hamiltonian,
+         lotka_volterra_2d_u, lotka_volterra_2d_g, lotka_volterra_2d_ϕ,
+         q₀, p₀, λ₀; parameters=params, invariants=(h=hamiltonian_iode,),
          v̄=lotka_volterra_2d_v, f̄=lotka_volterra_2d_f)
 end
 
 "Creates a partitioned DAE object for the Lotka-Volterra 2D model."
 function lotka_volterra_2d_pdae(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀), λ₀=zero(q₀); params=parameters)
     PDAE(lotka_volterra_2d_v_ham, lotka_volterra_2d_f_ham,
-         lotka_volterra_2d_u, lotka_volterra_2d_g,
-         lotka_volterra_2d_ϕ, q₀, p₀, λ₀;
-         parameters=params, h=hamiltonian,
+         lotka_volterra_2d_u, lotka_volterra_2d_g, lotka_volterra_2d_ϕ,
+         q₀, p₀, λ₀; parameters=params, invariants=(h=hamiltonian_pode,),
          v̄=lotka_volterra_2d_v, f̄=lotka_volterra_2d_f)
 end
 
 "Creates a variational DAE object for the Lotka-Volterra 2D model."
 function lotka_volterra_2d_ldae(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀), λ₀=zero(q₀); params=parameters)
     LDAE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f_ham,
-         lotka_volterra_2d_g, lotka_volterra_2d_g̅,
-         lotka_volterra_2d_ϕ, lotka_volterra_2d_ψ,
-         q₀, p₀, λ₀; parameters=params, h=hamiltonian,
+         lotka_volterra_2d_u, lotka_volterra_2d_g, lotka_volterra_2d_ϕ,
+         lotka_volterra_2d_u̅, lotka_volterra_2d_g̅, lotka_volterra_2d_ψ_lode,
+         lagrangian, lotka_volterra_2d_ω,
+         q₀, p₀, λ₀; parameters=params, invariants=(h=hamiltonian_iode,),
          v̄=lotka_volterra_2d_v, f̄=lotka_volterra_2d_f)
 end
 
 "Creates a variational DAE object for the Lotka-Volterra 2D model for use with SLRK integrators."
 function lotka_volterra_2d_slrk(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀), λ₀=zero(q₀); params=parameters)
     LDAE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f,
-            lotka_volterra_2d_g, lotka_volterra_2d_g̅,
-            lotka_volterra_2d_ϕ, lotka_volterra_2d_ψ,
-            q₀, p₀, λ₀; parameters=params, h=hamiltonian,
+            lotka_volterra_2d_u, lotka_volterra_2d_g, lotka_volterra_2d_ϕ,
+            lotka_volterra_2d_u̅, lotka_volterra_2d_g̅, lotka_volterra_2d_ψ,
+            lagrangian, lotka_volterra_2d_ω,
+            q₀, p₀, λ₀; parameters=params, invariants=(h=hamiltonian_iode,),
             v̄=lotka_volterra_2d_v, f̄=lotka_volterra_2d_f)
 end
 
 "Creates an implicit ODE object for the Lotka-Volterra 2D model for use with DG integrators."
 function lotka_volterra_2d_dg(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀); params=parameters)
-    IODE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f,
-            lotka_volterra_2d_g, q₀, p₀;
-            parameters=params, h=hamiltonian, v̄=lotka_volterra_2d_v)
+    IODE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f, lotka_volterra_2d_g,
+         q₀, p₀; parameters=params, invariants=(h=hamiltonian_iode,), v̄=lotka_volterra_2d_v)
 end
 
 

src/lotka_volterra_3d.jl

@@ -86,6 +86,8 @@ module LotkaVolterra3d
         A1*q[1] + A2*q[2] + A3*q[3] - B1*log(q[1]) - B2*log(q[2]) - B3*log(q[3])
     end
 
+    hamiltonian_iode(t, q, v) = hamiltonian(t, q)
+
     function casimir(t, q)
         log(q[1]) + log(q[2]) + log(q[3])
     end
@@ -100,7 +102,7 @@ module LotkaVolterra3d
 
 
     function lotka_volterra_3d_ode(q₀=q₀)
-        ODE(lotka_volterra_3d_v, q₀; h=hamiltonian)
+        ODE(lotka_volterra_3d_v, q₀; invariants=(h=hamiltonian,))
     end
 
 

src/lotka_volterra_4d.jl

@@ -228,9 +228,9 @@ module LotkaVolterra4d
         a₁*q[1] + a₂*q[2] + a₃*q[3] + a₄*q[4] + b₁*log(q[1]) + b₂*log(q[2]) + b₃*log(q[3]) + b₄*log(q[4])
     end
 
-    function hamiltonian(t, q, p, params)
-        hamiltonian(t, q, params)
-    end
+    hamiltonian_iode(t, q, v, params) = hamiltonian(t, q, params)
+
+    hamiltonian_pode(t, q, p, params) = hamiltonian(t, q, params)
 
 
     function dHd₁(t, q, params)
@@ -356,55 +356,55 @@ module LotkaVolterra4d
 
 
     function lotka_volterra_4d_ode(q₀=q₀, params=p)
-        ODE(lotka_volterra_4d_v, q₀; parameters=params, h=hamiltonian)
+        ODE(lotka_volterra_4d_v, q₀; parameters=params, invariants=(h=hamiltonian,))
     end
 
 
     function lotka_volterra_4d_pode(q₀=q₀, p₀=ϑ(0, q₀), params=p)
         PODE(lotka_volterra_4d_v, lotka_volterra_4d_f,
-             q₀, p₀; parameters=params, h=hamiltonian)
+             q₀, p₀; parameters=params, invariants=(h=hamiltonian_pode,))
     end
 
     function lotka_volterra_4d_iode(q₀=q₀, p₀=ϑ(0, q₀), params=p)
         IODE(lotka_volterra_4d_ϑ, lotka_volterra_4d_f,
              lotka_volterra_4d_g, q₀, p₀;
-             parameters=params, h=hamiltonian, v̄=lotka_volterra_4d_v)
+             parameters=params, invariants=(h=hamiltonian_iode,), v̄=lotka_volterra_4d_v)
     end
 
     function lotka_volterra_4d_lode(q₀=q₀, p₀=ϑ(0, q₀), params=p)
         LODE(lotka_volterra_4d_ϑ, lotka_volterra_4d_f,
              lotka_volterra_4d_g, q₀, p₀;
-             parameters=params, h=hamiltonian, v̄=lotka_volterra_4d_v,
+             parameters=params, invariants=(h=hamiltonian_iode,), v̄=lotka_volterra_4d_v,
              Ω=lotka_volterra_4d_ω, ∇H=lotka_volterra_4d_dH)
     end
 
     function lotka_volterra_4d_idae(q₀=q₀, p₀=ϑ(0, q₀), λ₀=zero(q₀), params=p)
         IDAE(lotka_volterra_4d_ϑ, lotka_volterra_4d_f,
              lotka_volterra_4d_u, lotka_volterra_4d_g,
              lotka_volterra_4d_ϕ, q₀, p₀, λ₀;
-             parameters=params, h=hamiltonian, v̄=lotka_volterra_4d_v)
+             parameters=params, invariants=(h=hamiltonian_iode,), v̄=lotka_volterra_4d_v)
     end
 
     function lotka_volterra_4d_pdae(q₀=q₀, p₀=ϑ(0, q₀), λ₀=zero(q₀), params=p)
         PDAE(lotka_volterra_4d_v_ham, lotka_volterra_4d_f_ham,
              lotka_volterra_4d_u, lotka_volterra_4d_g,
              lotka_volterra_4d_ϕ, q₀, p₀, λ₀;
              v̄=lotka_volterra_4d_v, f̄=lotka_volterra_4d_f,
-             parameters=params, h=hamiltonian)
+             parameters=params, invariants=(h=hamiltonian_pode,))
     end
 
     function lotka_volterra_4d_ldae(q₀=q₀, p₀=ϑ(0, q₀), λ₀=zero(q₀), params=p)
         LDAE(lotka_volterra_4d_ϑ, lotka_volterra_4d_f_ham,
              lotka_volterra_4d_g, lotka_volterra_4d_g̅,
              lotka_volterra_4d_ϕ, lotka_volterra_4d_ψ,
-             q₀, p₀, λ₀; parameters=params, h=hamiltonian,
+             q₀, p₀, λ₀; parameters=params, invariants=(h=hamiltonian_iode,),
              v̄=lotka_volterra_4d_v, f̄=lotka_volterra_4d_f,)
     end
 
     function lotka_volterra_4d_dg(q₀=q₀, p₀=ϑ(0, q₀), params=p)
         IODE(lotka_volterra_4d_ϑ, lotka_volterra_4d_f,
              lotka_volterra_4d_g, q₀, p₀;
-             parameters=params, h=hamiltonian, v̄=lotka_volterra_4d_v)
+             parameters=params, invariants=(h=hamiltonian_iode,), v̄=lotka_volterra_4d_v)
     end
 
 end