From e064d8345006502b7d7fe2f4dc40e4c92401af48 Mon Sep 17 00:00:00 2001 From: Luke Adams Date: Wed, 20 Mar 2024 14:48:42 -0600 Subject: [PATCH 1/7] add five pt stencil poisson solve --- src/poisson.jl | 61 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/src/poisson.jl b/src/poisson.jl index cdb265b..960b141 100644 --- a/src/poisson.jl +++ b/src/poisson.jl @@ -51,3 +51,64 @@ function step!(step::PoissonSolveFFT) inv(step.fft_plan) * step.ft_vector view(step.phi.values, eachindex(step.phi)) .= real.(step.ft_vector) end + +struct PoissonSolveFFT5{T,D,G,P,F<:AbstractField} <: AbstractSimulationStep + rho::F + phi::F + + Ksq_inv::Array{T,D} + ft_vector::Array{Complex{T},D} + + fft_plan::P + + function PoissonSolveFFT5(rho::F, phi::F) where {T,D,G,F<:AbstractField{T,D,G}} + # This restriction could possibly be relaxed to just require compatible grids... + @assert rho.grid === phi.grid + # Currently only supports periodic boundary conditions... + @assert all(rho.grid.periodic) + @assert num_elements(rho) == 1 && num_elements(phi) == 1 + + epsilon_0 = 8.8541878128e-12 + # Ksq_inv = reshape(copy(rho.values), size(rho.values)[1:end-1]...) + Ksq_inv = zeros(T, rho.grid.num_cells...) + ft_vector = zeros(Complex{T}, rho.grid.num_cells...) + + grid = rho.grid + # TODO: Make this use the cell_lengths logic in Grid + sim_lengths = grid.upper_bounds .- grid.lower_bounds + cell_lengths = sim_lengths ./ grid.num_cells + for I in eachindex(Ksq_inv) + It = Tuple(I) + + if any(x -> x == 1, It) + Ksq_inv[I] = 0 + continue + end + + # ks = 2π .* (It .- 1) ./ sim_lengths + ks = + 2π .* + ifelse.(It .<= size(Ksq_inv) ./ 2, It .- 1, It .- size(Ksq_inv) .- 1) ./ + sim_lengths + inv_Ksqs = + ( + ks .^ 2 .* sinc.(ks .* cell_lengths ./ 2 ./ pi) .^ 2 .* + (2 .+ cos.(ks .* cell_lengths)) ./ 3 + ) .^ (-1) ./ epsilon_0 + + Ksq_inv[I] = prod(inv_Ksqs) + end + + fft_plan = plan_fft!(ft_vector) + + new{T,D,G,typeof(fft_plan),F}(rho, phi, Ksq_inv, ft_vector, fft_plan) + end +end + +function step!(step::PoissonSolveFFT5) + step.ft_vector .= view(step.rho.values, eachindex(step.rho)) + step.fft_plan * step.ft_vector + step.ft_vector .= step.ft_vector .* step.Ksq_inv + inv(step.fft_plan) * step.ft_vector + view(step.phi.values, eachindex(step.phi)) .= real.(step.ft_vector) +end From 55dfca3c8342dd60afbd16885452bbccadc8b2b3 Mon Sep 17 00:00:00 2001 From: Luke Adams Date: Thu, 29 Aug 2024 16:15:24 -0600 Subject: [PATCH 2/7] add MultiplyField step --- src/ParticleInCell.jl | 2 +- src/field_utils.jl | 7 +++++++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/src/ParticleInCell.jl b/src/ParticleInCell.jl index b7d3040..814de14 100644 --- a/src/ParticleInCell.jl +++ b/src/ParticleInCell.jl @@ -44,7 +44,7 @@ end export step! include("field_utils.jl") -export FiniteDifferenceToEdges, AverageEdgesToNodes +export FiniteDifferenceToEdges, AverageEdgesToNodes, ZeroField, MultiplyField include("poisson.jl") export PoissonSolveFFT diff --git a/src/field_utils.jl b/src/field_utils.jl index 2a81e28..d747bdf 100644 --- a/src/field_utils.jl +++ b/src/field_utils.jl @@ -68,3 +68,10 @@ struct ZeroField{F} <: AbstractSimulationStep end step!(step::ZeroField) = step.field.values .= 0 + +struct MultiplyField{F,T} <: AbstractSimulationStep + field::F + x::T +end + +step!(step::MultiplyField) = step.field.values .*= step.x From 38b9b9fddc7519dcd7436487ce936485cd918114 Mon Sep 17 00:00:00 2001 From: Luke Adams Date: Thu, 29 Aug 2024 16:16:07 -0600 Subject: [PATCH 3/7] reduce number of guard cells in simulation creation helper function --- src/simulation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation.jl b/src/simulation.jl index f0c2aa6..900b425 100644 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -17,7 +17,7 @@ function create_electrostatic_simulation( sim = Simulation(AbstractSimulationStep[]) # Create fields - lower_guard_cells = div(interpolation_order, 2) + 1 + lower_guard_cells = div(interpolation_order, 2) rho = Field(grid, NodeOffset(), V, lower_guard_cells) phi = Field(grid, NodeOffset(), V, lower_guard_cells) Eedge = Field(grid, EdgeOffset(), V, lower_guard_cells) From e0ab4b7200dbfe1380838a78a7811fdaadb489a7 Mon Sep 17 00:00:00 2001 From: Luke Adams Date: Thu, 29 Aug 2024 16:17:01 -0600 Subject: [PATCH 4/7] add ability to create interpolation with an arbitrary shape function --- src/interpolation.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/interpolation.jl b/src/interpolation.jl index 174c92b..865c30d 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -163,6 +163,16 @@ struct BSplineChargeInterpolation{S,F,IF} <: AbstractSimulationStep interp_func, ) end + + function BSplineChargeInterpolation( + species::S, + rho::Field{T,N,NodeOffset}, + interp_width, + interp_func::F, + ) where {T,N,S,F} + + new{S,typeof(rho),F}(species, rho, -1, interp_width, interp_func) + end end function step!(step::BSplineChargeInterpolation) From edaff763912839a080a2a245933f4671d17f9247 Mon Sep 17 00:00:00 2001 From: Luke Adams Date: Thu, 29 Aug 2024 16:17:57 -0600 Subject: [PATCH 5/7] refactor poisson solve to allow for arbitrary eigenvalues --- src/poisson.jl | 72 ++++++++------------------------------------------ 1 file changed, 11 insertions(+), 61 deletions(-) diff --git a/src/poisson.jl b/src/poisson.jl index 960b141..e2ce48d 100644 --- a/src/poisson.jl +++ b/src/poisson.jl @@ -1,58 +1,8 @@ -struct PoissonSolveFFT{T,D,G,P,F<:AbstractField} <: AbstractSimulationStep - rho::F - phi::F - - Ksq_inv::Array{T,D} - ft_vector::Array{Complex{T},D} - - fft_plan::P - - function PoissonSolveFFT(rho::F, phi::F) where {T,N,O,D,G,F<:AbstractField{T,N,O,D,G}} - # This restriction could possibly be relaxed to just require compatible grids... - @assert rho.grid === phi.grid - # Currently only supports periodic boundary conditions... - @assert all(rho.grid.periodic) - @assert num_elements(rho) == 1 && num_elements(phi) == 1 - - epsilon_0 = 8.8541878128e-12 - # Ksq_inv = reshape(copy(rho.values), size(rho.values)[1:end-1]...) - Ksq_inv = zeros(T, rho.grid.num_cells...) - ft_vector = zeros(Complex{T}, rho.grid.num_cells...) - - grid = rho.grid - # TODO: Make this use the cell_lengths logic in Grid - sim_lengths = grid.upper_bounds .- grid.lower_bounds - cell_lengths = sim_lengths ./ grid.num_cells - for I in eachindex(Ksq_inv) - It = Tuple(I) +field_solve_three_pt(k, dx) = k^2 * sinc(k * dx / 2 / pi)^2 +field_solve_lagrange(k, dx) = k^2 * sinc(k * dx / 2 / pi)^2 * (2 + cos(k * dx)) / 3 +field_solve_five_pt(k, dx) = -1 * k^2 * sinc(k * dx / 2 / pi)^2 * (-7 + cos(k * dx)) / 6 - if any(x -> x == 1, It) - Ksq_inv[I] = 0 - continue - end - - ks = 2π .* (It .- 1) ./ sim_lengths - grid_angles = ks .* cell_lengths ./ 2 - inv_Ksqs = (cell_lengths ./ (2 .* sin.(grid_angles))) .^ 2 ./ epsilon_0 - - Ksq_inv[I] = prod(inv_Ksqs) - end - - fft_plan = plan_fft!(ft_vector) - - new{T,D,G,typeof(fft_plan),F}(rho, phi, Ksq_inv, ft_vector, fft_plan) - end -end - -function step!(step::PoissonSolveFFT) - step.ft_vector .= view(step.rho.values, eachindex(step.rho)) - step.fft_plan * step.ft_vector - step.ft_vector .= step.ft_vector .* step.Ksq_inv - inv(step.fft_plan) * step.ft_vector - view(step.phi.values, eachindex(step.phi)) .= real.(step.ft_vector) -end - -struct PoissonSolveFFT5{T,D,G,P,F<:AbstractField} <: AbstractSimulationStep +struct PoissonSolveFFT{T,D,G,P,F<:AbstractField} <: AbstractSimulationStep rho::F phi::F @@ -61,7 +11,11 @@ struct PoissonSolveFFT5{T,D,G,P,F<:AbstractField} <: AbstractSimulationStep fft_plan::P - function PoissonSolveFFT5(rho::F, phi::F) where {T,D,G,F<:AbstractField{T,D,G}} + function PoissonSolveFFT( + rho::F, + phi::F, + k2s = field_solve_three_pt, + ) where {T,D,G,F<:AbstractField{T,D,G}} # This restriction could possibly be relaxed to just require compatible grids... @assert rho.grid === phi.grid # Currently only supports periodic boundary conditions... @@ -90,11 +44,7 @@ struct PoissonSolveFFT5{T,D,G,P,F<:AbstractField} <: AbstractSimulationStep 2π .* ifelse.(It .<= size(Ksq_inv) ./ 2, It .- 1, It .- size(Ksq_inv) .- 1) ./ sim_lengths - inv_Ksqs = - ( - ks .^ 2 .* sinc.(ks .* cell_lengths ./ 2 ./ pi) .^ 2 .* - (2 .+ cos.(ks .* cell_lengths)) ./ 3 - ) .^ (-1) ./ epsilon_0 + inv_Ksqs = (k2s.(ks, cell_lengths)) .^ (-1) ./ epsilon_0 Ksq_inv[I] = prod(inv_Ksqs) end @@ -105,7 +55,7 @@ struct PoissonSolveFFT5{T,D,G,P,F<:AbstractField} <: AbstractSimulationStep end end -function step!(step::PoissonSolveFFT5) +function step!(step::PoissonSolveFFT) step.ft_vector .= view(step.rho.values, eachindex(step.rho)) step.fft_plan * step.ft_vector step.ft_vector .= step.ft_vector .* step.Ksq_inv From 0606706c1975fe3321fbde79adf44ce804e7fa7c Mon Sep 17 00:00:00 2001 From: Luke Adams Date: Thu, 29 Aug 2024 16:21:43 -0600 Subject: [PATCH 6/7] temporarily use first dimension for calculating E field location This is correct for node centered electric fields where all components are located at the same place in a cell, but wrong for edge electric fields. --- src/particle_updaters/electrostatic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/particle_updaters/electrostatic.jl b/src/particle_updaters/electrostatic.jl index 29bb943..9390798 100644 --- a/src/particle_updaters/electrostatic.jl +++ b/src/particle_updaters/electrostatic.jl @@ -56,7 +56,7 @@ function step!(step::ElectrostaticParticlePush) ) for I in Is - grid_cell_coord = cell_index_to_cell_coords(step.E, I) + grid_cell_coord = cell_index_to_cell_coords(step.E, I, 1) dist = Tuple(particle_cell_coord .- grid_cell_coord) interp_weights = step.interpolation_function.(dist) interp_weight = prod(interp_weights) From 666e228ceaa9a8ffed1aa857aa5196c504280066 Mon Sep 17 00:00:00 2001 From: Luke Adams Date: Thu, 29 Aug 2024 16:43:50 -0600 Subject: [PATCH 7/7] increase guard cells to account for averaging/finite differencing steps --- src/simulation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation.jl b/src/simulation.jl index 900b425..f0c2aa6 100644 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -17,7 +17,7 @@ function create_electrostatic_simulation( sim = Simulation(AbstractSimulationStep[]) # Create fields - lower_guard_cells = div(interpolation_order, 2) + lower_guard_cells = div(interpolation_order, 2) + 1 rho = Field(grid, NodeOffset(), V, lower_guard_cells) phi = Field(grid, NodeOffset(), V, lower_guard_cells) Eedge = Field(grid, EdgeOffset(), V, lower_guard_cells)