Use DispersionRelations package

Project.toml

@@ -4,6 +4,7 @@ authors = "JuliaVlasov organization"
 version = "0.1.0"
 
 [deps]
+DispersionRelations = "78f632dc-6214-4aca-89df-9f86f5ac7c9e"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

`

@@ -0,0 +1,121 @@
+# Quickstart - Landan damping simulation
+
+```@example 1
+using VlasovSolvers, Plots
+
+dev = CPU()
+```
+
+For now, only the `CPU` device is avalaible but we hope it will be possible 
+to run simulation also on `GPU`. 
+
+## Vlasov-Poisson 
+
+We consider the dimensionless Vlasov-Poisson equation for one species
+with a neutralizing background.
+
+```math
+\frac{\partial f}{\partial t}+ v\cdot \nabla_x f + E(t,x) \cdot \nabla_v f = 0,
+```
+
+```math
+- \Delta \phi = 1 - \rho, E = - \nabla \phi
+```
+
+```math
+\rho(t,x)  =  \int f(t,x,v)dv.
+```
+
+- [Vlasov Equation - Wikipedia](https://en.wikipedia.org/wiki/Vlasov_equation)
+
+
+## Landau damping
+
+- [Landau damping - Wikipedia](https://en.wikipedia.org/wiki/Landau_damping)
+
+We initialize the numerical domain by defining two grids, one for space, one for velocity:
+
+```@docs
+VlasovSolvers.OneDGrid
+```
+
+```@example 1
+nx, nv = 64, 64              # grid resolution
+xmin, xmax = 0, 4π           # X Domain length (m)
+vmin, vmax = -6, 6           # V Domain length (m)
+xgrid = OneDGrid(dev, nx, xmin, xmax)
+vgrid = OneDGrid(dev, nv, vmin, vmax)
+```
+
+We instantiate a distribution function type to store its values and the grids.
+
+```@docs
+VlasovSolvers.DistributionFunction
+```
+
+```@example 1
+f = DistributionFunction( xgrid, vgrid )
+```
+
+We initialize the distribution function for the Landau damping test case:
+
+```@docs
+VlasovSolvers.landau!
+```
+
+```@example 1
+α  = 0.001                   # Perturbation amplitude
+kx = 0.5                     # Wave number of perturbation
+landau!(f, α, kx)
+```
+
+To run the simulation we need to create a [`VlasovProblem`](@ref) and choose a 
+method to compute advections. We use the backward semi-lagrangian method with 
+5th order b-splines for interpolation [`BSLSpline`](@ref).
+
+```@docs
+VlasovSolvers.VlasovProblem
+```
+
+```@docs
+VlasovSolvers.BSLSpline
+```
+
+```@example 1
+prob = VlasovProblem(f, BSLSpline(5), dev)
+```
+
+```@example 1
+
+stepper = StrangSplitting()  # timestepper
+dt = 0.1                     # timestep
+nsteps = 1000                # total number of time-steps
+
+```
+
+
+```@example 1
+
+sol = solve(prob, stepper, dt, nsteps )
+
+t = sol.times
+
+plot(sol; yscale= :log10, label = "E")
+line, gamma = fit_complex_frequency(t, sol.energy)
+plot!(t, line; label = "growth = $(imag(gamma))")
+```
+
+## Simulation with Lagrange interpolation
+
+```@example 1
+landau!(f, α, kx)
+
+prob = VlasovProblem(f, BSLLagrange(9), dev)
+
+sol = solve(prob, stepper, dt, nsteps )
+
+t = sol.times
+
+plot(sol; label = "E")
+plot!(t, -0.1533*t.-5.50; label="-0.1533t.-5.5")
+```

docs/src/bump_on_tail.md

@@ -36,5 +36,5 @@ prob = VlasovProblem(f, BSLSpline(5), dev)
 
 sol = solve(prob, stepper, dt, nsteps )
 
-plot(sol, label=L"\frac{1}{2} \log(\int e^2dx)")
+plot(sol, yscale = :log10, label=L"\sqrt{\int e^2dx}")
 ```

docs/src/problem.md

