Skip to content

RZ: conserve charge and current under binomial filtering - #7059

Open
clarkse-he wants to merge 5 commits into
BLAST-WarpX:developmentfrom
Helion-Energy:fix-rz-filter-volume-weighting
Open

RZ: conserve charge and current under binomial filtering#7059
clarkse-he wants to merge 5 commits into
BLAST-WarpX:developmentfrom
Helion-Energy:fix-rz-filter-volume-weighting

Conversation

@clarkse-he

@clarkse-he clarkse-he commented Jul 15, 2026

Copy link
Copy Markdown
Contributor

Summary

The bilinear filter smooths rho and J as densities with an index-space stencil. On the unequal cell volumes of radial geometry (RZ / RCYLINDER / RSPHERE) this does not conserve the volume integral: charge and current are created or destroyed by every filter pass, with the bias concentrated near the axis where adjacent cells differ in volume by O(1) factors.

A naive extensive-quantity variant (filter rho*V, divide V back out) is conservative but wrong in the other direction: the index-symmetric stencil moves a quarter of the first ring's (large) charge into the (small) axis cells, spiking the axis density by more than 2x in our tests.

This PR applies the binomial passes to conserved densities in flux form: each pass is the divergence of a two-point diffusive flux with face-averaged volume factors, using the same volume conventions as the inverse-volume-scaling routines. This

  • conserves the volume integral of rho and J exactly in the interior,
  • leaves constant fields untouched,
  • reduces to the standard (1/4, 1/2, 1/4) stencil where volumes are uniform (Cartesian behavior unchanged),
  • passes zero flux through the axis face by construction (no below-axis guard dependence and no azimuthal-mode parity bookkeeping), and
  • passes zero flux through non-periodic domain boundaries, so nothing leaks into guard cells that are never folded back.

WarpX::ApplyFilterJ dispatches current-density filtering to the flux form in radial geometries and to ApplyFilterMF elsewhere; ApplyFilterandSumBoundaryRho uses the flux form for rho in radial geometries. Field smoothing (gather field, fluid velocities, temperature diagnostics) intentionally keeps the plain stencil: those are not conserved densities.

Validation

On an RZ uniform-plasma deck with charge-weighted diagnostics (single filter pass on identical deposits, filtered vs unfiltered):

  • near-axis rows are unbiased after filtering (ratios 0.99-1.00; the previous behavior distorted them),
  • the filtered/unfiltered total-charge ratio is a one-time redistribution bounded below 0.5% (outermost guard-layer bookkeeping) and does not accumulate over steps,
  • the profile response is confined to smoothing genuine deposit structure.

Notes for review (draft)

  1. RZ benchmark checksums of filtered tests will shift and need resetting once the approach is agreed on.
  2. Physics tests with tight tolerances tuned against the old filter should be re-examined: on one RZ deck whose axis deposit structure the two filters smooth differently, an electron-energy adiabat metric moved from 0.6% to 2.6% median.
  3. The flux-form sweeps fill one guard layer fewer per pass than the stencil form; the outermost guard layer keeps the raw deposit values and the guard sum is clamped to the well-defined region.
  4. Face weights are arithmetic means of the point volume factors; exact annulus areas are a possible refinement.
  5. Compensation: warpx.use_filter_compensation is only consumed by the PSATD k-space filter; the real-space BilinearFilter is purely binomial, so nothing else needs changing here. If a real-space compensator is ever added, it folds directly into this design: the (-a, 1+2a, -a) sharpening pass is one more flux-form sweep with coefficient -npass/4, and conservation/constant-preservation are properties of the flux structure independent of the sweep coefficient.

🤖 Generated with Claude Code

https://claude.ai/code/session_01FyFeNNpr5jir3ygSakbZCm

The bilinear filter smooths rho and J as densities with an index-space
stencil. On the unequal cell volumes of radial geometry this does not
conserve the volume integral: charge and current are created or
destroyed by every filter pass, with the bias concentrated near the
axis where adjacent cells differ in volume by O(1) factors. (An
extensive-quantity variant -- filter rho*V, divide V back out -- is
conservative but piles the first ring's charge onto the small axis
cells instead.)

Apply the binomial passes to conserved densities in flux form: each
pass becomes the divergence of a two-point diffusive flux with
face-averaged volume factors, using the same volume conventions as the
inverse-volume-scaling routines. This
 - conserves the volume integral of rho and J exactly in the interior,
 - leaves constant fields untouched,
 - reduces to the standard (1/4, 1/2, 1/4) stencil where volumes are
   uniform (Cartesian behavior unchanged, bitwise),
 - passes zero flux through the axis face by construction (no
   below-axis guard dependence, no mode-parity bookkeeping), and
 - passes zero flux through non-periodic domain boundaries, so nothing
   leaks into guard cells that are never folded back.

