Fix charge conservation in charge exchange and add automated test by oshapoval · Pull Request #6796 · BLAST-WarpX/warpx

oshapoval · 2026-04-21T14:27:52Z

Description updated by @RemiLehe

Motivation

We recently observed that combining charge exchange with the electrostatic solver (in the development branch) leads to unphysical heating (i.e. energy conservation breaks). After investigation, it seems that this is because the charge "jumps" from one macroparticle to another when a charge exchange event happens.

More specifically, consider the reaction
$$A^+ + B \rightarrow A + B^+$$
between one A+ macroparticle and B macroparticle that are in the same cell but not necessarily at the exact same location. By default (in the development branch), B+ is created at the position of B and A is created at the position of A+. But that means that the charge suddenly moves from one position to another in one timestep. When recomputing the electrostatic fields at the following timestep (with the electrostatic solver), the field configuration is thus substantially different. It seems that this provides a sudden "kick" to the particles and is what leads to a lack of energy conservation.

The solution implemented in this PR is to create the B+ particle at the position of A+ (and B at the position of A). This appears to fix energy conservation, for a uniform thermal plasma (see more details below):

What this PR does

Fix: Reintroduces a dedicated CHARGE_EXCHANGE reaction type (previously conflated with TWO_PRODUCT_REACTION) so that, upon a charge exchange event, the positions of the two product particles are swapped. This ensures that charge does not jump across space in a single timestep, restoring charge conservation.
New test: Adds a test analogous to test_2d_collisions_split_momentum_push_electrostatic (which verifies energy conservation when Coulomb scattering is combined with the electrostatic solver). A new input script is introduced that replaces Coulomb scattering with charge exchange and uses three species — neutrals, ions, and neutralizing electrons — instead of two. The existing analysis script (which checks energy conservation and energy equipartition) is reused but generalized to remove its hardcoded two-species assumption. The test is designed to fail if the position-swap fix is not applied.

…harge exchange CI test

…hange_2__ci_analysis

…change CI test

oshapoval · 2026-04-25T00:29:25Z

Energy conservation checks:

Without charge exchange fix (without swapping the positions):

With charge exchange fix (with swapping the positions):

RemiLehe · 2026-04-25T01:35:16Z

    # verify the total energy conservation
    total_energy_error_norm = np.max(np.abs(total_energy_error))
-    relative_tolerance = 1e-5 if electrostatic else 6e-5
+    relative_tolerance = 1.1e-5 if electrostatic else 6e-5


The new test requires the tolerance to be slightly relaxed.

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

for more information, see https://pre-commit.ci

RemiLehe · 2026-04-26T14:44:04Z


    # compare the final value of the field energy variation to the reference equipartition value
-    relative_tolerance = 1e-1 if electrostatic else 5e-1
+    relative_tolerance = 1.1e-1 if electrostatic else 5e-1


For the new test to pass, the tolerance needs to be very slightly increased.

roelof-groenewald · 2026-04-28T20:11:37Z

+                // Position components occupy indices 0 through PIdx::w-1 in m_rdata.
+                if (mask[i] == int(ScatteringProcessType::CHARGE_EXCHANGE)) {
+                    for (int idim = 0; idim < PIdx::w; ++idim) {
+                        soa_products_data[2].m_rdata[idim][slot2_idx] = soa_1.m_rdata[idim][src1_idx];
+                        soa_products_data[3].m_rdata[idim][slot3_idx] = soa_0.m_rdata[idim][src0_idx];
+                    }
+                }
+


What do you guys think about making this more future proof by using the GetPositionAsStored and SetPositionAsStored structs rather than relying on the position attributes being first in the order?

In retrospect, since the particle iterators are not themselves easily accessible from this location, I think reverting to this original approach is better (more efficient and more readable). Sorry for the flip-flopping!

Oh man, @oshapoval, I'm so sorry but I just realized there might be a problem with this approach due to the fact that theta is actually not before w in the attribute list for RZ, RCYLINDER, and RSPHERE - see

warpx/Source/Particles/WarpXParticleContainer.H

Line 74 in fcd1be2

static constexpr auto names = {"r", "z", "w", "ux", "uy", "uz", "theta"};

.

I think the balance of things then is to just have some code duplication here that will keep things easy to read. So making this something like:

