From 12ec4f6616f2d3a863e68d1c7daf23768d3f212c Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 13 May 2024 14:49:38 +0200 Subject: [PATCH] Some cleanup, refactoring, and updating for compatibility with the latest GeometricIntegrators. --- scripts/lenard_bernstein.jl | 6 ++--- ...in.jl => lenard_bernstein_conservative.jl} | 6 ++--- src/VlasovMethods.jl | 11 +++++----- src/methods/geometric_integrator.jl | 4 ++-- src/methods/splitting.jl | 2 +- src/models/collision_operator.jl | 2 ++ src/models/lenard_bernstein.jl | 2 +- ...lb.jl => lenard_bernstein_conservative.jl} | 6 ++--- src/models/model.jl | 2 +- src/models/vlasov_model.jl | 2 ++ src/models/vlasov_poisson.jl | 22 ++++++++++++------- 11 files changed, 38 insertions(+), 27 deletions(-) rename scripts/{cons_lenard_bernstein.jl => lenard_bernstein_conservative.jl} (93%) create mode 100644 src/models/collision_operator.jl rename src/models/{conservative_lb.jl => lenard_bernstein_conservative.jl} (92%) create mode 100644 src/models/vlasov_model.jl diff --git a/scripts/lenard_bernstein.jl b/scripts/lenard_bernstein.jl index 100e0a3..84c9003 100644 --- a/scripts/lenard_bernstein.jl +++ b/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)) diff --git a/scripts/cons_lenard_bernstein.jl b/scripts/lenard_bernstein_conservative.jl similarity index 93% rename from scripts/cons_lenard_bernstein.jl rename to scripts/lenard_bernstein_conservative.jl index e19c2a6..1687128 100644 --- a/scripts/cons_lenard_bernstein.jl +++ b/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) diff --git a/src/VlasovMethods.jl b/src/VlasovMethods.jl index 8503677..4b4afc7 100644 --- a/src/VlasovMethods.jl +++ b/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 diff --git a/src/methods/geometric_integrator.jl b/src/methods/geometric_integrator.jl index 8946bca..d520702 100644 --- a/src/methods/geometric_integrator.jl +++ b/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 diff --git a/src/methods/splitting.jl b/src/methods/splitting.jl index 97535d2..1ad351e 100644 --- a/src/methods/splitting.jl +++ b/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 diff --git a/src/models/collision_operator.jl b/src/models/collision_operator.jl new file mode 100644 index 0000000..4c905b9 --- /dev/null +++ b/src/models/collision_operator.jl @@ -0,0 +1,2 @@ + +abstract type CollisionOperator <: AbstractVlasovModel end diff --git a/src/models/lenard_bernstein.jl b/src/models/lenard_bernstein.jl index a6fa008..a1412c3 100644 --- a/src/models/lenard_bernstein.jl +++ b/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 diff --git a/src/models/conservative_lb.jl b/src/models/lenard_bernstein_conservative.jl similarity index 92% rename from src/models/conservative_lb.jl rename to src/models/lenard_bernstein_conservative.jl index 0c578e1..5970e6d 100644 --- a/src/models/conservative_lb.jl +++ b/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) diff --git a/src/models/model.jl b/src/models/model.jl index 146e546..26543fa 100644 --- a/src/models/model.jl +++ b/src/models/model.jl @@ -1,2 +1,2 @@ -abstract type VlasovModel end +abstract type AbstractVlasovModel end diff --git a/src/models/vlasov_model.jl b/src/models/vlasov_model.jl new file mode 100644 index 0000000..98b26c3 --- /dev/null +++ b/src/models/vlasov_model.jl @@ -0,0 +1,2 @@ + +abstract type VlasovModel <: AbstractVlasovModel end diff --git a/src/models/vlasov_poisson.jl b/src/models/vlasov_poisson.jl index e928477..414225b 100644 --- a/src/models/vlasov_poisson.jl +++ b/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,7 +30,11 @@ 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] @@ -43,6 +42,7 @@ function v_advection!(ż, t, z, params) end end +# Vector field for Lorentz force function v_lorentz_force!(ż, t, z, params) update_potential!(params.model) for i in axes(ż, 2) @@ -51,6 +51,7 @@ 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] @@ -58,6 +59,7 @@ function s_advection!(z, t, z̄, t̄, params) 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)