avoid NaN in maxParticleVelocity() (#6831)

This PR is an attempt to fix a random job hang that I've experienced
recently on Tuolumne.

I've tracked the hang down to `mypc->maxParticleVelocity()`, which is
called from `UpdateDtFromParticleSpeeds()` when using dynamic time
stepping. The code there does not guard against boxes with zero
particles, which eventually results in taking the sqrt of an
uninitialized value resulting in local `max_v = nan`. This PR avoids NaN
values.

With this fix, I have not yet observed any random job hanging.

Author: JustinRayAngus

Source/Particles/WarpXParticleContainer.cpp

@@ -2331,25 +2331,32 @@ amrex::ParticleReal WarpXParticleContainer::maxParticleVelocity(bool local) {
     ReduceOps<ReduceOpMax> reduce_op;
     ReduceData<ParticleReal> reduce_data(reduce_op);
 
+    amrex::Long np_total = 0;
+
     const int nLevels = finestLevel();
 
 #ifdef AMREX_USE_OMP
-#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#pragma omp parallel reduction(+:np_total) if (amrex::Gpu::notInLaunchRegion())
 #endif
     for (int lev = 0; lev <= nLevels; ++lev) {
         for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
         {
+            const amrex::Long np = pti.numParticles();
+            if (np == 0) { continue; }
+
+            np_total += np;
+
             auto *const ux = pti.GetAttribs(PIdx::ux).data();
             auto *const uy = pti.GetAttribs(PIdx::uy).data();
             auto *const uz = pti.GetAttribs(PIdx::uz).data();
 
-            reduce_op.eval(pti.numParticles(), reduce_data,
+            reduce_op.eval(np, reduce_data,
                 [=] AMREX_GPU_DEVICE (int ip)
                 { return (ux[ip]*ux[ip] + uy[ip]*uy[ip] + uz[ip]*uz[ip]) * inv_c2; });
         }
     }
 
-    const amrex::ParticleReal max_usq = amrex::get<0>(reduce_data.value());
+    const amrex::ParticleReal max_usq = (np_total > 0 ? amrex::get<0>(reduce_data.value()) : 0._prt);
 
     const amrex::ParticleReal gaminv = 1.0_prt/std::sqrt(1.0_prt + max_usq);
     amrex::ParticleReal max_v = gaminv * std::sqrt(max_usq) * PhysConst::c;