add do zerouperpB as flag

RevathiJambunathan · Oct 25, 2023 · b657b99 · b657b99
1 parent b471e06
commit b657b99
Show file tree

Hide file tree

Showing 6 changed files with 93 additions and 62 deletions.
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
@@ -2855,6 +2855,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
     amrex::Real gather_buffer_boxmin = Pulsar::m_gatherbuffer_min;
     amrex::Real gather_buffer_boxmax = Pulsar::m_gatherbuffer_max;
     amrex::Real re_scaledratio = Pulsar::m_re_scaledratio;
+    int do_zero_uperpB_driftframe = Pulsar::m_do_zero_uperpB_driftframe;
 #endif
 
 #ifdef AMREX_USE_OMP
@@ -3019,15 +3020,19 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
                                                              Exp, Eyp, Ezp, Bxp,
                                                              Byp, Bzp, qp, m,
 #ifdef PULSAR
-                                                             re_scaledratio,
+                                                             re_scaledratio, do_zero_uperpB_driftframe,
 #endif
                                                              dt);
                 } else if (pusher_algo == ParticlePusherAlgo::Boris) {
                     amrex::Real qp = q;
                     if (ion_lev) { qp *= ion_lev[ip]; }
                     UpdateMomentumBoris( ux[ip], uy[ip], uz[ip],
                                          Exp, Eyp, Ezp, Bxp,
-                                         Byp, Bzp, qp, m, dt);
+                                         Byp, Bzp, qp, m,
+#ifdef PULSAR
+                                         do_zero_uperpB_driftframe,
+#endif
+                                         dt);
                 } else if (pusher_algo == ParticlePusherAlgo::Vay) {
                     amrex::Real qp = q;
                     if (ion_lev){ qp *= ion_lev[ip]; }
@@ -3164,6 +3169,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti,
     amrex::Real gather_buffer_boxmin = Pulsar::m_gatherbuffer_min;
     amrex::Real gather_buffer_boxmax = Pulsar::m_gatherbuffer_max;
     amrex::Real re_scaledratio = Pulsar::m_re_scaledratio;
+    int do_zero_uperpB_driftframe = Pulsar::m_do_zero_uperpB_driftframe;
 #endif
 
 
@@ -3415,7 +3421,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti,
                               t_chi_max,
 #endif
 #ifdef PULSAR
-                              re_scaledratio,
+                              re_scaledratio, do_zero_uperpB_driftframe,
 #endif
                               dt);
 	   // PulsarPartDiagData);
@@ -3455,7 +3461,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti,
                                   m, q, pusher_algo, do_crr, do_copy,
                                   t_chi_max,
 #ifdef PULSAR
-                                  re_scaledratio,
+                                  re_scaledratio, do_zero_uperpB_driftframe,
 #endif
                                   dt);
             }

diff --git a/Source/Particles/PulsarParameters.cpp b/Source/Particles/PulsarParameters.cpp
@@ -153,7 +153,7 @@ amrex::Real Pulsar::m_gammarad_RR;
 amrex::Real Pulsar::m_BLC;
 amrex::Real Pulsar::m_beta_rec_RR;
 amrex::Real Pulsar::m_re_scaledratio = 1.;
-int Pulsar::m_do_zero_uperpB_driftframe;
+int Pulsar::m_do_zero_uperpB_driftframe = 0;
 
 Pulsar::Pulsar ()
 {

diff --git a/Source/Particles/Pusher/PushSelector.H b/Source/Particles/Pusher/PushSelector.H
@@ -67,6 +67,7 @@ void doParticlePush(const GetParticlePosition& GetPosition,
 #endif
 #ifdef PULSAR
                     const amrex::Real re_scaledratio,
+                    const int do_zero_uperpB_driftframe,
 #endif
                     const amrex::Real dt, amrex::Real * PulsarDiag = nullptr)
 {
@@ -86,14 +87,18 @@ void doParticlePush(const GetParticlePosition& GetPosition,
                                                          Ex, Ey, Ez, Bx,
                                                          By, Bz, qp, m,
 #ifdef PULSAR
-                                                         re_scaledratio,
+                                                         re_scaledratio, do_zero_uperpB_driftframe,
 #endif
                                                          dt);
             }
             else {
                 UpdateMomentumBoris( ux, uy, uz,
                                      Ex, Ey, Ez, Bx,
-                                     By, Bz, qp, m, dt);
+                                     By, Bz, qp, m,
+#ifdef PULSAR
+                                     do_zero_uperpB_driftframe,
+#endif
+                                     dt);
             }
             amrex::ParticleReal x, y, z;
             GetPosition(i, x, y, z);
@@ -107,7 +112,7 @@ void doParticlePush(const GetParticlePosition& GetPosition,
                                                      Ex, Ey, Ez, Bx,
                                                      By, Bz, qp, m,
 #ifdef PULSAR
-                                                     re_scaledratio,
+                                                     re_scaledratio, do_zero_uperpB_driftframe,
 #endif
                                                      dt);
             amrex::ParticleReal x, y, z;
