Generalize PIC deposition to arbitrary B-spline order (P=1,2,3)

src/Models/WaveGrowthModels1D.jl

@@ -50,7 +50,8 @@ function WaveGrowth1D(; grid::CartesianGridMesh1D,
     currents=nothing,
     periodic_boundary=true,
     boundary_type="same",
-    CBsets=nothing)
+    CBsets=nothing,
+    spline_order::Int=1)
 
     return WaveGrowth2D(;
         grid=grid,
@@ -65,7 +66,8 @@ function WaveGrowth1D(; grid::CartesianGridMesh1D,
         currents=currents,
         periodic_boundary=periodic_boundary,
         boundary_type=boundary_type,
-        CBsets=CBsets)
+        CBsets=CBsets,
+        spline_order=spline_order)
 end
 
 # end of module

src/Models/WaveGrowthModels2D.jl

@@ -88,6 +88,8 @@ mutable struct WaveGrowth2D{Grid<:AbstractGrid,
 
     MovieState::Mstat     # state of of the model. Only used for producing movieframes
 
+    spline_order::Base.Int  # B-spline order for particle->grid deposition (1=CIC, 2=TSC, 3=cubic)
+
 end
 
 function Base.getproperty(ow::WaveGrowth2D, s::Symbol)
@@ -214,12 +216,22 @@ function WaveGrowth2D(; grid::GG,
     periodic_boundary=true,
     boundary_type="same", # or "minimal", "same", default is same, only used if periodic_boundary is false
     CBsets=nothing,
+    spline_order::Int=1,  # B-spline deposition order: 1=CIC (default), 2=TSC, 3=cubic
     movie=false) where {PP<:Union{ParticleDefaults2D,String},GG<:AbstractGrid}
 
     if !isnothing(winds) && !(haskey(winds, :u) && haskey(winds, :v))
         error("`winds` must provide fields `u` and `v`")
     end
 
+    if !(spline_order in (1, 2, 3))
+        error("spline_order must be 1, 2, or 3 (got $spline_order)")
+    end
+    # Higher-order deposition is not yet supported on tripolar grids (multi-point seam remap is a
+    # separate task); guard rather than silently deposit at the wrong order.
+    if spline_order > 1 && typeof(grid) <: MeshGrids && occursin("Tripolar", string(nameof(typeof(grid.stats.Ny))))
+        error("spline_order > 1 is not yet supported on tripolar grids; use spline_order=1.")
+    end
+
     # initialize state {SharedArray} given grid and layers
     # Number of state variables 
     Nstate = 3
@@ -356,7 +368,8 @@ function WaveGrowth2D(; grid::GG,
         ocean_points, boundary_points,
         winds,
         currents,
-        Mstat)
+        Mstat,
+        spline_order)
 end
 
 

src/Operators/TimeSteppers.jl

