Enforce particle orbit bounds for implicit solvers

[JustinRayAngus]

This PR adds particle bounds checks prior to field gather and current deposition for the implicit solvers. The checks are included for each suborbit when particle suborbits are used.

This PR also adds checks that the number of cell crossings in each direction are within the bounds determined by particles.max_grid_crossings when using the mass matrices for the Jacobian and using the Villasenor deposition. For that configuration, particles.max_grid_crossings determines that size of the mass matrix containers.

Note that the particle bounds check prior to gather/deposition is not sufficient to ensure that each particle orbit remains within the maximum allowed value set by particles.max_grid_crossings. This is because a particle with a large orbit can be located anywhere within the tilebox, whereas the bounds check only catches particles near the tilebox boundaries.

When a particle is found to be out of bounds, or has too many cell crossings, the simulation aborts. In the future, the abort can be replaced with a trigger to subcycle the implicit step (see PR #7000).

This PR replaces PR #6156.

[JustinRayAngus]

@dpgrote @atmyers

I may be able to simplify some of the checks I've added here for particles being within range using:

AMREX_ASSERT_WITH_MESSAGE(
    amrex::numParticlesOutOfRange(pti, range) == 0,
    "Particles shape does not fit within tile (CPU) or guard cells (GPU) used for current deposition");

But it's not clear to me what this actually checks. I presume it uses the stored position of the particle and checks that it is within the specified range of the valid region of the box/tile that it lives in. Is that correct?

If yes, then I would need to convert the position from time-centered to time n+1 prior to calling this.

How does it work for curvilinear geometries? Does it convert xp and yp to rp = sqrt(xp^2 + yp^2)?

[dpgrote]

I don't know if it matters, but this is duplicating the calculation of xp_np1 here and also of x_new below (excepting lo) when called from the deposition routines. Is there a nice way of handling this to avoid this duplication? This might make this routine simpler, avoiding the code for different geometries. Though maybe not since the pusher routine which calls this doesn't calculate x_new.

Though a related question is why is the pusher routine calling this and not the gather where it matters?

[JustinRayAngus]

Yes, there is duplicate code for computing xp_np1. It's not just here though, it's in a lot of places. I've thought about this previously and I think it would be optimal to change the paradigm that the particle position lives at time n + 1/2 during the implicit solve to living at time n + 1. This should be done in a separate PR, as it will be a much bigger change.

Here is what CodeX has to say about putting the check inside doGatherShapeNImplicit.

I would keep the check in PushXPSingleStep, not move it into doGatherShapeNImplicit().

  Reasons:

  - doGatherShapeNImplicit() does not currently know nodal_lo, nodal_hi, or the error counter. Adding them
    would couple a low-level gather routine to pusher error handling.

  - Gathering can be disabled with do_gather == false; the particle-position invariant should still be
    enforced.

  - The final Picard position update may not be followed by another gather, so a gather-only check could miss
    the final invalid position.

  - Deposition also needs the invariant independently, since it is another memory-access boundary.
  - The current placement returns before an unsafe gather and validates every updated position, including the
    final one.

  A check inside the gather could provide extra defensive protection, but it could not replace the pusher and
  deposition checks. It would mainly duplicate them while requiring additional arguments or a changed return
  type.

  The clean ownership is:

  - PushXPSingleStep: validate every trial orbit and stop before gathering.
  - Implicit deposition routines: validate before depositing.
  - Host wrappers: inspect device counters and abort.

[JustinRayAngus]

@dpgrote The current calls to ParticleUtils::isImplicitParticlePositionInBounds() in ImplicitPushXP.cpp are placed after the particle position updates. This ensures that each particle is in bounds for both the gather and deposition operations, the latter of which occurs locally for suborbit particles.

[JustinRayAngus]

@dpgrote I could add a boolean flag to ParticleUtils::isImplicitParticlePositionInBounds() that says whether the passed position is to be interpreted at time n + 1/2 or at time n + 1, but I'm not sure that is a good design, or that it really saves that much.

[dpgrote]

The code is all good.