2.2 Simulation

function vlasov_ampere(nx, nv, xmin, xmax, vmin, vmax, tf, nt)
    """
    Full simulation of the Vlasov–Ampere system
    
    Parameters:
    - nx, nv: Number of spatial and velocity grid points
    - xmin, xmax: Spatial domain boundaries
    - vmin, vmax: Velocity domain boundaries  
    - tf: Final time
    - nt: Number of time steps
    
    Returns:
    - f: Final distribution function
    - nrj: Array of log-energy at each time step (for diagnostics)
    """

    # Create mesh objects
    meshx = UniformMesh(xmin, xmax, nx)
    meshv = UniformMesh(vmin, vmax, nv)
            
    # Prepare grid points
    x = meshx.points
    v = meshv.points
    
    # Landau damping parameters
    ϵ, kx = 0.001, 0.5   # Small amplitude perturbation at wavenumber kx
    
    # Allocate arrays
    f = zeros(Complex{Float64}, (nx, nv))    # Distribution: f[x_index, v_index]
    fᵀ = zeros(Complex{Float64}, (nv, nx))   # Transposed: fᵀ[v_index, x_index]
    
    # Use complex arithmetic for FFT efficiency
    # (We'll take real parts when needed)
    
    # Initialize distribution function
    f .= landau(ϵ, kx, x, v)
    
    # Transpose for velocity advection step
    transpose!(fᵀ, f)
    
    # Compute initial charge density and electric field
    ρ  = compute_rho(meshv, f)
    e  = zeros(ComplexF64, nx)
    e .= compute_e(meshx, ρ)
    
    # Storage for energy history (for analysis of Landau damping)
    nrj = Float64[]
    
    # Time step size
    dt = tf / nt
    
    # Create advection operators (stores Fourier wavenumbers)
    advection_x! = AmpereAdvection(meshx)
    advection_v! = AmpereAdvection(meshv)
    
    # Time stepping loop (Strang splitting)
    for i in 1:nt
        
        # Half-step velocity advection (using transposed array for cache efficiency)
        advection_v!(fᵀ, e, 0.5dt)
        
        # Transpose back for space advection
        transpose!(f, fᵀ)
        
        # Full-step space advection and field update
        advection_x!(f, e, v, dt)
        
        # Record energy: E_field = (1/2) ∫ E² dx
        # log(E_field) should decay linearly in time for Landau damping
        push!(nrj, log(sqrt((sum(e.^2)) * meshx.step)))
        
        # Transpose for next velocity advection
        transpose!(fᵀ, f)
        
        # Second half-step velocity advection
        advection_v!(fᵀ, e, 0.5dt)
        
    end
    
    # Return final distribution and energy history
    real(f), nrj
    
end

Time Integration Scheme (Strang Splitting):

The algorithm uses second-order Strang splitting to solve the coupled equations:

\[\frac{df}{dt} = -v\frac{\partial f}{\partial x} + E(x)\frac{\partial f}{\partial v}\]

At each time step:

1. Half velocity advection: \(0.5 \Delta t\)
2. Full space advection: \(\Delta t\) (with simultaneous field update)
3. Half velocity advection: \(0.5 \Delta t\)

This is second-order accurate and better than sequential splitting.

2.2.1 Running the Simulation

# Grid and domain parameters
nx, nv = 256, 256           # 256 points in space and velocity
xmin, xmax = 0., 4*π        # Spatial domain: [0, 4π]
vmin, vmax = -6., 6.        # Velocity domain: [-6, 6]

# Time parameters
tf = 60                     # Final time
nt = 600                    # Number of time steps
t = range(0, stop=tf, length=nt)  # Time array

# Theoretical prediction: Landau damping gives exponential decay
# E(t) ∝ exp(-γt) where γ ≈ 0.1533 (for our parameters)
# Plot the theoretical decay line
plot(t, -0.1533*t .- 5.48, label="Theory")

# Run the simulation
f, nrj = vlasov_ampere(nx, nv, xmin, xmax, vmin, vmax, tf, nt)

# Plot the simulated energy decay and compare with theory
plot!(t, nrj, label="Simulation")

• The log-energy should decay linearly with time: \(\log E(t) \approx -\gamma t + \text{const}\)
• The slope \(\gamma \approx 0.1533\) is the Landau damping rate
• When plotted on a log-linear scale, Landau damping appears as a straight line
• Deviations from linearity indicate numerical errors or nonlinear effects