From 7b9aa6dd3108f10ce088614abccb239645a5dfe0 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Sat, 3 Aug 2019 17:27:40 -0700 Subject: add ionization_level component --- Source/Particles/WarpXParticleContainer.cpp | 41 ++++++++++++++++------------- 1 file changed, 23 insertions(+), 18 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 89f233b2c..bce86a925 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -237,12 +237,14 @@ WarpXParticleContainer::AddNParticles (int lev, p.pos(1) = z[i]; #endif - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - auto& ptile = DefineAndReturnParticleTile(0, 0, 0); - ptile.push_back_real(particle_comps["xold"], x[i]); - ptile.push_back_real(particle_comps["yold"], y[i]); - ptile.push_back_real(particle_comps["zold"], z[i]); + if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ){ + auto& ptile = DefineAndReturnParticleTile(0, 0, 0); + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + { + ptile.push_back_real(particle_comps["xold"], x[i]); + ptile.push_back_real(particle_comps["yold"], y[i]); + ptile.push_back_real(particle_comps["zold"], z[i]); + } } particle_tile.push_back(p); @@ -255,12 +257,14 @@ WarpXParticleContainer::AddNParticles (int lev, particle_tile.push_back_real(PIdx::uy, vy + ibegin, vy + iend); particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend); - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - auto& ptile = DefineAndReturnParticleTile(0, 0, 0); - ptile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); - ptile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); - ptile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); + if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ){ + auto& ptile = DefineAndReturnParticleTile(0, 0, 0); + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + { + ptile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); + ptile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); + ptile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); + } } for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) @@ -432,6 +436,7 @@ void WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, RealVector& wp, RealVector& uxp, RealVector& uyp, RealVector& uzp, + Real* ion_lev, MultiFab* jx, MultiFab* jy, MultiFab* jz, const long offset, const long np_to_depose, int thread_num, int lev, int depos_lev, @@ -514,34 +519,34 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Esirkepov) { if (WarpX::nox == 1){ doEsirkepovDepositionShapeN<1>(xp, yp, zp, wp.dataPtr(), uxp.dataPtr(), - uyp.dataPtr(), uzp.dataPtr(), jx_arr, jy_arr, + uyp.dataPtr(), uzp.dataPtr(), ion_lev, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q); } else if (WarpX::nox == 2){ doEsirkepovDepositionShapeN<2>(xp, yp, zp, wp.dataPtr(), uxp.dataPtr(), - uyp.dataPtr(), uzp.dataPtr(), jx_arr, jy_arr, + uyp.dataPtr(), uzp.dataPtr(), ion_lev, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q); } else if (WarpX::nox == 3){ doEsirkepovDepositionShapeN<3>(xp, yp, zp, wp.dataPtr(), uxp.dataPtr(), - uyp.dataPtr(), uzp.dataPtr(), jx_arr, jy_arr, + uyp.dataPtr(), uzp.dataPtr(), ion_lev, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q); } } else { if (WarpX::nox == 1){ doDepositionShapeN<1>(xp, yp, zp, wp.dataPtr(), uxp.dataPtr(), - uyp.dataPtr(), uzp.dataPtr(), jx_arr, jy_arr, + uyp.dataPtr(), uzp.dataPtr(), ion_lev, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, stagger_shift, q); } else if (WarpX::nox == 2){ doDepositionShapeN<2>(xp, yp, zp, wp.dataPtr(), uxp.dataPtr(), - uyp.dataPtr(), uzp.dataPtr(), jx_arr, jy_arr, + uyp.dataPtr(), uzp.dataPtr(), ion_lev, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, stagger_shift, q); } else if (WarpX::nox == 3){ doDepositionShapeN<3>(xp, yp, zp, wp.dataPtr(), uxp.dataPtr(), - uyp.dataPtr(), uzp.dataPtr(), jx_arr, jy_arr, + uyp.dataPtr(), uzp.dataPtr(), ion_lev, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, stagger_shift, q); } -- cgit v1.2.3 From c7a56a07ee33f5d1cce29719a946036e36b5897c Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 8 Aug 2019 11:00:08 -0700 Subject: change ionization level from Real to Int --- Source/Diagnostics/ParticleIO.cpp | 3 ++- Source/Laser/LaserParticleContainer.cpp | 2 +- Source/Particles/Deposition/CurrentDeposition.H | 4 ++-- Source/Particles/PhysicalParticleContainer.cpp | 29 +++++++++++++------------ Source/Particles/WarpXParticleContainer.H | 9 ++++++-- Source/Particles/WarpXParticleContainer.cpp | 2 +- 6 files changed, 28 insertions(+), 21 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index 2b70c9379..618ecc6c5 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -114,7 +114,8 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const } if(pc->do_field_ionization){ - real_names.push_back("ionization_level"); + int_names.push_back("ionization_level"); + int_flags.resize(1, 1); } // Convert momentum to SI diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 37e8f8bf2..e2aebdd5e 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -504,7 +504,7 @@ LaserParticleContainer::Evolve (int lev, // Current Deposition // // Deposit inside domains - Real* ion_lev = nullptr; + int* ion_lev = nullptr; DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz, 0, np_current, thread_num, lev, lev, dt); diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index deff65424..86cf89cd3 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -26,7 +26,7 @@ void doDepositionShapeN(const amrex::Real * const xp, const amrex::Real * const uxp, const amrex::Real * const uyp, const amrex::Real * const uzp, - const amrex::Real * const ion_lev, + const int * const ion_lev, const amrex::Array4& jx_arr, const amrex::Array4& jy_arr, const amrex::Array4& jz_arr, @@ -184,7 +184,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, const amrex::Real * const uxp, const amrex::Real * const uyp, const amrex::Real * const uzp, - const amrex::Real * ion_lev, + const int * ion_lev, const amrex::Array4& Jx_arr, const amrex::Array4& Jy_arr, const amrex::Array4& Jz_arr, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 7fdb783b7..98d9e5d18 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -452,9 +452,9 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) pb[4] = soa.GetRealData(particle_comps["uyold"]).data() + old_size; pb[5] = soa.GetRealData(particle_comps["uzold"]).data() + old_size; } - Real* pi; + int* pi; if (do_field_ionization) { - pi = soa.GetRealData(particle_comps[ "ionization_level"]).data() + old_size; + pi = soa.GetIntData(particle_icomps["ionization_level"]).data() + old_size; } const GpuArray overlap_corner @@ -1293,9 +1293,9 @@ PhysicalParticleContainer::Evolve (int lev, lev, lev-1, dt); } } else { - Real* AMREX_RESTRICT ion_lev; + int* AMREX_RESTRICT ion_lev; if (do_field_ionization){ - ion_lev = pti.GetAttribs(particle_comps["ionization_level"]).dataPtr(); + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); } else { ion_lev = nullptr; } @@ -1523,9 +1523,9 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, copy_attribs(pti, x, y, z); } - Real* AMREX_RESTRICT ion_lev = nullptr; + int* AMREX_RESTRICT ion_lev = nullptr; if (do_field_ionization){ - ion_lev = pti.GetAttribs(particle_comps["ionization_level"]).dataPtr(); + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); } // Loop over the particles and update their momentum @@ -1972,7 +1972,7 @@ void PhysicalParticleContainer::InitIonizationModule() pp.get("ionization_product", ionization_product_name); pp.get("physical_element", physical_element); // Add Real component for ionization level - AddRealComp("ionization_level"); + AddIntComp("ionization_level"); plot_flags.resize(PIdx::nattribs + 1, 1); // Get atomic number and ionization energies from file int ion_element_id = ion_map_ids[physical_element]; @@ -2051,7 +2051,8 @@ PhysicalParticleContainer::buildIonizationMask(const amrex::MFIter& mfi, const i const Real * const AMREX_RESTRICT bx = soa.GetRealData(PIdx::Bx).data(); const Real * const AMREX_RESTRICT by = soa.GetRealData(PIdx::By).data(); const Real * const AMREX_RESTRICT bz = soa.GetRealData(PIdx::Bz).data(); - Real* ilev_real = soa.GetRealData(particle_comps["ionization_level"]).data(); + // Real* ilev_real = soa.GetIntData(particle_icomps["ionization_level"]).data(); + int* ion_lev = soa.GetIntData(particle_icomps["ionization_level"]).data(); Real c = PhysConst::c; Real c2_inv = 1./c/c; @@ -2063,9 +2064,9 @@ PhysicalParticleContainer::buildIonizationMask(const amrex::MFIter& mfi, const i np, [=] AMREX_GPU_DEVICE (long i) { // Get index of ionization_level - int ilev = (int) round(ilev_real[i]); + // int ilev = (int) round(ilev_real[i]); p_ionization_mask[i] = 0; - if ( ilev particle_comps; + std::map particle_icomps; int species_id; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 2e418aae0..a7c30764d 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -428,7 +428,7 @@ void WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, RealVector& wp, RealVector& uxp, RealVector& uyp, RealVector& uzp, - Real* ion_lev, + int* ion_lev, MultiFab* jx, MultiFab* jy, MultiFab* jz, const long offset, const long np_to_depose, int thread_num, int lev, int depos_lev, -- cgit v1.2.3 From 8dcbb28a52f429e22a6643ffa33e81c6224942e0 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 8 Aug 2019 11:31:22 -0700 Subject: add some const, and remove unnecessary includes --- Source/Particles/MultiParticleContainer.H | 1 - Source/Particles/MultiParticleContainer.cpp | 5 ----- Source/Particles/WarpXParticleContainer.H | 2 +- Source/Particles/WarpXParticleContainer.cpp | 2 +- 4 files changed, 2 insertions(+), 8 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 8f7fd56e3..fd9a33e86 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -7,7 +7,6 @@ #include #include -#include #include #include #include diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 6d977e16c..bff369edc 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -1,15 +1,10 @@ #include #include #include -#include -#include -#include -#include #include #include #include -#include using namespace amrex; diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 2613be167..dd5a150c4 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -184,7 +184,7 @@ public: RealVector& uxp, RealVector& uyp, RealVector& uzp, - int* ion_lev, + const int * const ion_lev, amrex::MultiFab* jx, amrex::MultiFab* jy, amrex::MultiFab* jz, diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index a7c30764d..eba072572 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -428,7 +428,7 @@ void WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, RealVector& wp, RealVector& uxp, RealVector& uyp, RealVector& uzp, - int* ion_lev, + const int * const ion_lev, MultiFab* jx, MultiFab* jy, MultiFab* jz, const long offset, const long np_to_depose, int thread_num, int lev, int depos_lev, -- cgit v1.2.3 From 16dd4d0491b32fb669a1ee07ed2866f68aafd6da Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 21 Aug 2019 09:14:21 -0700 Subject: charge deposition depends on ionization level --- Source/Laser/LaserParticleContainer.cpp | 12 ++++++++---- Source/Particles/Deposition/ChargeDeposition.H | 13 ++++++++++++- Source/Particles/PhysicalParticleContainer.cpp | 25 +++++++++++++++++++++---- Source/Particles/WarpXParticleContainer.H | 1 + Source/Particles/WarpXParticleContainer.cpp | 23 ++++++++++++++++------- 5 files changed, 58 insertions(+), 16 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index e9a2d42de..9fac71894 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -454,9 +454,10 @@ LaserParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_copy); if (rho) { - DepositCharge(pti, wp, rho, 0, 0, np_current, thread_num, lev, lev); + int* AMREX_RESTRICT ion_lev = nullptr; if (crho) { - DepositCharge(pti, wp, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); + DepositCharge(pti, wp, ion_lev, crho, 0, np_current, + np-np_current, thread_num, lev, lev-1); } } @@ -529,9 +530,12 @@ LaserParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_copy); if (rho) { - DepositCharge(pti, wp, rho, 1, 0, np_current, thread_num, lev, lev); + int* AMREX_RESTRICT ion_lev = nullptr; + DepositCharge(pti, wp, ion_lev, rho, 1, 0, + np_current, thread_num, lev, lev); if (crho) { - DepositCharge(pti, wp, crho, 1, np_current, np-np_current, thread_num, lev, lev-1); + DepositCharge(pti, wp, ion_lev, crho, 1, np_current, + np-np_current, thread_num, lev, lev-1); } } diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index 26d42370e..4253f5e76 100755 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -6,6 +6,10 @@ /* \brief Charge Deposition for thread thread_num * /param xp, yp, zp : Pointer to arrays of particle positions. * \param wp : Pointer to array of particle weights. + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. * \param rho_arr : Array4 of charge density, either full array or tile. * \param np_to_depose : Number of particles for which current is deposited. * \param dx : 3D cell size @@ -18,6 +22,7 @@ void doChargeDepositionShapeN(const amrex::Real * const xp, const amrex::Real * const yp, const amrex::Real * const zp, const amrex::Real * const wp, + const int * const ion_lev, const amrex::Array4& rho_arr, const long np_to_depose, const std::array& dx, @@ -25,6 +30,9 @@ void doChargeDepositionShapeN(const amrex::Real * const xp, const amrex::Dim3 lo, const amrex::Real q) { + // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer + // (do_ionization=1) + const bool do_ionization = ion_lev; const amrex::Real dxi = 1.0/dx[0]; const amrex::Real dzi = 1.0/dx[2]; #if (AMREX_SPACEDIM == 2) @@ -43,7 +51,10 @@ void doChargeDepositionShapeN(const amrex::Real * const xp, np_to_depose, [=] AMREX_GPU_DEVICE (long ip) { // --- Get particle quantities - const amrex::Real wq = q*wp[ip]*invvol; + amrex::Real wq = q*wp[ip]*invvol; + if (do_ionization){ + wq *= ion_lev[ip]; + } // --- Compute shape factors // x direction diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 4562b5c0a..3db413104 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1179,9 +1179,18 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_copy); if (rho) { - DepositCharge(pti, wp, rho, 0, 0, np_current, thread_num, lev, lev); + + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } + DepositCharge(pti, wp, ion_lev, rho, 0, 0, + np_current, thread_num, lev, lev); if (has_buffer){ - DepositCharge(pti, wp, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); + DepositCharge(pti, wp, ion_lev, crho, 0, np_current, + np-np_current, thread_num, lev, lev-1); } } @@ -1327,9 +1336,17 @@ PhysicalParticleContainer::Evolve (int lev, } if (rho) { - DepositCharge(pti, wp, rho, 1, 0, np_current, thread_num, lev, lev); + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } + DepositCharge(pti, wp, ion_lev, rho, 1, 0, + np_current, thread_num, lev, lev); if (has_buffer){ - DepositCharge(pti, wp, crho, 1, np_current, np-np_current, thread_num, lev, lev-1); + DepositCharge(pti, wp, ion_lev, crho, 1, np_current, + np-np_current, thread_num, lev, lev-1); } } diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 65eec3c2d..c252eee6b 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -171,6 +171,7 @@ public: virtual void DepositCharge(WarpXParIter& pti, RealVector& wp, + const int * const ion_lev, amrex::MultiFab* rho, int icomp, const long offset, diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 87c2d5e90..20f8c3e00 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -562,6 +562,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, void WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, + const int * const ion_lev, MultiFab* rho, int icomp, const long offset, const long np_to_depose, int thread_num, int lev, int depos_lev) @@ -623,14 +624,14 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, BL_PROFILE_VAR_START(blp_ppc_chd); if (WarpX::nox == 1){ - doChargeDepositionShapeN<1>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, - np_to_depose, dx, xyzmin, lo, q); + doChargeDepositionShapeN<1>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + rho_arr, np_to_depose, dx, xyzmin, lo, q); } else if (WarpX::nox == 2){ - doChargeDepositionShapeN<2>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, - np_to_depose, dx, xyzmin, lo, q); + doChargeDepositionShapeN<2>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + rho_arr, np_to_depose, dx, xyzmin, lo, q); } else if (WarpX::nox == 3){ - doChargeDepositionShapeN<3>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, - np_to_depose, dx, xyzmin, lo, q); + doChargeDepositionShapeN<3>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + rho_arr, np_to_depose, dx, xyzmin, lo, q); } BL_PROFILE_VAR_STOP(blp_ppc_chd); @@ -740,7 +741,15 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - DepositCharge(pti, wp, rho.get(), 0, 0, np, thread_num, lev, lev); + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } + + DepositCharge(pti, wp, ion_lev, rho.get(), 0, 0, np, + thread_num, lev, lev); } #ifdef _OPENMP } -- cgit v1.2.3 From d871e6f89e33aa5e9d5b308d19f84e86c5383d9a Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 21 Aug 2019 09:23:45 -0700 Subject: add comments for charge deposition --- Source/Particles/PhysicalParticleContainer.cpp | 3 ++- Source/Particles/WarpXParticleContainer.cpp | 22 ++++++++++++++++++++++ 2 files changed, 24 insertions(+), 1 deletion(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 3db413104..a82c97945 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1179,7 +1179,7 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_copy); if (rho) { - + // Deposit charge before particle push, in component 0 of MultiFab rho. int* AMREX_RESTRICT ion_lev; if (do_field_ionization){ ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); @@ -1336,6 +1336,7 @@ PhysicalParticleContainer::Evolve (int lev, } if (rho) { + // Deposit charge after particle push, in component 1 of MultiFab rho. int* AMREX_RESTRICT ion_lev; if (do_field_ionization){ ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 20f8c3e00..af56c58fb 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -416,6 +416,10 @@ WarpXParticleContainer::DepositCurrentFortran(WarpXParIter& pti, * \param pti : Particle iterator * \param wp : Array of particle weights * \param uxp uyp uzp : Array of particle + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. * \param jx jy jz : Full array of current density * \param offset : Index of first particle for which current is deposited * \param np_to_depose: Number of particles for which current is deposited. @@ -560,6 +564,24 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, #endif } +/* \brief Charge Deposition for thread thread_num + * \param pti : Particle iterator + * \param wp : Array of particle weights + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. + * \param rho : Full array of charge density + * \param icomp : Component of rho into which charge is deposited. + 0: old value (before particle push). + 1: new value (after particle push). + * \param offset : Index of first particle for which charge is deposited + * \param np_to_depose: Number of particles for which charge is deposited. + Particles [offset,offset+np_tp_depose] deposit charge + * \param thread_num : Thread number (if tiling) + * \param lev : Level of box that contains particles + * \param depos_lev : Level on which particles deposit (if buffers are used) + */ void WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, const int * const ion_lev, -- cgit v1.2.3 From eb4080b05342edd63c09091f393638c9053b7b19 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Wed, 21 Aug 2019 11:31:10 -0700 Subject: reimplement the boosted frame diagnostics without the runtime components --- Source/Diagnostics/ParticleIO.cpp | 11 ---- Source/Evolve/WarpXEvolveEM.cpp | 16 +++--- Source/Particles/MultiParticleContainer.cpp | 20 ------- Source/Particles/PhysicalParticleContainer.cpp | 65 ++++++---------------- .../Particles/RigidInjectedParticleContainer.cpp | 43 +++++--------- Source/Particles/WarpXParticleContainer.H | 2 + Source/Particles/WarpXParticleContainer.cpp | 27 --------- 7 files changed, 41 insertions(+), 143 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index f159e5302..4ac5cb6a1 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -102,17 +102,6 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const real_names.push_back("theta"); #endif - if (WarpX::do_boosted_frame_diagnostic && pc->DoBoostedFrameDiags()) - { - real_names.push_back("xold"); - real_names.push_back("yold"); - real_names.push_back("zold"); - - real_names.push_back("uxold"); - real_names.push_back("uyold"); - real_names.push_back("uzold"); - } - // Convert momentum to SI pc->ConvertUnits(ConvertDirection::WarpX_to_SI); // real_names contains a list of all particle attributes. diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 16b5905d1..ff55724f8 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -140,6 +140,14 @@ WarpX::EvolveEM (int numsteps) bool do_insitu = ((step+1) >= insitu_start) && (insitu_int > 0) && ((step+1) % insitu_int == 0); + if (do_boosted_frame_diagnostic) { + std::unique_ptr cell_centered_data = nullptr; + if (WarpX::do_boosted_frame_fields) { + cell_centered_data = GetCellCenteredData(); + } + myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); + } + bool move_j = is_synchronized || to_make_plot || do_insitu; // If is_synchronized we need to shift j too so that next step we can evolve E by dt/2. // We might need to move j because we are going to make a plotfile. @@ -173,14 +181,6 @@ WarpX::EvolveEM (int numsteps) t_new[i] = cur_time; } - if (do_boosted_frame_diagnostic) { - std::unique_ptr cell_centered_data = nullptr; - if (WarpX::do_boosted_frame_fields) { - cell_centered_data = GetCellCenteredData(); - } - myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); - } - // slice gen // if (to_make_plot || do_insitu || to_make_slice_plot) { diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 982e04e39..4980bc62a 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -42,26 +42,6 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) nspecies_boosted_frame_diags += 1; } } - - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - for (int i = 0; i < nspecies_boosted_frame_diags; ++i) - { - int is = map_species_boosted_frame_diags[i]; - allcontainers[is]->AddRealComp("xold"); - allcontainers[is]->AddRealComp("yold"); - allcontainers[is]->AddRealComp("zold"); - allcontainers[is]->AddRealComp("uxold"); - allcontainers[is]->AddRealComp("uyold"); - allcontainers[is]->AddRealComp("uzold"); - } - pc_tmp->AddRealComp("xold"); - pc_tmp->AddRealComp("yold"); - pc_tmp->AddRealComp("zold"); - pc_tmp->AddRealComp("uxold"); - pc_tmp->AddRealComp("uyold"); - pc_tmp->AddRealComp("uzold"); - } } void diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index d10390204..14807d0d9 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -196,18 +196,6 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, attribs[PIdx::uz] = u[2]; attribs[PIdx::w ] = weight; - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - // need to create old values - auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - particle_tile.push_back_real(particle_comps["xold"], x); - particle_tile.push_back_real(particle_comps["yold"], y); - particle_tile.push_back_real(particle_comps["zold"], z); - - particle_tile.push_back_real(particle_comps["uxold"], u[0]); - particle_tile.push_back_real(particle_comps["uyold"], u[1]); - particle_tile.push_back_real(particle_comps["uzold"], u[2]); - } // add particle AddOneParticle(0, 0, 0, x, y, z, attribs); } @@ -430,15 +418,6 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) for (int ia = 0; ia < PIdx::nattribs; ++ia) { pa[ia] = soa.GetRealData(ia).data() + old_size; } - GpuArray pb; - if (do_boosted) { - pb[0] = soa.GetRealData(particle_comps[ "xold"]).data() + old_size; - pb[1] = soa.GetRealData(particle_comps[ "yold"]).data() + old_size; - pb[2] = soa.GetRealData(particle_comps[ "zold"]).data() + old_size; - pb[3] = soa.GetRealData(particle_comps["uxold"]).data() + old_size; - pb[4] = soa.GetRealData(particle_comps["uyold"]).data() + old_size; - pb[5] = soa.GetRealData(particle_comps["uzold"]).data() + old_size; - } const GpuArray overlap_corner {AMREX_D_DECL(overlap_realbox.lo(0), @@ -589,15 +568,6 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) pa[PIdx::uy][ip] = u.y; pa[PIdx::uz][ip] = u.z; - if (do_boosted) { - pb[0][ip] = x; - pb[1][ip] = y; - pb[2][ip] = z; - pb[3][ip] = u.x; - pb[4][ip] = u.y; - pb[5][ip] = u.z; - } - #if (AMREX_SPACEDIM == 3) p.pos(0) = x; p.pos(1) = y; @@ -1633,21 +1603,22 @@ PhysicalParticleContainer::PushP (int lev, Real dt, void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, const Real* yp, const Real* zp) { - auto& attribs = pti.GetAttribs(); - Real* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); Real* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); - Real* AMREX_RESTRICT xpold = pti.GetAttribs(particle_comps["xold"]).dataPtr(); - Real* AMREX_RESTRICT ypold = pti.GetAttribs(particle_comps["yold"]).dataPtr(); - Real* AMREX_RESTRICT zpold = pti.GetAttribs(particle_comps["zold"]).dataPtr(); - Real* AMREX_RESTRICT uxpold = pti.GetAttribs(particle_comps["uxold"]).dataPtr(); - Real* AMREX_RESTRICT uypold = pti.GetAttribs(particle_comps["uyold"]).dataPtr(); - Real* AMREX_RESTRICT uzpold = pti.GetAttribs(particle_comps["uzold"]).dataPtr(); - - const long np = pti.numParticles(); + const auto np = pti.numParticles(); + const auto lev = pti.GetLevel(); + const auto index = pti.GetPairIndex(); + tmp_particle_data.resize(finestLevel()+1); + for (int i = 0; i < 6; ++i) tmp_particle_data[lev][index][i].resize(np); + Real* AMREX_RESTRICT xpold = tmp_particle_data[lev][index][0].dataPtr(); + Real* AMREX_RESTRICT ypold = tmp_particle_data[lev][index][1].dataPtr(); + Real* AMREX_RESTRICT zpold = tmp_particle_data[lev][index][2].dataPtr(); + Real* AMREX_RESTRICT uxpold = tmp_particle_data[lev][index][3].dataPtr(); + Real* AMREX_RESTRICT uypold = tmp_particle_data[lev][index][4].dataPtr(); + Real* AMREX_RESTRICT uzpold = tmp_particle_data[lev][index][5].dataPtr(); ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { @@ -1733,15 +1704,15 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real auto& uyp_new = attribs[PIdx::uy ]; auto& uzp_new = attribs[PIdx::uz ]; - auto& xp_old = pti.GetAttribs(particle_comps["xold"]); - auto& yp_old = pti.GetAttribs(particle_comps["yold"]); - auto& zp_old = pti.GetAttribs(particle_comps["zold"]); - auto& uxp_old = pti.GetAttribs(particle_comps["uxold"]); - auto& uyp_old = pti.GetAttribs(particle_comps["uyold"]); - auto& uzp_old = pti.GetAttribs(particle_comps["uzold"]); + auto& xp_old = tmp_particle_data[lev][index][0]; + auto& yp_old = tmp_particle_data[lev][index][1]; + auto& zp_old = tmp_particle_data[lev][index][2]; + auto& uxp_old = tmp_particle_data[lev][index][3]; + auto& uyp_old = tmp_particle_data[lev][index][4]; + auto& uzp_old = tmp_particle_data[lev][index][5]; const long np = pti.numParticles(); - + Real uzfrm = -WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c; Real inv_c2 = 1.0/PhysConst::c/PhysConst::c; diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 36cb9d224..0c94d1e33 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -239,15 +239,13 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, Real* const AMREX_RESTRICT Bzp = attribs[PIdx::Bz].dataPtr(); if (!done_injecting_lev) { - if (!(WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)) { - // If the old values are not already saved, create copies here. - xp_save = xp; - yp_save = yp; - zp_save = zp; - uxp_save = uxp; - uyp_save = uyp; - uzp_save = uzp; - } + // If the old values are not already saved, create copies here. + xp_save = xp; + yp_save = yp; + zp_save = zp; + uxp_save = uxp; + uyp_save = uyp; + uzp_save = uzp; // Scale the fields of particles about to cross the injection plane. // This only approximates what should be happening. The particles @@ -275,27 +273,12 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, if (!done_injecting_lev) { - Real* AMREX_RESTRICT x_save; - Real* AMREX_RESTRICT y_save; - Real* AMREX_RESTRICT z_save; - Real* AMREX_RESTRICT ux_save; - Real* AMREX_RESTRICT uy_save; - Real* AMREX_RESTRICT uz_save; - if (!(WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)) { - x_save = xp_save.dataPtr(); - y_save = yp_save.dataPtr(); - z_save = zp_save.dataPtr(); - ux_save = uxp_save.dataPtr(); - uy_save = uyp_save.dataPtr(); - uz_save = uzp_save.dataPtr(); - } else { - x_save = pti.GetAttribs(particle_comps["xold"]).dataPtr(); - y_save = pti.GetAttribs(particle_comps["yold"]).dataPtr(); - z_save = pti.GetAttribs(particle_comps["zold"]).dataPtr(); - ux_save = pti.GetAttribs(particle_comps["uxold"]).dataPtr(); - uy_save = pti.GetAttribs(particle_comps["uyold"]).dataPtr(); - uz_save = pti.GetAttribs(particle_comps["uzold"]).dataPtr(); - } + Real* AMREX_RESTRICT x_save = xp_save.dataPtr(); + Real* AMREX_RESTRICT y_save = yp_save.dataPtr(); + Real* AMREX_RESTRICT z_save = zp_save.dataPtr(); + Real* AMREX_RESTRICT ux_save = uxp_save.dataPtr(); + Real* AMREX_RESTRICT uy_save = uyp_save.dataPtr(); + Real* AMREX_RESTRICT uz_save = uzp_save.dataPtr(); // Undo the push for particles not injected yet. // The zp are advanced a fixed amount. diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index ac5b47ada..cbdf0cdbc 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -308,6 +308,8 @@ protected: // list of names of attributes to dump. amrex::Vector plot_vars; + amrex::Vector, std::array, 6> > > tmp_particle_data; + private: virtual void particlePostLocate(ParticleType& p, const amrex::ParticleLocData& pld, const int lev) override; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index befa5cfed..db0f683c7 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -85,17 +85,6 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) particle_comps["theta"] = PIdx::theta; #endif - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - particle_comps["xold"] = PIdx::nattribs; - particle_comps["yold"] = PIdx::nattribs+1; - particle_comps["zold"] = PIdx::nattribs+2; - particle_comps["uxold"] = PIdx::nattribs+3; - particle_comps["uyold"] = PIdx::nattribs+4; - particle_comps["uzold"] = PIdx::nattribs+5; - - } - // Initialize temporary local arrays for charge/current deposition int num_threads = 1; #ifdef _OPENMP @@ -238,14 +227,6 @@ WarpXParticleContainer::AddNParticles (int lev, p.pos(1) = z[i]; #endif - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - auto& ptile = DefineAndReturnParticleTile(0, 0, 0); - ptile.push_back_real(particle_comps["xold"], x[i]); - ptile.push_back_real(particle_comps["yold"], y[i]); - ptile.push_back_real(particle_comps["zold"], z[i]); - } - particle_tile.push_back(p); } @@ -256,14 +237,6 @@ WarpXParticleContainer::AddNParticles (int lev, particle_tile.push_back_real(PIdx::uy, vy + ibegin, vy + iend); particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend); - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - auto& ptile = DefineAndReturnParticleTile(0, 0, 0); - ptile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); - ptile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); - ptile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); - } - for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) { #ifdef WARPX_DIM_RZ -- cgit v1.2.3