refactor the code

Project.toml

@@ -4,6 +4,7 @@ authors = "Julia Vlasov"
 version = "0.1.0"
 
 [deps]
+DispersionRelations = "78f632dc-6214-4aca-89df-9f86f5ac7c9e"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 GEMPIC = "b6d65c3a-4a4e-11e9-25d0-d309dc85ddeb"
@@ -14,6 +15,7 @@ Sobol = "ed01d8cd-4d21-5b2a-85b4-cc3bdc58bad4"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]
 FFTW = "1"

src/ParticleInCell.jl

@@ -14,21 +14,19 @@ include("particle_1d1v.jl")
 include("landau_damping.jl")
 include("poisson_1d.jl")
 include("poisson_2d.jl")
-include("particle_mesh.jl")
+include("cloud_in_cell.jl")
 include("pushers.jl")
 include("fdtd.jl")
 include("pic.jl")
 include("ua_type.jl")
 include("meshfields.jl")
 include("read_particles.jl")
 include("plasma.jl")
-include("compute_rho.jl")
-include("interpolation.jl")
 include("integrate.jl")
 include("gnuplot.jl")
-include("poisson.jl")
 include("ua_steps.jl")
-include("landau.jl")
 include("compute_rho_cic.jl")
+include("compute_rho_m6.jl")
+include("interpolation_m6.jl")
 
 end

src/compute_rho.jl → src/compute_rho_m6.jl

@@ -26,7 +26,7 @@ end
 
 export compute_rho_m6!
 
