Skip to content

Commit

Permalink
limiting drift to c
Browse files Browse the repository at this point in the history
  • Loading branch information
@RevathiJambunathan
RevathiJambunathan committed Oct 26, 2023
1 parent b03ffa1 commit 6040ddc
Show file tree
Hide file tree
Showing 2 changed files with 52 additions and 43 deletions.
23 changes: 15 additions & 8 deletions Source/Particles/PulsarParameters.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2329,6 +2329,7 @@ Pulsar::FlagCellsForInjectionWithPcounts ()
amrex::Real sum_magnetization = m_sum_inj_magnetization;
amrex::Real GJdensitythreshold = m_injection_GJdensitythreshold;
// fill pcounts and injected cell flag
amrex::Real BulkVel = m_part_bulkVelocity;
amrex::Real chi = m_Chi;
for (amrex::MFIter mfi(*m_injection_flag[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
Expand Down Expand Up @@ -2424,6 +2425,7 @@ Pulsar::FlagCellsForInjectionWithPcounts ()
//amrex::Real num_part_cell = num_ppc_modified_real * factor;
if (use_injection_rate == 1) {
amrex::Real GJ_inj_rate = n_GJ * dx_lev[0] * dx_lev[0] * 3.e8;
if (BulkVel > 0.05) GJ_inj_rate = n_GJ * dx_lev[0] * dx_lev[0] * 3.e8 * BulkVel;
num_part_cell = GJ_inj_rate * injection_fac * dt / particle_wt;
}
int numpart_int = static_cast<int>(num_part_cell);
Expand Down Expand Up @@ -2465,12 +2467,12 @@ Pulsar::FlagCellsForInjectionWithPcounts ()
// pcount(i,j,k) = num_part + 1;
// }
//}
if (density_thresholdfactor >= 0 ) {
if ( amrex::Math::abs(rho_GJ) < (density_thresholdfactor * rho_GJ_fac) ) {
pcount(i,j,k) = 0;
injected_cell(i,j,k) = 0;
}
}
//if (density_thresholdfactor >= 0 ) {
// if ( amrex::Math::abs(rho_GJ) < (density_thresholdfactor * rho_GJ_fac) ) {
// pcount(i,j,k) = 0;
// injected_cell(i,j,k) = 0;
// }
//}
//if (density_thresholdfactor < 0 ) {
// if ( amrex::Math::abs(rho_GJ) < ( amrex::Math::abs(density_thresholdfactor) * rho_GJ_fac) ) {
// rho_GJ = rho_GJ_fac * amrex::Math::abs(density_thresholdfactor);
Expand Down Expand Up @@ -2540,8 +2542,13 @@ Pulsar::FlagCellsForInjectionWithPcounts ()
ConvertCartesianToSphericalCoord(x, y, z, xc,
r, theta, phi);
amrex::Real q = 1.609e-19;
amrex::Real rho_GJ = rho_GJ_fac * (1. - 3. * std::cos(theta) * std::cos(theta) );
amrex::Real n_GJ = amrex::Math::abs(rho_GJ)/q;
amrex::Real shifted_theta = theta - chi;
amrex::Real GJ_factor = amrex::Math::abs( 1. - 3. * std::cos(shifted_theta)*std::cos(shifted_theta));
if (GJ_factor < GJdensitythreshold) GJ_factor = GJdensitythreshold;
amrex::Real rho_GJ = rho_GJ_fac * GJ_factor;
amrex::Real n_GJ = amrex::Math::abs(rho_GJ)/q;
//amrex::Real rho_GJ = rho_GJ_fac * (1. - 3. * std::cos(theta) * std::cos(theta) );
//amrex::Real n_GJ = amrex::Math::abs(rho_GJ)/q;
//amrex::Real n = amrex::Math::abs(rho(i,j,k))/q;
amrex::Real n = ndens(i,j,k,0)+ndens(i,j,k,1);
if ( (n > (limit_GJ_factor * n_GJ) ) && limit_injection == 1 ) {
Expand Down
72 changes: 37 additions & 35 deletions Source/Particles/Pusher/UpdateMomentumBoris.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,44 +81,46 @@ void UpdateMomentumBoris(
/ ( (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
if ( (B_sq + Eprime_sq) > 0. ) {
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);
if ( (B_sq > 0. ) && ( std::sqrt(Eprime_sq) < std::sqrt(B_sq) ) ) {
amrex::Real vx_drift = (Ey*Bz - Ez*By)/(B_sq);
amrex::Real vy_drift = (Ez*Bx - Ex*Bz)/(B_sq);
amrex::Real vz_drift = (Ex*By - Ey*Bx)/(B_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 (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);
}
}
}
#endif
Expand Down

0 comments on commit 6040ddc

Please sign in to comment.