Some cleanup, refactoring, and updating for compatibility with the la…

…test GeometricIntegrators.

scripts/lenard_bernstein.jl

@@ -46,7 +46,7 @@ params = (ν = model.ν, idist = model.dist, fdist = model.ent.dist, model = mod
 
 
 f = projection(sol[:,1], dist, sdist)
-v = LB_rhs(collect(vgrid), params, f)
+v = VlasovMethods.LB_rhs(collect(vgrid), params, f)
 mom = [mapreduce(p -> p[1], +, sol[:,i]) for i in 1:step:length(sol)]
 enr = [mapreduce(p -> p[1].^2, +, sol[:,i]) for i in 1:step:length(sol)]
 
@@ -56,8 +56,8 @@ anim = @animate for i in 1:step:length(sol)
     # compute quantities for plotting
     f = projection(sol[:,i], dist, sdist)
     df = Derivative(1) * f
-    v = LB_rhs(collect(vgrid), params, f)
-    # v = CLB_rhs(collect(vgrid), params, f)
+    v = VlasovMethods.LB_rhs(collect(vgrid), params, f)
+    # v = VlasovMethods.CLB_rhs(collect(vgrid), params, f)
 
     plot(xlims = [-8, +8], ylims = [-0.5, +0.5], size=(1200,800))
 

scripts/cons_lenard_bernstein.jl → scripts/lenard_bernstein_conservative.jl

@@ -3,7 +3,7 @@ using BSplineKit
 using VlasovMethods
 
 # output file
-h5file = "lenard_bernstein.hdf5"
+h5file = "lenard_bernstein_conservative.hdf5"
 
 # params
 # parameters
@@ -67,7 +67,7 @@ anim = @animate for n in 1:step:size(z,2)
     f = projection(z[:,n], dist, sdist)
     df = Derivative(1) * f
     # v = LB_rhs(collect(vgrid), params, f)
-    v = CLB_rhs(collect(vgrid), params, f)
+    v = VlasovMethods.CLB_rhs(collect(vgrid), params, f)
 
     plot(xlims = [-8, +8], ylims = [-0.5, +0.5], size=(1200,800))
 
@@ -79,4 +79,4 @@ anim = @animate for n in 1:step:size(z,2)
 end
 
 # save animation to file
-gif(anim, "lenard_bernstein_anim.gif", fps=10)
+gif(anim, "lenard_bernstein_conservative_anim.gif", fps=10)

src/VlasovMethods.jl

@@ -18,7 +18,7 @@ using LinearAlgebra
 
 import GeometricEquations
 import GeometricEquations: ntime
-import GeometricIntegrators.Integrators
+import GeometricIntegrators
 import DifferentialEquations
 
 
@@ -72,15 +72,16 @@ export GeometricIntegrator
 
 # Vlasov models
 
-include("models/vlasov_poisson.jl")
+include("models/collision_operator.jl")
 include("models/lenard_bernstein.jl")
-include("models/conservative_lb.jl")
+include("models/lenard_bernstein_conservative.jl")
+
+include("models/vlasov_model.jl")
+include("models/vlasov_poisson.jl")
 
 export VlasovPoisson
 export LenardBernstein
 export ConservativeLenardBernstein
-export LB_rhs
-export CLB_rhs
 
 # Example Problems
 

src/methods/geometric_integrator.jl

@@ -24,12 +24,12 @@ function run!(method::GeometricIntegrator, h5file)
     h5z[:,1] = z₀
     h5t[1] = t₀
 
-    Integrators.initialize!(method.integrator)
+    GeometricIntegrators.Integrators.initialize!(method.integrator)
 
     # loop over time steps showing progress bar
     try
         @showprogress 5 for n in 1:ntime(method.equation)
-            Integrators.integrate!(method.integrator)
+            GeometricIntegrators.integrate!(method.integrator)
             h5z[:,n+1] = method.integrator.solstep.q
             h5t[n+1] = method.integrator.solstep.t
         end

src/methods/splitting.jl

@@ -4,7 +4,7 @@ struct SplittingMethod{MT, ET, IT} <: ParticleMethod
     equation::ET
     integrator::IT
 
-    function SplittingMethod(model::MT, equation::ET, integrator::IT) where {MT <: VlasovModel, ET <: GeometricEquations.GeometricProblem, IT}
+    function SplittingMethod(model::MT, equation::ET, integrator::IT) where {MT <: AbstractVlasovModel, ET <: GeometricEquations.GeometricProblem, IT}
         new{MT,ET,IT}(model, equation, integrator)
     end
 end

src/models/collision_operator.jl

@@ -0,0 +1,2 @@
+
+abstract type CollisionOperator <: AbstractVlasovModel end

src/models/lenard_bernstein.jl

@@ -1,4 +1,4 @@
-struct LenardBernstein{XD, VD, DT <: DistributionFunction{XD,VD}, ET <: Entropy, T} <: VlasovModel
+struct LenardBernstein{XD, VD, DT <: DistributionFunction{XD,VD}, ET <: Entropy, T} <: CollisionOperator
     dist::DT    # distribution function
     ent::ET     # entropy 
     ν::T        # collision frequency 

src/models/conservative_lb.jl → src/models/lenard_bernstein_conservative.jl

@@ -1,4 +1,4 @@
-struct ConservativeLenardBernstein{XD, VD, DT <: DistributionFunction{XD,VD}, ET <: Entropy, T} <: VlasovModel
+struct ConservativeLenardBernstein{XD, VD, DT <: DistributionFunction{XD,VD}, ET <: Entropy, T} <: CollisionOperator
     dist::DT    # distribution function
     ent::ET     # entropy 
     ν::T        # collision frequency 
@@ -96,8 +96,8 @@ function GeometricIntegrator(model::ConservativeLenardBernstein{1,1}, tspan::Tup
             parameters = params)
 
     # create integrator
-    int = Integrators.Integrator(equ, Integrators.RK438())
-    # int = Integrators.Integrator(equ, Integrators.CrankNicolson())
+    int = GeometricIntegrators.GeometricIntegrator(equ, GeometricIntegrators.RK438())
+    # int = GeometricIntegrators.GeometricIntegrator(equ, GeometricIntegrators.CrankNicolson())
 
     # put together splitting method
     GeometricIntegrator(model, equ, int)

src/models/model.jl

@@ -1,2 +1,2 @@
 
-abstract type VlasovModel end
+abstract type AbstractVlasovModel end

src/models/vlasov_model.jl

@@ -0,0 +1,2 @@
+
+abstract type VlasovModel <: AbstractVlasovModel end

src/models/vlasov_poisson.jl

@@ -11,17 +11,12 @@ struct VlasovPoisson{XD, VD, DT <: DistributionFunction{XD,VD}, PT <: Potential,
 end
 
 
-
 function update_potential!(model::VlasovPoisson)
     projection!(model.rhs, model.potential, model.distribution)
     PoissonSolvers.update!(model.potential, model.rhs)
 end
 
 
-
-
-
-
 ####################################################
 # Define Splitting Method for Vlasov-Poisson Model #
 ####################################################
@@ -35,14 +30,19 @@ function lorentz_force!(ż, t, z, params)
     end
 end
 
-# splitting fields
+###########################################################
+# Vlasov-Poisson 1D1V splitting fields for particles      #
+###########################################################
+
+# Vector field for advection
 function v_advection!(ż, t, z, params)
     for i in axes(ż, 2)
         ż[1,i] = z[2,i]
         ż[2,i] = 0
     end
 end
 
+# Vector field for Lorentz force
 function v_lorentz_force!(ż, t, z, params)
     update_potential!(params.model)
     for i in axes(ż, 2)
@@ -51,13 +51,15 @@ function v_lorentz_force!(ż, t, z, params)
     end
 end
 
+# Solution for advection
 function s_advection!(z, t, z̄, t̄, params)
     for i in axes(z, 2)
         z[1,i] = z̄[1,i] + (t-t̄) * z̄[2,i]
         z[2,i] = z̄[2,i]
     end
 end
 
+# Solution for Lorentz force
 function s_lorentz_force!(z, t, z̄, t̄, params)
     update_potential!(params.model)
     for i in axes(z, 2)
@@ -66,8 +68,12 @@ function s_lorentz_force!(z, t, z̄, t̄, params)
     end
 end
 
-
-function SplittingMethod(model::VlasovPoisson{1,1}, tspan::Tuple, tstep::Real)
+# Constructor for a splitting method from GeometricIntegrators
+# The problem is setup such that one solution step pushes all particles.
+# While this allows for a simple implementation, it is not well-suited
+# for parallelisation.
+# Have a look at the CollisionalVlasovPoisson model for an alternative approach.
+function SplittingMethod(model::VlasovPoisson{1, 1, <: ParticleDistribution}, tspan::Tuple, tstep::Real)
     # collect parameters
     params = (ϕ = model.potential, model = model)