Simulation

using Plots
using VlasovPoissonTwoSpecies
import SpecialFunctions: ellipk

function run(coef, data)

     mesh_x = Mesh(data.x_min, data.x_max, data.nx)
     mesh_v = Mesh(data.v_min, data.v_max, data.nv)

     x = mesh_x.x
     v = mesh_v.x

     tf = data.T_final
     nt = data.nb_time_steps
     dt = tf / nt

     eq_manager = EquilibriumManager(coef, mesh_x, mesh_v)
     output = OutputManager(data, eq_manager)

     fe = perturbate(x, v, eq_manager.fe, data.perturbation_init)
     fi = perturbate(x, v, eq_manager.fi, data.perturbation_init)

     scheme = WellBalanced( fe, fi, eq_manager)

     rho = zeros(mesh_x.nx)
     e = zeros(mesh_x.nx)

     for it = 1:nt

         if it % data.freq_save == 0
             if data.output
                 save(output, scheme, it * dt)
             end
         end

         compute_source(scheme, 0.5dt)

         advect(scheme.advection_x, scheme.ge, v, 0.5dt)
         advect(scheme.advection_x, scheme.gi, v, 0.5dt)

         rho .= compute_rho(scheme)
         e .= compute_e(mesh_x, rho)

         e .+= scheme.e_eq

         advect(scheme.advection_v, transpose(scheme.ge), -e, dt)
         advect(scheme.advection_v, transpose(scheme.gi), e, dt)

         advect(scheme.advection_x, scheme.ge, v, 0.5dt)
         advect(scheme.advection_x, scheme.gi, v, 0.5dt)

         scheme.fe .= scheme.fe_eq .+ scheme.ge
         scheme.fi .= scheme.fi_eq .+ scheme.gi

         compute_source(scheme, 0.5dt)

     end

     output

end

run (generic function with 1 method)

refine_dt              = 1
refine_factor          = 1
eps_perturbation       = 0.001

coef = Coef()

data                   = Data()
data.projection        = false
data.projection_type   = :coefficients   # :BGK
data.T_final           = 200
data.nb_time_steps     = 1000 * refine_factor * refine_dt
data.nx                = 64 * refine_factor
data.nv                = 64 * refine_factor
data.x_min             = 0
data.x_max             = 4 * ellipk(coef.m)
data.v_min             = -10
data.v_max             = 10
data.perturbation_init = (x, v) -> eps_perturbation * (cos(2π * x / (data.x_max - data.x_min) + 1))
data.output            = true
data.vtk               = true
data.freq_save         = 5 * refine_factor * refine_dt
data.freq_output       = 5 * refine_factor * refine_dt
data.freq_projection   = 5 * 5 * refine_factor * refine_dt

@time output = run(coef, data)

OutputManager(Data(200, 1000, 64, 64, 0.0, 12.623499791567363, -10.0, 10.0, true, true, 5, 5, false, Main.var"#1#2"(), false, :coefficients, 25), Mesh([0.0, 0.19724218424324003, 0.39448436848648005, 0.5917265527297201, 0.7889687369729601, 0.9862109212162002, 1.1834531054594402, 1.3806952897026803, 1.5779374739459202, 1.7751796581891603  …  10.651077949134962, 10.848320133378202, 11.045562317621442, 11.242804501864683, 11.440046686107923, 11.637288870351163, 11.834531054594402, 12.031773238837642, 12.229015423080883, 12.426257607324123], 0.0, 12.623499791567363, 64, 0.19724218424324003), Mesh([-10.0, -9.6875, -9.375, -9.0625, -8.75, -8.437500000000002, -8.125, -7.8125, -7.500000000000001, -7.187500000000002  …  6.875, 7.187500000000001, 7.499999999999999, 7.8125, 8.124999999999998, 8.4375, 8.75, 9.0625, 9.374999999999998, 9.6875], -10.0, 10.0, 64, 0.3125), 0, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0  …  191.0, 192.0, 193.0, 194.0, 195.0, 196.0, 197.0, 198.0, 199.0, 200.0], false, 101.32205662204053, 135.44972697017792, [2.4055603796007123e-8, 8.854259683771559e-8, 1.979418883445158e-7, 2.427082740840795e-7, 2.836705481777226e-7, 2.718626057737916e-7, 2.768050662470789e-7, 2.757363270342309e-7, 2.959553660210037e-7, 3.3104049759782334e-7  …  0.0008535675680169586, 0.0010048472880817595, 0.0011259653045098859, 0.0012188815470815185, 0.0013099942253739253, 0.001415177018088093, 0.0015454681047745796, 0.0016998576213672707, 0.0018665911929881358, 0.002015754658543447], [2.1723222132145855e-8, 6.300740091430097e-8, 1.5208588427493454e-7, 1.7985149259705158e-7, 2.0428028268838097e-7, 2.0811352017597217e-7, 1.9706718011122757e-7, 2.0958704158318425e-7, 2.1358690112736694e-7, 2.4878538436073813e-7  …  0.00018712975129299756, 0.00028585749948714266, 0.00035285159218202184, 0.00040163117899417214, 0.0004568637778072643, 0.0005399351134201287, 0.0006516877634131419, 0.0007909511420676961, 0.000944788074523153, 0.0010797024208329543])

plot(output.t, output.energy_fe, label="electrons")
plot!(output.t, output.energy_fi, label="ions", legend=:topleft)