From 13ed05a26fc08b8508bf87d316d11a63d3fbb64b Mon Sep 17 00:00:00 2001
From: RevathiJambunathan <revanathan@gmail.com>
Date: Tue, 14 Nov 2023 11:51:18 -0800
Subject: [PATCH]  print to debug

---
 Source/Particles/Pusher/PushSelector.H        | 27 +++++++++++++++++++
 Source/Particles/Pusher/UpdateMomentumBoris.H | 13 ++++++++-
 2 files changed, 39 insertions(+), 1 deletion(-)

diff --git a/Source/Particles/Pusher/PushSelector.H b/Source/Particles/Pusher/PushSelector.H
index 4b85ab68351..cfd02c77029 100644
--- a/Source/Particles/Pusher/PushSelector.H
+++ b/Source/Particles/Pusher/PushSelector.H
@@ -73,6 +73,11 @@ void doParticlePush(const GetParticlePosition& GetPosition,
 #endif
                     const amrex::Real dt, amrex::Real * PulsarDiag = nullptr)
 {
+#ifdef PULSAR
+    amrex::Real ux_pre = ux;
+    amrex::Real uy_pre = uy;
+    amrex::Real uz_pre = uz;
+#endif
     amrex::ParticleReal qp = a_q;
     if (ion_lev) { qp *= ion_lev; }
 
@@ -106,8 +111,19 @@ void doParticlePush(const GetParticlePosition& GetPosition,
             }
             amrex::ParticleReal x, y, z;
             GetPosition(i, x, y, z);
+#ifdef PULSAR
+            amrex::Real x_pre = x;
+            amrex::Real y_pre = y;
+            amrex::Real z_pre = z;
+#endif
             UpdatePosition(x, y, z, ux, uy, uz, dt );
             SetPosition(i, x, y, z);
+#ifdef PULSAR
+            if ( std::sqrt(x*x + y*y + z*z) > 125200 )
+            {
+                printf(" ux_pre = %e  uy_pre = %e  uz_pre = %e ux_post = %e uy_post = %e uz_post= %e x_pre = %f y_pre=%f z_pre=%f x=%f y=%f z=%f \n", ux_pre, uy_pre, uz_pre, ux, uy, uz, x_pre, y_pre, z_pre, x, y, z );
+            }
+#endif
         } else
 #endif
         {
@@ -136,8 +152,19 @@ void doParticlePush(const GetParticlePosition& GetPosition,
                              dt);
         amrex::ParticleReal x, y, z;
         GetPosition(i, x, y, z);
+#ifdef PULSAR
+        amrex::Real x_pre = x;
+        amrex::Real y_pre = y;
+        amrex::Real z_pre = z;
+#endif
         UpdatePosition(x, y, z, ux, uy, uz, dt );
         SetPosition(i, x, y, z);
+#ifdef PULSAR
+        if ( std::sqrt(x*x + y*y + z*z) > 125200 )
+        {
+            printf(" ux_pre = %e  uy_pre = %e  uz_pre = %e ux_post = %e uy_post = %e uz_post= %e x_pre = %f y_pre=%f z_pre=%f x=%f y=%f z=%f \n", ux_pre, uy_pre, uz_pre, ux, uy, uz, x_pre, y_pre, z_pre, x, y, z );
+        }
+#endif
     } else if (pusher_algo == ParticlePusherAlgo::Vay) {
         UpdateMomentumVay( ux, uy, uz,
                            Ex, Ey, Ez, Bx,
diff --git a/Source/Particles/Pusher/UpdateMomentumBoris.H b/Source/Particles/Pusher/UpdateMomentumBoris.H
index a4ab3ff4a3f..33e2bfbb3cd 100644
--- a/Source/Particles/Pusher/UpdateMomentumBoris.H
+++ b/Source/Particles/Pusher/UpdateMomentumBoris.H
@@ -29,6 +29,10 @@ void UpdateMomentumBoris(
 
     const amrex::ParticleReal econst = 0.5_prt*q*dt/m;
 
+#ifdef PULSAR
+    amrex::Real abs_u_pre = std::sqrt(ux*ux + uy*uy + uz*uz);
+#endif
+
     // First half-push for E
     ux += econst*Ex;
     uy += econst*Ey;
@@ -92,7 +96,7 @@ void UpdateMomentumBoris(
 
             // 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);
-            if (drift_mag < (0.99 * PhysConst::c) ) {
+            if (drift_mag < (PhysConst::c) ) {
 	        amrex::Real vxd_hat = vx_drift/drift_mag;
                 amrex::Real vyd_hat = vy_drift/drift_mag;
                 amrex::Real vzd_hat = vz_drift/drift_mag;
@@ -124,6 +128,13 @@ void UpdateMomentumBoris(
                 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);
+
+                // debugging
+                amrex::Real abs_u_post = std::sqrt(ux*ux + uy*uy + uz*uz);
+                if ( abs_u_post > abs_u_pre || (drift_mag > 0.99999999*PhysConst::c))
+                {
+                    printf(" premom = %e  postmom = %e  gammadrift = %e gammanewprime = %e B_sq = %e Eprime_sq = %e drift_mag = %e \n", abs_u_pre, abs_u_post, gamma_drift, gamma_new_prime, B_sq, Eprime_sq, drift_mag );
+                }
             }
         }
     }