replace floor by trunc

.gitignore

@@ -18,3 +18,7 @@ notebooks/quietstart.ipynb
 *.jld2
 notebooks/maxwell_1d_fem.ipynb
 notebooks/strong_landau_damping_spin.ipynb
+charge_density.png
+electric_field.png
+histograms.png
+potential_energy.png

src/simulation.jl → notebooks/simulation.jl

@@ -1,16 +1,15 @@
-
-include("mesh.jl")
-include("low_level_bsplines.jl")
-include("splinepp.jl")
-include("distributions.jl")
-include("maxwell_1d_fem.jl")
-include("particle_group.jl")
-include("particle_mesh_coupling.jl")
-include("particle_sampling.jl")
-include("landau_damping.jl")
-include("hamiltonian_splitting.jl")
-include("hamiltonian_splitting_boris.jl")
-include("diagnostics.jl")
+include("../src/mesh.jl")
+include("../src/low_level_bsplines.jl")
+include("../src/splinepp.jl")
+include("../src/distributions.jl")
+include("../src/maxwell_1d_fem.jl")
+include("../src/particle_group.jl")
+include("../src/particle_mesh_coupling.jl")
+include("../src/particle_sampling.jl")
+include("../src/landau_damping.jl")
+include("../src/hamiltonian_splitting.jl")
+include("../src/hamiltonian_splitting_boris.jl")
+include("../src/diagnostics.jl")
 
 const β = 0.0001
 const k  = 1.25
@@ -105,7 +104,4 @@ function run_simulation( steps )
 end
 
 run_simulation(1)
-using Profile
-Profile.clear()
-@profile run_simulation(500)
-Profile.print()
+@time run_simulation(500)

notebooks/strong_landau_damping.jl

@@ -59,28 +59,25 @@ histogram!(p[2,1], vp[1,:], weights=wp, normalize=true, bins = 100, lab = "")
 plot!(p[2,1], v-> exp( - v^2 / 2) * 4 / π^2 , -6, 6, lab="")
 histogram!(p[3,1], vp[2,:], weights=wp, normalize=true, bins = 100, lab = "")
 plot!(p[3,1], v-> exp( - v^2 / 2) * 4 / π^2 , -6, 6, lab="")
+savefig("histograms.png")
 
 
 kernel_smoother1 = ParticleMeshCoupling( domain, [nx], n_particles, spline_degree-1, :galerkin)    
 kernel_smoother0 = ParticleMeshCoupling( domain, [nx], n_particles, spline_degree, :galerkin)
 rho = zeros(Float64, nx)
-ex = zeros(Float64, nx)
-maxwell_solver = Maxwell1DFEM(domain, nx, spline_degree)
-solve_poisson!( ex, particle_group, 
-                kernel_smoother0, maxwell_solver, rho)
 xg = LinRange(xmin, xmax, nx)
 sval = eval_uniform_periodic_spline_curve(spline_degree-1, rho)
 plot( xg, sval )
+savefig("charge_density.png")
 
 efield_poisson = zeros(Float64, nx)
-# Init!ialize the field solver
+# Initialize the field solver
 maxwell_solver = Maxwell1DFEM(domain, nx, spline_degree)
 # efield by Poisson
-solve_poisson!( efield_poisson, particle_group, kernel_smoother0, maxwell_solver, rho )
-
 solve_poisson!( efield_poisson, particle_group, kernel_smoother0, maxwell_solver, rho )
 sval = eval_uniform_periodic_spline_curve(spline_degree-1, efield_poisson)
 plot( xg, sval )       
+savefig("electric_field.png")
 
 # Initialize the arrays for the spline coefficients of the fields
 
@@ -132,9 +129,10 @@ end
 
 # +
 
-@time results = run(1) # change number of steps
+@time results = run(10) # change number of steps
 
 plot(results[!,:Time], log.(results[!,:PotentialEnergyE1]))
+savefig("potential_energy.png")
 # -
 
 

src/hamiltonian_splitting.jl

@@ -209,6 +209,7 @@ function operatorHp1(h :: HamiltonianSplitting, dt :: Float64)
                     wi[1])
 
        x_new[1] = mod(x_new[1], h.Lx)
+
        set_x(h.particle_group, i_part, x_new)
        set_v(h.particle_group, i_part, vi)
 

src/particle_mesh_coupling.jl

@@ -101,7 +101,7 @@ function add_charge_pp!(rho_dofs :: Vector{Float64},
                         marker_charge)
 
     xi    = (position[1] - p.domain[1])/p.delta_x
-    index = floor(Int64, xi)+1
+    index = trunc(Int, xi)+1
     xi    = xi    - (index-1)
     index = index - p.spline_degree
 
@@ -133,7 +133,7 @@ function add_charge!( rho_dofs      :: Vector{Float64},
                       marker_charge :: Float64) 
 
     xi :: Float64 = (position[1] - p.domain[1])/p.delta_x[1]
-    index :: Int = floor(Int64, xi)+1
+    index :: Int = trunc(Int, xi)+1
     xi = xi - (index-1)
     index = index - p.spline_degree
 
