From 94ae11900131846e9f8b3704194673f7f02d8959 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Sun, 28 Jan 2024 10:37:19 -0800 Subject: [PATCH 1/9] Update Particle Container to Pure SoA (#3850) * Update Particle Container to Pure SoA Transition particle containers to pure SoA layouts. * Python: Pure SoA Particle * Update Particle Container to Pure SoA Transition particle containers to pure SoA layouts. * Python: Pure SoA Particle * fix a couple of places where we underflow the particle idcpu * fixed bad merge * fix bad merge again * Update Particle Container to Pure SoA Transition particle containers to pure SoA layouts. * Python: Pure SoA Particle * skip positions for pure SoA plotfiles * Update Particle Container to Pure SoA Transition particle containers to pure SoA layouts. Co-authored-by: Andrew Myers * Python: Pure SoA Particle * correctly account for positions in restart --------- Co-authored-by: Andrew Myers Co-authored-by: Andrew Myers --- Docs/source/usage/workflows/python_extend.rst | 3 +- .../particle_data_python/PICMI_inputs_2d.py | 6 +- .../PICMI_inputs_prev_pos_2d.py | 6 +- .../PICMI_inputs_runtime_component_analyze.py | 6 +- Python/pywarpx/particle_containers.py | 255 ++++++++---------- .../BackTransformParticleFunctor.H | 16 +- .../FlushFormats/FlushFormatAscent.cpp | 3 - .../FlushFormats/FlushFormatCheckpoint.cpp | 9 +- .../FlushFormats/FlushFormatPlotfile.cpp | 13 +- Source/Diagnostics/ParticleIO.cpp | 2 +- .../Diagnostics/ReducedDiags/FieldProbe.cpp | 6 +- .../FieldProbeParticleContainer.H | 20 +- .../FieldProbeParticleContainer.cpp | 36 +-- .../ReducedDiags/LoadBalanceCosts.cpp | 3 +- Source/Diagnostics/WarpXOpenPMD.H | 4 +- Source/Diagnostics/WarpXOpenPMD.cpp | 139 ++-------- .../ParticleBoundaryProcess.H | 5 +- Source/EmbeddedBoundary/ParticleScraper.H | 3 +- .../BinaryCollision/BinaryCollision.H | 16 +- .../Coulomb/PairWiseCoulombCollisionFunc.H | 7 +- .../Collision/BinaryCollision/DSMC/DSMC.H | 3 +- .../DSMC/SplitAndScatterFunc.H | 30 ++- .../NuclearFusion/NuclearFusionFunc.H | 14 +- .../ProtonBoronFusionInitializeMomentum.H | 10 +- .../TwoProductFusionInitializeMomentum.H | 4 +- .../BinaryCollision/ParticleCreationFunc.H | 60 +++-- .../BinaryCollision/ShuffleFisherYates.H | 2 +- .../Particles/Deposition/ChargeDeposition.H | 2 +- .../Particles/Deposition/CurrentDeposition.H | 2 +- .../ElementaryProcess/QEDPairGeneration.H | 2 +- .../ElementaryProcess/QEDPhotonEmission.H | 16 +- Source/Particles/LaserParticleContainer.H | 7 +- .../NamedComponentParticleContainer.H | 47 +++- Source/Particles/ParticleBoundaryBuffer.cpp | 4 +- Source/Particles/ParticleCreation/SmartCopy.H | 5 +- .../Particles/ParticleCreation/SmartCreate.H | 22 +- .../Particles/ParticleCreation/SmartUtils.H | 9 +- Source/Particles/PhysicalParticleContainer.H | 8 +- .../Particles/PhysicalParticleContainer.cpp | 129 ++++----- Source/Particles/Pusher/GetAndSetPosition.H | 152 +++++++---- .../Particles/Resampling/LevelingThinning.cpp | 5 +- Source/Particles/Sorting/Partition.cpp | 4 +- Source/Particles/Sorting/SortingUtils.H | 41 ++- Source/Particles/Sorting/SortingUtils.cpp | 2 +- Source/Particles/WarpXParticleContainer.H | 14 +- Source/Particles/WarpXParticleContainer.cpp | 89 +++--- .../Particles/ParticleBoundaryBuffer.cpp | 10 +- .../PinnedMemoryParticleContainer.cpp | 2 +- .../Particles/WarpXParticleContainer.cpp | 18 +- Source/Utils/ParticleUtils.H | 7 +- Source/Utils/ParticleUtils.cpp | 26 +- Source/ablastr/particles/IndexHandling.H | 41 --- Source/ablastr/particles/ParticleMoments.H | 25 +- 53 files changed, 661 insertions(+), 709 deletions(-) delete mode 100644 Source/ablastr/particles/IndexHandling.H diff --git a/Docs/source/usage/workflows/python_extend.rst b/Docs/source/usage/workflows/python_extend.rst index 6c9286c02ce..7a3ef0e7af9 100644 --- a/Docs/source/usage/workflows/python_extend.rst +++ b/Docs/source/usage/workflows/python_extend.rst @@ -252,7 +252,8 @@ Particles can be added to the simulation at specific positions and with specific .. autoclass:: pywarpx.particle_containers.ParticleContainerWrapper :members: -The ``get_particle_structs()`` and ``get_particle_arrays()`` functions are called +The ``get_particle_real_arrays()``, ``get_particle_int_arrays()`` and +``get_particle_idcpu_arrays()`` functions are called by several utility functions of the form ``get_particle_{comp_name}`` where ``comp_name`` is one of ``x``, ``y``, ``z``, ``r``, ``theta``, ``id``, ``cpu``, ``weight``, ``ux``, ``uy`` or ``uz``. diff --git a/Examples/Tests/particle_data_python/PICMI_inputs_2d.py b/Examples/Tests/particle_data_python/PICMI_inputs_2d.py index a4b7d9e134e..572871b8ed5 100755 --- a/Examples/Tests/particle_data_python/PICMI_inputs_2d.py +++ b/Examples/Tests/particle_data_python/PICMI_inputs_2d.py @@ -153,10 +153,10 @@ def add_particles(): ########################## assert (elec_wrapper.nps == 270 / (2 - args.unique)) -assert (elec_wrapper.particle_container.get_comp_index('w') == 0) -assert (elec_wrapper.particle_container.get_comp_index('newPid') == 4) +assert (elec_wrapper.particle_container.get_comp_index('w') == 2) +assert (elec_wrapper.particle_container.get_comp_index('newPid') == 6) -new_pid_vals = elec_wrapper.get_particle_arrays('newPid', 0) +new_pid_vals = elec_wrapper.get_particle_real_arrays('newPid', 0) for vals in new_pid_vals: assert np.allclose(vals, 5) diff --git a/Examples/Tests/particle_data_python/PICMI_inputs_prev_pos_2d.py b/Examples/Tests/particle_data_python/PICMI_inputs_prev_pos_2d.py index 1becd4464e7..5de9879f0f8 100755 --- a/Examples/Tests/particle_data_python/PICMI_inputs_prev_pos_2d.py +++ b/Examples/Tests/particle_data_python/PICMI_inputs_prev_pos_2d.py @@ -120,12 +120,12 @@ elec_count = elec_wrapper.nps # check that the runtime attributes have the right indices -assert (elec_wrapper.particle_container.get_comp_index('prev_x') == 4) -assert (elec_wrapper.particle_container.get_comp_index('prev_z') == 5) +assert (elec_wrapper.particle_container.get_comp_index('prev_x') == 6) +assert (elec_wrapper.particle_container.get_comp_index('prev_z') == 7) # sanity check that the prev_z values are reasonable and # that the correct number of values are returned -prev_z_vals = elec_wrapper.get_particle_arrays('prev_z', 0) +prev_z_vals = elec_wrapper.get_particle_real_arrays('prev_z', 0) running_count = 0 for z_vals in prev_z_vals: diff --git a/Examples/Tests/restart/PICMI_inputs_runtime_component_analyze.py b/Examples/Tests/restart/PICMI_inputs_runtime_component_analyze.py index 706dedb6959..32c9f4e5808 100755 --- a/Examples/Tests/restart/PICMI_inputs_runtime_component_analyze.py +++ b/Examples/Tests/restart/PICMI_inputs_runtime_component_analyze.py @@ -158,10 +158,10 @@ def add_particles(): ########################## assert electron_wrapper.nps == 90 -assert electron_wrapper.particle_container.get_comp_index("w") == 0 -assert electron_wrapper.particle_container.get_comp_index("newPid") == 4 +assert electron_wrapper.particle_container.get_comp_index("w") == 2 +assert electron_wrapper.particle_container.get_comp_index("newPid") == 6 -new_pid_vals = electron_wrapper.get_particle_arrays("newPid", 0) +new_pid_vals = electron_wrapper.get_particle_real_arrays("newPid", 0) for vals in new_pid_vals: assert np.allclose(vals, 5) diff --git a/Python/pywarpx/particle_containers.py b/Python/pywarpx/particle_containers.py index 273a981f4bd..72a675ec43c 100644 --- a/Python/pywarpx/particle_containers.py +++ b/Python/pywarpx/particle_containers.py @@ -117,8 +117,10 @@ def add_particles(self, x=None, y=None, z=None, ux=None, uy=None, kwargs[key] = np.full(maxlen, val) # --- The number of built in attributes + # --- The positions + built_in_attrs = libwarpx.dim # --- The three velocities - built_in_attrs = 3 + built_in_attrs += 3 if libwarpx.geometry_dim == 'rz': # --- With RZ, there is also theta built_in_attrs += 1 @@ -188,28 +190,25 @@ def add_real_comp(self, pid_name, comm=True): self.particle_container.add_real_comp(pid_name, comm) - def get_particle_structs(self, level, copy_to_host=False): + def get_particle_real_arrays(self, comp_name, level, copy_to_host=False): ''' - This returns a list of numpy or cupy arrays containing the particle struct data - on each tile for this process. The particle data is represented as a structured - array and contains the particle 'x', 'y', 'z', and 'idcpu'. + This returns a list of numpy or cupy arrays containing the particle real array data + on each tile for this process. - Unless copy_to_host is specified, the data for the structs are - not copied, but share the underlying memory buffer with WarpX. The + Unless copy_to_host is specified, the data for the arrays are not + copied, but share the underlying memory buffer with WarpX. The arrays are fully writeable. - Note that cupy does not support structs: - https://github.com/cupy/cupy/issues/2031 - and will return arrays of binary blobs for the AoS (DP: ``"|V24"``). If copied - to host via copy_to_host, we correct for the right numpy AoS type. - Parameters ---------- - level : int + comp_name : str + The component of the array data that will be returned + + level : int The refinement level to reference (default=0) - copy_to_host : bool + copy_to_host : bool For GPU-enabled runs, one can either return the GPU arrays (the default) or force a device-to-host copy. @@ -217,30 +216,29 @@ def get_particle_structs(self, level, copy_to_host=False): ------- List of arrays - The requested particle struct data + The requested particle array data ''' - particle_data = [] + comp_idx = self.particle_container.get_comp_index(comp_name) + + data_array = [] for pti in libwarpx.libwarpx_so.WarpXParIter(self.particle_container, level): + soa = pti.soa() + idx = soa.GetRealData(comp_idx) if copy_to_host: - particle_data.append(pti.aos().to_numpy(copy=True)) + data_array.append(idx.to_numpy(copy=True)) else: - if libwarpx.amr.Config.have_gpu: - libwarpx.amr.Print( - "get_particle_structs: cupy does not yet support structs. " - "https://github.com/cupy/cupy/issues/2031" - "Did you mean copy_to_host=True?" - ) xp, cupy_status = load_cupy() if cupy_status is not None: libwarpx.amr.Print(cupy_status) - aos_arr = xp.array(pti.aos(), copy=False) # void blobs for cupy - particle_data.append(aos_arr) - return particle_data + + data_array.append(xp.array(idx, copy=False)) + + return data_array - def get_particle_arrays(self, comp_name, level, copy_to_host=False): + def get_particle_int_arrays(self, comp_name, level, copy_to_host=False): ''' - This returns a list of numpy or cupy arrays containing the particle array data + This returns a list of numpy or cupy arrays containing the particle int array data on each tile for this process. Unless copy_to_host is specified, the data for the arrays are not @@ -266,12 +264,52 @@ def get_particle_arrays(self, comp_name, level, copy_to_host=False): List of arrays The requested particle array data ''' - comp_idx = self.particle_container.get_comp_index(comp_name) + comp_idx = self.particle_container.get_icomp_index(comp_name) data_array = [] for pti in libwarpx.libwarpx_so.WarpXParIter(self.particle_container, level): soa = pti.soa() - idx = soa.GetRealData(comp_idx) + idx = soa.GetIntData(comp_idx) + if copy_to_host: + data_array.append(idx.to_numpy(copy=True)) + else: + xp, cupy_status = load_cupy() + if cupy_status is not None: + libwarpx.amr.Print(cupy_status) + + data_array.append(xp.array(idx, copy=False)) + + return data_array + + + def get_particle_idcpu_arrays(self, level, copy_to_host=False): + ''' + This returns a list of numpy or cupy arrays containing the particle idcpu data + on each tile for this process. + + Unless copy_to_host is specified, the data for the arrays are not + copied, but share the underlying memory buffer with WarpX. The + arrays are fully writeable. + + Parameters + ---------- + level : int + The refinement level to reference (default=0) + + copy_to_host : bool + For GPU-enabled runs, one can either return the GPU + arrays (the default) or force a device-to-host copy. + + Returns + ------- + + List of arrays + The requested particle array data + ''' + data_array = [] + for pti in libwarpx.libwarpx_so.WarpXParIter(self.particle_container, level): + soa = pti.soa() + idx = soa.GetIdCPUData() if copy_to_host: data_array.append(idx.to_numpy(copy=True)) else: @@ -284,6 +322,31 @@ def get_particle_arrays(self, comp_name, level, copy_to_host=False): return data_array + def get_particle_idcpu(self, level=0, copy_to_host=False): + ''' + Return a list of numpy or cupy arrays containing the particle 'idcpu' + numbers on each tile. + + Parameters + ---------- + + level : int + The refinement level to reference (default=0) + + copy_to_host : bool + For GPU-enabled runs, one can either return the GPU + arrays (the default) or force a device-to-host copy. + + Returns + ------- + + List of arrays + The requested particle idcpu + ''' + return self.get_particle_idcpu_arrays(level, copy_to_host=copy_to_host) + idcpu = property(get_particle_idcpu) + + def get_particle_id(self, level=0, copy_to_host=False): ''' Return a list of numpy or cupy arrays containing the particle 'id' @@ -305,8 +368,8 @@ def get_particle_id(self, level=0, copy_to_host=False): List of arrays The requested particle ids ''' - structs = self.get_particle_structs(level, copy_to_host) - return [libwarpx.amr.unpack_ids(struct['cpuid']) for struct in structs] + idcpu = self.get_particle_idcpu(level, copy_to_host) + return [libwarpx.amr.unpack_ids(tile) for tile in idcpu] def get_particle_cpu(self, level=0, copy_to_host=False): @@ -330,8 +393,8 @@ def get_particle_cpu(self, level=0, copy_to_host=False): List of arrays The requested particle cpus ''' - structs = self.get_particle_structs(level, copy_to_host) - return [libwarpx.amr.unpack_cpus(struct['cpuid']) for struct in structs] + idcpu = self.get_particle_idcpu(level, copy_to_host) + return [libwarpx.amr.unpack_cpus(tile) for tile in idcpu] def get_particle_x(self, level=0, copy_to_host=False): @@ -355,20 +418,7 @@ def get_particle_x(self, level=0, copy_to_host=False): List of arrays The requested particle x position ''' - if copy_to_host: - # Use the numpy version of cosine - xp = np - else: - xp, cupy_status = load_cupy() - - structs = self.get_particle_structs(level, copy_to_host) - if libwarpx.geometry_dim == '3d' or libwarpx.geometry_dim == '2d': - return [struct['x'] for struct in structs] - elif libwarpx.geometry_dim == 'rz': - theta = self.get_particle_theta(level, copy_to_host) - return [struct['x']*xp.cos(theta) for struct, theta in zip(structs, theta)] - elif libwarpx.geometry_dim == '1d': - raise Exception('get_particle_x: There is no x coordinate with 1D Cartesian') + return self.get_particle_real_arrays('x', level, copy_to_host=copy_to_host) xp = property(get_particle_x) @@ -393,20 +443,7 @@ def get_particle_y(self, level=0, copy_to_host=False): List of arrays The requested particle y position ''' - if copy_to_host: - # Use the numpy version of sine - xp = np - else: - xp, cupy_status = load_cupy() - - structs = self.get_particle_structs(level, copy_to_host) - if libwarpx.geometry_dim == '3d': - return [struct['y'] for struct in structs] - elif libwarpx.geometry_dim == 'rz': - theta = self.get_particle_theta(level, copy_to_host) - return [struct['x']*xp.sin(theta) for struct, theta in zip(structs, theta)] - elif libwarpx.geometry_dim == '1d' or libwarpx.geometry_dim == '2d': - raise Exception('get_particle_y: There is no y coordinate with 1D or 2D Cartesian') + return self.get_particle_real_arrays('y', level, copy_to_host=copy_to_host) yp = property(get_particle_y) @@ -433,11 +470,12 @@ def get_particle_r(self, level=0, copy_to_host=False): ''' xp, cupy_status = load_cupy() - structs = self.get_particle_structs(level, copy_to_host) if libwarpx.geometry_dim == 'rz': - return [struct['x'] for struct in structs] + return self.get_particle_x(level, copy_to_host) elif libwarpx.geometry_dim == '3d': - return [xp.sqrt(struct['x']**2 + struct['y']**2) for struct in structs] + x = self.get_particle_x(level, copy_to_host) + y = self.get_particle_y(level, copy_to_host) + return xp.sqrt(x**2 + y**2) elif libwarpx.geometry_dim == '2d' or libwarpx.geometry_dim == '1d': raise Exception('get_particle_r: There is no r coordinate with 1D or 2D Cartesian') rp = property(get_particle_r) @@ -467,10 +505,11 @@ def get_particle_theta(self, level=0, copy_to_host=False): xp, cupy_status = load_cupy() if libwarpx.geometry_dim == 'rz': - return self.get_particle_arrays('theta', level, copy_to_host) + return self.get_particle_real_arrays('theta', level, copy_to_host) elif libwarpx.geometry_dim == '3d': - structs = self.get_particle_structs(level, copy_to_host) - return [xp.arctan2(struct['y'], struct['x']) for struct in structs] + x = self.get_particle_x(level, copy_to_host) + y = self.get_particle_y(level, copy_to_host) + return xp.arctan2(y, x) elif libwarpx.geometry_dim == '2d' or libwarpx.geometry_dim == '1d': raise Exception('get_particle_theta: There is no theta coordinate with 1D or 2D Cartesian') thetap = property(get_particle_theta) @@ -497,13 +536,7 @@ def get_particle_z(self, level=0, copy_to_host=False): List of arrays The requested particle z position ''' - structs = self.get_particle_structs(level, copy_to_host) - if libwarpx.geometry_dim == '3d': - return [struct['z'] for struct in structs] - elif libwarpx.geometry_dim == 'rz' or libwarpx.geometry_dim == '2d': - return [struct['y'] for struct in structs] - elif libwarpx.geometry_dim == '1d': - return [struct['x'] for struct in structs] + return self.get_particle_real_arrays('z', level, copy_to_host=copy_to_host) zp = property(get_particle_z) @@ -528,7 +561,7 @@ def get_particle_weight(self, level=0, copy_to_host=False): List of arrays The requested particle weight ''' - return self.get_particle_arrays('w', level, copy_to_host=copy_to_host) + return self.get_particle_real_arrays('w', level, copy_to_host=copy_to_host) wp = property(get_particle_weight) @@ -553,7 +586,7 @@ def get_particle_ux(self, level=0, copy_to_host=False): List of arrays The requested particle x momentum ''' - return self.get_particle_arrays('ux', level, copy_to_host=copy_to_host) + return self.get_particle_real_arrays('ux', level, copy_to_host=copy_to_host) uxp = property(get_particle_ux) @@ -578,7 +611,7 @@ def get_particle_uy(self, level=0, copy_to_host=False): List of arrays The requested particle y momentum ''' - return self.get_particle_arrays('uy', level, copy_to_host=copy_to_host) + return self.get_particle_real_arrays('uy', level, copy_to_host=copy_to_host) uyp = property(get_particle_uy) @@ -604,7 +637,7 @@ def get_particle_uz(self, level=0, copy_to_host=False): The requested particle z momentum ''' - return self.get_particle_arrays('uz', level, copy_to_host=copy_to_host) + return self.get_particle_real_arrays('uz', level, copy_to_host=copy_to_host) uzp = property(get_particle_uz) @@ -720,70 +753,6 @@ def get_particle_boundary_buffer_size(self, species_name, boundary, local=False) ) - def get_particle_boundary_buffer_structs( - self, species_name, boundary, level, copy_to_host=False - ): - ''' - This returns a list of numpy or cupy arrays containing the particle struct data - for a species that has been scraped by a specific simulation boundary. The - particle data is represented as a structured array and contains the - particle 'x', 'y', 'z', and 'idcpu'. - - Unless copy_to_host is specified, the data for the structs are - not copied, but share the underlying memory buffer with WarpX. The - arrays are fully writeable. - - Note that cupy does not support structs: - https://github.com/cupy/cupy/issues/2031 - and will return arrays of binary blobs for the AoS (DP: ``"|V24"``). If copied - to host via copy_to_host, we correct for the right numpy AoS type. - - Parameters - ---------- - - species_name : str - The species name that the data will be returned for - - boundary : str - The boundary from which to get the scraped particle data in the - form x/y/z_hi/lo or eb. - - level : int - The refinement level to reference (default=0) - - copy_to_host : bool - For GPU-enabled runs, one can either return the GPU - arrays (the default) or force a device-to-host copy. - - Returns - ------- - - List of arrays - The requested particle struct data - ''' - particle_container = self.particle_buffer.get_particle_container( - species_name, self._get_boundary_number(boundary) - ) - - particle_data = [] - for pti in libwarpx.libwarpx_so.BoundaryBufferParIter(particle_container, level): - if copy_to_host: - particle_data.append(pti.aos().to_numpy(copy=True)) - else: - if libwarpx.amr.Config.have_gpu: - libwarpx.amr.Print( - "get_particle_structs: cupy does not yet support structs. " - "https://github.com/cupy/cupy/issues/2031" - "Did you mean copy_to_host=True?" - ) - xp, cupy_status = load_cupy() - if cupy_status is not None: - libwarpx.amr.Print(cupy_status) - aos_arr = xp.array(pti.aos(), copy=False) # void blobs for cupy - particle_data.append(aos_arr) - return particle_data - - def get_particle_boundary_buffer(self, species_name, boundary, comp_name, level): ''' This returns a list of numpy or cupy arrays containing the particle array data diff --git a/Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.H b/Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.H index e88b522f148..fa00f6288f9 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.H +++ b/Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.H @@ -143,19 +143,19 @@ struct LorentzTransformParticles const amrex::ParticleReal uzp = uz_old_p * weight_old + uz_new_p * weight_new; #if defined (WARPX_DIM_3D) - dst.m_aos[i_dst].pos(0) = xp; - dst.m_aos[i_dst].pos(1) = yp; - dst.m_aos[i_dst].pos(2) = zp; + dst.m_rdata[PIdx::x][i_dst] = xp; + dst.m_rdata[PIdx::y][i_dst] = yp; + dst.m_rdata[PIdx::z][i_dst] = zp; #elif defined (WARPX_DIM_RZ) - dst.m_aos[i_dst].pos(0) = std::sqrt(xp*xp + yp*yp); - dst.m_aos[i_dst].pos(1) = zp; + dst.m_rdata[PIdx::x][i_dst] = std::sqrt(xp*xp + yp*yp); + dst.m_rdata[PIdx::z][i_dst] = zp; dst.m_rdata[PIdx::theta][i_dst] = std::atan2(yp, xp); #elif defined (WARPX_DIM_XZ) - dst.m_aos[i_dst].pos(0) = xp; - dst.m_aos[i_dst].pos(1) = zp; + dst.m_rdata[PIdx::x][i_dst] = xp; + dst.m_rdata[PIdx::z][i_dst] = zp; amrex::ignore_unused(yp); #elif defined (WARPX_DIM_1D_Z) - dst.m_aos[i_dst].pos(0) = zp; + dst.m_rdata[PIdx::z][i_dst] = zp; amrex::ignore_unused(xp, yp); #else amrex::ignore_unused(xp, yp, zp); diff --git a/Source/Diagnostics/FlushFormats/FlushFormatAscent.cpp b/Source/Diagnostics/FlushFormats/FlushFormatAscent.cpp index 980047e3b46..c224bc4f871 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatAscent.cpp +++ b/Source/Diagnostics/FlushFormats/FlushFormatAscent.cpp @@ -94,9 +94,6 @@ FlushFormatAscent::WriteParticles(const amrex::Vector& particle_di // get names of real comps std::map real_comps_map = pc->getParticleComps(); - // WarpXParticleContainer compile-time extra AoS attributes (Real): 0 - // WarpXParticleContainer compile-time extra AoS attributes (int): 0 - // WarpXParticleContainer compile-time extra SoA attributes (Real): PIdx::nattribs // not an efficient search, but N is small... for(int j = 0; j < PIdx::nattribs; ++j) diff --git a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp index 5f59cd723da..00c28eead51 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp +++ b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp @@ -178,8 +178,8 @@ FlushFormatCheckpoint::CheckpointParticles ( Vector real_names; Vector int_names; + // note: positions skipped here, since we reconstruct a plotfile SoA from them real_names.push_back("weight"); - real_names.push_back("momentum_x"); real_names.push_back("momentum_y"); real_names.push_back("momentum_z"); @@ -189,9 +189,12 @@ FlushFormatCheckpoint::CheckpointParticles ( #endif // get the names of the real comps - real_names.resize(pc->NumRealComps()); + // note: skips the mandatory AMREX_SPACEDIM positions for pure SoA + real_names.resize(pc->NumRealComps() - AMREX_SPACEDIM); auto runtime_rnames = pc->getParticleRuntimeComps(); - for (auto const& x : runtime_rnames) { real_names[x.second+PIdx::nattribs] = x.first; } + for (auto const& x : runtime_rnames) { + real_names[x.second + PIdx::nattribs - AMREX_SPACEDIM] = x.first; + } // and the int comps int_names.resize(pc->NumIntComps()); diff --git a/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp b/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp index df73ed34c94..5ec067a6eac 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp +++ b/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp @@ -355,8 +355,8 @@ FlushFormatPlotfile::WriteParticles(const std::string& dir, Vector int_flags; Vector real_flags; + // note: positions skipped here, since we reconstruct a plotfile SoA from them real_names.push_back("weight"); - real_names.push_back("momentum_x"); real_names.push_back("momentum_y"); real_names.push_back("momentum_z"); @@ -366,14 +366,21 @@ FlushFormatPlotfile::WriteParticles(const std::string& dir, #endif // get the names of the real comps - real_names.resize(tmp.NumRealComps()); + + // note: skips the mandatory AMREX_SPACEDIM positions for pure SoA + real_names.resize(tmp.NumRealComps() - AMREX_SPACEDIM); auto runtime_rnames = tmp.getParticleRuntimeComps(); - for (auto const& x : runtime_rnames) { real_names[x.second+PIdx::nattribs] = x.first; } + for (auto const& x : runtime_rnames) { + real_names[x.second + PIdx::nattribs - AMREX_SPACEDIM] = x.first; + } // plot any "extra" fields by default real_flags = part_diag.m_plot_flags; real_flags.resize(tmp.NumRealComps(), 1); + // note: skip the mandatory AMREX_SPACEDIM positions for pure SoA + real_flags.erase(real_flags.begin(), real_flags.begin() + AMREX_SPACEDIM); + // and the names int_names.resize(tmp.NumIntComps()); auto runtime_inames = tmp.getParticleRuntimeiComps(); diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index 7ca5e6541d7..a8bb9303fe1 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -160,7 +160,7 @@ MultiParticleContainer::Restart (const std::string& dir) ); } - for (int j = PIdx::nattribs; j < nr; ++j) { + for (int j = PIdx::nattribs-AMREX_SPACEDIM; j < nr; ++j) { const auto& comp_name = real_comp_names[j]; auto current_comp_names = pc->getParticleComps(); auto search = current_comp_names.find(comp_name); diff --git a/Source/Diagnostics/ReducedDiags/FieldProbe.cpp b/Source/Diagnostics/ReducedDiags/FieldProbe.cpp index 9f45392bb0a..24ad0e64ea8 100644 --- a/Source/Diagnostics/ReducedDiags/FieldProbe.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldProbe.cpp @@ -431,8 +431,6 @@ void FieldProbe::ComputeDiags (int step) { const auto getPosition = GetParticlePosition(pti); auto setPosition = SetParticlePosition(pti); - const auto& aos = pti.GetArrayOfStructs(); - const auto* AMREX_RESTRICT m_structs = aos().dataPtr(); auto const np = pti.numParticles(); if (update_particles_moving_window) @@ -482,6 +480,8 @@ void FieldProbe::ComputeDiags (int step) ParticleReal* const AMREX_RESTRICT part_Bz = attribs[FieldProbePIdx::Bz].dataPtr(); ParticleReal* const AMREX_RESTRICT part_S = attribs[FieldProbePIdx::S].dataPtr(); + auto * const AMREX_RESTRICT idcpu = pti.GetStructOfArrays().GetIdCPUData().data(); + const auto &xyzmin = WarpX::LowerCorner(box, lev, 0._rt); const std::array &dx = WarpX::CellSize(lev); @@ -556,7 +556,7 @@ void FieldProbe::ComputeDiags (int step) amrex::ParticleReal xp, yp, zp; getPosition(ip, xp, yp, zp); long idx = ip*noutputs; - dvp[idx++] = m_structs[ip].id(); + dvp[idx++] = amrex::ParticleIDWrapper{idcpu[ip]}; // all particles created on IO cpu dvp[idx++] = xp; dvp[idx++] = yp; dvp[idx++] = zp; diff --git a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H index c85bf8fd541..7d59ade5dc6 100644 --- a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H +++ b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H @@ -24,7 +24,14 @@ struct FieldProbePIdx { enum { - Ex = 0, Ey, Ez, +#if !defined (WARPX_DIM_1D_Z) + x, +#endif +#if defined (WARPX_DIM_3D) + y, +#endif + z, + Ex, Ey, Ez, Bx, By, Bz, S, //!< the Poynting vector #ifdef WARPX_DIM_RZ @@ -40,9 +47,14 @@ struct FieldProbePIdx * nattribs tells the particle container to allot 7 SOA values. */ class FieldProbeParticleContainer - : public amrex::ParticleContainer<0, 0, FieldProbePIdx::nattribs> + : public amrex::ParticleContainerPureSoA { public: + static constexpr int NStructReal = 0; + static constexpr int NStructInt = 0; + static constexpr int NReal = FieldProbePIdx::nattribs; + static constexpr int NInt = 0; + FieldProbeParticleContainer (amrex::AmrCore* amr_core); ~FieldProbeParticleContainer() override = default; @@ -52,9 +64,9 @@ public: FieldProbeParticleContainer& operator= ( FieldProbeParticleContainer&& ) = default; //! amrex iterator for our number of attributes - using iterator = amrex::ParIter<0, 0, FieldProbePIdx::nattribs, 0>; + using iterator = amrex::ParIterSoA; //! amrex iterator for our number of attributes (read-only) - using const_iterator = amrex::ParConstIter<0, 0, FieldProbePIdx::nattribs, 0>; + using const_iterator = amrex::ParConstIterSoA; //! similar to WarpXParticleContainer::AddNParticles but does not include u(x,y,z) void AddNParticles (int lev, amrex::Vector const & x, amrex::Vector const & y, amrex::Vector const & z); diff --git a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp index 1fd741ddc47..18137fe9b2c 100644 --- a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp @@ -59,7 +59,7 @@ using namespace amrex; FieldProbeParticleContainer::FieldProbeParticleContainer (AmrCore* amr_core) - : ParticleContainer<0, 0, FieldProbePIdx::nattribs>(amr_core->GetParGDB()) + : ParticleContainerPureSoA(amr_core->GetParGDB()) { SetParticleSize(); } @@ -89,33 +89,17 @@ FieldProbeParticleContainer::AddNParticles (int lev, * is then coppied to the permament tile which is stored on the particle * (particle_tile). */ + using PinnedTile = typename ContainerLike::ParticleTileType; - using PinnedTile = ParticleTile, - NArrayReal, NArrayInt, - amrex::PinnedArenaAllocator>; PinnedTile pinned_tile; pinned_tile.define(NumRuntimeRealComps(), NumRuntimeIntComps()); for (int i = 0; i < np; i++) { - ParticleType p; - p.id() = ParticleType::NextID(); - p.cpu() = ParallelDescriptor::MyProc(); -#if defined(WARPX_DIM_3D) - p.pos(0) = x[i]; - p.pos(1) = y[i]; - p.pos(2) = z[i]; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::ignore_unused(y); - p.pos(0) = x[i]; - p.pos(1) = z[i]; -#elif defined(WARPX_DIM_1D_Z) - amrex::ignore_unused(x, y); - p.pos(0) = z[i]; -#endif - - // write position, cpu id, and particle id to particle - pinned_tile.push_back(p); + auto & idcpu_data = pinned_tile.GetStructOfArrays().GetIdCPUData(); + idcpu_data.push_back(0); + amrex::ParticleIDWrapper{idcpu_data.back()} = ParticleType::NextID(); + amrex::ParticleCPUWrapper(idcpu_data.back()) = ParallelDescriptor::MyProc(); } // write Real attributes (SoA) to particle initialized zero @@ -125,7 +109,13 @@ FieldProbeParticleContainer::AddNParticles (int lev, #ifdef WARPX_DIM_RZ pinned_tile.push_back_real(FieldProbePIdx::theta, np, 0.0); #endif - +#if !defined (WARPX_DIM_1D_Z) + pinned_tile.push_back_real(FieldProbePIdx::x, x); +#endif +#if defined (WARPX_DIM_3D) + pinned_tile.push_back_real(FieldProbePIdx::y, y); +#endif + pinned_tile.push_back_real(FieldProbePIdx::z, z); pinned_tile.push_back_real(FieldProbePIdx::Ex, np, 0.0); pinned_tile.push_back_real(FieldProbePIdx::Ey, np, 0.0); pinned_tile.push_back_real(FieldProbePIdx::Ez, np, 0.0); diff --git a/Source/Diagnostics/ReducedDiags/LoadBalanceCosts.cpp b/Source/Diagnostics/ReducedDiags/LoadBalanceCosts.cpp index 893b00a5f00..b4e07b51982 100644 --- a/Source/Diagnostics/ReducedDiags/LoadBalanceCosts.cpp +++ b/Source/Diagnostics/ReducedDiags/LoadBalanceCosts.cpp @@ -56,8 +56,7 @@ namespace auto const & plev = pc.GetParticles(lev); auto const & ptile = plev.at(box_index); - auto const & aos = ptile.GetArrayOfStructs(); - auto const np = aos.numParticles(); + auto const np = ptile.numParticles(); num_macro_particles += np; } diff --git a/Source/Diagnostics/WarpXOpenPMD.H b/Source/Diagnostics/WarpXOpenPMD.H index e3b7b893d0a..110f1ada649 100644 --- a/Source/Diagnostics/WarpXOpenPMD.H +++ b/Source/Diagnostics/WarpXOpenPMD.H @@ -41,7 +41,7 @@ class WarpXParticleCounter { public: using ParticleContainer = typename WarpXParticleContainer::ContainerLike; - using ParticleIter = typename amrex::ParIter<0, 0, PIdx::nattribs, 0, amrex::PinnedArenaAllocator>; + using ParticleIter = typename amrex::ParIterSoA; WarpXParticleCounter (ParticleContainer* pc); [[nodiscard]] unsigned long GetTotalNumParticles () const {return m_Total;} @@ -77,7 +77,7 @@ class WarpXOpenPMDPlot { public: using ParticleContainer = typename WarpXParticleContainer::ContainerLike; - using ParticleIter = typename amrex::ParConstIter<0, 0, PIdx::nattribs, 0, amrex::PinnedArenaAllocator>; + using ParticleIter = typename amrex::ParConstIterSoA; /** Initialize openPMD I/O routines * diff --git a/Source/Diagnostics/WarpXOpenPMD.cpp b/Source/Diagnostics/WarpXOpenPMD.cpp index 71d96a47927..d9b70a2e18d 100644 --- a/Source/Diagnostics/WarpXOpenPMD.cpp +++ b/Source/Diagnostics/WarpXOpenPMD.cpp @@ -18,11 +18,9 @@ #include "WarpX.H" #include "OpenPMDHelpFunction.H" -#include #include #include -#include #include #include #include @@ -548,6 +546,13 @@ for (unsigned i = 0, n = particle_diags.size(); i < n; ++i) { // see openPMD ED-PIC extension for namings // note: an underscore separates the record name from its component // for non-scalar records +#if !defined (WARPX_DIM_1D_Z) + real_names.push_back("position_x"); +#endif +#if defined (WARPX_DIM_3D) + real_names.push_back("position_y"); +#endif + real_names.push_back("position_z"); real_names.push_back("weighting"); real_names.push_back("momentum_x"); real_names.push_back("momentum_y"); @@ -732,77 +737,7 @@ WarpXOpenPMDPlot::DumpToFile (ParticleContainer* pc, contributed_particles = true; - // get position and particle ID from aos - // note: this implementation iterates the AoS 4x... - // if we flush late as we do now, we can also copy out the data in one go - const auto &aos = pti.GetArrayOfStructs(); // size = numParticlesOnTile - { - // Save positions -#if defined(WARPX_DIM_RZ) - { - const std::shared_ptr z( - new amrex::ParticleReal[numParticleOnTile], - [](amrex::ParticleReal const *p) { delete[] p; } - ); - for (auto i = 0; i < numParticleOnTile; i++) { - z.get()[i] = aos[i].pos(1); // {0: "r", 1: "z"} - } - std::string const positionComponent = "z"; - currSpecies["position"]["z"].storeChunk(z, {offset}, {numParticleOnTile64}); - } - - // reconstruct x and y from polar coordinates r, theta - auto const& soa = pti.GetStructOfArrays(); - amrex::ParticleReal const* theta = soa.GetRealData(PIdx::theta).dataPtr(); - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(theta != nullptr, "openPMD: invalid theta pointer."); - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(int(soa.GetRealData(PIdx::theta).size()) == numParticleOnTile, - "openPMD: theta and tile size do not match"); - { - const std::shared_ptr< amrex::ParticleReal > x( - new amrex::ParticleReal[numParticleOnTile], - [](amrex::ParticleReal const *p){ delete[] p; } - ); - const std::shared_ptr< amrex::ParticleReal > y( - new amrex::ParticleReal[numParticleOnTile], - [](amrex::ParticleReal const *p){ delete[] p; } - ); - for (auto i=0; i 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]; - currSpecies["position"][positionComponent].storeChunk(curr, {offset}, - {numParticleOnTile64}); - } -#endif - - // save particle ID after converting it to a globally unique ID - const std::shared_ptr ids( - new uint64_t[numParticleOnTile], - [](uint64_t const *p) { delete[] p; } - ); - for (auto i = 0; i < numParticleOnTile; i++) { - ids.get()[i] = ablastr::particles::localIDtoGlobal(static_cast(aos[i].id()), static_cast(aos[i].cpu())); - } - const auto *const scalar = openPMD::RecordComponent::SCALAR; - currSpecies["id"][scalar].storeChunk(ids, {offset}, {numParticleOnTile64}); - - } - // save "extra" particle properties in AoS and SoA + // save particle properties SaveRealProperty(pti, currSpecies, offset, @@ -903,10 +838,9 @@ WarpXOpenPMDPlot::SetupRealProperties (ParticleContainer const * pc, std::set< std::string > addedRecords; // add meta-data per record only once for (auto idx=0; idxNumRealComps(); idx++) { - auto ii = ParticleContainer::NStructReal + idx; // jump over extra AoS names - if (write_real_comp[ii]) { + if (write_real_comp[idx]) { // handle scalar and non-scalar records by name - const auto [record_name, component_name] = detail::name2openPMD(real_comp_names[ii]); + const auto [record_name, component_name] = detail::name2openPMD(real_comp_names[idx]); auto currRecord = currSpecies[record_name]; // meta data for ED-PIC extension @@ -927,10 +861,9 @@ WarpXOpenPMDPlot::SetupRealProperties (ParticleContainer const * pc, } } for (auto idx=0; idx( numParticleOnTile ); - auto const& aos = pti.GetArrayOfStructs(); // size = numParticlesOnTile + auto const numParticleOnTile64 = static_cast(numParticleOnTile); auto const& soa = pti.GetStructOfArrays(); - // first we concatenate the AoS into contiguous arrays - { - // note: WarpX does not yet use extra AoS Real attributes - for( auto idx=0; idx d( - new amrex::ParticleReal[numParticleOnTile], - [](amrex::ParticleReal const *p){ delete[] p; } - ); - - for( auto kk=0; kk AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - void operator() (const PData& ptd, int i, + void operator() (PData& ptd, int i, const amrex::RealVect& /*pos*/, const amrex::RealVect& /*normal*/, amrex::RandomEngine const& /*engine*/) const noexcept { - auto& p = ptd.m_aos[i]; - p.id() = -p.id(); + ptd.id(i) = -ptd.id(i); } }; } diff --git a/Source/EmbeddedBoundary/ParticleScraper.H b/Source/EmbeddedBoundary/ParticleScraper.H index d6196c35f44..312b3762a7e 100644 --- a/Source/EmbeddedBoundary/ParticleScraper.H +++ b/Source/EmbeddedBoundary/ParticleScraper.H @@ -170,13 +170,12 @@ scrapeParticles (PC& pc, const amrex::Vector& distance_t auto& tile = pti.GetParticleTile(); auto ptd = tile.getParticleTileData(); const auto np = tile.numParticles(); - amrex::Particle<0,0> * const particles = tile.GetArrayOfStructs()().data(); auto phi = (*distance_to_eb[lev])[pti].array(); // signed distance function amrex::ParallelForRNG( np, [=] AMREX_GPU_DEVICE (const int ip, amrex::RandomEngine const& engine) noexcept { // skip particles that are already flagged for removal - if (particles[ip].id() < 0) return; + if (ptd.id(ip) < 0) return; amrex::ParticleReal xp, yp, zp; getPosition(ip, xp, yp, zp); diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H index 5c90dab25e6..c69f07acdb2 100644 --- a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H +++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H @@ -72,7 +72,8 @@ class BinaryCollision final // Define shortcuts for frequently-used type names using ParticleType = WarpXParticleContainer::ParticleType; using ParticleTileType = WarpXParticleContainer::ParticleTileType; - using ParticleBins = amrex::DenseBins; + using ParticleTileDataType = ParticleTileType::ParticleTileDataType; + using ParticleBins = amrex::DenseBins; using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; using index_type = ParticleBins::index_type; @@ -261,9 +262,6 @@ public: const amrex::ParticleReal q1 = species_1.getCharge(); const amrex::ParticleReal m1 = species_1.getMass(); auto get_position_1 = GetParticlePosition(ptile_1, getpos_offset); - // Needed to access the particle id - ParticleType * AMREX_RESTRICT - particle_ptr_1 = ptile_1.GetArrayOfStructs()().data(); amrex::Geometry const& geom = WarpX::GetInstance().Geom(lev); #if defined WARPX_DIM_1D_Z @@ -371,7 +369,7 @@ public: soa_1, soa_1, product_species_vector, tile_products_data, - particle_ptr_1, particle_ptr_1, m1, m1, + m1, m1, products_mass, p_mask, products_np, copy_species1, copy_species2, p_pair_indices_1, p_pair_indices_2, @@ -403,9 +401,6 @@ public: const amrex::ParticleReal q1 = species_1.getCharge(); const amrex::ParticleReal m1 = species_1.getMass(); auto get_position_1 = GetParticlePosition(ptile_1, getpos_offset); - // Needed to access the particle id - ParticleType * AMREX_RESTRICT - particle_ptr_1 = ptile_1.GetArrayOfStructs()().data(); // - Species 2 const auto soa_2 = ptile_2.getParticleTileData(); index_type* AMREX_RESTRICT indices_2 = bins_2.permutationPtr(); @@ -413,9 +408,6 @@ public: const amrex::ParticleReal q2 = species_2.getCharge(); const amrex::ParticleReal m2 = species_2.getMass(); auto get_position_2 = GetParticlePosition(ptile_2, getpos_offset); - // Needed to access the particle id - ParticleType * AMREX_RESTRICT - particle_ptr_2 = ptile_2.GetArrayOfStructs()().data(); amrex::Geometry const& geom = WarpX::GetInstance().Geom(lev); #if defined WARPX_DIM_1D_Z @@ -535,7 +527,7 @@ public: soa_1, soa_2, product_species_vector, tile_products_data, - particle_ptr_1, particle_ptr_2, m1, m2, + m1, m2, products_mass, p_mask, products_np, copy_species1, copy_species2, p_pair_indices_1, p_pair_indices_2, diff --git a/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H b/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H index feb7acf81d3..cfdc36d3c50 100644 --- a/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H @@ -23,10 +23,13 @@ * \brief This functor performs pairwise Coulomb collision on a single cell by calling the function * ElasticCollisionPerez. It also reads and contains the Coulomb logarithm. */ -class PairWiseCoulombCollisionFunc{ +class PairWiseCoulombCollisionFunc +{ // Define shortcuts for frequently-used type names using ParticleType = WarpXParticleContainer::ParticleType; - using ParticleBins = amrex::DenseBins; + using ParticleTileType = WarpXParticleContainer::ParticleTileType; + using ParticleTileDataType = ParticleTileType::ParticleTileDataType; + using ParticleBins = amrex::DenseBins; using index_type = ParticleBins::index_type; using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMC.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMC.H index c1be307b811..ab01eba2c81 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMC.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMC.H @@ -38,7 +38,8 @@ class DSMC final // Define shortcuts for frequently-used type names using ParticleType = WarpXParticleContainer::ParticleType; using ParticleTileType = WarpXParticleContainer::ParticleTileType; - using ParticleBins = amrex::DenseBins; + using ParticleTileDataType = ParticleTileType::ParticleTileDataType; + using ParticleBins = amrex::DenseBins; using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; using index_type = ParticleBins::index_type; diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H index c1fb7ee7e38..5a7a4f3237a 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H @@ -10,6 +10,9 @@ #define SPLIT_AND_SCATTER_FUNC_H_ #include "Particles/Collision/ScatteringProcess.H" +#include "Particles/NamedComponentParticleContainer.H" + +#include /** * \brief Function that performs the particle scattering and injection due @@ -55,8 +58,6 @@ int splitScatteringParticles ( const auto ptile1_data = ptile1.getParticleTileData(); const auto ptile2_data = ptile2.getParticleTileData(); - const Long minus_one_long = -1; - ParallelForRNG(n_total_pairs, [=] AMREX_GPU_DEVICE (int i, RandomEngine const& engine) noexcept { @@ -70,20 +71,37 @@ int splitScatteringParticles ( // starting with the parent particles auto& w1 = ptile1_data.m_rdata[PIdx::w][p_pair_indices_1[i]]; auto& w2 = ptile2_data.m_rdata[PIdx::w][p_pair_indices_2[i]]; + uint64_t* AMREX_RESTRICT idcpu1 = ptile1_data.m_idcpu; + uint64_t* AMREX_RESTRICT idcpu2 = ptile2_data.m_idcpu; + + // Note: Particle::atomicSetID should also be provided as a standalone helper function in AMReX + // to replace the following lambda. + auto const atomicSetIdMinus = [] AMREX_GPU_DEVICE (uint64_t & idcpu) + { + constexpr amrex::Long minus_one_long = -1; + uint64_t tmp = 0; + amrex::ParticleIDWrapper wrapper(tmp); + wrapper = minus_one_long; +#if defined(AMREX_USE_OMP) +#pragma omp atomic write + idcpu = wrapper.m_idata; +#else + auto *old_ptr = reinterpret_cast(&idcpu); + amrex::Gpu::Atomic::Exch(old_ptr, (unsigned long long) wrapper.m_idata); +#endif + }; // Remove p_pair_reaction_weight[i] from the colliding particles' weights. // If the colliding particle weight decreases to zero, remove particle by // setting its id to -1. Gpu::Atomic::AddNoRet(&w1, -p_pair_reaction_weight[i]); if (w1 <= 0._prt) { - auto& p = ptile1_data.m_aos[p_pair_indices_1[i]]; - p.atomicSetID(minus_one_long); + atomicSetIdMinus(idcpu1[p_pair_indices_1[i]]); } Gpu::Atomic::AddNoRet(&w2, -p_pair_reaction_weight[i]); if (w2 <= 0._prt) { - auto& p = ptile2_data.m_aos[p_pair_indices_2[i]]; - p.atomicSetID(minus_one_long); + atomicSetIdMinus(idcpu2[p_pair_indices_2[i]]); } // Set the child particle properties appropriately diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H index 397536b67bf..b2a2112ca68 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H @@ -33,10 +33,13 @@ * creation functor. * This functor also reads and contains the fusion multiplier. */ -class NuclearFusionFunc{ +class NuclearFusionFunc +{ // Define shortcuts for frequently-used type names using ParticleType = WarpXParticleContainer::ParticleType; - using ParticleBins = amrex::DenseBins; + using ParticleTileType = WarpXParticleContainer::ParticleTileType; + using ParticleTileDataType = ParticleTileType::ParticleTileDataType; + using ParticleBins = amrex::DenseBins; using index_type = ParticleBins::index_type; using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; @@ -154,12 +157,13 @@ public: // other species and we need to decrease their weight accordingly. // c1 corresponds to the minimum number of times a particle of species 1 will be paired // with a particle of species 2. Same for c2. - const index_type c1 = amrex::max(NI2/NI1,1u); - const index_type c2 = amrex::max(NI1/NI2,1u); + // index_type(1): https://github.com/AMReX-Codes/amrex/pull/3684 + const index_type c1 = amrex::max(NI2/NI1, index_type(1)); + const index_type c2 = amrex::max(NI1/NI2, index_type(1)); // multiplier ratio to take into account unsampled pairs const auto multiplier_ratio = static_cast( - (m_isSameSpecies)?(2u*max_N - 1):(max_N)); + m_isSameSpecies ? 2*max_N - 1 : max_N); #if (defined WARPX_DIM_RZ) amrex::ParticleReal * const AMREX_RESTRICT theta1 = soa_1.m_rdata[PIdx::theta]; diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H index 7b29267ec32..0b51d6b4b61 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H @@ -22,10 +22,12 @@ namespace { // Define shortcuts for frequently-used type names - using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; - using ParticleType = WarpXParticleContainer::ParticleType; - using ParticleBins = amrex::DenseBins; - using index_type = ParticleBins::index_type; + using SoaData_type = typename WarpXParticleContainer::ParticleTileType::ParticleTileDataType; + using ParticleType = typename WarpXParticleContainer::ParticleType; + using ParticleTileType = typename WarpXParticleContainer::ParticleTileType; + using ParticleTileDataType = typename ParticleTileType::ParticleTileDataType; + using ParticleBins = amrex::DenseBins; + using index_type = typename ParticleBins::index_type; /** * \brief This function initializes the momentum of the alpha particles produced from diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionInitializeMomentum.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionInitializeMomentum.H index be3f5b2d957..52e9db8aa94 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionInitializeMomentum.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionInitializeMomentum.H @@ -24,7 +24,9 @@ namespace { // Define shortcuts for frequently-used type names using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; using ParticleType = WarpXParticleContainer::ParticleType; - using ParticleBins = amrex::DenseBins; + using ParticleTileType = WarpXParticleContainer::ParticleTileType; + using ParticleTileDataType = ParticleTileType::ParticleTileDataType; + using ParticleBins = amrex::DenseBins; using index_type = ParticleBins::index_type; /** diff --git a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H index dc830b477df..3928c74223d 100644 --- a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H +++ b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H @@ -30,13 +30,15 @@ * \brief This functor creates particles produced from a binary collision and sets their initial * properties (position, momentum, weight). */ -class ParticleCreationFunc{ +class ParticleCreationFunc +{ // Define shortcuts for frequently-used type names - using ParticleType = WarpXParticleContainer::ParticleType; - using ParticleTileType = WarpXParticleContainer::ParticleTileType; - using ParticleBins = amrex::DenseBins; - using index_type = ParticleBins::index_type; - using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; + using ParticleType = typename WarpXParticleContainer::ParticleType; + using ParticleTileType = typename WarpXParticleContainer::ParticleTileType; + using ParticleTileDataType = typename ParticleTileType::ParticleTileDataType; + using ParticleBins = amrex::DenseBins; + using index_type = typename ParticleBins::index_type; + using SoaData_type = typename WarpXParticleContainer::ParticleTileType::ParticleTileDataType; public: /** @@ -69,12 +71,6 @@ public: * @param[in, out] soa_1 struct of array data of the first colliding particle species * @param[in, out] soa_2 struct of array data of the second colliding particle species * @param[out] tile_products array containing tile data of the product particles. - * @param[out] particle_ptr_1 pointer to data of the first colliding particle species. Is - * needed to set the id of a particle to -1 in order to delete it when its weight - * reaches 0. - * @param[out] particle_ptr_2 pointer to data of the second colliding particle species. Is - * needed to set the id of a particle to -1 in order to delete it when its weight - * reaches 0. * @param[in] m1 mass of the first colliding particle species * @param[in] m2 mass of the second colliding particle species * @param[in] products_mass array storing the mass of product particles @@ -102,7 +98,6 @@ public: const SoaData_type& soa_1, const SoaData_type& soa_2, const amrex::Vector& pc_products, ParticleTileType** AMREX_RESTRICT tile_products, - ParticleType* particle_ptr_1, ParticleType* particle_ptr_2, const amrex::ParticleReal& m1, const amrex::ParticleReal& m2, const amrex::Vector& products_mass, const index_type* AMREX_RESTRICT p_mask, @@ -137,6 +132,8 @@ public: amrex::ParticleReal* AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w]; amrex::ParticleReal* AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w]; + uint64_t* AMREX_RESTRICT idcpu1 = soa_1.m_idcpu; + uint64_t* AMREX_RESTRICT idcpu2 = soa_2.m_idcpu; // Create necessary GPU vectors, that will be used in the kernel below amrex::Vector soa_products; @@ -205,16 +202,32 @@ public: amrex::Gpu::Atomic::AddNoRet(&w2[p_pair_indices_2[i]], -p_pair_reaction_weight[i]); + // Note: Particle::atomicSetID should also be provided as a standalone helper function in AMReX + // to replace the following lambda. + auto const atomicSetIdMinus = [] AMREX_GPU_DEVICE (uint64_t & idcpu) + { + constexpr amrex::Long minus_one_long = -1; + uint64_t tmp = 0; + amrex::ParticleIDWrapper wrapper(tmp); + wrapper = minus_one_long; +#if defined(AMREX_USE_OMP) +#pragma omp atomic write + idcpu = wrapper.m_idata; +#else + auto *old_ptr = reinterpret_cast(&idcpu); + amrex::Gpu::Atomic::Exch(old_ptr, (unsigned long long) wrapper.m_idata); +#endif + }; + // If the colliding particle weight decreases to zero, remove particle by // setting its id to -1 - constexpr amrex::Long minus_one_long = -1; if (w1[p_pair_indices_1[i]] <= amrex::ParticleReal(0.)) { - particle_ptr_1[p_pair_indices_1[i]].atomicSetID(minus_one_long); + atomicSetIdMinus(idcpu1[p_pair_indices_1[i]]); } if (w2[p_pair_indices_2[i]] <= amrex::ParticleReal(0.)) { - particle_ptr_2[p_pair_indices_2[i]].atomicSetID(minus_one_long); + atomicSetIdMinus(idcpu2[p_pair_indices_2[i]]); } // Initialize the product particles' momentum, using a function depending on the @@ -294,12 +307,14 @@ private: * \brief This class does nothing and is used as second template parameter for binary collisions * that do not create particles. */ -class NoParticleCreationFunc{ - using ParticleType = WarpXParticleContainer::ParticleType; - using ParticleTileType = WarpXParticleContainer::ParticleTileType; - using ParticleBins = amrex::DenseBins; - using index_type = ParticleBins::index_type; - using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; +class NoParticleCreationFunc +{ + using ParticleType = typename WarpXParticleContainer::ParticleType; + using ParticleTileType = typename WarpXParticleContainer::ParticleTileType; + using ParticleTileDataType = typename ParticleTileType::ParticleTileDataType; + using ParticleBins = amrex::DenseBins; + using index_type = typename ParticleBins::index_type; + using SoaData_type = typename WarpXParticleContainer::ParticleTileType::ParticleTileDataType; public: NoParticleCreationFunc () = default; @@ -313,7 +328,6 @@ public: const SoaData_type& /*soa_1*/, const SoaData_type& /*soa_2*/, amrex::Vector& /*pc_products*/, ParticleTileType** /*tile_products*/, - ParticleType* /*particle_ptr_1*/, ParticleType* /*particle_ptr_2*/, const amrex::ParticleReal& /*m1*/, const amrex::ParticleReal& /*m2*/, const amrex::Vector& /*products_mass*/, const index_type* /*p_mask*/, const amrex::Vector& /*products_np*/, diff --git a/Source/Particles/Collision/BinaryCollision/ShuffleFisherYates.H b/Source/Particles/Collision/BinaryCollision/ShuffleFisherYates.H index 42259512b0d..3b8f72f4b84 100644 --- a/Source/Particles/Collision/BinaryCollision/ShuffleFisherYates.H +++ b/Source/Particles/Collision/BinaryCollision/ShuffleFisherYates.H @@ -12,7 +12,7 @@ /* \brief Shuffle array according to Fisher-Yates algorithm. * Only shuffle the part between is <= i < ie, n = ie-is. * T_index shall be - * amrex::DenseBins::index_type + * amrex::DenseBins::index_type */ template diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index d0db678dfda..d0822789015 100644 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -252,7 +252,7 @@ void doChargeDepositionSharedShapeN (const GetParticlePosition& GetPositio const int n_rz_azimuthal_modes, amrex::Real* cost, const long load_balance_costs_update_algo, - const amrex::DenseBins& a_bins, + const amrex::DenseBins& a_bins, const amrex::Box& box, const amrex::Geometry& geom, const amrex::IntVect& a_tbox_max_size, diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index efe6efcc788..7343a8c5626 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -592,7 +592,7 @@ void doDepositionSharedShapeN (const GetParticlePosition& GetPosition, int n_rz_azimuthal_modes, amrex::Real* cost, long load_balance_costs_update_algo, - const amrex::DenseBins& a_bins, + const amrex::DenseBins& a_bins, const amrex::Box& box, const amrex::Geometry& geom, const amrex::IntVect& a_tbox_max_size) diff --git a/Source/Particles/ElementaryProcess/QEDPairGeneration.H b/Source/Particles/ElementaryProcess/QEDPairGeneration.H index 5abc9282d4f..ca00b56323a 100644 --- a/Source/Particles/ElementaryProcess/QEDPairGeneration.H +++ b/Source/Particles/ElementaryProcess/QEDPairGeneration.H @@ -167,7 +167,7 @@ public: p_ux, p_uy, p_uz, engine); - src.m_aos[i_src].id() = -1; //destroy photon after pair generation + amrex::ParticleIDWrapper{src.m_idcpu[i_src]} = -1; // destroy photon after pair generation } private: diff --git a/Source/Particles/ElementaryProcess/QEDPhotonEmission.H b/Source/Particles/ElementaryProcess/QEDPhotonEmission.H index 8ba5c63ad57..f3099bf997f 100644 --- a/Source/Particles/ElementaryProcess/QEDPhotonEmission.H +++ b/Source/Particles/ElementaryProcess/QEDPhotonEmission.H @@ -22,6 +22,7 @@ #include #include #include +#include #include #include @@ -237,12 +238,11 @@ void cleanLowEnergyPhotons( const int old_size, const int num_added, const amrex::ParticleReal energy_threshold) { - auto pp = ptile.GetArrayOfStructs()().data() + old_size; - - const auto& soa = ptile.GetStructOfArrays(); + auto& soa = ptile.GetStructOfArrays(); + auto p_idcpu = soa.GetIdCPUData().data() + old_size; const auto p_ux = soa.GetRealData(PIdx::ux).data() + old_size; - const auto p_uy = soa.GetRealData(PIdx::uy).data() + old_size; - const auto p_uz = soa.GetRealData(PIdx::uz).data() + old_size; + const auto p_uy = soa.GetRealData(PIdx::uy).data() + old_size; + const auto p_uz = soa.GetRealData(PIdx::uz).data() + old_size; //The square of the energy threshold const auto energy_threshold2 = std::max( @@ -251,8 +251,6 @@ void cleanLowEnergyPhotons( amrex::ParallelFor(num_added, [=] AMREX_GPU_DEVICE (int ip) noexcept { - auto& p = pp[ip]; - const auto ux = p_ux[ip]; const auto uy = p_uy[ip]; const auto uz = p_uz[ip]; @@ -262,8 +260,8 @@ void cleanLowEnergyPhotons( constexpr amrex::ParticleReal me_c = PhysConst::m_e*PhysConst::c; const auto phot_energy2 = (ux*ux + uy*uy + uz*uz)*me_c*me_c; - if (phot_energy2 < energy_threshold2){ - p.id() = - 1; + if (phot_energy2 < energy_threshold2) { + amrex::ParticleIDWrapper{p_idcpu[ip]} = -1; } }); } diff --git a/Source/Particles/LaserParticleContainer.H b/Source/Particles/LaserParticleContainer.H index e6fa308431c..fac94ff20a3 100644 --- a/Source/Particles/LaserParticleContainer.H +++ b/Source/Particles/LaserParticleContainer.H @@ -56,10 +56,9 @@ public: * \brief Method to initialize runtime attributes. Does nothing for LaserParticleContainer. */ void DefaultInitializeRuntimeAttributes ( - amrex::ParticleTile, - NArrayReal, NArrayInt, amrex::PinnedArenaAllocator>& /*pinned_tile*/, - const int /*n_external_attr_real*/, - const int /*n_external_attr_int*/) final {} + typename ContainerLike::ParticleTileType& /*pinned_tile*/, + int /*n_external_attr_real*/, + int /*n_external_attr_int*/) final {} void ReadHeader (std::istream& is) final; diff --git a/Source/Particles/NamedComponentParticleContainer.H b/Source/Particles/NamedComponentParticleContainer.H index 3be0886425d..e04e7dba6df 100644 --- a/Source/Particles/NamedComponentParticleContainer.H +++ b/Source/Particles/NamedComponentParticleContainer.H @@ -18,12 +18,19 @@ #include -/** Particle Attributes stored in amrex::ParticleContainer's struct of array +/** Real Particle Attributes stored in amrex::ParticleContainer's struct of array */ struct PIdx { enum { - w = 0, ///< weight +#if !defined (WARPX_DIM_1D_Z) + x, +#endif +#if defined (WARPX_DIM_3D) + y, +#endif + z, + w, ///< weight ux, uy, uz, #ifdef WARPX_DIM_RZ theta, ///< RZ needs all three position components @@ -32,6 +39,15 @@ struct PIdx }; }; +/** Integer Particle Attributes stored in amrex::ParticleContainer's struct of array + */ +struct PIdxInt +{ + enum { + nattribs ///< number of attributes + }; +}; + /** Particle Container class that allows to add/access particle components * with a name (string) instead of doing so with an integer index. * (The "components" are all the particle quantities - except those @@ -45,11 +61,11 @@ struct PIdx */ template class T_Allocator=amrex::DefaultAllocator> class NamedComponentParticleContainer : -public amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator> +public amrex::ParticleContainerPureSoA { public: /** Construct an empty NamedComponentParticleContainer **/ - NamedComponentParticleContainer () : amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>() {} + NamedComponentParticleContainer () : amrex::ParticleContainerPureSoA() {} /** Construct a NamedComponentParticleContainer from an AmrParGDB object * @@ -61,8 +77,15 @@ public: * AMR hierarchy. Usually, this is generated by an AmrCore or AmrLevel object. */ NamedComponentParticleContainer (amrex::AmrParGDB* amr_pgdb) - : amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>(amr_pgdb) { + : amrex::ParticleContainerPureSoA(amr_pgdb) { // build up the map of string names to particle component numbers +#if !defined (WARPX_DIM_1D_Z) + particle_comps["x"] = PIdx::x; +#endif +#if defined (WARPX_DIM_3D) + particle_comps["y"] = PIdx::y; +#endif + particle_comps["z"] = PIdx::z; particle_comps["w"] = PIdx::w; particle_comps["ux"] = PIdx::ux; particle_comps["uy"] = PIdx::uy; @@ -85,12 +108,12 @@ public: * @param p_ricomps name-to-index map for run-time integer components */ NamedComponentParticleContainer( - amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator> && pc, + amrex::ParticleContainerPureSoA && pc, std::map p_comps, std::map p_icomps, std::map p_rcomps, std::map p_ricomps) - : amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>(std::move(pc)), + : amrex::ParticleContainerPureSoA(std::move(pc)), particle_comps(std::move(p_comps)), particle_icomps(std::move(p_icomps)), particle_runtime_comps(std::move(p_rcomps)), @@ -118,7 +141,7 @@ public: NamedComponentParticleContainer make_alike () const { auto tmp = NamedComponentParticleContainer( - amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>::template make_alike(), + amrex::ParticleContainerPureSoA::template make_alike(), particle_comps, particle_icomps, particle_runtime_comps, @@ -127,10 +150,10 @@ public: return tmp; } - using amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>::NumRealComps; - using amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>::NumIntComps; - using amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>::AddRealComp; - using amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>::AddIntComp; + using amrex::ParticleContainerPureSoA::NumRealComps; + using amrex::ParticleContainerPureSoA::NumIntComps; + using amrex::ParticleContainerPureSoA::AddRealComp; + using amrex::ParticleContainerPureSoA::AddIntComp; /** Allocate a new run-time real component * diff --git a/Source/Particles/ParticleBoundaryBuffer.cpp b/Source/Particles/ParticleBoundaryBuffer.cpp index 54c4396379d..88304bd8a9c 100644 --- a/Source/Particles/ParticleBoundaryBuffer.cpp +++ b/Source/Particles/ParticleBoundaryBuffer.cpp @@ -50,7 +50,7 @@ struct CopyAndTimestamp { void operator() (const DstData& dst, const SrcData& src, int src_i, int dst_i) const noexcept { - dst.m_aos[dst_i] = src.m_aos[src_i]; + dst.m_idcpu[dst_i] = src.m_idcpu[src_i]; for (int j = 0; j < SrcData::NAR; ++j) { dst.m_rdata[j][dst_i] = src.m_rdata[j][src_i]; } @@ -222,7 +222,7 @@ void ParticleBoundaryBuffer::gatherParticles (MultiParticleContainer& mypc, { WARPX_PROFILE("ParticleBoundaryBuffer::gatherParticles"); - using PIter = amrex::ParConstIter<0,0,PIdx::nattribs>; + using PIter = amrex::ParConstIterSoA; const auto& warpx_instance = WarpX::GetInstance(); const amrex::Geometry& geom = warpx_instance.Geom(0); auto plo = geom.ProbLoArray(); diff --git a/Source/Particles/ParticleCreation/SmartCopy.H b/Source/Particles/ParticleCreation/SmartCopy.H index 2c04baa18bb..6a6ceb3d290 100644 --- a/Source/Particles/ParticleCreation/SmartCopy.H +++ b/Source/Particles/ParticleCreation/SmartCopy.H @@ -26,7 +26,7 @@ * type. Second, if a given component name is found in both the src * and the dst, then the src value is copied. * - * Particle structs - positions and id numbers - are always copied. + * Particle positions and id numbers are always copied. * * You don't create this directly - use the SmartCopyFactory object below. */ @@ -48,9 +48,6 @@ struct SmartCopy void operator() (DstData& dst, const SrcData& src, int i_src, int i_dst, amrex::RandomEngine const& engine) const noexcept { - // the particle struct is always copied over - dst.m_aos[i_dst] = src.m_aos[i_src]; - // initialize the real components for (int j = 0; j < DstData::NAR; ++j) { dst.m_rdata[j][i_dst] = initializeRealValue(m_policy_real[j], engine); diff --git a/Source/Particles/ParticleCreation/SmartCreate.H b/Source/Particles/ParticleCreation/SmartCreate.H index 67d7767a5d3..7fb854ee083 100644 --- a/Source/Particles/ParticleCreation/SmartCreate.H +++ b/Source/Particles/ParticleCreation/SmartCreate.H @@ -14,6 +14,8 @@ #include #include #include +#include +#include /** * \brief This is a functor for performing a "smart create" that works @@ -47,23 +49,23 @@ struct SmartCreate const int id = 0) const noexcept { #if defined(WARPX_DIM_3D) - prt.m_aos[i_prt].pos(0) = x; - prt.m_aos[i_prt].pos(1) = y; - prt.m_aos[i_prt].pos(2) = z; + prt.m_rdata[PIdx::x][i_prt] = x; + prt.m_rdata[PIdx::y][i_prt] = y; + prt.m_rdata[PIdx::z][i_prt] = z; #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - prt.m_aos[i_prt].pos(0) = x; - prt.m_aos[i_prt].pos(1) = z; + prt.m_rdata[PIdx::x][i_prt] = x; + prt.m_rdata[PIdx::z][i_prt] = z; amrex::ignore_unused(y); #else - prt.m_aos[i_prt].pos(0) = z; + prt.m_rdata[PIdx::z][i_prt] = z; amrex::ignore_unused(x,y); #endif - prt.m_aos[i_prt].cpu() = cpu; - prt.m_aos[i_prt].id() = id; + amrex::ParticleIDWrapper{prt.m_idcpu[i_prt]} = id; + amrex::ParticleCPUWrapper{prt.m_idcpu[i_prt]} = cpu; - // initialize the real components - for (int j = 0; j < PartData::NAR; ++j) { + // initialize the real components after position + for (int j = AMREX_SPACEDIM; j < PartData::NAR; ++j) { prt.m_rdata[j][i_prt] = initializeRealValue(m_policy_real[j], engine); } for (int j = 0; j < prt.m_num_runtime_real; ++j) { diff --git a/Source/Particles/ParticleCreation/SmartUtils.H b/Source/Particles/ParticleCreation/SmartUtils.H index 732a12bb729..6b3d396900d 100644 --- a/Source/Particles/ParticleCreation/SmartUtils.H +++ b/Source/Particles/ParticleCreation/SmartUtils.H @@ -17,6 +17,7 @@ #include #include #include +#include #include #include @@ -60,12 +61,12 @@ void setNewParticleIDs (PTile& ptile, int old_size, int num_added) } const int cpuid = amrex::ParallelDescriptor::MyProc(); - auto pp = ptile.GetArrayOfStructs()().data() + old_size; + auto ptd = ptile.getParticleTileData(); amrex::ParallelFor(num_added, [=] AMREX_GPU_DEVICE (int ip) noexcept { - auto& p = pp[ip]; - p.id() = pid+ip; - p.cpu() = cpuid; + auto const new_id = ip + old_size; + amrex::ParticleIDWrapper{ptd.m_idcpu[new_id]} = pid+ip; + amrex::ParticleCPUWrapper{ptd.m_idcpu[new_id]} = cpuid; }); } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index a12ae75f629..edf91a84526 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -268,11 +268,9 @@ public: * @param[in] engine the random engine, used in initialization of QED optical depths */ void DefaultInitializeRuntimeAttributes ( - amrex::ParticleTile, - NArrayReal, NArrayInt, - amrex::PinnedArenaAllocator>& pinned_tile, - int n_external_attr_real, - int n_external_attr_int) final; + typename ContainerLike::ParticleTileType& pinned_tile, + int n_external_attr_real, + int n_external_attr_int) final; /** * \brief Apply NCI Godfrey filter to all components of E and B before gather diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 0d9180c090e..fd55e369040 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -194,8 +194,8 @@ namespace * and avoid any possible undefined behavior before the next call to redistribute) and sets * the particle id to -1 so that it can be effectively deleted. * - * \param p particle aos data - * \param pa particle soa data + * \param idcpu particle id soa data + * \param pa particle real soa data * \param ip index for soa data * \param do_field_ionization whether species has ionization * \param pi ionization level data @@ -206,20 +206,21 @@ namespace */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ZeroInitializeAndSetNegativeID ( - ParticleType& p, const GpuArray& pa, long& ip, + uint64_t * AMREX_RESTRICT idcpu, + const GpuArray& pa, long& ip, const bool& do_field_ionization, int* pi #ifdef WARPX_QED - ,const bool& has_quantum_sync, amrex::ParticleReal* p_optical_depth_QSR - ,const bool& has_breit_wheeler, amrex::ParticleReal* p_optical_depth_BW + ,const bool& has_quantum_sync, amrex::ParticleReal* AMREX_RESTRICT p_optical_depth_QSR + ,const bool& has_breit_wheeler, amrex::ParticleReal* AMREX_RESTRICT p_optical_depth_BW #endif ) noexcept { - p.pos(0) = 0._rt; + pa[PIdx::z][ip] = 0._rt; #if (AMREX_SPACEDIM >= 2) - p.pos(1) = 0._rt; + pa[PIdx::x][ip] = 0._rt; #endif #if defined(WARPX_DIM_3D) - p.pos(2) = 0._rt; + pa[PIdx::y][ip] = 0._rt; #endif pa[PIdx::w ][ip] = 0._rt; pa[PIdx::ux][ip] = 0._rt; @@ -234,7 +235,7 @@ namespace if (has_breit_wheeler) {p_optical_depth_BW[ip] = 0._rt;} #endif - p.id() = -1; + amrex::ParticleIDWrapper{idcpu[ip]} = -1; } } @@ -776,11 +777,9 @@ PhysicalParticleContainer::AddPlasmaFromFile(PlasmaInjector & plasma_injector, void PhysicalParticleContainer::DefaultInitializeRuntimeAttributes ( - amrex::ParticleTile, - NArrayReal, NArrayInt, - amrex::PinnedArenaAllocator>& pinned_tile, - const int n_external_attr_real, - const int n_external_attr_int) + typename ContainerLike::ParticleTileType& pinned_tile, + int n_external_attr_real, + int n_external_attr_int) { ParticleCreation::DefaultInitializeRuntimeAttributes(pinned_tile, n_external_attr_real, n_external_attr_int, @@ -1080,7 +1079,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int const int max_new_particles = Scan::ExclusiveSum(counts.size(), counts.data(), offset.data()); // Update NextID to include particles created in this function - Long pid; + int pid; #ifdef AMREX_USE_OMP #pragma omp critical (add_plasma_nextid) #endif @@ -1089,7 +1088,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int ParticleType::NextID(pid+max_new_particles); } WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - static_cast(pid + max_new_particles) < LastParticleID, + static_cast(pid) + static_cast(max_new_particles) < LongParticleIds::LastParticleID, "ERROR: overflow on particle id numbers"); const int cpuid = ParallelDescriptor::MyProc(); @@ -1100,16 +1099,16 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int DefineAndReturnParticleTile(lev, grid_id, tile_id); } - auto old_size = particle_tile.GetArrayOfStructs().size(); + auto old_size = particle_tile.size(); auto new_size = old_size + max_new_particles; particle_tile.resize(new_size); - ParticleType* pp = particle_tile.GetArrayOfStructs()().data() + old_size; auto& soa = particle_tile.GetStructOfArrays(); GpuArray pa; for (int ia = 0; ia < PIdx::nattribs; ++ia) { pa[ia] = soa.GetRealData(ia).data() + old_size; } + uint64_t * AMREX_RESTRICT pa_idcpu = soa.GetIdCPUData().data() + old_size; // user-defined integer and real attributes const auto n_user_int_attribs = static_cast(m_user_int_attribs.size()); const auto n_user_real_attribs = static_cast(m_user_real_attribs.size()); @@ -1222,9 +1221,8 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int for (int i_part = 0; i_part < pcounts[index]; ++i_part) { long ip = poffset[index] + i_part; - ParticleType& p = pp[ip]; - p.id() = pid+ip; - p.cpu() = cpuid; + amrex::ParticleIDWrapper{pa_idcpu[ip]} = pid+ip; + amrex::ParticleCPUWrapper{pa_idcpu[ip]} = cpuid; const XDim3 r = (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) ? // In the refined injection region: use refinement ratio `lrrfac` inj_pos->getPositionUnitBox(i_part, lrrfac, engine) : @@ -1234,7 +1232,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int #if defined(WARPX_DIM_3D) if (!tile_realbox.contains(XDim3{pos.x,pos.y,pos.z})) { - ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi + ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED ,loc_has_quantum_sync, p_optical_depth_QSR ,loc_has_breit_wheeler, p_optical_depth_BW @@ -1245,7 +1243,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) amrex::ignore_unused(k); if (!tile_realbox.contains(XDim3{pos.x,pos.z,0.0_rt})) { - ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi + ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED ,loc_has_quantum_sync, p_optical_depth_QSR ,loc_has_breit_wheeler, p_optical_depth_BW @@ -1256,7 +1254,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int #else amrex::ignore_unused(j,k); if (!tile_realbox.contains(XDim3{pos.z,0.0_rt,0.0_rt})) { - ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi + ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED ,loc_has_quantum_sync, p_optical_depth_QSR ,loc_has_breit_wheeler, p_optical_depth_BW @@ -1295,7 +1293,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int const Real z0 = applyBallisticCorrection(pos, inj_mom, gamma_boost, beta_boost, t); if (!inj_pos->insideBounds(xb, yb, z0)) { - ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi + ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED ,loc_has_quantum_sync, p_optical_depth_QSR ,loc_has_breit_wheeler, p_optical_depth_BW @@ -1309,7 +1307,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int // Remove particle if density below threshold if ( dens < density_min ){ - ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi + ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED ,loc_has_quantum_sync, p_optical_depth_QSR ,loc_has_breit_wheeler, p_optical_depth_BW @@ -1327,7 +1325,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int // If the particle is not within the lab-frame zmin, zmax, etc. // go to the next generated particle. if (!inj_pos->insideBounds(xb, yb, z0_lab)) { - ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi + ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED ,loc_has_quantum_sync, p_optical_depth_QSR ,loc_has_breit_wheeler, p_optical_depth_BW @@ -1339,7 +1337,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int dens = inj_rho->getDensity(pos.x, pos.y, z0_lab); // Remove particle if density below threshold if ( dens < density_min ){ - ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi + ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED ,loc_has_quantum_sync, p_optical_depth_QSR ,loc_has_breit_wheeler, p_optical_depth_BW @@ -1406,17 +1404,17 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int pa[PIdx::uz][ip] = u.z; #if defined(WARPX_DIM_3D) - p.pos(0) = pos.x; - p.pos(1) = pos.y; - p.pos(2) = pos.z; + pa[PIdx::x][ip] = pos.x; + pa[PIdx::y][ip] = pos.y; + pa[PIdx::z][ip] = pos.z; #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) #ifdef WARPX_DIM_RZ pa[PIdx::theta][ip] = theta; #endif - p.pos(0) = xb; - p.pos(1) = pos.z; + pa[PIdx::x][ip] = xb; + pa[PIdx::z][ip] = pos.z; #else - p.pos(0) = pos.z; + pa[PIdx::z][ip] = pos.z; #endif } }); @@ -1635,7 +1633,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, const int max_new_particles = Scan::ExclusiveSum(counts.size(), counts.data(), offset.data()); // Update NextID to include particles created in this function - Long pid; + int pid; #ifdef AMREX_USE_OMP #pragma omp critical (add_plasma_nextid) #endif @@ -1644,23 +1642,23 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, ParticleType::NextID(pid+max_new_particles); } WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - static_cast(pid + max_new_particles) < LastParticleID, + static_cast(pid) + static_cast(max_new_particles) < LongParticleIds::LastParticleID, "overflow on particle id numbers"); const int cpuid = ParallelDescriptor::MyProc(); auto& particle_tile = tmp_pc.DefineAndReturnParticleTile(0, grid_id, tile_id); - auto old_size = particle_tile.GetArrayOfStructs().size(); + auto old_size = particle_tile.size(); auto new_size = old_size + max_new_particles; particle_tile.resize(new_size); - ParticleType* pp = particle_tile.GetArrayOfStructs()().data() + old_size; auto& soa = particle_tile.GetStructOfArrays(); GpuArray pa; for (int ia = 0; ia < PIdx::nattribs; ++ia) { pa[ia] = soa.GetRealData(ia).data() + old_size; } + uint64_t * AMREX_RESTRICT pa_idcpu = soa.GetIdCPUData().data() + old_size; // user-defined integer and real attributes const auto n_user_int_attribs = static_cast(m_user_int_attribs.size()); @@ -1758,9 +1756,8 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, for (int i_part = 0; i_part < pcounts[index]; ++i_part) { const long ip = poffset[index] + i_part; - ParticleType& p = pp[ip]; - p.id() = pid+ip; - p.cpu() = cpuid; + amrex::ParticleIDWrapper{pa_idcpu[ip]} = pid+ip; + amrex::ParticleCPUWrapper{pa_idcpu[ip]} = cpuid; // This assumes the flux_pos is of type InjectorPositionRandomPlane const XDim3 r = (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) ? @@ -1785,19 +1782,19 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, // the particles will be within the domain. #if defined(WARPX_DIM_3D) if (!ParticleUtils::containsInclusive(tile_realbox, XDim3{ppos.x,ppos.y,ppos.z})) { - p.id() = -1; + amrex::ParticleIDWrapper{pa_idcpu[ip]} = -1; continue; } #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) amrex::ignore_unused(k); if (!ParticleUtils::containsInclusive(tile_realbox, XDim3{ppos.x,ppos.z,0.0_prt})) { - p.id() = -1; + amrex::ParticleIDWrapper{pa_idcpu[ip]} = -1; continue; } #else amrex::ignore_unused(j,k); if (!ParticleUtils::containsInclusive(tile_realbox, XDim3{ppos.z,0.0_prt,0.0_prt})) { - p.id() = -1; + amrex::ParticleIDWrapper{pa_idcpu[ip]} = -1; continue; } #endif @@ -1805,7 +1802,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, // If the particle's initial position is not within or on the species's // xmin, xmax, ymin, ymax, zmin, zmax, go to the next generated particle. if (!flux_pos->insideBoundsInclusive(ppos.x, ppos.y, ppos.z)) { - p.id() = -1; + amrex::ParticleIDWrapper{pa_idcpu[ip]} = -1; continue; } @@ -1839,7 +1836,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, Real flux = inj_flux->getFlux(ppos.x, ppos.y, ppos.z, t); // Remove particle if flux is negative or 0 if ( flux <=0 ){ - p.id() = -1; + amrex::ParticleIDWrapper{pa_idcpu[ip]} = -1; continue; } @@ -1898,18 +1895,18 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, UpdatePosition(ppos.x, ppos.y, ppos.z, pu.x, pu.y, pu.z, t_fract); #if defined(WARPX_DIM_3D) - p.pos(0) = ppos.x; - p.pos(1) = ppos.y; - p.pos(2) = ppos.z; + pa[PIdx::x][ip] = ppos.x; + pa[PIdx::y][ip] = ppos.y; + pa[PIdx::z][ip] = ppos.z; #elif defined(WARPX_DIM_RZ) pa[PIdx::theta][ip] = std::atan2(ppos.y, ppos.x); - p.pos(0) = std::sqrt(ppos.x*ppos.x + ppos.y*ppos.y); - p.pos(1) = ppos.z; + pa[PIdx::x][ip] = std::sqrt(ppos.x*ppos.x + ppos.y*ppos.y); + pa[PIdx::z][ip] = ppos.z; #elif defined(WARPX_DIM_XZ) - p.pos(0) = ppos.x; - p.pos(1) = ppos.z; + pa[PIdx::x][ip] = ppos.x; + pa[PIdx::z][ip] = ppos.z; #else - p.pos(0) = ppos.z; + pa[PIdx::z][ip] = ppos.z; #endif } }); @@ -2326,20 +2323,22 @@ PhysicalParticleContainer::SplitParticles (int lev) split_offset[1] /= ppc_nd[1]; split_offset[2] /= ppc_nd[2]; } - // particle Array Of Structs data - auto& particles = pti.GetArrayOfStructs(); // particle Struct Of Arrays data auto& attribs = pti.GetAttribs(); auto& wp = attribs[PIdx::w ]; auto& uxp = attribs[PIdx::ux]; auto& uyp = attribs[PIdx::uy]; auto& uzp = attribs[PIdx::uz]; + + ParticleTileType& ptile = ParticlesAt(lev, pti); + auto& soa = ptile.GetStructOfArrays(); + uint64_t * const AMREX_RESTRICT idcpu = soa.GetIdCPUData().data(); + const long np = pti.numParticles(); for(int i=0; i> attr_int; pctmp_split.AddNParticles(lev, np_split_to_add, - xp, yp, zp, uxp, uyp, uzp, - 1, attr, + xp, + yp, + zp, + uxp, + uyp, + uzp, + 1, + attr, 0, attr_int, - 1, NoSplitParticleID); + 1, LongParticleIds::NoSplitParticleID); // Copy particles from tmp to current particle container constexpr bool local_flag = true; addParticles(pctmp_split,local_flag); diff --git a/Source/Particles/Pusher/GetAndSetPosition.H b/Source/Particles/Pusher/GetAndSetPosition.H index e4477a2a60d..44641557756 100644 --- a/Source/Particles/Pusher/GetAndSetPosition.H +++ b/Source/Particles/Pusher/GetAndSetPosition.H @@ -30,24 +30,26 @@ void get_particle_position (const WarpXParticleContainer::SuperParticleType& p, amrex::ParticleReal& y, amrex::ParticleReal& z) noexcept { -#ifdef WARPX_DIM_RZ - const amrex::ParticleReal theta = p.rdata(T_PIdx::theta); - const amrex::ParticleReal r = p.pos(0); + using namespace amrex::literals; + +#if defined(WARPX_DIM_RZ) + amrex::ParticleReal const theta = p.rdata(T_PIdx::theta); + amrex::ParticleReal const r = p.pos(T_PIdx::x); x = r*std::cos(theta); y = r*std::sin(theta); - z = p.pos(1); -#elif WARPX_DIM_3D - x = p.pos(0); - y = p.pos(1); - z = p.pos(2); -#elif WARPX_DIM_XZ - x = p.pos(0); - y = amrex::ParticleReal(0.0); - z = p.pos(1); + z = p.pos(PIdx::z); +#elif defined(WARPX_DIM_3D) + x = p.pos(PIdx::x); + y = p.pos(PIdx::y); + z = p.pos(PIdx::z); +#elif defined(WARPX_DIM_XZ) + x = p.pos(PIdx::x); + y = 0_prt; + z = p.pos(PIdx::z); #else - x = amrex::ParticleReal(0.0); - y = amrex::ParticleReal(0.0); - z = p.pos(0); + x = 0_prt; + y = 0_prt; + z = p.pos(PIdx::z); #endif } @@ -59,10 +61,19 @@ void get_particle_position (const WarpXParticleContainer::SuperParticleType& p, template struct GetParticlePosition { - using PType = WarpXParticleContainer::ParticleType; using RType = amrex::ParticleReal; - const PType* AMREX_RESTRICT m_structs = nullptr; +#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) + const RType* AMREX_RESTRICT m_x = nullptr; + const RType* AMREX_RESTRICT m_z = nullptr; +#elif defined(WARPX_DIM_3D) + const RType* AMREX_RESTRICT m_x = nullptr; + const RType* AMREX_RESTRICT m_y = nullptr; + const RType* AMREX_RESTRICT m_z = nullptr; +#elif defined(WARPX_DIM_1D_Z) + const RType* AMREX_RESTRICT m_z = nullptr; +#endif + #if defined(WARPX_DIM_RZ) const RType* m_theta = nullptr; #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) @@ -84,10 +95,19 @@ struct GetParticlePosition template GetParticlePosition (const ptiType& a_pti, long a_offset = 0) noexcept { - const auto& aos = a_pti.GetArrayOfStructs(); - m_structs = aos().dataPtr() + a_offset; -#if defined(WARPX_DIM_RZ) const auto& soa = a_pti.GetStructOfArrays(); + +#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) + m_x = soa.GetRealData(PIdx::x).dataPtr() + a_offset; + m_z = soa.GetRealData(PIdx::z).dataPtr() + a_offset; +#elif defined(WARPX_DIM_3D) + m_x = soa.GetRealData(PIdx::x).dataPtr() + a_offset; + m_y = soa.GetRealData(PIdx::y).dataPtr() + a_offset; + m_z = soa.GetRealData(PIdx::z).dataPtr() + a_offset; +#elif defined(WARPX_DIM_1D_Z) + m_z = soa.GetRealData(PIdx::z).dataPtr() + a_offset; +#endif +#if defined(WARPX_DIM_RZ) m_theta = soa.GetRealData(T_PIdx::theta).dataPtr() + a_offset; #endif } @@ -98,24 +118,23 @@ struct GetParticlePosition AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() (const long i, RType& x, RType& y, RType& z) const noexcept { - const PType& p = m_structs[i]; #ifdef WARPX_DIM_RZ - const RType r = p.pos(0); + RType const r = m_x[i]; x = r*std::cos(m_theta[i]); y = r*std::sin(m_theta[i]); - z = p.pos(1); + z = m_z[i]; #elif WARPX_DIM_3D - x = p.pos(0); - y = p.pos(1); - z = p.pos(2); + x = m_x[i]; + y = m_y[i]; + z = m_z[i]; #elif WARPX_DIM_XZ - x = p.pos(0); + x = m_x[i]; y = m_y_default; - z = p.pos(1); + z = m_z[i]; #else x = m_x_default; y = m_y_default; - z = p.pos(0); + z = m_z[i]; #endif } @@ -127,23 +146,22 @@ struct GetParticlePosition AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void AsStored (const long i, RType& x, RType& y, RType& z) const noexcept { - const PType& p = m_structs[i]; #ifdef WARPX_DIM_RZ - x = p.pos(0); + x = m_x[i]; y = m_theta[i]; - z = p.pos(1); + z = m_z[i]; #elif WARPX_DIM_3D - x = p.pos(0); - y = p.pos(1); - z = p.pos(2); + x = m_x[i]; + y = m_y[i]; + z = m_z[i]; #elif WARPX_DIM_XZ - x = p.pos(0); + x = m_x[i]; y = m_y_default; - z = p.pos(1); + z = m_z[i]; #else x = m_x_default; y = m_y_default; - z = p.pos(0); + z = m_z[i]; #endif } }; @@ -158,10 +176,18 @@ struct GetParticlePosition template struct SetParticlePosition { - using PType = WarpXParticleContainer::ParticleType; using RType = amrex::ParticleReal; - PType* AMREX_RESTRICT m_structs; +#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) + RType* AMREX_RESTRICT m_x; + RType* AMREX_RESTRICT m_z; +#elif defined(WARPX_DIM_3D) + RType* AMREX_RESTRICT m_x; + RType* AMREX_RESTRICT m_y; + RType* AMREX_RESTRICT m_z; +#elif defined(WARPX_DIM_1D_Z) + RType* AMREX_RESTRICT m_z; +#endif #if defined(WARPX_DIM_RZ) RType* AMREX_RESTRICT m_theta; #endif @@ -169,10 +195,18 @@ struct SetParticlePosition template SetParticlePosition (const ptiType& a_pti, long a_offset = 0) noexcept { - auto& aos = a_pti.GetArrayOfStructs(); - m_structs = aos().dataPtr() + a_offset; -#if defined(WARPX_DIM_RZ) auto& soa = a_pti.GetStructOfArrays(); +#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) + m_x = soa.GetRealData(PIdx::x).dataPtr() + a_offset; + m_z = soa.GetRealData(PIdx::z).dataPtr() + a_offset; +#elif defined(WARPX_DIM_3D) + m_x = soa.GetRealData(PIdx::x).dataPtr() + a_offset; + m_y = soa.GetRealData(PIdx::y).dataPtr() + a_offset; + m_z = soa.GetRealData(PIdx::z).dataPtr() + a_offset; +#elif defined(WARPX_DIM_1D_Z) + m_z = soa.GetRealData(PIdx::z).dataPtr() + a_offset; +#endif +#if defined(WARPX_DIM_RZ) m_theta = soa.GetRealData(T_PIdx::theta).dataPtr() + a_offset; #endif } @@ -190,17 +224,17 @@ struct SetParticlePosition #endif #ifdef WARPX_DIM_RZ m_theta[i] = std::atan2(y, x); - m_structs[i].pos(0) = std::sqrt(x*x + y*y); - m_structs[i].pos(1) = z; + m_x[i] = std::sqrt(x*x + y*y); + m_z[i] = z; #elif WARPX_DIM_3D - m_structs[i].pos(0) = x; - m_structs[i].pos(1) = y; - m_structs[i].pos(2) = z; + m_x[i] = x; + m_y[i] = y; + m_z[i] = z; #elif WARPX_DIM_XZ - m_structs[i].pos(0) = x; - m_structs[i].pos(1) = z; + m_x[i] = x; + m_z[i] = z; #else - m_structs[i].pos(0) = z; + m_z[i] = z; #endif } @@ -218,18 +252,18 @@ struct SetParticlePosition amrex::ignore_unused(x,y); #endif #ifdef WARPX_DIM_RZ - m_structs[i].pos(0) = x; + m_x[i] = x; m_theta[i] = y; - m_structs[i].pos(1) = z; + m_z[i] = z; #elif WARPX_DIM_3D - m_structs[i].pos(0) = x; - m_structs[i].pos(1) = y; - m_structs[i].pos(2) = z; + m_x[i] = x; + m_y[i] = y; + m_z[i] = z; #elif WARPX_DIM_XZ - m_structs[i].pos(0) = x; - m_structs[i].pos(1) = z; + m_x[i] = x; + m_z[i] = z; #else - m_structs[i].pos(0) = z; + m_z[i] = z; #endif } }; diff --git a/Source/Particles/Resampling/LevelingThinning.cpp b/Source/Particles/Resampling/LevelingThinning.cpp index 680e33ebe6a..40f75c76eab 100644 --- a/Source/Particles/Resampling/LevelingThinning.cpp +++ b/Source/Particles/Resampling/LevelingThinning.cpp @@ -60,8 +60,7 @@ void LevelingThinning::operator() (WarpXParIter& pti, const int lev, auto& ptile = pc->ParticlesAt(lev, pti); auto& soa = ptile.GetStructOfArrays(); amrex::ParticleReal * const AMREX_RESTRICT w = soa.GetRealData(PIdx::w).data(); - WarpXParticleContainer::ParticleType * const AMREX_RESTRICT - particle_ptr = ptile.GetArrayOfStructs()().data(); + auto * const AMREX_RESTRICT idcpu = soa.GetIdCPUData().data(); // Using this function means that we must loop over the cells in the ParallelFor. In the case // of the leveling thinning algorithm, it would have possibly been more natural and more @@ -114,7 +113,7 @@ void LevelingThinning::operator() (WarpXParIter& pti, const int lev, // Remove particle with probability 1 - particle_weight/level_weight if (random_number > w[indices[i]]/level_weight) { - particle_ptr[indices[i]].id() = -1; + amrex::ParticleIDWrapper{idcpu[indices[i]]} = -1; } // Set particle weight to level weight otherwise else diff --git a/Source/Particles/Sorting/Partition.cpp b/Source/Particles/Sorting/Partition.cpp index 58511cfd5e7..58e3450f47d 100644 --- a/Source/Particles/Sorting/Partition.cpp +++ b/Source/Particles/Sorting/Partition.cpp @@ -61,7 +61,7 @@ PhysicalParticleContainer::PartitionParticlesInBuffers( // Initialize temporary arrays Gpu::DeviceVector inexflag; inexflag.resize(np); - Gpu::DeviceVector pid; + Gpu::DeviceVector pid; pid.resize(np); // First, partition particles into the larger buffer @@ -109,7 +109,7 @@ PhysicalParticleContainer::PartitionParticlesInBuffers( // - For each particle in the large buffer, find whether it is in // the smaller buffer, by looking up the mask. Store the answer in `inexflag`. amrex::ParallelFor( np - n_fine, - fillBufferFlagRemainingParticles(pti, bmasks, inexflag, Geom(lev), pid, n_fine) ); + fillBufferFlagRemainingParticles(pti, bmasks, inexflag, Geom(lev), pid, int(n_fine)) ); auto *const sep2 = stablePartition( sep, pid.end(), inexflag ); if (bmasks == gather_masks) { diff --git a/Source/Particles/Sorting/SortingUtils.H b/Source/Particles/Sorting/SortingUtils.H index ac2c63e88f8..ba7761bf48a 100644 --- a/Source/Particles/Sorting/SortingUtils.H +++ b/Source/Particles/Sorting/SortingUtils.H @@ -12,6 +12,7 @@ #include #include +#include /** \brief Fill the elements of the input vector with consecutive integer, @@ -19,7 +20,7 @@ * * \param[inout] v Vector of integers, to be filled by this routine */ -void fillWithConsecutiveIntegers( amrex::Gpu::DeviceVector& v ); +void fillWithConsecutiveIntegers( amrex::Gpu::DeviceVector& v ); /** \brief Find the indices that would reorder the elements of `predicate` * so that the elements with non-zero value precede the other elements @@ -41,7 +42,7 @@ ForwardIterator stablePartition(ForwardIterator const index_begin, int const* AMREX_RESTRICT predicate_ptr = predicate.dataPtr(); int N = static_cast(std::distance(index_begin, index_end)); auto num_true = amrex::StablePartition(&(*index_begin), N, - [predicate_ptr] AMREX_GPU_DEVICE (long i) { return predicate_ptr[i]; }); + [predicate_ptr] AMREX_GPU_DEVICE (int i) { return predicate_ptr[i]; }); ForwardIterator sep = index_begin; std::advance(sep, num_true); @@ -49,7 +50,7 @@ ForwardIterator stablePartition(ForwardIterator const index_begin, // On CPU: Use std library ForwardIterator const sep = std::stable_partition( index_begin, index_end, - [&predicate](long i) { return predicate[i]; } + [&predicate](int i) { return predicate[i]; } ); #endif return sep; @@ -88,7 +89,7 @@ class fillBufferFlag // Extract simple structure that can be used directly on the GPU m_domain{geom.Domain()}, m_inexflag_ptr{inexflag.dataPtr()}, - m_particles{pti.GetArrayOfStructs().data()}, + m_ptd{pti.GetParticleTile().getConstParticleTileData()}, m_buffer_mask{(*bmasks)[pti].array()} { for (int idim=0; idim m_buffer_mask; amrex::GpuArray m_prob_lo; amrex::GpuArray m_inv_cell_size; @@ -141,12 +140,12 @@ class fillBufferFlagRemainingParticles amrex::iMultiFab const* bmasks, amrex::Gpu::DeviceVector& inexflag, amrex::Geometry const& geom, - amrex::Gpu::DeviceVector const& particle_indices, - long const start_index ) : + amrex::Gpu::DeviceVector const& particle_indices, + int start_index ) : m_domain{geom.Domain()}, // Extract simple structure that can be used directly on the GPU m_inexflag_ptr{inexflag.dataPtr()}, - m_particles{pti.GetArrayOfStructs().data()}, + m_ptd{pti.GetParticleTile().getConstParticleTileData()}, m_buffer_mask{(*bmasks)[pti].array()}, m_start_index{start_index}, m_indices_ptr{particle_indices.dataPtr()} @@ -159,11 +158,11 @@ class fillBufferFlagRemainingParticles AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - void operator()( const long i ) const { + void operator()( const int i ) const { // Select a particle - auto const& p = m_particles[m_indices_ptr[i+m_start_index]]; + auto const j = m_indices_ptr[i+m_start_index]; // Find the index of the cell where this particle is located - amrex::IntVect const iv = amrex::getParticleCell( p, + amrex::IntVect const iv = amrex::getParticleCell( m_ptd, j, m_prob_lo, m_inv_cell_size, m_domain ); // Find the value of the buffer flag in this cell and // store it at the corresponding particle position in the array `inexflag` @@ -175,10 +174,10 @@ class fillBufferFlagRemainingParticles amrex::GpuArray m_inv_cell_size; amrex::Box m_domain; int* m_inexflag_ptr; - WarpXParticleContainer::ParticleType const* m_particles; + WarpXParticleContainer::ParticleTileType::ConstParticleTileDataType const m_ptd; amrex::Array4 m_buffer_mask; - long const m_start_index; - long const* m_indices_ptr; + int const m_start_index; + int const* m_indices_ptr; }; /** \brief Functor that copies the elements of `src` into `dst`, @@ -195,7 +194,7 @@ class copyAndReorder copyAndReorder( amrex::Gpu::DeviceVector const& src, amrex::Gpu::DeviceVector& dst, - amrex::Gpu::DeviceVector const& indices ): + amrex::Gpu::DeviceVector const& indices ): // Extract simple structure that can be used directly on the GPU m_src_ptr{src.dataPtr()}, m_dst_ptr{dst.dataPtr()}, @@ -203,14 +202,14 @@ class copyAndReorder {} AMREX_GPU_DEVICE AMREX_FORCE_INLINE - void operator()( const long ip ) const { + void operator()( const int ip ) const { m_dst_ptr[ip] = m_src_ptr[ m_indices_ptr[ip] ]; } private: T const* m_src_ptr; T* m_dst_ptr; - long const* m_indices_ptr; + int const* m_indices_ptr; }; #endif // WARPX_PARTICLES_SORTING_SORTINGUTILS_H_ diff --git a/Source/Particles/Sorting/SortingUtils.cpp b/Source/Particles/Sorting/SortingUtils.cpp index 699119e8e18..cd4b6a13c76 100644 --- a/Source/Particles/Sorting/SortingUtils.cpp +++ b/Source/Particles/Sorting/SortingUtils.cpp @@ -8,7 +8,7 @@ #include "SortingUtils.H" -void fillWithConsecutiveIntegers( amrex::Gpu::DeviceVector& v ) +void fillWithConsecutiveIntegers( amrex::Gpu::DeviceVector& v ) { #ifdef AMREX_USE_GPU // On GPU: Use amrex diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 33aa71d1c7d..a244afc2389 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -49,10 +49,10 @@ class WarpXParIter - : public amrex::ParIter<0,0,PIdx::nattribs> + : public amrex::ParIterSoA { public: - using amrex::ParIter<0,0,PIdx::nattribs>::ParIter; + using amrex::ParIterSoA::ParIterSoA; WarpXParIter (ContainerType& pc, int level); @@ -89,7 +89,7 @@ public: * particle container classes (that store a collection of particles) derive. Derived * classes can be used for plasma particles, photon particles, or non-physical * particles (e.g., for the laser antenna). - * It derives from amrex::ParticleContainer<0,0,PIdx::nattribs>, where the + * It derives from amrex::ParticleContainerPureSoA, where the * template arguments stand for the number of int and amrex::Real SoA and AoS * data in amrex::Particle. * - AoS amrex::Real: x, y, z (default), 0 additional (first template @@ -164,11 +164,9 @@ public: * class. */ virtual void DefaultInitializeRuntimeAttributes ( - amrex::ParticleTile, - NArrayReal, NArrayInt, - amrex::PinnedArenaAllocator>& pinned_tile, - int n_external_attr_real, - int n_external_attr_int) = 0; + typename ContainerLike::ParticleTileType& pinned_tile, + int n_external_attr_real, + int n_external_attr_int) = 0; /// /// This pushes the particle positions by one half time step. diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index e94aae5d573..48344b20fef 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -72,13 +72,13 @@ using namespace amrex; WarpXParIter::WarpXParIter (ContainerType& pc, int level) - : amrex::ParIter<0,0,PIdx::nattribs>(pc, level, + : amrex::ParIterSoA(pc, level, MFItInfo().SetDynamic(WarpX::do_dynamic_scheduling)) { } WarpXParIter::WarpXParIter (ContainerType& pc, int level, MFItInfo& info) - : amrex::ParIter<0,0,PIdx::nattribs>(pc, level, + : amrex::ParIterSoA(pc, level, info.SetDynamic(WarpX::do_dynamic_scheduling)) { } @@ -195,52 +195,54 @@ WarpXParticleContainer::AddNParticles (int /*lev*/, long n, // Redistribute() will move them to proper places. auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - using PinnedTile = amrex::ParticleTile, - NArrayReal, NArrayInt, - amrex::PinnedArenaAllocator>; + using PinnedTile = typename ContainerLike::ParticleTileType; PinnedTile pinned_tile; pinned_tile.define(NumRuntimeRealComps(), NumRuntimeIntComps()); const std::size_t np = iend-ibegin; #ifdef WARPX_DIM_RZ + amrex::Vector r(np); amrex::Vector theta(np); #endif for (auto i = ibegin; i < iend; ++i) { - ParticleType p; - if (id==-1) - { - p.id() = ParticleType::NextID(); + auto & idcpu_data = pinned_tile.GetStructOfArrays().GetIdCPUData(); + idcpu_data.push_back(0); + if (id==-1) { + amrex::ParticleIDWrapper{idcpu_data.back()} = ParticleType::NextID(); } else { - p.id() = id; + amrex::ParticleIDWrapper{idcpu_data.back()} = id; } - p.cpu() = amrex::ParallelDescriptor::MyProc(); + amrex::ParticleCPUWrapper(idcpu_data.back()) = ParallelDescriptor::MyProc(); + +#ifdef WARPX_DIM_RZ + r[i-ibegin] = std::sqrt(x[i]*x[i] + y[i]*y[i]); + theta[i-ibegin] = std::atan2(y[i], x[i]); +#endif + } + + if (np > 0) + { #if defined(WARPX_DIM_3D) - p.pos(0) = x[i]; - p.pos(1) = y[i]; - p.pos(2) = z[i]; + pinned_tile.push_back_real(PIdx::x, x.data() + ibegin, x.data() + iend); + pinned_tile.push_back_real(PIdx::y, y.data() + ibegin, y.data() + iend); + pinned_tile.push_back_real(PIdx::z, z.data() + ibegin, z.data() + iend); #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) amrex::ignore_unused(y); #ifdef WARPX_DIM_RZ - theta[i-ibegin] = std::atan2(y[i], x[i]); - p.pos(0) = std::sqrt(x[i]*x[i] + y[i]*y[i]); + pinned_tile.push_back_real(PIdx::x, r.data(), r.data() + np); #else - p.pos(0) = x[i]; + pinned_tile.push_back_real(PIdx::x, x.data() + ibegin, x.data() + iend); #endif - p.pos(1) = z[i]; + pinned_tile.push_back_real(PIdx::z, z.data() + ibegin, z.data() + iend); #else //AMREX_SPACEDIM == 1 amrex::ignore_unused(x,y); - p.pos(0) = z[i]; + pinned_tile.push_back_real(PIdx::z, z.data() + ibegin, z.data() + iend); #endif - pinned_tile.push_back(p); - } - - if (np > 0) - { - pinned_tile.push_back_real(PIdx::w , attr_real[0].data() + ibegin, attr_real[0].data() + iend); + pinned_tile.push_back_real(PIdx::w, attr_real[0].data() + ibegin, attr_real[0].data() + iend); pinned_tile.push_back_real(PIdx::ux, ux.data() + ibegin, ux.data() + iend); pinned_tile.push_back_real(PIdx::uy, uy.data() + ibegin, uy.data() + iend); pinned_tile.push_back_real(PIdx::uz, uz.data() + ibegin, uz.data() + iend); @@ -466,15 +468,14 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, //sort particles by bin WARPX_PROFILE_VAR_START(blp_sort); - amrex::DenseBins bins; + amrex::DenseBins bins; { auto& ptile = ParticlesAt(lev, pti); - auto& aos = ptile.GetArrayOfStructs(); - auto *pstruct_ptr = aos().dataPtr(); + auto ptd = ptile.getParticleTileData(); const int ntiles = numTilesInBox(box, true, bin_size); - bins.build(ptile.numParticles(), pstruct_ptr, ntiles, + bins.build(ptile.numParticles(), ptd, ntiles, [=] AMREX_GPU_HOST_DEVICE (const ParticleType& p) -> unsigned int { Box tbox; @@ -874,7 +875,7 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector const& wp, // HACK - sort particles by bin here. WARPX_PROFILE_VAR_START(blp_sort); - amrex::DenseBins bins; + amrex::DenseBins bins; { const Geometry& geom = Geom(lev); const auto dxi = geom.InvCellSizeArray(); @@ -882,16 +883,15 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector const& wp, const auto domain = geom.Domain(); auto& ptile = ParticlesAt(lev, pti); - auto& aos = ptile.GetArrayOfStructs(); - auto *pstruct_ptr = aos().dataPtr(); + auto ptd = ptile.getParticleTileData(); Box box = pti.validbox(); box.grow(ng_rho); const amrex::IntVect bin_size = WarpX::shared_tilesize; const int ntiles = numTilesInBox(box, true, bin_size); - bins.build(ptile.numParticles(), pstruct_ptr, ntiles, - [=] AMREX_GPU_HOST_DEVICE (const ParticleType& p) -> unsigned int + bins.build(ptile.numParticles(), ptd, ntiles, + [=] AMREX_GPU_HOST_DEVICE (ParticleType const & p) -> unsigned int { Box tbx; auto iv = getParticleCell(p, plo, dxi, domain); @@ -911,8 +911,7 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector const& wp, const auto domain = geom.Domain(); auto& ptile = ParticlesAt(lev, pti); - auto& aos = ptile.GetArrayOfStructs(); - auto *pstruct_ptr = aos().dataPtr(); + auto ptd = ptile.getParticleTileData(); Box box = pti.validbox(); box.grow(ng_rho); @@ -926,9 +925,10 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector const& wp, const auto bin_start = offsets_ptr[ibin]; const auto bin_stop = offsets_ptr[ibin+1]; if (bin_start < bin_stop) { - auto p = pstruct_ptr[permutation[bin_start]]; + // static_cast until https://github.com/AMReX-Codes/amrex/pull/3684 + auto const i = static_cast(permutation[bin_start]); Box tbx; - auto iv = getParticleCell(p, plo, dxi, domain); + auto iv = getParticleCell(ptd, i, plo, dxi, domain); AMREX_ASSERT(box.contains(iv)); [[maybe_unused]] auto tid = getTileIndex(iv, box, true, bin_size, tbx); AMREX_ASSERT(tid == ibin); @@ -1417,10 +1417,10 @@ WarpXParticleContainer::particlePostLocate(ParticleType& p, // Tag particle if goes to higher level. // It will be split later in the loop if (pld.m_lev == lev+1 - and p.id() != NoSplitParticleID + and p.id() != amrex::LongParticleIds::NoSplitParticleID and p.id() >= 0) { - p.id() = DoSplitParticleID; + p.id() = amrex::LongParticleIds::DoSplitParticleID; } if (pld.m_lev == lev-1){ @@ -1459,9 +1459,9 @@ WarpXParticleContainer::ApplyBoundaryConditions (){ const Real zmax = Geom(lev).ProbHi(WARPX_ZINDEX); ParticleTileType& ptile = ParticlesAt(lev, pti); - ParticleType * const pp = ptile.GetArrayOfStructs()().data(); auto& soa = ptile.GetStructOfArrays(); + uint64_t * const AMREX_RESTRICT idcpu = soa.GetIdCPUData().data(); amrex::ParticleReal * const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data(); amrex::ParticleReal * const AMREX_RESTRICT uy = soa.GetRealData(PIdx::uy).data(); amrex::ParticleReal * const AMREX_RESTRICT uz = soa.GetRealData(PIdx::uz).data(); @@ -1470,10 +1470,9 @@ WarpXParticleContainer::ApplyBoundaryConditions (){ amrex::ParallelForRNG( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i, amrex::RandomEngine const& engine) { - ParticleType& p = pp[i]; - // skip particles that are already flagged for removal - if (p.id() < 0) { return; } + auto const id = amrex::ParticleIDWrapper{idcpu[i]}; + if (id < 0) { return; } ParticleReal x, y, z; GetPosition.AsStored(i, x, y, z); @@ -1495,7 +1494,7 @@ WarpXParticleContainer::ApplyBoundaryConditions (){ boundary_conditions, engine); if (particle_lost) { - p.id() = -p.id(); + amrex::ParticleIDWrapper{idcpu[i]} = -id; } else { SetPosition.AsStored(i, x, y, z); } diff --git a/Source/Python/Particles/ParticleBoundaryBuffer.cpp b/Source/Python/Particles/ParticleBoundaryBuffer.cpp index 2a35faece9b..b04ac75e600 100644 --- a/Source/Python/Particles/ParticleBoundaryBuffer.cpp +++ b/Source/Python/Particles/ParticleBoundaryBuffer.cpp @@ -10,13 +10,13 @@ namespace warpx { class BoundaryBufferParIter - : public amrex::ParIter<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator> + : public amrex::ParIterSoA { public: - using amrex::ParIter<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator>::ParIter; + using amrex::ParIterSoA::ParIterSoA; BoundaryBufferParIter(ContainerType& pc, int level) : - amrex::ParIter<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator>(pc, level) {} + amrex::ParIterSoA(pc, level) {} }; } @@ -24,9 +24,9 @@ void init_BoundaryBufferParIter (py::module& m) { py::class_< warpx::BoundaryBufferParIter, - amrex::ParIter<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator> + amrex::ParIterSoA >(m, "BoundaryBufferParIter") - .def(py::init::ContainerType&, int>(), + .def(py::init::ContainerType&, int>(), py::arg("particle_container"), py::arg("level") ) ; diff --git a/Source/Python/Particles/PinnedMemoryParticleContainer.cpp b/Source/Python/Particles/PinnedMemoryParticleContainer.cpp index 600d56a62c9..d4f6a422dbe 100644 --- a/Source/Python/Particles/PinnedMemoryParticleContainer.cpp +++ b/Source/Python/Particles/PinnedMemoryParticleContainer.cpp @@ -13,6 +13,6 @@ void init_PinnedMemoryParticleContainer (py::module& m) { py::class_< PinnedMemoryParticleContainer, - amrex::ParticleContainer<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator> + amrex::ParticleContainerPureSoA > pmpc (m, "PinnedMemoryParticleContainer"); } diff --git a/Source/Python/Particles/WarpXParticleContainer.cpp b/Source/Python/Particles/WarpXParticleContainer.cpp index 1473a750941..07793a373f3 100644 --- a/Source/Python/Particles/WarpXParticleContainer.cpp +++ b/Source/Python/Particles/WarpXParticleContainer.cpp @@ -12,11 +12,11 @@ void init_WarpXParIter (py::module& m) { py::class_< - WarpXParIter, amrex::ParIter<0,0,PIdx::nattribs> + WarpXParIter, amrex::ParIterSoA >(m, "WarpXParIter") - .def(py::init::ContainerType&, int>(), + .def(py::init::ContainerType&, int>(), py::arg("particle_container"), py::arg("level")) - .def(py::init::ContainerType&, int, amrex::MFItInfo&>(), + .def(py::init::ContainerType&, int, amrex::MFItInfo&>(), py::arg("particle_container"), py::arg("level"), py::arg("info")) ; @@ -26,11 +26,11 @@ void init_WarpXParticleContainer (py::module& m) { py::class_< WarpXParticleContainer, - amrex::ParticleContainer<0, 0, PIdx::nattribs, 0> + amrex::ParticleContainerPureSoA > wpc (m, "WarpXParticleContainer"); wpc .def("add_real_comp", - [](WarpXParticleContainer& pc, const std::string& name, bool const comm) { pc.AddRealComp(name, comm); }, + [](WarpXParticleContainer& pc, const std::string& name, bool comm) { pc.AddRealComp(name, comm); }, py::arg("name"), py::arg("comm") ) .def("add_n_particles", @@ -93,6 +93,14 @@ void init_WarpXParticleContainer (py::module& m) }, py::arg("comp_name") ) + .def("get_icomp_index", + [](WarpXParticleContainer& pc, std::string comp_name) + { + auto particle_comps = pc.getParticleiComps(); + return particle_comps.at(comp_name); + }, + py::arg("comp_name") + ) .def("num_local_tiles_at_level", &WarpXParticleContainer::numLocalTilesAtLevel, py::arg("level") diff --git a/Source/Utils/ParticleUtils.H b/Source/Utils/ParticleUtils.H index b04176d4d83..7e3c89228ea 100644 --- a/Source/Utils/ParticleUtils.H +++ b/Source/Utils/ParticleUtils.H @@ -28,9 +28,10 @@ namespace ParticleUtils { * @param[in] mfi the MultiFAB iterator. * @param[in] ptile the particle tile. */ - amrex::DenseBins - findParticlesInEachCell(int lev, amrex::MFIter const& mfi, - WarpXParticleContainer::ParticleTileType const& ptile); + amrex::DenseBins + findParticlesInEachCell (int lev, + amrex::MFIter const & mfi, + WarpXParticleContainer::ParticleTileType & ptile); /** * \brief Return (relativistic) particle energy given velocity and mass. diff --git a/Source/Utils/ParticleUtils.cpp b/Source/Utils/ParticleUtils.cpp index 60e04f12b86..b8207b61fa0 100644 --- a/Source/Utils/ParticleUtils.cpp +++ b/Source/Utils/ParticleUtils.cpp @@ -22,24 +22,28 @@ #include #include -namespace ParticleUtils { +namespace ParticleUtils +{ using namespace amrex; + // Define shortcuts for frequently-used type names - using ParticleType = WarpXParticleContainer::ParticleType; - using ParticleTileType = WarpXParticleContainer::ParticleTileType; - using ParticleBins = DenseBins; - using index_type = ParticleBins::index_type; + using ParticleType = typename WarpXParticleContainer::ParticleType; + using ParticleTileType = typename WarpXParticleContainer::ParticleTileType; + using ParticleTileDataType = typename ParticleTileType::ParticleTileDataType; + using ParticleBins = DenseBins; + using index_type = typename ParticleBins::index_type; /* Find the particles and count the particles that are in each cell. Note that this does *not* rearrange particle arrays */ ParticleBins - findParticlesInEachCell( int const lev, MFIter const& mfi, - ParticleTileType const& ptile) { + findParticlesInEachCell (int lev, + MFIter const & mfi, + ParticleTileType & ptile) { // Extract particle structures for this tile int const np = ptile.numParticles(); - ParticleType const* particle_ptr = ptile.GetArrayOfStructs()().data(); + auto ptd = ptile.getParticleTileData(); // Extract box properties Geometry const& geom = WarpX::GetInstance().Geom(lev); @@ -51,9 +55,9 @@ namespace ParticleUtils { // Find particles that are in each cell; // results are stored in the object `bins`. ParticleBins bins; - bins.build(np, particle_ptr, cbx, + bins.build(np, ptd, cbx, // Pass lambda function that returns the cell index - [=] AMREX_GPU_DEVICE (const ParticleType& p) noexcept + [=] AMREX_GPU_DEVICE (ParticleType const & p) noexcept -> amrex::IntVect { return IntVect{AMREX_D_DECL( static_cast((p.pos(0)-plo[0])*dxi[0] - lo.x), @@ -64,4 +68,4 @@ namespace ParticleUtils { return bins; } -} +} // namespace ParticleUtils diff --git a/Source/ablastr/particles/IndexHandling.H b/Source/ablastr/particles/IndexHandling.H deleted file mode 100644 index 0ad5ca60446..00000000000 --- a/Source/ablastr/particles/IndexHandling.H +++ /dev/null @@ -1,41 +0,0 @@ -/* Copyright 2019-2022 Axel Huebl - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ -#ifndef ABLASTR_INDEX_HANDLING_H -#define ABLASTR_INDEX_HANDLING_H - -#include - - -namespace ablastr::particles { - - /** A helper function to derive a globally unique particle ID - * - * @param[in] id AMReX particle ID (on local cpu/rank), AoS .id - * @param[in] cpu AMReX particle CPU (rank) at creation of the particle, AoS .cpu - * @return global particle ID that is unique and permanent in the whole simulation - */ - constexpr uint64_t - localIDtoGlobal (int const id, int const cpu) - { - static_assert(sizeof(int) * 2u <= sizeof(uint64_t), - "int size might cause collisions in global IDs"); - // implementation: - // - we cast both 32-bit (or smaller) ints to a 64bit unsigned int - // - this will leave half of the "upper range" bits in the 64bit unsigned int zeroed out - // because the corresponding (extended) value range was not part of the value range in - // the int representation - // - we bit-shift the cpu into the upper half of zero bits in the 64 bit unsigned int - // (imagine this step as "adding a per-cpu/rank offset to the local integers") - // - then we add this offset - // note: the add is expressed as bitwise OR (|) since this saves us from writing - // brackets for operator precedence between + and << - return uint64_t(id) | uint64_t(cpu) << 32u; - } - -} // namespace ablastr::particles - -#endif // ABLASTR_INDEX_HANDLING_H diff --git a/Source/ablastr/particles/ParticleMoments.H b/Source/ablastr/particles/ParticleMoments.H index e45fb574cce..b648ccb28aa 100644 --- a/Source/ablastr/particles/ParticleMoments.H +++ b/Source/ablastr/particles/ParticleMoments.H @@ -35,7 +35,7 @@ namespace particles { amrex::ParticleReal, amrex::ParticleReal> MinAndMaxPositions (T_PC const & pc) { - using PType = typename T_PC::SuperParticleType; + using ConstParticleTileDataType = typename T_PC::ParticleTileType::ConstParticleTileDataType; // Get min and max for the local rank amrex::ReduceOps< @@ -46,11 +46,11 @@ namespace particles { amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal> >( pc, - [=] AMREX_GPU_DEVICE(PType const & p) noexcept + [=] AMREX_GPU_DEVICE(const ConstParticleTileDataType& ptd, const int i) noexcept { - amrex::ParticleReal const x = p.pos(0); - amrex::ParticleReal const y = p.pos(1); - amrex::ParticleReal const z = p.pos(2); + const amrex::ParticleReal x = ptd.rdata(0)[i]; + const amrex::ParticleReal y = ptd.rdata(1)[i]; + const amrex::ParticleReal z = ptd.rdata(2)[i]; return amrex::makeTuple(x, y, z, x, y, z); }, @@ -90,7 +90,8 @@ namespace particles { amrex::ParticleReal, amrex::ParticleReal> MeanAndStdPositions (T_PC const & pc) { - using PType = typename T_PC::SuperParticleType; + + using ConstParticleTileDataType = typename T_PC::ParticleTileType::ConstParticleTileDataType; amrex::ReduceOps< amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum, @@ -103,12 +104,14 @@ namespace particles { amrex::ParticleReal> >( pc, - [=] AMREX_GPU_DEVICE(const PType& p) noexcept + [=] AMREX_GPU_DEVICE(const ConstParticleTileDataType& ptd, const int i) noexcept { - amrex::ParticleReal const x = p.pos(0); - amrex::ParticleReal const y = p.pos(1); - amrex::ParticleReal const z = p.pos(2); - amrex::ParticleReal const w = p.rdata(T_RealSoAWeight); + + const amrex::ParticleReal x = ptd.rdata(0)[i]; + const amrex::ParticleReal y = ptd.rdata(1)[i]; + const amrex::ParticleReal z = ptd.rdata(2)[i]; + + const amrex::ParticleReal w = ptd.rdata(T_RealSoAWeight)[i]; return amrex::makeTuple(x, x*x, y, y*y, z, z*z, w); }, From 67c9aa2b9ff9f8a611b4a98ae1332e4794c4ac3a Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Mon, 29 Jan 2024 11:36:15 -0800 Subject: [PATCH 2/9] Remove particles that are initialized in the EB (#4585) --- Source/Particles/PhysicalParticleContainer.cpp | 16 ++++++++++++++++ Source/Particles/WarpXParticleContainer.cpp | 11 ++++++++++- Source/WarpX.H | 4 +++- 3 files changed, 29 insertions(+), 2 deletions(-) diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index fd55e369040..17c238e1121 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -39,6 +39,10 @@ #include "Utils/WarpXAlgorithmSelection.H" #include "Utils/WarpXConst.H" #include "Utils/WarpXProfilerWrapper.H" +#ifdef AMREX_USE_EB +# include "EmbeddedBoundary/ParticleBoundaryProcess.H" +# include "EmbeddedBoundary/ParticleScraper.H" +#endif #include "WarpX.H" #include @@ -1428,6 +1432,12 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int } } + // Remove particles that are inside the embedded boundaries +#ifdef AMREX_USE_EB + auto & distance_to_eb = WarpX::GetInstance().GetDistanceToEB(); + scrapeParticles( *this, amrex::GetVecOfConstPtrs(distance_to_eb), ParticleBoundaryProcess::Absorb()); +#endif + // The function that calls this is responsible for redistributing particles. } @@ -1920,6 +1930,12 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, } } + // Remove particles that are inside the embedded boundaries +#ifdef AMREX_USE_EB + auto & distance_to_eb = WarpX::GetInstance().GetDistanceToEB(); + scrapeParticles(tmp_pc, amrex::GetVecOfConstPtrs(distance_to_eb), ParticleBoundaryProcess::Absorb()); +#endif + // Redistribute the new particles that were added to the temporary container. // (This eliminates invalid particles, and makes sure that particles // are in the right tile.) diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 48344b20fef..0d927caacc8 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -60,7 +60,10 @@ #include #include #include - +#ifdef AMREX_USE_EB +# include "EmbeddedBoundary/ParticleBoundaryProcess.H" +# include "EmbeddedBoundary/ParticleScraper.H" +#endif #ifdef AMREX_USE_OMP # include @@ -293,6 +296,12 @@ WarpXParticleContainer::AddNParticles (int /*lev*/, long n, ); } + // Remove particles that are inside the embedded boundaries +#ifdef AMREX_USE_EB + auto & distance_to_eb = WarpX::GetInstance().GetDistanceToEB(); + scrapeParticles( *this, amrex::GetVecOfConstPtrs(distance_to_eb), ParticleBoundaryProcess::Absorb()); +#endif + Redistribute(); } diff --git a/Source/WarpX.H b/Source/WarpX.H index c9aa24d3973..e8c7ae79f7e 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -132,7 +132,9 @@ public: MacroscopicProperties& GetMacroscopicProperties () { return *m_macroscopic_properties; } HybridPICModel& GetHybridPICModel () { return *m_hybrid_pic_model; } MultiDiagnostics& GetMultiDiags () {return *multi_diags;} - +#ifdef AMREX_USE_EB + amrex::Vector >& GetDistanceToEB () {return m_distance_to_eb;} +#endif ParticleBoundaryBuffer& GetParticleBoundaryBuffer () { return *m_particle_boundary_buffer; } static void shiftMF (amrex::MultiFab& mf, const amrex::Geometry& geom, From c9df4485fc29e1ec5111906975d359efcace269a Mon Sep 17 00:00:00 2001 From: Roelof Groenewald <40245517+roelof-groenewald@users.noreply.github.com> Date: Mon, 29 Jan 2024 11:43:50 -0800 Subject: [PATCH 3/9] Update DSMC picmi class initialization function (#4646) --- Python/pywarpx/picmi.py | 2 +- Regression/WarpX-tests.ini | 36 ++++++++++++++++++------------------ 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index 857f857a554..89bd5af2eab 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -1575,7 +1575,7 @@ def __init__(self, name, species, scattering_processes, ndt=None, **kw): self.handle_init(kw) - def initialize_inputs(self): + def collision_initialize_inputs(self): collision = pywarpx.Collisions.newcollision(self.name) collision.type = 'dsmc' collision.species = [species.name for species in self.species] diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 75939c1fd12..7af5eb4b002 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -3079,24 +3079,24 @@ doVis = 0 compareParticles = 0 analysisRoutine = Examples/Tests/electrostatic_dirichlet_bc/analysis.py -# [Python_dsmc_1d] -# buildDir = . -# inputFile = Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py -# runtime_params = -# customRunCmd = python3 PICMI_inputs_1d.py --test --dsmc -# dim = 1 -# addToCompileString = USE_PYTHON_MAIN=TRUE USE_OPENPMD=TRUE QED=FALSE -# cmakeSetupOpts = -DWarpX_DIMS=1 -DWarpX_APP=OFF -DWarpX_PYTHON=ON -DWarpX_OPENPMD=ON -DWarpX_QED=OFF -# target = pip_install -# restartTest = 0 -# useMPI = 1 -# numprocs = 2 -# useOMP = 1 -# numthreads = 1 -# compileTest = 0 -# doVis = 0 -# compareParticles = 0 -# analysisRoutine = Examples/Physics_applications/capacitive_discharge/analysis_dsmc.py +[Python_dsmc_1d] +buildDir = . +inputFile = Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py +runtime_params = +customRunCmd = python3 PICMI_inputs_1d.py --test --dsmc +dim = 1 +addToCompileString = USE_PYTHON_MAIN=TRUE USE_OPENPMD=TRUE QED=FALSE +cmakeSetupOpts = -DWarpX_DIMS=1 -DWarpX_APP=OFF -DWarpX_PYTHON=ON -DWarpX_OPENPMD=ON -DWarpX_QED=OFF +target = pip_install +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 0 +analysisRoutine = Examples/Physics_applications/capacitive_discharge/analysis_dsmc.py [Python_ElectrostaticSphereEB] buildDir = . From 0284b5731ebba43dbe571b7c53c4692c0271f4c3 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Mon, 29 Jan 2024 16:03:09 -0800 Subject: [PATCH 4/9] AMReX: Weekly Update (#4647) --- .github/workflows/cuda.yml | 2 +- Regression/WarpX-GPU-tests.ini | 2 +- Regression/WarpX-tests.ini | 2 +- cmake/dependencies/AMReX.cmake | 2 +- run_test.sh | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/.github/workflows/cuda.yml b/.github/workflows/cuda.yml index 09f7f9e80fa..eca39369a08 100644 --- a/.github/workflows/cuda.yml +++ b/.github/workflows/cuda.yml @@ -115,7 +115,7 @@ jobs: which nvcc || echo "nvcc not in PATH!" git clone https://github.com/AMReX-Codes/amrex.git ../amrex - cd ../amrex && git checkout --detach 73b215557c0e842c3e829b683939bbb7a7e12373 && cd - + cd ../amrex && git checkout --detach ab99ea69089f4fffbfd727ed965d5ceb3d905baa && cd - make COMP=gcc QED=FALSE USE_MPI=TRUE USE_GPU=TRUE USE_OMP=FALSE USE_PSATD=TRUE USE_CCACHE=TRUE -j 2 ccache -s diff --git a/Regression/WarpX-GPU-tests.ini b/Regression/WarpX-GPU-tests.ini index f0557e05fc4..51ae03b60fa 100644 --- a/Regression/WarpX-GPU-tests.ini +++ b/Regression/WarpX-GPU-tests.ini @@ -60,7 +60,7 @@ emailBody = Check https://ccse.lbl.gov/pub/GpuRegressionTesting/WarpX/ for more [AMReX] dir = /home/regtester/git/amrex/ -branch = 73b215557c0e842c3e829b683939bbb7a7e12373 +branch = ab99ea69089f4fffbfd727ed965d5ceb3d905baa [source] dir = /home/regtester/git/WarpX diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 7af5eb4b002..ba5ae86c28e 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -59,7 +59,7 @@ emailBody = Check https://ccse.lbl.gov/pub/RegressionTesting/WarpX/ for more det [AMReX] dir = /home/regtester/AMReX_RegTesting/amrex/ -branch = 73b215557c0e842c3e829b683939bbb7a7e12373 +branch = ab99ea69089f4fffbfd727ed965d5ceb3d905baa [source] dir = /home/regtester/AMReX_RegTesting/warpx diff --git a/cmake/dependencies/AMReX.cmake b/cmake/dependencies/AMReX.cmake index 31c8919e38f..72173b2acb0 100644 --- a/cmake/dependencies/AMReX.cmake +++ b/cmake/dependencies/AMReX.cmake @@ -269,7 +269,7 @@ set(WarpX_amrex_src "" set(WarpX_amrex_repo "https://github.com/AMReX-Codes/amrex.git" CACHE STRING "Repository URI to pull and build AMReX from if(WarpX_amrex_internal)") -set(WarpX_amrex_branch "73b215557c0e842c3e829b683939bbb7a7e12373" +set(WarpX_amrex_branch "ab99ea69089f4fffbfd727ed965d5ceb3d905baa" CACHE STRING "Repository branch for WarpX_amrex_repo if(WarpX_amrex_internal)") diff --git a/run_test.sh b/run_test.sh index 77980bb00b0..20c6638f4d0 100755 --- a/run_test.sh +++ b/run_test.sh @@ -68,7 +68,7 @@ python3 -m pip install --upgrade -r warpx/Regression/requirements.txt # Clone AMReX and warpx-data git clone https://github.com/AMReX-Codes/amrex.git -cd amrex && git checkout --detach 73b215557c0e842c3e829b683939bbb7a7e12373 && cd - +cd amrex && git checkout --detach ab99ea69089f4fffbfd727ed965d5ceb3d905baa && cd - # warpx-data contains various required data sets git clone --depth 1 https://github.com/ECP-WarpX/warpx-data.git # openPMD-example-datasets contains various required data sets From 3e551e0a0b559114b91a341a5f5b7e0a85faeedd Mon Sep 17 00:00:00 2001 From: GivenChen <52523185+GivenChen@users.noreply.github.com> Date: Mon, 29 Jan 2024 22:18:33 -0500 Subject: [PATCH 5/9] Fix neutral particle got pushed (#4643) * Update PushSelector.H add the pointer of the array of ion_lev to the pusher function, and use this pointer to determine whether to change the charge instead of the value of ion_lev * Update PhysicalParticleContainer.cpp Change the variable in the targeted function * Alternative approach --------- Co-authored-by: Axel Huebl --- Source/Particles/PhysicalParticleContainer.cpp | 8 ++++---- Source/Particles/Pusher/PushSelector.H | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 17c238e1121..a8ed42ce419 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2880,7 +2880,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, doParticleMomentumPush<0>(ux[ip], uy[ip], uz[ip], Exp, Eyp, Ezp, Bxp, Byp, Bzp, - ion_lev ? ion_lev[ip] : 0, + ion_lev ? ion_lev[ip] : 1, m, q, pusher_algo, do_crr, #ifdef WARPX_QED t_chi_max, @@ -2900,7 +2900,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, doParticleMomentumPush<1>(ux[ip], uy[ip], uz[ip], Exp, Eyp, Ezp, Bxp, Byp, Bzp, - ion_lev ? ion_lev[ip] : 0, + ion_lev ? ion_lev[ip] : 1, m, q, pusher_algo, do_crr, t_chi_max, dt); @@ -3140,7 +3140,7 @@ PhysicalParticleContainer::ImplicitPushXP (WarpXParIter& pti, { doParticleMomentumPush<0>(ux[ip], uy[ip], uz[ip], Exp, Eyp, Ezp, Bxp, Byp, Bzp, - ion_lev ? ion_lev[ip] : 0, + ion_lev ? ion_lev[ip] : 1, m, q, pusher_algo, do_crr, #ifdef WARPX_QED t_chi_max, @@ -3152,7 +3152,7 @@ PhysicalParticleContainer::ImplicitPushXP (WarpXParIter& pti, if constexpr (qed_control == has_qed) { doParticleMomentumPush<1>(ux[ip], uy[ip], uz[ip], Exp, Eyp, Ezp, Bxp, Byp, Bzp, - ion_lev ? ion_lev[ip] : 0, + ion_lev ? ion_lev[ip] : 1, m, q, pusher_algo, do_crr, t_chi_max, dt); diff --git a/Source/Particles/Pusher/PushSelector.H b/Source/Particles/Pusher/PushSelector.H index 2a6ac8a7d6c..4a82e582bfb 100644 --- a/Source/Particles/Pusher/PushSelector.H +++ b/Source/Particles/Pusher/PushSelector.H @@ -57,7 +57,7 @@ void doParticleMomentumPush(amrex::ParticleReal& ux, const amrex::Real dt) { amrex::ParticleReal qp = a_q; - if (ion_lev) { qp *= ion_lev; } + qp *= ion_lev; if (do_crr) { #ifdef WARPX_QED From 33d7b496d11ca41ce9025c6df6feb56a2ecb8b30 Mon Sep 17 00:00:00 2001 From: Roelof Groenewald <40245517+roelof-groenewald@users.noreply.github.com> Date: Tue, 30 Jan 2024 06:50:50 -0800 Subject: [PATCH 6/9] reactivate CI test (Python_background_mcc_1d_tridiag) (#4650) --- Regression/WarpX-tests.ini | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index ba5ae86c28e..3489e8ad100 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -3021,24 +3021,24 @@ doVis = 0 compareParticles = 0 analysisRoutine = Examples/Physics_applications/capacitive_discharge/analysis_1d.py -# [Python_background_mcc_1d_tridiag] -# buildDir = . -# inputFile = Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py -# runtime_params = -# customRunCmd = python3 PICMI_inputs_1d.py --test -# dim = 1 -# addToCompileString = USE_PYTHON_MAIN=TRUE USE_OPENPMD=TRUE QED=FALSE -# cmakeSetupOpts = -DWarpX_DIMS=1 -DWarpX_APP=OFF -DWarpX_PYTHON=ON -DWarpX_OPENPMD=ON -DWarpX_QED=OFF -# target = pip_install -# restartTest = 0 -# useMPI = 1 -# numprocs = 2 -# useOMP = 1 -# numthreads = 1 -# compileTest = 0 -# doVis = 0 -# compareParticles = 0 -# analysisRoutine = Examples/Physics_applications/capacitive_discharge/analysis_1d.py +[Python_background_mcc_1d_tridiag] +buildDir = . +inputFile = Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py +runtime_params = +customRunCmd = python3 PICMI_inputs_1d.py --test +dim = 1 +addToCompileString = USE_PYTHON_MAIN=TRUE USE_OPENPMD=TRUE QED=FALSE +cmakeSetupOpts = -DWarpX_DIMS=1 -DWarpX_APP=OFF -DWarpX_PYTHON=ON -DWarpX_OPENPMD=ON -DWarpX_QED=OFF +target = pip_install +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 0 +analysisRoutine = Examples/Physics_applications/capacitive_discharge/analysis_1d.py [Python_collisionXZ] buildDir = . From 2708b4fb8fc16c847a49d5f7f367a9d0287d6042 Mon Sep 17 00:00:00 2001 From: Sarah Vickers <136610602+sevicky@users.noreply.github.com> Date: Tue, 30 Jan 2024 06:51:53 -0800 Subject: [PATCH 7/9] Grammar Change in Theory Intro (#4617) Split a run-on sentence into two sentences in the first paragraph of the introduction. --- Docs/source/theory/intro.rst | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Docs/source/theory/intro.rst b/Docs/source/theory/intro.rst index 10fbba679c2..2101d9d81cb 100644 --- a/Docs/source/theory/intro.rst +++ b/Docs/source/theory/intro.rst @@ -8,8 +8,7 @@ Introduction Plasma laser-driven (top) and charged-particles-driven (bottom) acceleration (rendering from 3-D Particle-In-Cell simulations). A laser beam (red and blue disks in top picture) or a charged particle beam (red dots in bottom picture) propagating (from left to right) through an under-dense plasma (not represented) displaces electrons, creating a plasma wakefield that supports very high electric fields (pale blue and yellow). These electric fields, which can be orders of magnitude larger than with conventional techniques, can be used to accelerate a short charged particle beam (white) to high-energy over a very short distance. -Computer simulations have had a profound impact on the design and understanding of past and present plasma acceleration experiments :cite:p:`i-Tsungpop06,i-Geddesjp08,i-Geddesscidac09,i-Geddespac09,i-Huangscidac09`, with -accurate modeling of wake formation, electron self-trapping and acceleration requiring fully kinetic methods (usually Particle-In-Cell) using large computational resources due to the wide range of space and time scales involved. Numerical modeling complements and guides the design and analysis of advanced accelerators, and can reduce development costs significantly. Despite the major recent experimental successes :cite:p:`i-LeemansPRL2014,i-Blumenfeld2007,i-BulanovSV2014,i-Steinke2016`, the various advanced acceleration concepts need significant progress to fulfill their potential. To this end, large-scale simulations will continue to be a key component toward reaching a detailed understanding of the complex interrelated physics phenomena at play. +Computer simulations have had a profound impact on the design and understanding of past and present plasma acceleration experiments :cite:p:`i-Tsungpop06,i-Geddesjp08,i-Geddesscidac09,i-Geddespac09,i-Huangscidac09`. Accurate modeling of wake formation, electron self-trapping and acceleration require fully kinetic methods (usually Particle-In-Cell) using large computational resources due to the wide range of space and time scales involved. Numerical modeling complements and guides the design and analysis of advanced accelerators, and can reduce development costs significantly. Despite the major recent experimental successes :cite:p:`i-LeemansPRL2014,i-Blumenfeld2007,i-BulanovSV2014,i-Steinke2016`, the various advanced acceleration concepts need significant progress to fulfill their potential. To this end, large-scale simulations will continue to be a key component toward reaching a detailed understanding of the complex interrelated physics phenomena at play. For such simulations, the most popular algorithm is the Particle-In-Cell (or PIC) technique, From 56e4a2a6e180c4f79ef1822493e11cc7561b59ed Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Tue, 30 Jan 2024 14:00:27 -0800 Subject: [PATCH 8/9] Revert "Update Particle Container to Pure SoA (#3850)" (#4652) This reverts commit 94ae11900131846e9f8b3704194673f7f02d8959. --- Docs/source/usage/workflows/python_extend.rst | 3 +- .../particle_data_python/PICMI_inputs_2d.py | 6 +- .../PICMI_inputs_prev_pos_2d.py | 6 +- .../PICMI_inputs_runtime_component_analyze.py | 6 +- Python/pywarpx/particle_containers.py | 255 ++++++++++-------- .../BackTransformParticleFunctor.H | 16 +- .../FlushFormats/FlushFormatAscent.cpp | 3 + .../FlushFormats/FlushFormatCheckpoint.cpp | 9 +- .../FlushFormats/FlushFormatPlotfile.cpp | 13 +- Source/Diagnostics/ParticleIO.cpp | 2 +- .../Diagnostics/ReducedDiags/FieldProbe.cpp | 6 +- .../FieldProbeParticleContainer.H | 20 +- .../FieldProbeParticleContainer.cpp | 36 ++- .../ReducedDiags/LoadBalanceCosts.cpp | 3 +- Source/Diagnostics/WarpXOpenPMD.H | 4 +- Source/Diagnostics/WarpXOpenPMD.cpp | 139 ++++++++-- .../ParticleBoundaryProcess.H | 5 +- Source/EmbeddedBoundary/ParticleScraper.H | 3 +- .../BinaryCollision/BinaryCollision.H | 16 +- .../Coulomb/PairWiseCoulombCollisionFunc.H | 7 +- .../Collision/BinaryCollision/DSMC/DSMC.H | 3 +- .../DSMC/SplitAndScatterFunc.H | 30 +-- .../NuclearFusion/NuclearFusionFunc.H | 14 +- .../ProtonBoronFusionInitializeMomentum.H | 10 +- .../TwoProductFusionInitializeMomentum.H | 4 +- .../BinaryCollision/ParticleCreationFunc.H | 60 ++--- .../BinaryCollision/ShuffleFisherYates.H | 2 +- .../Particles/Deposition/ChargeDeposition.H | 2 +- .../Particles/Deposition/CurrentDeposition.H | 2 +- .../ElementaryProcess/QEDPairGeneration.H | 2 +- .../ElementaryProcess/QEDPhotonEmission.H | 16 +- Source/Particles/LaserParticleContainer.H | 7 +- .../NamedComponentParticleContainer.H | 47 +--- Source/Particles/ParticleBoundaryBuffer.cpp | 4 +- Source/Particles/ParticleCreation/SmartCopy.H | 5 +- .../Particles/ParticleCreation/SmartCreate.H | 22 +- .../Particles/ParticleCreation/SmartUtils.H | 9 +- Source/Particles/PhysicalParticleContainer.H | 8 +- .../Particles/PhysicalParticleContainer.cpp | 129 +++++---- Source/Particles/Pusher/GetAndSetPosition.H | 152 ++++------- .../Particles/Resampling/LevelingThinning.cpp | 5 +- Source/Particles/Sorting/Partition.cpp | 4 +- Source/Particles/Sorting/SortingUtils.H | 41 +-- Source/Particles/Sorting/SortingUtils.cpp | 2 +- Source/Particles/WarpXParticleContainer.H | 14 +- Source/Particles/WarpXParticleContainer.cpp | 89 +++--- .../Particles/ParticleBoundaryBuffer.cpp | 10 +- .../PinnedMemoryParticleContainer.cpp | 2 +- .../Particles/WarpXParticleContainer.cpp | 18 +- Source/Utils/ParticleUtils.H | 7 +- Source/Utils/ParticleUtils.cpp | 26 +- Source/ablastr/particles/IndexHandling.H | 41 +++ Source/ablastr/particles/ParticleMoments.H | 25 +- 53 files changed, 709 insertions(+), 661 deletions(-) create mode 100644 Source/ablastr/particles/IndexHandling.H diff --git a/Docs/source/usage/workflows/python_extend.rst b/Docs/source/usage/workflows/python_extend.rst index 7a3ef0e7af9..6c9286c02ce 100644 --- a/Docs/source/usage/workflows/python_extend.rst +++ b/Docs/source/usage/workflows/python_extend.rst @@ -252,8 +252,7 @@ Particles can be added to the simulation at specific positions and with specific .. autoclass:: pywarpx.particle_containers.ParticleContainerWrapper :members: -The ``get_particle_real_arrays()``, ``get_particle_int_arrays()`` and -``get_particle_idcpu_arrays()`` functions are called +The ``get_particle_structs()`` and ``get_particle_arrays()`` functions are called by several utility functions of the form ``get_particle_{comp_name}`` where ``comp_name`` is one of ``x``, ``y``, ``z``, ``r``, ``theta``, ``id``, ``cpu``, ``weight``, ``ux``, ``uy`` or ``uz``. diff --git a/Examples/Tests/particle_data_python/PICMI_inputs_2d.py b/Examples/Tests/particle_data_python/PICMI_inputs_2d.py index 572871b8ed5..a4b7d9e134e 100755 --- a/Examples/Tests/particle_data_python/PICMI_inputs_2d.py +++ b/Examples/Tests/particle_data_python/PICMI_inputs_2d.py @@ -153,10 +153,10 @@ def add_particles(): ########################## assert (elec_wrapper.nps == 270 / (2 - args.unique)) -assert (elec_wrapper.particle_container.get_comp_index('w') == 2) -assert (elec_wrapper.particle_container.get_comp_index('newPid') == 6) +assert (elec_wrapper.particle_container.get_comp_index('w') == 0) +assert (elec_wrapper.particle_container.get_comp_index('newPid') == 4) -new_pid_vals = elec_wrapper.get_particle_real_arrays('newPid', 0) +new_pid_vals = elec_wrapper.get_particle_arrays('newPid', 0) for vals in new_pid_vals: assert np.allclose(vals, 5) diff --git a/Examples/Tests/particle_data_python/PICMI_inputs_prev_pos_2d.py b/Examples/Tests/particle_data_python/PICMI_inputs_prev_pos_2d.py index 5de9879f0f8..1becd4464e7 100755 --- a/Examples/Tests/particle_data_python/PICMI_inputs_prev_pos_2d.py +++ b/Examples/Tests/particle_data_python/PICMI_inputs_prev_pos_2d.py @@ -120,12 +120,12 @@ elec_count = elec_wrapper.nps # check that the runtime attributes have the right indices -assert (elec_wrapper.particle_container.get_comp_index('prev_x') == 6) -assert (elec_wrapper.particle_container.get_comp_index('prev_z') == 7) +assert (elec_wrapper.particle_container.get_comp_index('prev_x') == 4) +assert (elec_wrapper.particle_container.get_comp_index('prev_z') == 5) # sanity check that the prev_z values are reasonable and # that the correct number of values are returned -prev_z_vals = elec_wrapper.get_particle_real_arrays('prev_z', 0) +prev_z_vals = elec_wrapper.get_particle_arrays('prev_z', 0) running_count = 0 for z_vals in prev_z_vals: diff --git a/Examples/Tests/restart/PICMI_inputs_runtime_component_analyze.py b/Examples/Tests/restart/PICMI_inputs_runtime_component_analyze.py index 32c9f4e5808..706dedb6959 100755 --- a/Examples/Tests/restart/PICMI_inputs_runtime_component_analyze.py +++ b/Examples/Tests/restart/PICMI_inputs_runtime_component_analyze.py @@ -158,10 +158,10 @@ def add_particles(): ########################## assert electron_wrapper.nps == 90 -assert electron_wrapper.particle_container.get_comp_index("w") == 2 -assert electron_wrapper.particle_container.get_comp_index("newPid") == 6 +assert electron_wrapper.particle_container.get_comp_index("w") == 0 +assert electron_wrapper.particle_container.get_comp_index("newPid") == 4 -new_pid_vals = electron_wrapper.get_particle_real_arrays("newPid", 0) +new_pid_vals = electron_wrapper.get_particle_arrays("newPid", 0) for vals in new_pid_vals: assert np.allclose(vals, 5) diff --git a/Python/pywarpx/particle_containers.py b/Python/pywarpx/particle_containers.py index 72a675ec43c..273a981f4bd 100644 --- a/Python/pywarpx/particle_containers.py +++ b/Python/pywarpx/particle_containers.py @@ -117,10 +117,8 @@ def add_particles(self, x=None, y=None, z=None, ux=None, uy=None, kwargs[key] = np.full(maxlen, val) # --- The number of built in attributes - # --- The positions - built_in_attrs = libwarpx.dim # --- The three velocities - built_in_attrs += 3 + built_in_attrs = 3 if libwarpx.geometry_dim == 'rz': # --- With RZ, there is also theta built_in_attrs += 1 @@ -190,25 +188,28 @@ def add_real_comp(self, pid_name, comm=True): self.particle_container.add_real_comp(pid_name, comm) - def get_particle_real_arrays(self, comp_name, level, copy_to_host=False): + def get_particle_structs(self, level, copy_to_host=False): ''' - This returns a list of numpy or cupy arrays containing the particle real array data - on each tile for this process. + This returns a list of numpy or cupy arrays containing the particle struct data + on each tile for this process. The particle data is represented as a structured + array and contains the particle 'x', 'y', 'z', and 'idcpu'. - Unless copy_to_host is specified, the data for the arrays are not - copied, but share the underlying memory buffer with WarpX. The + Unless copy_to_host is specified, the data for the structs are + not copied, but share the underlying memory buffer with WarpX. The arrays are fully writeable. + Note that cupy does not support structs: + https://github.com/cupy/cupy/issues/2031 + and will return arrays of binary blobs for the AoS (DP: ``"|V24"``). If copied + to host via copy_to_host, we correct for the right numpy AoS type. + Parameters ---------- - comp_name : str - The component of the array data that will be returned - - level : int + level : int The refinement level to reference (default=0) - copy_to_host : bool + copy_to_host : bool For GPU-enabled runs, one can either return the GPU arrays (the default) or force a device-to-host copy. @@ -216,29 +217,30 @@ def get_particle_real_arrays(self, comp_name, level, copy_to_host=False): ------- List of arrays - The requested particle array data + The requested particle struct data ''' - comp_idx = self.particle_container.get_comp_index(comp_name) - - data_array = [] + particle_data = [] for pti in libwarpx.libwarpx_so.WarpXParIter(self.particle_container, level): - soa = pti.soa() - idx = soa.GetRealData(comp_idx) if copy_to_host: - data_array.append(idx.to_numpy(copy=True)) + particle_data.append(pti.aos().to_numpy(copy=True)) else: + if libwarpx.amr.Config.have_gpu: + libwarpx.amr.Print( + "get_particle_structs: cupy does not yet support structs. " + "https://github.com/cupy/cupy/issues/2031" + "Did you mean copy_to_host=True?" + ) xp, cupy_status = load_cupy() if cupy_status is not None: libwarpx.amr.Print(cupy_status) - - data_array.append(xp.array(idx, copy=False)) - - return data_array + aos_arr = xp.array(pti.aos(), copy=False) # void blobs for cupy + particle_data.append(aos_arr) + return particle_data - def get_particle_int_arrays(self, comp_name, level, copy_to_host=False): + def get_particle_arrays(self, comp_name, level, copy_to_host=False): ''' - This returns a list of numpy or cupy arrays containing the particle int array data + This returns a list of numpy or cupy arrays containing the particle array data on each tile for this process. Unless copy_to_host is specified, the data for the arrays are not @@ -264,52 +266,12 @@ def get_particle_int_arrays(self, comp_name, level, copy_to_host=False): List of arrays The requested particle array data ''' - comp_idx = self.particle_container.get_icomp_index(comp_name) - - data_array = [] - for pti in libwarpx.libwarpx_so.WarpXParIter(self.particle_container, level): - soa = pti.soa() - idx = soa.GetIntData(comp_idx) - if copy_to_host: - data_array.append(idx.to_numpy(copy=True)) - else: - xp, cupy_status = load_cupy() - if cupy_status is not None: - libwarpx.amr.Print(cupy_status) - - data_array.append(xp.array(idx, copy=False)) - - return data_array - - - def get_particle_idcpu_arrays(self, level, copy_to_host=False): - ''' - This returns a list of numpy or cupy arrays containing the particle idcpu data - on each tile for this process. - - Unless copy_to_host is specified, the data for the arrays are not - copied, but share the underlying memory buffer with WarpX. The - arrays are fully writeable. - - Parameters - ---------- - level : int - The refinement level to reference (default=0) - - copy_to_host : bool - For GPU-enabled runs, one can either return the GPU - arrays (the default) or force a device-to-host copy. - - Returns - ------- + comp_idx = self.particle_container.get_comp_index(comp_name) - List of arrays - The requested particle array data - ''' data_array = [] for pti in libwarpx.libwarpx_so.WarpXParIter(self.particle_container, level): soa = pti.soa() - idx = soa.GetIdCPUData() + idx = soa.GetRealData(comp_idx) if copy_to_host: data_array.append(idx.to_numpy(copy=True)) else: @@ -322,31 +284,6 @@ def get_particle_idcpu_arrays(self, level, copy_to_host=False): return data_array - def get_particle_idcpu(self, level=0, copy_to_host=False): - ''' - Return a list of numpy or cupy arrays containing the particle 'idcpu' - numbers on each tile. - - Parameters - ---------- - - level : int - The refinement level to reference (default=0) - - copy_to_host : bool - For GPU-enabled runs, one can either return the GPU - arrays (the default) or force a device-to-host copy. - - Returns - ------- - - List of arrays - The requested particle idcpu - ''' - return self.get_particle_idcpu_arrays(level, copy_to_host=copy_to_host) - idcpu = property(get_particle_idcpu) - - def get_particle_id(self, level=0, copy_to_host=False): ''' Return a list of numpy or cupy arrays containing the particle 'id' @@ -368,8 +305,8 @@ def get_particle_id(self, level=0, copy_to_host=False): List of arrays The requested particle ids ''' - idcpu = self.get_particle_idcpu(level, copy_to_host) - return [libwarpx.amr.unpack_ids(tile) for tile in idcpu] + structs = self.get_particle_structs(level, copy_to_host) + return [libwarpx.amr.unpack_ids(struct['cpuid']) for struct in structs] def get_particle_cpu(self, level=0, copy_to_host=False): @@ -393,8 +330,8 @@ def get_particle_cpu(self, level=0, copy_to_host=False): List of arrays The requested particle cpus ''' - idcpu = self.get_particle_idcpu(level, copy_to_host) - return [libwarpx.amr.unpack_cpus(tile) for tile in idcpu] + structs = self.get_particle_structs(level, copy_to_host) + return [libwarpx.amr.unpack_cpus(struct['cpuid']) for struct in structs] def get_particle_x(self, level=0, copy_to_host=False): @@ -418,7 +355,20 @@ def get_particle_x(self, level=0, copy_to_host=False): List of arrays The requested particle x position ''' - return self.get_particle_real_arrays('x', level, copy_to_host=copy_to_host) + if copy_to_host: + # Use the numpy version of cosine + xp = np + else: + xp, cupy_status = load_cupy() + + structs = self.get_particle_structs(level, copy_to_host) + if libwarpx.geometry_dim == '3d' or libwarpx.geometry_dim == '2d': + return [struct['x'] for struct in structs] + elif libwarpx.geometry_dim == 'rz': + theta = self.get_particle_theta(level, copy_to_host) + return [struct['x']*xp.cos(theta) for struct, theta in zip(structs, theta)] + elif libwarpx.geometry_dim == '1d': + raise Exception('get_particle_x: There is no x coordinate with 1D Cartesian') xp = property(get_particle_x) @@ -443,7 +393,20 @@ def get_particle_y(self, level=0, copy_to_host=False): List of arrays The requested particle y position ''' - return self.get_particle_real_arrays('y', level, copy_to_host=copy_to_host) + if copy_to_host: + # Use the numpy version of sine + xp = np + else: + xp, cupy_status = load_cupy() + + structs = self.get_particle_structs(level, copy_to_host) + if libwarpx.geometry_dim == '3d': + return [struct['y'] for struct in structs] + elif libwarpx.geometry_dim == 'rz': + theta = self.get_particle_theta(level, copy_to_host) + return [struct['x']*xp.sin(theta) for struct, theta in zip(structs, theta)] + elif libwarpx.geometry_dim == '1d' or libwarpx.geometry_dim == '2d': + raise Exception('get_particle_y: There is no y coordinate with 1D or 2D Cartesian') yp = property(get_particle_y) @@ -470,12 +433,11 @@ def get_particle_r(self, level=0, copy_to_host=False): ''' xp, cupy_status = load_cupy() + structs = self.get_particle_structs(level, copy_to_host) if libwarpx.geometry_dim == 'rz': - return self.get_particle_x(level, copy_to_host) + return [struct['x'] for struct in structs] elif libwarpx.geometry_dim == '3d': - x = self.get_particle_x(level, copy_to_host) - y = self.get_particle_y(level, copy_to_host) - return xp.sqrt(x**2 + y**2) + return [xp.sqrt(struct['x']**2 + struct['y']**2) for struct in structs] elif libwarpx.geometry_dim == '2d' or libwarpx.geometry_dim == '1d': raise Exception('get_particle_r: There is no r coordinate with 1D or 2D Cartesian') rp = property(get_particle_r) @@ -505,11 +467,10 @@ def get_particle_theta(self, level=0, copy_to_host=False): xp, cupy_status = load_cupy() if libwarpx.geometry_dim == 'rz': - return self.get_particle_real_arrays('theta', level, copy_to_host) + return self.get_particle_arrays('theta', level, copy_to_host) elif libwarpx.geometry_dim == '3d': - x = self.get_particle_x(level, copy_to_host) - y = self.get_particle_y(level, copy_to_host) - return xp.arctan2(y, x) + structs = self.get_particle_structs(level, copy_to_host) + return [xp.arctan2(struct['y'], struct['x']) for struct in structs] elif libwarpx.geometry_dim == '2d' or libwarpx.geometry_dim == '1d': raise Exception('get_particle_theta: There is no theta coordinate with 1D or 2D Cartesian') thetap = property(get_particle_theta) @@ -536,7 +497,13 @@ def get_particle_z(self, level=0, copy_to_host=False): List of arrays The requested particle z position ''' - return self.get_particle_real_arrays('z', level, copy_to_host=copy_to_host) + structs = self.get_particle_structs(level, copy_to_host) + if libwarpx.geometry_dim == '3d': + return [struct['z'] for struct in structs] + elif libwarpx.geometry_dim == 'rz' or libwarpx.geometry_dim == '2d': + return [struct['y'] for struct in structs] + elif libwarpx.geometry_dim == '1d': + return [struct['x'] for struct in structs] zp = property(get_particle_z) @@ -561,7 +528,7 @@ def get_particle_weight(self, level=0, copy_to_host=False): List of arrays The requested particle weight ''' - return self.get_particle_real_arrays('w', level, copy_to_host=copy_to_host) + return self.get_particle_arrays('w', level, copy_to_host=copy_to_host) wp = property(get_particle_weight) @@ -586,7 +553,7 @@ def get_particle_ux(self, level=0, copy_to_host=False): List of arrays The requested particle x momentum ''' - return self.get_particle_real_arrays('ux', level, copy_to_host=copy_to_host) + return self.get_particle_arrays('ux', level, copy_to_host=copy_to_host) uxp = property(get_particle_ux) @@ -611,7 +578,7 @@ def get_particle_uy(self, level=0, copy_to_host=False): List of arrays The requested particle y momentum ''' - return self.get_particle_real_arrays('uy', level, copy_to_host=copy_to_host) + return self.get_particle_arrays('uy', level, copy_to_host=copy_to_host) uyp = property(get_particle_uy) @@ -637,7 +604,7 @@ def get_particle_uz(self, level=0, copy_to_host=False): The requested particle z momentum ''' - return self.get_particle_real_arrays('uz', level, copy_to_host=copy_to_host) + return self.get_particle_arrays('uz', level, copy_to_host=copy_to_host) uzp = property(get_particle_uz) @@ -753,6 +720,70 @@ def get_particle_boundary_buffer_size(self, species_name, boundary, local=False) ) + def get_particle_boundary_buffer_structs( + self, species_name, boundary, level, copy_to_host=False + ): + ''' + This returns a list of numpy or cupy arrays containing the particle struct data + for a species that has been scraped by a specific simulation boundary. The + particle data is represented as a structured array and contains the + particle 'x', 'y', 'z', and 'idcpu'. + + Unless copy_to_host is specified, the data for the structs are + not copied, but share the underlying memory buffer with WarpX. The + arrays are fully writeable. + + Note that cupy does not support structs: + https://github.com/cupy/cupy/issues/2031 + and will return arrays of binary blobs for the AoS (DP: ``"|V24"``). If copied + to host via copy_to_host, we correct for the right numpy AoS type. + + Parameters + ---------- + + species_name : str + The species name that the data will be returned for + + boundary : str + The boundary from which to get the scraped particle data in the + form x/y/z_hi/lo or eb. + + level : int + The refinement level to reference (default=0) + + copy_to_host : bool + For GPU-enabled runs, one can either return the GPU + arrays (the default) or force a device-to-host copy. + + Returns + ------- + + List of arrays + The requested particle struct data + ''' + particle_container = self.particle_buffer.get_particle_container( + species_name, self._get_boundary_number(boundary) + ) + + particle_data = [] + for pti in libwarpx.libwarpx_so.BoundaryBufferParIter(particle_container, level): + if copy_to_host: + particle_data.append(pti.aos().to_numpy(copy=True)) + else: + if libwarpx.amr.Config.have_gpu: + libwarpx.amr.Print( + "get_particle_structs: cupy does not yet support structs. " + "https://github.com/cupy/cupy/issues/2031" + "Did you mean copy_to_host=True?" + ) + xp, cupy_status = load_cupy() + if cupy_status is not None: + libwarpx.amr.Print(cupy_status) + aos_arr = xp.array(pti.aos(), copy=False) # void blobs for cupy + particle_data.append(aos_arr) + return particle_data + + def get_particle_boundary_buffer(self, species_name, boundary, comp_name, level): ''' This returns a list of numpy or cupy arrays containing the particle array data diff --git a/Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.H b/Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.H index fa00f6288f9..e88b522f148 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.H +++ b/Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.H @@ -143,19 +143,19 @@ struct LorentzTransformParticles const amrex::ParticleReal uzp = uz_old_p * weight_old + uz_new_p * weight_new; #if defined (WARPX_DIM_3D) - dst.m_rdata[PIdx::x][i_dst] = xp; - dst.m_rdata[PIdx::y][i_dst] = yp; - dst.m_rdata[PIdx::z][i_dst] = zp; + dst.m_aos[i_dst].pos(0) = xp; + dst.m_aos[i_dst].pos(1) = yp; + dst.m_aos[i_dst].pos(2) = zp; #elif defined (WARPX_DIM_RZ) - dst.m_rdata[PIdx::x][i_dst] = std::sqrt(xp*xp + yp*yp); - dst.m_rdata[PIdx::z][i_dst] = zp; + dst.m_aos[i_dst].pos(0) = std::sqrt(xp*xp + yp*yp); + dst.m_aos[i_dst].pos(1) = zp; dst.m_rdata[PIdx::theta][i_dst] = std::atan2(yp, xp); #elif defined (WARPX_DIM_XZ) - dst.m_rdata[PIdx::x][i_dst] = xp; - dst.m_rdata[PIdx::z][i_dst] = zp; + dst.m_aos[i_dst].pos(0) = xp; + dst.m_aos[i_dst].pos(1) = zp; amrex::ignore_unused(yp); #elif defined (WARPX_DIM_1D_Z) - dst.m_rdata[PIdx::z][i_dst] = zp; + dst.m_aos[i_dst].pos(0) = zp; amrex::ignore_unused(xp, yp); #else amrex::ignore_unused(xp, yp, zp); diff --git a/Source/Diagnostics/FlushFormats/FlushFormatAscent.cpp b/Source/Diagnostics/FlushFormats/FlushFormatAscent.cpp index c224bc4f871..980047e3b46 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatAscent.cpp +++ b/Source/Diagnostics/FlushFormats/FlushFormatAscent.cpp @@ -94,6 +94,9 @@ FlushFormatAscent::WriteParticles(const amrex::Vector& particle_di // get names of real comps std::map real_comps_map = pc->getParticleComps(); + // WarpXParticleContainer compile-time extra AoS attributes (Real): 0 + // WarpXParticleContainer compile-time extra AoS attributes (int): 0 + // WarpXParticleContainer compile-time extra SoA attributes (Real): PIdx::nattribs // not an efficient search, but N is small... for(int j = 0; j < PIdx::nattribs; ++j) diff --git a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp index 00c28eead51..5f59cd723da 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp +++ b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp @@ -178,8 +178,8 @@ FlushFormatCheckpoint::CheckpointParticles ( Vector real_names; Vector int_names; - // note: positions skipped here, since we reconstruct a plotfile SoA from them real_names.push_back("weight"); + real_names.push_back("momentum_x"); real_names.push_back("momentum_y"); real_names.push_back("momentum_z"); @@ -189,12 +189,9 @@ FlushFormatCheckpoint::CheckpointParticles ( #endif // get the names of the real comps - // note: skips the mandatory AMREX_SPACEDIM positions for pure SoA - real_names.resize(pc->NumRealComps() - AMREX_SPACEDIM); + real_names.resize(pc->NumRealComps()); auto runtime_rnames = pc->getParticleRuntimeComps(); - for (auto const& x : runtime_rnames) { - real_names[x.second + PIdx::nattribs - AMREX_SPACEDIM] = x.first; - } + for (auto const& x : runtime_rnames) { real_names[x.second+PIdx::nattribs] = x.first; } // and the int comps int_names.resize(pc->NumIntComps()); diff --git a/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp b/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp index 5ec067a6eac..df73ed34c94 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp +++ b/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp @@ -355,8 +355,8 @@ FlushFormatPlotfile::WriteParticles(const std::string& dir, Vector int_flags; Vector real_flags; - // note: positions skipped here, since we reconstruct a plotfile SoA from them real_names.push_back("weight"); + real_names.push_back("momentum_x"); real_names.push_back("momentum_y"); real_names.push_back("momentum_z"); @@ -366,21 +366,14 @@ FlushFormatPlotfile::WriteParticles(const std::string& dir, #endif // get the names of the real comps - - // note: skips the mandatory AMREX_SPACEDIM positions for pure SoA - real_names.resize(tmp.NumRealComps() - AMREX_SPACEDIM); + real_names.resize(tmp.NumRealComps()); auto runtime_rnames = tmp.getParticleRuntimeComps(); - for (auto const& x : runtime_rnames) { - real_names[x.second + PIdx::nattribs - AMREX_SPACEDIM] = x.first; - } + for (auto const& x : runtime_rnames) { real_names[x.second+PIdx::nattribs] = x.first; } // plot any "extra" fields by default real_flags = part_diag.m_plot_flags; real_flags.resize(tmp.NumRealComps(), 1); - // note: skip the mandatory AMREX_SPACEDIM positions for pure SoA - real_flags.erase(real_flags.begin(), real_flags.begin() + AMREX_SPACEDIM); - // and the names int_names.resize(tmp.NumIntComps()); auto runtime_inames = tmp.getParticleRuntimeiComps(); diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index a8bb9303fe1..7ca5e6541d7 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -160,7 +160,7 @@ MultiParticleContainer::Restart (const std::string& dir) ); } - for (int j = PIdx::nattribs-AMREX_SPACEDIM; j < nr; ++j) { + for (int j = PIdx::nattribs; j < nr; ++j) { const auto& comp_name = real_comp_names[j]; auto current_comp_names = pc->getParticleComps(); auto search = current_comp_names.find(comp_name); diff --git a/Source/Diagnostics/ReducedDiags/FieldProbe.cpp b/Source/Diagnostics/ReducedDiags/FieldProbe.cpp index 24ad0e64ea8..9f45392bb0a 100644 --- a/Source/Diagnostics/ReducedDiags/FieldProbe.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldProbe.cpp @@ -431,6 +431,8 @@ void FieldProbe::ComputeDiags (int step) { const auto getPosition = GetParticlePosition(pti); auto setPosition = SetParticlePosition(pti); + const auto& aos = pti.GetArrayOfStructs(); + const auto* AMREX_RESTRICT m_structs = aos().dataPtr(); auto const np = pti.numParticles(); if (update_particles_moving_window) @@ -480,8 +482,6 @@ void FieldProbe::ComputeDiags (int step) ParticleReal* const AMREX_RESTRICT part_Bz = attribs[FieldProbePIdx::Bz].dataPtr(); ParticleReal* const AMREX_RESTRICT part_S = attribs[FieldProbePIdx::S].dataPtr(); - auto * const AMREX_RESTRICT idcpu = pti.GetStructOfArrays().GetIdCPUData().data(); - const auto &xyzmin = WarpX::LowerCorner(box, lev, 0._rt); const std::array &dx = WarpX::CellSize(lev); @@ -556,7 +556,7 @@ void FieldProbe::ComputeDiags (int step) amrex::ParticleReal xp, yp, zp; getPosition(ip, xp, yp, zp); long idx = ip*noutputs; - dvp[idx++] = amrex::ParticleIDWrapper{idcpu[ip]}; // all particles created on IO cpu + dvp[idx++] = m_structs[ip].id(); dvp[idx++] = xp; dvp[idx++] = yp; dvp[idx++] = zp; diff --git a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H index 7d59ade5dc6..c85bf8fd541 100644 --- a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H +++ b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H @@ -24,14 +24,7 @@ struct FieldProbePIdx { enum { -#if !defined (WARPX_DIM_1D_Z) - x, -#endif -#if defined (WARPX_DIM_3D) - y, -#endif - z, - Ex, Ey, Ez, + Ex = 0, Ey, Ez, Bx, By, Bz, S, //!< the Poynting vector #ifdef WARPX_DIM_RZ @@ -47,14 +40,9 @@ struct FieldProbePIdx * nattribs tells the particle container to allot 7 SOA values. */ class FieldProbeParticleContainer - : public amrex::ParticleContainerPureSoA + : public amrex::ParticleContainer<0, 0, FieldProbePIdx::nattribs> { public: - static constexpr int NStructReal = 0; - static constexpr int NStructInt = 0; - static constexpr int NReal = FieldProbePIdx::nattribs; - static constexpr int NInt = 0; - FieldProbeParticleContainer (amrex::AmrCore* amr_core); ~FieldProbeParticleContainer() override = default; @@ -64,9 +52,9 @@ public: FieldProbeParticleContainer& operator= ( FieldProbeParticleContainer&& ) = default; //! amrex iterator for our number of attributes - using iterator = amrex::ParIterSoA; + using iterator = amrex::ParIter<0, 0, FieldProbePIdx::nattribs, 0>; //! amrex iterator for our number of attributes (read-only) - using const_iterator = amrex::ParConstIterSoA; + using const_iterator = amrex::ParConstIter<0, 0, FieldProbePIdx::nattribs, 0>; //! similar to WarpXParticleContainer::AddNParticles but does not include u(x,y,z) void AddNParticles (int lev, amrex::Vector const & x, amrex::Vector const & y, amrex::Vector const & z); diff --git a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp index 18137fe9b2c..1fd741ddc47 100644 --- a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp @@ -59,7 +59,7 @@ using namespace amrex; FieldProbeParticleContainer::FieldProbeParticleContainer (AmrCore* amr_core) - : ParticleContainerPureSoA(amr_core->GetParGDB()) + : ParticleContainer<0, 0, FieldProbePIdx::nattribs>(amr_core->GetParGDB()) { SetParticleSize(); } @@ -89,17 +89,33 @@ FieldProbeParticleContainer::AddNParticles (int lev, * is then coppied to the permament tile which is stored on the particle * (particle_tile). */ - using PinnedTile = typename ContainerLike::ParticleTileType; + using PinnedTile = ParticleTile, + NArrayReal, NArrayInt, + amrex::PinnedArenaAllocator>; PinnedTile pinned_tile; pinned_tile.define(NumRuntimeRealComps(), NumRuntimeIntComps()); for (int i = 0; i < np; i++) { - auto & idcpu_data = pinned_tile.GetStructOfArrays().GetIdCPUData(); - idcpu_data.push_back(0); - amrex::ParticleIDWrapper{idcpu_data.back()} = ParticleType::NextID(); - amrex::ParticleCPUWrapper(idcpu_data.back()) = ParallelDescriptor::MyProc(); + ParticleType p; + p.id() = ParticleType::NextID(); + p.cpu() = ParallelDescriptor::MyProc(); +#if defined(WARPX_DIM_3D) + p.pos(0) = x[i]; + p.pos(1) = y[i]; + p.pos(2) = z[i]; +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + amrex::ignore_unused(y); + p.pos(0) = x[i]; + p.pos(1) = z[i]; +#elif defined(WARPX_DIM_1D_Z) + amrex::ignore_unused(x, y); + p.pos(0) = z[i]; +#endif + + // write position, cpu id, and particle id to particle + pinned_tile.push_back(p); } // write Real attributes (SoA) to particle initialized zero @@ -109,13 +125,7 @@ FieldProbeParticleContainer::AddNParticles (int lev, #ifdef WARPX_DIM_RZ pinned_tile.push_back_real(FieldProbePIdx::theta, np, 0.0); #endif -#if !defined (WARPX_DIM_1D_Z) - pinned_tile.push_back_real(FieldProbePIdx::x, x); -#endif -#if defined (WARPX_DIM_3D) - pinned_tile.push_back_real(FieldProbePIdx::y, y); -#endif - pinned_tile.push_back_real(FieldProbePIdx::z, z); + pinned_tile.push_back_real(FieldProbePIdx::Ex, np, 0.0); pinned_tile.push_back_real(FieldProbePIdx::Ey, np, 0.0); pinned_tile.push_back_real(FieldProbePIdx::Ez, np, 0.0); diff --git a/Source/Diagnostics/ReducedDiags/LoadBalanceCosts.cpp b/Source/Diagnostics/ReducedDiags/LoadBalanceCosts.cpp index b4e07b51982..893b00a5f00 100644 --- a/Source/Diagnostics/ReducedDiags/LoadBalanceCosts.cpp +++ b/Source/Diagnostics/ReducedDiags/LoadBalanceCosts.cpp @@ -56,7 +56,8 @@ namespace auto const & plev = pc.GetParticles(lev); auto const & ptile = plev.at(box_index); - auto const np = ptile.numParticles(); + auto const & aos = ptile.GetArrayOfStructs(); + auto const np = aos.numParticles(); num_macro_particles += np; } diff --git a/Source/Diagnostics/WarpXOpenPMD.H b/Source/Diagnostics/WarpXOpenPMD.H index 110f1ada649..e3b7b893d0a 100644 --- a/Source/Diagnostics/WarpXOpenPMD.H +++ b/Source/Diagnostics/WarpXOpenPMD.H @@ -41,7 +41,7 @@ class WarpXParticleCounter { public: using ParticleContainer = typename WarpXParticleContainer::ContainerLike; - using ParticleIter = typename amrex::ParIterSoA; + using ParticleIter = typename amrex::ParIter<0, 0, PIdx::nattribs, 0, amrex::PinnedArenaAllocator>; WarpXParticleCounter (ParticleContainer* pc); [[nodiscard]] unsigned long GetTotalNumParticles () const {return m_Total;} @@ -77,7 +77,7 @@ class WarpXOpenPMDPlot { public: using ParticleContainer = typename WarpXParticleContainer::ContainerLike; - using ParticleIter = typename amrex::ParConstIterSoA; + using ParticleIter = typename amrex::ParConstIter<0, 0, PIdx::nattribs, 0, amrex::PinnedArenaAllocator>; /** Initialize openPMD I/O routines * diff --git a/Source/Diagnostics/WarpXOpenPMD.cpp b/Source/Diagnostics/WarpXOpenPMD.cpp index d9b70a2e18d..71d96a47927 100644 --- a/Source/Diagnostics/WarpXOpenPMD.cpp +++ b/Source/Diagnostics/WarpXOpenPMD.cpp @@ -18,9 +18,11 @@ #include "WarpX.H" #include "OpenPMDHelpFunction.H" +#include #include #include +#include #include #include #include @@ -546,13 +548,6 @@ for (unsigned i = 0, n = particle_diags.size(); i < n; ++i) { // see openPMD ED-PIC extension for namings // note: an underscore separates the record name from its component // for non-scalar records -#if !defined (WARPX_DIM_1D_Z) - real_names.push_back("position_x"); -#endif -#if defined (WARPX_DIM_3D) - real_names.push_back("position_y"); -#endif - real_names.push_back("position_z"); real_names.push_back("weighting"); real_names.push_back("momentum_x"); real_names.push_back("momentum_y"); @@ -737,7 +732,77 @@ WarpXOpenPMDPlot::DumpToFile (ParticleContainer* pc, contributed_particles = true; - // save particle properties + // get position and particle ID from aos + // note: this implementation iterates the AoS 4x... + // if we flush late as we do now, we can also copy out the data in one go + const auto &aos = pti.GetArrayOfStructs(); // size = numParticlesOnTile + { + // Save positions +#if defined(WARPX_DIM_RZ) + { + const std::shared_ptr z( + new amrex::ParticleReal[numParticleOnTile], + [](amrex::ParticleReal const *p) { delete[] p; } + ); + for (auto i = 0; i < numParticleOnTile; i++) { + z.get()[i] = aos[i].pos(1); // {0: "r", 1: "z"} + } + std::string const positionComponent = "z"; + currSpecies["position"]["z"].storeChunk(z, {offset}, {numParticleOnTile64}); + } + + // reconstruct x and y from polar coordinates r, theta + auto const& soa = pti.GetStructOfArrays(); + amrex::ParticleReal const* theta = soa.GetRealData(PIdx::theta).dataPtr(); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(theta != nullptr, "openPMD: invalid theta pointer."); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(int(soa.GetRealData(PIdx::theta).size()) == numParticleOnTile, + "openPMD: theta and tile size do not match"); + { + const std::shared_ptr< amrex::ParticleReal > x( + new amrex::ParticleReal[numParticleOnTile], + [](amrex::ParticleReal const *p){ delete[] p; } + ); + const std::shared_ptr< amrex::ParticleReal > y( + new amrex::ParticleReal[numParticleOnTile], + [](amrex::ParticleReal const *p){ delete[] p; } + ); + for (auto i=0; i 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]; + currSpecies["position"][positionComponent].storeChunk(curr, {offset}, + {numParticleOnTile64}); + } +#endif + + // save particle ID after converting it to a globally unique ID + const std::shared_ptr ids( + new uint64_t[numParticleOnTile], + [](uint64_t const *p) { delete[] p; } + ); + for (auto i = 0; i < numParticleOnTile; i++) { + ids.get()[i] = ablastr::particles::localIDtoGlobal(static_cast(aos[i].id()), static_cast(aos[i].cpu())); + } + const auto *const scalar = openPMD::RecordComponent::SCALAR; + currSpecies["id"][scalar].storeChunk(ids, {offset}, {numParticleOnTile64}); + + } + // save "extra" particle properties in AoS and SoA SaveRealProperty(pti, currSpecies, offset, @@ -838,9 +903,10 @@ WarpXOpenPMDPlot::SetupRealProperties (ParticleContainer const * pc, std::set< std::string > addedRecords; // add meta-data per record only once for (auto idx=0; idxNumRealComps(); idx++) { - if (write_real_comp[idx]) { + auto ii = ParticleContainer::NStructReal + idx; // jump over extra AoS names + if (write_real_comp[ii]) { // handle scalar and non-scalar records by name - const auto [record_name, component_name] = detail::name2openPMD(real_comp_names[idx]); + const auto [record_name, component_name] = detail::name2openPMD(real_comp_names[ii]); auto currRecord = currSpecies[record_name]; // meta data for ED-PIC extension @@ -861,9 +927,10 @@ WarpXOpenPMDPlot::SetupRealProperties (ParticleContainer const * pc, } } for (auto idx=0; idx(numParticleOnTile); + auto const numParticleOnTile64 = static_cast( numParticleOnTile ); + auto const& aos = pti.GetArrayOfStructs(); // size = numParticlesOnTile auto const& soa = pti.GetStructOfArrays(); + // first we concatenate the AoS into contiguous arrays + { + // note: WarpX does not yet use extra AoS Real attributes + for( auto idx=0; idx d( + new amrex::ParticleReal[numParticleOnTile], + [](amrex::ParticleReal const *p){ delete[] p; } + ); + + for( auto kk=0; kk AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - void operator() (PData& ptd, int i, + void operator() (const PData& ptd, int i, const amrex::RealVect& /*pos*/, const amrex::RealVect& /*normal*/, amrex::RandomEngine const& /*engine*/) const noexcept { - ptd.id(i) = -ptd.id(i); + auto& p = ptd.m_aos[i]; + p.id() = -p.id(); } }; } diff --git a/Source/EmbeddedBoundary/ParticleScraper.H b/Source/EmbeddedBoundary/ParticleScraper.H index 312b3762a7e..d6196c35f44 100644 --- a/Source/EmbeddedBoundary/ParticleScraper.H +++ b/Source/EmbeddedBoundary/ParticleScraper.H @@ -170,12 +170,13 @@ scrapeParticles (PC& pc, const amrex::Vector& distance_t auto& tile = pti.GetParticleTile(); auto ptd = tile.getParticleTileData(); const auto np = tile.numParticles(); + amrex::Particle<0,0> * const particles = tile.GetArrayOfStructs()().data(); auto phi = (*distance_to_eb[lev])[pti].array(); // signed distance function amrex::ParallelForRNG( np, [=] AMREX_GPU_DEVICE (const int ip, amrex::RandomEngine const& engine) noexcept { // skip particles that are already flagged for removal - if (ptd.id(ip) < 0) return; + if (particles[ip].id() < 0) return; amrex::ParticleReal xp, yp, zp; getPosition(ip, xp, yp, zp); diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H index c69f07acdb2..5c90dab25e6 100644 --- a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H +++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H @@ -72,8 +72,7 @@ class BinaryCollision final // Define shortcuts for frequently-used type names using ParticleType = WarpXParticleContainer::ParticleType; using ParticleTileType = WarpXParticleContainer::ParticleTileType; - using ParticleTileDataType = ParticleTileType::ParticleTileDataType; - using ParticleBins = amrex::DenseBins; + using ParticleBins = amrex::DenseBins; using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; using index_type = ParticleBins::index_type; @@ -262,6 +261,9 @@ public: const amrex::ParticleReal q1 = species_1.getCharge(); const amrex::ParticleReal m1 = species_1.getMass(); auto get_position_1 = GetParticlePosition(ptile_1, getpos_offset); + // Needed to access the particle id + ParticleType * AMREX_RESTRICT + particle_ptr_1 = ptile_1.GetArrayOfStructs()().data(); amrex::Geometry const& geom = WarpX::GetInstance().Geom(lev); #if defined WARPX_DIM_1D_Z @@ -369,7 +371,7 @@ public: soa_1, soa_1, product_species_vector, tile_products_data, - m1, m1, + particle_ptr_1, particle_ptr_1, m1, m1, products_mass, p_mask, products_np, copy_species1, copy_species2, p_pair_indices_1, p_pair_indices_2, @@ -401,6 +403,9 @@ public: const amrex::ParticleReal q1 = species_1.getCharge(); const amrex::ParticleReal m1 = species_1.getMass(); auto get_position_1 = GetParticlePosition(ptile_1, getpos_offset); + // Needed to access the particle id + ParticleType * AMREX_RESTRICT + particle_ptr_1 = ptile_1.GetArrayOfStructs()().data(); // - Species 2 const auto soa_2 = ptile_2.getParticleTileData(); index_type* AMREX_RESTRICT indices_2 = bins_2.permutationPtr(); @@ -408,6 +413,9 @@ public: const amrex::ParticleReal q2 = species_2.getCharge(); const amrex::ParticleReal m2 = species_2.getMass(); auto get_position_2 = GetParticlePosition(ptile_2, getpos_offset); + // Needed to access the particle id + ParticleType * AMREX_RESTRICT + particle_ptr_2 = ptile_2.GetArrayOfStructs()().data(); amrex::Geometry const& geom = WarpX::GetInstance().Geom(lev); #if defined WARPX_DIM_1D_Z @@ -527,7 +535,7 @@ public: soa_1, soa_2, product_species_vector, tile_products_data, - m1, m2, + particle_ptr_1, particle_ptr_2, m1, m2, products_mass, p_mask, products_np, copy_species1, copy_species2, p_pair_indices_1, p_pair_indices_2, diff --git a/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H b/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H index cfdc36d3c50..feb7acf81d3 100644 --- a/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H @@ -23,13 +23,10 @@ * \brief This functor performs pairwise Coulomb collision on a single cell by calling the function * ElasticCollisionPerez. It also reads and contains the Coulomb logarithm. */ -class PairWiseCoulombCollisionFunc -{ +class PairWiseCoulombCollisionFunc{ // Define shortcuts for frequently-used type names using ParticleType = WarpXParticleContainer::ParticleType; - using ParticleTileType = WarpXParticleContainer::ParticleTileType; - using ParticleTileDataType = ParticleTileType::ParticleTileDataType; - using ParticleBins = amrex::DenseBins; + using ParticleBins = amrex::DenseBins; using index_type = ParticleBins::index_type; using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMC.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMC.H index ab01eba2c81..c1be307b811 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMC.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMC.H @@ -38,8 +38,7 @@ class DSMC final // Define shortcuts for frequently-used type names using ParticleType = WarpXParticleContainer::ParticleType; using ParticleTileType = WarpXParticleContainer::ParticleTileType; - using ParticleTileDataType = ParticleTileType::ParticleTileDataType; - using ParticleBins = amrex::DenseBins; + using ParticleBins = amrex::DenseBins; using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; using index_type = ParticleBins::index_type; diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H index 5a7a4f3237a..c1fb7ee7e38 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H @@ -10,9 +10,6 @@ #define SPLIT_AND_SCATTER_FUNC_H_ #include "Particles/Collision/ScatteringProcess.H" -#include "Particles/NamedComponentParticleContainer.H" - -#include /** * \brief Function that performs the particle scattering and injection due @@ -58,6 +55,8 @@ int splitScatteringParticles ( const auto ptile1_data = ptile1.getParticleTileData(); const auto ptile2_data = ptile2.getParticleTileData(); + const Long minus_one_long = -1; + ParallelForRNG(n_total_pairs, [=] AMREX_GPU_DEVICE (int i, RandomEngine const& engine) noexcept { @@ -71,37 +70,20 @@ int splitScatteringParticles ( // starting with the parent particles auto& w1 = ptile1_data.m_rdata[PIdx::w][p_pair_indices_1[i]]; auto& w2 = ptile2_data.m_rdata[PIdx::w][p_pair_indices_2[i]]; - uint64_t* AMREX_RESTRICT idcpu1 = ptile1_data.m_idcpu; - uint64_t* AMREX_RESTRICT idcpu2 = ptile2_data.m_idcpu; - - // Note: Particle::atomicSetID should also be provided as a standalone helper function in AMReX - // to replace the following lambda. - auto const atomicSetIdMinus = [] AMREX_GPU_DEVICE (uint64_t & idcpu) - { - constexpr amrex::Long minus_one_long = -1; - uint64_t tmp = 0; - amrex::ParticleIDWrapper wrapper(tmp); - wrapper = minus_one_long; -#if defined(AMREX_USE_OMP) -#pragma omp atomic write - idcpu = wrapper.m_idata; -#else - auto *old_ptr = reinterpret_cast(&idcpu); - amrex::Gpu::Atomic::Exch(old_ptr, (unsigned long long) wrapper.m_idata); -#endif - }; // Remove p_pair_reaction_weight[i] from the colliding particles' weights. // If the colliding particle weight decreases to zero, remove particle by // setting its id to -1. Gpu::Atomic::AddNoRet(&w1, -p_pair_reaction_weight[i]); if (w1 <= 0._prt) { - atomicSetIdMinus(idcpu1[p_pair_indices_1[i]]); + auto& p = ptile1_data.m_aos[p_pair_indices_1[i]]; + p.atomicSetID(minus_one_long); } Gpu::Atomic::AddNoRet(&w2, -p_pair_reaction_weight[i]); if (w2 <= 0._prt) { - atomicSetIdMinus(idcpu2[p_pair_indices_2[i]]); + auto& p = ptile2_data.m_aos[p_pair_indices_2[i]]; + p.atomicSetID(minus_one_long); } // Set the child particle properties appropriately diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H index b2a2112ca68..397536b67bf 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H @@ -33,13 +33,10 @@ * creation functor. * This functor also reads and contains the fusion multiplier. */ -class NuclearFusionFunc -{ +class NuclearFusionFunc{ // Define shortcuts for frequently-used type names using ParticleType = WarpXParticleContainer::ParticleType; - using ParticleTileType = WarpXParticleContainer::ParticleTileType; - using ParticleTileDataType = ParticleTileType::ParticleTileDataType; - using ParticleBins = amrex::DenseBins; + using ParticleBins = amrex::DenseBins; using index_type = ParticleBins::index_type; using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; @@ -157,13 +154,12 @@ public: // other species and we need to decrease their weight accordingly. // c1 corresponds to the minimum number of times a particle of species 1 will be paired // with a particle of species 2. Same for c2. - // index_type(1): https://github.com/AMReX-Codes/amrex/pull/3684 - const index_type c1 = amrex::max(NI2/NI1, index_type(1)); - const index_type c2 = amrex::max(NI1/NI2, index_type(1)); + const index_type c1 = amrex::max(NI2/NI1,1u); + const index_type c2 = amrex::max(NI1/NI2,1u); // multiplier ratio to take into account unsampled pairs const auto multiplier_ratio = static_cast( - m_isSameSpecies ? 2*max_N - 1 : max_N); + (m_isSameSpecies)?(2u*max_N - 1):(max_N)); #if (defined WARPX_DIM_RZ) amrex::ParticleReal * const AMREX_RESTRICT theta1 = soa_1.m_rdata[PIdx::theta]; diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H index 0b51d6b4b61..7b29267ec32 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H @@ -22,12 +22,10 @@ namespace { // Define shortcuts for frequently-used type names - using SoaData_type = typename WarpXParticleContainer::ParticleTileType::ParticleTileDataType; - using ParticleType = typename WarpXParticleContainer::ParticleType; - using ParticleTileType = typename WarpXParticleContainer::ParticleTileType; - using ParticleTileDataType = typename ParticleTileType::ParticleTileDataType; - using ParticleBins = amrex::DenseBins; - using index_type = typename ParticleBins::index_type; + using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; + using ParticleType = WarpXParticleContainer::ParticleType; + using ParticleBins = amrex::DenseBins; + using index_type = ParticleBins::index_type; /** * \brief This function initializes the momentum of the alpha particles produced from diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionInitializeMomentum.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionInitializeMomentum.H index 52e9db8aa94..be3f5b2d957 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionInitializeMomentum.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionInitializeMomentum.H @@ -24,9 +24,7 @@ namespace { // Define shortcuts for frequently-used type names using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; using ParticleType = WarpXParticleContainer::ParticleType; - using ParticleTileType = WarpXParticleContainer::ParticleTileType; - using ParticleTileDataType = ParticleTileType::ParticleTileDataType; - using ParticleBins = amrex::DenseBins; + using ParticleBins = amrex::DenseBins; using index_type = ParticleBins::index_type; /** diff --git a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H index 3928c74223d..dc830b477df 100644 --- a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H +++ b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H @@ -30,15 +30,13 @@ * \brief This functor creates particles produced from a binary collision and sets their initial * properties (position, momentum, weight). */ -class ParticleCreationFunc -{ +class ParticleCreationFunc{ // Define shortcuts for frequently-used type names - using ParticleType = typename WarpXParticleContainer::ParticleType; - using ParticleTileType = typename WarpXParticleContainer::ParticleTileType; - using ParticleTileDataType = typename ParticleTileType::ParticleTileDataType; - using ParticleBins = amrex::DenseBins; - using index_type = typename ParticleBins::index_type; - using SoaData_type = typename WarpXParticleContainer::ParticleTileType::ParticleTileDataType; + using ParticleType = WarpXParticleContainer::ParticleType; + using ParticleTileType = WarpXParticleContainer::ParticleTileType; + using ParticleBins = amrex::DenseBins; + using index_type = ParticleBins::index_type; + using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; public: /** @@ -71,6 +69,12 @@ public: * @param[in, out] soa_1 struct of array data of the first colliding particle species * @param[in, out] soa_2 struct of array data of the second colliding particle species * @param[out] tile_products array containing tile data of the product particles. + * @param[out] particle_ptr_1 pointer to data of the first colliding particle species. Is + * needed to set the id of a particle to -1 in order to delete it when its weight + * reaches 0. + * @param[out] particle_ptr_2 pointer to data of the second colliding particle species. Is + * needed to set the id of a particle to -1 in order to delete it when its weight + * reaches 0. * @param[in] m1 mass of the first colliding particle species * @param[in] m2 mass of the second colliding particle species * @param[in] products_mass array storing the mass of product particles @@ -98,6 +102,7 @@ public: const SoaData_type& soa_1, const SoaData_type& soa_2, const amrex::Vector& pc_products, ParticleTileType** AMREX_RESTRICT tile_products, + ParticleType* particle_ptr_1, ParticleType* particle_ptr_2, const amrex::ParticleReal& m1, const amrex::ParticleReal& m2, const amrex::Vector& products_mass, const index_type* AMREX_RESTRICT p_mask, @@ -132,8 +137,6 @@ public: amrex::ParticleReal* AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w]; amrex::ParticleReal* AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w]; - uint64_t* AMREX_RESTRICT idcpu1 = soa_1.m_idcpu; - uint64_t* AMREX_RESTRICT idcpu2 = soa_2.m_idcpu; // Create necessary GPU vectors, that will be used in the kernel below amrex::Vector soa_products; @@ -202,32 +205,16 @@ public: amrex::Gpu::Atomic::AddNoRet(&w2[p_pair_indices_2[i]], -p_pair_reaction_weight[i]); - // Note: Particle::atomicSetID should also be provided as a standalone helper function in AMReX - // to replace the following lambda. - auto const atomicSetIdMinus = [] AMREX_GPU_DEVICE (uint64_t & idcpu) - { - constexpr amrex::Long minus_one_long = -1; - uint64_t tmp = 0; - amrex::ParticleIDWrapper wrapper(tmp); - wrapper = minus_one_long; -#if defined(AMREX_USE_OMP) -#pragma omp atomic write - idcpu = wrapper.m_idata; -#else - auto *old_ptr = reinterpret_cast(&idcpu); - amrex::Gpu::Atomic::Exch(old_ptr, (unsigned long long) wrapper.m_idata); -#endif - }; - // If the colliding particle weight decreases to zero, remove particle by // setting its id to -1 + constexpr amrex::Long minus_one_long = -1; if (w1[p_pair_indices_1[i]] <= amrex::ParticleReal(0.)) { - atomicSetIdMinus(idcpu1[p_pair_indices_1[i]]); + particle_ptr_1[p_pair_indices_1[i]].atomicSetID(minus_one_long); } if (w2[p_pair_indices_2[i]] <= amrex::ParticleReal(0.)) { - atomicSetIdMinus(idcpu2[p_pair_indices_2[i]]); + particle_ptr_2[p_pair_indices_2[i]].atomicSetID(minus_one_long); } // Initialize the product particles' momentum, using a function depending on the @@ -307,14 +294,12 @@ private: * \brief This class does nothing and is used as second template parameter for binary collisions * that do not create particles. */ -class NoParticleCreationFunc -{ - using ParticleType = typename WarpXParticleContainer::ParticleType; - using ParticleTileType = typename WarpXParticleContainer::ParticleTileType; - using ParticleTileDataType = typename ParticleTileType::ParticleTileDataType; - using ParticleBins = amrex::DenseBins; - using index_type = typename ParticleBins::index_type; - using SoaData_type = typename WarpXParticleContainer::ParticleTileType::ParticleTileDataType; +class NoParticleCreationFunc{ + using ParticleType = WarpXParticleContainer::ParticleType; + using ParticleTileType = WarpXParticleContainer::ParticleTileType; + using ParticleBins = amrex::DenseBins; + using index_type = ParticleBins::index_type; + using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; public: NoParticleCreationFunc () = default; @@ -328,6 +313,7 @@ public: const SoaData_type& /*soa_1*/, const SoaData_type& /*soa_2*/, amrex::Vector& /*pc_products*/, ParticleTileType** /*tile_products*/, + ParticleType* /*particle_ptr_1*/, ParticleType* /*particle_ptr_2*/, const amrex::ParticleReal& /*m1*/, const amrex::ParticleReal& /*m2*/, const amrex::Vector& /*products_mass*/, const index_type* /*p_mask*/, const amrex::Vector& /*products_np*/, diff --git a/Source/Particles/Collision/BinaryCollision/ShuffleFisherYates.H b/Source/Particles/Collision/BinaryCollision/ShuffleFisherYates.H index 3b8f72f4b84..42259512b0d 100644 --- a/Source/Particles/Collision/BinaryCollision/ShuffleFisherYates.H +++ b/Source/Particles/Collision/BinaryCollision/ShuffleFisherYates.H @@ -12,7 +12,7 @@ /* \brief Shuffle array according to Fisher-Yates algorithm. * Only shuffle the part between is <= i < ie, n = ie-is. * T_index shall be - * amrex::DenseBins::index_type + * amrex::DenseBins::index_type */ template diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index d0822789015..d0db678dfda 100644 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -252,7 +252,7 @@ void doChargeDepositionSharedShapeN (const GetParticlePosition& GetPositio const int n_rz_azimuthal_modes, amrex::Real* cost, const long load_balance_costs_update_algo, - const amrex::DenseBins& a_bins, + const amrex::DenseBins& a_bins, const amrex::Box& box, const amrex::Geometry& geom, const amrex::IntVect& a_tbox_max_size, diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index 7343a8c5626..efe6efcc788 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -592,7 +592,7 @@ void doDepositionSharedShapeN (const GetParticlePosition& GetPosition, int n_rz_azimuthal_modes, amrex::Real* cost, long load_balance_costs_update_algo, - const amrex::DenseBins& a_bins, + const amrex::DenseBins& a_bins, const amrex::Box& box, const amrex::Geometry& geom, const amrex::IntVect& a_tbox_max_size) diff --git a/Source/Particles/ElementaryProcess/QEDPairGeneration.H b/Source/Particles/ElementaryProcess/QEDPairGeneration.H index ca00b56323a..5abc9282d4f 100644 --- a/Source/Particles/ElementaryProcess/QEDPairGeneration.H +++ b/Source/Particles/ElementaryProcess/QEDPairGeneration.H @@ -167,7 +167,7 @@ public: p_ux, p_uy, p_uz, engine); - amrex::ParticleIDWrapper{src.m_idcpu[i_src]} = -1; // destroy photon after pair generation + src.m_aos[i_src].id() = -1; //destroy photon after pair generation } private: diff --git a/Source/Particles/ElementaryProcess/QEDPhotonEmission.H b/Source/Particles/ElementaryProcess/QEDPhotonEmission.H index f3099bf997f..8ba5c63ad57 100644 --- a/Source/Particles/ElementaryProcess/QEDPhotonEmission.H +++ b/Source/Particles/ElementaryProcess/QEDPhotonEmission.H @@ -22,7 +22,6 @@ #include #include #include -#include #include #include @@ -238,11 +237,12 @@ void cleanLowEnergyPhotons( const int old_size, const int num_added, const amrex::ParticleReal energy_threshold) { - auto& soa = ptile.GetStructOfArrays(); - auto p_idcpu = soa.GetIdCPUData().data() + old_size; + auto pp = ptile.GetArrayOfStructs()().data() + old_size; + + const auto& soa = ptile.GetStructOfArrays(); const auto p_ux = soa.GetRealData(PIdx::ux).data() + old_size; - const auto p_uy = soa.GetRealData(PIdx::uy).data() + old_size; - const auto p_uz = soa.GetRealData(PIdx::uz).data() + old_size; + const auto p_uy = soa.GetRealData(PIdx::uy).data() + old_size; + const auto p_uz = soa.GetRealData(PIdx::uz).data() + old_size; //The square of the energy threshold const auto energy_threshold2 = std::max( @@ -251,6 +251,8 @@ void cleanLowEnergyPhotons( amrex::ParallelFor(num_added, [=] AMREX_GPU_DEVICE (int ip) noexcept { + auto& p = pp[ip]; + const auto ux = p_ux[ip]; const auto uy = p_uy[ip]; const auto uz = p_uz[ip]; @@ -260,8 +262,8 @@ void cleanLowEnergyPhotons( constexpr amrex::ParticleReal me_c = PhysConst::m_e*PhysConst::c; const auto phot_energy2 = (ux*ux + uy*uy + uz*uz)*me_c*me_c; - if (phot_energy2 < energy_threshold2) { - amrex::ParticleIDWrapper{p_idcpu[ip]} = -1; + if (phot_energy2 < energy_threshold2){ + p.id() = - 1; } }); } diff --git a/Source/Particles/LaserParticleContainer.H b/Source/Particles/LaserParticleContainer.H index fac94ff20a3..e6fa308431c 100644 --- a/Source/Particles/LaserParticleContainer.H +++ b/Source/Particles/LaserParticleContainer.H @@ -56,9 +56,10 @@ public: * \brief Method to initialize runtime attributes. Does nothing for LaserParticleContainer. */ void DefaultInitializeRuntimeAttributes ( - typename ContainerLike::ParticleTileType& /*pinned_tile*/, - int /*n_external_attr_real*/, - int /*n_external_attr_int*/) final {} + amrex::ParticleTile, + NArrayReal, NArrayInt, amrex::PinnedArenaAllocator>& /*pinned_tile*/, + const int /*n_external_attr_real*/, + const int /*n_external_attr_int*/) final {} void ReadHeader (std::istream& is) final; diff --git a/Source/Particles/NamedComponentParticleContainer.H b/Source/Particles/NamedComponentParticleContainer.H index e04e7dba6df..3be0886425d 100644 --- a/Source/Particles/NamedComponentParticleContainer.H +++ b/Source/Particles/NamedComponentParticleContainer.H @@ -18,19 +18,12 @@ #include -/** Real Particle Attributes stored in amrex::ParticleContainer's struct of array +/** Particle Attributes stored in amrex::ParticleContainer's struct of array */ struct PIdx { enum { -#if !defined (WARPX_DIM_1D_Z) - x, -#endif -#if defined (WARPX_DIM_3D) - y, -#endif - z, - w, ///< weight + w = 0, ///< weight ux, uy, uz, #ifdef WARPX_DIM_RZ theta, ///< RZ needs all three position components @@ -39,15 +32,6 @@ struct PIdx }; }; -/** Integer Particle Attributes stored in amrex::ParticleContainer's struct of array - */ -struct PIdxInt -{ - enum { - nattribs ///< number of attributes - }; -}; - /** Particle Container class that allows to add/access particle components * with a name (string) instead of doing so with an integer index. * (The "components" are all the particle quantities - except those @@ -61,11 +45,11 @@ struct PIdxInt */ template class T_Allocator=amrex::DefaultAllocator> class NamedComponentParticleContainer : -public amrex::ParticleContainerPureSoA +public amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator> { public: /** Construct an empty NamedComponentParticleContainer **/ - NamedComponentParticleContainer () : amrex::ParticleContainerPureSoA() {} + NamedComponentParticleContainer () : amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>() {} /** Construct a NamedComponentParticleContainer from an AmrParGDB object * @@ -77,15 +61,8 @@ public: * AMR hierarchy. Usually, this is generated by an AmrCore or AmrLevel object. */ NamedComponentParticleContainer (amrex::AmrParGDB* amr_pgdb) - : amrex::ParticleContainerPureSoA(amr_pgdb) { + : amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>(amr_pgdb) { // build up the map of string names to particle component numbers -#if !defined (WARPX_DIM_1D_Z) - particle_comps["x"] = PIdx::x; -#endif -#if defined (WARPX_DIM_3D) - particle_comps["y"] = PIdx::y; -#endif - particle_comps["z"] = PIdx::z; particle_comps["w"] = PIdx::w; particle_comps["ux"] = PIdx::ux; particle_comps["uy"] = PIdx::uy; @@ -108,12 +85,12 @@ public: * @param p_ricomps name-to-index map for run-time integer components */ NamedComponentParticleContainer( - amrex::ParticleContainerPureSoA && pc, + amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator> && pc, std::map p_comps, std::map p_icomps, std::map p_rcomps, std::map p_ricomps) - : amrex::ParticleContainerPureSoA(std::move(pc)), + : amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>(std::move(pc)), particle_comps(std::move(p_comps)), particle_icomps(std::move(p_icomps)), particle_runtime_comps(std::move(p_rcomps)), @@ -141,7 +118,7 @@ public: NamedComponentParticleContainer make_alike () const { auto tmp = NamedComponentParticleContainer( - amrex::ParticleContainerPureSoA::template make_alike(), + amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>::template make_alike(), particle_comps, particle_icomps, particle_runtime_comps, @@ -150,10 +127,10 @@ public: return tmp; } - using amrex::ParticleContainerPureSoA::NumRealComps; - using amrex::ParticleContainerPureSoA::NumIntComps; - using amrex::ParticleContainerPureSoA::AddRealComp; - using amrex::ParticleContainerPureSoA::AddIntComp; + using amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>::NumRealComps; + using amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>::NumIntComps; + using amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>::AddRealComp; + using amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>::AddIntComp; /** Allocate a new run-time real component * diff --git a/Source/Particles/ParticleBoundaryBuffer.cpp b/Source/Particles/ParticleBoundaryBuffer.cpp index 88304bd8a9c..54c4396379d 100644 --- a/Source/Particles/ParticleBoundaryBuffer.cpp +++ b/Source/Particles/ParticleBoundaryBuffer.cpp @@ -50,7 +50,7 @@ struct CopyAndTimestamp { void operator() (const DstData& dst, const SrcData& src, int src_i, int dst_i) const noexcept { - dst.m_idcpu[dst_i] = src.m_idcpu[src_i]; + dst.m_aos[dst_i] = src.m_aos[src_i]; for (int j = 0; j < SrcData::NAR; ++j) { dst.m_rdata[j][dst_i] = src.m_rdata[j][src_i]; } @@ -222,7 +222,7 @@ void ParticleBoundaryBuffer::gatherParticles (MultiParticleContainer& mypc, { WARPX_PROFILE("ParticleBoundaryBuffer::gatherParticles"); - using PIter = amrex::ParConstIterSoA; + using PIter = amrex::ParConstIter<0,0,PIdx::nattribs>; const auto& warpx_instance = WarpX::GetInstance(); const amrex::Geometry& geom = warpx_instance.Geom(0); auto plo = geom.ProbLoArray(); diff --git a/Source/Particles/ParticleCreation/SmartCopy.H b/Source/Particles/ParticleCreation/SmartCopy.H index 6a6ceb3d290..2c04baa18bb 100644 --- a/Source/Particles/ParticleCreation/SmartCopy.H +++ b/Source/Particles/ParticleCreation/SmartCopy.H @@ -26,7 +26,7 @@ * type. Second, if a given component name is found in both the src * and the dst, then the src value is copied. * - * Particle positions and id numbers are always copied. + * Particle structs - positions and id numbers - are always copied. * * You don't create this directly - use the SmartCopyFactory object below. */ @@ -48,6 +48,9 @@ struct SmartCopy void operator() (DstData& dst, const SrcData& src, int i_src, int i_dst, amrex::RandomEngine const& engine) const noexcept { + // the particle struct is always copied over + dst.m_aos[i_dst] = src.m_aos[i_src]; + // initialize the real components for (int j = 0; j < DstData::NAR; ++j) { dst.m_rdata[j][i_dst] = initializeRealValue(m_policy_real[j], engine); diff --git a/Source/Particles/ParticleCreation/SmartCreate.H b/Source/Particles/ParticleCreation/SmartCreate.H index 7fb854ee083..67d7767a5d3 100644 --- a/Source/Particles/ParticleCreation/SmartCreate.H +++ b/Source/Particles/ParticleCreation/SmartCreate.H @@ -14,8 +14,6 @@ #include #include #include -#include -#include /** * \brief This is a functor for performing a "smart create" that works @@ -49,23 +47,23 @@ struct SmartCreate const int id = 0) const noexcept { #if defined(WARPX_DIM_3D) - prt.m_rdata[PIdx::x][i_prt] = x; - prt.m_rdata[PIdx::y][i_prt] = y; - prt.m_rdata[PIdx::z][i_prt] = z; + prt.m_aos[i_prt].pos(0) = x; + prt.m_aos[i_prt].pos(1) = y; + prt.m_aos[i_prt].pos(2) = z; #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - prt.m_rdata[PIdx::x][i_prt] = x; - prt.m_rdata[PIdx::z][i_prt] = z; + prt.m_aos[i_prt].pos(0) = x; + prt.m_aos[i_prt].pos(1) = z; amrex::ignore_unused(y); #else - prt.m_rdata[PIdx::z][i_prt] = z; + prt.m_aos[i_prt].pos(0) = z; amrex::ignore_unused(x,y); #endif - amrex::ParticleIDWrapper{prt.m_idcpu[i_prt]} = id; - amrex::ParticleCPUWrapper{prt.m_idcpu[i_prt]} = cpu; + prt.m_aos[i_prt].cpu() = cpu; + prt.m_aos[i_prt].id() = id; - // initialize the real components after position - for (int j = AMREX_SPACEDIM; j < PartData::NAR; ++j) { + // initialize the real components + for (int j = 0; j < PartData::NAR; ++j) { prt.m_rdata[j][i_prt] = initializeRealValue(m_policy_real[j], engine); } for (int j = 0; j < prt.m_num_runtime_real; ++j) { diff --git a/Source/Particles/ParticleCreation/SmartUtils.H b/Source/Particles/ParticleCreation/SmartUtils.H index 6b3d396900d..732a12bb729 100644 --- a/Source/Particles/ParticleCreation/SmartUtils.H +++ b/Source/Particles/ParticleCreation/SmartUtils.H @@ -17,7 +17,6 @@ #include #include #include -#include #include #include @@ -61,12 +60,12 @@ void setNewParticleIDs (PTile& ptile, int old_size, int num_added) } const int cpuid = amrex::ParallelDescriptor::MyProc(); - auto ptd = ptile.getParticleTileData(); + auto pp = ptile.GetArrayOfStructs()().data() + old_size; amrex::ParallelFor(num_added, [=] AMREX_GPU_DEVICE (int ip) noexcept { - auto const new_id = ip + old_size; - amrex::ParticleIDWrapper{ptd.m_idcpu[new_id]} = pid+ip; - amrex::ParticleCPUWrapper{ptd.m_idcpu[new_id]} = cpuid; + auto& p = pp[ip]; + p.id() = pid+ip; + p.cpu() = cpuid; }); } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index edf91a84526..a12ae75f629 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -268,9 +268,11 @@ public: * @param[in] engine the random engine, used in initialization of QED optical depths */ void DefaultInitializeRuntimeAttributes ( - typename ContainerLike::ParticleTileType& pinned_tile, - int n_external_attr_real, - int n_external_attr_int) final; + amrex::ParticleTile, + NArrayReal, NArrayInt, + amrex::PinnedArenaAllocator>& pinned_tile, + int n_external_attr_real, + int n_external_attr_int) final; /** * \brief Apply NCI Godfrey filter to all components of E and B before gather diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index a8ed42ce419..e39cd79b55d 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -198,8 +198,8 @@ namespace * and avoid any possible undefined behavior before the next call to redistribute) and sets * the particle id to -1 so that it can be effectively deleted. * - * \param idcpu particle id soa data - * \param pa particle real soa data + * \param p particle aos data + * \param pa particle soa data * \param ip index for soa data * \param do_field_ionization whether species has ionization * \param pi ionization level data @@ -210,21 +210,20 @@ namespace */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ZeroInitializeAndSetNegativeID ( - uint64_t * AMREX_RESTRICT idcpu, - const GpuArray& pa, long& ip, + ParticleType& p, const GpuArray& pa, long& ip, const bool& do_field_ionization, int* pi #ifdef WARPX_QED - ,const bool& has_quantum_sync, amrex::ParticleReal* AMREX_RESTRICT p_optical_depth_QSR - ,const bool& has_breit_wheeler, amrex::ParticleReal* AMREX_RESTRICT p_optical_depth_BW + ,const bool& has_quantum_sync, amrex::ParticleReal* p_optical_depth_QSR + ,const bool& has_breit_wheeler, amrex::ParticleReal* p_optical_depth_BW #endif ) noexcept { - pa[PIdx::z][ip] = 0._rt; + p.pos(0) = 0._rt; #if (AMREX_SPACEDIM >= 2) - pa[PIdx::x][ip] = 0._rt; + p.pos(1) = 0._rt; #endif #if defined(WARPX_DIM_3D) - pa[PIdx::y][ip] = 0._rt; + p.pos(2) = 0._rt; #endif pa[PIdx::w ][ip] = 0._rt; pa[PIdx::ux][ip] = 0._rt; @@ -239,7 +238,7 @@ namespace if (has_breit_wheeler) {p_optical_depth_BW[ip] = 0._rt;} #endif - amrex::ParticleIDWrapper{idcpu[ip]} = -1; + p.id() = -1; } } @@ -781,9 +780,11 @@ PhysicalParticleContainer::AddPlasmaFromFile(PlasmaInjector & plasma_injector, void PhysicalParticleContainer::DefaultInitializeRuntimeAttributes ( - typename ContainerLike::ParticleTileType& pinned_tile, - int n_external_attr_real, - int n_external_attr_int) + amrex::ParticleTile, + NArrayReal, NArrayInt, + amrex::PinnedArenaAllocator>& pinned_tile, + const int n_external_attr_real, + const int n_external_attr_int) { ParticleCreation::DefaultInitializeRuntimeAttributes(pinned_tile, n_external_attr_real, n_external_attr_int, @@ -1083,7 +1084,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int const int max_new_particles = Scan::ExclusiveSum(counts.size(), counts.data(), offset.data()); // Update NextID to include particles created in this function - int pid; + Long pid; #ifdef AMREX_USE_OMP #pragma omp critical (add_plasma_nextid) #endif @@ -1092,7 +1093,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int ParticleType::NextID(pid+max_new_particles); } WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - static_cast(pid) + static_cast(max_new_particles) < LongParticleIds::LastParticleID, + static_cast(pid + max_new_particles) < LastParticleID, "ERROR: overflow on particle id numbers"); const int cpuid = ParallelDescriptor::MyProc(); @@ -1103,16 +1104,16 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int DefineAndReturnParticleTile(lev, grid_id, tile_id); } - auto old_size = particle_tile.size(); + auto old_size = particle_tile.GetArrayOfStructs().size(); auto new_size = old_size + max_new_particles; particle_tile.resize(new_size); + ParticleType* pp = particle_tile.GetArrayOfStructs()().data() + old_size; auto& soa = particle_tile.GetStructOfArrays(); GpuArray pa; for (int ia = 0; ia < PIdx::nattribs; ++ia) { pa[ia] = soa.GetRealData(ia).data() + old_size; } - uint64_t * AMREX_RESTRICT pa_idcpu = soa.GetIdCPUData().data() + old_size; // user-defined integer and real attributes const auto n_user_int_attribs = static_cast(m_user_int_attribs.size()); const auto n_user_real_attribs = static_cast(m_user_real_attribs.size()); @@ -1225,8 +1226,9 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int for (int i_part = 0; i_part < pcounts[index]; ++i_part) { long ip = poffset[index] + i_part; - amrex::ParticleIDWrapper{pa_idcpu[ip]} = pid+ip; - amrex::ParticleCPUWrapper{pa_idcpu[ip]} = cpuid; + ParticleType& p = pp[ip]; + p.id() = pid+ip; + p.cpu() = cpuid; const XDim3 r = (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) ? // In the refined injection region: use refinement ratio `lrrfac` inj_pos->getPositionUnitBox(i_part, lrrfac, engine) : @@ -1236,7 +1238,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int #if defined(WARPX_DIM_3D) if (!tile_realbox.contains(XDim3{pos.x,pos.y,pos.z})) { - ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi + ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED ,loc_has_quantum_sync, p_optical_depth_QSR ,loc_has_breit_wheeler, p_optical_depth_BW @@ -1247,7 +1249,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) amrex::ignore_unused(k); if (!tile_realbox.contains(XDim3{pos.x,pos.z,0.0_rt})) { - ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi + ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED ,loc_has_quantum_sync, p_optical_depth_QSR ,loc_has_breit_wheeler, p_optical_depth_BW @@ -1258,7 +1260,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int #else amrex::ignore_unused(j,k); if (!tile_realbox.contains(XDim3{pos.z,0.0_rt,0.0_rt})) { - ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi + ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED ,loc_has_quantum_sync, p_optical_depth_QSR ,loc_has_breit_wheeler, p_optical_depth_BW @@ -1297,7 +1299,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int const Real z0 = applyBallisticCorrection(pos, inj_mom, gamma_boost, beta_boost, t); if (!inj_pos->insideBounds(xb, yb, z0)) { - ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi + ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED ,loc_has_quantum_sync, p_optical_depth_QSR ,loc_has_breit_wheeler, p_optical_depth_BW @@ -1311,7 +1313,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int // Remove particle if density below threshold if ( dens < density_min ){ - ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi + ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED ,loc_has_quantum_sync, p_optical_depth_QSR ,loc_has_breit_wheeler, p_optical_depth_BW @@ -1329,7 +1331,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int // If the particle is not within the lab-frame zmin, zmax, etc. // go to the next generated particle. if (!inj_pos->insideBounds(xb, yb, z0_lab)) { - ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi + ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED ,loc_has_quantum_sync, p_optical_depth_QSR ,loc_has_breit_wheeler, p_optical_depth_BW @@ -1341,7 +1343,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int dens = inj_rho->getDensity(pos.x, pos.y, z0_lab); // Remove particle if density below threshold if ( dens < density_min ){ - ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi + ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED ,loc_has_quantum_sync, p_optical_depth_QSR ,loc_has_breit_wheeler, p_optical_depth_BW @@ -1408,17 +1410,17 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int pa[PIdx::uz][ip] = u.z; #if defined(WARPX_DIM_3D) - pa[PIdx::x][ip] = pos.x; - pa[PIdx::y][ip] = pos.y; - pa[PIdx::z][ip] = pos.z; + p.pos(0) = pos.x; + p.pos(1) = pos.y; + p.pos(2) = pos.z; #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) #ifdef WARPX_DIM_RZ pa[PIdx::theta][ip] = theta; #endif - pa[PIdx::x][ip] = xb; - pa[PIdx::z][ip] = pos.z; + p.pos(0) = xb; + p.pos(1) = pos.z; #else - pa[PIdx::z][ip] = pos.z; + p.pos(0) = pos.z; #endif } }); @@ -1643,7 +1645,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, const int max_new_particles = Scan::ExclusiveSum(counts.size(), counts.data(), offset.data()); // Update NextID to include particles created in this function - int pid; + Long pid; #ifdef AMREX_USE_OMP #pragma omp critical (add_plasma_nextid) #endif @@ -1652,23 +1654,23 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, ParticleType::NextID(pid+max_new_particles); } WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - static_cast(pid) + static_cast(max_new_particles) < LongParticleIds::LastParticleID, + static_cast(pid + max_new_particles) < LastParticleID, "overflow on particle id numbers"); const int cpuid = ParallelDescriptor::MyProc(); auto& particle_tile = tmp_pc.DefineAndReturnParticleTile(0, grid_id, tile_id); - auto old_size = particle_tile.size(); + auto old_size = particle_tile.GetArrayOfStructs().size(); auto new_size = old_size + max_new_particles; particle_tile.resize(new_size); + ParticleType* pp = particle_tile.GetArrayOfStructs()().data() + old_size; auto& soa = particle_tile.GetStructOfArrays(); GpuArray pa; for (int ia = 0; ia < PIdx::nattribs; ++ia) { pa[ia] = soa.GetRealData(ia).data() + old_size; } - uint64_t * AMREX_RESTRICT pa_idcpu = soa.GetIdCPUData().data() + old_size; // user-defined integer and real attributes const auto n_user_int_attribs = static_cast(m_user_int_attribs.size()); @@ -1766,8 +1768,9 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, for (int i_part = 0; i_part < pcounts[index]; ++i_part) { const long ip = poffset[index] + i_part; - amrex::ParticleIDWrapper{pa_idcpu[ip]} = pid+ip; - amrex::ParticleCPUWrapper{pa_idcpu[ip]} = cpuid; + ParticleType& p = pp[ip]; + p.id() = pid+ip; + p.cpu() = cpuid; // This assumes the flux_pos is of type InjectorPositionRandomPlane const XDim3 r = (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) ? @@ -1792,19 +1795,19 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, // the particles will be within the domain. #if defined(WARPX_DIM_3D) if (!ParticleUtils::containsInclusive(tile_realbox, XDim3{ppos.x,ppos.y,ppos.z})) { - amrex::ParticleIDWrapper{pa_idcpu[ip]} = -1; + p.id() = -1; continue; } #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) amrex::ignore_unused(k); if (!ParticleUtils::containsInclusive(tile_realbox, XDim3{ppos.x,ppos.z,0.0_prt})) { - amrex::ParticleIDWrapper{pa_idcpu[ip]} = -1; + p.id() = -1; continue; } #else amrex::ignore_unused(j,k); if (!ParticleUtils::containsInclusive(tile_realbox, XDim3{ppos.z,0.0_prt,0.0_prt})) { - amrex::ParticleIDWrapper{pa_idcpu[ip]} = -1; + p.id() = -1; continue; } #endif @@ -1812,7 +1815,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, // If the particle's initial position is not within or on the species's // xmin, xmax, ymin, ymax, zmin, zmax, go to the next generated particle. if (!flux_pos->insideBoundsInclusive(ppos.x, ppos.y, ppos.z)) { - amrex::ParticleIDWrapper{pa_idcpu[ip]} = -1; + p.id() = -1; continue; } @@ -1846,7 +1849,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, Real flux = inj_flux->getFlux(ppos.x, ppos.y, ppos.z, t); // Remove particle if flux is negative or 0 if ( flux <=0 ){ - amrex::ParticleIDWrapper{pa_idcpu[ip]} = -1; + p.id() = -1; continue; } @@ -1905,18 +1908,18 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, UpdatePosition(ppos.x, ppos.y, ppos.z, pu.x, pu.y, pu.z, t_fract); #if defined(WARPX_DIM_3D) - pa[PIdx::x][ip] = ppos.x; - pa[PIdx::y][ip] = ppos.y; - pa[PIdx::z][ip] = ppos.z; + p.pos(0) = ppos.x; + p.pos(1) = ppos.y; + p.pos(2) = ppos.z; #elif defined(WARPX_DIM_RZ) pa[PIdx::theta][ip] = std::atan2(ppos.y, ppos.x); - pa[PIdx::x][ip] = std::sqrt(ppos.x*ppos.x + ppos.y*ppos.y); - pa[PIdx::z][ip] = ppos.z; + p.pos(0) = std::sqrt(ppos.x*ppos.x + ppos.y*ppos.y); + p.pos(1) = ppos.z; #elif defined(WARPX_DIM_XZ) - pa[PIdx::x][ip] = ppos.x; - pa[PIdx::z][ip] = ppos.z; + p.pos(0) = ppos.x; + p.pos(1) = ppos.z; #else - pa[PIdx::z][ip] = ppos.z; + p.pos(0) = ppos.z; #endif } }); @@ -2339,22 +2342,20 @@ PhysicalParticleContainer::SplitParticles (int lev) split_offset[1] /= ppc_nd[1]; split_offset[2] /= ppc_nd[2]; } + // particle Array Of Structs data + auto& particles = pti.GetArrayOfStructs(); // particle Struct Of Arrays data auto& attribs = pti.GetAttribs(); auto& wp = attribs[PIdx::w ]; auto& uxp = attribs[PIdx::ux]; auto& uyp = attribs[PIdx::uy]; auto& uzp = attribs[PIdx::uz]; - - ParticleTileType& ptile = ParticlesAt(lev, pti); - auto& soa = ptile.GetStructOfArrays(); - uint64_t * const AMREX_RESTRICT idcpu = soa.GetIdCPUData().data(); - const long np = pti.numParticles(); for(int i=0; i> attr_int; pctmp_split.AddNParticles(lev, np_split_to_add, - xp, - yp, - zp, - uxp, - uyp, - uzp, - 1, - attr, + xp, yp, zp, uxp, uyp, uzp, + 1, attr, 0, attr_int, - 1, LongParticleIds::NoSplitParticleID); + 1, NoSplitParticleID); // Copy particles from tmp to current particle container constexpr bool local_flag = true; addParticles(pctmp_split,local_flag); diff --git a/Source/Particles/Pusher/GetAndSetPosition.H b/Source/Particles/Pusher/GetAndSetPosition.H index 44641557756..e4477a2a60d 100644 --- a/Source/Particles/Pusher/GetAndSetPosition.H +++ b/Source/Particles/Pusher/GetAndSetPosition.H @@ -30,26 +30,24 @@ void get_particle_position (const WarpXParticleContainer::SuperParticleType& p, amrex::ParticleReal& y, amrex::ParticleReal& z) noexcept { - using namespace amrex::literals; - -#if defined(WARPX_DIM_RZ) - amrex::ParticleReal const theta = p.rdata(T_PIdx::theta); - amrex::ParticleReal const r = p.pos(T_PIdx::x); +#ifdef WARPX_DIM_RZ + const amrex::ParticleReal theta = p.rdata(T_PIdx::theta); + const amrex::ParticleReal r = p.pos(0); x = r*std::cos(theta); y = r*std::sin(theta); - z = p.pos(PIdx::z); -#elif defined(WARPX_DIM_3D) - x = p.pos(PIdx::x); - y = p.pos(PIdx::y); - z = p.pos(PIdx::z); -#elif defined(WARPX_DIM_XZ) - x = p.pos(PIdx::x); - y = 0_prt; - z = p.pos(PIdx::z); + z = p.pos(1); +#elif WARPX_DIM_3D + x = p.pos(0); + y = p.pos(1); + z = p.pos(2); +#elif WARPX_DIM_XZ + x = p.pos(0); + y = amrex::ParticleReal(0.0); + z = p.pos(1); #else - x = 0_prt; - y = 0_prt; - z = p.pos(PIdx::z); + x = amrex::ParticleReal(0.0); + y = amrex::ParticleReal(0.0); + z = p.pos(0); #endif } @@ -61,19 +59,10 @@ void get_particle_position (const WarpXParticleContainer::SuperParticleType& p, template struct GetParticlePosition { + using PType = WarpXParticleContainer::ParticleType; using RType = amrex::ParticleReal; -#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) - const RType* AMREX_RESTRICT m_x = nullptr; - const RType* AMREX_RESTRICT m_z = nullptr; -#elif defined(WARPX_DIM_3D) - const RType* AMREX_RESTRICT m_x = nullptr; - const RType* AMREX_RESTRICT m_y = nullptr; - const RType* AMREX_RESTRICT m_z = nullptr; -#elif defined(WARPX_DIM_1D_Z) - const RType* AMREX_RESTRICT m_z = nullptr; -#endif - + const PType* AMREX_RESTRICT m_structs = nullptr; #if defined(WARPX_DIM_RZ) const RType* m_theta = nullptr; #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) @@ -95,19 +84,10 @@ struct GetParticlePosition template GetParticlePosition (const ptiType& a_pti, long a_offset = 0) noexcept { - const auto& soa = a_pti.GetStructOfArrays(); - -#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) - m_x = soa.GetRealData(PIdx::x).dataPtr() + a_offset; - m_z = soa.GetRealData(PIdx::z).dataPtr() + a_offset; -#elif defined(WARPX_DIM_3D) - m_x = soa.GetRealData(PIdx::x).dataPtr() + a_offset; - m_y = soa.GetRealData(PIdx::y).dataPtr() + a_offset; - m_z = soa.GetRealData(PIdx::z).dataPtr() + a_offset; -#elif defined(WARPX_DIM_1D_Z) - m_z = soa.GetRealData(PIdx::z).dataPtr() + a_offset; -#endif + const auto& aos = a_pti.GetArrayOfStructs(); + m_structs = aos().dataPtr() + a_offset; #if defined(WARPX_DIM_RZ) + const auto& soa = a_pti.GetStructOfArrays(); m_theta = soa.GetRealData(T_PIdx::theta).dataPtr() + a_offset; #endif } @@ -118,23 +98,24 @@ struct GetParticlePosition AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() (const long i, RType& x, RType& y, RType& z) const noexcept { + const PType& p = m_structs[i]; #ifdef WARPX_DIM_RZ - RType const r = m_x[i]; + const RType r = p.pos(0); x = r*std::cos(m_theta[i]); y = r*std::sin(m_theta[i]); - z = m_z[i]; + z = p.pos(1); #elif WARPX_DIM_3D - x = m_x[i]; - y = m_y[i]; - z = m_z[i]; + x = p.pos(0); + y = p.pos(1); + z = p.pos(2); #elif WARPX_DIM_XZ - x = m_x[i]; + x = p.pos(0); y = m_y_default; - z = m_z[i]; + z = p.pos(1); #else x = m_x_default; y = m_y_default; - z = m_z[i]; + z = p.pos(0); #endif } @@ -146,22 +127,23 @@ struct GetParticlePosition AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void AsStored (const long i, RType& x, RType& y, RType& z) const noexcept { + const PType& p = m_structs[i]; #ifdef WARPX_DIM_RZ - x = m_x[i]; + x = p.pos(0); y = m_theta[i]; - z = m_z[i]; + z = p.pos(1); #elif WARPX_DIM_3D - x = m_x[i]; - y = m_y[i]; - z = m_z[i]; + x = p.pos(0); + y = p.pos(1); + z = p.pos(2); #elif WARPX_DIM_XZ - x = m_x[i]; + x = p.pos(0); y = m_y_default; - z = m_z[i]; + z = p.pos(1); #else x = m_x_default; y = m_y_default; - z = m_z[i]; + z = p.pos(0); #endif } }; @@ -176,18 +158,10 @@ struct GetParticlePosition template struct SetParticlePosition { + using PType = WarpXParticleContainer::ParticleType; using RType = amrex::ParticleReal; -#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) - RType* AMREX_RESTRICT m_x; - RType* AMREX_RESTRICT m_z; -#elif defined(WARPX_DIM_3D) - RType* AMREX_RESTRICT m_x; - RType* AMREX_RESTRICT m_y; - RType* AMREX_RESTRICT m_z; -#elif defined(WARPX_DIM_1D_Z) - RType* AMREX_RESTRICT m_z; -#endif + PType* AMREX_RESTRICT m_structs; #if defined(WARPX_DIM_RZ) RType* AMREX_RESTRICT m_theta; #endif @@ -195,18 +169,10 @@ struct SetParticlePosition template SetParticlePosition (const ptiType& a_pti, long a_offset = 0) noexcept { - auto& soa = a_pti.GetStructOfArrays(); -#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) - m_x = soa.GetRealData(PIdx::x).dataPtr() + a_offset; - m_z = soa.GetRealData(PIdx::z).dataPtr() + a_offset; -#elif defined(WARPX_DIM_3D) - m_x = soa.GetRealData(PIdx::x).dataPtr() + a_offset; - m_y = soa.GetRealData(PIdx::y).dataPtr() + a_offset; - m_z = soa.GetRealData(PIdx::z).dataPtr() + a_offset; -#elif defined(WARPX_DIM_1D_Z) - m_z = soa.GetRealData(PIdx::z).dataPtr() + a_offset; -#endif + auto& aos = a_pti.GetArrayOfStructs(); + m_structs = aos().dataPtr() + a_offset; #if defined(WARPX_DIM_RZ) + auto& soa = a_pti.GetStructOfArrays(); m_theta = soa.GetRealData(T_PIdx::theta).dataPtr() + a_offset; #endif } @@ -224,17 +190,17 @@ struct SetParticlePosition #endif #ifdef WARPX_DIM_RZ m_theta[i] = std::atan2(y, x); - m_x[i] = std::sqrt(x*x + y*y); - m_z[i] = z; + m_structs[i].pos(0) = std::sqrt(x*x + y*y); + m_structs[i].pos(1) = z; #elif WARPX_DIM_3D - m_x[i] = x; - m_y[i] = y; - m_z[i] = z; + m_structs[i].pos(0) = x; + m_structs[i].pos(1) = y; + m_structs[i].pos(2) = z; #elif WARPX_DIM_XZ - m_x[i] = x; - m_z[i] = z; + m_structs[i].pos(0) = x; + m_structs[i].pos(1) = z; #else - m_z[i] = z; + m_structs[i].pos(0) = z; #endif } @@ -252,18 +218,18 @@ struct SetParticlePosition amrex::ignore_unused(x,y); #endif #ifdef WARPX_DIM_RZ - m_x[i] = x; + m_structs[i].pos(0) = x; m_theta[i] = y; - m_z[i] = z; + m_structs[i].pos(1) = z; #elif WARPX_DIM_3D - m_x[i] = x; - m_y[i] = y; - m_z[i] = z; + m_structs[i].pos(0) = x; + m_structs[i].pos(1) = y; + m_structs[i].pos(2) = z; #elif WARPX_DIM_XZ - m_x[i] = x; - m_z[i] = z; + m_structs[i].pos(0) = x; + m_structs[i].pos(1) = z; #else - m_z[i] = z; + m_structs[i].pos(0) = z; #endif } }; diff --git a/Source/Particles/Resampling/LevelingThinning.cpp b/Source/Particles/Resampling/LevelingThinning.cpp index 40f75c76eab..680e33ebe6a 100644 --- a/Source/Particles/Resampling/LevelingThinning.cpp +++ b/Source/Particles/Resampling/LevelingThinning.cpp @@ -60,7 +60,8 @@ void LevelingThinning::operator() (WarpXParIter& pti, const int lev, auto& ptile = pc->ParticlesAt(lev, pti); auto& soa = ptile.GetStructOfArrays(); amrex::ParticleReal * const AMREX_RESTRICT w = soa.GetRealData(PIdx::w).data(); - auto * const AMREX_RESTRICT idcpu = soa.GetIdCPUData().data(); + WarpXParticleContainer::ParticleType * const AMREX_RESTRICT + particle_ptr = ptile.GetArrayOfStructs()().data(); // Using this function means that we must loop over the cells in the ParallelFor. In the case // of the leveling thinning algorithm, it would have possibly been more natural and more @@ -113,7 +114,7 @@ void LevelingThinning::operator() (WarpXParIter& pti, const int lev, // Remove particle with probability 1 - particle_weight/level_weight if (random_number > w[indices[i]]/level_weight) { - amrex::ParticleIDWrapper{idcpu[indices[i]]} = -1; + particle_ptr[indices[i]].id() = -1; } // Set particle weight to level weight otherwise else diff --git a/Source/Particles/Sorting/Partition.cpp b/Source/Particles/Sorting/Partition.cpp index 58e3450f47d..58511cfd5e7 100644 --- a/Source/Particles/Sorting/Partition.cpp +++ b/Source/Particles/Sorting/Partition.cpp @@ -61,7 +61,7 @@ PhysicalParticleContainer::PartitionParticlesInBuffers( // Initialize temporary arrays Gpu::DeviceVector inexflag; inexflag.resize(np); - Gpu::DeviceVector pid; + Gpu::DeviceVector pid; pid.resize(np); // First, partition particles into the larger buffer @@ -109,7 +109,7 @@ PhysicalParticleContainer::PartitionParticlesInBuffers( // - For each particle in the large buffer, find whether it is in // the smaller buffer, by looking up the mask. Store the answer in `inexflag`. amrex::ParallelFor( np - n_fine, - fillBufferFlagRemainingParticles(pti, bmasks, inexflag, Geom(lev), pid, int(n_fine)) ); + fillBufferFlagRemainingParticles(pti, bmasks, inexflag, Geom(lev), pid, n_fine) ); auto *const sep2 = stablePartition( sep, pid.end(), inexflag ); if (bmasks == gather_masks) { diff --git a/Source/Particles/Sorting/SortingUtils.H b/Source/Particles/Sorting/SortingUtils.H index ba7761bf48a..ac2c63e88f8 100644 --- a/Source/Particles/Sorting/SortingUtils.H +++ b/Source/Particles/Sorting/SortingUtils.H @@ -12,7 +12,6 @@ #include #include -#include /** \brief Fill the elements of the input vector with consecutive integer, @@ -20,7 +19,7 @@ * * \param[inout] v Vector of integers, to be filled by this routine */ -void fillWithConsecutiveIntegers( amrex::Gpu::DeviceVector& v ); +void fillWithConsecutiveIntegers( amrex::Gpu::DeviceVector& v ); /** \brief Find the indices that would reorder the elements of `predicate` * so that the elements with non-zero value precede the other elements @@ -42,7 +41,7 @@ ForwardIterator stablePartition(ForwardIterator const index_begin, int const* AMREX_RESTRICT predicate_ptr = predicate.dataPtr(); int N = static_cast(std::distance(index_begin, index_end)); auto num_true = amrex::StablePartition(&(*index_begin), N, - [predicate_ptr] AMREX_GPU_DEVICE (int i) { return predicate_ptr[i]; }); + [predicate_ptr] AMREX_GPU_DEVICE (long i) { return predicate_ptr[i]; }); ForwardIterator sep = index_begin; std::advance(sep, num_true); @@ -50,7 +49,7 @@ ForwardIterator stablePartition(ForwardIterator const index_begin, // On CPU: Use std library ForwardIterator const sep = std::stable_partition( index_begin, index_end, - [&predicate](int i) { return predicate[i]; } + [&predicate](long i) { return predicate[i]; } ); #endif return sep; @@ -89,7 +88,7 @@ class fillBufferFlag // Extract simple structure that can be used directly on the GPU m_domain{geom.Domain()}, m_inexflag_ptr{inexflag.dataPtr()}, - m_ptd{pti.GetParticleTile().getConstParticleTileData()}, + m_particles{pti.GetArrayOfStructs().data()}, m_buffer_mask{(*bmasks)[pti].array()} { for (int idim=0; idim m_buffer_mask; amrex::GpuArray m_prob_lo; amrex::GpuArray m_inv_cell_size; @@ -140,12 +141,12 @@ class fillBufferFlagRemainingParticles amrex::iMultiFab const* bmasks, amrex::Gpu::DeviceVector& inexflag, amrex::Geometry const& geom, - amrex::Gpu::DeviceVector const& particle_indices, - int start_index ) : + amrex::Gpu::DeviceVector const& particle_indices, + long const start_index ) : m_domain{geom.Domain()}, // Extract simple structure that can be used directly on the GPU m_inexflag_ptr{inexflag.dataPtr()}, - m_ptd{pti.GetParticleTile().getConstParticleTileData()}, + m_particles{pti.GetArrayOfStructs().data()}, m_buffer_mask{(*bmasks)[pti].array()}, m_start_index{start_index}, m_indices_ptr{particle_indices.dataPtr()} @@ -158,11 +159,11 @@ class fillBufferFlagRemainingParticles AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - void operator()( const int i ) const { + void operator()( const long i ) const { // Select a particle - auto const j = m_indices_ptr[i+m_start_index]; + auto const& p = m_particles[m_indices_ptr[i+m_start_index]]; // Find the index of the cell where this particle is located - amrex::IntVect const iv = amrex::getParticleCell( m_ptd, j, + amrex::IntVect const iv = amrex::getParticleCell( p, m_prob_lo, m_inv_cell_size, m_domain ); // Find the value of the buffer flag in this cell and // store it at the corresponding particle position in the array `inexflag` @@ -174,10 +175,10 @@ class fillBufferFlagRemainingParticles amrex::GpuArray m_inv_cell_size; amrex::Box m_domain; int* m_inexflag_ptr; - WarpXParticleContainer::ParticleTileType::ConstParticleTileDataType const m_ptd; + WarpXParticleContainer::ParticleType const* m_particles; amrex::Array4 m_buffer_mask; - int const m_start_index; - int const* m_indices_ptr; + long const m_start_index; + long const* m_indices_ptr; }; /** \brief Functor that copies the elements of `src` into `dst`, @@ -194,7 +195,7 @@ class copyAndReorder copyAndReorder( amrex::Gpu::DeviceVector const& src, amrex::Gpu::DeviceVector& dst, - amrex::Gpu::DeviceVector const& indices ): + amrex::Gpu::DeviceVector const& indices ): // Extract simple structure that can be used directly on the GPU m_src_ptr{src.dataPtr()}, m_dst_ptr{dst.dataPtr()}, @@ -202,14 +203,14 @@ class copyAndReorder {} AMREX_GPU_DEVICE AMREX_FORCE_INLINE - void operator()( const int ip ) const { + void operator()( const long ip ) const { m_dst_ptr[ip] = m_src_ptr[ m_indices_ptr[ip] ]; } private: T const* m_src_ptr; T* m_dst_ptr; - int const* m_indices_ptr; + long const* m_indices_ptr; }; #endif // WARPX_PARTICLES_SORTING_SORTINGUTILS_H_ diff --git a/Source/Particles/Sorting/SortingUtils.cpp b/Source/Particles/Sorting/SortingUtils.cpp index cd4b6a13c76..699119e8e18 100644 --- a/Source/Particles/Sorting/SortingUtils.cpp +++ b/Source/Particles/Sorting/SortingUtils.cpp @@ -8,7 +8,7 @@ #include "SortingUtils.H" -void fillWithConsecutiveIntegers( amrex::Gpu::DeviceVector& v ) +void fillWithConsecutiveIntegers( amrex::Gpu::DeviceVector& v ) { #ifdef AMREX_USE_GPU // On GPU: Use amrex diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index a244afc2389..33aa71d1c7d 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -49,10 +49,10 @@ class WarpXParIter - : public amrex::ParIterSoA + : public amrex::ParIter<0,0,PIdx::nattribs> { public: - using amrex::ParIterSoA::ParIterSoA; + using amrex::ParIter<0,0,PIdx::nattribs>::ParIter; WarpXParIter (ContainerType& pc, int level); @@ -89,7 +89,7 @@ public: * particle container classes (that store a collection of particles) derive. Derived * classes can be used for plasma particles, photon particles, or non-physical * particles (e.g., for the laser antenna). - * It derives from amrex::ParticleContainerPureSoA, where the + * It derives from amrex::ParticleContainer<0,0,PIdx::nattribs>, where the * template arguments stand for the number of int and amrex::Real SoA and AoS * data in amrex::Particle. * - AoS amrex::Real: x, y, z (default), 0 additional (first template @@ -164,9 +164,11 @@ public: * class. */ virtual void DefaultInitializeRuntimeAttributes ( - typename ContainerLike::ParticleTileType& pinned_tile, - int n_external_attr_real, - int n_external_attr_int) = 0; + amrex::ParticleTile, + NArrayReal, NArrayInt, + amrex::PinnedArenaAllocator>& pinned_tile, + int n_external_attr_real, + int n_external_attr_int) = 0; /// /// This pushes the particle positions by one half time step. diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 0d927caacc8..85bb1c3f4b8 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -75,13 +75,13 @@ using namespace amrex; WarpXParIter::WarpXParIter (ContainerType& pc, int level) - : amrex::ParIterSoA(pc, level, + : amrex::ParIter<0,0,PIdx::nattribs>(pc, level, MFItInfo().SetDynamic(WarpX::do_dynamic_scheduling)) { } WarpXParIter::WarpXParIter (ContainerType& pc, int level, MFItInfo& info) - : amrex::ParIterSoA(pc, level, + : amrex::ParIter<0,0,PIdx::nattribs>(pc, level, info.SetDynamic(WarpX::do_dynamic_scheduling)) { } @@ -198,54 +198,52 @@ WarpXParticleContainer::AddNParticles (int /*lev*/, long n, // Redistribute() will move them to proper places. auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - using PinnedTile = typename ContainerLike::ParticleTileType; + using PinnedTile = amrex::ParticleTile, + NArrayReal, NArrayInt, + amrex::PinnedArenaAllocator>; PinnedTile pinned_tile; pinned_tile.define(NumRuntimeRealComps(), NumRuntimeIntComps()); const std::size_t np = iend-ibegin; #ifdef WARPX_DIM_RZ - amrex::Vector r(np); amrex::Vector theta(np); #endif for (auto i = ibegin; i < iend; ++i) { - auto & idcpu_data = pinned_tile.GetStructOfArrays().GetIdCPUData(); - idcpu_data.push_back(0); - if (id==-1) { - amrex::ParticleIDWrapper{idcpu_data.back()} = ParticleType::NextID(); + ParticleType p; + if (id==-1) + { + p.id() = ParticleType::NextID(); } else { - amrex::ParticleIDWrapper{idcpu_data.back()} = id; + p.id() = id; } - amrex::ParticleCPUWrapper(idcpu_data.back()) = ParallelDescriptor::MyProc(); - -#ifdef WARPX_DIM_RZ - r[i-ibegin] = std::sqrt(x[i]*x[i] + y[i]*y[i]); - theta[i-ibegin] = std::atan2(y[i], x[i]); -#endif - } - - if (np > 0) - { + p.cpu() = amrex::ParallelDescriptor::MyProc(); #if defined(WARPX_DIM_3D) - pinned_tile.push_back_real(PIdx::x, x.data() + ibegin, x.data() + iend); - pinned_tile.push_back_real(PIdx::y, y.data() + ibegin, y.data() + iend); - pinned_tile.push_back_real(PIdx::z, z.data() + ibegin, z.data() + iend); + p.pos(0) = x[i]; + p.pos(1) = y[i]; + p.pos(2) = z[i]; #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) amrex::ignore_unused(y); #ifdef WARPX_DIM_RZ - pinned_tile.push_back_real(PIdx::x, r.data(), r.data() + np); + theta[i-ibegin] = std::atan2(y[i], x[i]); + p.pos(0) = std::sqrt(x[i]*x[i] + y[i]*y[i]); #else - pinned_tile.push_back_real(PIdx::x, x.data() + ibegin, x.data() + iend); + p.pos(0) = x[i]; #endif - pinned_tile.push_back_real(PIdx::z, z.data() + ibegin, z.data() + iend); + p.pos(1) = z[i]; #else //AMREX_SPACEDIM == 1 amrex::ignore_unused(x,y); - pinned_tile.push_back_real(PIdx::z, z.data() + ibegin, z.data() + iend); + p.pos(0) = z[i]; #endif - pinned_tile.push_back_real(PIdx::w, attr_real[0].data() + ibegin, attr_real[0].data() + iend); + pinned_tile.push_back(p); + } + + if (np > 0) + { + pinned_tile.push_back_real(PIdx::w , attr_real[0].data() + ibegin, attr_real[0].data() + iend); pinned_tile.push_back_real(PIdx::ux, ux.data() + ibegin, ux.data() + iend); pinned_tile.push_back_real(PIdx::uy, uy.data() + ibegin, uy.data() + iend); pinned_tile.push_back_real(PIdx::uz, uz.data() + ibegin, uz.data() + iend); @@ -477,14 +475,15 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, //sort particles by bin WARPX_PROFILE_VAR_START(blp_sort); - amrex::DenseBins bins; + amrex::DenseBins bins; { auto& ptile = ParticlesAt(lev, pti); - auto ptd = ptile.getParticleTileData(); + auto& aos = ptile.GetArrayOfStructs(); + auto *pstruct_ptr = aos().dataPtr(); const int ntiles = numTilesInBox(box, true, bin_size); - bins.build(ptile.numParticles(), ptd, ntiles, + bins.build(ptile.numParticles(), pstruct_ptr, ntiles, [=] AMREX_GPU_HOST_DEVICE (const ParticleType& p) -> unsigned int { Box tbox; @@ -884,7 +883,7 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector const& wp, // HACK - sort particles by bin here. WARPX_PROFILE_VAR_START(blp_sort); - amrex::DenseBins bins; + amrex::DenseBins bins; { const Geometry& geom = Geom(lev); const auto dxi = geom.InvCellSizeArray(); @@ -892,15 +891,16 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector const& wp, const auto domain = geom.Domain(); auto& ptile = ParticlesAt(lev, pti); - auto ptd = ptile.getParticleTileData(); + auto& aos = ptile.GetArrayOfStructs(); + auto *pstruct_ptr = aos().dataPtr(); Box box = pti.validbox(); box.grow(ng_rho); const amrex::IntVect bin_size = WarpX::shared_tilesize; const int ntiles = numTilesInBox(box, true, bin_size); - bins.build(ptile.numParticles(), ptd, ntiles, - [=] AMREX_GPU_HOST_DEVICE (ParticleType const & p) -> unsigned int + bins.build(ptile.numParticles(), pstruct_ptr, ntiles, + [=] AMREX_GPU_HOST_DEVICE (const ParticleType& p) -> unsigned int { Box tbx; auto iv = getParticleCell(p, plo, dxi, domain); @@ -920,7 +920,8 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector const& wp, const auto domain = geom.Domain(); auto& ptile = ParticlesAt(lev, pti); - auto ptd = ptile.getParticleTileData(); + auto& aos = ptile.GetArrayOfStructs(); + auto *pstruct_ptr = aos().dataPtr(); Box box = pti.validbox(); box.grow(ng_rho); @@ -934,10 +935,9 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector const& wp, const auto bin_start = offsets_ptr[ibin]; const auto bin_stop = offsets_ptr[ibin+1]; if (bin_start < bin_stop) { - // static_cast until https://github.com/AMReX-Codes/amrex/pull/3684 - auto const i = static_cast(permutation[bin_start]); + auto p = pstruct_ptr[permutation[bin_start]]; Box tbx; - auto iv = getParticleCell(ptd, i, plo, dxi, domain); + auto iv = getParticleCell(p, plo, dxi, domain); AMREX_ASSERT(box.contains(iv)); [[maybe_unused]] auto tid = getTileIndex(iv, box, true, bin_size, tbx); AMREX_ASSERT(tid == ibin); @@ -1426,10 +1426,10 @@ WarpXParticleContainer::particlePostLocate(ParticleType& p, // Tag particle if goes to higher level. // It will be split later in the loop if (pld.m_lev == lev+1 - and p.id() != amrex::LongParticleIds::NoSplitParticleID + and p.id() != NoSplitParticleID and p.id() >= 0) { - p.id() = amrex::LongParticleIds::DoSplitParticleID; + p.id() = DoSplitParticleID; } if (pld.m_lev == lev-1){ @@ -1468,9 +1468,9 @@ WarpXParticleContainer::ApplyBoundaryConditions (){ const Real zmax = Geom(lev).ProbHi(WARPX_ZINDEX); ParticleTileType& ptile = ParticlesAt(lev, pti); + ParticleType * const pp = ptile.GetArrayOfStructs()().data(); auto& soa = ptile.GetStructOfArrays(); - uint64_t * const AMREX_RESTRICT idcpu = soa.GetIdCPUData().data(); amrex::ParticleReal * const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data(); amrex::ParticleReal * const AMREX_RESTRICT uy = soa.GetRealData(PIdx::uy).data(); amrex::ParticleReal * const AMREX_RESTRICT uz = soa.GetRealData(PIdx::uz).data(); @@ -1479,9 +1479,10 @@ WarpXParticleContainer::ApplyBoundaryConditions (){ amrex::ParallelForRNG( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i, amrex::RandomEngine const& engine) { + ParticleType& p = pp[i]; + // skip particles that are already flagged for removal - auto const id = amrex::ParticleIDWrapper{idcpu[i]}; - if (id < 0) { return; } + if (p.id() < 0) { return; } ParticleReal x, y, z; GetPosition.AsStored(i, x, y, z); @@ -1503,7 +1504,7 @@ WarpXParticleContainer::ApplyBoundaryConditions (){ boundary_conditions, engine); if (particle_lost) { - amrex::ParticleIDWrapper{idcpu[i]} = -id; + p.id() = -p.id(); } else { SetPosition.AsStored(i, x, y, z); } diff --git a/Source/Python/Particles/ParticleBoundaryBuffer.cpp b/Source/Python/Particles/ParticleBoundaryBuffer.cpp index b04ac75e600..2a35faece9b 100644 --- a/Source/Python/Particles/ParticleBoundaryBuffer.cpp +++ b/Source/Python/Particles/ParticleBoundaryBuffer.cpp @@ -10,13 +10,13 @@ namespace warpx { class BoundaryBufferParIter - : public amrex::ParIterSoA + : public amrex::ParIter<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator> { public: - using amrex::ParIterSoA::ParIterSoA; + using amrex::ParIter<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator>::ParIter; BoundaryBufferParIter(ContainerType& pc, int level) : - amrex::ParIterSoA(pc, level) {} + amrex::ParIter<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator>(pc, level) {} }; } @@ -24,9 +24,9 @@ void init_BoundaryBufferParIter (py::module& m) { py::class_< warpx::BoundaryBufferParIter, - amrex::ParIterSoA + amrex::ParIter<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator> >(m, "BoundaryBufferParIter") - .def(py::init::ContainerType&, int>(), + .def(py::init::ContainerType&, int>(), py::arg("particle_container"), py::arg("level") ) ; diff --git a/Source/Python/Particles/PinnedMemoryParticleContainer.cpp b/Source/Python/Particles/PinnedMemoryParticleContainer.cpp index d4f6a422dbe..600d56a62c9 100644 --- a/Source/Python/Particles/PinnedMemoryParticleContainer.cpp +++ b/Source/Python/Particles/PinnedMemoryParticleContainer.cpp @@ -13,6 +13,6 @@ void init_PinnedMemoryParticleContainer (py::module& m) { py::class_< PinnedMemoryParticleContainer, - amrex::ParticleContainerPureSoA + amrex::ParticleContainer<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator> > pmpc (m, "PinnedMemoryParticleContainer"); } diff --git a/Source/Python/Particles/WarpXParticleContainer.cpp b/Source/Python/Particles/WarpXParticleContainer.cpp index 07793a373f3..1473a750941 100644 --- a/Source/Python/Particles/WarpXParticleContainer.cpp +++ b/Source/Python/Particles/WarpXParticleContainer.cpp @@ -12,11 +12,11 @@ void init_WarpXParIter (py::module& m) { py::class_< - WarpXParIter, amrex::ParIterSoA + WarpXParIter, amrex::ParIter<0,0,PIdx::nattribs> >(m, "WarpXParIter") - .def(py::init::ContainerType&, int>(), + .def(py::init::ContainerType&, int>(), py::arg("particle_container"), py::arg("level")) - .def(py::init::ContainerType&, int, amrex::MFItInfo&>(), + .def(py::init::ContainerType&, int, amrex::MFItInfo&>(), py::arg("particle_container"), py::arg("level"), py::arg("info")) ; @@ -26,11 +26,11 @@ void init_WarpXParticleContainer (py::module& m) { py::class_< WarpXParticleContainer, - amrex::ParticleContainerPureSoA + amrex::ParticleContainer<0, 0, PIdx::nattribs, 0> > wpc (m, "WarpXParticleContainer"); wpc .def("add_real_comp", - [](WarpXParticleContainer& pc, const std::string& name, bool comm) { pc.AddRealComp(name, comm); }, + [](WarpXParticleContainer& pc, const std::string& name, bool const comm) { pc.AddRealComp(name, comm); }, py::arg("name"), py::arg("comm") ) .def("add_n_particles", @@ -93,14 +93,6 @@ void init_WarpXParticleContainer (py::module& m) }, py::arg("comp_name") ) - .def("get_icomp_index", - [](WarpXParticleContainer& pc, std::string comp_name) - { - auto particle_comps = pc.getParticleiComps(); - return particle_comps.at(comp_name); - }, - py::arg("comp_name") - ) .def("num_local_tiles_at_level", &WarpXParticleContainer::numLocalTilesAtLevel, py::arg("level") diff --git a/Source/Utils/ParticleUtils.H b/Source/Utils/ParticleUtils.H index 7e3c89228ea..b04176d4d83 100644 --- a/Source/Utils/ParticleUtils.H +++ b/Source/Utils/ParticleUtils.H @@ -28,10 +28,9 @@ namespace ParticleUtils { * @param[in] mfi the MultiFAB iterator. * @param[in] ptile the particle tile. */ - amrex::DenseBins - findParticlesInEachCell (int lev, - amrex::MFIter const & mfi, - WarpXParticleContainer::ParticleTileType & ptile); + amrex::DenseBins + findParticlesInEachCell(int lev, amrex::MFIter const& mfi, + WarpXParticleContainer::ParticleTileType const& ptile); /** * \brief Return (relativistic) particle energy given velocity and mass. diff --git a/Source/Utils/ParticleUtils.cpp b/Source/Utils/ParticleUtils.cpp index b8207b61fa0..60e04f12b86 100644 --- a/Source/Utils/ParticleUtils.cpp +++ b/Source/Utils/ParticleUtils.cpp @@ -22,28 +22,24 @@ #include #include -namespace ParticleUtils -{ +namespace ParticleUtils { using namespace amrex; - // Define shortcuts for frequently-used type names - using ParticleType = typename WarpXParticleContainer::ParticleType; - using ParticleTileType = typename WarpXParticleContainer::ParticleTileType; - using ParticleTileDataType = typename ParticleTileType::ParticleTileDataType; - using ParticleBins = DenseBins; - using index_type = typename ParticleBins::index_type; + using ParticleType = WarpXParticleContainer::ParticleType; + using ParticleTileType = WarpXParticleContainer::ParticleTileType; + using ParticleBins = DenseBins; + using index_type = ParticleBins::index_type; /* Find the particles and count the particles that are in each cell. Note that this does *not* rearrange particle arrays */ ParticleBins - findParticlesInEachCell (int lev, - MFIter const & mfi, - ParticleTileType & ptile) { + findParticlesInEachCell( int const lev, MFIter const& mfi, + ParticleTileType const& ptile) { // Extract particle structures for this tile int const np = ptile.numParticles(); - auto ptd = ptile.getParticleTileData(); + ParticleType const* particle_ptr = ptile.GetArrayOfStructs()().data(); // Extract box properties Geometry const& geom = WarpX::GetInstance().Geom(lev); @@ -55,9 +51,9 @@ namespace ParticleUtils // Find particles that are in each cell; // results are stored in the object `bins`. ParticleBins bins; - bins.build(np, ptd, cbx, + bins.build(np, particle_ptr, cbx, // Pass lambda function that returns the cell index - [=] AMREX_GPU_DEVICE (ParticleType const & p) noexcept -> amrex::IntVect + [=] AMREX_GPU_DEVICE (const ParticleType& p) noexcept { return IntVect{AMREX_D_DECL( static_cast((p.pos(0)-plo[0])*dxi[0] - lo.x), @@ -68,4 +64,4 @@ namespace ParticleUtils return bins; } -} // namespace ParticleUtils +} diff --git a/Source/ablastr/particles/IndexHandling.H b/Source/ablastr/particles/IndexHandling.H new file mode 100644 index 00000000000..0ad5ca60446 --- /dev/null +++ b/Source/ablastr/particles/IndexHandling.H @@ -0,0 +1,41 @@ +/* Copyright 2019-2022 Axel Huebl + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef ABLASTR_INDEX_HANDLING_H +#define ABLASTR_INDEX_HANDLING_H + +#include + + +namespace ablastr::particles { + + /** A helper function to derive a globally unique particle ID + * + * @param[in] id AMReX particle ID (on local cpu/rank), AoS .id + * @param[in] cpu AMReX particle CPU (rank) at creation of the particle, AoS .cpu + * @return global particle ID that is unique and permanent in the whole simulation + */ + constexpr uint64_t + localIDtoGlobal (int const id, int const cpu) + { + static_assert(sizeof(int) * 2u <= sizeof(uint64_t), + "int size might cause collisions in global IDs"); + // implementation: + // - we cast both 32-bit (or smaller) ints to a 64bit unsigned int + // - this will leave half of the "upper range" bits in the 64bit unsigned int zeroed out + // because the corresponding (extended) value range was not part of the value range in + // the int representation + // - we bit-shift the cpu into the upper half of zero bits in the 64 bit unsigned int + // (imagine this step as "adding a per-cpu/rank offset to the local integers") + // - then we add this offset + // note: the add is expressed as bitwise OR (|) since this saves us from writing + // brackets for operator precedence between + and << + return uint64_t(id) | uint64_t(cpu) << 32u; + } + +} // namespace ablastr::particles + +#endif // ABLASTR_INDEX_HANDLING_H diff --git a/Source/ablastr/particles/ParticleMoments.H b/Source/ablastr/particles/ParticleMoments.H index b648ccb28aa..e45fb574cce 100644 --- a/Source/ablastr/particles/ParticleMoments.H +++ b/Source/ablastr/particles/ParticleMoments.H @@ -35,7 +35,7 @@ namespace particles { amrex::ParticleReal, amrex::ParticleReal> MinAndMaxPositions (T_PC const & pc) { - using ConstParticleTileDataType = typename T_PC::ParticleTileType::ConstParticleTileDataType; + using PType = typename T_PC::SuperParticleType; // Get min and max for the local rank amrex::ReduceOps< @@ -46,11 +46,11 @@ namespace particles { amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal> >( pc, - [=] AMREX_GPU_DEVICE(const ConstParticleTileDataType& ptd, const int i) noexcept + [=] AMREX_GPU_DEVICE(PType const & p) noexcept { - const amrex::ParticleReal x = ptd.rdata(0)[i]; - const amrex::ParticleReal y = ptd.rdata(1)[i]; - const amrex::ParticleReal z = ptd.rdata(2)[i]; + amrex::ParticleReal const x = p.pos(0); + amrex::ParticleReal const y = p.pos(1); + amrex::ParticleReal const z = p.pos(2); return amrex::makeTuple(x, y, z, x, y, z); }, @@ -90,8 +90,7 @@ namespace particles { amrex::ParticleReal, amrex::ParticleReal> MeanAndStdPositions (T_PC const & pc) { - - using ConstParticleTileDataType = typename T_PC::ParticleTileType::ConstParticleTileDataType; + using PType = typename T_PC::SuperParticleType; amrex::ReduceOps< amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum, @@ -104,14 +103,12 @@ namespace particles { amrex::ParticleReal> >( pc, - [=] AMREX_GPU_DEVICE(const ConstParticleTileDataType& ptd, const int i) noexcept + [=] AMREX_GPU_DEVICE(const PType& p) noexcept { - - const amrex::ParticleReal x = ptd.rdata(0)[i]; - const amrex::ParticleReal y = ptd.rdata(1)[i]; - const amrex::ParticleReal z = ptd.rdata(2)[i]; - - const amrex::ParticleReal w = ptd.rdata(T_RealSoAWeight)[i]; + amrex::ParticleReal const x = p.pos(0); + amrex::ParticleReal const y = p.pos(1); + amrex::ParticleReal const z = p.pos(2); + amrex::ParticleReal const w = p.rdata(T_RealSoAWeight); return amrex::makeTuple(x, x*x, y, y*y, z, z*z, w); }, From fa422b501814a8beb3e3a7ac2e13bc176e010fac Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Wed, 31 Jan 2024 11:25:37 -0800 Subject: [PATCH 9/9] AMReX: Latest `development` (#4654) --- .github/workflows/cuda.yml | 2 +- Regression/WarpX-GPU-tests.ini | 2 +- Regression/WarpX-tests.ini | 2 +- cmake/dependencies/AMReX.cmake | 2 +- run_test.sh | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/.github/workflows/cuda.yml b/.github/workflows/cuda.yml index eca39369a08..9960cdfbb29 100644 --- a/.github/workflows/cuda.yml +++ b/.github/workflows/cuda.yml @@ -115,7 +115,7 @@ jobs: which nvcc || echo "nvcc not in PATH!" git clone https://github.com/AMReX-Codes/amrex.git ../amrex - cd ../amrex && git checkout --detach ab99ea69089f4fffbfd727ed965d5ceb3d905baa && cd - + cd ../amrex && git checkout --detach 689144d157a0106faf3d0ae89f8d90b0250cf975 && cd - make COMP=gcc QED=FALSE USE_MPI=TRUE USE_GPU=TRUE USE_OMP=FALSE USE_PSATD=TRUE USE_CCACHE=TRUE -j 2 ccache -s diff --git a/Regression/WarpX-GPU-tests.ini b/Regression/WarpX-GPU-tests.ini index 51ae03b60fa..70c62b190fb 100644 --- a/Regression/WarpX-GPU-tests.ini +++ b/Regression/WarpX-GPU-tests.ini @@ -60,7 +60,7 @@ emailBody = Check https://ccse.lbl.gov/pub/GpuRegressionTesting/WarpX/ for more [AMReX] dir = /home/regtester/git/amrex/ -branch = ab99ea69089f4fffbfd727ed965d5ceb3d905baa +branch = 689144d157a0106faf3d0ae89f8d90b0250cf975 [source] dir = /home/regtester/git/WarpX diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 3489e8ad100..ab11d70dfcc 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -59,7 +59,7 @@ emailBody = Check https://ccse.lbl.gov/pub/RegressionTesting/WarpX/ for more det [AMReX] dir = /home/regtester/AMReX_RegTesting/amrex/ -branch = ab99ea69089f4fffbfd727ed965d5ceb3d905baa +branch = 689144d157a0106faf3d0ae89f8d90b0250cf975 [source] dir = /home/regtester/AMReX_RegTesting/warpx diff --git a/cmake/dependencies/AMReX.cmake b/cmake/dependencies/AMReX.cmake index 72173b2acb0..9b74c8db7fd 100644 --- a/cmake/dependencies/AMReX.cmake +++ b/cmake/dependencies/AMReX.cmake @@ -269,7 +269,7 @@ set(WarpX_amrex_src "" set(WarpX_amrex_repo "https://github.com/AMReX-Codes/amrex.git" CACHE STRING "Repository URI to pull and build AMReX from if(WarpX_amrex_internal)") -set(WarpX_amrex_branch "ab99ea69089f4fffbfd727ed965d5ceb3d905baa" +set(WarpX_amrex_branch "689144d157a0106faf3d0ae89f8d90b0250cf975" CACHE STRING "Repository branch for WarpX_amrex_repo if(WarpX_amrex_internal)") diff --git a/run_test.sh b/run_test.sh index 20c6638f4d0..48857d264cb 100755 --- a/run_test.sh +++ b/run_test.sh @@ -68,7 +68,7 @@ python3 -m pip install --upgrade -r warpx/Regression/requirements.txt # Clone AMReX and warpx-data git clone https://github.com/AMReX-Codes/amrex.git -cd amrex && git checkout --detach ab99ea69089f4fffbfd727ed965d5ceb3d905baa && cd - +cd amrex && git checkout --detach 689144d157a0106faf3d0ae89f8d90b0250cf975 && cd - # warpx-data contains various required data sets git clone --depth 1 https://github.com/ECP-WarpX/warpx-data.git # openPMD-example-datasets contains various required data sets