@@ -100,8 +100,9 @@ sol = solve(prob, stepper, dt, nsteps )
 
 t = sol.times
 
-plot(sol; label = "E")
-plot!(t, -0.1533*t.-5.50; label="-0.1533t.-5.5")
+plot(sol; yscale= :log10, label = "E")
+line, gamma = fit_complex_frequency(t, sol.energy)
+plot!(t, line; label = "growth = $(imag(gamma))")
 ```
 
 ## Simulation with Lagrange interpolation
@@ -115,6 +116,7 @@ sol = solve(prob, stepper, dt, nsteps )
 
 t = sol.times
 
-plot(sol; label = "E")
-plot!(t, -0.1533*t.-5.50; label="-0.1533t.-5.5")
+plot(sol; yscale= :log10, label = "E")
+line, gamma = fit_complex_frequency(t, sol.energy)
+plot!(t, line; label = "growth = $(imag(gamma))")
 ```

docs/src/vlasov-ampere.md

@@ -62,6 +62,7 @@ dt = 0.1
 
 sol = solve(prob, stepper, dt, nsteps )
 
-plot(sol.times, -0.1533*sol.times .- 5.48)
-plot!(sol; label="ampere" )
+plot(sol; yscale= :log10, label = "E")
+line, gamma = fit_complex_frequency(sol.times, sol.energy)
+plot!(sol.times, line; label = "growth = $(imag(gamma))")
 ```

src/VlasovSolvers.jl

@@ -1,10 +1,12 @@
 module VlasovSolvers
 
+  using DispersionRelations
   using DocStringExtensions
   using FFTW
   using LinearAlgebra
   using Statistics
 
+  export fit_complex_frequency
   export solve
 
   include("devices.jl")

src/advection.jl

@@ -144,7 +144,7 @@ function advection_v!( f :: DistributionFunction, dt, method::BSLSpline)
     dx = f.xgrid.step
     nx = f.xgrid.len
     e  = compute_e(f)
-    nrj = 0.5*log(sum(e.*e)*dx)
+    nrj = sqrt(sum(e.*e)*dx)
     fᵗ = collect(transpose(f.values))
     advection!(fᵗ, f.vgrid, e, dt; method.p)
     f.values .= transpose(fᵗ)
@@ -174,7 +174,7 @@ function advection_v!( f :: DistributionFunction, dt, method::BSLLagrange)
     dx = f.xgrid.step
     nx = f.xgrid.len
     e  = compute_e(f)
-    nrj = 0.5*log(sum(e.*e)*dx)
+    nrj = sqrt(sum(e.*e)*dx)
 
     nv = f.vgrid.len
     fi = zeros(nv)

src/fourier.jl

@@ -15,12 +15,12 @@ struct Fourier <: AbstractMethod
         nx  = meshx.len
         dx  = meshx.step
         Lx  = meshx.stop - meshx.start
-        kx  = 2π/Lx .* [0:nx÷2-1;-nx÷2:-1]
+        kx  = 2π/Lx .* fftfreq(nx, nx)
 
         nv  = meshv.len
         dv  = meshv.step
         Lv  = meshv.stop - meshv.start
-        kv  = 2π/Lv .* [0:nv÷2-1;-nv÷2:-1]
+        kv  = 2π/Lv .* fftfreq(nx, nx)
 
         new( kx, kv)
 

src/problems.jl

@@ -80,15 +80,15 @@ function solve( problem::VlasovProblem{Fourier}, stepper::StrangSplitting, dt, n
     sol = VlasovSolution1D1V() 
     time = 0.0
     push!(sol.times, time)
-    energy = log(sqrt((sum(e.^2)) * dx))
+    energy = sqrt((sum(e.^2)) * dx)
     push!(sol.energy, energy)
 
     for it in 1:nsteps
 
         advection_v!(fᵀ, problem.method, e, 0.5dt)
         transpose!(f,fᵀ)
         advection_x!( f, problem.method, e, v, dt)
-        energy = log(sqrt((sum(e.^2)) * dx))
+        energy = sqrt((sum(e.^2)) * dx)
         time += dt
         push!(sol.times, time)
         push!(sol.energy, energy)