Generalize PIC deposition to arbitrary B-spline order (P=1,2,3)#62
Open
mochell wants to merge 3 commits into
Open
Generalize PIC deposition to arbitrary B-spline order (P=1,2,3)#62mochell wants to merge 3 commits into
mochell wants to merge 3 commits into
Conversation
Replace the hard-coded linear (CIC) particle->grid weighting with a selectable B-spline order P=1,2,3, set via `spline_order` at model setup, to reduce the grid-axis energy focusing / axisymmetry break documented in issues #59 and #60. - wni becomes parametric wni{N} (N=P+1 points/axis); explicit constructor. - Add closed-form bspline_1d Val{P} kernels (CIC/TSC/cubic), odd orders centered on floor, even on round; partition of unity => exact conservation. P=1 is bit-identical to the previous linear kernel (regression guard). - construct_loop generalized to the N^2 stencil via ntuple+Val. - Thread the order from WaveGrowth2D.spline_order through a Val function barrier in time_step! -> advance! -> ParticleToNode! (compiles once per order; no per-step recompile). WaveGrowth1D forwards it. - Guard spline_order>1 on tripolar grids (multi-point seam remap is a separate task). Tests: new test_pic_bspline_deposition.jl (partition of unity, P=1 bit identity); P=1,2,3 sweeps added to the 2D energy, 1D charge (redirected to the live 2D deposit), and 1D wave-growth periodic tests (re-enabled). Benchmark for (P+1)^2 cost scaling. Full suite green. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Collect per-deposit time and allocations for P=1,2,3 and save a 2-panel scaling plot (time + allocations vs (P+1)^2 stencil points, with ideal proportional reference) to plots/bspline_deposition_scaling.png. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
movie_time_step! (the entry point used by the comparison harness and movie
output) called time_step!_advance without the Val{P} spline barrier, so it
silently deposited at P=1 regardless of model.spline_order. Resolve the order
once at the barrier here too, mirroring time_step!.
Add a regression test asserting movie_time_step! produces order-dependent
fields (P=1 vs P=3 differ) while conserving total energy.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #60
Closes #59
What
Replaces the hard-coded linear (CIC, order 1) particle→grid weighting with a selectable B-spline order
P = 1, 2, 3, set via aspline_orderkeyword at model setup. Higher orders (TSC, cubic) are radially smoother, which removes the Cartesian grid-axis energy focusing that breaks axisymmetry (#59).Why
The advect → deposit-to-4-nodes → remesh cycle discretizes axis-aligned propagation differently from diagonal propagation, funnelling energy onto the grid axes (#59). The deposit is purely additive and B-spline weights satisfy partition of unity (Σw = 1) at any order, so total energy/momentum are conserved exactly regardless of
P— the fix is a smoother kernel, not a tunable hack.Changes
wnibecomes parametricwni{N}(N = P+1 points/axis) with an explicit constructor.bspline_1dkernels viaVal{P}(CIC/TSC/cubic); odd orders centered onfloor, even onround. P=1 is bit-identical to the previous linear kernel (regression guard).construct_loopgeneralized to the N² stencil viantuple+Val(type-stable, verified with@code_warntype).WaveGrowth2D.spline_orderthrough aValfunction barrier intime_step!→advance!→ParticleToNode!— compiles once per order (cached), no per-step recompile, no per-particle dispatch.WaveGrowth1Dforwards it.spline_order > 1on tripolar grids (multi-point seam remap deferred to Support higher-order (P>1) B-spline PIC deposition on tripolar grids #61).Tests & benchmark
test_pic_bspline_deposition.jl: partition of unity + P=1 bit-identity.runtests.jl).(P+1)²cost scaling with no per-order overhead.Result (pure-propagation ring, dx=10 km, 30 h)
The grid-axis spike ratio (cardinal/diagonal) drops from 11.17× (P=1) → 1.05× (P=3) with identical energy conservation (188 → 376) — the #59 cross pattern is removed at higher order.
Follow-up
🤖 Generated with Claude Code