diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 7bda9bc2ac0..836af67d4c2 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2856,6 +2856,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, 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; + amrex::Real rmax_zero_uperp_driftframe = Pulsar::m_rmax_zero_uperpB_driftframe; #endif #ifdef AMREX_USE_OMP @@ -3020,7 +3021,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, Exp, Eyp, Ezp, Bxp, Byp, Bzp, qp, m, #ifdef PULSAR - re_scaledratio, do_zero_uperpB_driftframe, + re_scaledratio, do_zero_uperpB_driftframe, r_p, rmax_zero_uperp_driftframe, #endif dt); } else if (pusher_algo == ParticlePusherAlgo::Boris) { @@ -3030,7 +3031,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, Exp, Eyp, Ezp, Bxp, Byp, Bzp, qp, m, #ifdef PULSAR - do_zero_uperpB_driftframe, + do_zero_uperpB_driftframe, r_p, rmax_zero_uperp_driftframe, #endif dt); } else if (pusher_algo == ParticlePusherAlgo::Vay) { @@ -3170,6 +3171,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, 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; + amrex::Real rmax_zero_uperp_driftframe = Pulsar::m_rmax_zero_uperpB_driftframe; #endif @@ -3421,7 +3423,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, t_chi_max, #endif #ifdef PULSAR - re_scaledratio, do_zero_uperpB_driftframe, + re_scaledratio, do_zero_uperpB_driftframe, r_p, rmax_zero_uperp_driftframe, #endif dt); // PulsarPartDiagData); @@ -3461,7 +3463,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, m, q, pusher_algo, do_crr, do_copy, t_chi_max, #ifdef PULSAR - re_scaledratio, do_zero_uperpB_driftframe, + re_scaledratio, do_zero_uperpB_driftframe, r_p, rmax_zero_uperp_driftframe, #endif dt); } diff --git a/Source/Particles/PulsarParameters.H b/Source/Particles/PulsarParameters.H index 86a03f62d9c..dc0b10292a0 100644 --- a/Source/Particles/PulsarParameters.H +++ b/Source/Particles/PulsarParameters.H @@ -714,6 +714,7 @@ public: static amrex::Real m_beta_rec_RR; static amrex::Real m_re_scaledratio; static int m_do_zero_uperpB_driftframe; + static amrex::Real m_rmax_zero_uperpB_driftframe; private: }; diff --git a/Source/Particles/PulsarParameters.cpp b/Source/Particles/PulsarParameters.cpp index c28d5f86753..1b4e44b8eb0 100644 --- a/Source/Particles/PulsarParameters.cpp +++ b/Source/Particles/PulsarParameters.cpp @@ -154,6 +154,7 @@ 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 = 0; +amrex::Real Pulsar::m_rmax_zero_uperpB_driftframe; Pulsar::Pulsar () { @@ -406,6 +407,7 @@ Pulsar::ReadParameters () { } pp.get("do_zero_uperpB_driftframe",m_do_zero_uperpB_driftframe); amrex::Print() << " do zero uperp B " << m_do_zero_uperpB_driftframe << "\n"; + pp.get("rmax_zero_uperpB_driftframe",m_rmax_zero_uperpB_driftframe); } diff --git a/Source/Particles/Pusher/PushSelector.H b/Source/Particles/Pusher/PushSelector.H index 8808c1462eb..4b85ab68351 100644 --- a/Source/Particles/Pusher/PushSelector.H +++ b/Source/Particles/Pusher/PushSelector.H @@ -68,6 +68,8 @@ void doParticlePush(const GetParticlePosition& GetPosition, #ifdef PULSAR const amrex::Real re_scaledratio, const int do_zero_uperpB_driftframe, + const amrex::ParticleReal r_p, + const amrex::Real rmax_zero_uperp_driftframe, #endif const amrex::Real dt, amrex::Real * PulsarDiag = nullptr) { @@ -87,7 +89,8 @@ void doParticlePush(const GetParticlePosition& GetPosition, Ex, Ey, Ez, Bx, By, Bz, qp, m, #ifdef PULSAR - re_scaledratio, do_zero_uperpB_driftframe, + re_scaledratio, do_zero_uperpB_driftframe, r_p, + rmax_zero_uperp_driftframe, #endif dt); } @@ -96,7 +99,8 @@ void doParticlePush(const GetParticlePosition& GetPosition, Ex, Ey, Ez, Bx, By, Bz, qp, m, #ifdef PULSAR - do_zero_uperpB_driftframe, + do_zero_uperpB_driftframe, r_p, + rmax_zero_uperp_driftframe, #endif dt); } @@ -112,7 +116,8 @@ void doParticlePush(const GetParticlePosition& GetPosition, Ex, Ey, Ez, Bx, By, Bz, qp, m, #ifdef PULSAR - re_scaledratio, do_zero_uperpB_driftframe, + re_scaledratio, do_zero_uperpB_driftframe, r_p, + rmax_zero_uperp_driftframe, #endif dt); amrex::ParticleReal x, y, z; @@ -125,7 +130,8 @@ void doParticlePush(const GetParticlePosition& GetPosition, Ex, Ey, Ez, Bx, By, Bz, qp, m, #ifdef PULSAR - do_zero_uperpB_driftframe, + do_zero_uperpB_driftframe, r_p, + rmax_zero_uperp_driftframe, #endif dt); amrex::ParticleReal x, y, z; diff --git a/Source/Particles/Pusher/UpdateMomentumBoris.H b/Source/Particles/Pusher/UpdateMomentumBoris.H index 36b3b756bc9..a4ab3ff4a3f 100644 --- a/Source/Particles/Pusher/UpdateMomentumBoris.H +++ b/Source/Particles/Pusher/UpdateMomentumBoris.H @@ -20,6 +20,8 @@ void UpdateMomentumBoris( const amrex::ParticleReal q, const amrex::ParticleReal m, #ifdef PULSAR const int do_zero_uperpB_driftframe, + const amrex::ParticleReal r_p, + const amrex::Real rmax_zero_uperpB_driftframe, #endif const amrex::Real dt, amrex::Real * PulsarDiag = nullptr ) { @@ -68,7 +70,7 @@ void UpdateMomentumBoris( #ifdef PULSAR - if (do_zero_uperpB_driftframe) { + if (do_zero_uperpB_driftframe && (r_p < rmax_zero_uperpB_driftframe) ) { // 1(a) Find drift velocity amrex::Real B_sq = Bx*Bx + By*By + Bz*Bz; @@ -83,14 +85,14 @@ void UpdateMomentumBoris( + std::sqrt( ( (B_sq*c2) - E_sq) * ( (B_sq*c2) - E_sq ) + 4._prt * EdotB_sq * c2 ) ); // only if E << B -- else fix with thesis - if ( (B_sq + Eprime_sq/c2) > 0 ) { + if ( (B_sq + Eprime_sq/c2) > 1.e-30 ) { amrex::Real vx_drift = (Ey*Bz - Ez*By)/(B_sq + Eprime_sq/c2); amrex::Real vy_drift = (Ez*Bx - Ex*Bz)/(B_sq + Eprime_sq/c2); amrex::Real vz_drift = (Ex*By - Ey*Bx)/(B_sq + Eprime_sq/c2); // 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 < PhysConst::c) { + if (drift_mag < (0.99 * 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; diff --git a/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H b/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H index 7b36b83e5a4..497177a8e66 100644 --- a/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H +++ b/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H @@ -29,6 +29,8 @@ void UpdateMomentumBorisWithRadiationReaction( #ifdef PULSAR const amrex::Real re_scaledratio, const int do_zero_uperpB_driftframe, + const amrex::ParticleReal r_p, + const amrex::Real rmax_zero_uperpB_driftframe, #endif const amrex::Real dt ) { @@ -49,7 +51,7 @@ void UpdateMomentumBorisWithRadiationReaction( Bx, By, Bz, q, m, #ifdef PULSAR - do_zero_uperpB_driftframe, + do_zero_uperpB_driftframe, r_p, rmax_zero_uperpB_driftframe, #endif dt ); diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 56552c77965..4b898a2b902 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -411,6 +411,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; + amrex::Real rmax_zero_uperp_driftframe = Pulsar::m_rmax_zero_uperpB_driftframe; #endif amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (long ip) @@ -422,6 +423,9 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, amrex::ParticleReal xp, yp, zp; getPosition(ip, xp, yp, zp); +#ifdef PULSAR + amrex::ParticleReal r_p = std::sqrt(xp*xp + yp*yp + zp*zp); +#endif amrex::ParticleReal Exp = 0._prt, Eyp = 0._prt, Ezp = 0._prt; amrex::ParticleReal Bxp = 0._prt, Byp = 0._prt, Bzp = 0._prt; @@ -442,7 +446,8 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, Byp, Bzp, qp, m, #ifdef PULSAR re_scaledratio, - do_zero_uperpB_driftframe, + do_zero_uperpB_driftframe, r_p, + rmax_zero_uperp_driftframe, #endif dt); } else if (pusher_algo == ParticlePusherAlgo::Boris) { @@ -450,7 +455,8 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, Exp, Eyp, Ezp, Bxp, Byp, Bzp, qp, m, #ifdef PULSAR - do_zero_uperpB_driftframe, + do_zero_uperpB_driftframe, r_p, + rmax_zero_uperp_driftframe, #endif dt); } else if (pusher_algo == ParticlePusherAlgo::Vay) {