diff --git a/Source/Particles/PulsarParameters.cpp b/Source/Particles/PulsarParameters.cpp index 364f6f0732a..02650dc3ec9 100644 --- a/Source/Particles/PulsarParameters.cpp +++ b/Source/Particles/PulsarParameters.cpp @@ -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) { @@ -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(num_part_cell); @@ -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); @@ -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 ) { diff --git a/Source/Particles/Pusher/UpdateMomentumBoris.H b/Source/Particles/Pusher/UpdateMomentumBoris.H index acf6a6d90fa..7c47ed32af7 100644 --- a/Source/Particles/Pusher/UpdateMomentumBoris.H +++ b/Source/Particles/Pusher/UpdateMomentumBoris.H @@ -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