From 0b4925aa4738a61d1dca1cc521144ab6880b3260 Mon Sep 17 00:00:00 2001 From: Pierre Navaro Date: Sat, 11 Apr 2020 22:14:33 +0200 Subject: [PATCH] replace floor by trunc --- .gitignore | 4 ++++ {src => notebooks}/simulation.jl | 30 +++++++++++++----------------- notebooks/strong_landau_damping.jl | 14 ++++++-------- src/hamiltonian_splitting.jl | 1 + src/particle_mesh_coupling.jl | 16 ++++++++-------- src/particle_mesh_coupling_spin.jl | 16 ++++++++-------- 6 files changed, 40 insertions(+), 41 deletions(-) rename {src => notebooks}/simulation.jl (86%) diff --git a/.gitignore b/.gitignore index 4526a57..7b8c310 100644 --- a/.gitignore +++ b/.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 diff --git a/src/simulation.jl b/notebooks/simulation.jl similarity index 86% rename from src/simulation.jl rename to notebooks/simulation.jl index 34c7a59..72a0e98 100644 --- a/src/simulation.jl +++ b/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) diff --git a/notebooks/strong_landau_damping.jl b/notebooks/strong_landau_damping.jl index a6a6201..467a815 100644 --- a/notebooks/strong_landau_damping.jl +++ b/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") # - diff --git a/src/hamiltonian_splitting.jl b/src/hamiltonian_splitting.jl index 2cca0c9..f4a4870 100644 --- a/src/hamiltonian_splitting.jl +++ b/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) diff --git a/src/particle_mesh_coupling.jl b/src/particle_mesh_coupling.jl index 8d7997d..c7a44aa 100644 --- a/src/particle_mesh_coupling.jl +++ b/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) diff --git a/src/particle_mesh_coupling_spin.jl b/src/particle_mesh_coupling_spin.jl index c4575a5..38da7f4 100755 --- a/src/particle_mesh_coupling_spin.jl +++ b/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)