From e80ca983aa27b023b0528c0ed2875a5745ea6599 Mon Sep 17 00:00:00 2001 From: David Grote Date: Fri, 7 Jun 2024 17:54:08 -0700 Subject: [PATCH] Cleanup diagnostic dimension macros (#4973) * Clean up BeamRelavent diagnostics, using get_particle_position * In Field diagnostics, simplify calculation of dV * Clean up of ParticleExtrema * In WarpXOpenPMD, removed duplicate particle coordinate transformation with RZ * Added more amrex prefixes in ParticleExtrema * Fix const * Small cleanup in WarpXOpenPMD * Clean up calculation of dV * Cleanup call to CellSize * Fix in FieldEnergy --- .../Diagnostics/ReducedDiags/BeamRelevant.cpp | 54 +--- .../Diagnostics/ReducedDiags/FieldEnergy.cpp | 14 +- .../ReducedDiags/FieldMomentum.cpp | 12 +- .../ReducedDiags/ParticleExtrema.cpp | 271 ++++++------------ Source/Diagnostics/WarpXOpenPMD.cpp | 43 +-- 5 files changed, 119 insertions(+), 275 deletions(-) diff --git a/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp b/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp index ae7d5230b4c..9b5a28a3516 100644 --- a/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp +++ b/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp @@ -175,15 +175,6 @@ void BeamRelevant::ComputeDiags (int step) // inverse of speed of light squared Real constexpr inv_c2 = 1.0_rt / (PhysConst::c * PhysConst::c); - // If 2D-XZ, p.pos(1) is z, rather than p.pos(2). -#if (defined WARPX_DIM_3D) - int const index_z = 2; -#elif (defined WARPX_DIM_XZ || defined WARPX_DIM_RZ) - int const index_z = 1; -#elif (defined WARPX_DIM_1D_Z) - int const index_z = 0; -#endif - // loop over species for (int i_s = 0; i_s < nSpecies; ++i_s) { @@ -212,26 +203,14 @@ void BeamRelevant::ComputeDiags (int step) const ParticleReal p_uy = p.rdata(PIdx::uy); const ParticleReal p_uz = p.rdata(PIdx::uz); const ParticleReal p_us = p_ux*p_ux + p_uy*p_uy + p_uz*p_uz; - const ParticleReal p_pos0 = p.pos(0); const ParticleReal p_w = p.rdata(PIdx::w); -#if defined(WARPX_DIM_3D) - const ParticleReal p_pos1 = p.pos(1); - const ParticleReal p_x_mean = p_pos0*p_w; - const ParticleReal p_y_mean = p_pos1*p_w; -#elif defined(WARPX_DIM_RZ) - const ParticleReal p_theta = p.rdata(PIdx::theta); - const ParticleReal p_x_mean = p_pos0*std::cos(p_theta)*p_w; - const ParticleReal p_y_mean = p_pos0*std::sin(p_theta)*p_w; -#elif defined(WARPX_DIM_XZ) - const ParticleReal p_x_mean = p_pos0*p_w; - const ParticleReal p_y_mean = 0; -#elif defined(WARPX_DIM_1D_Z) - amrex::ignore_unused(p_pos0); - const ParticleReal p_x_mean = 0; - const ParticleReal p_y_mean = 0; -#endif - const ParticleReal p_z_mean = p.pos(index_z)*p_w; + ParticleReal p_x, p_y, p_z; + get_particle_position(p, p_x, p_y, p_z); + + const ParticleReal p_x_mean = p_x*p_w; + const ParticleReal p_y_mean = p_y*p_w; + const ParticleReal p_z_mean = p_z*p_w; const ParticleReal p_ux_mean = p_ux*p_w; const ParticleReal p_uy_mean = p_uy*p_w; @@ -292,25 +271,8 @@ void BeamRelevant::ComputeDiags (int step) const ParticleReal p_gm = std::sqrt(1.0_rt+p_us*inv_c2); const ParticleReal p_w = p.rdata(PIdx::w); -#if (defined WARPX_DIM_1D_Z) - const ParticleReal p_x = 0.0; - const ParticleReal p_y = 0.0; -#elif (defined WARPX_DIM_RZ) - const ParticleReal p_pos0 = p.pos(0); - const ParticleReal p_theta = p.rdata(PIdx::theta); - const ParticleReal p_x = p_pos0*std::cos(p_theta); - const ParticleReal p_y = p_pos0*std::sin(p_theta); -#elif (defined WARPX_DIM_XZ) - const ParticleReal p_pos0 = p.pos(0); - const ParticleReal p_x = p_pos0; - const ParticleReal p_y = 0.0; -#else - const ParticleReal p_pos0 = p.pos(0); - const ParticleReal p_pos1 = p.pos(1); - const ParticleReal p_x = p_pos0; - const ParticleReal p_y = p_pos1; -#endif - const ParticleReal p_z = p.pos(index_z); + ParticleReal p_x, p_y, p_z; + get_particle_position(p, p_x, p_y, p_z); const ParticleReal p_x_ms = (p_x-x_mean)*(p_x-x_mean)*p_w; const ParticleReal p_y_ms = (p_y-y_mean)*(p_y-y_mean)*p_w; diff --git a/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp b/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp index f4b4e2a39a1..40ef1a088e6 100644 --- a/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp @@ -98,15 +98,9 @@ void FieldEnergy::ComputeDiags (int step) const MultiFab & By = warpx.getField(FieldType::Bfield_aux, lev,1); const MultiFab & Bz = warpx.getField(FieldType::Bfield_aux, lev,2); - // get cell size - Geometry const & geom = warpx.Geom(lev); -#if defined(WARPX_DIM_1D_Z) - auto dV = geom.CellSize(0); -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - auto dV = geom.CellSize(0) * geom.CellSize(1); -#elif defined(WARPX_DIM_3D) - auto dV = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2); -#endif + // get cell volume + const std::array &dx = WarpX::CellSize(lev); + const amrex::Real dV = dx[0]*dx[1]*dx[2]; #if defined(WARPX_DIM_RZ) amrex::Real const tmpEx = ComputeNorm2RZ(Ex, lev); @@ -119,6 +113,8 @@ void FieldEnergy::ComputeDiags (int step) amrex::Real const tmpBz = ComputeNorm2RZ(Bz, lev); amrex::Real const Bs = tmpBx + tmpBy + tmpBz; #else + Geometry const & geom = warpx.Geom(lev); + // compute E squared Real const tmpEx = Ex.norm2(0,geom.periodicity()); Real const tmpEy = Ey.norm2(0,geom.periodicity()); diff --git a/Source/Diagnostics/ReducedDiags/FieldMomentum.cpp b/Source/Diagnostics/ReducedDiags/FieldMomentum.cpp index f182acd5ba2..7eb16efecff 100644 --- a/Source/Diagnostics/ReducedDiags/FieldMomentum.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldMomentum.cpp @@ -183,15 +183,9 @@ void FieldMomentum::ComputeDiags (int step) amrex::Real ExB_z = amrex::get<2>(r); amrex::ParallelDescriptor::ReduceRealSum({ExB_x,ExB_y,ExB_z}); - // Get cell size - amrex::Geometry const & geom = warpx.Geom(lev); -#if defined(WARPX_DIM_1D_Z) - auto dV = geom.CellSize(0); -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - auto dV = geom.CellSize(0) * geom.CellSize(1); -#elif defined(WARPX_DIM_3D) - auto dV = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2); -#endif + // Get cell volume + const std::array &dx = WarpX::CellSize(lev); + const amrex::Real dV = dx[0]*dx[1]*dx[2]; // Save data (offset: 3 values for each refinement level) const int offset = lev*3; diff --git a/Source/Diagnostics/ReducedDiags/ParticleExtrema.cpp b/Source/Diagnostics/ReducedDiags/ParticleExtrema.cpp index 0cc5429be7a..9adfd3f238c 100644 --- a/Source/Diagnostics/ReducedDiags/ParticleExtrema.cpp +++ b/Source/Diagnostics/ReducedDiags/ParticleExtrema.cpp @@ -51,7 +51,7 @@ #include #include -using namespace amrex; +using namespace amrex::literals; using namespace warpx::fields; // constructor @@ -59,7 +59,7 @@ ParticleExtrema::ParticleExtrema (const std::string& rd_name) : ReducedDiags{rd_name} { // read species name - const ParmParse pp_rd_name(rd_name); + const amrex::ParmParse pp_rd_name(rd_name); pp_rd_name.get("species",m_species_name); // get WarpX class object @@ -122,7 +122,7 @@ ParticleExtrema::ParticleExtrema (const std::string& rd_name) m_data.resize(all_diag_names.size()); - if (ParallelDescriptor::IOProcessor()) + if (amrex::ParallelDescriptor::IOProcessor()) { if ( m_write_header ) { @@ -165,16 +165,7 @@ void ParticleExtrema::ComputeDiags (int step) const auto species_names = mypc.GetSpeciesNames(); // inverse of speed of light squared - Real constexpr inv_c2 = 1.0_rt / (PhysConst::c * PhysConst::c); - - // If 2D-XZ, p.pos(1) is z, rather than p.pos(2). -#if (defined WARPX_DIM_3D) - int const index_z = 2; -#elif (defined WARPX_DIM_XZ || defined WARPX_DIM_RZ) - int const index_z = 1; -#elif (defined WARPX_DIM_1D_Z) - int const index_z = 0; -#endif + amrex::Real constexpr inv_c2 = 1.0_rt / (PhysConst::c * PhysConst::c); // loop over species for (int i_s = 0; i_s < nSpecies; ++i_s) @@ -193,172 +184,72 @@ void ParticleExtrema::ComputeDiags (int step) } using PType = typename WarpXParticleContainer::SuperParticleType; - - // xmin -#if (defined WARPX_DIM_RZ) - Real xmin = ReduceMin( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { return p.pos(0)*std::cos(p.rdata(PIdx::theta)); }); -#elif (defined WARPX_DIM_1D_Z) - Real xmin = 0.0_rt; -#else - Real xmin = ReduceMin( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { return p.pos(0); }); -#endif - - // xmax -#if (defined WARPX_DIM_RZ) - Real xmax = ReduceMax( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { return p.pos(0)*std::cos(p.rdata(PIdx::theta)); }); -#elif (defined WARPX_DIM_1D_Z) - Real xmax = 0.0_rt; -#else - Real xmax = ReduceMax( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { return p.pos(0); }); -#endif - - // ymin -#if (defined WARPX_DIM_RZ) - Real ymin = ReduceMin( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { return p.pos(0)*std::sin(p.rdata(PIdx::theta)); }); -#elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z) - Real ymin = 0.0_rt; -#else - Real ymin = ReduceMin( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { return p.pos(1); }); -#endif - - // ymax -#if (defined WARPX_DIM_RZ) - Real ymax = ReduceMax( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { return p.pos(0)*std::sin(p.rdata(PIdx::theta)); }); -#elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z) - Real ymax = 0.0_rt; -#else - Real ymax = ReduceMax( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { return p.pos(1); }); -#endif - - // zmin - Real zmin = ReduceMin( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { return p.pos(index_z); }); - - // zmax - Real zmax = ReduceMax( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { return p.pos(index_z); }); - - // uxmin - Real uxmin = ReduceMin( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { return p.rdata(PIdx::ux); }); - - // uxmax - Real uxmax = ReduceMax( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { return p.rdata(PIdx::ux); }); - - // uymin - Real uymin = ReduceMin( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { return p.rdata(PIdx::uy); }); - - // uymax - Real uymax = ReduceMax( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { return p.rdata(PIdx::uy); }); - - // uzmin - Real uzmin = ReduceMin( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { return p.rdata(PIdx::uz); }); - - // uzmax - Real uzmax = ReduceMax( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { return p.rdata(PIdx::uz); }); - - // gmin - Real gmin = 0.0_rt; - if ( is_photon ) { - gmin = ReduceMin( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) + using OpMin = amrex::ReduceOpMin; + using OpMax = amrex::ReduceOpMax; + + amrex::ReduceOps reduce_ops; + auto posminmax = amrex::ParticleReduce>( + myspc, + [=] AMREX_GPU_DEVICE(const PType& p) noexcept -> amrex::GpuTuple { - const Real ux = p.rdata(PIdx::ux); - const Real uy = p.rdata(PIdx::uy); - const Real uz = p.rdata(PIdx::uz); - const Real us = ux*ux + uy*uy + uz*uz; - return std::sqrt(us*inv_c2); - }); - } else { - gmin = ReduceMin( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) + amrex::ParticleReal x, y, z; + get_particle_position(p, x, y, z); + amrex::Real const w = p.rdata(PIdx::w); + return {w, x, y, z, w, x, y, z}; + }, + reduce_ops); + + amrex::Real wmin = amrex::get<0>(posminmax); + amrex::Real xmin = amrex::get<1>(posminmax); + amrex::Real ymin = amrex::get<2>(posminmax); + amrex::Real zmin = amrex::get<3>(posminmax); + amrex::Real wmax = amrex::get<4>(posminmax); + amrex::Real xmax = amrex::get<5>(posminmax); + amrex::Real ymax = amrex::get<6>(posminmax); + amrex::Real zmax = amrex::get<7>(posminmax); + + amrex::Real const gfactor = (is_photon ? 0._rt : 1._rt); + auto uminmax = amrex::ParticleReduce>( + myspc, + [=] AMREX_GPU_DEVICE(const PType& p) noexcept -> amrex::GpuTuple { - const Real ux = p.rdata(PIdx::ux); - const Real uy = p.rdata(PIdx::uy); - const Real uz = p.rdata(PIdx::uz); - const Real us = ux*ux + uy*uy + uz*uz; - return std::sqrt(1.0_rt + us*inv_c2); - }); - } - - // gmax - Real gmax = 0.0_rt; - if ( is_photon ) { - gmax = ReduceMax( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { - const Real ux = p.rdata(PIdx::ux); - const Real uy = p.rdata(PIdx::uy); - const Real uz = p.rdata(PIdx::uz); - const Real us = ux*ux + uy*uy + uz*uz; - return std::sqrt(us*inv_c2); - }); - } else { - gmax = ReduceMax( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { - const Real ux = p.rdata(PIdx::ux); - const Real uy = p.rdata(PIdx::uy); - const Real uz = p.rdata(PIdx::uz); - const Real us = ux*ux + uy*uy + uz*uz; - return std::sqrt(1.0_rt + us*inv_c2); - }); - } - - // wmin - Real wmin = ReduceMin( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { return p.rdata(PIdx::w); }); - - // wmax - Real wmax = ReduceMax( myspc, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) - { return p.rdata(PIdx::w); }); - - ParallelDescriptor::ReduceRealMin({xmin,ymin,zmin,uxmin,uymin,uzmin,gmin,wmin}); - ParallelDescriptor::ReduceRealMax({xmax,ymax,zmax,uxmax,uymax,uzmax,gmax,wmax}); + amrex::Real const ux = p.rdata(PIdx::ux); + amrex::Real const uy = p.rdata(PIdx::uy); + amrex::Real const uz = p.rdata(PIdx::uz); + amrex::Real const g = std::sqrt(gfactor + (ux*ux + uy*uy + uz*uz)*inv_c2); + return {g, ux, uy, uz, g, ux, uy, uz}; + }, + reduce_ops); + + amrex::Real gmin = amrex::get<0>(uminmax); + amrex::Real uxmin = amrex::get<1>(uminmax); + amrex::Real uymin = amrex::get<2>(uminmax); + amrex::Real uzmin = amrex::get<3>(uminmax); + amrex::Real gmax = amrex::get<4>(uminmax); + amrex::Real uxmax = amrex::get<5>(uminmax); + amrex::Real uymax = amrex::get<6>(uminmax); + amrex::Real uzmax = amrex::get<7>(uminmax); + + amrex::ParallelDescriptor::ReduceRealMin({xmin,ymin,zmin,uxmin,uymin,uzmin,gmin,wmin}); + amrex::ParallelDescriptor::ReduceRealMax({xmax,ymax,zmax,uxmax,uymax,uzmax,gmax,wmax}); #if (defined WARPX_QED) // get number of level (int) const auto level_number = WarpX::GetInstance().finestLevel(); // compute chimin and chimax - Real chimin_f = 0.0_rt; - Real chimax_f = 0.0_rt; + amrex::Real chimin_f = 0.0_rt; + amrex::Real chimax_f = 0.0_rt; if (myspc.DoQED()) { // declare chi arrays - std::vector chimin, chimax; + std::vector chimin, chimax; chimin.resize(level_number+1,0.0_rt); chimax.resize(level_number+1,0.0_rt); @@ -374,17 +265,17 @@ void ParticleExtrema::ComputeDiags (int step) { // define variables in preparation for field gathering const std::array& dx = WarpX::CellSize(std::max(lev, 0)); - const GpuArray dx_arr = {dx[0], dx[1], dx[2]}; - const MultiFab & Ex = warpx.getField(FieldType::Efield_aux, lev,0); - const MultiFab & Ey = warpx.getField(FieldType::Efield_aux, lev,1); - const MultiFab & Ez = warpx.getField(FieldType::Efield_aux, lev,2); - const MultiFab & Bx = warpx.getField(FieldType::Bfield_aux, lev,0); - const MultiFab & By = warpx.getField(FieldType::Bfield_aux, lev,1); - const MultiFab & Bz = warpx.getField(FieldType::Bfield_aux, lev,2); + const amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; + const amrex::MultiFab & Ex = warpx.getField(FieldType::Efield_aux, lev,0); + const amrex::MultiFab & Ey = warpx.getField(FieldType::Efield_aux, lev,1); + const amrex::MultiFab & Ez = warpx.getField(FieldType::Efield_aux, lev,2); + const amrex::MultiFab & Bx = warpx.getField(FieldType::Bfield_aux, lev,0); + const amrex::MultiFab & By = warpx.getField(FieldType::Bfield_aux, lev,1); + const amrex::MultiFab & Bz = warpx.getField(FieldType::Bfield_aux, lev,2); // declare reduce_op - ReduceOps reduce_op; - ReduceData reduce_data(reduce_op); + amrex::ReduceOps reduce_op; + amrex::ReduceData reduce_data(reduce_op); using ReduceTuple = typename decltype(reduce_data)::Type; // Loop over boxes @@ -408,28 +299,28 @@ void ParticleExtrema::ComputeDiags (int step) // define variables in preparation for field gathering amrex::Box box = pti.tilebox(); box.grow(ngEB); - const Dim3 lo = amrex::lbound(box); + const amrex::Dim3 lo = amrex::lbound(box); const std::array& xyzmin = WarpX::LowerCorner(box, lev, 0._rt); - const GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; + const amrex::GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; const auto& ex_arr = Ex[pti].array(); const auto& ey_arr = Ey[pti].array(); const auto& ez_arr = Ez[pti].array(); const auto& bx_arr = Bx[pti].array(); const auto& by_arr = By[pti].array(); const auto& bz_arr = Bz[pti].array(); - const IndexType ex_type = Ex[pti].box().ixType(); - const IndexType ey_type = Ey[pti].box().ixType(); - const IndexType ez_type = Ez[pti].box().ixType(); - const IndexType bx_type = Bx[pti].box().ixType(); - const IndexType by_type = By[pti].box().ixType(); - const IndexType bz_type = Bz[pti].box().ixType(); + const amrex::IndexType ex_type = Ex[pti].box().ixType(); + const amrex::IndexType ey_type = Ey[pti].box().ixType(); + const amrex::IndexType ez_type = Ez[pti].box().ixType(); + const amrex::IndexType bx_type = Bx[pti].box().ixType(); + const amrex::IndexType by_type = By[pti].box().ixType(); + const amrex::IndexType bz_type = Bz[pti].box().ixType(); // evaluate reduce_op reduce_op.eval(pti.numParticles(), reduce_data, [=] AMREX_GPU_DEVICE (int i) -> ReduceTuple { // get external fields - ParticleReal xp, yp, zp; + amrex::ParticleReal xp, yp, zp; GetPosition(i, xp, yp, zp); amrex::ParticleReal ex = Ex_external_particle; amrex::ParticleReal ey = Ey_external_particle; @@ -449,7 +340,7 @@ void ParticleExtrema::ComputeDiags (int step) dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, nox, galerkin_interpolation); // compute chi - Real chi = 0.0_rt; + amrex::Real chi = 0.0_rt; if ( is_photon ) { chi = QedUtils::chi_photon(ux[i]*m, uy[i]*m, uz[i]*m, ex, ey, ez, bx, by, bz); @@ -461,13 +352,13 @@ void ParticleExtrema::ComputeDiags (int step) }); } auto val = reduce_data.value(); - chimin[lev] = get<0>(val); - chimax[lev] = get<1>(val); + chimin[lev] = amrex::get<0>(val); + chimax[lev] = amrex::get<1>(val); } chimin_f = *std::min_element(chimin.begin(), chimin.end()); chimax_f = *std::max_element(chimax.begin(), chimax.end()); - ParallelDescriptor::ReduceRealMin(chimin_f, ParallelDescriptor::IOProcessorNumber()); - ParallelDescriptor::ReduceRealMax(chimax_f, ParallelDescriptor::IOProcessorNumber()); + amrex::ParallelDescriptor::ReduceRealMin(chimin_f, amrex::ParallelDescriptor::IOProcessorNumber()); + amrex::ParallelDescriptor::ReduceRealMax(chimax_f, amrex::ParallelDescriptor::IOProcessorNumber()); } #endif diff --git a/Source/Diagnostics/WarpXOpenPMD.cpp b/Source/Diagnostics/WarpXOpenPMD.cpp index 45d5dc1cdee..36827bd316a 100644 --- a/Source/Diagnostics/WarpXOpenPMD.cpp +++ b/Source/Diagnostics/WarpXOpenPMD.cpp @@ -911,33 +911,34 @@ WarpXOpenPMDPlot::SaveRealProperty (ParticleIter& pti, { auto const real_counter = std::min(write_real_comp.size(), real_comp_names.size()); +#if defined(WARPX_DIM_RZ) // reconstruct Cartesian positions for RZ simulations // r,z,theta -> x,y,z -#if defined(WARPX_DIM_RZ) - auto const * const r = soa.GetRealData(PIdx::x).data(); - auto const * const theta = soa.GetRealData(PIdx::theta).data(); + // If each comp is being written, create a temporary array, otherwise create an empty array. + std::shared_ptr const x( + new amrex::ParticleReal[(write_real_comp[0] ? numParticleOnTile : 0)], + [](amrex::ParticleReal const *p) { delete[] p; } + ); + std::shared_ptr const y( + new amrex::ParticleReal[(write_real_comp[1] ? numParticleOnTile : 0)], + [](amrex::ParticleReal const *p) { delete[] p; } + ); + const auto& tile = pti.GetParticleTile(); + const auto& ptd = tile.getConstParticleTileData(); + + for (int i = 0; i < numParticleOnTile; ++i) { + const auto& p = ptd.getSuperParticle(i); + amrex::ParticleReal xp, yp, zp; + get_particle_position(p, xp, yp, zp); + if (write_real_comp[0]) { x.get()[i] = xp; } + if (write_real_comp[1]) { y.get()[i] = yp; } + } if (write_real_comp[0]) { - std::shared_ptr const x( - new amrex::ParticleReal[numParticleOnTile], - [](amrex::ParticleReal const *p) { delete[] p; } - ); - for (int i = 0; i < numParticleOnTile; ++i) { - x.get()[i] = r[i] * std::cos(theta[i]); - } - getComponentRecord(real_comp_names[0]).storeChunk( - x, {offset}, {numParticleOnTile64}); + getComponentRecord(real_comp_names[0]).storeChunk(x, {offset}, {numParticleOnTile64}); } if (write_real_comp[1]) { - std::shared_ptr const y( - new amrex::ParticleReal[numParticleOnTile], - [](amrex::ParticleReal const *p) { delete[] p; } - ); - for (int i = 0; i < numParticleOnTile; ++i) { - y.get()[i] = r[i] * std::sin(theta[i]); - } - getComponentRecord(real_comp_names[1]).storeChunk( - y, {offset}, {numParticleOnTile64}); + getComponentRecord(real_comp_names[1]).storeChunk(y, {offset}, {numParticleOnTile64}); } #endif