@@ -118,7 +123,11 @@ void doParticlePush(const GetParticlePosition& GetPosition,
     } else if (pusher_algo == ParticlePusherAlgo::Boris) {
         UpdateMomentumBoris( ux, uy, uz,
                              Ex, Ey, Ez, Bx,
-                             By, Bz, qp, m, dt, PulsarDiag);
+                             By, Bz, qp, m,
+#ifdef PULSAR
+                             do_zero_uperpB_driftframe,
+#endif
+                             dt);
         amrex::ParticleReal x, y, z;
         GetPosition(i, x, y, z);
         UpdatePosition(x, y, z, ux, uy, uz, dt );

diff --git a/Source/Particles/Pusher/UpdateMomentumBoris.H b/Source/Particles/Pusher/UpdateMomentumBoris.H
@@ -17,7 +17,11 @@ void UpdateMomentumBoris(
     amrex::ParticleReal& ux, amrex::ParticleReal& uy, amrex::ParticleReal& uz,
     const amrex::ParticleReal Ex, const amrex::ParticleReal Ey, const amrex::ParticleReal Ez,
     const amrex::ParticleReal Bx, const amrex::ParticleReal By, const amrex::ParticleReal Bz,
-    const amrex::ParticleReal q, const amrex::ParticleReal m, const amrex::Real dt, amrex::Real * PulsarDiag = nullptr )
+    const amrex::ParticleReal q, const amrex::ParticleReal m,
+#ifdef PULSAR
+    const int do_zero_uperpB_driftframe,
+#endif
+    const amrex::Real dt, amrex::Real * PulsarDiag = nullptr )
 {
     using namespace amrex::literals;
 
@@ -64,56 +68,57 @@ void UpdateMomentumBoris(
 
 
 #ifdef PULSAR
-    // 1(a) Find drift velocity
-    amrex::Real B_sq = Bx*Bx + By*By + Bz*Bz;
-
-    amrex::Real E_sq = Ex*Ex + Ey*Ey + Ez*Ez;
-
-    amrex::Real EdotB = Ex*Bx + Ey*By + Ez*Bz;
-    amrex::Real EdotB_sq = EdotB*EdotB;
-
-    amrex::Real Eprime_sq = 2._prt * EdotB_sq
-              / ( (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 (do_zero_uperpB_driftframe) {
+        // 1(a) Find drift velocity
+        amrex::Real B_sq = Bx*Bx + By*By + Bz*Bz;
+
+        amrex::Real E_sq = Ex*Ex + Ey*Ey + Ez*Ez;
+
+        amrex::Real EdotB = Ex*Bx + Ey*By + Ez*Bz;
+        amrex::Real EdotB_sq = EdotB*EdotB;
+
+        amrex::Real Eprime_sq = 2._prt * EdotB_sq
+                  / ( (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);
+    }
 #endif
 
 }

diff --git a/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H b/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H
@@ -28,6 +28,7 @@ void UpdateMomentumBorisWithRadiationReaction(
     const amrex::ParticleReal q, const amrex::ParticleReal m,
 #ifdef PULSAR
     const amrex::Real re_scaledratio,
+    const int do_zero_uperpB_driftframe,
 #endif
     const amrex::Real dt )
 {
@@ -46,7 +47,11 @@ void UpdateMomentumBorisWithRadiationReaction(
         ux, uy, uz,
         Ex, Ey, Ez,
         Bx, By, Bz,
-        q, m, dt );
+        q, m,
+#ifdef PULSAR
+        do_zero_uperpB_driftframe,
+#endif
+        dt );
 
     //Estimation of the normalized momentum at intermediate (integer) time
     const amrex::ParticleReal ux_n = (ux+ux_old)*0.5_prt;

diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp
@@ -410,6 +410,7 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt,
 
 #ifdef PULSAR
             amrex::Real re_scaledratio = Pulsar::m_re_scaledratio;
+            int do_zero_uperpB_driftframe = Pulsar::m_do_zero_uperpB_driftframe;
 #endif
 
             amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (long ip)
@@ -441,12 +442,17 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt,
                                                              Byp, Bzp, qp, m,
 #ifdef PULSAR
                                                              re_scaledratio,
+                                                             do_zero_uperpB_driftframe,
 #endif
                                                              dt);
                 } else if (pusher_algo == ParticlePusherAlgo::Boris) {
                     UpdateMomentumBoris( uxpp[ip], uypp[ip], uzpp[ip],
                                          Exp, Eyp, Ezp, Bxp,
-                                         Byp, Bzp, qp, m, dt);
+                                         Byp, Bzp, qp, m,
+#ifdef PULSAR
+                                         do_zero_uperpB_driftframe,
+#endif
+                                         dt);
                 } else if (pusher_algo == ParticlePusherAlgo::Vay) {
                     UpdateMomentumVay( uxpp[ip], uypp[ip], uzpp[ip],
                                        Exp, Eyp, Ezp, Bxp,