Skip to content

Commit

Permalink
add theta min max for zerouperpB=0
Browse files Browse the repository at this point in the history
  • Loading branch information
@RevathiJambunathan
RevathiJambunathan committed Dec 21, 2023
1 parent 0a5ceac commit 2874189
Show file tree
Hide file tree
Showing 7 changed files with 95 additions and 61 deletions.
11 changes: 8 additions & 3 deletions Source/Particles/PhysicalParticleContainer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2861,6 +2861,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
}
int do_zero_uperpB_driftframe = 0;
amrex::Real rmax_zero_uperp_driftframe = 0.;
amrex::Real dtheta = Pulsar::m_zerouperp_dtheta;
if (cur_time > Pulsar::m_zerouperp_start_time) {
do_zero_uperpB_driftframe = Pulsar::m_do_zero_uperpB_driftframe;
rmax_zero_uperp_driftframe = Pulsar::m_rmax_zero_uperpB_driftframe;
Expand Down Expand Up @@ -3031,6 +3032,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
#ifdef PULSAR
re_scaledratio, scaledRR_rmax,
do_zero_uperpB_driftframe, r_p,
theta_p, dtheta, chi_data,
rmax_zero_uperp_driftframe,
#endif
dt);
Expand All @@ -3041,7 +3043,8 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
Exp, Eyp, Ezp, Bxp,
Byp, Bzp, qp, m,
#ifdef PULSAR
do_zero_uperpB_driftframe, r_p, rmax_zero_uperp_driftframe,
do_zero_uperpB_driftframe, r_p,
theta_p, dtheta, chi_data, rmax_zero_uperp_driftframe,
#endif
dt);
} else if (pusher_algo == ParticlePusherAlgo::Vay) {
Expand Down Expand Up @@ -3184,6 +3187,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti,
if (cur_time > Pulsar::m_RR_start_time) { re_scaledratio = Pulsar::m_re_scaledratio;}
int do_zero_uperpB_driftframe = 0;
amrex::Real rmax_zero_uperp_driftframe = 0.;
amrex::Real dtheta = Pulsar::m_zerouperp_dtheta;
if (cur_time > Pulsar::m_zerouperp_start_time) {
do_zero_uperpB_driftframe = Pulsar::m_do_zero_uperpB_driftframe;
rmax_zero_uperp_driftframe = Pulsar::m_rmax_zero_uperpB_driftframe;
Expand Down Expand Up @@ -3440,7 +3444,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti,
#endif
#ifdef PULSAR
re_scaledratio, scaledRR_rmax,
do_zero_uperpB_driftframe, r_p, rmax_zero_uperp_driftframe,
do_zero_uperpB_driftframe, r_p, theta_p, dtheta, chi_data, rmax_zero_uperp_driftframe,
#endif
dt);
// PulsarPartDiagData);
Expand Down Expand Up @@ -3481,7 +3485,8 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti,
t_chi_max,
#ifdef PULSAR
re_scaledratio, scaledRR_rmax,
do_zero_uperpB_driftframe, r_p, rmax_zero_uperp_driftframe,
do_zero_uperpB_driftframe, r_p,
theta_p, dtheta, chi_data, rmax_zero_uperp_driftframe,
#endif
dt);
}
Expand Down
1 change: 1 addition & 0 deletions Source/Particles/PulsarParameters.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -718,6 +718,7 @@ public:
static amrex::Real m_rmax_zero_uperpB_driftframe;
static amrex::Real m_RR_start_time;
static amrex::Real m_zerouperp_start_time;
static amrex::Real m_zerouperp_dtheta;
private:
};

Expand Down
2 changes: 2 additions & 0 deletions Source/Particles/PulsarParameters.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,6 +158,7 @@ amrex::Real Pulsar::m_rmax_zero_uperpB_driftframe;
amrex::Real Pulsar::m_RR_start_time;
amrex::Real Pulsar::m_zerouperp_start_time;
amrex::Real Pulsar::m_scaledRR_rmax;
amrex::Real Pulsar::m_zerouperp_dtheta;

Pulsar::Pulsar ()
{
Expand Down Expand Up @@ -414,6 +415,7 @@ Pulsar::ReadParameters () {
pp.get("RR_start_time",m_RR_start_time);
pp.get("zerouperp_start_time",m_zerouperp_start_time);
pp.get("scaledRR_rmax",m_scaledRR_rmax);
pp.get("zerouperp_dtheta",m_zerouperp_dtheta);
}


Expand Down
7 changes: 7 additions & 0 deletions Source/Particles/Pusher/PushSelector.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,9 @@ void doParticlePush(const GetParticlePosition& GetPosition,
const amrex::Real scaledRR_rmax,
const int do_zero_uperpB_driftframe,
const amrex::ParticleReal r_p,
const amrex::ParticleReal theta_p,
const amrex::Real dtheta,
const amrex::Real chi_data,
const amrex::Real rmax_zero_uperp_driftframe,
#endif
const amrex::Real dt, amrex::Real * PulsarDiag = nullptr)
Expand Down Expand Up @@ -97,6 +100,7 @@ void doParticlePush(const GetParticlePosition& GetPosition,
#ifdef PULSAR
re_scaledratio, scaledRR_rmax,
do_zero_uperpB_driftframe, r_p,
theta_p, dtheta, chi_data,
rmax_zero_uperp_driftframe,
#endif
dt);
Expand All @@ -107,6 +111,7 @@ void doParticlePush(const GetParticlePosition& GetPosition,
By, Bz, qp, m,
#ifdef PULSAR
do_zero_uperpB_driftframe, r_p,
theta_p, dtheta, chi_data,
rmax_zero_uperp_driftframe,
#endif
dt);
Expand All @@ -130,6 +135,7 @@ void doParticlePush(const GetParticlePosition& GetPosition,
#ifdef PULSAR
re_scaledratio, scaledRR_rmax,
do_zero_uperpB_driftframe, r_p,
theta_p, dtheta, chi_data,
rmax_zero_uperp_driftframe,
#endif
dt);
Expand All @@ -144,6 +150,7 @@ void doParticlePush(const GetParticlePosition& GetPosition,
By, Bz, qp, m,
#ifdef PULSAR
do_zero_uperpB_driftframe, r_p,
theta_p, dtheta, chi_data,
rmax_zero_uperp_driftframe,
#endif
dt);
Expand Down
122 changes: 65 additions & 57 deletions Source/Particles/Pusher/UpdateMomentumBoris.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,9 @@ void UpdateMomentumBoris(
#ifdef PULSAR
const int do_zero_uperpB_driftframe,
const amrex::ParticleReal r_p,
const amrex::ParticleReal theta_p,
const amrex::Real dtheta,
const amrex::Real chi_data,
const amrex::Real rmax_zero_uperpB_driftframe,
#endif
const amrex::Real dt, amrex::Real * PulsarDiag = nullptr )
Expand Down Expand Up @@ -75,64 +78,69 @@ void UpdateMomentumBoris(

#ifdef PULSAR
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;

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 c2 = PhysConst::c*PhysConst::c;
amrex::Real Eprime_sq = 2._prt * EdotB_sq * c2
/ ( ( (B_sq * c2) - E_sq)
+ 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) > 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) ) {
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);

// debugging
amrex::Real abs_u_post = std::sqrt(ux*ux + uy*uy + uz*uz);
// equator is at 90, chi_data accounts for obliquity
amrex::Real theta_min = (90. + chi_data - dtheta) * MathConst::pi / 180.;
amrex::Real theta_max = (90. + chi_data + dtheta) * MathConst::pi / 180.;
if ( (theta_p > theta_min) && (theta_p < theta_max) ) {
// 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 c2 = PhysConst::c*PhysConst::c;
amrex::Real Eprime_sq = 2._prt * EdotB_sq * c2
/ ( ( (B_sq * c2) - E_sq)
+ 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) > 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) ) {
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);

// debugging
amrex::Real abs_u_post = std::sqrt(ux*ux + uy*uy + uz*uz);
}
}
}
} //theta condition
}
#endif

Expand Down
7 changes: 6 additions & 1 deletion Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,9 @@ void UpdateMomentumBorisWithRadiationReaction(
const amrex::Real scaledRR_rmax,
const int do_zero_uperpB_driftframe,
const amrex::ParticleReal r_p,
const amrex::ParticleReal theta_p,
const amrex::Real dtheta,
const amrex::Real chi_data,
const amrex::Real rmax_zero_uperpB_driftframe,
#endif
const amrex::Real dt )
Expand All @@ -52,7 +55,9 @@ void UpdateMomentumBorisWithRadiationReaction(
Bx, By, Bz,
q, m,
#ifdef PULSAR
do_zero_uperpB_driftframe, r_p, rmax_zero_uperpB_driftframe,
do_zero_uperpB_driftframe, r_p,
theta_p, dtheta, chi_data,
rmax_zero_uperpB_driftframe,
#endif
dt );

Expand Down
6 changes: 6 additions & 0 deletions Source/Particles/RigidInjectedParticleContainer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -413,6 +413,8 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt,
amrex::Real scaledRR_rmax = Pulsar::m_scaledRR_rmax;
int do_zero_uperpB_driftframe = Pulsar::m_do_zero_uperpB_driftframe;
amrex::Real rmax_zero_uperp_driftframe = Pulsar::m_rmax_zero_uperpB_driftframe;
amrex::Real dtheta = Pulsar::m_zerouperp_dtheta;
amrex::Real chi_data = Pulsar::m_Chi;
#endif

amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (long ip)
Expand All @@ -425,7 +427,9 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt,
getPosition(ip, xp, yp, zp);

#ifdef PULSAR
// RigidParticleContainer not used for Pulsars.
amrex::ParticleReal r_p = std::sqrt(xp*xp + yp*yp + zp*zp);
amrex::ParticleReal theta_p = 0.;
#endif
amrex::ParticleReal Exp = 0._prt, Eyp = 0._prt, Ezp = 0._prt;
amrex::ParticleReal Bxp = 0._prt, Byp = 0._prt, Bzp = 0._prt;
Expand All @@ -448,6 +452,7 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt,
#ifdef PULSAR
re_scaledratio, scaledRR_rmax,
do_zero_uperpB_driftframe, r_p,
theta_p, dtheta, chi_data,
rmax_zero_uperp_driftframe,
#endif
dt);
Expand All @@ -457,6 +462,7 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt,
Byp, Bzp, qp, m,
#ifdef PULSAR
do_zero_uperpB_driftframe, r_p,
theta_p, dtheta, chi_data,
rmax_zero_uperp_driftframe,
#endif
dt);
Expand Down

0 comments on commit 2874189

Please sign in to comment.