Skip to content

Enforce particle orbit bounds for implicit solvers #7225

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
RemiLehe merged 23 commits into from
Sep 1, 2026
Merged

Enforce particle orbit bounds for implicit solvers #7225

Show file tree
Hide file tree
Changes from 7 commits
Commits
Show all changes
23 commits
Select commit Hold shift + click to select a range
f03f38f
refactor grid_crossing checks when using MM.
JustinRayAngus Aug 30, 2026
bc331e7
add cell crossings check to villasenor deposition.
JustinRayAngus Aug 30, 2026
776ddcf
check that suborbit position is within range.
JustinRayAngus Aug 30, 2026
a525170
catch all positions in suborbits.
JustinRayAngus Aug 30, 2026
abf3673
clean up message.
JustinRayAngus Aug 31, 2026
095b026
add param to brief.
JustinRayAngus Aug 31, 2026
6c9ef26
mark maybe unused.
JustinRayAngus Aug 31, 2026
bdafff1
refactored.
JustinRayAngus Sep 1, 2026
dce19bc
spacing.
JustinRayAngus Sep 1, 2026
ebe9e31
update briefs.
JustinRayAngus Sep 1, 2026
ecd076e
maybe unusd.
JustinRayAngus Sep 1, 2026
d9e6870
add guards.
JustinRayAngus Sep 1, 2026
77c9a2e
capture lambda issue fixed.
JustinRayAngus Sep 1, 2026
88918cb
use const.
JustinRayAngus Sep 1, 2026
ca1c507
nodal box bug fix.
JustinRayAngus Sep 1, 2026
5b1b4bc
bug fix.
JustinRayAngus Sep 1, 2026
4d35717
remove unneeded check in DepositMassMatrices.
JustinRayAngus Sep 1, 2026
ea1e8bf
add bounds check to implicit Esirkepov deposition.
JustinRayAngus Sep 1, 2026
ed73963
use device-safe check.
JustinRayAngus Sep 1, 2026
1f7b312
device-safe position handling for pusher.
JustinRayAngus Sep 1, 2026
eb28e12
device-safe position handling for deposit.
JustinRayAngus Sep 1, 2026
7b00cf7
PushXPSingleStep now discerns between unconverged and out-of-bounds.
JustinRayAngus Sep 1, 2026
9744790
higher-level error handling for ALL implicit depositions.
JustinRayAngus Sep 1, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
84 changes: 43 additions & 41 deletions Source/FieldSolver/ImplicitSolvers/ImplicitSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -752,13 +752,15 @@ void ImplicitSolver::InitializeMassMatrices ()
}
else if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Villasenor) {
#ifndef WARPX_DIM_3D
const int max_crossings = ngJ[0] - shape + 1;
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(max_crossings > 0,
const int max_grid_crossings = ngJ[0] - shape + 1;
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(max_grid_crossings > 0,
"Mass Matrices for Jacobian with Villasenor deposition requires particles.max_grid_crossings > 0.");
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(max_crossings == m_WarpX->particle_max_grid_crossings,
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(max_grid_crossings == WarpX::particle_max_grid_crossings,
"Guard cells for J are not consistent with particle_max_grid_crossings.");
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(max_crossings <= 2,
"Mass Matrices for Jacobian with Villasenor deposition requires particles.max_grid_crossings <= 2.");
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
max_grid_crossings <= WarpX::villasenor_mass_matrices_max_grid_crossings,
"Mass matrices for the Jacobian with Villasenor deposition support "
"particles.max_grid_crossings <= WarpX::villasenor_mass_matrices_max_grid_crossings.");
#endif
// Comment on direction-dependent number of mass matrices components
// set below for charge-conserving Villasenor deposition:
Expand All @@ -767,47 +769,47 @@ void ImplicitSolver::InitializeMassMatrices ()
// 1 + 2*shape (both comps nodal)
#if defined(WARPX_DIM_1D_Z)
// x and y are nodal, z is centered
m_ncomp_xx[0] = 1 + 2*shape + 2*max_crossings;
m_ncomp_xy[0] = 1 + 2*shape + 2*max_crossings;
m_ncomp_xz[0] = 0 + 2*shape + 2*max_crossings;
m_ncomp_yx[0] = 1 + 2*shape + 2*max_crossings;
m_ncomp_yy[0] = 1 + 2*shape + 2*max_crossings;
m_ncomp_yz[0] = 0 + 2*shape + 2*max_crossings;
m_ncomp_zx[0] = 0 + 2*shape + 2*max_crossings;
m_ncomp_zy[0] = 0 + 2*shape + 2*max_crossings;
m_ncomp_zz[0] = 1 + 2*(shape-1) + 2*max_crossings;
m_ncomp_xx[0] = 1 + 2*shape + 2*max_grid_crossings;
m_ncomp_xy[0] = 1 + 2*shape + 2*max_grid_crossings;
m_ncomp_xz[0] = 0 + 2*shape + 2*max_grid_crossings;
m_ncomp_yx[0] = 1 + 2*shape + 2*max_grid_crossings;
m_ncomp_yy[0] = 1 + 2*shape + 2*max_grid_crossings;
m_ncomp_yz[0] = 0 + 2*shape + 2*max_grid_crossings;
m_ncomp_zx[0] = 0 + 2*shape + 2*max_grid_crossings;
m_ncomp_zy[0] = 0 + 2*shape + 2*max_grid_crossings;
m_ncomp_zz[0] = 1 + 2*(shape-1) + 2*max_grid_crossings;
#elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
// x is centered, y and z are nodal
m_ncomp_xx[0] = 1 + 2*(shape-1) + 2*max_crossings;
m_ncomp_xy[0] = 0 + 2*shape + 2*max_crossings;
m_ncomp_xz[0] = 0 + 2*shape + 2*max_crossings;
m_ncomp_yx[0] = 0 + 2*shape + 2*max_crossings;
m_ncomp_yy[0] = 1 + 2*shape + 2*max_crossings;
m_ncomp_yz[0] = 1 + 2*shape + 2*max_crossings;
m_ncomp_zx[0] = 0 + 2*shape + 2*max_crossings;
m_ncomp_zy[0] = 1 + 2*shape + 2*max_crossings;
m_ncomp_zz[0] = 1 + 2*shape + 2*max_crossings;
m_ncomp_xx[0] = 1 + 2*(shape-1) + 2*max_grid_crossings;
m_ncomp_xy[0] = 0 + 2*shape + 2*max_grid_crossings;
m_ncomp_xz[0] = 0 + 2*shape + 2*max_grid_crossings;
m_ncomp_yx[0] = 0 + 2*shape + 2*max_grid_crossings;
m_ncomp_yy[0] = 1 + 2*shape + 2*max_grid_crossings;
m_ncomp_yz[0] = 1 + 2*shape + 2*max_grid_crossings;
m_ncomp_zx[0] = 0 + 2*shape + 2*max_grid_crossings;
m_ncomp_zy[0] = 1 + 2*shape + 2*max_grid_crossings;
m_ncomp_zz[0] = 1 + 2*shape + 2*max_grid_crossings;
#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
// dir = 0: x is centered, y and z are nodal
m_ncomp_xx[0] = 1 + 2*(shape-1) + 2*max_crossings;
m_ncomp_xy[0] = 0 + 2*shape + 2*max_crossings;
m_ncomp_xz[0] = 0 + 2*shape + 2*max_crossings;
m_ncomp_yx[0] = 0 + 2*shape + 2*max_crossings;
m_ncomp_yy[0] = 1 + 2*shape + 2*max_crossings;
m_ncomp_yz[0] = 1 + 2*shape + 2*max_crossings;
m_ncomp_zx[0] = 0 + 2*shape + 2*max_crossings;
m_ncomp_zy[0] = 1 + 2*shape + 2*max_crossings;
m_ncomp_zz[0] = 1 + 2*shape + 2*max_crossings;
m_ncomp_xx[0] = 1 + 2*(shape-1) + 2*max_grid_crossings;
m_ncomp_xy[0] = 0 + 2*shape + 2*max_grid_crossings;
m_ncomp_xz[0] = 0 + 2*shape + 2*max_grid_crossings;
m_ncomp_yx[0] = 0 + 2*shape + 2*max_grid_crossings;
m_ncomp_yy[0] = 1 + 2*shape + 2*max_grid_crossings;
m_ncomp_yz[0] = 1 + 2*shape + 2*max_grid_crossings;
m_ncomp_zx[0] = 0 + 2*shape + 2*max_grid_crossings;
m_ncomp_zy[0] = 1 + 2*shape + 2*max_grid_crossings;
m_ncomp_zz[0] = 1 + 2*shape + 2*max_grid_crossings;
// dir = 1: x and y are nodal, z is centered
m_ncomp_xx[1] = 1 + 2*shape + 2*max_crossings;
m_ncomp_xy[1] = 1 + 2*shape + 2*max_crossings;
m_ncomp_xz[1] = 0 + 2*shape + 2*max_crossings;
m_ncomp_yx[1] = 1 + 2*shape + 2*max_crossings;
m_ncomp_yy[1] = 1 + 2*shape + 2*max_crossings;
m_ncomp_yz[1] = 0 + 2*shape + 2*max_crossings;
m_ncomp_zx[1] = 0 + 2*shape + 2*max_crossings;
m_ncomp_zy[1] = 0 + 2*shape + 2*max_crossings;
m_ncomp_zz[1] = 1 + 2*(shape-1) + 2*max_crossings;
m_ncomp_xx[1] = 1 + 2*shape + 2*max_grid_crossings;
m_ncomp_xy[1] = 1 + 2*shape + 2*max_grid_crossings;
m_ncomp_xz[1] = 0 + 2*shape + 2*max_grid_crossings;
m_ncomp_yx[1] = 1 + 2*shape + 2*max_grid_crossings;
m_ncomp_yy[1] = 1 + 2*shape + 2*max_grid_crossings;
m_ncomp_yz[1] = 0 + 2*shape + 2*max_grid_crossings;
m_ncomp_zx[1] = 0 + 2*shape + 2*max_grid_crossings;
m_ncomp_zy[1] = 0 + 2*shape + 2*max_grid_crossings;
m_ncomp_zz[1] = 1 + 2*(shape-1) + 2*max_grid_crossings;
#endif
for (int dir=0; dir<AMREX_SPACEDIM; dir++) {
Nc_tot_xx *= m_ncomp_xx[dir];
Expand Down
45 changes: 28 additions & 17 deletions Source/Particles/Deposition/CurrentDeposition.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1497,6 +1497,7 @@ void doChargeConservingDepositionShapeNImplicit ([[maybe_unused]]const amrex::Pa
* \param lo Index lower bounds of domain
* \param invvol One over cell volume
* \param n_rz_azimuthal_modes Number of azimuthal modes when using RZ geometry
* \param max_grid_crossings Maximum allowed number of grid crossings in each direction.
*/
template <int depos_order>
AMREX_GPU_HOST_DEVICE AMREX_INLINE
Expand All @@ -1521,7 +1522,8 @@ void VillasenorDepositionShapeNKernel ([[maybe_unused]] amrex::ParticleReal cons
const amrex::GpuArray<amrex::GpuArray<bool,2>, AMREX_SPACEDIM> & do_cropping,
amrex::Dim3 const lo,
amrex::Real const invvol,
[[maybe_unused]] int const n_rz_azimuthal_modes)
[[maybe_unused]] int const n_rz_azimuthal_modes,
int const max_grid_crossings)
{

using namespace amrex::literals;
Expand Down Expand Up @@ -1664,9 +1666,14 @@ void VillasenorDepositionShapeNKernel ([[maybe_unused]] amrex::ParticleReal cons
const int cell_crossings_z = std::abs(k_new-k_old);
num_segments += cell_crossings_z;

// need to assert that the number of cell crossings in each direction
// is within the range permitted by the number of guard cells
// e.g., if (num_segments > 7) ...
AMREX_ALWAYS_ASSERT_WITH_MESSAGE( cell_crossings_x <= max_grid_crossings,
"Villasenor deposition error: x-cell_crossings exceeds max_grid_crossings.");

AMREX_ALWAYS_ASSERT_WITH_MESSAGE( cell_crossings_y <= max_grid_crossings,
"Villasenor deposition error: y-cell_crossings exceeds max_grid_crossings.");

AMREX_ALWAYS_ASSERT_WITH_MESSAGE( cell_crossings_z <= max_grid_crossings,
"Villasenor deposition error: z-cell_crossings exceeds max_grid_crossings.");

// compute the initial cell location used to find the cell crossings.
const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
Expand Down Expand Up @@ -1845,9 +1852,11 @@ void VillasenorDepositionShapeNKernel ([[maybe_unused]] amrex::ParticleReal cons
const int cell_crossings_z = std::abs(k_new-k_old);
num_segments += cell_crossings_z;

// need to assert that the number of cell crossings in each direction
// is within the range permitted by the number of guard cells
// e.g., if (num_segments > 5) ...
AMREX_ALWAYS_ASSERT_WITH_MESSAGE( cell_crossings_x <= max_grid_crossings,
"Villasenor deposition error: x-cell_crossings exceeds max_grid_crossings.");

AMREX_ALWAYS_ASSERT_WITH_MESSAGE( cell_crossings_z <= max_grid_crossings,
"Villasenor deposition error: z-cell_crossings exceeds max_grid_crossings.");

// compute the initial cell location used to find the cell crossings.
const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
Expand Down Expand Up @@ -2009,9 +2018,8 @@ void VillasenorDepositionShapeNKernel ([[maybe_unused]] amrex::ParticleReal cons
const int cell_crossings_z = std::abs(k_new-k_old);
num_segments += cell_crossings_z;

// need to assert that the number of cell crossings in each direction
// is within the range permitted by the number of guard cells
// e.g., if (num_segments > 3) ...
AMREX_ALWAYS_ASSERT_WITH_MESSAGE( cell_crossings_z <= max_grid_crossings,
"Villasenor deposition error: z-cell_crossings exceeds max_grid_crossings.");

// compute the initial cell location used to find the cell crossings.
const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
Expand Down Expand Up @@ -2088,9 +2096,8 @@ void VillasenorDepositionShapeNKernel ([[maybe_unused]] amrex::ParticleReal cons
const int cell_crossings_x = std::abs(i_new-i_old);
num_segments += cell_crossings_x;

// need to assert that the number of cell crossings in each direction
// is within the range permitted by the number of guard cells
// e.g., if (num_segments > 3) ...
AMREX_ALWAYS_ASSERT_WITH_MESSAGE( cell_crossings_x <= max_grid_crossings,
"Villasenor deposition error: x-cell_crossings exceeds max_grid_crossings.");

// compute the initial cell location used to find the cell crossings.
const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
Expand Down Expand Up @@ -2188,6 +2195,7 @@ void VillasenorDepositionShapeNKernel ([[maybe_unused]] amrex::ParticleReal cons
* \param lo Index lower bounds of domain.
* \param q species charge.
* \param n_rz_azimuthal_modes Number of azimuthal modes when using RZ geometry.
* \param max_grid_crossings Maximum allowed number of grid crossings in each direction.
*/
template <int depos_order>
void doVillasenorDepositionShapeNExplicit (const GetParticlePosition<PIdx>& GetPosition,
Expand All @@ -2208,7 +2216,8 @@ void doVillasenorDepositionShapeNExplicit (const GetParticlePosition<PIdx>& GetP
const amrex::GpuArray<amrex::GpuArray<bool,2>, AMREX_SPACEDIM> & do_cropping,
const amrex::Dim3 lo,
const amrex::Real q,
[[maybe_unused]] const int n_rz_azimuthal_modes)
[[maybe_unused]] const int n_rz_azimuthal_modes,
const int max_grid_crossings)
{
using namespace amrex::literals;

Expand Down Expand Up @@ -2249,7 +2258,7 @@ void doVillasenorDepositionShapeNExplicit (const GetParticlePosition<PIdx>& GetP
uxp[ip], uyp[ip], uzp[ip], gaminv,
Jx_arr, Jy_arr, Jz_arr,
dt, dinv, xyzmin, domain_double, do_cropping,
lo, invvol, n_rz_azimuthal_modes);
lo, invvol, n_rz_azimuthal_modes, max_grid_crossings);

});
}
Expand Down Expand Up @@ -2281,6 +2290,7 @@ void doVillasenorDepositionShapeNExplicit (const GetParticlePosition<PIdx>& GetP
* \param lo Index lower bounds of domain.
* \param q species charge.
* \param n_rz_azimuthal_modes Number of azimuthal modes when using RZ geometry.
* \param max_grid_crossings Maximum allowed number of grid crossings in each direction.
*/
template <int depos_order>
void doVillasenorDepositionShapeNImplicit ([[maybe_unused]]const amrex::ParticleReal * const xp_n_data,
Expand All @@ -2306,7 +2316,8 @@ void doVillasenorDepositionShapeNImplicit ([[maybe_unused]]const amrex::Particle
const amrex::GpuArray<amrex::GpuArray<bool,2>, AMREX_SPACEDIM> & do_cropping,
const amrex::Dim3 lo,
const amrex::Real q,
[[maybe_unused]] const int n_rz_azimuthal_modes)
[[maybe_unused]] const int n_rz_azimuthal_modes,
const int max_grid_crossings)
{
using namespace amrex::literals;

Expand Down Expand Up @@ -2357,7 +2368,7 @@ void doVillasenorDepositionShapeNImplicit ([[maybe_unused]]const amrex::Particle
uxp_nph[ip], uyp_nph[ip], uzp_nph[ip], gaminv,
Jx_arr, Jy_arr, Jz_arr,
dt, dinv, xyzmin, domain_double, do_cropping,
lo, invvol, n_rz_azimuthal_modes);
lo, invvol, n_rz_azimuthal_modes, max_grid_crossings);

});
}
Expand Down
Loading
Loading