Paper example

Test corresponding to Fig. 4 of the paper

using Plots
using Random
using SpinGEMPIC

import GEMPIC: OneDGrid, Maxwell1DFEM
import GEMPIC: l2projection!

function run_simulation( steps, Δt)

    σ, μ = 0.17, 0.0
    kx, α = 1.22, 0.02

    xmin, xmax = 0, 4pi/kx
    domain = [xmin, xmax, xmax - xmin]
    nx = 128
    n_particles = 20000
    mesh = OneDGrid( xmin, xmax, nx)
    spline_degree = 3

    df = CosGaussian(kx, α, σ, μ )

    rng = MersenneTwister(123)
    mass, charge = 1.0, 1.0

    particle_group = ParticleGroup( n_particles, mass, charge, 1)
    sample!(rng, particle_group, df, mesh, method = :quietstart)

    kernel_smoother2 = ParticleMeshCoupling( mesh, n_particles, spline_degree-2)
    kernel_smoother1 = ParticleMeshCoupling( mesh, n_particles, spline_degree-1)
    kernel_smoother0 = ParticleMeshCoupling( mesh, n_particles, spline_degree)

    maxwell_solver = Maxwell1DFEM(mesh, spline_degree)

    rho = zeros(nx)
    efield_poisson = zeros(nx)

    solve_poisson!( efield_poisson, particle_group, kernel_smoother0, maxwell_solver, rho )

    k0 = 2*kx
    E0 = 0.325
    ww = 2.63
    Ey(x) = E0*cos(k0*x)
    Ez(x) = E0*sin(k0*x)
    Ay(x) = -E0/ww*sin(k0*x)
    Az(x) = E0/ww*cos(k0*x)

    efield_dofs = [ zeros(nx), zeros(nx), zeros(nx)]
    efield_dofs[1] .= efield_poisson
    afield_dofs = [zeros(nx), zeros(nx)]

    l2projection!( efield_dofs[2], maxwell_solver, Ey, spline_degree)
    l2projection!( efield_dofs[3], maxwell_solver, Ez, spline_degree)
    l2projection!( afield_dofs[1], maxwell_solver, Ay, spline_degree)
    l2projection!( afield_dofs[2], maxwell_solver, Az, spline_degree)

    propagator = HamiltonianSplitting( maxwell_solver,
                                       kernel_smoother0,
                                       kernel_smoother1,
                                       kernel_smoother2,
                                       efield_dofs,
                                       afield_dofs,
                                       domain);

    thdiag = TimeHistoryDiagnostics( maxwell_solver,
                            kernel_smoother0, kernel_smoother1 )

    write_step!(thdiag, 0.0, spline_degree,
                        efield_dofs,  afield_dofs,
                        efield_poisson,
                        propagator, particle_group)

    for j = 1:steps # loop over time

        strang_splitting!(propagator, particle_group, Δt, 1)

        write_step!(thdiag, j * Δt, spline_degree,
                        efield_dofs,  afield_dofs,
                        efield_poisson,
                        propagator, particle_group)

    end

    return thdiag

end

steps, Δt = 200, 0.05

thdiag = run_simulation(steps, Δt)

time = thdiag.data[!, :Time]
energy = thdiag.data[!, :PotentialEnergyE1]
plot( time, 0.5*log.(energy.^2), xlabel = "time", ylabel = "Ex energy", label = "")