Example

Set physical and numerical parameters

using DynamicElectricSheath
using Plots

physics = Physics()

ν = physics.ν
λ = physics.λ
T = physics.T
μ = physics.μ
xmin, xmax = physics.xmin, physics.xmax
vmin, vmax = physics.vmin, physics.vmax

numerics = Discretization(physics)

Nt = numerics.Nt
dt = numerics.dt
CFL_x = numerics.CFL_x
CFL_v = numerics.CFL_v
Nx = numerics.Nx
dx = numerics.dx
Nv = numerics.Nv
dv = numerics.dv

println("Parameters : ")
println("ν = $ν, λ = $λ, μ = $μ")
println("T = $T, Nt = $Nt, dt = $dt (CFL_x = $CFL_x, CFL_v = $CFL_v)")
println("[xmin,xmax] = [$xmin,$xmax], Nx = $Nx, dx = $dx")
println("[vmin,vmax] = [$vmin,$vmax], Nv = $Nv, dv = $dv")

Parameters :
ν = 20.0, λ = 0.5, μ = 0.5
T = 10.0, Nt = 60000, dt = 0.00016666666666666666 (CFL_x = 0.0001666666666666667, CFL_v = 0.000999000999000999)
[xmin,xmax] = [-1.0,1.0], Nx = 600, dx = 0.0033333333333333335
[vmin,vmax] = [-10.0,10.0], Nv = 1001, dv = 0.01998001998001998

Set mesh grid

xx = LinRange(xmin, xmax, Nx + 1)
vv = LinRange(vmin, vmax, Nv + 1)
vv_plus = vv .* (vv .> 0.0)
vv_minus = vv .* (vv .< 0.0)

1002-element Vector{Float64}:
 -10.0
  -9.98001998001998
  -9.96003996003996
  -9.940059940059939
  -9.920079920079921
  -9.9000999000999
  -9.88011988011988
  -9.86013986013986
  -9.84015984015984
  -9.82017982017982
   ⋮
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0

Initialize fields

EE = E0.(xx)
fi = fi_0.(xx, vv')
fe = fe_0.(xx, vv')

601×1002 Matrix{Float64}:
 5.85739e-20  6.47213e-20  7.14995e-20  …  6.47213e-20  5.85739e-20
 6.26119e-20  6.91831e-20  7.64286e-20     6.91831e-20  6.26119e-20
 6.69284e-20  7.39525e-20  8.16975e-20     7.39525e-20  6.69284e-20
 7.15423e-20  7.90507e-20  8.73297e-20     7.90507e-20  7.15423e-20
 7.64744e-20  8.45004e-20  9.33501e-20     8.45004e-20  7.64744e-20
 8.17465e-20  9.03258e-20  9.97856e-20  …  9.03258e-20  8.17465e-20
 8.7382e-20   9.65528e-20  1.06665e-19     9.65528e-20  8.7382e-20
 9.3406e-20   1.03209e-19  1.14018e-19     1.03209e-19  9.3406e-20
 9.98454e-20  1.10324e-19  1.21878e-19     1.10324e-19  9.98454e-20
 1.06729e-19  1.1793e-19   1.30281e-19     1.1793e-19   1.06729e-19
 ⋮                                      ⋱  ⋮            
 9.98454e-20  1.10324e-19  1.21878e-19     1.10324e-19  9.98454e-20
 9.3406e-20   1.03209e-19  1.14018e-19     1.03209e-19  9.3406e-20
 8.7382e-20   9.65528e-20  1.06665e-19     9.65528e-20  8.7382e-20
 8.17465e-20  9.03258e-20  9.97856e-20  …  9.03258e-20  8.17465e-20
 7.64744e-20  8.45004e-20  9.33501e-20     8.45004e-20  7.64744e-20
 7.15423e-20  7.90507e-20  8.73297e-20     7.90507e-20  7.15423e-20
 6.69284e-20  7.39525e-20  8.16975e-20     7.39525e-20  6.69284e-20
 6.26119e-20  6.91831e-20  7.64286e-20     6.91831e-20  6.26119e-20
 5.85739e-20  6.47213e-20  7.14995e-20  …  6.47213e-20  5.85739e-20

p = plot(layout=(1,2), size=(800,300))
contourf!(p[1], xx, vv, fi', label="ions")
contourf!(p[2], xx, vv, fe', label="electrons")

ρ = zeros(Nx + 1)       # charge density
ρi = zeros(Nx + 1)        # ion charge density
ρe = zeros(Nx + 1)        # electron charge density
compute_charge!(ρi, fi, dv)
compute_charge!(ρe, fe, dv)
compute_charge!(ρ, fi .- fe, dv)
plot(xx, ρ)

# Boundary conditions (all 0, never updated afterwards)
fi[begin, :] .= 0.0
fi[end, :] .= 0.0 # speed distribution is almost 0 
fi[:, begin] .= 0.0
fi[:, end] .= 0.0 # non-emmiting wall
fe[begin, :] .= 0.0
fe[end, :] .= 0.0 # speed distribution is almost 0 
fe[:, begin] .= 0.0
fe[:, end] .= 0.0 # non-emmiting wall


for n = 1:Nt # loop over time

    compute_charge!(ρi, fi, dv)
    compute_charge!(ρe, fe, dv)
    compute_charge!(ρ, fi .- fe, dv)

    J_l, J_r = compute_current(fi, fe, vv, dv)

    EE_minus, EE_plus = compute_e!(EE, ρ, λ, J_l, J_r, dx, dt)

    advection!(fi, fe, vv_plus, vv_minus, EE_plus, EE_minus, ν, μ, dx, dv, dt)

end