The new WarpX::ApplyFilterJ dispatches current-density filtering to the
flux form in radial geometries and to ApplyFilterMF elsewhere;
ApplyFilterandSumBoundaryRho uses the flux form for rho in radial
geometries. Field smoothing (the gather field, fluid velocities,
temperature diagnostics) intentionally keeps the plain stencil: those
are not conserved densities.

Verified on an RZ uniform-plasma deck (charge-weighted diagnostics):
the near-axis rows are unbiased after filtering (the previous behavior
distorted them), the profile response is confined to smoothing genuine
deposit structure, and the filtered/unfiltered total-charge ratio is a
one-time redistribution bounded below 0.5 percent (outermost guard-layer
bookkeeping) rather than a secular loss. RZ benchmark checksums of
filtered tests will shift and need resetting; results on tests with
tight physics tolerances that were tuned against the old filter should
be re-examined (observed: an RZ electron-energy adiabat metric moved
from 0.6 to 2.6 percent median on a deck whose axis deposit structure
the two filters smooth differently).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FyFeNNpr5jir3ygSakbZCm
@clarkse-he
clarkse-he marked this pull request as ready for review July 16, 2026 00:25
clarkse added a commit to Alfvenic/WarpXC that referenced this pull request Jul 25, 2026
Brings upstream/development plus three of Helion PRs from the
integration base: BLAST-WarpX#7040 isotropic Ohm-law operators, BLAST-WarpX#7038 divB/divE
from the solver fp fields, BLAST-WarpX#7059 RZ charge/current conservation under
binomial filtering.

Assembled and tested on warpxc-base first, where the vendor branches
reach the kernels they were written against. See the follow-up commit
for the metric re-homing this fork requires.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FtmsmZbwzoVwHgPzCrxzQA
clarkse added a commit to Alfvenic/WarpXC that referenced this pull request Jul 26, 2026
The azimuthal factor, the Verboncoeur axis fraction and the on- and
off-axis cell-volume weights were written out four times: in
ApplyInverseVolumeScalingToCurrentDensity, ...ToChargeDensity,
...ToMassMatricesPC and ApplyVolumeWeightedFilter. They had already
drifted -- the filter still carried the pi*dr*avf on-axis form that was
just corrected in the other three, and the spherical Verboncoeur fraction
existed in two of the four. Move all of it into
Utils/RadialVolumeScaling.H so there is one definition to be right or
wrong, and have every site call it.

Behaviour is unchanged: the on-axis probe still reads 0.995 of a uniform
density after the refactor, and the full local gate is 22/22.

Also parse boundary.verboncoeur_axis_correction in the RTZ and 3D builds.
It was guarded to RZ/RCYLINDER/RSPHERE only, so on a runtime-cylindrical
deck the option was silently ignored and pinned to its default with no
warning. Now measurable: setting it to 0 moves the axis node by exactly
4/3, the ratio of the two fractions.

The helper makes the 3D spherical origin an explicit NaN rather than a
silent fall-through to the cylindrical weight. That case is unreachable
today (r_lo > 0 is asserted), but the r == 0 branches had no is_spherical
test at all while the r != 0 branches did, so the fall-through was there
waiting. The origin is a real coordinate singularity -- every (theta, phi)
collapses to one point -- so regularising it means treating the innermost
shell as a single degree of freedom, not picking a different constant.

