diff --git a/cmake/dependencies/ABLASTR.cmake b/cmake/dependencies/ABLASTR.cmake index f2df7a9c3..ff22e6441 100644 --- a/cmake/dependencies/ABLASTR.cmake +++ b/cmake/dependencies/ABLASTR.cmake @@ -178,7 +178,7 @@ set(ImpactX_openpmd_src "" set(ImpactX_ablastr_repo "https://github.com/ECP-WarpX/WarpX.git" CACHE STRING "Repository URI to pull and build ABLASTR from if(ImpactX_ablastr_internal)") -set(ImpactX_ablastr_branch "24.01" +set(ImpactX_ablastr_branch "94ae11900131846e9f8b3704194673f7f02d8959" CACHE STRING "Repository branch for ImpactX_ablastr_repo if(ImpactX_ablastr_internal)") diff --git a/examples/epac2004_benchmarks/input_fodo_rf_SC.in b/examples/epac2004_benchmarks/input_fodo_rf_SC.in index ade26631c..8008089e1 100644 --- a/examples/epac2004_benchmarks/input_fodo_rf_SC.in +++ b/examples/epac2004_benchmarks/input_fodo_rf_SC.in @@ -125,4 +125,4 @@ geometry.prob_relative = 4.0 ############################################################################### # Diagnostics ############################################################################### -diag.slice_step_diagnostics = true +diag.slice_step_diagnostics = false diff --git a/examples/fodo/run_fodo_programmable.py b/examples/fodo/run_fodo_programmable.py index 1af61bbdb..2d42b51c9 100755 --- a/examples/fodo/run_fodo_programmable.py +++ b/examples/fodo/run_fodo_programmable.py @@ -78,16 +78,16 @@ def my_drift(pge, pti, refpart): else: array = np.array - # access AoS data such as positions and cpu/id - aos = pti.aos() - aos_arr = array(aos, copy=False) - # access SoA data such as momentum + # access particle attributes soa = pti.soa() real_arrays = soa.GetRealData() - px = array(real_arrays[0], copy=False) - py = array(real_arrays[1], copy=False) - pt = array(real_arrays[2], copy=False) + x = array(real_arrays[0], copy=False) + y = array(real_arrays[1], copy=False) + t = array(real_arrays[2], copy=False) + px = array(real_arrays[3], copy=False) + py = array(real_arrays[4], copy=False) + pt = array(real_arrays[5], copy=False) # length of the current slice slice_ds = pge.ds / pge.nslice @@ -97,9 +97,9 @@ def my_drift(pge, pti, refpart): betgam2 = pt_ref**2 - 1.0 # advance position and momentum (drift) - aos_arr[:]["x"] += slice_ds * px[:] - aos_arr[:]["y"] += slice_ds * py[:] - aos_arr[:]["z"] += (slice_ds / betgam2) * pt[:] + x[:] += slice_ds * px[:] + y[:] += slice_ds * py[:] + t[:] += (slice_ds / betgam2) * pt[:] def my_ref_drift(pge, refpart): diff --git a/examples/pytorch_surrogate_model/run_ml_surrogate.py b/examples/pytorch_surrogate_model/run_ml_surrogate.py index 1de19c005..df09fc28c 100644 --- a/examples/pytorch_surrogate_model/run_ml_surrogate.py +++ b/examples/pytorch_surrogate_model/run_ml_surrogate.py @@ -11,6 +11,14 @@ from urllib import request import numpy as np + +try: + import cupy as cp + + cupy_available = True +except ImportError: + cupy_available = False + from surrogate_model_definitions import surrogate_model try: @@ -20,10 +28,11 @@ sys.exit(0) from impactx import ( + Config, + CoordSystem, ImpactX, ImpactXParIter, RefPart, - TransformationDirection, coordinate_transformation, distribution, elements, @@ -79,7 +88,18 @@ def __init__(self, stage_i, surrogate_model, surrogate_length, stage_start): self.ds = surrogate_length def surrogate_push(self, pc, step): - array = np.array + # CPU/GPU logic + if Config.have_gpu: + if cupy_available: + array = cp.array + else: + print("Warning: GPU found but cupy not available! Try managed...") + array = np.array + if Config.gpu_backend == "SYCL": + print("Warning: SYCL GPU backend not yet implemented for Python") + + else: + array = np.array ref_part = pc.ref_particle() ref_z_i = ref_part.z @@ -100,26 +120,28 @@ def surrogate_push(self, pc, step): ref_part_final = torch.tensor([0, 0, ref_z_f, 0, 0, ref_uz_f]) # transform - coordinate_transformation(pc, TransformationDirection.to_fixed_t) + coordinate_transformation(pc, direction=CoordSystem.t) for lvl in range(pc.finest_level + 1): for pti in ImpactXParIter(pc, level=lvl): - aos = pti.aos() - aos_arr = array(aos, copy=False) - soa = pti.soa() real_arrays = soa.GetRealData() - px = array(real_arrays[0], copy=False) - py = array(real_arrays[1], copy=False) - pt = array(real_arrays[2], copy=False) - data_arr = ( - torch.tensor( - np.vstack( - [aos_arr["x"], aos_arr["y"], aos_arr["z"], real_arrays[:3]] - ) - ) - .float() - .T + x = array(real_arrays[0], copy=False) + y = array(real_arrays[1], copy=False) + t = array(real_arrays[2], copy=False) + px = array(real_arrays[3], copy=False) + py = array(real_arrays[4], copy=False) + pt = array(real_arrays[5], copy=False) + data_arr = torch.stack( + ( + x, + y, + t, + px, + py, + py, + ), + dim=0, ) data_arr[:, 0] += ref_part.x @@ -147,9 +169,9 @@ def surrogate_push(self, pc, step): data_arr_post_model[:, 3 + ii] -= ref_part_final[3 + ii] data_arr_post_model[:, 3 + ii] /= ref_beta_gamma_final - aos_arr["x"] = data_arr_post_model[:, 0] - aos_arr["y"] = data_arr_post_model[:, 1] - aos_arr["z"] = data_arr_post_model[:, 2] + x[:] = data_arr_post_model[:, 0] + y[:] = data_arr_post_model[:, 1] + t[:] = data_arr_post_model[:, 2] px[:] = data_arr_post_model[:, 3] py[:] = data_arr_post_model[:, 4] pt[:] = data_arr_post_model[:, 5] @@ -174,7 +196,7 @@ def surrogate_push(self, pc, step): # ref_part.s += pge1.ds # ref_part.t += pge1.ds / ref_beta - coordinate_transformation(pc, TransformationDirection.to_fixed_s) + coordinate_transformation(pc, direction=CoordSystem.s) ## Done! diff --git a/src/particles/CollectLost.cpp b/src/particles/CollectLost.cpp index a8f8b8ede..5d4d27ba6 100644 --- a/src/particles/CollectLost.cpp +++ b/src/particles/CollectLost.cpp @@ -27,9 +27,9 @@ namespace impactx using DstData = ImpactXParticleContainer::ParticleTileType::ParticleTileDataType; AMREX_GPU_HOST_DEVICE - void operator() (DstData const &dst, SrcData const &src, int src_ip, int dst_ip) const noexcept { - dst.m_aos[dst_ip] = src.m_aos[src_ip]; - + void operator() (DstData const &dst, SrcData const &src, int src_ip, int dst_ip) const noexcept + { + dst.m_idcpu[dst_ip] = src.m_idcpu[src_ip]; for (int j = 0; j < SrcData::NAR; ++j) dst.m_rdata[j][dst_ip] = src.m_rdata[j][src_ip]; for (int j = 0; j < src.m_num_runtime_real; ++j) @@ -141,8 +141,7 @@ namespace impactx // move down int const new_index = ip - n_removed; - ptile_src_data.m_aos[new_index] = ptile_src_data.m_aos[ip]; - + ptile_src_data.m_idcpu[new_index] = ptile_src_data.m_idcpu[ip]; for (int j = 0; j < SrcData::NAR; ++j) ptile_src_data.m_rdata[j][new_index] = ptile_src_data.m_rdata[j][ip]; for (int j = 0; j < ptile_src_data.m_num_runtime_real; ++j) diff --git a/src/particles/ImpactXParticleContainer.H b/src/particles/ImpactXParticleContainer.H index fb0deccc0..d5a68e805 100644 --- a/src/particles/ImpactXParticleContainer.H +++ b/src/particles/ImpactXParticleContainer.H @@ -18,6 +18,7 @@ #include #include #include +#include #include #include @@ -35,43 +36,16 @@ namespace impactx t ///< fixed t as the independent variable }; - /** AMReX pre-defined Real attributes - * - * These are the AMReX pre-defined struct indexes for the Real attributes - * stored in an AoS in ImpactXParticleContainer. We document this here, - * because we change the meaning of these "positions" depending on the - * coordinate system we are currently in. - */ - struct RealAoS - { - enum - { - x, ///< position in x [m] (at fixed s OR fixed t) - y, ///< position in y [m] (at fixed s OR fixed t) - t, ///< c * time-of-flight [m] (at fixed s) - nattribs ///< the number of attributes above (always last) - }; - - // at fixed t, the third component represents the position z - enum { - z = t ///< position in z [m] (at fixed t) - }; - - //! named labels for fixed s - static constexpr auto names_s = { "position_x", "position_y", "position_t" }; - //! named labels for fixed t - static constexpr auto names_t = { "position_x", "position_y", "position_z" }; - static_assert(names_s.size() == nattribs); - static_assert(names_t.size() == nattribs); - }; - - /** This struct indexes the additional Real attributes + /** This struct indexes the Real attributes * stored in an SoA in ImpactXParticleContainer */ struct RealSoA { enum { + x, ///< position in x [m] (at fixed s or t) + y, ///< position in y [m] (at fixed s or t) + t, ///< time-of-flight ct [m] (at fixed s) px, ///< momentum in x, scaled by the magnitude of the reference momentum [unitless] (at fixed s or t) py, ///< momentum in y, scaled by the magnitude of the reference momentum [unitless] (at fixed s or t) pt, ///< energy deviation, scaled by speed of light * the magnitude of the reference momentum [unitless] (at fixed s) @@ -80,27 +54,28 @@ namespace impactx nattribs ///< the number of attributes above (always last) }; - // at fixed t, the third component represents the momentum in z + // at fixed t, the third component represents the position z, the 6th component represents the momentum in z enum { + z = t, ///< position in z [m] (at fixed t) pz = pt ///< momentum in z, scaled by the magnitude of the reference momentum [unitless] (at fixed t) }; //! named labels for fixed s - static constexpr auto names_s = { "momentum_x", "momentum_y", "momentum_t", "qm", "weighting" }; + static constexpr auto names_s = { "position_x", "position_y", "position_t", "momentum_x", "momentum_y", "momentum_t", "qm", "weighting" }; //! named labels for fixed t - static constexpr auto names_t = { "momentum_x", "momentum_y", "momentum_z", "qm", "weighting" }; + static constexpr auto names_t = { "position_x", "position_y", "position_z", "momentum_x", "momentum_y", "momentum_z", "qm", "weighting" }; static_assert(names_s.size() == nattribs); static_assert(names_t.size() == nattribs); }; - /** This struct indexes the additional Integer attributes + /** This struct indexes the Integer attributes * stored in an SoA in ImpactXParticleContainer */ struct IntSoA { enum { - nattribs ///< the number of particles above (always last) + nattribs ///< the number of attributes above (always last) }; }; @@ -109,15 +84,15 @@ namespace impactx * We subclass here to change the default threading strategy, which is * `static` in AMReX, to `dynamic` in ImpactX. */ - class ParIter - : public amrex::ParIter<0, 0, RealSoA::nattribs, IntSoA::nattribs> + class ParIterSoA + : public amrex::ParIterSoA { public: - using amrex::ParIter<0, 0, RealSoA::nattribs, IntSoA::nattribs>::ParIter; + using amrex::ParIterSoA::ParIterSoA; - ParIter (ContainerType& pc, int level); + ParIterSoA (ContainerType& pc, int level); - ParIter (ContainerType& pc, int level, amrex::MFItInfo& info); + ParIterSoA (ContainerType& pc, int level, amrex::MFItInfo& info); }; /** Const AMReX iterator for particle boxes - data is read only. @@ -125,15 +100,15 @@ namespace impactx * We subclass here to change the default threading strategy, which is * `static` in AMReX, to `dynamic` in ImpactX. */ - class ParConstIter - : public amrex::ParConstIter<0, 0, RealSoA::nattribs, IntSoA::nattribs> + class ParConstIterSoA + : public amrex::ParConstIterSoA { public: - using amrex::ParConstIter<0, 0, RealSoA::nattribs, IntSoA::nattribs>::ParConstIter; + using amrex::ParConstIterSoA::ParConstIterSoA; - ParConstIter (ContainerType& pc, int level); + ParConstIterSoA (ContainerType& pc, int level); - ParConstIter (ContainerType& pc, int level, amrex::MFItInfo& info); + ParConstIterSoA (ContainerType& pc, int level, amrex::MFItInfo& info); }; /** Beam Particles in ImpactX @@ -141,14 +116,14 @@ namespace impactx * This class stores particles, distributed over MPI ranks. */ class ImpactXParticleContainer - : public amrex::ParticleContainer<0, 0, RealSoA::nattribs, IntSoA::nattribs> + : public amrex::ParticleContainerPureSoA { public: //! amrex iterator for particle boxes - using iterator = impactx::ParIter; + using iterator = impactx::ParIterSoA; //! amrex constant iterator for particle boxes (read-only) - using const_iterator = impactx::ParConstIter; + using const_iterator = impactx::ParConstIterSoA; //! Construct a new particle container ImpactXParticleContainer (initialization::AmrCoreData* amr_core); @@ -276,10 +251,6 @@ namespace impactx DepositCharge (std::unordered_map & rho, amrex::Vector const & ref_ratio); - /** Get the name of each Real AoS component */ - std::vector - RealAoS_names () const; - /** Get the name of each Real SoA component */ std::vector RealSoA_names () const; @@ -311,10 +282,6 @@ namespace impactx }; // ImpactXParticleContainer - /** Get the name of each Real AoS component */ - std::vector - get_RealAoS_names (); - /** Get the name of each Real SoA component * * @param num_real_comps number of compile-time + runtime arrays diff --git a/src/particles/ImpactXParticleContainer.cpp b/src/particles/ImpactXParticleContainer.cpp index a28e33525..94f33fcdc 100644 --- a/src/particles/ImpactXParticleContainer.cpp +++ b/src/particles/ImpactXParticleContainer.cpp @@ -19,7 +19,7 @@ #include #include #include -#include +#include #include #include @@ -38,24 +38,24 @@ namespace namespace impactx { - ParIter::ParIter (ContainerType& pc, int level) - : amrex::ParIter<0, 0, RealSoA::nattribs, IntSoA::nattribs>(pc, level, + ParIterSoA::ParIterSoA (ContainerType& pc, int level) + : amrex::ParIterSoA(pc, level, amrex::MFItInfo().SetDynamic(do_omp_dynamic())) {} - ParIter::ParIter (ContainerType& pc, int level, amrex::MFItInfo& info) - : amrex::ParIter<0, 0, RealSoA::nattribs, IntSoA::nattribs>(pc, level, + ParIterSoA::ParIterSoA (ContainerType& pc, int level, amrex::MFItInfo& info) + : amrex::ParIterSoA(pc, level, info.SetDynamic(do_omp_dynamic())) {} - ParConstIter::ParConstIter (ContainerType& pc, int level) - : amrex::ParConstIter<0, 0, RealSoA::nattribs, IntSoA::nattribs>(pc, level, + ParConstIterSoA::ParConstIterSoA (ContainerType& pc, int level) + : amrex::ParConstIterSoA(pc, level, amrex::MFItInfo().SetDynamic(do_omp_dynamic())) {} - ParConstIter::ParConstIter (ContainerType& pc, int level, amrex::MFItInfo& info) - : amrex::ParConstIter<0, 0, RealSoA::nattribs, IntSoA::nattribs>(pc, level, + ParConstIterSoA::ParConstIterSoA (ContainerType& pc, int level, amrex::MFItInfo& info) + : amrex::ParConstIterSoA(pc, level, info.SetDynamic(do_omp_dynamic())) {} ImpactXParticleContainer::ImpactXParticleContainer (initialization::AmrCoreData* amr_core) - : amrex::ParticleContainer<0, 0, RealSoA::nattribs, IntSoA::nattribs>(amr_core->GetParGDB()) + : amrex::ParticleContainerPureSoA(amr_core->GetParGDB()) { SetParticleSize(); } @@ -157,14 +157,18 @@ namespace impactx const int cpuid = amrex::ParallelDescriptor::MyProc(); - auto * AMREX_RESTRICT pstructs = particle_tile.GetArrayOfStructs()().dataPtr(); auto & soa = particle_tile.GetStructOfArrays().GetRealData(); + amrex::ParticleReal * const AMREX_RESTRICT x_arr = soa[RealSoA::x].dataPtr(); + amrex::ParticleReal * const AMREX_RESTRICT y_arr = soa[RealSoA::y].dataPtr(); + amrex::ParticleReal * const AMREX_RESTRICT t_arr = soa[RealSoA::t].dataPtr(); amrex::ParticleReal * const AMREX_RESTRICT px_arr = soa[RealSoA::px].dataPtr(); amrex::ParticleReal * const AMREX_RESTRICT py_arr = soa[RealSoA::py].dataPtr(); amrex::ParticleReal * const AMREX_RESTRICT pt_arr = soa[RealSoA::pt].dataPtr(); amrex::ParticleReal * const AMREX_RESTRICT qm_arr = soa[RealSoA::qm].dataPtr(); amrex::ParticleReal * const AMREX_RESTRICT w_arr = soa[RealSoA::w ].dataPtr(); + uint64_t * const AMREX_RESTRICT idcpu_arr = particle_tile.GetStructOfArrays().GetIdCPUData().dataPtr(); + amrex::ParticleReal const * const AMREX_RESTRICT x_ptr = x.data(); amrex::ParticleReal const * const AMREX_RESTRICT y_ptr = y.data(); amrex::ParticleReal const * const AMREX_RESTRICT t_ptr = t.data(); @@ -175,12 +179,12 @@ namespace impactx amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE (int i) noexcept { - ParticleType& p = pstructs[old_np + i]; - p.id() = pid + i; - p.cpu() = cpuid; - p.pos(RealAoS::x) = x_ptr[i]; - p.pos(RealAoS::y) = y_ptr[i]; - p.pos(RealAoS::t) = t_ptr[i]; + amrex::ParticleIDWrapper{idcpu_arr[old_np+i]} = pid + i; + amrex::ParticleCPUWrapper{idcpu_arr[old_np+i]} = cpuid; + + x_arr[old_np+i] = x_ptr[i]; + y_arr[old_np+i] = y_ptr[i]; + t_arr[old_np+i] = t_ptr[i]; px_arr[old_np+i] = px_ptr[i]; py_arr[old_np+i] = py_ptr[i]; @@ -240,12 +244,6 @@ namespace impactx >(*this); } - std::vector - ImpactXParticleContainer::RealAoS_names () const - { - return get_RealAoS_names(); - } - std::vector ImpactXParticleContainer::RealSoA_names () const { @@ -264,17 +262,6 @@ namespace impactx m_coordsystem = coord_system; } - std::vector - get_RealAoS_names () - { - std::vector real_aos_names(RealAoS::names_s.size()); - - // compile-time attributes - std::copy(RealAoS::names_s.begin(), RealAoS::names_s.end(), real_aos_names.begin()); - - return real_aos_names; - } - std::vector get_RealSoA_names (int num_real_comps) { diff --git a/src/particles/diagnostics/DiagnosticOutput.cpp b/src/particles/diagnostics/DiagnosticOutput.cpp index 978208fc3..f9304afd1 100644 --- a/src/particles/diagnostics/DiagnosticOutput.cpp +++ b/src/particles/diagnostics/DiagnosticOutput.cpp @@ -11,13 +11,12 @@ #include "NonlinearLensInvariants.H" #include "ReducedBeamCharacteristics.H" -#include - #include // for BL_PROFILE #include // for AMREX_RESTRICT #include // for ParmParse #include // for ParticleReal #include // for PrintToFile +#include // for constructor of SoAParticle #include #include @@ -115,67 +114,58 @@ namespace impactx::diagnostics int const nLevel = tmp.finestLevel(); for (int lev = 0; lev <= nLevel; ++lev) { // loop over all particle boxes - using ParIt = typename decltype(tmp)::ParConstIterType; - for (ParIt pti(tmp, lev); pti.isValid(); ++pti) { + using MyPinnedParIter = amrex::ParIterSoA; + for (MyPinnedParIter pti(tmp, lev); pti.isValid(); ++pti) { const int np = pti.numParticles(); - // preparing access to particle data: AoS - using PType = ImpactXParticleContainer::ParticleType; - auto const &aos = pti.GetArrayOfStructs(); - PType const *const AMREX_RESTRICT aos_ptr = aos().dataPtr(); - // preparing access to particle data: SoA of Reals - auto &soa_real = pti.GetStructOfArrays().GetRealData(); - amrex::ParticleReal const *const AMREX_RESTRICT part_px = soa_real[RealSoA::px].dataPtr(); - amrex::ParticleReal const *const AMREX_RESTRICT part_py = soa_real[RealSoA::py].dataPtr(); - - // TODO: Refactor since this condition is always true here - if (otype == OutputType::PrintNonlinearLensInvariants) { - using namespace amrex::literals; + auto const& soa = pti.GetStructOfArrays(); + auto const& part_x = soa.GetRealData(RealSoA::x); + auto const& part_y = soa.GetRealData(RealSoA::y); + auto const& part_px = soa.GetRealData(RealSoA::px); + auto const& part_py = soa.GetRealData(RealSoA::py); - // Parse the diagnostic parameters - amrex::ParmParse pp_diag("diag"); + auto const& part_idcpu = soa.GetIdCPUData().dataPtr(); - amrex::ParticleReal alpha = 0.0; - pp_diag.queryAdd("alpha", alpha); + // Parse the diagnostic parameters + amrex::ParmParse pp_diag("diag"); - amrex::ParticleReal beta = 1.0; - pp_diag.queryAdd("beta", beta); + amrex::ParticleReal alpha = 0.0; + pp_diag.queryAdd("alpha", alpha); - amrex::ParticleReal tn = 0.4; - pp_diag.queryAdd("tn", tn); + amrex::ParticleReal beta = 1.0; + pp_diag.queryAdd("beta", beta); - amrex::ParticleReal cn = 0.01; - pp_diag.queryAdd("cn", cn); + amrex::ParticleReal tn = 0.4; + pp_diag.queryAdd("tn", tn); - NonlinearLensInvariants const nonlinear_lens_invariants(alpha, beta, tn, cn); + amrex::ParticleReal cn = 0.01; + pp_diag.queryAdd("cn", cn); - // print out particles - for (int i = 0; i < np; ++i) { + NonlinearLensInvariants const nonlinear_lens_invariants(alpha, beta, tn, cn); - // access AoS data such as positions and cpu/id - PType const &p = aos_ptr[i]; - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - uint64_t const global_id = ablastr::particles::localIDtoGlobal(p.id(), p.cpu()); + // print out particles (this hack works only on CPU and on GPUs with + // unified memory access) + for (int i = 0; i < np; ++i) { + amrex::ParticleReal const x = part_x[i]; + amrex::ParticleReal const y = part_y[i]; + uint64_t const global_id = part_idcpu[i]; - // access SoA Real data - amrex::ParticleReal const px = part_px[i]; - amrex::ParticleReal const py = part_py[i]; + amrex::ParticleReal const px = part_px[i]; + amrex::ParticleReal const py = part_py[i]; - // calculate invariants of motion - NonlinearLensInvariants::Data const HI_out = - nonlinear_lens_invariants(x, y, px, py); + // calculate invariants of motion + NonlinearLensInvariants::Data const HI_out = + nonlinear_lens_invariants(x, y, px, py); - // write particle invariant data to file - file_handler - << global_id << " " - << HI_out.H << " " << HI_out.I << "\n"; + // write particle invariant data to file + file_handler + << global_id << " " + << HI_out.H << " " << HI_out.I << "\n"; - } // i=0...np - } // if( otype == OutputType::PrintInvariants) + } // i=0...np } // end loop over all particle boxes - } // env mesh-refinement level loop + } // end mesh-refinement level loop } } diff --git a/src/particles/diagnostics/ReducedBeamCharacteristics.cpp b/src/particles/diagnostics/ReducedBeamCharacteristics.cpp index 86f79d7be..063c7c3c2 100644 --- a/src/particles/diagnostics/ReducedBeamCharacteristics.cpp +++ b/src/particles/diagnostics/ReducedBeamCharacteristics.cpp @@ -33,7 +33,7 @@ namespace impactx::diagnostics // reference particle charge in C amrex::ParticleReal const q_C = ref_part.charge; - // preparing access to particle data: AoS and SoA + // preparing access to particle data: SoA using PType = typename ImpactXParticleContainer::SuperParticleType; // prepare reduction operations for calculation of mean and min/max values in 6D phase space @@ -78,10 +78,10 @@ namespace impactx::diagnostics amrex::ParticleReal, amrex::ParticleReal > { - // access AoS particle position data - const amrex::ParticleReal p_pos0 = p.pos(0); - const amrex::ParticleReal p_pos1 = p.pos(1); - const amrex::ParticleReal p_pos2 = p.pos(2); + // access particle position data + const amrex::ParticleReal p_x = p.rdata(RealSoA::x); + const amrex::ParticleReal p_y = p.rdata(RealSoA::y); + const amrex::ParticleReal p_t = p.rdata(RealSoA::t); // access SoA particle momentum data and weighting const amrex::ParticleReal p_w = p.rdata(RealSoA::w); @@ -90,20 +90,20 @@ namespace impactx::diagnostics const amrex::ParticleReal p_pt = p.rdata(RealSoA::pt); // prepare mean position values - const amrex::ParticleReal p_x_mean = p_pos0*p_w; - const amrex::ParticleReal p_y_mean = p_pos1*p_w; - const amrex::ParticleReal p_t_mean = p_pos2*p_w; + const amrex::ParticleReal p_x_mean = p_x * p_w; + const amrex::ParticleReal p_y_mean = p_y * p_w; + const amrex::ParticleReal p_t_mean = p_t * p_w; - const amrex::ParticleReal p_px_mean = p_px*p_w; - const amrex::ParticleReal p_py_mean = p_py*p_w; - const amrex::ParticleReal p_pt_mean = p_pt*p_w; + const amrex::ParticleReal p_px_mean = p_px * p_w; + const amrex::ParticleReal p_py_mean = p_py * p_w; + const amrex::ParticleReal p_pt_mean = p_pt * p_w; return {p_w, p_x_mean, p_y_mean, p_t_mean, p_px_mean, p_py_mean, p_pt_mean, - p_pos0, p_pos0, - p_pos1, p_pos1, - p_pos2, p_pos2, + p_x, p_x, + p_y, p_y, + p_t, p_t, p_px, p_px, p_py, p_py, p_pt, p_pt}; @@ -199,13 +199,10 @@ namespace impactx::diagnostics const amrex::ParticleReal p_px = p.rdata(RealSoA::px); const amrex::ParticleReal p_py = p.rdata(RealSoA::py); const amrex::ParticleReal p_pt = p.rdata(RealSoA::pt); - // access AoS particle position data - const amrex::ParticleReal p_pos0 = p.pos(0); - const amrex::ParticleReal p_pos1 = p.pos(1); - const amrex::ParticleReal p_pos2 = p.pos(2); - const amrex::ParticleReal p_x = p_pos0; - const amrex::ParticleReal p_y = p_pos1; - const amrex::ParticleReal p_t = p_pos2; + // access position data + const amrex::ParticleReal p_x = p.rdata(RealSoA::x); + const amrex::ParticleReal p_y = p.rdata(RealSoA::y); + const amrex::ParticleReal p_t = p.rdata(RealSoA::t); // prepare mean square for positions const amrex::ParticleReal p_x_ms = (p_x-x_mean)*(p_x-x_mean)*p_w; const amrex::ParticleReal p_y_ms = (p_y-y_mean)*(p_y-y_mean)*p_w; diff --git a/src/particles/elements/Aperture.H b/src/particles/elements/Aperture.H index d42c79f82..b8d7a8f78 100644 --- a/src/particles/elements/Aperture.H +++ b/src/particles/elements/Aperture.H @@ -17,7 +17,6 @@ #include "mixin/thin.H" #include "mixin/nofinalize.H" -#include #include #include @@ -69,29 +68,33 @@ namespace impactx /** This is an aperture functor, so that a variable of this type can be used like an * aperture function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t (unused) * @param px particle momentum in x * @param py particle momentum in y - * @param pt particle momentum in t + * @param pt particle momentum in t (unused) + * @param idcpu particle global index * @param refpart reference particle (unused) */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, - [[maybe_unused]] amrex::ParticleReal & AMREX_RESTRICT px, - [[maybe_unused]] amrex::ParticleReal & AMREX_RESTRICT py, + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + [[maybe_unused]] amrex::ParticleReal & AMREX_RESTRICT t, + amrex::ParticleReal & AMREX_RESTRICT px, + amrex::ParticleReal & AMREX_RESTRICT py, [[maybe_unused]] amrex::ParticleReal & AMREX_RESTRICT pt, - [[maybe_unused]] RefPart const & refpart) const + uint64_t & AMREX_RESTRICT idcpu, + [[maybe_unused]] RefPart const & refpart + ) const { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_in(x, y, px, py); - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - auto const id = p.id(); + auto const id = amrex::ParticleIDWrapper{idcpu}; // scale horizontal and vertical coordinates amrex::ParticleReal const u = x / m_xmax; @@ -102,19 +105,19 @@ namespace impactx { case Shape::rectangular : // default if (pow(u,2)>1 || pow(v,2) > 1_prt) { - p.id() = -id; + amrex::ParticleIDWrapper{idcpu} = -id; } break; case Shape::elliptical : if (pow(u,2)+pow(v,2) > 1_prt) { - p.id() = -id; + amrex::ParticleIDWrapper{idcpu} = -id; } break; } // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. */ diff --git a/src/particles/elements/Buncher.H b/src/particles/elements/Buncher.H index b8dda736b..e23b307d6 100644 --- a/src/particles/elements/Buncher.H +++ b/src/particles/elements/Buncher.H @@ -59,29 +59,30 @@ namespace impactx /** This is a buncher functor, so that a variable of this type can be used like a * buncher function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, amrex::ParticleReal & AMREX_RESTRICT px, amrex::ParticleReal & AMREX_RESTRICT py, amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, RefPart const & refpart) const { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); - - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + shift_in(x, y, px, py); // access reference particle values to find (beta*gamma)^2 amrex::ParticleReal const pt_ref = refpart.pt; @@ -93,13 +94,8 @@ namespace impactx amrex::ParticleReal ptout = pt; // advance position and momentum - p.pos(RealAoS::x) = x; pxout = px + m_k*m_V/(2.0_prt*betgam2)*x; - - p.pos(RealAoS::y) = y; pyout = py + m_k*m_V/(2.0_prt*betgam2)*y; - - p.pos(RealAoS::t) = t; ptout = pt - m_k*m_V*t; // assign updated momenta @@ -108,7 +104,7 @@ namespace impactx pt = ptout; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. */ diff --git a/src/particles/elements/CFbend.H b/src/particles/elements/CFbend.H index b925a5ab6..069e86b21 100644 --- a/src/particles/elements/CFbend.H +++ b/src/particles/elements/CFbend.H @@ -68,29 +68,36 @@ namespace impactx /** This is a cfbend functor, so that a variable of this type can be used like a cfbend function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, - amrex::ParticleReal & AMREX_RESTRICT px, - amrex::ParticleReal & AMREX_RESTRICT py, - amrex::ParticleReal & AMREX_RESTRICT pt, - RefPart const & refpart) const + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, + amrex::ParticleReal & AMREX_RESTRICT px, + amrex::ParticleReal & AMREX_RESTRICT py, + amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, + RefPart const & refpart + ) const { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_in(x, y, px, py); - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + // initialize output values + amrex::ParticleReal xout = x; + amrex::ParticleReal yout = y; + amrex::ParticleReal tout = t; // initialize output values of momenta amrex::ParticleReal pxout = px; @@ -116,11 +123,11 @@ namespace impactx + (sinx - omegax*slice_ds)/(gx*omegax*pow(bet,2)*pow(m_rc,2)); // advance position and momentum (focusing) - p.pos(RealAoS::x) = cosx*x + sinx/omegax*px - (1.0_prt - cosx)/(gx*bet*m_rc)*pt; - pxout = -omegax*sinx*x + cosx*px - sinx/(omegax*bet*m_rc)*pt; + x = cosx*xout + sinx/omegax*px - (1.0_prt - cosx)/(gx*bet*m_rc)*pt; + pxout = -omegax*sinx*xout + cosx*px - sinx/(omegax*bet*m_rc)*pt; - p.pos(RealAoS::t) = sinx/(omegax*bet*m_rc)*x + (1.0_prt - cosx)/(gx*bet*m_rc)*px - + t + r56*pt; + y = sinx/(omegax*bet*m_rc)*xout + (1.0_prt - cosx)/(gx*bet*m_rc)*px + + tout + r56*pt; ptout = pt; } else { // calculate expensive terms once @@ -130,11 +137,11 @@ namespace impactx + (sinhx - omegax*slice_ds)/(gx*omegax*pow(bet,2)*pow(m_rc,2)); // advance position and momentum (defocusing) - p.pos(RealAoS::x) = coshx*x + sinhx/omegax*px - (1.0_prt - coshx)/(gx*bet*m_rc)*pt; - pxout = omegax*sinhx*x + coshx*px - sinhx/(omegax*bet*m_rc)*pt; + x = coshx*xout + sinhx/omegax*px - (1.0_prt - coshx)/(gx*bet*m_rc)*pt; + pxout = omegax*sinhx*xout + coshx*px - sinhx/(omegax*bet*m_rc)*pt; - p.pos(RealAoS::t) = sinhx/(omegax*bet*m_rc)*x + (1.0_prt - coshx)/(gx*bet*m_rc)*px - + t + r56*pt; + t = sinhx/(omegax*bet*m_rc)*xout + (1.0_prt - coshx)/(gx*bet*m_rc)*px + + tout + r56*pt; ptout = pt; } @@ -147,8 +154,8 @@ namespace impactx auto const [siny, cosy] = amrex::Math::sincos(omegay * slice_ds); // advance position and momentum (focusing) - p.pos(RealAoS::y) = cosy*y + siny/omegay*py; - pyout = -omegay*siny*y + cosy*py; + y = cosy*yout + siny/omegay*py; + pyout = -omegay*siny*yout + cosy*py; } else { // calculate expensive terms once @@ -156,8 +163,8 @@ namespace impactx amrex::ParticleReal const coshy = cosh(omegay * slice_ds); // advance position and momentum (defocusing) - p.pos(RealAoS::y) = coshy*y + sinhy/omegay*py; - pyout = omegay*sinhy*y + coshy*py; + y = coshy*yout + sinhy/omegay*py; + pyout = omegay*sinhy*yout + coshy*py; } // assign updated momenta @@ -166,7 +173,7 @@ namespace impactx pt = ptout; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. diff --git a/src/particles/elements/ChrDrift.H b/src/particles/elements/ChrDrift.H index 92f940ce3..337a3eca7 100644 --- a/src/particles/elements/ChrDrift.H +++ b/src/particles/elements/ChrDrift.H @@ -61,29 +61,36 @@ namespace impactx /** This is a chrdrift functor, so that a variable of this type can be used like a chrdrift function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, - amrex::ParticleReal & AMREX_RESTRICT px, - amrex::ParticleReal & AMREX_RESTRICT py, - amrex::ParticleReal & AMREX_RESTRICT pt, - RefPart const & refpart) const { - + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, + amrex::ParticleReal & AMREX_RESTRICT px, + amrex::ParticleReal & AMREX_RESTRICT py, + amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, + RefPart const & refpart + ) const + { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_in(x, y, px, py); - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + // access position data + amrex::ParticleReal const xout = x; + amrex::ParticleReal const yout = y; + amrex::ParticleReal const tout = t; // initialize output values of momenta amrex::ParticleReal const pxout = px; @@ -102,9 +109,9 @@ namespace impactx delta1 = sqrt(1_prt - 2_prt*pt/bet + pow(pt,2)); // advance transverse position and momentum (drift) - p.pos(RealAoS::x) = x + slice_ds * px / delta1; + x = xout + slice_ds * px / delta1; // pxout = px; - p.pos(RealAoS::y) = y + slice_ds * py / delta1; + y = yout + slice_ds * py / delta1; // pyout = py; // the corresponding symplectic update to t @@ -113,7 +120,7 @@ namespace impactx term = -2_prt + pow(gam,2)*term; term = (-1_prt+bet*pt)*term; term = term/(2_prt*pow(bet,3)*pow(gam,2)); - p.pos(RealAoS::t) = t - slice_ds*(1_prt/bet + term/pow(delta1,3)); + t = tout - slice_ds * (1_prt / bet + term / pow(delta1, 3)); // ptout = pt; // assign updated momenta @@ -122,7 +129,7 @@ namespace impactx pt = ptout; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. @@ -130,8 +137,8 @@ namespace impactx * @param[in,out] refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - void operator() (RefPart & AMREX_RESTRICT refpart) const { - + void operator() (RefPart & AMREX_RESTRICT refpart) const + { using namespace amrex::literals; // for _rt and _prt // assign input reference particle values @@ -159,7 +166,6 @@ namespace impactx // advance integrated path length refpart.s = s + slice_ds; - } }; diff --git a/src/particles/elements/ChrQuad.H b/src/particles/elements/ChrQuad.H index 88901aabb..02463f46e 100644 --- a/src/particles/elements/ChrQuad.H +++ b/src/particles/elements/ChrQuad.H @@ -73,30 +73,36 @@ namespace impactx /** This is a chrquad functor, so that a variable of this type can be used like a chrquad function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, amrex::ParticleReal & AMREX_RESTRICT px, amrex::ParticleReal & AMREX_RESTRICT py, amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, RefPart const & refpart ) const { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_in(x, y, px, py); - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + // access position data + amrex::ParticleReal const xout = x; + amrex::ParticleReal const yout = y; + amrex::ParticleReal const tout = t; // length of the current slice amrex::ParticleReal const slice_ds = m_ds / nslice(); @@ -125,33 +131,33 @@ namespace impactx amrex::ParticleReal const ptout = pt; // paceholder variables - amrex::ParticleReal q1 = x; - amrex::ParticleReal q2 = y; + amrex::ParticleReal q1 = xout; + amrex::ParticleReal q2 = yout; amrex::ParticleReal p1 = px; amrex::ParticleReal p2 = py; if(g > 0.0) { // advance transverse position and momentum (focusing quad) - p.pos(RealAoS::x) = cos(omega*slice_ds)*x + + x = cos(omega*slice_ds) * xout + sin(omega*slice_ds)/(omega*delta1)*px; - pxout = -omega*delta1*sin(omega*slice_ds)*x + cos(omega*slice_ds)*px; + pxout = -omega * delta1 * sin(omega*slice_ds) * xout + cos(omega * slice_ds) * px; - p.pos(RealAoS::y) = cosh(omega*slice_ds)*y + + y = cosh(omega*slice_ds) * yout + sinh(omega*slice_ds)/(omega*delta1)*py; - pyout = omega*delta1*sinh(omega*slice_ds)*y + cosh(omega*slice_ds)*py; + pyout = omega * delta1 * sinh(omega*slice_ds) * yout + cosh(omega * slice_ds) * py; } else { // advance transverse position and momentum (defocusing quad) - p.pos(RealAoS::x) = cosh(omega*slice_ds)*x + + x = cosh(omega*slice_ds) * xout + sinh(omega*slice_ds)/(omega*delta1)*px; - pxout = omega*delta1*sinh(omega*slice_ds)*x + cosh(omega*slice_ds)*px; + pxout = omega * delta1 * sinh(omega*slice_ds) * xout + cosh(omega * slice_ds) * px; - p.pos(RealAoS::y) = cos(omega*slice_ds)*y + + y = cos(omega*slice_ds) * yout + sin(omega*slice_ds)/(omega*delta1)*py; - pyout = -omega*delta1*sin(omega*slice_ds)*y + cos(omega*slice_ds)*py; + pyout = -omega * delta1 * sin(omega*slice_ds) * yout + cos(omega * slice_ds) * py; - q1 = y; - q2 = x; + q1 = yout; + q2 = xout; p1 = py; p2 = px; } @@ -160,7 +166,7 @@ namespace impactx // the corresponding symplectic update to t amrex::ParticleReal const term = pt + delta/bet; - amrex::ParticleReal const t0 = t - term*slice_ds/delta1; + amrex::ParticleReal const t0 = tout - term * slice_ds / delta1; amrex::ParticleReal const w = omega*delta1; amrex::ParticleReal const term1 = -(pow(p2,2)+pow(q2,2)*pow(w,2))*sinh(2_prt*slice_ds*omega); @@ -169,8 +175,8 @@ namespace impactx amrex::ParticleReal const term4 = -2_prt*q1*p1*w*cos(2_prt*slice_ds*omega); amrex::ParticleReal const term5 = 2_prt*omega*(q1*p1*delta1 + q2*p2*delta1 -(pow(p1,2)+pow(p2,2))*slice_ds - (pow(q1,2)-pow(q2,2))*pow(w,2)*slice_ds); - p.pos(RealAoS::t) = t0 + (-1_prt+bet*pt)/(8_prt*bet*pow(delta1,3)*omega) - *(term1+term2+term3+term4+term5); + t = t0 + (-1_prt+bet*pt)/(8_prt*bet*pow(delta1,3)*omega) + *(term1+term2+term3+term4+term5); // ptout = pt; @@ -180,7 +186,7 @@ namespace impactx pt = ptout; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. diff --git a/src/particles/elements/ChrUniformAcc.H b/src/particles/elements/ChrUniformAcc.H index f1734168f..dfbdd794f 100644 --- a/src/particles/elements/ChrUniformAcc.H +++ b/src/particles/elements/ChrUniformAcc.H @@ -68,30 +68,31 @@ namespace impactx /** This is a chracc functor, so that a variable of this type can be used like a chracc function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, amrex::ParticleReal & AMREX_RESTRICT px, amrex::ParticleReal & AMREX_RESTRICT py, amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, RefPart const & refpart ) const { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); - - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + shift_in(x, y, px, py); // length of the current slice amrex::ParticleReal const slice_ds = m_ds / nslice(); @@ -148,13 +149,13 @@ namespace impactx pt = ptout; // advance positions and momenta using map for rotation - p.pos(RealAoS::x) = cos(theta)*xout + sin(theta)*yout; + x = cos(theta)*xout + sin(theta)*yout; pxout = cos(theta)*px + sin(theta)*py; - p.pos(RealAoS::y) = -sin(theta)*xout + cos(theta)*yout; + y = -sin(theta)*xout + cos(theta)*yout; pyout = -sin(theta)*px + cos(theta)*py; - p.pos(RealAoS::t) = tout; + t = tout; ptout = pt; // assign updated momenta @@ -168,7 +169,7 @@ namespace impactx pt = pt/bgf; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. diff --git a/src/particles/elements/ConstF.H b/src/particles/elements/ConstF.H index c0dae3af1..41d390a38 100644 --- a/src/particles/elements/ConstF.H +++ b/src/particles/elements/ConstF.H @@ -65,35 +65,39 @@ namespace impactx /** This pushes a single particle, relative to the reference particle * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, amrex::ParticleReal & AMREX_RESTRICT px, amrex::ParticleReal & AMREX_RESTRICT py, amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, RefPart const & refpart) const { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); - - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + shift_in(x, y, px, py); // access reference particle values to find beta*gamma^2 amrex::ParticleReal const pt_ref = refpart.pt; amrex::ParticleReal const betgam2 = pow(pt_ref, 2) - 1.0_prt; - // intialize output values of momenta + // intialize output values + amrex::ParticleReal xout = x; + amrex::ParticleReal yout = y; + amrex::ParticleReal tout = t; amrex::ParticleReal pxout = px; amrex::ParticleReal pyout = py; amrex::ParticleReal ptout = pt; @@ -102,22 +106,25 @@ namespace impactx amrex::ParticleReal const slice_ds = m_ds / nslice(); // advance position and momentum - p.pos(RealAoS::x) = cos(m_kx*slice_ds)*x + sin(m_kx*slice_ds)/m_kx*px; + xout = cos(m_kx*slice_ds)*x + sin(m_kx*slice_ds)/m_kx*px; pxout = -m_kx*sin(m_kx*slice_ds)*x + cos(m_kx*slice_ds)*px; - p.pos(RealAoS::y) = cos(m_ky*slice_ds)*y + sin(m_ky*slice_ds)/m_ky*py; + yout = cos(m_ky*slice_ds)*y + sin(m_ky*slice_ds)/m_ky*py; pyout = -m_ky*sin(m_ky*slice_ds)*y + cos(m_ky*slice_ds)*py; - p.pos(RealAoS::t) = cos(m_kt*slice_ds)*t + sin(m_kt*slice_ds)/(betgam2*m_kt)*pt; + tout = cos(m_kt*slice_ds)*t + sin(m_kt*slice_ds)/(betgam2*m_kt)*pt; ptout = -(m_kt*betgam2)*sin(m_kt*slice_ds)*t + cos(m_kt*slice_ds)*pt; - // assign updated momenta + // assign updated values + x = xout; + y = yout; + t = tout; px = pxout; py = pyout; pt = ptout; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. diff --git a/src/particles/elements/DipEdge.H b/src/particles/elements/DipEdge.H index 9a297cf30..ba9410f4f 100644 --- a/src/particles/elements/DipEdge.H +++ b/src/particles/elements/DipEdge.H @@ -70,30 +70,30 @@ namespace impactx /** This is a dipedge functor, so that a variable of this type can be used like a * dipedge function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t (unchanged) + * @param idcpu particle global index (unused) * @param refpart reference particle (unused) */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + [[maybe_unused]] amrex::ParticleReal & AMREX_RESTRICT t, amrex::ParticleReal & AMREX_RESTRICT px, amrex::ParticleReal & AMREX_RESTRICT py, [[maybe_unused]] amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, [[maybe_unused]] RefPart const & refpart) const { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); - - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - - // access reference particle values if needed + shift_in(x, y, px, py); // edge focusing matrix elements (zero gap) amrex::ParticleReal const R21 = tan(m_psi)/m_rc; @@ -110,7 +110,7 @@ namespace impactx py = py + R43*y; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. */ diff --git a/src/particles/elements/Drift.H b/src/particles/elements/Drift.H index 607ac057f..d4ef04acb 100644 --- a/src/particles/elements/Drift.H +++ b/src/particles/elements/Drift.H @@ -58,34 +58,39 @@ namespace impactx /** This is a drift functor, so that a variable of this type can be used like a drift function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, - amrex::ParticleReal & AMREX_RESTRICT px, - amrex::ParticleReal & AMREX_RESTRICT py, - amrex::ParticleReal & AMREX_RESTRICT pt, - RefPart const & refpart) const { - + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, + amrex::ParticleReal & AMREX_RESTRICT px, + amrex::ParticleReal & AMREX_RESTRICT py, + amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, + RefPart const & refpart + ) const + { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); - - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + shift_in(x, y, px, py); - // initialize output values of momenta - amrex::ParticleReal const pxout = px; - amrex::ParticleReal const pyout = py; - amrex::ParticleReal const ptout = pt; + // intialize output values + amrex::ParticleReal xout = x; + amrex::ParticleReal yout = y; + amrex::ParticleReal tout = t; + amrex::ParticleReal pxout = px; + amrex::ParticleReal pyout = py; + amrex::ParticleReal ptout = pt; // length of the current slice amrex::ParticleReal const slice_ds = m_ds / nslice(); @@ -95,20 +100,23 @@ namespace impactx amrex::ParticleReal const betgam2 = pow(pt_ref, 2) - 1.0_prt; // advance position and momentum (drift) - p.pos(RealAoS::x) = x + slice_ds * px; + xout = x + slice_ds * px; // pxout = px; - p.pos(RealAoS::y) = y + slice_ds * py; + yout = y + slice_ds * py; // pyout = py; - p.pos(RealAoS::t) = t + (slice_ds/betgam2) * pt; + tout = t + (slice_ds/betgam2) * pt; // ptout = pt; - // assign updated momenta + // assign updated values + x = xout; + y = yout; + t = tout; px = pxout; py = pyout; pt = ptout; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. diff --git a/src/particles/elements/ExactDrift.H b/src/particles/elements/ExactDrift.H index 19f00abe0..2e732130f 100644 --- a/src/particles/elements/ExactDrift.H +++ b/src/particles/elements/ExactDrift.H @@ -59,29 +59,36 @@ namespace impactx /** This is an exactdrift functor, so that a variable of this type can be used like * an exactdrift function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, amrex::ParticleReal & AMREX_RESTRICT px, amrex::ParticleReal & AMREX_RESTRICT py, amrex::ParticleReal & AMREX_RESTRICT pt, - RefPart const & refpart) const { - + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, + RefPart const & refpart + ) const + { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_in(x, y, px, py); - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + // initialize output values + amrex::ParticleReal const xout = x; + amrex::ParticleReal const yout = y; + amrex::ParticleReal const tout = t; // initialize output values of momenta amrex::ParticleReal const pxout = px; @@ -100,12 +107,11 @@ namespace impactx 1_prt/pow(betgam,2) - pow(px,2) - pow(py,2)); // advance position and momentum (exact drift) - p.pos(RealAoS::x) = x + slice_ds * px / pzden; + x = xout + slice_ds * px / pzden; // pxout = px; - p.pos(RealAoS::y) = y + slice_ds * py / pzden; + y = yout + slice_ds * py / pzden; // pyout = py; - p.pos(RealAoS::t) = t - slice_ds * (1_prt/bet + - (pt-1_prt/bet)/pzden); + t = tout - slice_ds * (1_prt / bet + (pt-1_prt/bet)/pzden); // ptout = pt; // assign updated momenta @@ -114,7 +120,7 @@ namespace impactx pt = ptout; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. diff --git a/src/particles/elements/ExactSbend.H b/src/particles/elements/ExactSbend.H index 0b8ee014b..8d9501f11 100644 --- a/src/particles/elements/ExactSbend.H +++ b/src/particles/elements/ExactSbend.H @@ -76,29 +76,36 @@ namespace impactx /** This is an ExactSbend functor, so that a variable of this type can be used like an ExactSbend function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, - amrex::ParticleReal & AMREX_RESTRICT px, - amrex::ParticleReal & AMREX_RESTRICT py, - amrex::ParticleReal & AMREX_RESTRICT pt, - RefPart const & refpart) const { - + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, + amrex::ParticleReal & AMREX_RESTRICT px, + amrex::ParticleReal & AMREX_RESTRICT py, + amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, + RefPart const & refpart + ) const + { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_in(x, y, px, py); - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + // access position data + amrex::ParticleReal const xout = x; + amrex::ParticleReal const yout = y; + amrex::ParticleReal const tout = t; // angle of arc for the current slice amrex::ParticleReal const slice_phi = m_phi / nslice(); @@ -117,7 +124,7 @@ namespace impactx // assign intermediate quantities amrex::ParticleReal const pperp = sqrt(pow(pt,2)-2.0_prt/bet*pt-pow(py,2)+1.0_prt); amrex::ParticleReal const pzi = sqrt(pow(pperp,2)-pow(px,2)); - amrex::ParticleReal const rho = rc + x; + amrex::ParticleReal const rho = rc + xout; auto const [sin_phi, cos_phi] = amrex::Math::sincos(slice_phi); // update momenta @@ -130,9 +137,9 @@ namespace impactx amrex::ParticleReal const theta = slice_phi + asin(px/pperp) - asin(pxout/pperp); // update position coordinates - p.pos(RealAoS::x) = -rc + rho*cos_phi + rc*(pzf + px*sin_phi - pzi*cos_phi); - p.pos(RealAoS::y) = y + theta*rc*py; - p.pos(RealAoS::t) = t - theta*rc*(pt - 1.0_prt/bet) - m_phi*rc/bet; + x = -rc + rho*cos_phi + rc*(pzf + px*sin_phi - pzi*cos_phi); + y = yout + theta * rc * py; + t = tout - theta * rc * (pt - 1.0_prt / bet) - m_phi * rc / bet; // assign updated momenta px = pxout; @@ -140,7 +147,7 @@ namespace impactx pt = ptout; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. @@ -148,8 +155,8 @@ namespace impactx * @param[in,out] refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - void operator() (RefPart & AMREX_RESTRICT refpart) const { - + void operator() (RefPart & AMREX_RESTRICT refpart) const + { using namespace amrex::literals; // for _rt and _prt // assign input reference particle values diff --git a/src/particles/elements/Kicker.H b/src/particles/elements/Kicker.H index bd735747b..021298646 100644 --- a/src/particles/elements/Kicker.H +++ b/src/particles/elements/Kicker.H @@ -68,29 +68,36 @@ namespace impactx /** This is a transverse kicker functor, so that a variable of this type can be * used like a function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle (unused) */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, amrex::ParticleReal & AMREX_RESTRICT px, amrex::ParticleReal & AMREX_RESTRICT py, amrex::ParticleReal & AMREX_RESTRICT pt, - [[maybe_unused]] RefPart const & refpart) const + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, + [[maybe_unused]] RefPart const & refpart + ) const { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_in(x, y, px, py); - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + // access position data + amrex::ParticleReal const xout = x; + amrex::ParticleReal const yout = y; + amrex::ParticleReal const tout = t; // normalize quad units to MAD-X convention if needed amrex::ParticleReal dpx = m_xkick; @@ -106,13 +113,13 @@ namespace impactx amrex::ParticleReal ptout = pt; // advance position and momentum - p.pos(RealAoS::x) = x; + x = xout; pxout = px + dpx; - p.pos(RealAoS::y) = y; + y = yout; pyout = py + dpy; - p.pos(RealAoS::t) = t; + t = tout; ptout = pt; // assign updated momenta @@ -121,7 +128,7 @@ namespace impactx pt = ptout; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. */ diff --git a/src/particles/elements/Multipole.H b/src/particles/elements/Multipole.H index a198b8fa8..12389098d 100644 --- a/src/particles/elements/Multipole.H +++ b/src/particles/elements/Multipole.H @@ -67,38 +67,43 @@ namespace impactx /** This is a multipole functor, so that a variable of this type can be used like a * multipole function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle (unused) */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + [[maybe_unused]] amrex::ParticleReal & AMREX_RESTRICT t, amrex::ParticleReal & AMREX_RESTRICT px, amrex::ParticleReal & AMREX_RESTRICT py, amrex::ParticleReal & AMREX_RESTRICT pt, - [[maybe_unused]] RefPart const & refpart) const { - + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, + [[maybe_unused]] RefPart const & refpart + ) const + { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_in(x, y, px, py); // a complex type with two amrex::ParticleReal using Complex = amrex::GpuComplex; - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); - // access reference particle values to find (beta*gamma)^2 //amrex::ParticleReal const pt_ref = refpart.pt; //amrex::ParticleReal const betgam2 = pow(pt_ref, 2) - 1.0_prt; - // intialize output values of momenta + // intialize output values + amrex::ParticleReal xout = x; + amrex::ParticleReal yout = y; + amrex::ParticleReal tout = t; amrex::ParticleReal pxout = px; amrex::ParticleReal pyout = py; amrex::ParticleReal ptout = pt; @@ -115,22 +120,25 @@ namespace impactx amrex::ParticleReal const dpy = kick.m_imag/m_mfactorial; // advance position and momentum - p.pos(RealAoS::x) = x; + // xout = x; pxout = px + dpx; - p.pos(RealAoS::y) = y; + // yout = y; pyout = py + dpy; - p.pos(RealAoS::t) = t; + // tout = t; ptout = pt; - // assign updated momenta + // assign updated values + x = xout; + y = yout; + t = tout; px = pxout; py = pyout; pt = ptout; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. */ diff --git a/src/particles/elements/None.H b/src/particles/elements/None.H index ecad8223f..2f09d4935 100644 --- a/src/particles/elements/None.H +++ b/src/particles/elements/None.H @@ -51,19 +51,26 @@ namespace impactx /** Does nothing to a particle. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - [[maybe_unused]] PType& AMREX_RESTRICT p, + [[maybe_unused]] amrex::ParticleReal & AMREX_RESTRICT x, + [[maybe_unused]] amrex::ParticleReal & AMREX_RESTRICT y, + [[maybe_unused]] amrex::ParticleReal & AMREX_RESTRICT t, [[maybe_unused]] amrex::ParticleReal & AMREX_RESTRICT px, [[maybe_unused]] amrex::ParticleReal & AMREX_RESTRICT py, [[maybe_unused]] amrex::ParticleReal & AMREX_RESTRICT pt, - [[maybe_unused]] RefPart const & refpart) const + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, + [[maybe_unused]] RefPart const & refpart + ) const { // nothing to do } diff --git a/src/particles/elements/NonlinearLens.H b/src/particles/elements/NonlinearLens.H index 0bae32d3d..3a48a5430 100644 --- a/src/particles/elements/NonlinearLens.H +++ b/src/particles/elements/NonlinearLens.H @@ -64,38 +64,43 @@ namespace impactx /** This is a nonlinear lens functor, so that a variable of this type * can be used like a nonlinear lens function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle (unused) */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, - amrex::ParticleReal & AMREX_RESTRICT px, - amrex::ParticleReal & AMREX_RESTRICT py, - amrex::ParticleReal & AMREX_RESTRICT pt, - [[maybe_unused]] RefPart const & refpart) const { - + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + [[maybe_unused]] amrex::ParticleReal & AMREX_RESTRICT t, + amrex::ParticleReal & AMREX_RESTRICT px, + amrex::ParticleReal & AMREX_RESTRICT py, + amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, + [[maybe_unused]] RefPart const & refpart + ) const + { using namespace amrex::literals; // for _rt and _prt // a complex type with two amrex::ParticleReal using Complex = amrex::GpuComplex; // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); - - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + shift_in(x, y, px, py); // access reference particle values to find (beta*gamma)^2 //amrex::ParticleReal const pt_ref = refpart.pt; //amrex::ParticleReal const betgam2 = pow(pt_ref, 2) - 1.0_prt; - // intialize output values of momenta + // intialize output values + amrex::ParticleReal xout = x; + amrex::ParticleReal yout = y; + amrex::ParticleReal tout = t; amrex::ParticleReal pxout = px; amrex::ParticleReal pyout = py; amrex::ParticleReal ptout = pt; @@ -125,22 +130,25 @@ namespace impactx amrex::ParticleReal const dpy = -kick*dF.m_imag; // advance position and momentum - p.pos(RealAoS::x) = x; + // xout = x; pxout = px + dpx; - p.pos(RealAoS::y) = y; + // yout = y; pyout = py + dpy; - p.pos(RealAoS::t) = t; + // tout = t; ptout = pt; - // assign updated momenta + // assign updated values + x = xout; + y = yout; + t = tout; px = pxout; py = pyout; pt = ptout; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. */ diff --git a/src/particles/elements/PRot.H b/src/particles/elements/PRot.H index bffb4d479..5a51fbaf2 100644 --- a/src/particles/elements/PRot.H +++ b/src/particles/elements/PRot.H @@ -58,32 +58,36 @@ namespace impactx /** This is a prot functor, so that a variable of this type can be used like a * prot function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType & AMREX_RESTRICT p, + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, amrex::ParticleReal & AMREX_RESTRICT px, amrex::ParticleReal & AMREX_RESTRICT py, amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, RefPart const & refpart ) const { using namespace amrex::literals; // for _rt and _prt - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); - // access reference particle values to find beta: amrex::ParticleReal const beta = refpart.beta(); - // intialize output values of momenta + // intialize output values + amrex::ParticleReal xout = x; + amrex::ParticleReal yout = y; + amrex::ParticleReal tout = t; amrex::ParticleReal pxout = px; amrex::ParticleReal pyout = py; amrex::ParticleReal ptout = pt; @@ -98,16 +102,19 @@ namespace impactx amrex::ParticleReal const pzf = pz*cos_theta - (px + sin_phi_in)*sin_theta; // advance position and momentum - p.pos(RealAoS::x) = x*pz/pzf; + xout = x*pz/pzf; pxout = px*cos_theta + (pz - cos_phi_in)*sin_theta; - p.pos(RealAoS::y) = y + py*x*sin_theta/pzf; + yout = y + py*x*sin_theta/pzf; pyout = py; - p.pos(RealAoS::t) = t - (pt - 1.0_prt/beta)*x*sin_theta/pzf; + tout = t - (pt - 1.0_prt/beta)*x*sin_theta/pzf; ptout = pt; - // assign updated momenta + // assign updated values + x = xout; + y = yout; + t = tout; px = pxout; py = pyout; pt = ptout; diff --git a/src/particles/elements/Quad.H b/src/particles/elements/Quad.H index 75d82ae7b..d8c78c806 100644 --- a/src/particles/elements/Quad.H +++ b/src/particles/elements/Quad.H @@ -64,30 +64,31 @@ namespace impactx /** This is a quad functor, so that a variable of this type can be used like a quad function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, amrex::ParticleReal & AMREX_RESTRICT px, amrex::ParticleReal & AMREX_RESTRICT py, amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, RefPart const & refpart ) const { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); - - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + shift_in(x, y, px, py); // length of the current slice amrex::ParticleReal const slice_ds = m_ds / nslice(); @@ -99,40 +100,46 @@ namespace impactx // compute phase advance per unit length in s (in rad/m) amrex::ParticleReal const omega = sqrt(std::abs(m_k)); - // intialize output values of momenta + // intialize output values + amrex::ParticleReal xout = x; + amrex::ParticleReal yout = y; + amrex::ParticleReal tout = t; amrex::ParticleReal pxout = px; amrex::ParticleReal pyout = py; amrex::ParticleReal const ptout = pt; if(m_k > 0.0) { // advance position and momentum (focusing quad) - p.pos(RealAoS::x) = cos(omega*slice_ds)*x + sin(omega*slice_ds)/omega*px; + xout = cos(omega*slice_ds)*x + sin(omega*slice_ds)/omega*px; pxout = -omega*sin(omega*slice_ds)*x + cos(omega*slice_ds)*px; - p.pos(RealAoS::y) = cosh(omega*slice_ds)*y + sinh(omega*slice_ds)/omega*py; + yout = cosh(omega*slice_ds)*y + sinh(omega*slice_ds)/omega*py; pyout = omega*sinh(omega*slice_ds)*y + cosh(omega*slice_ds)*py; - p.pos(RealAoS::t) = t + (slice_ds/betgam2)*pt; + tout = t + (slice_ds/betgam2)*pt; // ptout = pt; } else { // advance position and momentum (defocusing quad) - p.pos(RealAoS::x) = cosh(omega*slice_ds)*x + sinh(omega*slice_ds)/omega*px; + xout = cosh(omega*slice_ds)*x + sinh(omega*slice_ds)/omega*px; pxout = omega*sinh(omega*slice_ds)*x + cosh(omega*slice_ds)*px; - p.pos(RealAoS::y) = cos(omega*slice_ds)*y + sin(omega*slice_ds)/omega*py; + yout = cos(omega*slice_ds)*y + sin(omega*slice_ds)/omega*py; pyout = -omega*sin(omega*slice_ds)*y + cos(omega*slice_ds)*py; - p.pos(RealAoS::t) = t + (slice_ds/betgam2)*pt; + tout = t + (slice_ds/betgam2)*pt; // ptout = pt; } - // assign updated momenta + // assign updated values + x = xout; + y = yout; + t = tout; px = pxout; py = pyout; pt = ptout; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. diff --git a/src/particles/elements/RFCavity.H b/src/particles/elements/RFCavity.H index b4e4025c8..5c955239d 100644 --- a/src/particles/elements/RFCavity.H +++ b/src/particles/elements/RFCavity.H @@ -177,30 +177,36 @@ namespace RFCavityData /** This is an RF cavity functor, so that a variable of this type can be used like * an RF cavity function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, amrex::ParticleReal & AMREX_RESTRICT px, amrex::ParticleReal & AMREX_RESTRICT py, amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, [[maybe_unused]] RefPart const & refpart ) const { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_in(x, y, px, py); - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + // intialize output values + amrex::ParticleReal xout = x; + amrex::ParticleReal yout = y; + amrex::ParticleReal tout = t; // initialize output values of momenta amrex::ParticleReal pxout = px; @@ -217,26 +223,29 @@ namespace RFCavityData // so that, e.g., R(3,4) = dyf/dpyi. // push particles using the linear map - p.pos(RealAoS::x) = R(1,1)*x + R(1,2)*px + R(1,3)*y + xout = R(1,1)*x + R(1,2)*px + R(1,3)*y + R(1,4)*py + R(1,5)*t + R(1,6)*pt; pxout = R(2,1)*x + R(2,2)*px + R(2,3)*y + R(2,4)*py + R(2,5)*t + R(2,6)*pt; - p.pos(RealAoS::y) = R(3,1)*x + R(3,2)*px + R(3,3)*y + yout = R(3,1)*x + R(3,2)*px + R(3,3)*y + R(3,4)*py + R(3,5)*t + R(3,6)*pt; pyout = R(4,1)*x + R(4,2)*px + R(4,3)*y + R(4,4)*py + R(4,5)*t + R(4,6)*pt; - p.pos(RealAoS::t) = R(5,1)*x + R(5,2)*px + R(5,3)*y + tout = R(5,1)*x + R(5,2)*px + R(5,3)*y + R(5,4)*py + R(5,5)*t + R(5,6)*pt; ptout = R(6,1)*x + R(6,2)*px + R(6,3)*y + R(6,4)*py + R(6,5)*t + R(6,6)*pt; - // assign updated momenta + // assign updated values + x = xout; + y = yout; + t = tout; px = pxout; py = pyout; pt = ptout; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. diff --git a/src/particles/elements/Sbend.H b/src/particles/elements/Sbend.H index f30a4723f..eb37fb0e3 100644 --- a/src/particles/elements/Sbend.H +++ b/src/particles/elements/Sbend.H @@ -62,29 +62,36 @@ namespace impactx /** This is a sbend functor, so that a variable of this type can be used like a sbend function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, - amrex::ParticleReal & AMREX_RESTRICT px, - amrex::ParticleReal & AMREX_RESTRICT py, - amrex::ParticleReal & AMREX_RESTRICT pt, - RefPart const & refpart) const { - + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, + amrex::ParticleReal & AMREX_RESTRICT px, + amrex::ParticleReal & AMREX_RESTRICT py, + amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, + RefPart const & refpart + ) const + { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_in(x, y, px, py); - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + // intialize output values + amrex::ParticleReal xout = x; + amrex::ParticleReal yout = y; + amrex::ParticleReal tout = t; // initialize output values of momenta amrex::ParticleReal pxout = px; @@ -104,27 +111,30 @@ namespace impactx auto const [sin_theta, cos_theta] = amrex::Math::sincos(theta); // advance position and momentum (sector bend) - p.pos(RealAoS::x) = cos_theta*x + m_rc*sin_theta*px + xout = cos_theta*x + m_rc*sin_theta*px - (m_rc/bet)*(1.0_prt - cos_theta)*pt; pxout = -sin_theta/m_rc*x + cos_theta*px - sin_theta/bet*pt; - p.pos(RealAoS::y) = y + m_rc*theta*py; + yout = y + m_rc*theta*py; // pyout = py; - p.pos(RealAoS::t) = sin_theta/bet*x + m_rc/bet*(1.0_prt - cos_theta)*px + t + tout = sin_theta/bet*x + m_rc/bet*(1.0_prt - cos_theta)*px + t + m_rc*(-theta+sin_theta/(bet*bet))*pt; // ptout = pt; - // assign updated momenta + // assign updated values + x = xout; + y = yout; + t = tout; px = pxout; py = pyout; pt = ptout; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. diff --git a/src/particles/elements/ShortRF.H b/src/particles/elements/ShortRF.H index a8b8c81cf..3c46ae8e2 100644 --- a/src/particles/elements/ShortRF.H +++ b/src/particles/elements/ShortRF.H @@ -64,29 +64,32 @@ namespace impactx /** This is a shortrf functor, so that a variable of this type can be used like a * shortrf function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, - amrex::ParticleReal & AMREX_RESTRICT px, - amrex::ParticleReal & AMREX_RESTRICT py, - amrex::ParticleReal & AMREX_RESTRICT pt, - RefPart const & refpart) const { + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, + amrex::ParticleReal & AMREX_RESTRICT px, + amrex::ParticleReal & AMREX_RESTRICT py, + amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, + RefPart const & refpart + ) const + { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); - - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + shift_in(x, y, px, py); // Define parameters and intermediate constants using ablastr::constant::math::pi; @@ -105,22 +108,28 @@ namespace impactx py = py*bgi; pt = pt*bgi; - // intialize output values of momenta + // intialize output values + amrex::ParticleReal xout = x; + amrex::ParticleReal yout = y; + amrex::ParticleReal tout = t; amrex::ParticleReal pxout = px; amrex::ParticleReal pyout = py; amrex::ParticleReal ptout = pt; // advance position and momentum in dynamic units - p.pos(RealAoS::x) = x; + // xout = x; pxout = px; - p.pos(RealAoS::y) = y; + // yout = y; pyout = py; - p.pos(RealAoS::t) = t; + // tout = t; ptout = pt - m_V*cos(k*t + phi) + m_V*cos(phi); - // assign updated momenta + // assign updated values + x = xout; + y = yout; + t = tout; px = pxout; py = pyout; pt = ptout; @@ -131,7 +140,7 @@ namespace impactx pt = pt/bgf; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. */ @@ -142,8 +151,8 @@ namespace impactx * @param[in,out] refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - void operator() (RefPart & AMREX_RESTRICT refpart) const { - + void operator() (RefPart & AMREX_RESTRICT refpart) const + { using namespace amrex::literals; // for _rt and _prt // assign input reference particle values diff --git a/src/particles/elements/SoftQuad.H b/src/particles/elements/SoftQuad.H index 6edf73e8f..1bc8faf8d 100644 --- a/src/particles/elements/SoftQuad.H +++ b/src/particles/elements/SoftQuad.H @@ -182,30 +182,36 @@ namespace SoftQuadrupoleData /** This is a soft-edge quadrupole functor, so that a variable of this type can be used * like a soft-edge quadrupole function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, amrex::ParticleReal & AMREX_RESTRICT px, amrex::ParticleReal & AMREX_RESTRICT py, amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, [[maybe_unused]] RefPart const & refpart ) const { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_in(x, y, px, py); - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + // intialize output values + amrex::ParticleReal xout = x; + amrex::ParticleReal yout = y; + amrex::ParticleReal tout = t; // initialize output values of momenta amrex::ParticleReal pxout = px; @@ -222,26 +228,29 @@ namespace SoftQuadrupoleData // so that, e.g., R(3,4) = dyf/dpyi. // push particles using the linear map - p.pos(RealAoS::x) = R(1,1)*x + R(1,2)*px + R(1,3)*y + xout = R(1,1)*x + R(1,2)*px + R(1,3)*y + R(1,4)*py + R(1,5)*t + R(1,6)*pt; pxout = R(2,1)*x + R(2,2)*px + R(2,3)*y + R(2,4)*py + R(2,5)*t + R(2,6)*pt; - p.pos(RealAoS::y) = R(3,1)*x + R(3,2)*px + R(3,3)*y + yout = R(3,1)*x + R(3,2)*px + R(3,3)*y + R(3,4)*py + R(3,5)*t + R(3,6)*pt; pyout = R(4,1)*x + R(4,2)*px + R(4,3)*y + R(4,4)*py + R(4,5)*t + R(4,6)*pt; - p.pos(RealAoS::t) = R(5,1)*x + R(5,2)*px + R(5,3)*y + tout = R(5,1)*x + R(5,2)*px + R(5,3)*y + R(5,4)*py + R(5,5)*t + R(5,6)*pt; ptout = R(6,1)*x + R(6,2)*px + R(6,3)*y + R(6,4)*py + R(6,5)*t + R(6,6)*pt; - // assign updated momenta + // assign updated values + x = xout; + y = yout; + t = tout; px = pxout; py = pyout; pt = ptout; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. diff --git a/src/particles/elements/SoftSol.H b/src/particles/elements/SoftSol.H index d7c0d25db..b4e05ce2c 100644 --- a/src/particles/elements/SoftSol.H +++ b/src/particles/elements/SoftSol.H @@ -187,30 +187,36 @@ namespace SoftSolenoidData /** This is a soft-edge solenoid functor, so that a variable of this type can be used * like a soft-edge solenoid function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, amrex::ParticleReal & AMREX_RESTRICT px, amrex::ParticleReal & AMREX_RESTRICT py, amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, [[maybe_unused]] RefPart const & refpart ) const { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_in(x, y, px, py); - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + // intialize output values + amrex::ParticleReal xout = x; + amrex::ParticleReal yout = y; + amrex::ParticleReal tout = t; // initialize output values of momenta amrex::ParticleReal pxout = px; @@ -227,26 +233,29 @@ namespace SoftSolenoidData // so that, e.g., R(3,4) = dyf/dpyi. // push particles using the linear map - p.pos(RealAoS::x) = R(1,1)*x + R(1,2)*px + R(1,3)*y + xout = R(1,1)*x + R(1,2)*px + R(1,3)*y + R(1,4)*py + R(1,5)*t + R(1,6)*pt; pxout = R(2,1)*x + R(2,2)*px + R(2,3)*y + R(2,4)*py + R(2,5)*t + R(2,6)*pt; - p.pos(RealAoS::y) = R(3,1)*x + R(3,2)*px + R(3,3)*y + yout = R(3,1)*x + R(3,2)*px + R(3,3)*y + R(3,4)*py + R(3,5)*t + R(3,6)*pt; pyout = R(4,1)*x + R(4,2)*px + R(4,3)*y + R(4,4)*py + R(4,5)*t + R(4,6)*pt; - p.pos(RealAoS::t) = R(5,1)*x + R(5,2)*px + R(5,3)*y + tout = R(5,1)*x + R(5,2)*px + R(5,3)*y + R(5,4)*py + R(5,5)*t + R(5,6)*pt; ptout = R(6,1)*x + R(6,2)*px + R(6,3)*y + R(6,4)*py + R(6,5)*t + R(6,6)*pt; - // assign updated momenta + // assign updated values + x = xout; + y = yout; + t = tout; px = pxout; py = pyout; pt = ptout; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. diff --git a/src/particles/elements/Sol.H b/src/particles/elements/Sol.H index 334e223a7..7404cce73 100644 --- a/src/particles/elements/Sol.H +++ b/src/particles/elements/Sol.H @@ -63,30 +63,31 @@ namespace impactx /** This is a sol functor, so that a variable of this type can be used like a sol function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType & AMREX_RESTRICT p, + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, amrex::ParticleReal & AMREX_RESTRICT px, amrex::ParticleReal & AMREX_RESTRICT py, amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, RefPart const & refpart ) const { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); - - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + shift_in(x, y, px, py); // length of the current slice amrex::ParticleReal const slice_ds = m_ds / nslice(); @@ -97,7 +98,7 @@ namespace impactx // compute phase advance per unit length (in rad/m) and // rotation angle (in rad) - amrex::ParticleReal const alpha = m_ks/2.0_prt; + amrex::ParticleReal const alpha = m_ks / 2.0_prt; amrex::ParticleReal const theta = alpha*slice_ds; // intialize output values @@ -125,22 +126,22 @@ namespace impactx pt = ptout; // advance positions and momenta using map for rotation - p.pos(RealAoS::x) = cos_theta*xout + sin_theta*yout; + x = cos_theta*xout + sin_theta*yout; pxout = cos_theta*px + sin_theta*py; - p.pos(RealAoS::y) = -sin_theta*xout + cos_theta*yout; + y = -sin_theta*xout + cos_theta*yout; pyout = -sin_theta*px + cos_theta*py; - p.pos(RealAoS::t) = tout; + t = tout; ptout = pt; - // assign updated momenta + // assign updated values px = pxout; py = pyout; pt = ptout; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. diff --git a/src/particles/elements/ThinDipole.H b/src/particles/elements/ThinDipole.H index a7a06e716..4e145920d 100644 --- a/src/particles/elements/ThinDipole.H +++ b/src/particles/elements/ThinDipole.H @@ -64,29 +64,36 @@ namespace impactx /** This is a multipole functor, so that a variable of this type can be used like a * multipole function. * - * @param p Particle AoS data for positions and cpu/id + * @param x particle position in x + * @param y particle position in y + * @param t particle position in t * @param px particle momentum in x * @param py particle momentum in y * @param pt particle momentum in t + * @param idcpu particle global index (unused) * @param refpart reference particle */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& AMREX_RESTRICT p, - amrex::ParticleReal & AMREX_RESTRICT px, - amrex::ParticleReal & AMREX_RESTRICT py, - amrex::ParticleReal & AMREX_RESTRICT pt, - RefPart const & refpart) const { - + amrex::ParticleReal & AMREX_RESTRICT x, + amrex::ParticleReal & AMREX_RESTRICT y, + amrex::ParticleReal & AMREX_RESTRICT t, + amrex::ParticleReal & AMREX_RESTRICT px, + amrex::ParticleReal & AMREX_RESTRICT py, + amrex::ParticleReal & AMREX_RESTRICT pt, + [[maybe_unused]] uint64_t & AMREX_RESTRICT idcpu, + RefPart const & refpart + ) const + { using namespace amrex::literals; // for _rt and _prt // shift due to alignment errors of the element - shift_in(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_in(x, y, px, py); - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); + // access position data + amrex::ParticleReal const xout = x; + amrex::ParticleReal const yout = y; + amrex::ParticleReal const tout = t; // access reference particle to find relativistic beta amrex::ParticleReal const beta_ref = refpart.beta(); @@ -105,13 +112,13 @@ namespace impactx amrex::ParticleReal kx = 1.0_prt/m_rc; // advance position and momentum - p.pos(RealAoS::x) = x; - pxout = px - pow(kx,2)*ds*x + kx*ds*f; //eq (3.2b) + x = xout; + pxout = px - pow(kx,2) * ds * xout + kx * ds * f; //eq (3.2b) - p.pos(RealAoS::y) = y; + y = yout; pyout = py; - p.pos(RealAoS::t) = t + kx*x*ds*fprime; //eq (3.2e) + t = tout + kx * xout * ds * fprime; //eq (3.2e) ptout = pt; // assign updated momenta @@ -120,7 +127,7 @@ namespace impactx pt = ptout; // undo shift due to alignment errors of the element - shift_out(p.pos(RealAoS::x), p.pos(RealAoS::y), px, py); + shift_out(x, y, px, py); } /** This pushes the reference particle. */ diff --git a/src/particles/elements/diagnostics/openPMD.cpp b/src/particles/elements/diagnostics/openPMD.cpp index b325e94e9..dba138d10 100644 --- a/src/particles/elements/diagnostics/openPMD.cpp +++ b/src/particles/elements/diagnostics/openPMD.cpp @@ -12,8 +12,6 @@ #include "ImpactXVersion.H" #include "particles/ImpactXParticleContainer.H" -#include - #include #include #include @@ -252,15 +250,6 @@ namespace detail auto d_fl = io::Dataset(dtype_fl, {np}); auto d_ui = io::Dataset(dtype_ui, {np}); - // AoS: Real - { - std::vector real_aos_names = get_RealAoS_names(); - for (auto real_idx=0; real_idx < RealAoS::nattribs; real_idx++) { - auto const component_name = real_aos_names.at(real_idx); - getComponentRecord(component_name).resetDataset(d_fl); - } - } - // reference particle information beam.setAttribute( "beta_ref", ref_part.beta() ); beam.setAttribute( "gamma_ref", ref_part.gamma() ); @@ -278,6 +267,7 @@ namespace detail // openPMD coarse position: for global coordinates { + beam["positionOffset"]["x"].resetDataset(d_fl); beam["positionOffset"]["x"].makeConstant(ref_part.x); beam["positionOffset"]["y"].resetDataset(d_fl); @@ -286,7 +276,7 @@ namespace detail beam["positionOffset"]["t"].makeConstant(ref_part.t); } - // AoS: Int + // unique, global particle index beam["id"][scalar].resetDataset(d_ui); // SoA: Real @@ -368,9 +358,6 @@ namespace detail auto & offset = m_offset.at(currentLevel); // ... - // preparing access to particle data: AoS - auto& aos = pti.GetArrayOfStructs(); - // series & iteration auto series = std::any_cast(m_series); io::WriteIterations iterations = series.writeIterations(); @@ -391,36 +378,12 @@ namespace detail return detail::get_component_record(beam, std::move(comp_name)); }; - // AoS: position and particle ID + // SoA + auto const& soa = pti.GetStructOfArrays(); + // particle id arrays { - using vs = std::vector; - vs const positionComponents{"x", "y", "t"}; // TODO: generalize - for (auto currDim = 0; currDim < AMREX_SPACEDIM; currDim++) { - std::shared_ptr const curr( - new amrex::ParticleReal[numParticleOnTile], - [](amrex::ParticleReal const *p) { delete[] p; } - ); - for (auto i = 0; i < numParticleOnTile; i++) { - curr.get()[i] = aos[i].pos(currDim); - } - std::string const& positionComponent = positionComponents[currDim]; - beam["position"][positionComponent].storeChunk(curr, {offset}, - {numParticleOnTile64}); - } - - // save particle ID after converting it to a globally unique ID - std::shared_ptr const ids( - new uint64_t[numParticleOnTile], - [](uint64_t const *p) { delete[] p; } - ); - for (auto i = 0; i < numParticleOnTile; i++) { - ids.get()[i] = ablastr::particles::localIDtoGlobal(aos[i].id(), aos[i].cpu()); - } - beam["id"][scalar].storeChunk(ids, {offset}, {numParticleOnTile64}); + beam["id"][scalar].storeChunkRaw(soa.GetIdCPUData().data(), {offset}, {numParticleOnTile64}); } - - // SoA: everything else - auto const& soa = pti.GetStructOfArrays(); // SoA floating point (ParticleReal) properties { std::vector real_soa_names = get_RealSoA_names(soa.NumRealComps()); diff --git a/src/particles/elements/mixin/beamoptic.H b/src/particles/elements/mixin/beamoptic.H index e4908463a..fd58ac600 100644 --- a/src/particles/elements/mixin/beamoptic.H +++ b/src/particles/elements/mixin/beamoptic.H @@ -40,25 +40,34 @@ namespace detail struct PushSingleParticle { using PType = ImpactXParticleContainer::ParticleType; + using ParticleTileType = amrex::ParticleTile,RealSoA::nattribs, IntSoA::nattribs>; /** Constructor taking in pointers to particle data * * @param element the beamline element to push through - * @param aos_ptr the array-of-struct with position and ids + * @param part_x the array to the particle position (x) + * @param part_y the array to the particle position (y) + * @param part_t the array to the particle position (t) * @param part_px the array to the particle momentum (x) * @param part_py the array to the particle momentum (y) * @param part_pt the array to the particle momentum (t) + * @param part_idcpu the array to the particle global index * @param ref_part the struct containing the reference particle */ PushSingleParticle (T_Element element, - PType* AMREX_RESTRICT aos_ptr, + amrex::ParticleReal* AMREX_RESTRICT part_x, + amrex::ParticleReal* AMREX_RESTRICT part_y, + amrex::ParticleReal* AMREX_RESTRICT part_t, amrex::ParticleReal* AMREX_RESTRICT part_px, amrex::ParticleReal* AMREX_RESTRICT part_py, amrex::ParticleReal* AMREX_RESTRICT part_pt, + uint64_t* AMREX_RESTRICT part_idcpu, RefPart ref_part) - : m_element(std::move(element)), m_aos_ptr(aos_ptr), - m_part_px(part_px), m_part_py(part_py), m_part_pt(part_pt), - m_ref_part(ref_part) + : m_element(std::move(element)), + m_part_x(part_x), m_part_y(part_y), m_part_t(part_t), + m_part_px(part_px), m_part_py(part_py), m_part_pt(part_pt), + m_part_idcpu(part_idcpu), + m_ref_part(std::move(ref_part)) { } @@ -75,25 +84,30 @@ namespace detail void operator() (long i) const { - // access AoS data such as positions and cpu/id - PType& AMREX_RESTRICT p = m_aos_ptr[i]; - // access SoA Real data + // note: an optimizing compiler will eliminate loads of unused parameters + amrex::ParticleReal & AMREX_RESTRICT x = m_part_x[i]; + amrex::ParticleReal & AMREX_RESTRICT y = m_part_y[i]; + amrex::ParticleReal & AMREX_RESTRICT t = m_part_t[i]; amrex::ParticleReal & AMREX_RESTRICT px = m_part_px[i]; amrex::ParticleReal & AMREX_RESTRICT py = m_part_py[i]; amrex::ParticleReal & AMREX_RESTRICT pt = m_part_pt[i]; + uint64_t & AMREX_RESTRICT idcpu = m_part_idcpu[i]; // push through element - m_element(p, px, py, pt, m_ref_part); + m_element(x, y, t, px, py, pt, idcpu, m_ref_part); } private: T_Element const m_element; - PType* const AMREX_RESTRICT m_aos_ptr; + amrex::ParticleReal* const AMREX_RESTRICT m_part_x; + amrex::ParticleReal* const AMREX_RESTRICT m_part_y; + amrex::ParticleReal* const AMREX_RESTRICT m_part_t; amrex::ParticleReal* const AMREX_RESTRICT m_part_px; amrex::ParticleReal* const AMREX_RESTRICT m_part_py; amrex::ParticleReal* const AMREX_RESTRICT m_part_pt; + uint64_t* const AMREX_RESTRICT m_part_idcpu; RefPart const m_ref_part; }; @@ -107,19 +121,19 @@ namespace detail ) { const int np = pti.numParticles(); - // preparing access to particle data: AoS - using PType = ImpactXParticleContainer::ParticleType; - auto& aos = pti.GetArrayOfStructs(); - PType* AMREX_RESTRICT aos_ptr = aos().dataPtr(); - // preparing access to particle data: SoA of Reals auto& soa_real = pti.GetStructOfArrays().GetRealData(); + amrex::ParticleReal* const AMREX_RESTRICT part_x = soa_real[RealSoA::x].dataPtr(); + amrex::ParticleReal* const AMREX_RESTRICT part_y = soa_real[RealSoA::y].dataPtr(); + amrex::ParticleReal* const AMREX_RESTRICT part_t = soa_real[RealSoA::t].dataPtr(); amrex::ParticleReal* const AMREX_RESTRICT part_px = soa_real[RealSoA::px].dataPtr(); amrex::ParticleReal* const AMREX_RESTRICT part_py = soa_real[RealSoA::py].dataPtr(); amrex::ParticleReal* const AMREX_RESTRICT part_pt = soa_real[RealSoA::pt].dataPtr(); + uint64_t* const AMREX_RESTRICT part_idcpu = pti.GetStructOfArrays().GetIdCPUData().dataPtr(); + detail::PushSingleParticle const pushSingleParticle( - element, aos_ptr, part_px, part_py, part_pt, ref_part); + element, part_x, part_y, part_t, part_px, part_py, part_pt, part_idcpu, ref_part); // loop over beam particles in the box amrex::ParallelFor(np, pushSingleParticle); } diff --git a/src/particles/integrators/Integrators.H b/src/particles/integrators/Integrators.H index c0ead2075..06381a5fc 100644 --- a/src/particles/integrators/Integrators.H +++ b/src/particles/integrators/Integrators.H @@ -10,8 +10,6 @@ #ifndef IMPACTX_INTEGRATORS_H_ #define IMPACTX_INTEGRATORS_H_ -#include - #include // for AMREX_RESTRICT #include // for ParticleReal diff --git a/src/particles/spacecharge/GatherAndPush.cpp b/src/particles/spacecharge/GatherAndPush.cpp index 2fc38d16e..fe0b7454f 100644 --- a/src/particles/spacecharge/GatherAndPush.cpp +++ b/src/particles/spacecharge/GatherAndPush.cpp @@ -63,13 +63,11 @@ namespace impactx::spacecharge amrex::ParticleReal const dt = slice_ds / pc.GetRefParticle().beta() / c0_SI; - // preparing access to particle data: AoS - using PType = ImpactXParticleContainer::ParticleType; - auto const & aos = pti.GetArrayOfStructs(); - PType const * const AMREX_RESTRICT aos_ptr = aos().dataPtr(); - // preparing access to particle data: SoA of Reals auto& soa_real = pti.GetStructOfArrays().GetRealData(); + amrex::ParticleReal* const AMREX_RESTRICT part_x = soa_real[RealSoA::x].dataPtr(); + amrex::ParticleReal* const AMREX_RESTRICT part_y = soa_real[RealSoA::y].dataPtr(); + amrex::ParticleReal* const AMREX_RESTRICT part_z = soa_real[RealSoA::z].dataPtr(); // note: currently for a fixed t amrex::ParticleReal* const AMREX_RESTRICT part_px = soa_real[RealSoA::px].dataPtr(); amrex::ParticleReal* const AMREX_RESTRICT part_py = soa_real[RealSoA::py].dataPtr(); amrex::ParticleReal* const AMREX_RESTRICT part_pz = soa_real[RealSoA::pz].dataPtr(); // note: currently for a fixed t @@ -79,10 +77,10 @@ namespace impactx::spacecharge // gather to each particle and push momentum amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE (int i) { - // access AoS data such as positions and cpu/id - PType const & AMREX_RESTRICT p = aos_ptr[i]; - // access SoA Real data + amrex::ParticleReal & AMREX_RESTRICT x = part_x[i]; + amrex::ParticleReal & AMREX_RESTRICT y = part_y[i]; + amrex::ParticleReal & AMREX_RESTRICT z = part_z[i]; amrex::ParticleReal & AMREX_RESTRICT px = part_px[i]; amrex::ParticleReal & AMREX_RESTRICT py = part_py[i]; amrex::ParticleReal & AMREX_RESTRICT pz = part_pz[i]; @@ -90,7 +88,7 @@ namespace impactx::spacecharge // force gather amrex::GpuArray const field_interp = ablastr::particles::doGatherVectorFieldNodal ( - p.pos(RealAoS::x), p.pos(RealAoS::y), p.pos(RealAoS::t), + x, y, z, scf_arr_x, scf_arr_y, scf_arr_z, invdr, prob_lo); diff --git a/src/particles/transformation/CoordinateTransformation.cpp b/src/particles/transformation/CoordinateTransformation.cpp index b52e86afe..cbf105533 100644 --- a/src/particles/transformation/CoordinateTransformation.cpp +++ b/src/particles/transformation/CoordinateTransformation.cpp @@ -48,19 +48,17 @@ namespace impactx::transformation for (ParIt pti(pc, lev); pti.isValid(); ++pti) { const int np = pti.numParticles(); - // preparing access to particle data: AoS - using PType = ImpactXParticleContainer::ParticleType; - auto &aos = pti.GetArrayOfStructs(); - PType *AMREX_RESTRICT aos_ptr = aos().dataPtr(); - // preparing access to particle data: SoA of Reals auto &soa_real = pti.GetStructOfArrays().GetRealData(); + amrex::ParticleReal *const AMREX_RESTRICT part_x = soa_real[RealSoA::x].dataPtr(); + amrex::ParticleReal *const AMREX_RESTRICT part_y = soa_real[RealSoA::y].dataPtr(); amrex::ParticleReal *const AMREX_RESTRICT part_px = soa_real[RealSoA::px].dataPtr(); amrex::ParticleReal *const AMREX_RESTRICT part_py = soa_real[RealSoA::py].dataPtr(); if (direction == CoordSystem::s) { BL_PROFILE("impactx::transformation::CoordinateTransformation::to_fixed_s"); + amrex::ParticleReal *const AMREX_RESTRICT part_z = soa_real[RealSoA::z].dataPtr(); amrex::ParticleReal *const AMREX_RESTRICT part_pz = soa_real[RealSoA::pz].dataPtr(); // Design value of pz/mc = beta*gamma @@ -68,33 +66,34 @@ namespace impactx::transformation ToFixedS const to_s(pzd); amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(long i) { - // access AoS data such as positions and cpu/id - PType &p = aos_ptr[i]; - - // access SoA Real data + amrex::ParticleReal &x = part_x[i]; + amrex::ParticleReal &y = part_y[i]; + amrex::ParticleReal &z = part_z[i]; amrex::ParticleReal &px = part_px[i]; amrex::ParticleReal &py = part_py[i]; amrex::ParticleReal &pz = part_pz[i]; - to_s(p, px, py, pz); + to_s(x, y, z, px, py, pz); }); } else { BL_PROFILE("impactx::transformation::CoordinateTransformation::to_fixed_t"); + amrex::ParticleReal *const AMREX_RESTRICT part_t = soa_real[RealSoA::t].dataPtr(); amrex::ParticleReal *const AMREX_RESTRICT part_pt = soa_real[RealSoA::pt].dataPtr(); - amrex::ParticleReal const ptd = pd; // Design value of pt/mc2 = -gamma. + // Design value of pt/mc2 = -gamma. + amrex::ParticleReal const ptd = pd; + ToFixedT const to_t(ptd); amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(long i) { - // access AoS data such as positions and cpu/id - PType &p = aos_ptr[i]; - - // access SoA Real data + amrex::ParticleReal &x = part_x[i]; + amrex::ParticleReal &y = part_y[i]; + amrex::ParticleReal &t = part_t[i]; amrex::ParticleReal &px = part_px[i]; amrex::ParticleReal &py = part_py[i]; amrex::ParticleReal &pt = part_pt[i]; - to_t(p, px, py, pt); + to_t(x, y, t, px, py, pt); }); } } // end loop over all particle boxes diff --git a/src/particles/transformation/ToFixedS.H b/src/particles/transformation/ToFixedS.H index c69fa0576..10cb1b297 100644 --- a/src/particles/transformation/ToFixedS.H +++ b/src/particles/transformation/ToFixedS.H @@ -40,25 +40,24 @@ namespace impactx::transformation /** This is a t-to-s map, so that a variable of this type can be used like a * t-to-s function. * - * @param[inout] p Particle AoS data for positions and cpu/id + * @param[inout] x particle position in x + * @param[inout] y particle position in y + * @param[inout] z particle position in z (in), in t (out) * @param[inout] px particle momentum in x * @param[inout] py particle momentum in y * @param[inout] pz particle momentum in z (in), in t (out) */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& p, + amrex::ParticleReal & x, + amrex::ParticleReal & y, + amrex::ParticleReal & z, amrex::ParticleReal & px, amrex::ParticleReal & py, amrex::ParticleReal & pz) const { using namespace amrex::literals; - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const z = p.pos(RealAoS::z); - // compute value of reference ptd = -gamma amrex::ParticleReal const argd = 1.0_prt + pow(m_pzd, 2); AMREX_ASSERT_WITH_MESSAGE(argd > 0.0_prt, "invalid ptd arg (<=0)"); @@ -76,11 +75,12 @@ namespace impactx::transformation amrex::ParticleReal const ptf = arg > 0.0_prt ? -sqrt(arg) : -1.0_prt; // transform position and momentum (from fixed t to fixed s) - p.pos(RealAoS::x) = x - px * z / (m_pzd + pz); + x = x - px * z / (m_pzd + pz); // px = px; - p.pos(RealAoS::y) = y - py * z / (m_pzd + pz); + y = y - py * z / (m_pzd + pz); // py = py; - p.pos(RealAoS::t) = ptf * z / (m_pzd + pz); + auto & t = z; // We store t in the same memory slot as z. + t = ptf * z / (m_pzd + pz); auto & pt = pz; // We store pt in the same memory slot as pz. pt = ptf - ptdf; diff --git a/src/particles/transformation/ToFixedT.H b/src/particles/transformation/ToFixedT.H index 669a3264d..577436ea3 100644 --- a/src/particles/transformation/ToFixedT.H +++ b/src/particles/transformation/ToFixedT.H @@ -40,25 +40,24 @@ namespace impactx::transformation /** This is a s-to-t map, so that a variable of this type can be used like a * s-to-t function. * - * @param[inout] p Particle AoS data for positions and cpu/id + * @param[inout] x particle position in x + * @param[inout] y particle position in y + * @param[inout] t particle position in t (in), in z (out) * @param[inout] px particle momentum in x * @param[inout] py particle momentum in y * @param[inout] pt particle momentum in t (in), in z (out) */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() ( - PType& p, + amrex::ParticleReal & x, + amrex::ParticleReal & y, + amrex::ParticleReal & t, amrex::ParticleReal & px, amrex::ParticleReal & py, amrex::ParticleReal & pt) const { using namespace amrex::literals; - // access AoS data such as positions and cpu/id - amrex::ParticleReal const x = p.pos(RealAoS::x); - amrex::ParticleReal const y = p.pos(RealAoS::y); - amrex::ParticleReal const t = p.pos(RealAoS::t); - // small tolerance to avoid NaN for pz<0: constexpr amrex::ParticleReal tol = 1.0e-8_prt; @@ -79,11 +78,12 @@ namespace impactx::transformation amrex::ParticleReal const pzf = arg > 0.0_prt ? sqrt(arg) : tol; // transform position and momentum (from fixed s to fixed t) - p.pos(RealAoS::x) = x + px*t/(m_ptd+pt); + x = x + px*t/(m_ptd+pt); // px = px; - p.pos(RealAoS::y) = y + py*t/(m_ptd+pt); + y = y + py*t/(m_ptd+pt); // py = py; - p.pos(RealAoS::z) = pzf * t / (m_ptd + pt); + auto & z = t; // We store z in the same memory slot as t. + z = pzf * t / (m_ptd + pt); auto & pz = pt; // We store pz in the same memory slot as pt. pz = pzf - pzdf; diff --git a/src/python/ImpactXParticleContainer.cpp b/src/python/ImpactXParticleContainer.cpp index 3b6fc3f57..93ef54fad 100644 --- a/src/python/ImpactXParticleContainer.cpp +++ b/src/python/ImpactXParticleContainer.cpp @@ -28,35 +28,31 @@ void init_impactxparticlecontainer(py::module& m) .export_values(); py::class_< - ParIter, - amrex::ParIter<0, 0, RealSoA::nattribs, IntSoA::nattribs> + ParIterSoA, + amrex::ParIterSoA >(m, "ImpactXParIter") - .def(py::init(), + .def(py::init(), py::arg("particle_container"), py::arg("level")) - .def(py::init(), + .def(py::init(), py::arg("particle_container"), py::arg("level"), py::arg("info")) ; py::class_< - ParConstIter, - amrex::ParConstIter<0, 0, RealSoA::nattribs, IntSoA::nattribs> + ParConstIterSoA, + amrex::ParConstIterSoA >(m, "ImpactXParConstIter") - .def(py::init(), + .def(py::init(), py::arg("particle_container"), py::arg("level")) - .def(py::init(), + .def(py::init(), py::arg("particle_container"), py::arg("level"), py::arg("info")) ; py::class_< ImpactXParticleContainer, - amrex::ParticleContainer<0, 0, RealSoA::nattribs, IntSoA::nattribs> + amrex::ParticleContainerPureSoA >(m, "ImpactXParticleContainer") //.def(py::init<>()) - .def_property_readonly_static("RealAoS", - [](py::object /* pc */){ return py::type::of(); }, - "RealAoS attribute name labels" - ) .def_property_readonly_static("RealSoA", [](py::object /* pc */){ return py::type::of(); }, "RealSoA attribute name labels" @@ -125,16 +121,10 @@ void init_impactxparticlecontainer(py::module& m) ) */ - .def_property_readonly("RealAoS_names", &ImpactXParticleContainer::RealAoS_names, - "Get the name of each Real AoS component") - .def_property_readonly("RealSoA_names", &ImpactXParticleContainer::RealSoA_names, "Get the name of each Real SoA component") ; - m.def("get_RealAoS_names", &get_RealAoS_names, - "Get the name of each Real AoS component"); - m.def("get_RealSoA_names", &get_RealSoA_names, py::arg("num_real_comps"), "Get the name of each Real SoA component\n\nnum_real_comps: pass number of compile-time + runtime arrays"); diff --git a/src/python/impactx/ImpactXParticleContainer.py b/src/python/impactx/ImpactXParticleContainer.py index 8547348db..82648cd94 100644 --- a/src/python/impactx/ImpactXParticleContainer.py +++ b/src/python/impactx/ImpactXParticleContainer.py @@ -36,9 +36,7 @@ def ix_pc_to_df(self, local=True, comm=None, root_rank=0): # todo: check if currently in fixed s or fixed t and pick name accordingly names = [] - for n in self.RealAoS_names: - names.append(n) - names.append("cpuid") + for n in self.RealSoA_names: names.append(n)