soa_products_data[2].m_rdata[PIdx::z][slot2_idx] = soa_1.m_rdata[PIdx::z][src1_idx]; soa_products_data[3].m_rdata[PIdx::z][slot3_idx] = soa_0.m_rdata[PIdx::z][src0_idx]; #if defined(WARPX_DIM_RZ) soa_products_data[2].m_rdata[PIdx::x][slot2_idx] = soa_1.m_rdata[PIdx::x][src1_idx]; soa_products_data[2].m_rdata[PIdx::theta][slot2_idx] = soa_1.m_rdata[PIdx::theta][src1_idx]; soa_products_data[3].m_rdata[PIdx::x][slot3_idx] = soa_0.m_rdata[PIdx::x][src0_idx]; soa_products_data[3].m_rdata[PIdx::theta][slot3_idx] = soa_0.m_rdata[PIdx::theta][src0_idx]; #elif defined(WARPX_DIM_RCYLINDER) soa_products_data[2].m_rdata[PIdx::theta][slot2_idx] = soa_1.m_rdata[PIdx::theta][src1_idx]; soa_products_data[3].m_rdata[PIdx::theta][slot3_idx] = soa_0.m_rdata[PIdx::theta][src0_idx]; #elif defined(WARPX_DIM_RSPHERE) soa_products_data[2].m_rdata[PIdx::theta][slot2_idx] = soa_1.m_rdata[PIdx::theta][src1_idx]; soa_products_data[3].m_rdata[PIdx::theta][slot3_idx] = soa_0.m_rdata[PIdx::theta][src0_idx]; soa_products_data[2].m_rdata[PIdx::phi][slot2_idx] = soa_1.m_rdata[PIdx::phi][src1_idx]; soa_products_data[3].m_rdata[PIdx::phi][slot3_idx] = soa_0.m_rdata[PIdx::phi][src0_idx]; #elif defined(WARPX_DIM_3D) soa_products_data[2].m_rdata[PIdx::x][slot2_idx] = soa_1.m_rdata[PIdx::x][src1_idx]; soa_products_data[2].m_rdata[PIdx::y][slot2_idx] = soa_1.m_rdata[PIdx::y][src1_idx]; soa_products_data[3].m_rdata[PIdx::x][slot3_idx] = soa_0.m_rdata[PIdx::x][src0_idx]; soa_products_data[3].m_rdata[PIdx::y][slot3_idx] = soa_0.m_rdata[PIdx::y][src0_idx]; #elif defined(WARPX_DIM_XZ) soa_products_data[2].m_rdata[PIdx::x][slot2_idx] = soa_1.m_rdata[PIdx::x][src1_idx]; soa_products_data[3].m_rdata[PIdx::x][slot3_idx] = soa_0.m_rdata[PIdx::x][src0_idx]; #endif

Thoughts?

Maybe I don't know what is happening here, but could this be done using SmartCopy?

@roelof-groenewald Good point, thanks for catching this.

@dpgrote SmartCopy is being used here to create two new particles (soa_products_data[2] particle at slot2_idx based on soa_0 parent particle at src0_idx and soa_products_data[3] particle at slot3_idx based on parent particle sao_1 particle at src1_idx). However, in order to keep charge conservation during the charge exchange process the location of the two new particles is being swapped so that their locations depend on the other one's parent particle.

The previous `for (int idim = 0; idim < PIdx::w; ++idim)` loop assumed position attributes were stored at indices 0..PIdx::w-1, which is true for 1D_Z, XZ, and 3D, but NOT for RZ, RCYLINDER, and RSPHERE — there theta and phi come after `w` in PIdx. Replace with a range-based for over a `static constexpr std::array<int, N>` whose contents are selected per-dimension via `#if defined(WARPX_DIM_*)`. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

A file-scope `static constexpr std::array` has no device address, so NVCC fails with "identifier 'position_indices' is undefined in device code" when the range-based for loop tries to call begin()/end() on it. Declaring the array as a local `constexpr` inside the device lambda bakes the values into the kernel as compile-time constants and works on both CPU and GPU. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

roelof-groenewald

Thanks! Looks good now!

RemiLehe and others added 8 commits April 3, 2026 07:38

Add energy conservation test with charge exchange

3ea212f

Add charge exchange type ; switch positions in that case

d9600be

Generalize analysis script to support energy conservation check for c…

58fd6aa

…harge exchange CI test

Add the corresponding checksums for test_2d_charge_exchange_dsmc CI

6592e8d

Merge remote-tracking branch 'origin/development' into fix_charge_exc…

23c92ad

…hange_2__ci_analysis

Update checksums

d80af6d

Merge remote-tracking branch 'origin/development' into fix_charge_exc…

aba1fd7

…hange_2__ci_analysis

