Skip to content

Conversation

@n01r
Copy link
Member

@n01r n01r commented Nov 14, 2025

Triggered by a badly formed momentum_expression, I encountered a segmentation fault that did not leave a message.
Compiling again with debug symbols, an abort message appeared that pointed to https://github.com/AMReX-Codes/amrex/blob/15c4a2c360549fdcfac7fd02c328efa44e476d59/Src/Base/AMReX_Array4.H#L46.

This PR implements a check for excessive particle momenta before particles are actually initialized.

It implements CheckInitialParticleVelocity() to detect when particles might have excessive velocities that cause them to cross multiple cells per timestep, possibly leading to physically incorrect simulation results.

The check runs once during species initialization, evaluating momentum functions at all cell centers where particles will be created. It issues warnings at three severity levels based on displacement per timestep: >1, >sqrt(2), or >sqrt(3) cells.

The check converts momentum (u=gamma*v/c) to velocity and accounts for moving window velocities. Users can disable it with <species>.do_initial_velocity_check = 0 for special cases (e.g., passive tracer particles).

Updated documentation in parameters.rst and Python PICMI bindings.

n01r and others added 3 commits November 14, 2025 13:53
Implement CheckInitialParticleVelocity() to detect when particles
have excessive velocities that cause them to cross multiple cells
per timestep, leading to physically incorrect simulation results.

The check runs once during species initialization, evaluating
momentum functions at all cell centers where particles will be
created. It issues warnings at three severity levels based on
displacement per timestep: >1, >sqrt(2), or >sqrt(3) cells.

The check correctly converts momentum (u=gamma*v/c) to velocity
and accounts for moving window velocities. Users can disable it
with <species>.do_initial_velocity_check = 0 for special cases.

Updated documentation in parameters.rst and Python PICMI bindings.
Add `ignore_unused(sy)` so unused-variable-warnings do not cause Github
checks to fail.
@dpgrote
Copy link
Member

dpgrote commented Nov 17, 2025

This check is a good idea. Though, could the check be done in AddPlasma to avoid the extra loop over the injection sources? This would also allow checking the actual particle velocities (including any thermal piece) instead of only checking the bulk momentum.

@n01r
Copy link
Member Author

n01r commented Nov 17, 2025

I wondered which one to choose but thought it would produce a more lightweight check with less overhead, in case there was particle injection every step.
But if we think it's worth it, I could look into checking momentum for every particle in AddPlasma.

@titoiride
Copy link
Member

I'm curious about this check and I have a question: which case do you have in mind? What I'm thinking is that in a "standard" explicit PIC, the CFL is what limits a particle motion, not the particle momentum. If a particle crosses more than one cell in a single step there might be a problem with the temporal step and you should probably lower the CFL. Implicit PIC schemes shouldn't have the cell limitation (AFAIK). In other words, for instance in a 1D case dt = dx and even if a particle were to move with c (infinite momentum), it would cross at most 1 cell per step.

I'm wondering if I'm missing some edge case here.

@n01r
Copy link
Member Author

n01r commented Nov 17, 2025

I'd be happy to flesh out the specific cases with different solvers in follow-ups. Currently, none of our tests throw a warning so we fulfill it everywhere so far. A possible case I have in mind is where you're setting up momentum functions that scale with your coordinate (e.g. in a rigid rotor ring current distribution), but then as you might change resolution to check scaling, you violate the condition.

Of course, not the particle momentum is the more fundamental but as you say, the time step, @titoiride. But in a typical workflow, you calculate your general resolution and then set up your sim but then you might add a beam or plasma current later in addition to the base setup and may forget about the actual speed of the particles. So yea, it's mostly a quality-of-life improvement and hence I kept it more lightweight and out of AddParticles.

@n01r
Copy link
Member Author

n01r commented Nov 17, 2025

I'm wondering if I'm missing some edge case here.

Oh, particularly the hybrid Ohm solver is semi-implicit. You can have (and actually should) have large cells and you can run with large time steps. But only the B-field is evolved with substeps. Ion motion defines the resolution, and the user is required to set const_dt, instead of a CFL (which is ignored).

Maybe we only need the check for the hybrid mode?

@RemiLehe RemiLehe self-requested a review November 25, 2025 23:34
@RemiLehe RemiLehe self-assigned this Nov 25, 2025
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants