Skip to content

Commit

Permalink
add partilce diags, and also fix AddPairinjection to remove particles…
Browse files Browse the repository at this point in the history 
… within rmin rmax
  • Loading branch information
@RevathiJambunathan
RevathiJambunathan committed Dec 6, 2023
1 parent a17c62e commit ee910ae
Show file tree
Hide file tree
Showing 3 changed files with 148 additions and 4 deletions.
17 changes: 17 additions & 0 deletions Source/Particles/MultiParticleContainer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2149,6 +2149,23 @@ MultiParticleContainer::PulsarPairInjection ()
,loc_has_breit_wheeler_sp2, p_optical_depth_BW_sp2
#endif
);
continue;
}
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
,loc_has_breit_wheeler_sp1, p_optical_depth_BW_sp1
#endif
);
ZeroInitializeAndSetNegativeID(p_sp2, pa_sp2, ip
#ifdef WARPX_QED
,loc_has_quantum_sync_sp2, p_optical_depth_QSR_sp2
,loc_has_breit_wheeler_sp2, p_optical_depth_BW_sp2
#endif
);
continue;
}
// Initialize particle velocity along the Bfield line
XDim3 u;
Expand Down
102 changes: 99 additions & 3 deletions Source/Particles/PhysicalParticleContainer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3172,6 +3172,18 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti,
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"];
const int i_TimeEnterInStar = particle_comps["TimeEnterInStar"];
const int i_TimeEnterInMag = particle_comps["TimeEnterInMag"];
const int i_TimeReenterInMag = particle_comps["TimeReenterInMag"];

amrex::ParticleReal* runtimeattr_EnterInStar = pti.GetAttribs(i_EnterInStar ).dataPtr() + offset;
amrex::ParticleReal* runtimeattr_EnterInMag = pti.GetAttribs(i_EnterInMag ).dataPtr() + offset;
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;
#endif


Expand Down Expand Up @@ -3320,6 +3332,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti,
amrex::ParticleReal r_p, theta_p, phi_p;
Pulsar::ConvertCartesianToSphericalCoord( xp, yp, zp, center_star_arr,
r_p, theta_p, phi_p);
amrex::ParticleReal rp_old = r_p;
Pulsar::getExternalEBOnParticle(r_p, theta_p, phi_p, chi_data,
AddExternalMonopoleOnly,
AddPulsarVacuumEFields,
Expand Down Expand Up @@ -3470,6 +3483,28 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti,
}
#endif

#ifdef PULSAR
getPosition(ip, xp, yp, zp);
Pulsar::ConvertCartesianToSphericalCoord( xp, yp, zp, center_star_arr,
r_p, theta_p, phi_p);
amrex::Real rp_new = r_p;
if ( (rp_old > Rstar_data) && (rp_new < Rstar_data) ) {
// Particle has entered Rstar. SetFlag , EnterInStar = 1, and note time
runtimeattr_EnterInStar[ip] = 1.;
runtimeattr_TimeEnterInStar[ip] = cur_time;
}
if ( (rp_old < Rstar_data) && (rp_new > Rstar_data) ) {
// Particle has entered magnetosphere. Set Flag to EnterInMag = 1, and note time, TimeEnterInMag
runtimeattr_EnterInMag[ip] = 1.;
runtimeattr_TimeEnterInMag[ip] = cur_time;
if ( (runtimeattr_EnterInStar[ip] < 1.1) && (runtimeattr_EnterInStar[ip] > 0.9) ) {
// check if particle has a flag from previous time, EnterInStar == 1,
// then this particle is reentering mag, note time for re-entering mag TimeReenterInMag
runtimeattr_TimeReenterInMag[ip] = cur_time;
}
}
#endif