-function compute_rho_m6!( fields    :: MeshFields, 
+function compute_rho_m6!( fields    :: MeshFields2D, 
                           particles :: Particles, 
                           xt        :: Array{ComplexF64,3},
                           ua        :: UA ) 
@@ -178,7 +178,7 @@ function compute_rho_m6!( fields    :: MeshFields,
 
 end 
 
-function compute_rho_m6!( fields  :: MeshFields, particles :: Particles) 
+function compute_rho_m6!( fields  :: MeshFields2D, particles :: Particles) 
 
     fill!(fields.ρ , 0.0)
     nx = fields.mesh.nx

src/gnuplot.jl

@@ -1,7 +1,7 @@
 export gnuplot
 export errors
 
-function gnuplot( filename :: String, fields :: MeshFields )
+function gnuplot( filename :: String, fields :: MeshFields2D )
 
     open(filename, "w") do f
 
@@ -26,7 +26,7 @@ function gnuplot( filename :: String, fields :: MeshFields )
 
 end
 
-function errors( computed :: MeshFields, reference :: MeshFields )
+function errors( computed :: MeshFields2D, reference :: MeshFields2D )
 
     err = maximum(abs.(computed.e .- reference.e))
 

src/interpolation.jl → src/interpolation_m6.jl

@@ -1,7 +1,7 @@
 export interpol_eb_m6!
 
 function interpol_eb_m6!( e      :: Array{Float64,3}, 
-                          fields :: MeshFields, 
+                          fields :: MeshFields2D, 
                           x      :: Array{ComplexF64,3},
                           nbpart :: Int64,
                           ntau   :: Int64)
@@ -122,7 +122,7 @@ function interpol_eb_m6!( e      :: Array{Float64,3},
 
 end 
 
-function interpol_eb_m6!( particles :: Particles, fields :: MeshFields )
+function interpol_eb_m6!( particles :: Particles, fields :: MeshFields2D )
 
     nx = fields.mesh.nx
     ny = fields.mesh.ny

src/landau_damping.jl

@@ -124,3 +124,4 @@ function sample!(pg::ParticleGroup{1,1}, mesh::OneDGrid, d::LandauDamping)
     end
 
 end
+

src/meshfields.jl

@@ -1,13 +1,13 @@
-export MeshFields
+export MeshFields2D
 
-struct MeshFields
+struct MeshFields2D
 
     mesh :: TwoDGrid
 
     e :: Array{Float64,3}
     ρ :: Array{Float64,2}
 
-    function MeshFields( mesh :: TwoDGrid )
+    function MeshFields2D( mesh :: TwoDGrid )
 
 	     nx, ny = mesh.nx, mesh.ny
 
@@ -20,3 +20,89 @@ struct MeshFields
 
 end
 
+using FFTW
+
+export Poisson2D
+
+struct Poisson2D
+
+    mesh   :: TwoDGrid
+    kx     :: Array{Float64, 2}
+    ky     :: Array{Float64, 2}
+    ρ̃      :: Array{ComplexF64, 2}
+    fft_ex :: Array{ComplexF64, 2}
+    fft_ey :: Array{ComplexF64, 2}
+
+    function Poisson2D( mesh :: TwoDGrid )
+
+        nx = mesh.nx
+        ny = mesh.ny
+
+        kx0 = 2π / (mesh.xmax - mesh.xmin)
+        ky0 = 2π / (mesh.ymax - mesh.ymin)
+
+        kx = zeros(Float64, (nx÷2+1,ny))
+        ky = zeros(Float64, (nx÷2+1,ny))
+
+        for ik=1:nx÷2+1
+           kx1 = (ik-1)*kx0
+           for jk = 1:ny÷2
+              kx[ik,jk] = kx1
+              ky[ik,jk] = (jk-1)*ky0
+           end
+           for jk = ny÷2+1:ny
+              kx[ik,jk] = kx1
+              ky[ik,jk] = (jk-1-ny)*ky0
+           end 
+        end
+
+        kx[1,1] = 1.0
+        k2  = kx .* kx .+ ky .* ky
+        kx .= kx ./ k2
+        ky .= ky ./ k2
+
+        ρ̃ = zeros(ComplexF64,(nx÷2+1,ny))
+
+        fft_ex = similar(ρ̃)
+        fft_ey = similar(ρ̃)
+
+        new( mesh, kx, ky, ρ̃, fft_ex, fft_ey)
+
+    end	    
+
+end
+
+
+"""
+
+   poisson!( fields )
+
+Solve the equation Δ Φ = - fields.ρ
+
+ fields.ex = ∂ Φ / ∂ x
+ fields.ey = ∂ Φ / ∂ y 
+
+
+"""
+function solve!( fields :: MeshFields2D, p :: Poisson2D )
+
+    nx, ny = p.mesh.nx, p.mesh.ny
+    dx, dy = p.mesh.dx, p.mesh.dy
+
+    p.ρ̃ .= rfft(view(fields.ρ,1:nx,1:ny))
+
+    p.fft_ex .= -1im .* p.kx .* p.ρ̃
+    p.fft_ey .= -1im .* p.ky .* p.ρ̃
+
+    fields.e[1,1:nx,1:ny] .= irfft(p.fft_ex, nx)
+    fields.e[2,1:nx,1:ny] .= irfft(p.fft_ey, nx)
+
+    fields.e[1,nx+1,:] .= view(fields.e,1,1,:)
+    fields.e[1,:,ny+1] .= view(fields.e,1,:,1)
+    fields.e[2,nx+1,:] .= view(fields.e,2,1,:)
+    fields.e[2,:,ny+1] .= view(fields.e,2,:,1)
+
+    @views sum(  fields.e[1,:,:] .* fields.e[1,:,:] 
+              .+ fields.e[2,:,:] .* fields.e[2,:,:]) * dx * dy
+
+end

src/particles.jl

@@ -34,4 +34,50 @@ mutable struct Particles
 
 end
 
+export landau_sampling
+
+function landau_sampling( mesh :: TwoDGrid, nbpart :: Int64 )
+
+    kx     = 0.5
+    alpha  = 0.05
+
+    nx, ny = mesh.nx, mesh.ny
+    dx, dy = mesh.dx, mesh.dy
+    dimx = mesh.xmax - mesh.xmin
+    dimy = mesh.ymax - mesh.ymin
+
+    weight = (dimx * dimy) / nbpart
+
+    particles = Particles(nbpart, weight )
+
+    function newton(r)
+        x0, x1 = 0.0, 1.0
+        r *= 2π / kx
+        while (abs(x1-x0) > 1e-12)
+            p = x0 + alpha * sin( kx * x0) / kx 
+            f = 1 + alpha * cos( kx * x0)
+            x0, x1 = x1, x0 - (p - r) / f
+        end
+        x1
+    end
+
+    s = Sobol.SobolSeq(3)
+    nbpart = pg.n_particles
+
+    for i=1:nbpart
+
+        v = sqrt(-2 * log( (i-0.5)/nbpart))
+        r1, r2, r3  = Sobol.next!(s)
+        θ = r1 * 2π
+        particles.x[1,i] = xmin + newton(r2) * ( xmax - xmin )
+        particles.x[2,i] = ymin + r3 * ( ymax - ymin )
+        particles.v[1,i] = v * cos(θ)
+        particles.v[2,i] = v * sin(θ)
+
+    end
+
+    particles
+
+end
+