Update for recent versions of GeometricIntegrators.

Project.toml

@@ -8,37 +8,36 @@ DecFP = "55939f99-70c6-5e9b-8bb0-5071ed7d61fd"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ElectromagneticFields = "d6c1ba6f-ee03-53af-b876-68cefeb88ec8"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-GeometricIntegrators = "dcce2d33-59f6-5b8d-9047-0defad88ae06"
+GeometricEquations = "c85262ba-a08a-430a-b926-d29770767bf2"
 GeometricProblems = "18cb22b4-ad41-5c80-9c5f-710df63fbdc9"
+GeometricSolutions = "7843afe4-64f4-4df4-9231-049495c56661"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PoincareInvariants = "26663084-47d3-540f-bd97-40ca743aafa4"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
-SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
-Showoff = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
+SimpleSolvers = "36b790f5-6434-4fbe-b711-1f64a1e2f6a2"
 
 [compat]
-DecFP = "^0.4.9, 1.0"
-Documenter = "0.23, 0.24, 0.25, 0.26, 0.27"
+DecFP = "1"
+Documenter = "1"
 ElectromagneticFields = "0.5"
 ForwardDiff = "0.10"
-GeometricIntegrators = "^0.3.1, 0.4, 0.6"
-GeometricProblems = "^0.1, 0.2"
-LaTeXStrings = "1.0, 1.1"
-Parameters = "0.12"
+GeometricEquations = "0.18"
+GeometricProblems = "0.6"
+GeometricSolutions = "0.3"
+LaTeXStrings = "1"
+Parameters = "0.12, 0.13"
 Plots = "1"
-PoincareInvariants = "0.2, 0.3"
 RecipesBase = "1"
-SafeTestsets = "0.0.1"
-Showoff = "^0.3.1, 1.0"
-julia = "^1.3"
+julia = "1.6"
 
 [extras]
+GeometricIntegrators = "dcce2d33-59f6-5b8d-9047-0defad88ae06"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["GeometricIntegrators", "SafeTestsets", "Test"]

docs/Project.toml

@@ -1,3 +1,6 @@
 [deps]
+ChargedParticleDynamics = "209b3593-7183-5600-9f8f-df5ce7fddd64"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-"ChargedParticleDynamics" = "209b3593-7183-5600-9f8f-df5ce7fddd64"
+ElectromagneticFields = "d6c1ba6f-ee03-53af-b876-68cefeb88ec8"
+GeometricIntegrators = "dcce2d33-59f6-5b8d-9047-0defad88ae06"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

docs/make.jl