@@ -167,12 +167,12 @@ function add_current_update_v_pp!( j_dofs        :: AbstractArray,
     # contributes to, and r_old, its position (normalized to cell size one).
 
     xi = (position_old[1] - p.domain[1]) / p.delta_x[1]
-    index_old = floor(Int64, xi)
+    index_old = trunc(Int, xi)
     r_old = xi - index_old
 
     # Compute the new box index index_new and normalized position r_old.
     xi = (position_new[1] - p.domain[1]) / p.delta_x[1]
-    index_new = floor(Int64, xi)
+    index_new = trunc(Int, xi)
     r_new = xi - index_new
 
     if index_old == index_new
@@ -275,13 +275,13 @@ function add_current_update_v!( j_dofs        :: AbstractArray,
 
 
     xi = (position_old[1] - p.domain[1]) / p.delta_x[1]
-    index_old = floor(Int64,xi)
+    index_old = trunc(Int,xi)
     r_old = xi - index_old
 
     # Compute the new box index index_new and normalized position r_old.
 
     xi = (position_new[1] - p.domain[1]) / p.delta_x[1]
-    index_new = floor(Int64, xi)
+    index_new = trunc(Int, xi)
     r_new = xi - index_new
 
     if index_old == index_new
@@ -390,7 +390,7 @@ function evaluate_pp(p             :: ParticleMeshCoupling,
                      field_dofs_pp :: Array{Float64,2})
 
     xi = (position[1] - p.domain[1])/p.delta_x[1]
-    index = floor(Int64, xi)+1
+    index = trunc(Int, xi)+1
     xi = xi - (index-1)
 
     horner_1d(p.spline_degree, field_dofs_pp, xi, index)
@@ -411,7 +411,7 @@ function evaluate(p          :: ParticleMeshCoupling,
                   field_dofs :: Vector{Float64})
 
     xi = (position[1] - p.domain[1])/p.delta_x[1]
-    index = floor(Int64, xi)+1
+    index = trunc(Int, xi)+1
     xi = xi - (index-1)
     index = index - p.spline_degree
     p.spline_val .= uniform_bsplines_eval_basis(p.spline_degree, xi)

src/particle_mesh_coupling_spin.jl

@@ -99,7 +99,7 @@ function add_charge_pp!(rho_dofs :: Vector{Float64},
                         marker_charge)
 
     xi    = (position[1] - p.domain[1])/p.delta_x
-    index = floor(Int64, xi)+1
+    index = trunc(Int64, xi)+1
     xi    = xi    - (index-1)
     index = index - p.spline_degree
 
@@ -134,7 +134,7 @@ function add_charge!( rho_dofs      :: Vector{Float64},
                       marker_charge :: Float64) 
 
     xi    = (position[1] - p.domain[1])/p.delta_x[1]
-    index = floor(Int64, xi)+1
+    index = trunc(Int64, xi)+1
     xi    = xi - (index-1)
     index = index - p.spline_degree
 
@@ -168,12 +168,12 @@ function add_current_update_v_pp!( j_dofs        :: AbstractArray,
     # contributes to, and r_old, its position (normalized to cell size one).
 
     xi = (position_old[1] - p.domain[1]) / p.delta_x[1]
-    index_old = floor(Int64, xi)
+    index_old = trunc(Int64, xi)
     r_old = xi - index_old
 
     # Compute the new box index index_new and normalized position r_old.
     xi = (position_new[1] - p.domain[1]) / p.delta_x[1]
-    index_new = floor(Int64, xi)
+    index_new = trunc(Int64, xi)
     r_new = xi - index_new
 
     if index_old == index_new
@@ -275,13 +275,13 @@ function add_current_update_v!( j_dofs        :: AbstractArray,
 
 
     xi = (position_old[1] - p.domain[1]) / p.delta_x[1]
-    index_old = floor(Int64,xi)
+    index_old = trunc(Int64,xi)
     r_old = xi - index_old
 
     # Compute the new box index index_new and normalized position r_old.
 
     xi = (position_new[1] - p.domain[1]) / p.delta_x[1]
-    index_new = floor(Int64, xi)
+    index_new = trunc(Int64, xi)
     r_new = xi - index_new
 
     if index_old == index_new
@@ -393,7 +393,7 @@ function evaluate_pp(p             :: SpinParticleMeshCoupling,
                      field_dofs_pp :: Array{Float64,2})
 
     xi = (position[1] - p.domain[1])/p.delta_x[1]
-    index = floor(Int64, xi)+1
+    index = trunc(Int64, xi)+1
     xi = xi - (index-1)
 
     horner_1d(p.spline_degree, field_dofs_pp, xi, index)
@@ -414,7 +414,7 @@ function evaluate(p          :: SpinParticleMeshCoupling,
                   field_dofs :: Vector{Float64})
 
     xi = (position[1] - p.domain[1])/p.delta_x[1]
-    index = floor(Int64, xi)+1
+    index = trunc(Int64, xi)+1
     xi = xi - (index-1)
     index = index - p.spline_degree
     p.spline_val .= uniform_bsplines_eval_basis(p.spline_degree, xi)