split_momentum_push silently enabled by default for EB runs after #6617, breaking energy conservation with SIPIC solver

[anecas1979-cyber]

Summary

Commit b603cbb (PR #6617) changed the default behavior of split_momentum_push for simulations using embedded boundaries (EB), causing catastrophic energy non-conservation when running the semi-implicit electrostatic solver (warpx_effective_potential=True) with DSMCCollisions.

Affected configuration

• Solver: ElectrostaticSolver, method='Multigrid', warpx_effective_potential=True (SIPIC)
• Collisions: DSMCCollisions
• Geometry: 3D Cartesian with embedded boundary (EB::enabled() == true)
• Evolve scheme: Explicit
• Build: 3D, CUDA, double precision, MPI

Root cause

In Source/WarpX.cpp::ReadParameters(), PR #6617 changed:

// Before (working)
if (evolve_scheme == EvolveScheme::Explicit && !EB::enabled()) {
    m_collisions_split_momentum_push = true;
}

// After (broken)
if (evolve_scheme == EvolveScheme::Explicit) {
    m_collisions_split_momentum_push = true;
}

Since warpx_effective_potential uses EvolveScheme::Explicit, the new condition evaluates to true for SIPIC+EB runs. Before this commit, !EB::enabled() blocked split_momentum_push from activating. After this commit, it is enabled by default for all explicit+EB simulations with collisions, without any user-facing warning.

The mid-v push is not compatible with the SIPIC effective potential solver: splitting the Boris push and scattering particle velocities mid-step via DSMC breaks the energy consistency that the semi-implicit field solve relies upon.

Bisection

Commit	Description	Energy conserved?
58f84cf	Replace WARPX_PROFILE with ABLASTR_PROFILE	✅ Yes
b603cbb	Enable mid-v push with embedded boundaries (#6617)	❌ No

Observed symptoms

1. Silent energy drain (0–15 µs): total energy decays ~4× faster than in the working commit.

2. Catastrophic energy spike at ~18 µs: total energy increases from ~4.5×10⁴ J to ~9.2×10⁵ J (20× the initial value), predominantly in electrons.

Part_energy.txt — working commit (58f84cf):

step    time(s)           total(J)          ions(J)           electrons(J)
0       0.000e+00         4.7603e+04        2.3983e+04        2.3620e+04
7490    7.490e-06         4.7243e+04        2.3928e+04        2.3314e+04
14980   1.498e-05         4.7105e+04        2.3821e+04        2.3283e+04
20972   2.097e-05         4.7020e+04        2.3752e+04        2.3267e+04

Part_energy.txt — broken commit (b603cbb):

step    time(s)           total(J)          ions(J)           electrons(J)
0       0.000e+00         4.7601e+04        2.3982e+04        2.3619e+04
7490    7.490e-06         4.6269e+04        2.3935e+04        2.2334e+04
14980   1.498e-05         4.5437e+04        2.3742e+04        2.1693e+04
17227   1.723e-05         4.7535e+04        2.3866e+04        2.3668e+04  <- sudden jump
17976   1.798e-05         9.1887e+05        1.7563e+05        7.4323e+05  <- catastrophic spike

Throughput also drops from 1.65 steps/s (58f84cf) to 1.43 steps/s (b603cbb), consistent with the extra particle push kernel being active.

Workaround

Add to simulation inputs before WarpX initialization:

sys.argv += ['collisions.split_momentum_push=0']

Suggested fix

Restore the !EB::enabled() guard until split_momentum_push is validated for EB, or add an explicit check blocking split_momentum_push when warpx_effective_potential=True:

if (evolve_scheme == EvolveScheme::Explicit && !warpx_effective_potential) {
    m_collisions_split_momentum_push = true;
}

Alternatively, emit a warning and force split_momentum_push=false when both EB and warpx_effective_potential are active.

Reported by

Ales Necas

[RemiLehe]

Thanks for this bug report. While your data clearly supports the fact that energy conservation is indeed broken, I am not sure I understand why.

The above text mentions that "The mid-v push is not compatible with the SIPIC effective potential solver: splitting the Boris push and scattering particle velocities mid-step via DSMC breaks the energy consistency that the semi-implicit field solve relies upon." However, I am not sure I agree with this statement. More specifically, I am not sure that the effective potential solver has any built-in energy conservation properties, esp. ones that requires the Boris push to be unsplit.

For more context, here is the paper that shows that the mid-v push has better energy-conserving properties for the explicit electrostatic solver:
https://journals.aps.org/pre/abstract/10.1103/PhysRevE.111.025306
While I agree that this paper only applies to the explicit electrostatic solver, I would be interested in understanding which specific part of the derivation fails for the effective potential solver, and whether one could explain, in mathematical terms, why the energy conservation is worse with the mid-v push, for the effective potential solver.

[anecas1979-cyber]

Thank you so much.

I agree there are many moving parts: SIPIC + DSMC + EB + mid-v push
Looking at inputs_test_2d_background_mcc_picmi.py — the CI test modified by PR #6617 does not include an embedded boundary. So the passing CI tests confirm only that mid-v push + MCC without EB is correct. The combination of mid-v push + any collision type + EB is untested. This suggests the regression may not be SIPIC or DSMC specific — it could affect any simulation using collisions + EB + mid-v push.

[RemiLehe]

Ah yes, good point.
To confirm this, in your use case, would you be able to run with the explicit electrostatic solver instead of the effective potential solver?

[anecas1979-cyber]

Thank you, I took the spacecraft charging example and running ES + explicit + DSMC + ionization + EB + RZ, so far I was unable to reproduce the spike. Running with WarpX (26.02-49-gb603cbb135e7).

[anecas1979-cyber]

I wrote a reproducer (below) that failed so far to reproduce. I will go back to the complex 3D run with FRC + SIPIC + EB + DSMC + biasing + ionization and see if I can track what could cause the spike:

#!/usr/bin/env python3

--- Minimal reproducer for WarpX issue #6747:

--- split_momentum_push silently enabled for EB runs after commit b603cbb (PR #6617),

--- causing energy non-conservation with DSMC ionization + EB + ElectrostaticSolver.

---

--- Components: ES Multigrid + SIPIC scheme + EB sphere + electrons + protons

--- + neutral hydrogen + DSMCCollisions (ionization)

---

--- split_momentum_push activates automatically (default=true for explicit+collisions

--- after b603cbb). Ionization creates newborn particles that skip the first

--- half-push but receive the second, injecting spurious energy each ionization event.

---

--- To confirm workaround add before sim.step():

--- import sys; sys.argv += ['collisions.split_momentum_push=0']

---

--- Related: #6747

import os

GPU assignment MUST happen before any MPI/WarpX import

slurm_id = os.environ.get('SLURM_LOCALID')
pbs_id = os.environ.get('PMI_LOCAL_RANK')

if slurm_id is not None:
ngpus = int(os.environ.get('SLURM_GPUS_PER_NODE', 4))
os.environ['CUDA_VISIBLE_DEVICES'] = str(ngpus - 1 - int(slurm_id))
elif pbs_id is not None:
ngpus = int(os.environ.get('PBS_NGPUS', 4))
os.environ['CUDA_VISIBLE_DEVICES'] = str(ngpus - 1 - int(pbs_id))

import numpy as np
from pywarpx import picmi
from mpi4py import MPI as mpi

##########################

parameters

##########################

dt = 1.27e-8
max_steps = 5000
diag_period = 10

nr, nz = 40, 80
rmin, rmax = 0.0, 3.0
zmin, zmax = -3.0, 3.0

n_plasma = 7.0e9 # m^-3 electrons + protons
n_neutrals = 1.0e19 # m^-3 neutral hydrogen — high enough for ionizations
Te = 85.0 # eV
Ti = 0.05 * Te
Tn = 0.1 # eV cold neutrals

qe = picmi.constants.q_e
m_e = picmi.constants.m_e
m_i = 1836.0 * m_e

v_eth = (qe * Te / m_e) ** 0.5
v_ith = (qe * Ti / m_i) ** 0.5
v_nth = (qe * Tn / m_i) ** 0.5

cross_sec_dir = "/home/anecas1979/src/koios-data-main/XS/"

##########################

species

##########################

electrons = picmi.Species(
particle_type="electron",
name="electrons",
initial_distribution=picmi.UniformDistribution(
density=n_plasma,
rms_velocity=[v_eth, v_eth, v_eth],
),
)

protons = picmi.Species(
particle_type="proton",
name="protons",
initial_distribution=picmi.UniformDistribution(
density=n_plasma,
rms_velocity=[v_ith, v_ith, v_ith],
),
)

neutrals = picmi.Species(
name="neutrals",
charge=0,
mass=m_i,
initial_distribution=picmi.UniformDistribution(
density=n_neutrals,
rms_velocity=[v_nth, v_nth, v_nth],
),
warpx_do_not_deposit=1,
warpx_do_not_gather=1,
)

##########################

collisions

##########################

dsmc_colls = picmi.DSMCCollisions(
name="coll_elec_dsmc",
species=[electrons, neutrals],
product_species=[electrons, protons],
ndt=1,
scattering_processes={
"ionization": {
"cross_section": cross_sec_dir + "hydrogen_electron_impact_ionization.dat",
"energy": 13.598,
"target_species": neutrals,
}
},
)

##########################

grid, solver, EB

##########################

grid = picmi.CylindricalGrid(
number_of_cells=[nr, nz],
n_azimuthal_modes=1,
lower_bound=[rmin, zmin],
upper_bound=[rmax, zmax],
lower_boundary_conditions=["none", "dirichlet"],
upper_boundary_conditions=["dirichlet", "dirichlet"],
lower_boundary_conditions_particles=["none", "absorbing"],
upper_boundary_conditions_particles=["absorbing", "absorbing"],
)

solver = picmi.ElectrostaticSolver(
grid=grid,
method="Multigrid",
warpx_absolute_tolerance=1e-5,
warpx_effective_potential=True,
warpx_effective_potential_factor=6
)

embedded_boundary = picmi.EmbeddedBoundary(
implicit_function="-(x2+y2+z2-radius2)",
potential=0.0,
radius=0.3277,
)

##########################

diagnostics

##########################

energy_diag = picmi.ReducedDiagnostic(
diag_type="ParticleEnergy",
name="part_energy",
period=diag_period,
)

##########################

simulation

##########################

sim = picmi.Simulation(
solver=solver,
time_step_size=dt,
max_steps=max_steps,
warpx_embedded_boundary=embedded_boundary,
warpx_random_seed=1,
)

layout = picmi.PseudoRandomLayout(n_macroparticles_per_cell=5, grid=grid)

sim.add_species(electrons, layout=layout)
sim.add_species(protons, layout=layout)
sim.add_species(neutrals, layout=layout)

sim.add_diagnostic(energy_diag)

sim.collisions = [dsmc_colls]

sim.step(max_steps)