More fixes and updates for new versions of PoissonSolvers, GeometricEquations, and GeometricIntegrators.

Author: michakraus

Project.toml

@@ -31,15 +31,15 @@ TerminalLoggers = "5d786b92-1e48-4d6f-9151-6b4477ca9bed"
 AdaptiveRejectionSampling = "0.1"
 BSplineKit = "0.14, 0.15, 0.16, 0.17"
 DifferentialEquations = "7"
-GeometricEquations = "0.14"
+GeometricEquations = "0.16"
 GeometricIntegrators = "0.13"
 HDF5 = "0.15, 0.16, 0.17"
 LaTeXStrings = "1.2, 1.3"
 OffsetArrays = "1"
 Parameters = "0.12"
 ParticleMethods = "0.1"
 Plots = "1"
-PoissonSolvers = "0.3"
+PoissonSolvers = "0.3.5"
 ProgressMeter = "1"
 Sobol = "1.5"
 SpecialFunctions = "2"

scripts/lenard_bernstein.jl

@@ -68,4 +68,3 @@ anim = @animate for i in 1:step:length(sol)
 end
 println("saving animation")
 gif(anim,"lenard_bernstein.gif",fps=2)
-

scripts/lenard_bernstein_conservative.jl

@@ -79,4 +79,4 @@ anim = @animate for n in 1:step:size(z,2)
 end
 
 # save animation to file
-gif(anim, "lenard_bernstein_conservative_anim.gif", fps=10)
+gif(anim, "lenard_bernstein_conservative.gif", fps=10)

scripts/vlasov_poisson.jl

@@ -1,4 +1,3 @@
-
 # import libraries
 using PoissonSolvers
 using VlasovMethods
@@ -15,7 +14,7 @@ domain = (0.0, 1.0)
 h5file = "vlasov_poisson.hdf5"
 
 # create and initialize particle distribution function
-dist = initialize!(ParticleDistribution(1, 1, npart), Normal())
+dist = initialize!(ParticleDistribution(1, 1, npart), NormalDistribution())
 # dist = initialize!(ParticleDistribution(1, 1, npart), BumpOnTail())
 
 # create electrostatic potential
@@ -70,5 +69,4 @@ anim = @animate for n in axes(z,3)
 end
 
 # save animation to file
-gif(anim, "vlasov_poisson_anim.gif", fps=10)
-
+gif(anim, "vlasov_poisson.gif", fps=10)

src/methods/splitting.jl

@@ -33,12 +33,12 @@ function run!(method::SplittingMethod, h5file)
     h5z = create_dataset(h5, "z", eltype(z₀), ((nd, np, ntime(method.equation)+1), (nd, np, -1)), chunk=(nd,np,1))
     copy_to_hdf5(h5z, z₀, 0)
 
-    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)
             copy_to_hdf5(h5z, method.integrator.solstep.q, n)
         end
     finally

src/models/lenard_bernstein.jl

@@ -76,8 +76,8 @@ function GeometricIntegrator(model::LenardBernstein{1,1}, tspan::Tuple, tstep::R
             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/vlasov_poisson.jl

@@ -1,19 +1,17 @@
 
-struct VlasovPoisson{XD, VD, DT <: DistributionFunction{XD,VD}, PT <: Potential, RT <: AbstractVector} <: VlasovModel
+struct VlasovPoisson{XD, VD, DT <: DistributionFunction{XD,VD}, PT <: Potential} <: VlasovModel
     distribution::DT
     potential::PT
-    rhs::RT
 
     function VlasovPoisson(dist::DistributionFunction{XD,VD}, potential) where {XD,VD}
-        rhs = zero(potential.coefficients)
-        new{XD, VD, typeof(dist), typeof(potential), typeof(rhs)}(dist, potential, rhs)
+        new{XD, VD, typeof(dist), typeof(potential)}(dist, potential)
     end
 end
 
 
 function update_potential!(model::VlasovPoisson)
-    projection!(model.rhs, model.potential, model.distribution)
-    PoissonSolvers.update!(model.potential, model.rhs)
+    projection!(model.potential, model.distribution)
+    PoissonSolvers.update!(model.potential)
 end
 
 
@@ -42,8 +40,8 @@ function v_advection!(ż, t, z, params)
     end
 end
 
-# Vector field for Lorentz force
-function v_lorentz_force!(ż, t, z, params)
+# Vector field for acceleration
+function v_acceleration!(ż, t, z, params)
     update_potential!(params.model)
     for i in axes(ż, 2)
         ż[1,i] = 0
@@ -60,7 +58,7 @@ function s_advection!(z, t, z̄, t̄, params)
 end
 
 # Solution for Lorentz force
-function s_lorentz_force!(z, t, z̄, t̄, params)
+function s_acceleration!(z, t, z̄, t̄, params)
     update_potential!(params.model)
     for i in axes(z, 2)
         z[1,i] = z̄[1,i]
@@ -72,20 +70,19 @@ end
 # 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)
 
     # create geometric problem
     equ = GeometricEquations.SODEProblem(
-            (v_advection!, v_lorentz_force!),
-            (s_advection!, s_lorentz_force!),
+            (v_advection!, v_acceleration!),
+            (s_advection!, s_acceleration!),
             tspan, tstep, copy(model.distribution.particles.z);
             parameters = params)
 
     # create integrator
-    int = Integrators.Integrator(equ, Integrators.Strang())
+    int = GeometricIntegrators.GeometricIntegrator(equ, GeometricIntegrators.Strang())
 
     # put together splitting method
     SplittingMethod(model, equ, int)

src/projections/potential.jl

@@ -1,6 +1,6 @@
 
-function projection!(out, potential::PoissonSolvers.Potential{<:PeriodicBSplineBasis}, distribution::ParticleDistribution)
-    b = Splines.PeriodicVector(out)
+function projection!(potential::PoissonSolvers.Potential{<:PeriodicBSplineBasis}, distribution::ParticleDistribution)
+    b = Splines.PeriodicVector(potential.rhs)
     b .= 0
 
     basis = potential.basis
@@ -18,5 +18,5 @@ function projection!(out, potential::PoissonSolvers.Potential{<:PeriodicBSplineB
         end
     end
 
-    return out
+    return potential
 end

test/projections_tests.jl

@@ -22,9 +22,9 @@ using VlasovMethods: projection!
     dist.particles.w .= (ones(npart) ./ npart)'
 
     rhs = zero(potential.coefficients)
-    projection!(rhs, potential, dist)
+    projection!(potential, dist)
 
-    ρ = Spline(potential.basis, potential.solver.Mfac \ rhs)
+    ρ = Spline(potential.basis, potential.solver.Mfac \ potential.rhs)
 
     x = domain[begin]:0.1:domain[end]