Add tolerance to ToFixedT to avoid pz<=0.

ECP-WarpX · Nov 10, 2023 · a3ed094 · a3ed094
1 parent 3733447
commit a3ed094
Showing 1 changed file with 7 additions and 4 deletions.
diff --git a/src/particles/transformation/ToFixedT.H b/src/particles/transformation/ToFixedT.H
@@ -61,10 +61,13 @@ namespace transformation
             amrex::ParticleReal const y = p.pos(RealAoS::y);
             amrex::ParticleReal const t = p.pos(RealAoS::t);
 
+            // small tolerance to avoid NaN for pz<0:
+            amrex::ParticleReal const tol = 1.0e-8_prt;
+
             // compute value of reference pzd = beta*gamma
             amrex::ParticleReal const argd = -1.0_prt + pow(m_ptd, 2);
-            AMREX_ASSERT_WITH_MESSAGE(argd > 0.0_prt, "invalid pzd arg (<=0)");
-            amrex::ParticleReal const pzdf = argd > 0.0_prt ? sqrt(argd) : 0.0_prt;
+            // AMREX_ASSERT_WITH_MESSAGE(argd > 0.0_prt, "invalid pzd arg (<=0)");
+            amrex::ParticleReal const pzdf = argd > 0.0_prt ? sqrt(argd) : tol;
 
             // transform momenta to dynamic units (eg, so that momenta are
             // normalized by mc):
@@ -74,8 +77,8 @@ namespace transformation
 
             // compute value of particle pz = beta*gamma
             amrex::ParticleReal const arg = -1.0_prt + pow(m_ptd+pt, 2) - pow(px, 2) - pow(py, 2);
-            AMREX_ASSERT_WITH_MESSAGE(arg > 0.0_prt, "invalid pz arg (<=0)");
-            amrex::ParticleReal const pzf = arg > 0.0_prt ? sqrt(arg) : 0.0_prt;
+            // AMREX_ASSERT_WITH_MESSAGE(arg > 0.0_prt, "invalid pz arg (<=0)");
+            amrex::ParticleReal const pzf = arg > 0.0_prt ? sqrt(arg) : tol;
 
             // transform position and momentum (from fixed s to fixed t)
             p.pos(RealAoS::x) = x + px*t/(m_ptd+pt);