GitHub

Quickstart - Landan damping simulation

using VlasovSolvers, Plots

dev = CPU()
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.

\[\frac{\partial f}{\partial t}+ v\cdot \nabla_x f + E(t,x) \cdot \nabla_v f = 0,\]

\[- \Delta \phi = 1 - \rho, E = - \nabla \phi\]

\[\rho(t,x) = \int f(t,x,v)dv.\]

Landau damping

We initialize the numerical domain by defining two grids, one for space, one for velocity:

VlasovSolvers.OneDGrid — Type
struct OneDGrid

  • dev::VlasovSolvers.AbstractDevice

  • len::Int64

  • start::Real

  • stop::Real

  • points::LinRange

  • step::Real

source
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)
OneDGrid(CPU(), 64, -6, 6, LinRange{Float64}(-6.0, 5.8125, 64), 0.1875)

We instantiate a distribution function type to store its values and the grids.

f = DistributionFunction( xgrid, vgrid )
DistributionFunction(OneDGrid(CPU(), 64, 0, 12.566370614359172, LinRange{Float64}(0.0, 12.370021073509811, 64), 0.19634954084936207), OneDGrid(CPU(), 64, -6, 6, LinRange{Float64}(-6.0, 5.8125, 64), 0.1875), AbstractFloat[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])

We initialize the distribution function for the Landau damping test case:

VlasovSolvers.landau! — Function
landau!(f, α, kx)

Initialize the distribution function for the Landau damping test case

\[f(x,v,t=0) = \frac{1}{\sqrt{2\pi}} ( 1 + \cos{kx \cdot x} ) \exp (-\frac{v^2}{2})\]

source
α  = 0.001                   # Perturbation amplitude
kx = 0.5                     # Wave number of perturbation
landau!(f, α, kx)
64×64 Matrix{AbstractFloat}:
 6.08196e-9  1.84073e-8  5.37861e-8  …  1.51733e-7  5.37861e-8  1.84073e-8
 6.08193e-9  1.84072e-8  5.37858e-8     1.51732e-7  5.37858e-8  1.84072e-8
 6.08184e-9  1.8407e-8   5.3785e-8      1.5173e-7   5.3785e-8   1.8407e-8
 6.0817e-9   1.84065e-8  5.37837e-8     1.51727e-7  5.37837e-8  1.84065e-8
 6.0815e-9   1.84059e-8  5.3782e-8      1.51722e-7  5.3782e-8   1.84059e-8
 6.08124e-9  1.84052e-8  5.37797e-8  …  1.51715e-7  5.37797e-8  1.84052e-8
 6.08093e-9  1.84042e-8  5.3777e-8      1.51708e-7  5.3777e-8   1.84042e-8
 6.08058e-9  1.84032e-8  5.37739e-8     1.51699e-7  5.37739e-8  1.84032e-8
 6.08018e-9  1.84019e-8  5.37703e-8     1.51689e-7  5.37703e-8  1.84019e-8
 6.07974e-9  1.84006e-8  5.37664e-8     1.51678e-7  5.37664e-8  1.84006e-8
 ⋮                                   ⋱                          
 6.07974e-9  1.84006e-8  5.37664e-8  …  1.51678e-7  5.37664e-8  1.84006e-8
 6.08018e-9  1.84019e-8  5.37703e-8     1.51689e-7  5.37703e-8  1.84019e-8
 6.08058e-9  1.84032e-8  5.37739e-8     1.51699e-7  5.37739e-8  1.84032e-8
 6.08093e-9  1.84042e-8  5.3777e-8      1.51708e-7  5.3777e-8   1.84042e-8
 6.08124e-9  1.84052e-8  5.37797e-8     1.51715e-7  5.37797e-8  1.84052e-8
 6.0815e-9   1.84059e-8  5.3782e-8   …  1.51722e-7  5.3782e-8   1.84059e-8
 6.0817e-9   1.84065e-8  5.37837e-8     1.51727e-7  5.37837e-8  1.84065e-8
 6.08184e-9  1.8407e-8   5.3785e-8      1.5173e-7   5.3785e-8   1.8407e-8
 6.08193e-9  1.84072e-8  5.37858e-8     1.51732e-7  5.37858e-8  1.84072e-8

To run the simulation we need to create a VlasovProblem and choose a method to compute advections. We use the backward semi-lagrangian method with 5th order b-splines for interpolation BSLSpline.

VlasovSolvers.VlasovProblem — Type
struct VlasovProblem{Method<:VlasovSolvers.AbstractMethod} <: VlasovSolvers.AbstractProblem

  • f::DistributionFunction

  • method::VlasovSolvers.AbstractMethod

  • dev::VlasovSolvers.AbstractDevice

source
prob = VlasovProblem(f, BSLSpline(5), dev)
VlasovProblem{BSLSpline}(DistributionFunction(OneDGrid(CPU(), 64, 0, 12.566370614359172, LinRange{Float64}(0.0, 12.370021073509811, 64), 0.19634954084936207), OneDGrid(CPU(), 64, -6, 6, LinRange{Float64}(-6.0, 5.8125, 64), 0.1875), AbstractFloat[6.0819587326731085e-9 1.84073249236779e-8 … 5.378605883219747e-8 1.84073249236779e-8; 6.0819294756364805e-9 1.8407236375924915e-8 … 5.378580009635882e-8 1.8407236375924915e-8; … ; 6.081841986287847e-9 1.8406971585429215e-8 … 5.378502638061039e-8 1.8406971585429215e-8; 6.0819294756364805e-9 1.8407236375924915e-8 … 5.378580009635882e-8 1.8407236375924915e-8]), BSLSpline(5), CPU())
stepper = StrangSplitting()  # timestepper
dt = 0.1                     # timestep
nsteps = 1000                # total number of time-steps
1000
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))")
Example block output

Simulation with Lagrange interpolation

landau!(f, α, kx)

prob = VlasovProblem(f, BSLLagrange(9), dev)

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))")
Example block output