safeguard against 0 denominator

RevathiJambunathan · Oct 25, 2023 · c74bc3c · c74bc3c
1 parent b657b99
commit c74bc3c
Showing 1 changed file with 39 additions and 37 deletions.
diff --git a/Source/Particles/Pusher/UpdateMomentumBoris.H b/Source/Particles/Pusher/UpdateMomentumBoris.H
@@ -81,43 +81,45 @@ void UpdateMomentumBoris(
                   / ( (B_sq - E_sq) + std::sqrt( (B_sq - E_sq) * (B_sq - E_sq ) + 4._prt * EdotB_sq ) );
 
         // only if E << B -- else fix with thesis
-        amrex::Real vx_drift = (Ey*Bz - Ez*By)/(B_sq + Eprime_sq);
-        amrex::Real vy_drift = (Ez*Bx - Ex*Bz)/(B_sq + Eprime_sq);
-        amrex::Real vz_drift = (Ex*By - Ey*Bx)/(B_sq + Eprime_sq);
-
-        // 1(b) unit vector along drift velocity
-        amrex::Real drift_mag = std::sqrt(vx_drift*vx_drift + vy_drift*vy_drift + vz_drift*vz_drift);
-        amrex::Real vxd_hat = vx_drift/drift_mag;
-        amrex::Real vyd_hat = vy_drift/drift_mag;
-        amrex::Real vzd_hat = vz_drift/drift_mag;
-
-        // 1(c) gamma_drift = gamma_o = 1 / (sqrt ( 1 - drift_mag^2/c^2) )
-        amrex::Real gamma_drift = 1._prt/std::sqrt( 1._prt - drift_mag*drift_mag*inv_c2);
-
-        // 2 Bfield transform
-        // By construction, only Bperp to drifting frame is non-zero
-        amrex::Real Bx_prime = gamma_drift * ( Bx - inv_c2 * (vy_drift*Ez - vz_drift*Ey) );
-        amrex::Real By_prime = gamma_drift * ( By - inv_c2 * (vz_drift*Ex - vx_drift*Ez) );
-        amrex::Real Bz_prime = gamma_drift * ( Bz - inv_c2 * (vx_drift*Ey - vy_drift*Ex) );
-
-        amrex::Real Bprime_mag = std::sqrt(Bx_prime*Bx_prime + By_prime*By_prime + Bz_prime*Bz_prime);
-        amrex::Real Bxprime_hat = Bx_prime/Bprime_mag;
-        amrex::Real Byprime_hat = By_prime/Bprime_mag;
-        amrex::Real Bzprime_hat = Bz_prime/Bprime_mag;
-
-        // New u momentum in drift frame
-        amrex::Real udotBprimehat = ux*Bxprime_hat + uy*Byprime_hat + uz*Bzprime_hat;
-        amrex::Real ux_new_prime = udotBprimehat * Bxprime_hat;
-        amrex::Real uy_new_prime = udotBprimehat * Byprime_hat;
-        amrex::Real uz_new_prime = udotBprimehat * Bzprime_hat;
-
-        // New u momentum in lab frame
-        amrex::Real gamma_new_prime = std::sqrt( 1._prt + inv_c2 * ( ux_new_prime*ux_new_prime
-                                                                   + uy_new_prime*uy_new_prime
-                                                                   + uz_new_prime*uz_new_prime ) );
-        ux = gamma_drift * ( ux_new_prime + gamma_new_prime * vx_drift);
-        uy = gamma_drift * ( uy_new_prime + gamma_new_prime * vy_drift);
-        uz = gamma_drift * ( uz_new_prime + gamma_new_prime * vz_drift);
+        if ( (B_sq + Eprime_sq) > 0. ) {
+            amrex::Real vx_drift = (Ey*Bz - Ez*By)/(B_sq + Eprime_sq);
+            amrex::Real vy_drift = (Ez*Bx - Ex*Bz)/(B_sq + Eprime_sq);
+            amrex::Real vz_drift = (Ex*By - Ey*Bx)/(B_sq + Eprime_sq);
+
+            // 1(b) unit vector along drift velocity
+            amrex::Real drift_mag = std::sqrt(vx_drift*vx_drift + vy_drift*vy_drift + vz_drift*vz_drift);
+            amrex::Real vxd_hat = vx_drift/drift_mag;
+            amrex::Real vyd_hat = vy_drift/drift_mag;
+            amrex::Real vzd_hat = vz_drift/drift_mag;
+
+            // 1(c) gamma_drift = gamma_o = 1 / (sqrt ( 1 - drift_mag^2/c^2) )
+            amrex::Real gamma_drift = 1._prt/std::sqrt( 1._prt - drift_mag*drift_mag*inv_c2);
+
+            // 2 Bfield transform
+            // By construction, only Bperp to drifting frame is non-zero
+            amrex::Real Bx_prime = gamma_drift * ( Bx - inv_c2 * (vy_drift*Ez - vz_drift*Ey) );
+            amrex::Real By_prime = gamma_drift * ( By - inv_c2 * (vz_drift*Ex - vx_drift*Ez) );
+            amrex::Real Bz_prime = gamma_drift * ( Bz - inv_c2 * (vx_drift*Ey - vy_drift*Ex) );
+
+            amrex::Real Bprime_mag = std::sqrt(Bx_prime*Bx_prime + By_prime*By_prime + Bz_prime*Bz_prime);
+            amrex::Real Bxprime_hat = Bx_prime/Bprime_mag;
+            amrex::Real Byprime_hat = By_prime/Bprime_mag;
+            amrex::Real Bzprime_hat = Bz_prime/Bprime_mag;
+
+            // New u momentum in drift frame
+            amrex::Real udotBprimehat = ux*Bxprime_hat + uy*Byprime_hat + uz*Bzprime_hat;
+            amrex::Real ux_new_prime = udotBprimehat * Bxprime_hat;
+            amrex::Real uy_new_prime = udotBprimehat * Byprime_hat;
+            amrex::Real uz_new_prime = udotBprimehat * Bzprime_hat;
+
+            // New u momentum in lab frame
+            amrex::Real gamma_new_prime = std::sqrt( 1._prt + inv_c2 * ( ux_new_prime*ux_new_prime
+                                                                       + uy_new_prime*uy_new_prime
+                                                                       + uz_new_prime*uz_new_prime ) );
+            ux = gamma_drift * ( ux_new_prime + gamma_new_prime * vx_drift);
+            uy = gamma_drift * ( uy_new_prime + gamma_new_prime * vy_drift);
+            uz = gamma_drift * ( uz_new_prime + gamma_new_prime * vz_drift);
+        }
     }
 #endif