AutoPas  3.0.0
Loading...
Searching...
No Matches
VLCCellPairNeighborList.h
Go to the documentation of this file.
1
7#pragma once
10
11namespace autopas {
12
18template <class ParticleCell>
19class TraversalSelector;
20
21template <class Particle_T>
22class VLCCellPairTraversalInterface;
29template <class Particle_T>
31 public:
36
40 using SoAPairOfParticleAndList = std::pair<size_t, std::vector<size_t, autopas::AlignedAllocator<size_t>>>;
41
46 using SoAListType = typename std::vector<std::vector<std::vector<SoAPairOfParticleAndList>>>;
47
48 [[nodiscard]] ContainerOption getContainerType() const override { return ContainerOption::pairwiseVerletLists; }
49
50 size_t getNumberOfPartners(const Particle_T *particle) const override {
51 size_t listSize{0};
52 bool particleFound{false};
53 for (auto &neighborListsForOneCell : _aosNeighborList) {
54 for (const auto &neighborListsForCellCellPair : neighborListsForOneCell) {
55 // Check for the desired particle in the first cell. We check all "partner" cells because the particle might not
56 // have any particles from one "partner" cell in its list.
57 for (size_t i{0}; i < neighborListsForCellCellPair.size(); ++i) {
58 if (neighborListsForCellCellPair[i].first == particle) {
59 particleFound = true;
60 // We accumulate the number of partners. Partners can be in all lists of neighboring cells (middle for loop)
61 listSize += neighborListsForCellCellPair[i].second.size();
62 // Since we found the particle, we can skip all the other particles in the current list
63 break;
64 }
65 }
66 }
67 if (particleFound) {
68 // We've found the particle in the cell with neighbor lists neighborListsForOneCell. So we can skip all the
69 // other cells (outer for loop)
70 return listSize;
71 }
72 }
73 return 0lu;
74 }
75
81 return _aosNeighborList;
82 }
83
88 auto &getSoANeighborList() { return _soaNeighborList; }
89
93 void buildAoSNeighborList(TraversalOption vlcTraversalOpt, LinkedCells<Particle_T> &linkedCells, bool useNewton3) {
94 using namespace utils::ArrayMath::literals;
95 // Sanity check.
96 if (linkedCells.getCellBlock().getCellsPerInteractionLength() > 1) {
98 "VLCCellPairNeighborList::buildAoSNeighborList() was called with a CSF < 1 but it only supports CSF>=1.");
99 }
100 // Define some aliases
101 auto &neighborLists = getAoSNeighborList();
102 auto &cells = linkedCells.getCells();
103 const auto interactionLength = linkedCells.getInteractionLength();
104 const auto interactionLengthSquared = interactionLength * interactionLength;
105 const auto boxSizeWithHalo = linkedCells.getBoxMax() - linkedCells.getBoxMin() +
106 std::array<double, 3>{interactionLength, interactionLength, interactionLength} * 2.;
107
108 // Helper lambda to compute the relative index from two cells within in a 3x3x3 cell-cube
109 auto relativeNeighborhoodIndex = [&](auto cellIndex1, auto cellIndex2) {
110 const auto cellsPerDimensionWithHalo = linkedCells.getCellBlock().getCellsPerDimensionWithHalo();
111 const auto threeDPosCell1 = utils::ThreeDimensionalMapping::oneToThreeD(static_cast<long unsigned>(cellIndex1),
112 cellsPerDimensionWithHalo);
113 const auto threeDPosCell2 = utils::ThreeDimensionalMapping::oneToThreeD(static_cast<long unsigned>(cellIndex2),
114 cellsPerDimensionWithHalo);
115 const auto offset = threeDPosCell2 - threeDPosCell1;
116 return (offset[0] + 1) * 9 + (offset[1] + 1) * 3 + (offset[2] + 1);
117 };
118
119 // This assumes homogeneous distribution and some overestimation.
120 const auto listLengthEstimate = VerletListsCellsHelpers::estimateListLength(
121 linkedCells.getNumberOfParticles(IteratorBehavior::ownedOrHalo), boxSizeWithHalo, interactionLength, 1.3);
122
123 // Reset lists. Don't free any memory, only mark as unused.
124 this->setLinkedCellsPointer(&linkedCells);
125 for (auto &neighborCellLists : neighborLists) {
126 for (auto &cellLists : neighborCellLists) {
127 for (auto &[particlePtr, neighbors] : cellLists) {
128 particlePtr = nullptr;
129 neighbors.clear();
130 }
131 }
132 }
133 neighborLists.resize(cells.size());
134 // each cell hast 27 neighbors including itself
135 for (auto &neighborList : neighborLists) {
136 neighborList.resize(27);
137 }
138
139 /* This must not be a doc comment (with two **) to not confuse doxygen.
140 * Helper function to insert a pointer into a list of the base cell.
141 * It considers the cases that neither particle is in the base cell
142 * and in that case finds or creates the appropriate list.
143 *
144 * @param p1 Reference to source particle.
145 * @param p2 Reference to target particle.
146 * @param neighborList Reference to the list where the particle pair should be stored.
147 */
148 auto insert = [&](auto &p1, auto &p2, auto &neighborList) {
149 // Check if the base cell already has a list for p1
150 auto iter = std::find_if(neighborList.begin(), neighborList.end(), [&](const auto &pair) {
151 const auto &[particlePtr, list] = pair;
152 return particlePtr == &p1;
153 });
154 // If yes, append p2 to it.
155 if (iter != neighborList.end()) {
156 iter->second.push_back(&p2);
157 } else {
158 // If no, create one (or reuse an empty pair), reserve space for the list and emplace p2
159 if (auto insertLoc = std::find_if(neighborList.begin(), neighborList.end(),
160 [&](const auto &pair) {
161 const auto &[particlePtr, list] = pair;
162 return particlePtr == nullptr;
163 });
164 insertLoc != neighborList.end()) {
165 auto &[particlePtr, neighbors] = *insertLoc;
166 particlePtr = &p1;
167 neighbors.reserve(listLengthEstimate);
168 neighbors.push_back(&p2);
169 } else {
170 neighborList.emplace_back(&p1, std::vector<Particle_T *>{});
171 neighborList.back().second.reserve(listLengthEstimate);
172 neighborList.back().second.push_back(&p2);
173 }
174 }
175 };
176
177 const auto &cellsPerDim =
178 utils::ArrayUtils::static_cast_copy_array<int>(linkedCells.getCellBlock().getCellsPerDimensionWithHalo());
179 // Vector of offsets from the base cell for the given base step
180 // and respective factors for the fraction of particles per cell that need neighbor lists in the base cell.
181 const auto offsets = VerletListsCellsHelpers::buildBaseStep(cellsPerDim, vlcTraversalOpt);
182
183 int xEnd{};
184 int yEnd{};
185 int zEnd{};
186
187 switch (vlcTraversalOpt) {
188 case TraversalOption::vlc_c08:
189 // Go over all cells except the very last layer and create lists per base step.
190 xEnd = cellsPerDim[0] - 1;
191 yEnd = cellsPerDim[1] - 1;
192 zEnd = cellsPerDim[2] - 1;
193 break;
194 default:
195 xEnd = cellsPerDim[0];
196 yEnd = cellsPerDim[1];
197 zEnd = cellsPerDim[2];
198 break;
199 }
200
201 // Since there are no loop dependencies merge all for loops and create 10 chunks per thread.
202 AUTOPAS_OPENMP(parallel for collapse(3) schedule(dynamic, std::max(cells.size() / (autopas::autopas_get_max_threads() * 10), 1ul)))
203 for (int z = 0; z < zEnd; ++z) {
204 for (int y = 0; y < yEnd; ++y) {
205 for (int x = 0; x < xEnd; ++x) {
206 // aliases
207 const auto cellIndexBase = utils::ThreeDimensionalMapping::threeToOneD(x, y, z, cellsPerDim);
208 auto &baseCell = cells[cellIndexBase];
209 auto &baseCellsLists = neighborLists[cellIndexBase];
210 auto threadNum = autopas_get_thread_num();
211
212 // Build lists for this base step according to predefined cell pairs
213 for (const auto &[offset1, offset2, _] : offsets) {
214 const auto cell1Index1D = cellIndexBase + offset1;
215 const auto cell2Index1D = cellIndexBase + offset2;
216 auto &cell1List = neighborLists[cell1Index1D];
217 auto &cell2List = neighborLists[cell2Index1D];
218
219 // For all traversals ensures the partner cell is not outside the boundary
220 const auto cell2Index3D = utils::ThreeDimensionalMapping::oneToThreeD(cell2Index1D, cellsPerDim);
221 if (cell2Index3D[0] >= cellsPerDim[0] or cell2Index3D[0] < 0 or cell2Index3D[1] >= cellsPerDim[1] or
222 cell2Index3D[1] < 0 or cell2Index3D[2] >= cellsPerDim[2] or cell2Index3D[2] < 0) {
223 continue;
224 }
225
226 // Skip if both cells only contain halos or dummys
227 if (not(cells[cell1Index1D].getPossibleParticleOwnerships() == OwnershipState::owned) and
228 not(cells[cell2Index1D].getPossibleParticleOwnerships() == OwnershipState::owned)) {
229 continue;
230 }
231
232 // Go over all particle pairs in the two cells and insert close pairs into their respective lists
233 for (size_t particleIndexCell1 = 0; particleIndexCell1 < cells[cell1Index1D].size(); ++particleIndexCell1) {
234 auto &p1 = cells[cell1Index1D][particleIndexCell1];
235
236 // Determine the starting index for the second particle in the pair.
237 // If both cells are the same and the traversal is newton3 compatible, this is the next particle in the
238 // cell. If the cells are different or the traversal is not newton3 compatible, this is index 0. For
239 // non-newton 3 compatible traversals (vlp_c01) we have to consider the interaction in both directions, so
240 // we always start from index 0.
241 size_t startIndexCell2 = 0;
242 if (cell1Index1D == cell2Index1D and vlcTraversalOpt != TraversalOption::vlp_c01) {
243 startIndexCell2 = particleIndexCell1 + 1;
244 }
245
246 for (size_t particleIndexCell2 = startIndexCell2; particleIndexCell2 < cells[cell2Index1D].size();
247 ++particleIndexCell2) {
248 auto &p2 = cells[cell2Index1D][particleIndexCell2];
249 // Ignore dummies and self interaction
250 if (&p1 == &p2 or p1.isDummy() or p2.isDummy()) {
251 continue;
252 }
253
254 // If the distance is less than interaction length add the pair to the list
255 const auto distVec = p2.getR() - p1.getR();
256 const auto distSquared = utils::ArrayMath::dot(distVec, distVec);
257 if (distSquared < interactionLengthSquared) {
258 {
259 const size_t secondCellIndexInFirst = relativeNeighborhoodIndex(cell1Index1D, cell2Index1D);
260 insert(p1, p2, cell1List[secondCellIndexInFirst]);
261 }
262 // If the traversal is Newton3 compatible, Newton3 is used for building regardless of if the actual
263 // interactions will use Newton3. In the case that Newton3 will not be used, the inverse interaction
264 // also needs to be stored in p2's list. If the traversal is Newton3 incompatible, the insertion into
265 // p2's list will occur when the p2 particle is p1. This is ensured above by the startIndexCell2.
266 if (not useNewton3 and not(vlcTraversalOpt == TraversalOption::vlp_c01)) {
267 {
268 const size_t secondCellIndexInFirst = relativeNeighborhoodIndex(cell2Index1D, cell1Index1D);
269 insert(p2, p1, cell2List[secondCellIndexInFirst]);
270 }
271 }
272 }
273 }
274 }
275 }
276 }
277 }
278 }
279
280 // Cleanup: Remove any unused ptr-list pairs to avoid accessing nullptr
281 for (auto &neighborCellLists : neighborLists) {
282 for (auto &cellLists : neighborCellLists) {
283 cellLists.erase(std::remove_if(cellLists.begin(), cellLists.end(),
284 [](const auto &pair) {
285 const auto &[particlePtr, neighbors] = pair;
286 return particlePtr == nullptr;
287 }),
288 cellLists.end());
289 }
290 }
291 }
292
293 void generateSoAFromAoS(LinkedCells<Particle_T> &linkedCells) override {
294 _soaNeighborList.clear();
295
296 // particle pointer to global index of particle
297 std::unordered_map<Particle_T *, size_t> particlePtrToIndex;
298 particlePtrToIndex.reserve(linkedCells.size());
299 size_t i = 0;
300 for (auto iter = linkedCells.begin(IteratorBehavior::ownedOrHaloOrDummy); iter.isValid(); ++iter, ++i) {
301 particlePtrToIndex[&(*iter)] = i;
302 }
303
304 _soaNeighborList.resize(linkedCells.getCells().size());
305
306 // iterate over cells and for each create the soa lists from the aos lists
307 for (size_t firstCellIndex = 0; firstCellIndex < _aosNeighborList.size(); ++firstCellIndex) {
308 const auto &aosLists = _aosNeighborList[firstCellIndex];
309 auto &soaLists = _soaNeighborList[firstCellIndex];
310 soaLists.resize(aosLists.size());
311
312 // iterate over each cell's neighboring cells
313 for (size_t secondCellIndex = 0; secondCellIndex < aosLists.size(); ++secondCellIndex) {
314 const auto &aosCellPairLists = aosLists[secondCellIndex];
315 auto &soaCellPairLists = soaLists[secondCellIndex];
316 soaCellPairLists.reserve(aosCellPairLists.capacity());
317
318 // iterate over pairs of particle and neighbor list
319 for (const auto &[particlePtr, neighbors] : aosCellPairLists) {
320 // global index of current particle
321 size_t currentParticleGlobalIndex = particlePtrToIndex.at(particlePtr);
322
323 // create SoA neighbor list for current particle
324 std::vector<size_t, autopas::AlignedAllocator<size_t>> currentSoANeighborList{};
325 currentSoANeighborList.reserve(neighbors.size());
326
327 // fill the SoA neighbor list with the indices of the particles from the corresponding AoS neighbor list
328 for (const auto &neighborOfCurrentParticle : neighbors) {
329 currentSoANeighborList.emplace_back(particlePtrToIndex.at(neighborOfCurrentParticle));
330 }
331
332 // add the newly constructed pair of particle index and SoA particle neighbor list to cell pair
333 soaCellPairLists.emplace_back(currentParticleGlobalIndex, currentSoANeighborList);
334 }
335 }
336 }
337 }
338
339 void setUpTraversal(TraversalInterface *traversal) override {
340 auto vTraversal = dynamic_cast<VLCCellPairTraversalInterface<Particle_T> *>(traversal);
341
342 if (vTraversal) {
343 vTraversal->setVerletList(*this);
344 } else {
345 auto traversal2 =
347 if (traversal2) {
348 traversal2->setVerletList(*this);
349 } else {
351 "Trying to use a traversal of wrong type in VerletListCells.h. TraversalID: {}",
352 traversal->getTraversalType());
353 }
354 }
355 }
356
357 private:
362 std::vector<std::vector<std::vector<std::pair<Particle_T *, std::vector<Particle_T *>>>>>();
363
368 SoAListType _soaNeighborList = std::vector<
369 std::vector<std::vector<std::pair<size_t, std::vector<size_t, autopas::AlignedAllocator<size_t>>>>>>();
370};
371} // namespace autopas
#define AUTOPAS_OPENMP(args)
Empty macro to throw away any arguments.
Definition: WrapOpenMP.h:126
double getInteractionLength() const final
Return the interaction length (cutoff+skin) of the container.
Definition: CellBasedParticleContainer.h:89
size_t getNumberOfParticles(IteratorBehavior behavior) const override
Get the number of particles with respect to the specified IteratorBehavior.
Definition: CellBasedParticleContainer.h:113
size_t size() const override
Get the total number of particles saved in the container (owned + halo + dummy).
Definition: CellBasedParticleContainer.h:131
const std::array< double, 3 > & getBoxMin() const final
Get the lower corner of the container without halo.
Definition: CellBasedParticleContainer.h:74
const std::array< double, 3 > & getBoxMax() const final
Get the upper corner of the container without halo.
Definition: CellBasedParticleContainer.h:69
LinkedCells class.
Definition: LinkedCells.h:40
internal::CellBlock3D< ParticleCell > & getCellBlock()
Get the cell block, not supposed to be used except by verlet lists.
Definition: LinkedCells.h:458
ContainerIterator< ParticleType, true, false > begin(IteratorBehavior behavior=autopas::IteratorBehavior::ownedOrHalo, typename ContainerIterator< ParticleType, true, false >::ParticleVecType *additionalVectors=nullptr) override
Iterate over all particles using for(auto iter = container.begin(); iter.isValid(); ++iter) .
Definition: LinkedCells.h:307
std::vector< ParticleCell > & getCells()
Returns a non-const reference to the cell data structure.
Definition: LinkedCells.h:470
This interface serves as a common parent class for all traversals.
Definition: TraversalInterface.h:18
virtual TraversalOption getTraversalType() const =0
Return a enum representing the name of the traversal class.
Neighbor list to be used with VerletListsCells container.
Definition: VLCCellPairNeighborList.h:30
auto & getSoANeighborList()
Returns the neighbor list in SoA layout.
Definition: VLCCellPairNeighborList.h:88
ContainerOption getContainerType() const override
Returns the container type of this neighbor list and the container it belongs to.
Definition: VLCCellPairNeighborList.h:48
typename std::vector< std::vector< std::vector< SoAPairOfParticleAndList > > > SoAListType
Helper type definition.
Definition: VLCCellPairNeighborList.h:46
void buildAoSNeighborList(TraversalOption vlcTraversalOpt, LinkedCells< Particle_T > &linkedCells, bool useNewton3)
Builds AoS neighbor list from underlying linked cells object.
Definition: VLCCellPairNeighborList.h:93
std::pair< size_t, std::vector< size_t, autopas::AlignedAllocator< size_t > > > SoAPairOfParticleAndList
Helper type definition.
Definition: VLCCellPairNeighborList.h:40
VerletListsCellsHelpers::PairwiseNeighborListsType< Particle_T > & getAoSNeighborList()
Returns the neighbor list in AoS layout.
Definition: VLCCellPairNeighborList.h:80
size_t getNumberOfPartners(const Particle_T *particle) const override
Gets the number of neighbors over all neighbor lists that belong to this particle.
Definition: VLCCellPairNeighborList.h:50
void setUpTraversal(TraversalInterface *traversal) override
Assigns the current traversal to the correct traversal interface.
Definition: VLCCellPairNeighborList.h:339
void generateSoAFromAoS(LinkedCells< Particle_T > &linkedCells) override
Generates neighbor list in SoA layout from available neighbor list in AoS layout.
Definition: VLCCellPairNeighborList.h:293
typename VerletListsCellsHelpers::PairwiseNeighborListsType< Particle_T > ListType
Type of the data structure used to save the neighbor lists.
Definition: VLCCellPairNeighborList.h:35
Interface for traversals used with VLCCellPairNeighborList.
Definition: VLCCellPairTraversalInterface.h:16
void setVerletList(VLCCellPairNeighborList< Particle_T > &verlet)
Sets the verlet list for the traversal to iterate over.
Definition: VLCCellPairTraversalInterface.h:22
Interface of neighbor lists to be used with VerletListsCells container.
Definition: VLCNeighborListInterface.h:19
void setLinkedCellsPointer(LinkedCells< Particle_T > *linkedCells)
Set the Linked Cells Pointer for this List.
Definition: VLCNeighborListInterface.h:110
This class provides the Traversal Interface for the verlet lists cells container.
Definition: VLCTraversalInterface.h:27
virtual void setVerletList(NeighborList &verlet)
Sets the verlet list for the traversal to iterate over.
Definition: VLCTraversalInterface.h:41
static void exception(const Exception e)
Handle an exception derived by std::exception.
Definition: ExceptionHandler.h:63
std::vector< std::vector< std::vector< std::pair< Particle_T *, std::vector< Particle_T * > > > > > PairwiseNeighborListsType
Pairwise verlet lists: For every cell a vector, for every neighboring cell a vector of particle-neigh...
Definition: VerletListsCellsHelpers.h:41
size_t estimateListLength(size_t numParticles, const std::array< double, 3 > &boxSize, double interactionLength, double correctionFactor)
Simple heuristic to calculate the average number of particles per verlet list assuming particles are ...
Definition: VerletListsCellsHelpers.cpp:16
std::vector< BaseStepOffsets > buildBaseStep(const std::array< int, 3 > &cellsPerDim, const TraversalOption traversal)
Builds the list of offsets from the base cell for the c01, c08, and c18 base step.
Definition: VerletListsCellsHelpers.cpp:26
constexpr T dot(const std::array< T, SIZE > &a, const std::array< T, SIZE > &b)
Generates the dot product of two arrays.
Definition: ArrayMath.h:233
constexpr std::array< T, 3 > oneToThreeD(T ind, const std::array< T, 3 > &dims)
Convert a 1d index to a 3d index.
Definition: ThreeDimensionalMapping.h:55
constexpr T threeToOneD(T x, T y, T z, const std::array< T, 3 > &dims)
Convert a 3d index to a 1d index.
Definition: ThreeDimensionalMapping.h:29
This is the main namespace of AutoPas.
Definition: AutoPasDecl.h:32
int autopas_get_max_threads()
Dummy for omp_get_max_threads() when no OpenMP is available.
Definition: WrapOpenMP.h:144
@ owned
Owned state, a particle with this state is an actual particle and owned by the current AutoPas object...
int autopas_get_thread_num()
Dummy for omp_set_lock() when no OpenMP is available.
Definition: WrapOpenMP.h:132