@@ -202,7 +202,11 @@ function time_step!(model::Abstract2DModel, Δt::Float64; forcing=nothing, callb
         @show model.ParticleCollection[:, 6].on
     end 
 
-    time_step!_advance(model, Δt, forcing_xy, FailedCollection)
+    # Resolve the B-spline deposition order Int -> Val{P} ONCE here (function barrier). Below this
+    # point everything specializes on Val{P} and compiles once per order (cached across steps);
+    # the only per-step cost is this single dispatch. See issues #59/#60.
+    spline_val = Val(hasproperty(model, :spline_order) ? model.spline_order : 1)
+    time_step!_advance(model, Δt, forcing_xy, FailedCollection, spline_val)
     # @threads for a_particle in model.ParticleCollection[model.ocean_points]
     #     #@info a_particle.position_ij
     #     mapping_2D.advance!(    a_particle, model.State, FailedCollection,
@@ -252,7 +256,7 @@ function time_step!(model::Abstract2DModel, Δt::Float64; forcing=nothing, callb
 
 end
 
-function time_step!_advance(model::Abstract2DModel, Δt::Float64, forcing_xy, FailedCollection::Vector{AbstractMarkedParticleInstance})
+function time_step!_advance(model::Abstract2DModel, Δt::Float64, forcing_xy, FailedCollection::Vector{AbstractMarkedParticleInstance}, spline::Val{P}=Val(1)) where {P}
 
     @threads for a_particle in model.ParticleCollection[model.ocean_points]
         #@info a_particle.position_ij
@@ -269,7 +273,8 @@ function time_step!_advance(model::Abstract2DModel, Δt::Float64, forcing_xy, Fa
                 model.ODEsettings.log_energy_maximum,
                 model.ODEsettings.wind_min_squared,
                 model.periodic_boundary,
-                model.ODEdefaults)
+                model.ODEdefaults,
+                spline)
 
         # if (a_particle.position_ij[2] == 6) & (particle_on != a_particle.on)            
         #     @info "after advance! outside: $(a_particle.position_ij) particle on change :$(particle_on) to $(a_particle.on)"
@@ -327,7 +332,9 @@ function movie_time_step!(model::Abstract2DModel, Δt; forcing=nothing, callback
     current_forcing === nothing && error("2D movie_time_step! requires forcing data (pass `forcing=...` to Simulation or set `model.winds`)")
     forcing_xy = time_slice_forcing(current_forcing, model.grid, model.clock.time)
 
-    time_step!_advance(model, Δt, forcing_xy, FailedCollection)
+    # resolve B-spline deposition order Int -> Val{P} once (function barrier); see time_step!
+    spline_val = Val(hasproperty(model, :spline_order) ? model.spline_order : 1)
+    time_step!_advance(model, Δt, forcing_xy, FailedCollection, spline_val)
 
     model.MovieState = copy(model.State)
 

src/Operators/mapping_2D.jl

@@ -52,8 +52,10 @@ S       Shared array where particles are stored
 G       (TwoDGrid) Grid that defines the nodepositions
 """
 
-function ParticleToNode!(PI::AbstractParticleInstance, S::StateTypeL1, G::TwoDGrid, periodic_boundary::Bool)
-
+function ParticleToNode!(PI::AbstractParticleInstance, S::StateTypeL1, G::TwoDGrid, periodic_boundary::Bool, spline::Val=Val(1))
+        # NOTE: TwoDGrid is deprecated (issue #43); higher-order B-spline deposition is only
+        # supported on the live MeshGrids/CartesianGrid path. `spline` is accepted for dispatch
+        # compatibility but ignored here — this path always deposits at order 1 (CIC).
         #u[4], u[5] are the x and y positions of the particle
         #index_positions, weights = PIC.compute_weights_and_index(G, PI.ODEIntegrator.u[4], PI.ODEIntegrator.u[5])
         weights_and_index = PIC.compute_weights_and_index_mininal(G, PI.ODEIntegrator.u[4], PI.ODEIntegrator.u[5])
@@ -68,10 +70,10 @@ function ParticleToNode!(PI::AbstractParticleInstance, S::StateTypeL1, G::TwoDGr
         nothing
 end
 
-function ParticleToNode!(PI::AbstractParticleInstance, S::StateTypeL1, G::MeshGrids, periodic_boundary::Bool)
-        
+function ParticleToNode!(PI::AbstractParticleInstance, S::StateTypeL1, G::MeshGrids, periodic_boundary::Bool, spline::Val{P}=Val(1)) where {P}
+
         #u[4], u[5] are the x and y positions of the particle. For the CartesianGrid2D these are cooridnates relative to the particle node
-        weights_and_index = PIC.compute_weights_and_index_mininal(PI.position_ij, PI.ODEIntegrator.u[4], PI.ODEIntegrator.u[5])
+        weights_and_index = PIC.compute_weights_and_index_mininal(PI.position_ij, PI.ODEIntegrator.u[4], PI.ODEIntegrator.u[5], spline)
         # @info PI.position_ij, weights_and_index
 
         #ui[1:2] .= PI.position_xy
@@ -185,9 +187,10 @@ function advance!(PI::AbstractParticleInstance,
                         DT::Float64, 
                         log_energy_maximum::Float64,
                         wind_min_squared::Float64,
-                        periodic_boundary::Bool, 
+                        periodic_boundary::Bool,
                         default_particle::PP,
-                        ) where {PP<:Union{ParticleDefaults,Nothing},FF<:Union{ForcingCollection,ForcingData,NamedTuple{(:u, :v)},Tuple{Float64,Float64}}}
+                        spline::Val{P}=Val(1),
+                        ) where {PP<:Union{ParticleDefaults,Nothing},FF<:Union{ForcingCollection,ForcingData,NamedTuple{(:u, :v)},Tuple{Float64,Float64}},P}
         #@show PI.position_ij
 
         #if ~PI.boundary # if point is not a 
@@ -280,8 +283,8 @@ function advance!(PI::AbstractParticleInstance,
         end
 
         #if PI.ODEIntegrator.u[1] > -13.0 #ODEs.log_energy_minimum # the minimum enerçy is distributed to 4 neighbouring particles
-        if PI.on 
-                ParticleToNode!(PI, S, Grid, periodic_boundary)
+        if PI.on
+                ParticleToNode!(PI, S, Grid, periodic_boundary, spline)
         end
 
         return PI

src/ParticleInCell.jl

@@ -88,9 +88,65 @@ function get_relative_i_and_w(zp_normed::Float64)
 end
 
 
+# ---------------------- B-spline deposition kernels ------------------
+# Generalize the linear (CIC, order 1) weighting above to arbitrary B-spline order P = 1, 2, 3,
+# selectable at model setup, to reduce grid-imprint anisotropy (see issues #59, #60).
+#
+# Centering convention (standard PIC, required for partition of unity Σw = 1):
+#   - odd order  (P = 1, 3): stencil centered on floor(z)
+#   - even order (P = 2):    stencil centered on round(z) (nearest node)
+# Stencil width is N = P + 1. Weights sum to 1 at every order, so the additive deposit conserves
+# total charge/energy exactly regardless of P.
+#
+# Two entry points mirror the linear functions above and keep P = 1 byte-identical:
+#   - get_absolute_i_and_w(z, ::Val{P})           absolute indices, δ = z - center (no rounding)
+#   - get_absolute_i_and_w(z, i_node, ::Val{P})   indices relative to particle node, δ rounded to 6 digits
+
+@inline _bspline_center(z::Float64, ::Val{1}) = floor(z)
+@inline _bspline_center(z::Float64, ::Val{2}) = round(z)
+@inline _bspline_center(z::Float64, ::Val{3}) = floor(z)
+
+@inline _bspline_offsets(::Val{1}) = SVector{2,Int64}(0, 1)
+@inline _bspline_offsets(::Val{2}) = SVector{3,Int64}(-1, 0, 1)
+@inline _bspline_offsets(::Val{3}) = SVector{4,Int64}(-1, 0, 1, 2)
+
+@inline _bspline_weights(δ::Float64, ::Val{1}) = SVector{2,Float64}(1.0 - δ, δ)
+@inline _bspline_weights(δ::Float64, ::Val{2}) =
+    SVector{3,Float64}(0.5 * (0.5 - δ)^2, 0.75 - δ^2, 0.5 * (0.5 + δ)^2)
+@inline function _bspline_weights(δ::Float64, ::Val{3})
+    ω = 1.0 / 6.0
+    return SVector{4,Float64}(ω * (1 - δ)^3, ω * (3δ^3 - 6δ^2 + 4), ω * (-3δ^3 + 3δ^2 + 3δ + 1), ω * δ^3)
+end
+
+"""
+get_absolute_i_and_w(zp_normed::Float64, spline::Val{P})
+B-spline generalization of the absolute-index variant. For P = 1 returns the same indices and
+weights as `get_absolute_i_and_w(zp_normed)`.
+"""
+@inline function get_absolute_i_and_w(zp_normed::Float64, spline::Val{P}) where {P}
+    center = _bspline_center(zp_normed, spline)
+    δ = zp_normed - center
+    i0 = Int(center) + 1
+    return i0 .+ _bspline_offsets(spline), _bspline_weights(δ, spline)
+end
+
+"""
+get_absolute_i_and_w(zp_normed::Float64, i_node::Int64, spline::Val{P})
+B-spline generalization of the particle-node-relative variant (used by MeshGrids/CartesianGrid).
+For P = 1 returns the same indices and weights as `get_absolute_i_and_w(zp_normed, i_node)`,
+including the `round(..., digits=6)` on the fractional offset.
+"""
+@inline function get_absolute_i_and_w(zp_normed::Float64, i_node::Int64, spline::Val{P}) where {P}
+    center = _bspline_center(zp_normed, spline)
+    δ = round(zp_normed - center, digits=6)
+    i0 = Int(center) + i_node
+    return i0 .+ _bspline_offsets(spline), _bspline_weights(δ, spline)
+end
+
+
 
 """
-norm_distance(xp::T, xmin::T, dx::T) 
+norm_distance(xp::T, xmin::T, dx::T)
 returns normalized distance
 """
 function norm_distance(xp::T, xmin::T, dx::T) where {T<:AbstractFloat}
@@ -156,6 +212,18 @@ function compute_weights_and_index_mininal(ij::II, xp::Float64, yp::Float64) whe
     return wni(xi, xw, yi, yw)
 end
 
+"""
+compute_weights_and_index_mininal(ij, xp, yp, spline::Val{P})
+B-spline (order P) variant of the particle-node-relative 2D wrapper. Returns a `wni{P+1}`.
+`spline = Val(1)` is bit-identical to the 3-argument method above.
+"""
+function compute_weights_and_index_mininal(ij::II, xp::Float64, yp::Float64, spline::Val{P}) where {II<:Union{Tuple{Int,Int},CartesianIndex},P}
+    xi, xw = get_absolute_i_and_w(xp, ij[1], spline)
+    yi, yw = get_absolute_i_and_w(yp, ij[2], spline)
+
+    return wni(xi, xw, yi, yw)
+end
+
 """
 compute_weights_and_index(g_pars::OneDGrid, xp::Float64 )
 returns indexes and weights for in 2D for single x point
@@ -322,7 +390,7 @@ function push_to_grid!(grid::StateTypeL1,
                             Nx::Int, Ny::Int,
                             periodic::Bool=true) where {CC<:Union{Vector{Float64},SVector{3,Float64},MVector{3,AbstractFloat}},
                                                             II<:Union{Tuple{Int,Int},SVector{2,Int64}},
-                                                            WW<:Union{Tuple{Float64,Float64},SVector{2,Float16}}}
+                                                            WW<:Tuple{Float64,Float64}}
     if periodic
         grid[ wrap_index!(index_pos[1], Nx) , wrap_index!(index_pos[2], Ny), : ] += weights[1] * weights[2] * charge
     else
@@ -345,7 +413,7 @@ function push_to_grid!(grid::StateTypeL1,
                             Nx::AbstractBoundary, 
                             Ny::AbstractBoundary) where {CC<:Union{Vector{Float64},SVector{3,Float64},MVector{3,AbstractFloat}},
                                                             II<:Union{Tuple{Int,Int},SVector{2,Int64}},
-                                                            WW<:Union{Tuple{Float64,Float64},SVector{2,Float16}}}
+                                                            WW<:Tuple{Float64,Float64}}
 
     # conditions where nothing should be returned
     if  (Nx isa N_NonPeriodic) & ~test_domain(index_pos[1], Nx.N) | # non-periodic in x and y-position is out of domain