Not done here: wiring the volume-weighted filter (BLAST-WarpX#7059) to the runtime
curvilinear metrics. Enabling it regressed the RTZ face-deposition volume
integrals by 4-5% low, i.e. it lost current. The suspected cause and the
work needed are recorded at the dispatch site; it deliberately still falls
back to the plain stencil.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FtmsmZbwzoVwHgPzCrxzQA
@EZoni EZoni added geometry: RZ axisymmetric 2D and quasi-3D geometry: cylinder/sphere labels Jul 28, 2026
@EZoni
EZoni requested a review from dpgrote July 28, 2026 22:47
Comment thread Source/Parallelization/WarpXComm.cpp Outdated
clarkse-he and others added 2 commits July 29, 2026 11:24
The radial direction is never periodic (enforced at setup), so zero the
domain-edge face weights unconditionally in the dir = 0 sweep.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011UeLDH36zCMcaDAfdjvQr1
Each flux-form pass consumes one defined guard layer of the working
arrays while pushing mass one layer outward. Sized to the source guard
width, npass passes left the outermost npass layers of the output
stale, and the subsequent guard-cell sum folded those stale layers,
dropping the interface flux: on a periodic-z hybrid deck filtered J_z
lost ~5e-6 of the total current per z pass at npass = 1, growing to
~3e-3 at npass = 8, and rho leaked for npass > 1 as well (the rho path
seeds and clamps, which only covers npass = 1 headroom).

Give the temporaries 2*npass extra guard layers per swept direction.
The filter runs before any guard-cell sum, so source guard layers hold
only local deposits bounded by the source guard width -- the extension
layers are genuinely zero. After npass passes the defined region ends
at src.ng + npass, covering the final mass reach, so every layer the
guard sum folds holds filtered data. Verified on ballistic hybrid RZ
twins (bit-identical particles): rho and J_z volume integrals conserved
to <= 4e-16 for radial-only, axial-only, and combined filtering at
npass = 1 and 8.

With the filter conservative in the radial direction, the HybridPIC
"use at your own risk" radial-filtering warning is obsolete; remove it.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011UeLDH36zCMcaDAfdjvQr1
@clarkse-he

Copy link
Copy Markdown
Contributor Author

Conservation evidence: flux-form filter vs. development

Setup: hybrid-PIC RZ (the solver that permits radial filter passes), 64×128, m=0,
axis-peaked Gaussian proton column (σ_r ≈ 8 Δr) with a small axial drift,
periodic z, filter applied to the fresh ρ/J deposits every step. Metric: the
code-native volume quadrature Σ w·u (2πr off axis, π Δr f_axis on axis — the
same convention as the inverse-volume-scaling routines) of the raw staggered
rho_fp/jz_fp, for each filtered run against an unfiltered twin. Deposition
conserves the twin's integral exactly, so the ratio isolates the filter.

conservation figure

  • Left — 100-step hybrid run, standard 1 pass/direction: development
    creates/destroys 2–4×10⁻³ of the total charge as a persistent systematic
    (bias concentrated near the axis); this branch sits at the measurement floor
    (≤ 2×10⁻¹⁶ early, ≤ 10⁻⁹ late, the late rise being the twin's own slow drift).
  • Right — ballistic twins (E = B = 0, so filtered/unfiltered particles stay
    bit-identical and the floor is pure roundoff): development's error grows
    linearly with pass count for both total charge and total axial current,
    reaching 2.4×10⁻² at 8 passes; this branch stays ≤ 4×10⁻¹⁶ for both
    quantities at every pass count.

The npass sweep is also what exposed the guard-budget issue fixed in
e012480: the filter's working arrays now carry 2·npass extra guard layers,
so the subsequent guard-cell sum only ever folds filtered data.

clarkse-he and others added 2 commits July 29, 2026 13:59
test_rz_laser_acceleration, test_rz_laser_acceleration_opmd, and
test_rz_ohm_solver_cylinder_compression_picmi are the RZ tests that
enable binomial filtering, and the flux-form filter intentionally
changes filtered rho and J in radial geometry. Benchmarks regenerated
at the current head (including the guard-headroom fix); all three
evaluate at zero relative error against the regenerating runs.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011UeLDH36zCMcaDAfdjvQr1
The three ion-temperature moments of the hybrid cylinder-compression
test are sensitive to the run environment (thread-dependent sampling)
at the few-1e-5 level, above the test's 1e-6 tolerance, so the locally
regenerated values do not transfer. All field and density keys matched
CI exactly; adopt CI's printed checksums for the three T keys.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011UeLDH36zCMcaDAfdjvQr1
clarkse added a commit to Alfvenic/WarpXC that referenced this pull request Jul 30, 2026
…his fork

Bisected with purpose-built RZ control builds. Clean boundary: pristine
upstream/development passes 5/5, upstream + BLAST-WarpX#7040 + BLAST-WarpX#7038 passes 5/5, and adding
BLAST-WarpX#7059 (db36e72, 'RZ: conserve charge and current under binomial filtering')
fails on rho keys at 2.14e-3 / 1.84e-3. The fork reproduces the base's numbers
exactly, so none of its ~49 commits nor vendor rounds 2-3 contribute.

Not a defect -- changing RZ filtered charge is that PR's purpose. It simply
re-blessed no benchmarks (zero benchmark files touched), so the shift shows up as
two red RZ checksums.

Corrects my own framing: I had recorded 1b as 'something on this FORK moved RZ
charge DEPOSITION'. Both halves were wrong -- it is the vendor stack, and it is
the filter rather than deposition. Three suspects I had nominated are eliminated
by measurement in the entry, including finding 16, whose IEEE-exactness claim
holds bitwise.

Leaves one thing unexplained rather than papering over it: why a z-only filter
sweep changes rho when BLAST-WarpX#7059 claims the volume-weighted filter reduces to the
standard stencil where volumes are uniform, as they are in z.

Flags the real action: BLAST-WarpX#7059 is Eric's own open upstream PR, so upstream CI
should fail these two tests on review. Better fixed there than re-blessed here;
our baselines stay un-blessed pending his call.
@clarkse-he
clarkse-he requested a review from dpgrote July 31, 2026 00:19
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

geometry: cylinder/sphere geometry: RZ axisymmetric 2D and quasi-3D

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants