Skip to content

Commit

Permalink
fix conflicts
Browse files Browse the repository at this point in the history
  • Loading branch information
@RevathiJambunathan
RevathiJambunathan committed Dec 21, 2023
2 parents 0d780b3 + 2874189 commit 565304c
Show file tree
Hide file tree
Showing 8 changed files with 141 additions and 91 deletions.
4 changes: 2 additions & 2 deletions Source/Particles/MultiParticleContainer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2151,8 +2151,8 @@ MultiParticleContainer::PulsarPairInjection ()
);
continue;
}
if ( !((rad > pulsar_particle_inject_rmin) &&
(rad < pulsar_particle_inject_rmax) ) ) {
if ( ! ( ( rad > pulsar_particle_inject_rmin) &&
( rad < pulsar_particle_inject_rmax) )) {
ZeroInitializeAndSetNegativeID(p_sp1, pa_sp1, ip
#ifdef WARPX_QED
,loc_has_quantum_sync_sp1, p_optical_depth_QSR_sp1
Expand Down
43 changes: 33 additions & 10 deletions Source/Particles/PhysicalParticleContainer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2854,9 +2854,18 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
amrex::Real BC_width = Pulsar::m_BC_width;
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;
amrex::Real rmax_zero_uperp_driftframe = Pulsar::m_rmax_zero_uperpB_driftframe;
amrex::Real re_scaledratio = 1.;
amrex::Real scaledRR_rmax = Pulsar::m_scaledRR_rmax;
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;
}
#endif

#ifdef AMREX_USE_OMP
Expand Down Expand Up @@ -3021,7 +3030,10 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
Exp, Eyp, Ezp, Bxp,
Byp, Bzp, qp, m,
#ifdef PULSAR
re_scaledratio, do_zero_uperpB_driftframe, r_p, rmax_zero_uperp_driftframe,
re_scaledratio, scaledRR_rmax,
do_zero_uperpB_driftframe, r_p,
theta_p, dtheta, chi_data,
rmax_zero_uperp_driftframe,
#endif
dt);
} else if (pusher_algo == ParticlePusherAlgo::Boris) {
Expand All @@ -3031,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 @@ -3169,9 +3182,6 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti,
amrex::Real BC_width = Pulsar::m_BC_width;
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;
amrex::Real rmax_zero_uperp_driftframe = Pulsar::m_rmax_zero_uperpB_driftframe;

const int i_EnterInStar = particle_comps["EnterInStar"];
const int i_EnterInMag = particle_comps["EnterInMag"];
Expand All @@ -3184,6 +3194,16 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti,
amrex::ParticleReal* runtimeattr_TimeEnterInStar = pti.GetAttribs(i_TimeEnterInStar ).dataPtr() + offset;
amrex::ParticleReal* runtimeattr_TimeEnterInMag = pti.GetAttribs(i_TimeEnterInMag ).dataPtr() + offset;
amrex::ParticleReal* runtimeattr_TimeReenterInMag = pti.GetAttribs(i_TimeReenterInMag).dataPtr() + offset;
amrex::Real re_scaledratio = 1.;
amrex::Real scaledRR_rmax = Pulsar::m_scaledRR_rmax;
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;
}
#endif


Expand Down Expand Up @@ -3436,7 +3456,8 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti,
t_chi_max,
#endif
#ifdef PULSAR
re_scaledratio, do_zero_uperpB_driftframe, r_p, rmax_zero_uperp_driftframe,
re_scaledratio, scaledRR_rmax,
do_zero_uperpB_driftframe, r_p, theta_p, dtheta, chi_data, rmax_zero_uperp_driftframe,
#endif
dt);
// PulsarPartDiagData);
Expand Down Expand Up @@ -3476,7 +3497,9 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti,
m, q, pusher_algo, do_crr, do_copy,
t_chi_max,
#ifdef PULSAR
re_scaledratio, do_zero_uperpB_driftframe, r_p, rmax_zero_uperp_driftframe,
re_scaledratio, scaledRR_rmax,
do_zero_uperpB_driftframe, r_p,
theta_p, dtheta, chi_data, rmax_zero_uperp_driftframe,
#endif
dt);
}
Expand Down
4 changes: 4 additions & 0 deletions Source/Particles/PulsarParameters.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -713,8 +713,12 @@ public:
static amrex::Real m_BLC;
static amrex::Real m_beta_rec_RR;
static amrex::Real m_re_scaledratio;
static amrex::Real m_scaledRR_rmax;
static int m_do_zero_uperpB_driftframe;
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
8 changes: 8 additions & 0 deletions Source/Particles/PulsarParameters.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -155,6 +155,10 @@ 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;
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 @@ -408,6 +412,10 @@ 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);
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
26 changes: 12 additions & 14 deletions Source/Particles/Pusher/PushSelector.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,8 +67,12 @@ void doParticlePush(const GetParticlePosition& GetPosition,
#endif
#ifdef PULSAR
const amrex::Real re_scaledratio,
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 All @@ -94,7 +98,9 @@ void doParticlePush(const GetParticlePosition& GetPosition,
Ex, Ey, Ez, Bx,
By, Bz, qp, m,
#ifdef PULSAR
re_scaledratio, do_zero_uperpB_driftframe, r_p,
re_scaledratio, scaledRR_rmax,
do_zero_uperpB_driftframe, r_p,
theta_p, dtheta, chi_data,
rmax_zero_uperp_driftframe,
#endif
dt);
Expand All @@ -105,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 @@ -118,12 +125,6 @@ void doParticlePush(const GetParticlePosition& GetPosition,
#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
{
Expand All @@ -132,7 +133,9 @@ void doParticlePush(const GetParticlePosition& GetPosition,
Ex, Ey, Ez, Bx,
By, Bz, qp, m,
#ifdef PULSAR
re_scaledratio, do_zero_uperpB_driftframe, r_p,
re_scaledratio, scaledRR_rmax,
do_zero_uperpB_driftframe, r_p,
theta_p, dtheta, chi_data,
rmax_zero_uperp_driftframe,
#endif
dt);
Expand All @@ -147,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 All @@ -159,12 +163,6 @@ void doParticlePush(const GetParticlePosition& GetPosition,
#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,
Expand Down
122 changes: 62 additions & 60 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,68 +78,67 @@ 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);
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 );
// 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);

}
}
}
} //theta condition
}
#endif

Expand Down
16 changes: 12 additions & 4 deletions Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,12 @@ void UpdateMomentumBorisWithRadiationReaction(
const amrex::ParticleReal q, const amrex::ParticleReal m,
#ifdef PULSAR
const amrex::Real re_scaledratio,
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 @@ -51,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 Expand Up @@ -108,9 +114,11 @@ void UpdateMomentumBorisWithRadiationReaction(
RRcoeff*(PhysConst::c*(flx_q*By - fly_q*Bx) + bdotE*Ez - coeff*bz_n);

//Update momentum using the RR force
ux += frx*dt;
uy += fry*dt;
uz += frz*dt;
if (r_p < scaledRR_rmax) {
ux += frx*dt;
uy += fry*dt;
uz += frz*dt;
}
}

#endif // WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_BORIS_WITHRR_H_
Loading

0 comments on commit 565304c

Please sign in to comment.