@@ -3,6 +3,7 @@ using ChargedParticleDynamics
 
 makedocs(
     sitename = "ChargedParticleDynamics.jl",
+    warnonly = Documenter.except(:autodocs_block, :cross_references, :docs_block, :doctest, :eval_block, :example_block, :footnote, :linkcheck_remotes, :linkcheck, :meta_block, :parse_error, :setup_block),
     format = Documenter.HTML(
                 prettyurls = get(ENV, "CI", nothing) == "true",
                 assets = [asset("assets/style.css", class=:css, islocal=true)]),
@@ -17,4 +18,7 @@ makedocs(
 
 deploydocs(
     repo   = "github.com/JuliaPlasma/ChargedParticleDynamics.jl",
+    devurl = "latest",
+    devbranch = "main",
+    push_preview = true,
 )

examples/Project.toml

@@ -4,7 +4,13 @@ DecFP = "55939f99-70c6-5e9b-8bb0-5071ed7d61fd"
 ElectromagneticFields = "d6c1ba6f-ee03-53af-b876-68cefeb88ec8"
 GeometricIntegrators = "dcce2d33-59f6-5b8d-9047-0defad88ae06"
 GeometricProblems = "18cb22b4-ad41-5c80-9c5f-710df63fbdc9"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+PGFPlotsX = "8314cec4-20b6-5062-9cdb-752b83310925"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+Plotly = "58dd65bb-95f3-509e-9936-c39a10fdeae7"
 PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PoincareInvariants = "26663084-47d3-540f-bd97-40ca743aafa4"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+SimpleSolvers = "36b790f5-6434-4fbe-b711-1f64a1e2f6a2"
 Weave = "44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9"

src/ChargedParticle3d.jl

@@ -1,5 +1,6 @@
 module ChargedParticle3d
 
+    include("charged_particle_3d/singular_field_canonical.jl")
     include("charged_particle_3d/singular_field.jl")
     include("charged_particle_3d/symmetric_field.jl")
     include("charged_particle_3d/solovev_iter.jl")

src/ChargedParticlePlots.jl

@@ -1,6 +1,6 @@
 module ChargedParticlePlots
 
-    using GeometricIntegrators
+    # using GeometricIntegrators
     using Plots
 
 

src/charged_particle_3d/charged_particle_3d_canonical.jl

@@ -1,12 +1,18 @@
+using GeometricEquations: IODEProblem, LODEProblem, PODEProblem
+using GeometricSolutions: GeometricSolution, DataSeries, TimeSeries
 
-using GeometricIntegrators.Equations
+import GeometricProblems.Diagnostics: compute_invariant, compute_invariant_error
 
 
+const Δt = 0.01
+const tᵢ = 0.0
+const tspan = (tᵢ, 10.0)
+
 ϑ₁(t, q, v) = g₁₁(t,q) * v[1] + A₁(t, q)
 ϑ₂(t, q, v) = g₂₂(t,q) * v[2] + A₂(t, q)
 ϑ₃(t, q, v) = g₃₃(t,q) * v[3] + A₃(t, q)
 
-function ϑ(t, q, v, θ)
+function ϑ(θ, t, q, v)
     θ[1] = ϑ₁(t,q,v)
     θ[2] = ϑ₂(t,q,v)
     θ[3] = ϑ₃(t,q,v)
@@ -17,32 +23,26 @@ v¹(t, q, p) = g¹¹(t, q) * (p[1] - A₁(t, q))
 v²(t, q, p) = g²²(t, q) * (p[2] - A₂(t, q))
 v³(t, q, p) = g³³(t, q) * (p[3] - A₃(t, q))
 
+lagrangian(t,q,v) = (g₁₁(t,q) * v[1]^2 + g₂₂(t,q) * v[2]^2 + g₃₃(t,q) * v[3]^2) / 2 - φ(t,q)
 hamiltonian(t,q,p) = (g₁₁(t,q) * v¹(t,q,p)^2 + g₂₂(t,q) * v²(t,q,p)^2 + g₃₃(t,q) * v³(t,q,p)^2) / 2 + φ(t,q)
 toroidal_momentum(t,q,p) = p[3]
 
 
-function charged_particle_3d_pᵢ(qᵢ, vᵢ, tᵢ=0)
+function charged_particle_3d_pᵢ(tᵢ, qᵢ, vᵢ)
     pᵢ = zero(qᵢ)
-
-    if ndims(qᵢ) == 1
-        ϑ(tᵢ, qᵢ, vᵢ, pᵢ)
-    else
-        for i in 1:size(qᵢ,2)
-            @views ϑ(tᵢ, qᵢ[:,i], vᵢ[:,i], pᵢ[:,i])
-        end
-    end
+    ϑ(pᵢ, tᵢ, qᵢ, vᵢ)
     pᵢ
 end
 
 
-function charged_particle_3d_pode_v(t, q, p, v)
+function charged_particle_3d_pode_v(v, t, q, p, params)
     v[1] = v¹(t,q,p)
     v[2] = v²(t,q,p)
     v[3] = v³(t,q,p)
     nothing
 end
 
-function charged_particle_3d_pode_f(t, q, p, f)
+function charged_particle_3d_pode_f(f, t, q, p, params)
     f[1] = dA₁dx₁(t,q) * v¹(t,q,p) + dA₂dx₁(t,q) * v²(t,q,p) + dA₃dx₁(t,q) * v³(t,q,p) + E₁(t,q) +
            (dg₁₁dx₁(t,q) * v¹(t,q,p)^2 + dg₂₂dx₁(t,q) * v²(t,q,p)^2 + dg₃₃dx₁(t,q) * v³(t,q,p)^2) / 2
     f[2] = dA₁dx₂(t,q) * v¹(t,q,p) + dA₂dx₂(t,q) * v²(t,q,p) + dA₃dx₂(t,q) * v³(t,q,p) + E₂(t,q) +
@@ -52,9 +52,9 @@ function charged_particle_3d_pode_f(t, q, p, f)
     nothing
 end
 
-charged_particle_3d_iode_ϑ(t, q, v, θ) = ϑ(t, q, v, θ)
+charged_particle_3d_iode_ϑ(θ, t, q, v, params) = ϑ(θ, t, q, v)
 
-function charged_particle_3d_iode_f(t, q, v, f)
+function charged_particle_3d_iode_f(f, t, q, v, params)
     f[1] = dA₁dx₁(t,q) * v[1] + dA₂dx₁(t,q) * v[2] + dA₃dx₁(t,q) * v[3] + E₁(t,q) +
            (dg₁₁dx₁(t,q) * v[1]^2 + dg₂₂dx₁(t,q) * v[2]^2 + dg₃₃dx₁(t,q) * v[3]^2) / 2
     f[2] = dA₁dx₂(t,q) * v[1] + dA₂dx₂(t,q) * v[2] + dA₃dx₂(t,q) * v[3] + E₂(t,q) +
@@ -64,27 +64,51 @@ function charged_particle_3d_iode_f(t, q, v, f)
     nothing
 end
 
-function charged_particle_3d_iode_g(t, q, v, f)
-    f[1] = dA₁dx₁(t,q) * v[1] + dA₂dx₁(t,q) * v[2] + dA₃dx₁(t,q) * v[3] +
-           (dg₁₁dx₁(t,q) * v[1]^2 + dg₂₂dx₁(t,q) * v[2]^2 + dg₃₃dx₁(t,q) * v[3]^2) / 2
-    f[2] = dA₁dx₂(t,q) * v[1] + dA₂dx₂(t,q) * v[2] + dA₃dx₂(t,q) * v[3] +
-           (dg₁₁dx₂(t,q) * v[1]^2 + dg₂₂dx₂(t,q) * v[2]^2 + dg₃₃dx₂(t,q) * v[3]^2) / 2
-    f[3] = dA₁dx₃(t,q) * v[1] + dA₂dx₃(t,q) * v[2] + dA₃dx₃(t,q) * v[3] +
-           (dg₁₁dx₃(t,q) * v[1]^2 + dg₂₂dx₃(t,q) * v[2]^2 + dg₃₃dx₃(t,q) * v[3]^2) / 2
+function charged_particle_3d_iode_g(g, t, q, v, λ, params)
+    g[1] = dA₁dx₁(t,q) * λ[1] + dA₂dx₁(t,q) * λ[2] + dA₃dx₁(t,q) * λ[3] +
+           (dg₁₁dx₁(t,q) * λ[1]^2 + dg₂₂dx₁(t,q) * λ[2]^2 + dg₃₃dx₁(t,q) * λ[3]^2) / 2
+    g[2] = dA₁dx₂(t,q) * λ[1] + dA₂dx₂(t,q) * λ[2] + dA₃dx₂(t,q) * λ[3] +
+           (dg₁₁dx₂(t,q) * λ[1]^2 + dg₂₂dx₂(t,q) * λ[2]^2 + dg₃₃dx₂(t,q) * λ[3]^2) / 2
+    g[3] = dA₁dx₃(t,q) * λ[1] + dA₂dx₃(t,q) * λ[2] + dA₃dx₃(t,q) * λ[3] +
+           (dg₁₁dx₃(t,q) * λ[1]^2 + dg₂₂dx₃(t,q) * λ[2]^2 + dg₃₃dx₃(t,q) * λ[3]^2) / 2
     nothing
 end
 
 
 # function charged_particle_3d_pode(q₀=qᵢ, v₀=vᵢ)
 #     PODE(charged_particle_3d_pode_v, charged_particle_3d_pode_f, q₀, charged_particle_3d_pᵢ(q₀, v₀))
 # end
-function charged_particle_3d_pode(q₀=qᵢ, p₀=pᵢ)
-    PODE(charged_particle_3d_pode_v, charged_particle_3d_pode_f, q₀, p₀)
+function charged_particle_3d_pode(q₀=qᵢ, p₀=pᵢ; tspan=tspan, tstep=Δt)
+    PODEProblem(
+        charged_particle_3d_pode_v,
+        charged_particle_3d_pode_f,
+        tspan, tstep, q₀, p₀;
+        invariants = (h = hamiltonian,))
 end
 
 
-function charged_particle_3d_iode(q₀=qᵢ, p₀=pᵢ)
-    IODE(charged_particle_3d_iode_ϑ, charged_particle_3d_iode_f, charged_particle_3d_iode_g,
-            q₀, p₀;
-            h=hamiltonian, v = (t, q, v) -> nothing)
+function charged_particle_3d_iode(q₀=qᵢ, p₀=pᵢ; tspan=tspan, tstep=Δt)
+    IODEProblem(
+        charged_particle_3d_iode_ϑ,
+        charged_particle_3d_iode_f,
+        charged_particle_3d_iode_g,
+        tspan, tstep, q₀, p₀;
+        invariants = (h = hamiltonian,))
 end
+
+function charged_particle_3d_lode(q₀=qᵢ, p₀=pᵢ; tspan=tspan, tstep=Δt)
+    LODEProblem(
+        charged_particle_3d_iode_ϑ,
+        charged_particle_3d_iode_f,
+        charged_particle_3d_iode_g,
+        ω, lagrangian,
+        tspan, tstep, q₀, p₀;
+        invariants = (h = hamiltonian,))
+end
+
+
+compute_energy(t::TimeSeries, q::DataSeries, p::DataSeries) = compute_invariant(t, q, p, hamiltonian)
+compute_energy(sol::GeometricSolution) = compute_energy(sol.t, sol.q, sol.p)
+
+compute_energy_error(t::TimeSeries, q::DataSeries, p::DataSeries) = compute_invariant_error(t, q, p, hamiltonian)
+compute_energy_error(sol::GeometricSolution) = compute_energy_error(sol.t, sol.q, sol.p)