diff --git a/Project.toml b/Project.toml index 349a50e..06dd673 100644 --- a/Project.toml +++ b/Project.toml @@ -31,7 +31,7 @@ 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" @@ -39,7 +39,7 @@ OffsetArrays = "1" Parameters = "0.12" ParticleMethods = "0.1" Plots = "1" -PoissonSolvers = "0.3" +PoissonSolvers = "0.3.5" ProgressMeter = "1" Sobol = "1.5" SpecialFunctions = "2" diff --git a/scripts/lenard_bernstein.jl b/scripts/lenard_bernstein.jl index 84c9003..36fe6ef 100644 --- a/scripts/lenard_bernstein.jl +++ b/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) - diff --git a/scripts/lenard_bernstein_conservative.jl b/scripts/lenard_bernstein_conservative.jl index 1687128..ca460d3 100644 --- a/scripts/lenard_bernstein_conservative.jl +++ b/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) diff --git a/scripts/vlasov_poisson.jl b/scripts/vlasov_poisson.jl index b8b6fdc..a1b4095 100644 --- a/scripts/vlasov_poisson.jl +++ b/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) diff --git a/src/methods/splitting.jl b/src/methods/splitting.jl index 1ad351e..304d939 100644 --- a/src/methods/splitting.jl +++ b/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 diff --git a/src/models/lenard_bernstein.jl b/src/models/lenard_bernstein.jl index a1412c3..bb17b75 100644 --- a/src/models/lenard_bernstein.jl +++ b/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) diff --git a/src/models/vlasov_poisson.jl b/src/models/vlasov_poisson.jl index 414225b..8f02fde 100644 --- a/src/models/vlasov_poisson.jl +++ b/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) diff --git a/src/projections/potential.jl b/src/projections/potential.jl index e7c7899..319e8f2 100644 --- a/src/projections/potential.jl +++ b/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 diff --git a/test/projections_tests.jl b/test/projections_tests.jl index 28122fd..0079f88 100644 --- a/test/projections_tests.jl +++ b/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]