@@ -501,10 +569,18 @@ function construct_loop(wni::FieldVector{4,SVector})
 #     return zip(idx, wtx)
 # end
 
-function construct_loop(wni::FieldVector{4,SVector})
-    idx = SVector{4,Tuple{Int,Int}}(                    (wni.xi[1], wni.yi[1]), (wni.xi[2], wni.yi[1]), (wni.xi[1], wni.yi[2]), (wni.xi[2], wni.yi[2]))
-    wtx = SVector{4,Tuple{AbstractFloat,AbstractFloat}}((wni.xw[1], wni.yw[1]), (wni.xw[2], wni.yw[1]), (wni.xw[1], wni.yw[2]), (wni.xw[2], wni.yw[2]))
-    return zip(idx, wtx)
+"""
+construct_loop(w::wni{N})
+Builds the full N×N separable stencil (N = spline order + 1) as a zip of (index, weight) pairs.
+Flatten is column-major (x fastest); for N = 2 this reproduces the original CIC ordering
+(x1,y1),(x2,y1),(x1,y2),(x2,y2) exactly. Weights are kept as the (xw, yw) 2-tuple so the
+single-point push_to_grid! arithmetic `weights[1]*weights[2]*charge` is unchanged.
+"""
+function construct_loop(w::wni{N}) where {N}
+    M = N * N
+    idx = ntuple(k -> (w.xi[(k - 1) % N + 1], w.yi[(k - 1) ÷ N + 1]), Val(M))
+    wtx = ntuple(k -> (w.xw[(k - 1) % N + 1], w.yw[(k - 1) ÷ N + 1]), Val(M))
+    return zip(SVector{M,Tuple{Int,Int}}(idx), SVector{M,Tuple{Float64,Float64}}(wtx))
 end
 
 
