From ee910aee73ace8d95a78ebeb4b613fe8fa1fb55b Mon Sep 17 00:00:00 2001 From: RevathiJambunathan Date: Wed, 6 Dec 2023 13:12:33 -0800 Subject: [PATCH] add partilce diags, and also fix AddPairinjection to remove particles within rmin rmax --- Source/Particles/MultiParticleContainer.cpp | 17 +++ .../Particles/PhysicalParticleContainer.cpp | 102 +++++++++++++++++- Source/Particles/PulsarParameters.cpp | 33 +++++- 3 files changed, 148 insertions(+), 4 deletions(-) diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index fed9af4562f..f09f6806ee1 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -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; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 836af67d4c2..1fdc7e583c0 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -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 @@ -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, @@ -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; @@ -3683,12 +3718,22 @@ void PhysicalParticleContainer::PulsarParticleRemoval () { amrex::Gpu::DeviceScalar sumParticles(0); amrex::Gpu::DeviceScalar sumWeight(0.0); + amrex::Gpu::DeviceScalar sumParticlesRemovedAfterEnterStar(0); + amrex::Gpu::DeviceScalar sumParticlesInStar(0); + amrex::Gpu::DeviceScalar sumParticlesEnterInStar(0); + amrex::Gpu::DeviceScalar sumParticlesEnterInMag(0); + amrex::Gpu::DeviceScalar sumParticlesReenterInMag(0); const amrex::Real q = this->charge; amrex::GpuArray 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++) { @@ -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) { @@ -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"; } } diff --git a/Source/Particles/PulsarParameters.cpp b/Source/Particles/PulsarParameters.cpp index 1b4e44b8eb0..0e47d26b295 100644 --- a/Source/Particles/PulsarParameters.cpp +++ b/Source/Particles/PulsarParameters.cpp @@ -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 @@ -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 reduce_ops; amrex::Real cur_time = warpx.gett_new(0); auto ws_r = amrex::ParticleReduce< @@ -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; @@ -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 @@ -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