Add examples and fix documentation

docs/make.jl

@@ -18,7 +18,6 @@ makedocs(;
             ),
         ),
     ),
-    modules=[GEMPIC],
     pages=[
         "Documentation" => [
             "index.md",

examples/distributions.jl

@@ -0,0 +1,62 @@
+# -*- coding: utf-8 -*-
+# ---
+# jupyter:
+#   jupytext:
+#     comment_magics: false
+#     formats: ipynb,jl:light
+#     text_representation:
+#       extension: .jl
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.4.2
+#   kernelspec:
+#     display_name: Julia 1.4.0
+#     language: julia
+#     name: julia-1.4
+# ---
+
+# ## Initialize a distribution
+
+using GEMPIC
+
+@doc SumCosGaussian
+
+# +
+using Plots
+
+k, α = 0.5, 0.1
+nx, nv = 64, 128
+xmin, xmax = 0, 2π/k
+vmin, vmax = -6, 6
+xg = LinRange(xmin, xmax, nx+1)[1:end-1] 
+vg = LinRange(vmin, vmax, nv+1)[1:end-1];
+
+# -
+
+df_1d1v = SumCosGaussian{1,1}( [[k]] , [α] , [[1.0]] , [[0.0]], [1.0])
+f = zeros(Float64, (nx,nv))
+for  j in eachindex(vg), i in eachindex(xg)
+    f[i,j] = df_1d1v( xg[i], vg[j])
+end
+surface(xg, vg, f')
+
+two_gaussians = SumCosGaussian{1,1}( [[k]] , [0.0] , [[0.2],[0.2]] , [[-1.0],[1.0]], [0.5, 0.5])
+f = zeros(Float64, (nx,nv))
+for  j in eachindex(vg), i in eachindex(xg)
+    f[i,j] = two_gaussians( xg[i], vg[j])
+end
+surface(xg, vg, f')
+
+# +
+fxv(x, v) = ( 1 + α * cos(k * x )) / sqrt(2π) * exp(- (v^2)/ 2)
+
+surface(xg, vg, fxv)
+# -
+
+@doc CosSumGaussian
+
+df = CosSumGaussian{1,1}( [[k]], [0.1], [[1.0]], [[0.0]], [1.0] )
+
+plot( xg, [GEMPIC.eval_x_density(df, x) for x in xg])
+
+plot( vg, [GEMPIC.eval_v_density(df, v) for v in vg])

examples/maxwell_1d_fem.jl

@@ -0,0 +1,150 @@
+# -*- coding: utf-8 -*-
+# ---
+# jupyter:
+#   jupytext:
+#     comment_magics: false
+#     formats: ipynb,jl:light
+#     text_representation:
+#       extension: .jl
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.7.1
+#   kernelspec:
+#     display_name: Julia 1.5.3
+#     language: julia
+#     name: julia-1.5
+# ---
+
+# # Maxwell solver FEM 1D
+
+using GEMPIC, LinearAlgebra, Plots
+
+const mode = 2
+
+eta1_min = .0
+eta1_max = 2π
+nc_eta1  = 128
+Lx = eta1_max - eta1_min
+delta_eta1 = Lx / nc_eta1
+deg = 3
+
+mesh = Mesh( eta1_min , eta1_max, nc_eta1)
+maxwell_1d = Maxwell1DFEM(mesh, deg);
+
+ex = zeros(Float64, nc_eta1)
+ey = zeros(Float64, nc_eta1)
+bz = zeros(Float64, nc_eta1);
+
+bz_exact = zeros(Float64, nc_eta1)
+ex_exact = zeros(Float64, nc_eta1)
+ey_exact = zeros(Float64, nc_eta1)
+rho      = zeros(Float64, nc_eta1)
+sval     = zeros(Float64, nc_eta1);
+
+cos_k(x) = cos(mode*2*pi*x/Lx) 
+sin_k(x) = sin(mode*2*pi*x/Lx) 
+
+# Test Poisson
+#-------------
+# Set exact solution
+for i = 1:nc_eta1
+   xi = eta1_min + (i-1)*delta_eta1
+   ex_exact[i] = sin_k(xi)/(2.0*mode*pi/Lx)
+end
+
+x = LinRange(eta1_min, eta1_max, nc_eta1+1)[1:end-1]
+
+compute_rhs_from_function!( rho, maxwell_1d, cos_k, deg)
+
+compute_e_from_rho!( ex, maxwell_1d, rho ) 
+
+plot(x, ex_exact)
+sval = eval_uniform_periodic_spline_curve(deg-1, ex)
+scatter!(x, sval)
+
+dt = .5 * delta_eta1
+
+for i = 1:nc_eta1
+   xi = eta1_min + (i-1)*delta_eta1
+   ex_exact[i] = - cos_k(xi)*dt
+end
+
+compute_rhs_from_function!(rho, maxwell_1d, cos_k, deg-1)
+fill!(ex, 0.0)
+compute_e_from_j!(ex, maxwell_1d, dt .* rho, 1 )
+sval =  eval_uniform_periodic_spline_curve(deg-1, ex);
+
+plot(x, ex_exact)
+scatter!(x, sval)
+
+# +
+mode     = 2
+eta1_min = .0
+eta1_max = 2π
+nc_eta1  = 256
+
+Lx = eta1_max - eta1_min
+
+delta_eta1 = Lx / nc_eta1
+
+deg = 3
+
+maxwell_1d = Maxwell1DFEM(mesh, deg)
+
+ex = zeros(Float64, nc_eta1)
+ey = zeros(Float64, nc_eta1)
+bz = zeros(Float64, nc_eta1)
+
+bz_exact = zeros(Float64, nc_eta1)
+ex_exact = zeros(Float64, nc_eta1)
+ey_exact = zeros(Float64, nc_eta1)
+rho      = zeros(Float64, nc_eta1)
+sval     = zeros(Float64, nc_eta1)
+
+cos_k(x) = cos(mode*2*pi*x/Lx) 
+sin_k(x) = sin(mode*2*pi*x/Lx) 
+
+# +
+# Test Maxwell on By and Ez 
+#--------------------------
+# Set time stepping parameters
+
+time  = 0.0
+dt    = .5 * delta_eta1
+nstep = 10
+
+# Compute initial fields 
+ex = 0.0 # 0-form -> splines of degree deg
+l2projection!(bz, maxwell_1d, cos_k, deg-1) # 0-form -> splines of degree deg-1
+
+#plot(bz)
+
+# +
+
+compute_b_from_e!( bz, maxwell_1d, 0.5*dt, ey)
+
+plot(bz)
+# -
+
+compute_e_from_b!( ey, maxwell_1d,     dt, bz)
+plot(ey)
+
+# +
+compute_b_from_e!( bz, maxwell_1d, 0.5*dt, ey)
+
+time = time + dt
+
+for i = 1:nc_eta1
+   xi = eta1_min + (i-1)*delta_eta1
+   ey_exact[i] = sin(mode * 2π * xi/Lx) * sin(mode * 2π * time/Lx)
+   bz_exact[i] = cos(mode * 2π * xi/Lx) * cos(mode * 2π * time/Lx)
+end
+
+sval = eval_uniform_periodic_spline_curve(deg, ey)
+@show err_ey = norm(sval .- ey_exact)
+
+sval = eval_uniform_periodic_spline_curve(deg-1, bz)
+@show err_bz = norm(sval .- bz_exact)
+# -
+
+

examples/vp1d1v.jl

@@ -0,0 +1,64 @@
+using ParticleInCell
+using Plots
+using Random
+
+dt = 0.1
+nsteps = 100
+alpha = 0.5
+kx = 0.5
+
+nx = 64
+xmin, xmax = 0.0, 2π / kx
+
+n_particles = 1000000
+degree_smoother = 3
+
+mesh = OneDGrid( xmin, xmax, nx)
+
+particles = ParticleGroup{1,1}( n_particles, charge=1.0, mass=1.0, n_weights=1)
+
+sampler = LandauDamping( alpha, kx )
+
+ParticleInCell.sample!( particles, mesh, sampler)
+
+particles.array[3,:] .= (xmax - xmin) ./ n_particles;
+
+poisson = OneDPoisson( mesh )
+kernel = ParticleMeshCoupling1D( mesh, n_particles, degree_smoother, :collocation)
+
+ex = zeros(nx)
+rho = zeros(nx)
+
+for i_part = 1:particles.n_particles
+    xi = particles.array[1, i_part]
+    wi = particles.array[3, i_part]
+    GEMPIC.add_charge!(rho, kernel, xi, wi)
+end
+
+compute_e_from_rho!(ex, poisson, rho)
+
+problem = OneDPoissonPIC( poisson, kernel )
+
+dt = 0.1
+nsteps = 100
+alpha = 0.1
+kx = 0.5
+
+propagator = OneDSplittingOperator( problem, particles ) 
+
+energy = Float64[]
+
+for j=1:nsteps
+
+    ParticleInCell.operator_t!(propagator, 0.5dt)
+    ParticleInCell.charge_deposition!(propagator)
+    ParticleInCell.solve_fields!(propagator)
+    ParticleInCell.operator_v!(propagator, dt)
+    ParticleInCell.operator_t!(propagator, 0.5dt)
+
+    push!(energy, compute_field_energy(problem))
+
+end
+
+t = collect(0:nsteps) .* dt
+plot(log.(energy))