#ifdef WARPX_QED
auto foo_local_has_quantum_sync = local_has_quantum_sync;
auto foo_podq = p_optical_depth_QSR;
Expand Down Expand Up @@ -3683,12 +3718,22 @@ void PhysicalParticleContainer::PulsarParticleRemoval ()
{
amrex::Gpu::DeviceScalar<int> sumParticles(0);
amrex::Gpu::DeviceScalar<amrex::Real> sumWeight(0.0);
amrex::Gpu::DeviceScalar<int> sumParticlesRemovedAfterEnterStar(0);
amrex::Gpu::DeviceScalar<int> sumParticlesInStar(0);
amrex::Gpu::DeviceScalar<int> sumParticlesEnterInStar(0);
amrex::Gpu::DeviceScalar<int> sumParticlesEnterInMag(0);
amrex::Gpu::DeviceScalar<int> sumParticlesReenterInMag(0);
const amrex::Real q = this->charge;
amrex::GpuArray<amrex::Real, 3> center_star_arr;
for (int i_dim = 0; i_dim < 3; ++i_dim) {
center_star_arr[i_dim] = Pulsar::m_center_star[i_dim];
}
amrex::Real max_particle_absorption_radius = Pulsar::m_max_particle_absorption_radius;
amrex::Real Rstar_data = Pulsar::m_R_star;
auto& warpx = WarpX::GetInstance();
amrex::Real cur_time = warpx.gett_new(0);
amrex::Real dt = warpx.getdt(0);

// Remove Particles From inside sphere
for (int lev = 0; lev <= finestLevel(); lev++)
{
Expand All @@ -3705,8 +3750,23 @@ void PhysicalParticleContainer::PulsarParticleRemoval ()
auto& attribs = pti.GetAttribs();
auto& wp = attribs[PIdx::w];
int* pSum_d = sumParticles.dataPtr();
int* pSumRemovedEnterInStar_d = sumParticlesRemovedAfterEnterStar.dataPtr();
int* pSumInStar_d = sumParticlesInStar.dataPtr();
int* pSumEnterInStar_d = sumParticlesEnterInStar.dataPtr();
int* pSumEnterInMag_d = sumParticlesEnterInMag.dataPtr();
int* pSumReenterInMag_d = sumParticlesReenterInMag.dataPtr();
amrex::Real* wSum_d = sumWeight.dataPtr();
amrex::Real* wp_d = wp.dataPtr();
const int i_EnterInStar = particle_comps["EnterInStar"];
const int i_EnterInMag = particle_comps["EnterInMag"];
const int i_TimeEnterInStar = particle_comps["TimeEnterInStar"];
const int i_TimeEnterInMag = particle_comps["TimeEnterInMag"];
const int i_TimeReenterInMag = particle_comps["TimeReenterInMag"];
amrex::ParticleReal* runtime_EnterInStar = pti.GetAttribs(i_EnterInStar).dataPtr();
amrex::ParticleReal* runtime_EnterInMag = pti.GetAttribs(i_EnterInMag).dataPtr();
amrex::ParticleReal* runtime_TimeReenterInMag = pti.GetAttribs(i_TimeReenterInMag).dataPtr();
amrex::ParticleReal* runtime_TimeEnterInStar = pti.GetAttribs(i_TimeEnterInStar).dataPtr();
amrex::ParticleReal* runtime_TimeEnterInMag = pti.GetAttribs(i_TimeEnterInMag).dataPtr();

amrex::ParallelFor(pti.numParticles(),
[=] AMREX_GPU_DEVICE (long i) {
Expand All @@ -3715,22 +3775,58 @@ void PhysicalParticleContainer::PulsarParticleRemoval ()
Real r = std::sqrt((x-xc)*(x-xc)
+ (y-yc)*(y-yc)
+ (z-zc)*(z-zc));
// flag particles to be deleted and count them
if (r <= (max_particle_absorption_radius)) {
pp[i].id() = -1;
// atomic add
int const unity = 1;
amrex::Gpu::Atomic::AddNoRet(pSum_d,unity);
amrex::Gpu::Atomic::AddNoRet(wSum_d,wp_d[i]);

// check if flag for EnterInStar==1, and count how many deleted particles are from entering
if ( (runtime_EnterInStar[i] < 1.1) && (runtime_EnterInStar[i] > 0.9)) {
int const one = 1;
amrex::Gpu::Atomic::AddNoRet(pSumRemovedEnterInStar_d, one);
}
}
// Total Particles In Star
if (r < Rstar_data) {
int const unity = 1;
amrex::Gpu::Atomic::AddNoRet(pSumInStar_d,unity);
}

//Count particles Entering Star in this timestep
if ( (runtime_EnterInStar[i] < 1.1) && (runtime_EnterInStar[i] > 0.9)) {
// Because cur time has been updated by the time we do this diagnostic
if ( (runtime_TimeEnterInStar[i] < (cur_time - dt + 0.1 * dt) ) &&
(runtime_TimeEnterInStar[i] > (cur_time - dt - 0.1 * dt) )) {
int const unity = 1;
amrex::Gpu::Atomic::AddNoRet(pSumEnterInStar_d,unity);
}
}

// Counter particles entering magnetosphere in this timestep
if ( (runtime_EnterInMag[i] < 1.1) && (runtime_EnterInMag[i] > 0.9) ) {
if ( (runtime_TimeEnterInMag[i] < (cur_time - dt + 0.1 * dt) ) &&
(runtime_TimeEnterInMag[i] > (cur_time - dt - 0.1 * dt) )) {
int const unity = 1;
amrex::Gpu::Atomic::AddNoRet(pSumEnterInMag_d,unity);
}
// How many particle entering the magnetosphere are re-entering
if ( (runtime_TimeReenterInMag[i] < (cur_time - dt + 0.1 * dt) ) &&
(runtime_TimeReenterInMag[i] > (cur_time - dt - 0.1 * dt) )) {
int const unity = 1;
amrex::Gpu::Atomic::AddNoRet(pSumReenterInMag_d,unity);
}
}
});
}
amrex::Gpu::synchronize();
}
auto& warpx = WarpX::GetInstance();
if (q > 0) {
amrex::AllPrintToFile("TotalPositronsEnteringStar") << warpx.gett_new(0) << " " << sumParticles.dataValue() << " " << sumWeight.dataValue() << "\n";
amrex::AllPrintToFile("TotalPositronsDeletedInStar") << warpx.getistep(0) << " " << warpx.gett_new(0) << " " << sumParticles.dataValue() << " " << sumWeight.dataValue() << " " << sumParticlesRemovedAfterEnterStar.dataValue() << " " << sumParticlesInStar.dataValue() << " " << sumParticlesEnterInStar.dataValue() << " " << sumParticlesEnterInMag.dataValue() << " " << sumParticlesReenterInMag.dataValue() << "\n";
} else {
amrex::AllPrintToFile("TotalElectronsEnteringStar") << warpx.gett_new(0) << " " << sumParticles.dataValue() << " " << sumWeight.dataValue() << "\n";
amrex::AllPrintToFile("TotalElectronsDeletedInStar") << warpx.getistep(0) << " "<< warpx.gett_new(0) << " " << sumParticles.dataValue() << " " << sumWeight.dataValue() << " " << sumParticlesRemovedAfterEnterStar.dataValue() << " " << sumParticlesInStar.dataValue() << " " << sumParticlesEnterInStar.dataValue() << " " << sumParticlesEnterInMag.dataValue() << " " << sumParticlesReenterInMag.dataValue() << "\n";
}

}
Expand Down
33 changes: 32 additions & 1 deletion Source/Particles/PulsarParameters.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1601,6 +1601,8 @@ Pulsar::TuneSigma0Threshold (const int step)
const int nspecies = species_names.size();
amrex::Real total_weight_allspecies = 0._rt;
amrex::Real dt = warpx.getdt(0);
amrex::Real NpSp0_injected = 0._rt;
amrex::Real NpSp1_injected = 0._rt;
// Total number of cells that have injected particles
// amrex::Real total_injected_cells = SumInjectionFlag();
// // for debugging print injected cell data
Expand All @@ -1614,6 +1616,8 @@ Pulsar::TuneSigma0Threshold (const int step)
for (int isp = 0; isp < nspecies; ++isp) {
amrex::Real ws_total = 0._rt;
auto& pc = warpx.GetPartContainer().GetParticleContainer(isp);
auto pcomps = pc.getParticleComps();
const int i_injectiontime = pcomps["injectiontime"] - PIdx::nattribs;
amrex::ReduceOps<amrex::ReduceOpSum> reduce_ops;
amrex::Real cur_time = warpx.gett_new(0);
auto ws_r = amrex::ParticleReduce<
Expand All @@ -1624,7 +1628,7 @@ Pulsar::TuneSigma0Threshold (const int step)
auto p = ptd.getSuperParticle(i);
amrex::ParticleReal wp = p.rdata(PIdx::w);
amrex::Real filter = 0._rt;
amrex::ParticleReal injectiontime = ptd.m_runtime_rdata[0][i];
amrex::ParticleReal injectiontime = ptd.m_runtime_rdata[i_injectiontime][i];
if ( injectiontime < cur_time + 0.1_rt*dt and
injectiontime > cur_time - 0.1_rt*dt ) {
filter = 1._rt;
Expand All @@ -1635,6 +1639,32 @@ Pulsar::TuneSigma0Threshold (const int step)
);
ws_total = amrex::get<0>(ws_r);
amrex::ParallelDescriptor::ReduceRealSum(ws_total);

amrex::Real Ntotal = 0._rt;
auto Ns_r = amrex::ParticleReduce<
amrex::ReduceData < amrex::ParticleReal> >
( pc,
[=] AMREX_GPU_DEVICE (const PTDType &ptd, const int i) noexcept
{
auto p = ptd.getSuperParticle(i);
amrex::ParticleReal wp = p.rdata(PIdx::w);
amrex::Real filter = 0._rt;
amrex::ParticleReal injectiontime = ptd.m_runtime_rdata[i_injectiontime][i];
if ( injectiontime < cur_time + 0.1_rt*dt and
injectiontime > cur_time - 0.1_rt*dt ) {
filter = 1._rt;
}
return (filter);
},
reduce_ops
);
Ntotal = amrex::get<0>(Ns_r);
amrex::ParallelDescriptor::ReduceRealSum(Ntotal);
if (isp == 0) {
NpSp0_injected = Ntotal;
} else {
NpSp1_injected = Ntotal;
}
total_weight_allspecies += ws_total;
}
// injection rate is sum of particle weight over all species per timestep
Expand Down Expand Up @@ -1825,6 +1855,7 @@ Pulsar::TuneSigma0Threshold (const int step)
// amrex::Real Btp = m_B_star * (c_chi * s_theta - s_chi * c_theta * c_psi);
// amrex::Real Bphip = m_B_star * s_chi * s_psi;
amrex::AllPrintToFile("ROI") << warpx.getistep(0) << " " << warpx.gett_new(0) << " " << dt << " " << specified_injection_rate << " " << current_injection_rate <<"\n";
amrex::AllPrintToFile("NpInjectedInStar") << warpx.getistep(0) << " " << warpx.gett_new(0) << " " << dt << " " << NpSp0_injected << " " << NpSp1_injected << "\n";
}

amrex::Real
Expand Down

0 comments on commit ee910ae

Please sign in to comment.