Generalize equipartition formula and add back its check for charge ex…

59b3358

…change CI test

RemiLehe mentioned this pull request Apr 25, 2026

[WIP] Make DSMC charge exchange locally charge conserving #6744

Closed

RemiLehe requested a review from roelof-groenewald April 25, 2026 01:08

RemiLehe changed the title ~~Fix charge exchange & generalizes the analysis script~~ Fix charge conservation in charge exchange and add automated test Apr 25, 2026

RemiLehe changed the title ~~Fix charge conservation in charge exchange and add automated test~~ [WIP] Fix charge conservation in charge exchange and add automated test Apr 25, 2026

RemiLehe assigned roelof-groenewald Apr 25, 2026

Update tolerance

1c80c0d

RemiLehe reviewed Apr 25, 2026

View reviewed changes

Comment thread Examples/Tests/collision/inputs_test_2d_charge_exchange_dsmc Outdated

RemiLehe changed the title ~~[WIP] Fix charge conservation in charge exchange and add automated test~~ Fix charge conservation in charge exchange and add automated test Apr 25, 2026

RemiLehe reviewed Apr 25, 2026

View reviewed changes

Comment thread Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H Outdated

RemiLehe and others added 3 commits April 25, 2026 10:52

Apply suggestions from code review

acd9cb2

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

[pre-commit.ci] auto fixes from pre-commit.com hooks

d7eecf3

for more information, see https://pre-commit.ci

Modify tolerance

2e3ad1d

RemiLehe reviewed Apr 26, 2026

View reviewed changes

roelof-groenewald reviewed Apr 28, 2026

View reviewed changes

Comment thread Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H

roelof-groenewald reviewed Apr 28, 2026

View reviewed changes

Address review suggestion: change bitwise to logical AND

39fcf1d

oshapoval force-pushed the fix_charge_exchange_2__ci_analysis branch from 4f6ff4e to 39fcf1d Compare April 29, 2026 22:44

RemiLehe force-pushed the fix_charge_exchange_2__ci_analysis branch from 7ccf3c6 to 31eaa1c Compare May 1, 2026 14:39

Fix const-ness

0f9b6f9

roelof-groenewald approved these changes May 1, 2026

View reviewed changes

roelof-groenewald enabled auto-merge (squash) May 1, 2026 20:42

roelof-groenewald merged commit 66006c9 into BLAST-WarpX:development May 1, 2026
47 of 48 checks passed

EZoni added bug Something isn't working component: electrostatic electrostatic solver bug: affects latest release Bug also exists in latest release version labels May 1, 2026

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix charge conservation in charge exchange and add automated test#6796

Fix charge conservation in charge exchange and add automated test#6796
roelof-groenewald merged 16 commits intoBLAST-WarpX:developmentfrom
oshapoval:fix_charge_exchange_2__ci_analysis

oshapoval commented Apr 21, 2026 •

edited by RemiLehe

Loading

Uh oh!

oshapoval commented Apr 25, 2026

Uh oh!

RemiLehe Apr 25, 2026

Uh oh!

Uh oh!

Uh oh!

RemiLehe Apr 26, 2026

Uh oh!

Uh oh!

roelof-groenewald Apr 28, 2026 •

edited

Loading

Uh oh!

roelof-groenewald Apr 29, 2026

Uh oh!

roelof-groenewald Apr 30, 2026

Uh oh!

dpgrote Apr 30, 2026

Uh oh!

oshapoval Apr 30, 2026

Uh oh!

roelof-groenewald Apr 30, 2026 •

edited

Loading

Uh oh!

roelof-groenewald left a comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

Conversation

oshapoval commented Apr 21, 2026 • edited by RemiLehe Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Motivation

What this PR does

Uh oh!

oshapoval commented Apr 25, 2026

Uh oh!

RemiLehe Apr 25, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

RemiLehe Apr 26, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

roelof-groenewald Apr 28, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

roelof-groenewald Apr 29, 2026

Choose a reason for hiding this comment

Uh oh!

roelof-groenewald Apr 30, 2026

Choose a reason for hiding this comment

Uh oh!

dpgrote Apr 30, 2026

Choose a reason for hiding this comment

Uh oh!

oshapoval Apr 30, 2026

Choose a reason for hiding this comment

Uh oh!

roelof-groenewald Apr 30, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

roelof-groenewald left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

oshapoval commented Apr 21, 2026 •

edited by RemiLehe

Loading

roelof-groenewald Apr 28, 2026 •

edited

Loading

roelof-groenewald Apr 30, 2026 •

edited

Loading