@@ -513,11 +589,11 @@ wrapper over FieldVector weight&index (wni),
 """
 function push_to_grid!(grid::StateTypeL1,
     charge::CC,
-    wni::FieldVector,
+    weights_and_index::wni,
     Nx::Int, Ny::Int,
     periodic::Bool=true) where CC <: Union{Vector{Float64}, SVector{3, Float64}}
     #@info "this is version D"
-    for (i, w) in construct_loop(wni)
+    for (i, w) in construct_loop(weights_and_index)
         push_to_grid!(grid, charge, i, w, Nx, Ny, periodic)
     end
 end
@@ -529,10 +605,10 @@ wrapper over FieldVector weight&index (wni),
 """
 function push_to_grid!(grid::StateTypeL1,
     charge::CC,
-    wni::FieldVector,
+    weights_and_index::wni,
     Nx::AbstractBoundary, Ny::AbstractBoundary) where {CC<:Union{Vector{Float64},SVector{3,Float64}}}
     #@info "this is version D"
-    for (i, w) in construct_loop(wni)
+    for (i, w) in construct_loop(weights_and_index)
         push_to_grid!(grid, charge, i, w, Nx, Ny)
     end
 end

src/custom_structures.jl

@@ -60,14 +60,23 @@ Base.copy(s::ParticleInstance1D) = ParticleInstance1D(s.position_ij, s.position_
 Base.copy(s::ParticleInstance2D) = ParticleInstance2D(s.position_ij, s.position_xy, s.ODEIntegrator, s.boundary, s.on)
 
 # Regridding types:
-"""Weights & Index (wni) FieldVector """
-struct wni{TI<:SVector,TF<:SVector} <: FieldVector{4,SVector}
+"""
+Weights & Index (wni) for a separable N-point-per-axis deposition stencil (N = spline order + 1).
+`xi`/`yi` are the N integer node indices per axis, `xw`/`yw` the N B-spline weights per axis.
+N is encoded in the type so `construct_loop` unrolls the N^2 stencil per order.
+"""
+struct wni{N,TI<:SVector{N,Int64},TF<:SVector{N,Float64}}
         xi::TI
         xw::TF
         yi::TI
         yw::TF
 end
 
+# Explicit outer constructor so N is inferred from the stencil width (avoids relying on
+# static-parameter solving of N through the SVector{N,...} field constraints).
+wni(xi::SVector{N,Int64}, xw::SVector{N,Float64}, yi::SVector{N,Int64}, yw::SVector{N,Float64}) where {N} =
+        wni{N,SVector{N,Int64},SVector{N,Float64}}(xi, xw, yi, yw)
+
 # Define Boundary Types:
 
 struct N_Periodic{T} <: AbstractBoundary

test/Project.toml

@@ -2,5 +2,6 @@
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
 SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"