Line data Source code
1 : #ifndef _MOLECULARDYNAMICS_MOLECULARCONTAINER_H_
2 : #define _MOLECULARDYNAMICS_MOLECULARCONTAINER_H_
3 :
4 : #include <Kokkos_Core.hpp>
5 :
6 : #include <simplemd/MolecularDynamicsDefinitions.h>
7 : #include <simplemd/Molecule.h>
8 : #include <simplemd/LinkedCell.h>
9 : #include <simplemd/services/ParallelTopologyService.h>
10 : #include <simplemd/services/MolecularPropertiesService.h>
11 :
12 : namespace simplemd {
13 : namespace services {
14 : // forward declarations to remove circular dependencies
15 : class ParallelTopologyService;
16 : } // namespace services
17 :
18 : class MoleculeContainer;
19 : } // namespace simplemd
20 :
21 : /**
22 : * @brief Contains molecules, creates linked cells and manages its own memory.
23 : *
24 : * The MoleculeContainer class serves as the sole container responsible for all molecule-related work.
25 : * Molecules are stored in a 2D Kokkos View, where each row corresponds to a linked cell (ghost included).
26 : * While creating a container, the maximum possible capacity (in molecules) of any cell is defined by cellCapacity.
27 : * Thus, each linked cell has a numMolecules() <= capacity.
28 : * All data with an index between numMolecules and capacity is treated as garbage data.
29 : */
30 : class simplemd::MoleculeContainer {
31 : public:
32 : /**
33 : * @brief Construct a new MoleculeContainer object
34 : *
35 : * Uses parallelTopologyService to store relevant information related to spatial layout,
36 : * so that 3D coordinated can be converted to 3D and 1D linked cell indices.
37 : *
38 : * @param parallelTopologyService Used to extract and store local number of cells, domain offset etc.
39 : * @param cellCapacity The maximum capacity of any cell. Cannot be changed, must provide ample room at compiletime.
40 : */
41 : MoleculeContainer(simplemd::services::ParallelTopologyService& parallelTopologyService, int cellCapacity);
42 :
43 : /**
44 : * @brief Inserts a molecule into a specific linked cell.
45 : *
46 : * @param cellIdx The one-dimensional index of the linked cell to insert the molecule into (ghost included).
47 : * @param molecule The molecule to be inserted.
48 : */
49 : void insert(unsigned int cellIdx, const simplemd::Molecule& molecule);
50 :
51 : /**
52 : * @brief Inserts a molecule into a specific linked cell.
53 : *
54 : * @param cellIdx The vector index of the linked cell to insert the molecule into (ghost included).
55 : * @param molecule The molecule to be inserted.
56 : */
57 : void insert(tarch::la::Vector<MD_DIM, unsigned int> cellIdx, const simplemd::Molecule& molecule) { insert(vectorIndexToLinear(cellIdx), molecule); }
58 :
59 : /**
60 : * @brief Inserts a molecule into the container.
61 : *
62 : * The function calculates the appropriate linked cell to insert the molecule into.
63 : *
64 : * @param molecule The molecule to be inserted.
65 : */
66 : void insert(const simplemd::Molecule& molecule);
67 :
68 : /**
69 : * @brief Removes a molecule, given its ID and linked cell index.
70 : *
71 : * This is done by swapping the last molecule of the linked cell into position moleculeIdx,
72 : * and decrementing the number of molecules in this cell.
73 : *
74 : * @param cellIdx 1D linked cell index of the molecule (ghost included).
75 : * @param moleculeIdx Index of the molecule within the linked cell.
76 : */
77 : void remove(unsigned int cellIdx, unsigned int moleculeIdx);
78 :
79 : /**
80 : * @brief Clears a linked cell.
81 : *
82 : * This function simply sets the number of molecules in the linked cell as 0.
83 : * All data in the linked cell is treated as garbage data, and insertion starts from index 0
84 : * within the linked cell again, overwriting previous data.
85 : *
86 : * @param cellIdx 1D index of the linked cell to clear (ghost included).
87 : */
88 : void clearLinkedCell(unsigned int cellIdx);
89 :
90 : /**
91 : * @brief Removes all outgoing molecules from a linked cell and moves them to the appropriate destination.
92 : *
93 : * @param cellIdx 1D index of the linked cell to sort (ghost included).
94 : */
95 : void sort(unsigned int cellIdx);
96 :
97 : /**
98 : * @brief Puts all molecules into their appropriate linked cells.
99 : *
100 : * Since a molecule can only travel to its immediate neighbours per timestep, sorting is done with a red-black
101 : * traversal of the whole domain (including ghost cells). This ensures that there are no concurrency issues since
102 : * no two linked cells should try to write to the same linked cell. Concurrency is implemented in quarter shells.
103 : * Since this includes ghost shells, they should be empty before sort() is called.
104 : */
105 : void sort();
106 :
107 : /**
108 : * @brief Get the Molecule at linked cell index i (ghost included), position j
109 : *
110 : * @param i
111 : * @param j
112 : * @return simplemd::Molecule&
113 : */
114 : KOKKOS_FUNCTION simplemd::Molecule& getMoleculeAt(size_t i, size_t j) const;
115 :
116 : /**
117 : * @brief Returns the linked cell at 1D index idx (ghost included)
118 : *
119 : * @param idx
120 : * @return simplemd::LinkedCell
121 : */
122 : KOKKOS_FUNCTION simplemd::LinkedCell operator[](const size_t idx) const;
123 :
124 : /**
125 : * @brief Returns the linked cell at 3D index idx (ghost included)
126 : *
127 : * @param idx
128 : * @return simplemd::LinkedCell
129 : */
130 : KOKKOS_FUNCTION simplemd::LinkedCell operator[](const tarch::la::Vector<MD_DIM, unsigned int> cellIdx) const;
131 :
132 : /**
133 : * @brief Get the total number of cells in the container
134 : *
135 : * @return size_t
136 : */
137 : KOKKOS_FUNCTION size_t getLocalNumberOfCellsScalarWithGhost() const;
138 :
139 : /**
140 : * @brief Returns the number of molecules in all cells
141 : *
142 : * @return const size_t
143 : */
144 : size_t getLocalNumberOfMoleculesWithGhost() const;
145 :
146 : /**
147 : * @brief returns the index of the first (non-ghost) cell along each dimension
148 : *
149 : * @return tarch::la::Vector
150 : */
151 : const tarch::la::Vector<MD_DIM, unsigned int>& getLocalIndexOfFirstCell() const;
152 :
153 : /**
154 : * @brief returns the number of (non-ghost) cells along each dimension
155 : *
156 : * @return tarch::la::Vector
157 : */
158 : const tarch::la::Vector<MD_DIM, unsigned int> getLocalNumberOfCells() const;
159 :
160 : /**
161 : * @brief returns the local cell index vector for the local cell index cellIndex
162 : *
163 : * @return tarch::la::Vector
164 : */
165 : KOKKOS_FUNCTION tarch::la::Vector<MD_DIM, unsigned int> getLocalCellIndexVector(const size_t cellIndex) const;
166 :
167 : /**
168 : * @brief can be used to apply a molecule-mapping which is iterated over all molecules of this process
169 : * uses static member in mapping class (A::IsParallel) to determine whether the parallel or serial iterator will be called
170 : */
171 : template <class A> void iterateMolecules(A& a);
172 :
173 : /**
174 : * @brief can be used to apply a molecule-with-cell-mapping which is iterated over all molecules of this process
175 : * uses static member in mapping class (A::IsParallel) to determine whether the parallel or serial iterator will be called
176 : */
177 : template <class A> void iterateMoleculesWithCell(A& a);
178 :
179 : /**
180 : * @brief iterates over all cells in the range defined by the lower left front
181 : * corner cell lowerLeftFrontCell and the size of the domain cellRange.
182 : * cellRange defines a number of cells in each spatial direction that the
183 : * class A shall be applied to. lowerLeftFrontCell needs to be given in local
184 : * coordinates.
185 : *
186 : * uses static member in mapping class (A::IsParallel) to determine whether the parallel or serial iterator will be called.
187 : */
188 : template <class A>
189 : void iterateCells(A& a, const tarch::la::Vector<MD_DIM, unsigned int>& lowerLeftFrontCell, const tarch::la::Vector<MD_DIM, unsigned int>& cellRange);
190 :
191 : /**
192 : * @brief iterates over all cells in the inner part (i.e. does not consider the
193 : * ghost layer)
194 : *
195 : * uses static member in mapping class (A::IsParallel) to determine whether the parallel or serial iterator will be called
196 : */
197 : template <class A> void iterateCells(A& a);
198 :
199 : /**
200 : * @brief iterates over all cell pairs for the cells in the inner region of each
201 : * local process
202 : *
203 : * uses static member in mapping class (A::IsParallel) to determine whether the parallel or serial iterator will be called
204 : */
205 : template <class A> void iterateCellPairs(A& a);
206 :
207 : /**
208 : * @brief iterates over all cell pairs cell1 and cell2 with cell1 in the range
209 : * described by lowerLeftFrontCell and cellRange; cell2 does not need to lie
210 : * within the range (example: iterate only over lowerLeftFrontCell=(1,1,1) and
211 : * cellRange=(1,1,1). Then, we will consider amongst others the pair
212 : * (0,0,0),(1,1,1)).
213 : *
214 : * uses static member in mapping class (A::IsParallel) to determine whether the parallel or serial iterator will be called
215 : */
216 : template <class A>
217 : void iterateCellPairs(A& a, const tarch::la::Vector<MD_DIM, unsigned int>& lowerLeftFrontCell, const tarch::la::Vector<MD_DIM, unsigned int>& cellRange);
218 :
219 : /**
220 : * @brief applies molecule mapping without any node-level parallelisation
221 : */
222 : template <class A> void iterateMoleculesSerial(A& a);
223 :
224 : /**
225 : * @brief applies molecule mapping while parallelising using Kokkos
226 : */
227 : template <class A> void iterateMoleculesParallel(A& a);
228 :
229 : /**
230 : * @brief applies molecule-with-cell mapping to all neighbors of cell
231 : */
232 : template <class A> KOKKOS_FUNCTION void handleCellNeighbors(A& a, Molecule& m, const LinkedCell& cell) const;
233 :
234 : /**
235 : * @brief applies molecule-with-cell mapping without any node-level parallelisation
236 : */
237 : template <class A> void iterateMoleculesWithCellSerial(A& a);
238 :
239 : /**
240 : * @brief applies molecule-with-cell mapping while parallelising using Kokkos
241 : */
242 : template <class A> void iterateMoleculesWithCellParallel(A& a);
243 :
244 : /**
245 : * @brief iterates over cells in parallel using Kokkos
246 : */
247 : template <class A>
248 : void iterateCellsParallel(A& a, const tarch::la::Vector<MD_DIM, unsigned int>& lowerLeftFrontCell, const tarch::la::Vector<MD_DIM, unsigned int>& cellRange);
249 :
250 : /**
251 : * @brief iterates over all cells in the inner part (i.e. does not consider the
252 : * ghost layer) in parallel using Kokkos
253 : */
254 : template <class A> void iterateCellsParallel(A& a);
255 :
256 : /**
257 : * @brief iterates over all cell pairs for the cells in the inner region of each
258 : * local process in parallel using Kokkos
259 : */
260 : template <class A> void iterateCellPairsParallel(A& a);
261 :
262 : /**
263 : * @brief iterates over all cell pairs cell1 and cell2 with cell1 in the range
264 : * described by lowerLeftFrontCell and cellRange in parallel using Kokkos; cell2 does not need to lie
265 : * within the range (example: iterate only over lowerLeftFrontCell=(1,1,1) and
266 : * cellRange=(1,1,1). Then, we will consider amongst others the pair
267 : * (0,0,0),(1,1,1)).
268 : */
269 : template <class A>
270 : void iterateCellPairsParallel(A& a, const tarch::la::Vector<MD_DIM, unsigned int>& lowerLeftFrontCell,
271 : const tarch::la::Vector<MD_DIM, unsigned int>& cellRange);
272 :
273 : /**
274 : * @brief iterates over cells without parallelization
275 : */
276 : template <class A>
277 : void iterateCellsSerial(A& a, const tarch::la::Vector<MD_DIM, unsigned int>& lowerLeftFrontCell, const tarch::la::Vector<MD_DIM, unsigned int>& cellRange);
278 :
279 : /**
280 : * @brief iterates over all cells in the inner part (i.e. does not consider the
281 : * ghost layer) without parallelisation */
282 : template <class A> void iterateCellsSerial(A& a);
283 :
284 : /**
285 : * @brief iterates over all cell pairs for the cells in the inner region of each
286 : * local process without parallelisation */
287 : template <class A> void iterateCellPairsSerial(A& a);
288 :
289 : /**
290 : * @brief iterates over all cell pairs cell1 and cell2 with cell1 in the range
291 : * described by lowerLeftFrontCell and cellRange without parallelisation; cell2 does not need to lie
292 : * within the range (example: iterate only over lowerLeftFrontCell=(1,1,1) and
293 : * cellRange=(1,1,1). Then, we will consider amongst others the pair
294 : * (0,0,0),(1,1,1)).
295 : */
296 : template <class A>
297 : void iterateCellPairsSerial(A& a, const tarch::la::Vector<MD_DIM, unsigned int>& lowerLeftFrontCell,
298 : const tarch::la::Vector<MD_DIM, unsigned int>& cellRange);
299 :
300 : /**
301 : * @brief Converts a global 3D spatial coordinate to a local 1D linked cell index.
302 : *
303 : * This can be used to find the linked cell that a molecule belongs to.
304 : *
305 : * @param position 3D spatial coordinate.
306 : * @return unsigned int
307 : */
308 : KOKKOS_FUNCTION unsigned int positionToCellIndex(const tarch::la::Vector<MD_DIM, double>& position) const;
309 :
310 : private:
311 : /**
312 : * @brief Converts a 3D local linked cell index into a 1D local linked cell index.
313 : *
314 : * @param vectorIndex 3D local index of the linked cell.
315 : * @return const unsigned int
316 : */
317 : KOKKOS_FUNCTION size_t vectorIndexToLinear(const tarch::la::Vector<MD_DIM, unsigned int>& vectorIndex) const;
318 :
319 : /**
320 : * @brief returns true if the local cell index cellIndex describes a linked cell
321 : * within the ghost layer
322 : *
323 : * @param cellIndex The linear cell index
324 : * @return bool
325 : */
326 : KOKKOS_INLINE_FUNCTION bool isGhostCell(const size_t cellIndex) const;
327 :
328 : KOKKOS_FUNCTION void printCellMolecules(size_t cellIndex) const;
329 : KOKKOS_FUNCTION void printNonGhostCells(bool shouldPrintCells, const char* const label) const;
330 :
331 : /** number of cells per direction in the local domain */
332 : const tarch::la::Vector<MD_DIM, unsigned int> _numCells;
333 : /** The number of ghost cells around the local domain along each axis on each side*/
334 : const tarch::la::Vector<MD_DIM, unsigned int> _ghostCellLayerThickness;
335 : const tarch::la::Vector<MD_DIM, unsigned int> _numLocalCellsNoGhost;
336 :
337 : /** maximum number of particles a cell (a row of the view) can contain
338 : * if this is exceeded when writing to cell, the simulation behaviour is undefined
339 : */
340 : size_t _cellCapacity;
341 :
342 : /** domain size */
343 : const tarch::la::Vector<MD_DIM, double> _domainSize;
344 :
345 : /** global domain offset */
346 : const tarch::la::Vector<MD_DIM, double> _domainOffset;
347 :
348 : /** mesh width of the linked cells */
349 : const tarch::la::Vector<MD_DIM, double> _meshWidth;
350 :
351 : /** global index of the first cell of this domain */
352 : const tarch::la::Vector<MD_DIM, unsigned int> _globalIndexOfFirstCell;
353 :
354 : /** local index of the first cell within this domain */
355 : const tarch::la::Vector<MD_DIM, unsigned int> _localIndexOfFirstCell;
356 :
357 : Kokkos::View<simplemd::Molecule**, Kokkos::LayoutRight, Kokkos::SharedSpace> _moleculeData;
358 : Kokkos::View<size_t*, Kokkos::LayoutRight, Kokkos::SharedSpace> _linkedCellNumMolecules;
359 : Kokkos::View<bool*> _linkedCellIsGhostCell;
360 : /** index offsets of all 26 neighbor cell directions */
361 : Kokkos::View<int*> _neighborOffsets;
362 : };
363 :
364 280 : template <class A> void simplemd::MoleculeContainer::iterateMolecules(A& a) {
365 : if constexpr (A::IsParallel) {
366 144 : iterateMoleculesParallel(a);
367 : } else {
368 136 : iterateMoleculesSerial(a);
369 : }
370 280 : }
371 :
372 136 : template <class A> void simplemd::MoleculeContainer::iterateMoleculesSerial(A& a) {
373 124 : a.beginMoleculeIteration();
374 1378728 : for (unsigned int i = 0; i < _linkedCellNumMolecules.size(); i++) {
375 1394692 : for (unsigned int j = 0; j < _linkedCellNumMolecules(i); j++) {
376 : #if (MD_DEBUG == MD_YES)
377 : std::cout << "Handle molecule " << j << " in cell #" << i << std::endl;
378 : #endif
379 16700 : a.handleMolecule(getMoleculeAt(i, j));
380 : }
381 : }
382 124 : a.endMoleculeIteration();
383 136 : }
384 :
385 144 : template <class A> void simplemd::MoleculeContainer::iterateMoleculesParallel(A& a) {
386 : #if (MD_DEBUG == MD_YES)
387 : iterateMoleculesSerial(a);
388 : #else
389 : a.beginMoleculeIteration();
390 144 : printNonGhostCells(true, "host start iterateMoleculesParallel");
391 144 : Kokkos::parallel_for(
392 432 : "simplemd::MoleculeContainer::iterateMoleculesParallel", Kokkos::RangePolicy<MainExecSpace>(0, _linkedCellIsGhostCell.size()),
393 724992 : KOKKOS_CLASS_LAMBDA(const unsigned int i) {
394 1400544 : printNonGhostCells(i == 0, "device start iterateMoleculesParallel");
395 1417644 : for (unsigned int j = 0; j < _linkedCellNumMolecules(i); j++) {
396 34200 : a.handleMolecule(getMoleculeAt(i, j));
397 : }
398 1400544 : printNonGhostCells(i == 0, "device end iterateMoleculesParallel");
399 : });
400 144 : Kokkos::fence(); // Ensure results are available on the host
401 144 : printNonGhostCells(true, "host end iterateMoleculesParallel");
402 : a.endMoleculeIteration();
403 : #endif
404 144 : }
405 :
406 : template <class A> void simplemd::MoleculeContainer::iterateMoleculesWithCell(A& a) {
407 : if constexpr (A::IsParallel) {
408 : iterateMoleculesWithCellParallel(a);
409 : } else {
410 : iterateMoleculesWithCellSerial(a);
411 : }
412 : }
413 :
414 : template <class A> void simplemd::MoleculeContainer::handleCellNeighbors(A& a, Molecule& m, const LinkedCell& cell) const {
415 : unsigned int index = cell.getIndex();
416 : for (unsigned int i = 0; i < 26; i++) {
417 : auto cell2 = (*this)[index + _neighborOffsets(i)];
418 : a.handleMolecule(m, cell2);
419 : }
420 : }
421 :
422 : template <class A> void simplemd::MoleculeContainer::iterateMoleculesWithCellSerial(A& a) {
423 : a.beginMoleculeIteration();
424 : for (unsigned int i = 0; i < _linkedCellNumMolecules.size(); i++) {
425 : for (unsigned int j = 0; j < _linkedCellNumMolecules(i); j++) {
426 : simplemd::LinkedCell cell = (*this)[i];
427 : Molecule& m = getMoleculeAt(i, j);
428 : a.handleMolecule(m, cell);
429 : handleCellNeighbors(a, m, cell);
430 : }
431 : }
432 : a.endMoleculeIteration();
433 : }
434 :
435 : template <class A> void simplemd::MoleculeContainer::iterateMoleculesWithCellParallel(A& a) {
436 : a.beginMoleculeIteration();
437 : const unsigned int threads_per_cell = 5;
438 : const unsigned int length = _linkedCellNumMolecules.size() * threads_per_cell;
439 : Kokkos::parallel_for(
440 : "simplemd::MoleculeContainer::iterateMoleculesWithCellParallel", Kokkos::RangePolicy<MainExecSpace>(0, length),
441 : KOKKOS_CLASS_LAMBDA(const unsigned int i) {
442 : const unsigned int cellIndex = i / threads_per_cell;
443 : simplemd::LinkedCell cell = (*this)[cellIndex];
444 :
445 : for (unsigned int j = i % threads_per_cell; j < _linkedCellNumMolecules(cellIndex); j += threads_per_cell) {
446 : Molecule& m = getMoleculeAt(cellIndex, j);
447 : a.handleMolecule(m, cell);
448 : handleCellNeighbors(a, m, cell);
449 : }
450 : });
451 : Kokkos::fence(); // Ensure results are available on the host
452 :
453 : a.endMoleculeIteration();
454 : }
455 :
456 : template <class A>
457 4 : void simplemd::MoleculeContainer::iterateCellsSerial(A& a, const tarch::la::Vector<MD_DIM, unsigned int>& lowerLeftFrontCell,
458 : const tarch::la::Vector<MD_DIM, unsigned int>& cellRange) {
459 4 : unsigned int index = 0;
460 : #if (MD_ERROR == MD_YES)
461 16 : for (unsigned int d = 0; d < MD_DIM; d++) {
462 12 : if (cellRange[d] == 0) {
463 0 : Kokkos::printf("cellRange(%d)==0!\n", d);
464 : Kokkos::abort("ERROR simplemd::MoleculeContainer::iterateCells: zero in cell range\n");
465 : }
466 12 : if (lowerLeftFrontCell[d] + cellRange[d] > 2 * _ghostCellLayerThickness[d] + _numLocalCellsNoGhost[d]) {
467 : Kokkos::abort("ERROR simplemd::MoleculeContainer::iterateCells(): "
468 : "defined Range does not fit into local sub-domain!\n");
469 : }
470 : }
471 : #endif
472 :
473 : // start iteration();
474 : a.beginCellIteration();
475 4 : tarch::la::Vector<MD_DIM, unsigned int> coords(0);
476 : // loop over domain
477 : #if (MD_DIM > 2)
478 204 : for (coords[2] = lowerLeftFrontCell[2]; coords[2] < lowerLeftFrontCell[2] + cellRange[2]; coords[2]++) {
479 : #endif
480 : #if (MD_DIM > 1)
481 6200 : for (coords[1] = lowerLeftFrontCell[1]; coords[1] < lowerLeftFrontCell[1] + cellRange[1]; coords[1]++) {
482 : #endif
483 306000 : for (coords[0] = lowerLeftFrontCell[0]; coords[0] < lowerLeftFrontCell[0] + cellRange[0]; coords[0]++) {
484 : #if (MD_DEBUG == MD_YES)
485 : std::cout << "Handle cell " << coords << std::endl;
486 : #endif
487 300000 : index = vectorIndexToLinear(coords);
488 300000 : simplemd::LinkedCell cell = (*this)[index];
489 300000 : a.handleCell(cell);
490 : }
491 : #if (MD_DIM > 1)
492 : }
493 : #endif
494 : #if (MD_DIM > 2)
495 : }
496 : #endif
497 :
498 : // end iteration();
499 : a.endCellIteration();
500 4 : }
501 :
502 : template <class A> void simplemd::MoleculeContainer::iterateCellsSerial(A& a) { iterateCellsSerial(a, _ghostCellLayerThickness, _numLocalCellsNoGhost); }
503 :
504 : template <class A>
505 4 : void simplemd::MoleculeContainer::iterateCellPairsSerial(A& a, const tarch::la::Vector<MD_DIM, unsigned int>& lowerLeftFrontCell,
506 : const tarch::la::Vector<MD_DIM, unsigned int>& cellRange) {
507 4 : tarch::la::Vector<MD_LINKED_CELL_NEIGHBOURS / 2, unsigned int> neighbourOffset;
508 4 : tarch::la::Vector<MD_LINKED_CELL_NEIGHBOURS / 2, unsigned int> indexOffset;
509 : #if (MD_DIM == 1)
510 : neighbourOffset[0] = 1;
511 : indexOffset[0] = 0;
512 : #elif (MD_DIM == 2)
513 : indexOffset[0] = 0;
514 : neighbourOffset[0] = 1;
515 : indexOffset[1] = 0;
516 : neighbourOffset[1] = _numLocalCellsNoGhost[0] + 2;
517 : indexOffset[2] = 0;
518 : neighbourOffset[2] = _numLocalCellsNoGhost[0] + 3;
519 : indexOffset[3] = 1;
520 : neighbourOffset[3] = _numLocalCellsNoGhost[0] + 2;
521 : #elif (MD_DIM == 3)
522 4 : indexOffset[0] = 0;
523 4 : neighbourOffset[0] = 1;
524 4 : indexOffset[1] = 0;
525 4 : neighbourOffset[1] = _numLocalCellsNoGhost[0] + 2;
526 4 : indexOffset[2] = 0;
527 4 : neighbourOffset[2] = _numLocalCellsNoGhost[0] + 3;
528 4 : indexOffset[3] = 0;
529 4 : neighbourOffset[3] = (_numLocalCellsNoGhost[0] + 2) * (_numLocalCellsNoGhost[1] + 2);
530 4 : indexOffset[4] = 0;
531 4 : neighbourOffset[4] = (_numLocalCellsNoGhost[0] + 2) * (_numLocalCellsNoGhost[1] + 2) + 1;
532 4 : indexOffset[5] = 0;
533 4 : neighbourOffset[5] = (_numLocalCellsNoGhost[0] + 2) * (_numLocalCellsNoGhost[1] + 2) + (_numLocalCellsNoGhost[0] + 2);
534 4 : indexOffset[6] = 0;
535 4 : neighbourOffset[6] = (_numLocalCellsNoGhost[0] + 2) * (_numLocalCellsNoGhost[1] + 2) + (_numLocalCellsNoGhost[0] + 2) + 1;
536 :
537 4 : indexOffset[7] = 1;
538 4 : neighbourOffset[7] = _numLocalCellsNoGhost[0] + 2;
539 4 : indexOffset[8] = 1;
540 4 : neighbourOffset[8] = (_numLocalCellsNoGhost[0] + 2) * (_numLocalCellsNoGhost[1] + 2);
541 4 : indexOffset[9] = 1;
542 4 : neighbourOffset[9] = (_numLocalCellsNoGhost[0] + 2) * (_numLocalCellsNoGhost[1] + 2) + (_numLocalCellsNoGhost[0] + 2);
543 :
544 4 : indexOffset[10] = _numLocalCellsNoGhost[0] + 2;
545 4 : neighbourOffset[10] = (_numLocalCellsNoGhost[0] + 2) * (_numLocalCellsNoGhost[1] + 2);
546 4 : indexOffset[11] = _numLocalCellsNoGhost[0] + 2;
547 4 : neighbourOffset[11] = (_numLocalCellsNoGhost[0] + 2) * (_numLocalCellsNoGhost[1] + 2) + 1;
548 :
549 4 : indexOffset[12] = (_numLocalCellsNoGhost[0] + 2) + 1;
550 4 : neighbourOffset[12] = (_numLocalCellsNoGhost[0] + 2) * (_numLocalCellsNoGhost[1] + 2);
551 : #endif
552 :
553 : #if (MD_ERROR == MD_YES)
554 16 : for (unsigned int d = 0; d < MD_DIM; d++) {
555 12 : if (cellRange[d] > 2 * _ghostCellLayerThickness[d] + _numLocalCellsNoGhost[d] - 1) {
556 0 : std::cout << "ERROR simplemd::MoleculeContainer::iterateCellPairs(): "
557 : "defined Range does not fit into local sub-domain!"
558 0 : << std::endl;
559 0 : exit(EXIT_FAILURE);
560 : }
561 : }
562 : #endif
563 :
564 : // start iteration();
565 : a.beginCellIteration();
566 4 : tarch::la::Vector<MD_DIM, unsigned int> coords(0);
567 4 : unsigned int coordsCell1Buffer(0);
568 4 : unsigned int coordsCell2Buffer(0);
569 : unsigned int index;
570 :
571 : // loop over domain
572 : #if (MD_DIM > 2)
573 208 : for (coords[2] = lowerLeftFrontCell[2]; coords[2] < lowerLeftFrontCell[2] + cellRange[2]; coords[2]++) {
574 : #endif
575 : #if (MD_DIM > 1)
576 6528 : for (coords[1] = lowerLeftFrontCell[1]; coords[1] < lowerLeftFrontCell[1] + cellRange[1]; coords[1]++) {
577 : #endif
578 328848 : for (coords[0] = lowerLeftFrontCell[0]; coords[0] < lowerLeftFrontCell[0] + cellRange[0]; coords[0]++) {
579 : // handle cell itself
580 322524 : index = vectorIndexToLinear(coords);
581 : #if (MD_DEBUG == MD_YES)
582 : std::cout << "iterateCellPairs: Single index " << index << std::endl;
583 : #endif
584 322524 : simplemd::LinkedCell cell = (*this)[index];
585 : a.handleCell(cell);
586 : // handle pairs (lower,left,back-oriented cells)
587 4515336 : for (unsigned int i = 0; i < MD_LINKED_CELL_NEIGHBOURS / 2; i++) {
588 : #if (MD_DEBUG == MD_YES)
589 : std::cout << "iterateCellPairs: Pair index " << index + indexOffset[i] << "," << index + neighbourOffset[i] << std::endl;
590 : #endif
591 :
592 4192812 : coordsCell1Buffer = index + indexOffset[i];
593 4192812 : coordsCell2Buffer = index + neighbourOffset[i];
594 4192812 : simplemd::LinkedCell cell1 = (*this)[coordsCell1Buffer];
595 4192812 : simplemd::LinkedCell cell2 = (*this)[coordsCell2Buffer];
596 4192812 : a.handleCellPair(cell1, cell2, coordsCell1Buffer, coordsCell2Buffer);
597 : }
598 : } // coords(0)
599 : #if (MD_DIM > 1)
600 : }
601 : #endif
602 : #if (MD_DIM > 2)
603 : }
604 : #endif
605 :
606 : // end iteration();
607 : a.endCellIteration();
608 4 : }
609 :
610 : template <class A> void simplemd::MoleculeContainer::iterateCellPairsSerial(A& a) {
611 : const tarch::la::Vector<MD_DIM, unsigned int> pairIterationStart(0);
612 : const tarch::la::Vector<MD_DIM, unsigned int> pairIterationLength(getLocalNumberOfCells() + getLocalIndexOfFirstCell());
613 : iterateCellPairsSerial(a, pairIterationStart, pairIterationLength);
614 : }
615 :
616 : template <class A>
617 4 : void simplemd::MoleculeContainer::iterateCellsParallel(A& a, const tarch::la::Vector<MD_DIM, unsigned int>& lowerLeftFrontCell,
618 : const tarch::la::Vector<MD_DIM, unsigned int>& cellRange) {
619 : #if (MD_ERROR == MD_YES)
620 16 : for (unsigned int d = 0; d < MD_DIM; d++) {
621 12 : if (cellRange[d] == 0) {
622 0 : Kokkos::printf("cellRange(%d)==0!\n", d);
623 : Kokkos::abort("ERROR simplemd::MoleculeContainer::iterateCells: zero in cell range\n");
624 : }
625 12 : if (lowerLeftFrontCell[d] + cellRange[d] > 2 * _ghostCellLayerThickness[d] + _numLocalCellsNoGhost[d]) {
626 : Kokkos::abort("ERROR simplemd::MoleculeContainer::iterateCells(): "
627 : "defined Range does not fit into local sub-domain!\n");
628 : }
629 : }
630 : #endif
631 :
632 : // start iteration();
633 : a.beginCellIteration();
634 : /**
635 : * The size<> Vector stores the number of cells, plus ghost layer.
636 : * If there are (1,1,1) ghost cells per dimension, getLocalIndexOfFirstCell will return (1,1,1)
637 : * Thus this is multiplied by 2 to account for ghost cells in both locations (begin, end) per axis
638 : * and then added to local number of cells
639 : */
640 4 : const tarch::la::Vector<MD_DIM, unsigned int> size(getLocalNumberOfCells() + 2u * getLocalIndexOfFirstCell());
641 4 : const int length = cellRange[0]
642 : #if (MD_DIM > 1)
643 4 : * cellRange[1]
644 : #endif
645 : #if (MD_DIM > 2)
646 4 : * cellRange[2]
647 : #endif
648 : ;
649 : // loop over domain, but with a single loop
650 4 : printNonGhostCells(true, "host start iterateCellsParallel");
651 4 : Kokkos::parallel_for(
652 24 : "simplemd::MoleculeContainer::iterateCellsParallel", Kokkos::RangePolicy<MainExecSpace>(0, length), KOKKOS_CLASS_LAMBDA(const unsigned int i) {
653 300000 : printNonGhostCells(i == 0, "device start iterateCellsParallel");
654 : // compute index of the current cell
655 : #if (MD_DIM > 1)
656 300000 : int helpIndex1 = i;
657 300000 : int helpIndex2 = 0;
658 : #endif
659 300000 : unsigned int index = 0;
660 :
661 : #if (MD_DIM > 2)
662 : // determine plane within traversed block
663 300000 : helpIndex2 = helpIndex1 / (cellRange[0] * cellRange[1]);
664 : // save rest of index in helpIndex1
665 300000 : helpIndex1 = helpIndex1 - helpIndex2 * (cellRange[0] * cellRange[1]);
666 : // compute contribution to index
667 300000 : index += (lowerLeftFrontCell[2] + helpIndex2) * size[0] * size[1];
668 : #endif
669 : #if (MD_DIM > 1)
670 : // determine plane within traversed block
671 300000 : helpIndex2 = helpIndex1 / cellRange[0];
672 : // save rest of index in helpIndex1
673 300000 : helpIndex1 = helpIndex1 - helpIndex2 * cellRange[0];
674 : // compute contribution to index
675 300000 : index += (lowerLeftFrontCell[1] + helpIndex2) * size[0];
676 : // compute contribution for last dimension
677 300000 : index += (lowerLeftFrontCell[0] + helpIndex1);
678 : #else
679 : index = lowerLeftFrontCell[0] + i;
680 : #endif
681 : #if (MD_DEBUG == MD_YES)
682 : Kokkos::printf("Handle cell %u\n", index);
683 : #endif
684 :
685 : // handle cell
686 300000 : auto cell = (*this)[index];
687 300000 : a.handleCell(cell);
688 300000 : printNonGhostCells(i == 0, "device end iterateCellsParallel");
689 : }); // Kokkos::parallel_for
690 4 : Kokkos::fence(); // Ensure results are available on the host
691 4 : printNonGhostCells(true, "host end iterateCellsParallel");
692 : // end iteration();
693 : a.endCellIteration();
694 4 : }
695 :
696 : template <class A> void simplemd::MoleculeContainer::iterateCellsParallel(A& a) { iterateCellsParallel(a, _ghostCellLayerThickness, _numLocalCellsNoGhost); }
697 :
698 : template <class A>
699 136 : void simplemd::MoleculeContainer::iterateCellPairsParallel(A& a, const tarch::la::Vector<MD_DIM, unsigned int>& lowerLeftFrontCell,
700 : const tarch::la::Vector<MD_DIM, unsigned int>& cellRange) {
701 136 : tarch::la::Vector<MD_LINKED_CELL_NEIGHBOURS / 2, unsigned int> neighbourOffset;
702 136 : tarch::la::Vector<MD_LINKED_CELL_NEIGHBOURS / 2, unsigned int> indexOffset;
703 : #if (MD_DIM == 1)
704 : neighbourOffset[0] = 1;
705 : indexOffset[0] = 0;
706 : #elif (MD_DIM == 2)
707 : indexOffset[0] = 0;
708 : neighbourOffset[0] = 1;
709 : indexOffset[1] = 0;
710 : neighbourOffset[1] = _numLocalCellsNoGhost[0] + 2;
711 : indexOffset[2] = 0;
712 : neighbourOffset[2] = _numLocalCellsNoGhost[0] + 3;
713 : indexOffset[3] = 1;
714 : neighbourOffset[3] = _numLocalCellsNoGhost[0] + 2;
715 : #elif (MD_DIM == 3)
716 136 : indexOffset[0] = 0;
717 136 : neighbourOffset[0] = 1;
718 136 : indexOffset[1] = 0;
719 136 : neighbourOffset[1] = _numLocalCellsNoGhost[0] + 2;
720 136 : indexOffset[2] = 0;
721 136 : neighbourOffset[2] = _numLocalCellsNoGhost[0] + 3;
722 136 : indexOffset[3] = 0;
723 136 : neighbourOffset[3] = (_numLocalCellsNoGhost[0] + 2) * (_numLocalCellsNoGhost[1] + 2);
724 136 : indexOffset[4] = 0;
725 136 : neighbourOffset[4] = (_numLocalCellsNoGhost[0] + 2) * (_numLocalCellsNoGhost[1] + 2) + 1;
726 136 : indexOffset[5] = 0;
727 136 : neighbourOffset[5] = (_numLocalCellsNoGhost[0] + 2) * (_numLocalCellsNoGhost[1] + 2) + (_numLocalCellsNoGhost[0] + 2);
728 136 : indexOffset[6] = 0;
729 136 : neighbourOffset[6] = (_numLocalCellsNoGhost[0] + 2) * (_numLocalCellsNoGhost[1] + 2) + (_numLocalCellsNoGhost[0] + 2) + 1;
730 :
731 136 : indexOffset[7] = 1;
732 136 : neighbourOffset[7] = _numLocalCellsNoGhost[0] + 2;
733 136 : indexOffset[8] = 1;
734 136 : neighbourOffset[8] = (_numLocalCellsNoGhost[0] + 2) * (_numLocalCellsNoGhost[1] + 2);
735 136 : indexOffset[9] = 1;
736 136 : neighbourOffset[9] = (_numLocalCellsNoGhost[0] + 2) * (_numLocalCellsNoGhost[1] + 2) + (_numLocalCellsNoGhost[0] + 2);
737 :
738 136 : indexOffset[10] = _numLocalCellsNoGhost[0] + 2;
739 136 : neighbourOffset[10] = (_numLocalCellsNoGhost[0] + 2) * (_numLocalCellsNoGhost[1] + 2);
740 136 : indexOffset[11] = _numLocalCellsNoGhost[0] + 2;
741 136 : neighbourOffset[11] = (_numLocalCellsNoGhost[0] + 2) * (_numLocalCellsNoGhost[1] + 2) + 1;
742 :
743 136 : indexOffset[12] = (_numLocalCellsNoGhost[0] + 2) + 1;
744 136 : neighbourOffset[12] = (_numLocalCellsNoGhost[0] + 2) * (_numLocalCellsNoGhost[1] + 2);
745 : #endif
746 :
747 : #if (MD_ERROR == MD_YES)
748 544 : for (unsigned int d = 0; d < MD_DIM; d++) {
749 408 : if (cellRange[d] > 2 * _ghostCellLayerThickness[d] + _numLocalCellsNoGhost[d] - 1) {
750 : Kokkos::abort("ERROR simplemd::MoleculeContainer::iterateCellPairs(): "
751 : "defined Range does not fit into local sub-domain!\n");
752 : }
753 : }
754 : #endif
755 :
756 : // start iteration();
757 132 : a.beginCellIteration();
758 : /**
759 : * The size<> Vector stores the number of cells, plus ghost layer.
760 : * If there are (1,1,1) ghost cells per dimension, getLocalIndexOfFirstCell will return (1,1,1)
761 : * Thus this is multiplied by 2 to account for ghost cells in both locations (begin, end) per axis
762 : * and then added to local number of cells
763 : */
764 136 : const tarch::la::Vector<MD_DIM, unsigned int> size(getLocalNumberOfCells() + 2u * getLocalIndexOfFirstCell());
765 :
766 : // iterate over the domain in a red-black manner
767 : #if (MD_DIM > 2)
768 408 : for (unsigned int z = 0; z < 2; z++) {
769 : #endif
770 : #if (MD_DIM > 1)
771 816 : for (unsigned int y = 0; y < 2; y++) {
772 : #endif
773 1632 : for (unsigned int x = 0; x < 2; x++) {
774 : // determine range/ length of blocks for red-black traversal.
775 : // For odd block sizes, we need to do some more work in the
776 : // x/y/z==0-traversals. The second x/y/z==1-traversals are reduced by
777 : // the normal integer-rounding in this case.
778 1088 : const tarch::la::Vector<MD_DIM, unsigned int> lengthVector((cellRange[0] + (cellRange[0] % 2) * (x == 0)) / 2
779 : #if (MD_DIM > 1)
780 : ,
781 1088 : (cellRange[1] + (cellRange[1] % 2) * (y == 0)) / 2
782 : #endif
783 : #if (MD_DIM > 2)
784 : ,
785 1088 : (cellRange[2] + (cellRange[2] % 2) * (z == 0)) / 2
786 : #endif
787 : );
788 1088 : const int length = lengthVector[0]
789 : #if (MD_DIM > 1)
790 1088 : * lengthVector[1]
791 : #endif
792 : #if (MD_DIM > 2)
793 1088 : * lengthVector[2]
794 : #endif
795 : ;
796 :
797 : // parallelise loop for all cells that are to be traversed in this way
798 1088 : printNonGhostCells(true, "host start iterateCellPairsParallel");
799 1088 : Kokkos::parallel_for(
800 587656 : "simplemd::MoleculeContainer::iterateCellPairsParallel", Kokkos::RangePolicy<MainExecSpace>(0, length), KOKKOS_CLASS_LAMBDA(const unsigned int j) {
801 612528 : printNonGhostCells(j == 0, "device start iterateCellPairsParallel");
802 : // compute index of the current cell
803 612528 : unsigned int index = 0;
804 : #if (MD_DIM > 1)
805 612528 : int helpIndex1 = j;
806 612528 : int helpIndex2 = 0;
807 : #endif
808 : unsigned int coordsCell1Buffer;
809 : unsigned int coordsCell2Buffer;
810 :
811 : #if (MD_DIM > 2)
812 : // determine plane within traversed block
813 612528 : helpIndex2 = helpIndex1 / (lengthVector[0] * lengthVector[1]);
814 : // save rest of index in helpIndex1
815 612528 : helpIndex1 = helpIndex1 - helpIndex2 * (lengthVector[0] * lengthVector[1]);
816 : // compute contribution to index
817 612528 : index += (lowerLeftFrontCell[2] + 2 * helpIndex2 + z) * size[0] * size[1];
818 : #endif
819 : #if (MD_DIM > 1)
820 : // determine plane within traversed block
821 612528 : helpIndex2 = helpIndex1 / lengthVector[0];
822 : // save rest of index in helpIndex1
823 612528 : helpIndex1 = helpIndex1 - helpIndex2 * lengthVector[0];
824 : // compute contribution to index
825 612528 : index += (lowerLeftFrontCell[1] + 2 * helpIndex2 + y) * size[0];
826 : // compute contribution for last dimension
827 612528 : index += (lowerLeftFrontCell[0] + 2 * helpIndex1 + x);
828 : #else
829 : index = lowerLeftFrontCell[0] + 2 * j + x;
830 : #endif
831 : #if (MD_DEBUG == MD_YES)
832 : Kokkos::printf("Handle cell %d\n", index);
833 : #endif
834 :
835 612528 : simplemd::LinkedCell cell = (*this)[index];
836 290004 : a.handleCell(cell);
837 : // handle pairs (lower,left,back-oriented cells)
838 8575392 : for (unsigned int i = 0; i < MD_LINKED_CELL_NEIGHBOURS / 2; i++) {
839 : #if (MD_DEBUG == MD_YES)
840 : Kokkos::printf("iterateCellPairs: Pair index %d, %d\n", index + indexOffset[i], index + neighbourOffset[i]);
841 : #endif
842 7962864 : coordsCell1Buffer = index + indexOffset[i];
843 7962864 : coordsCell2Buffer = index + neighbourOffset[i];
844 7962864 : auto cell1 = (*this)[coordsCell1Buffer];
845 7962864 : auto cell2 = (*this)[coordsCell2Buffer];
846 7962864 : a.handleCellPair(cell1, cell2, coordsCell1Buffer, coordsCell2Buffer);
847 : }
848 612528 : printNonGhostCells(j == 0, "device end iterateCellPairsParallel");
849 : }); // j, Kokkos::parallel_for
850 1088 : Kokkos::fence(); // Ensure results are available on the host
851 1088 : printNonGhostCells(true, "host end iterateCellPairsParallel");
852 : } // x
853 : #if (MD_DIM > 1)
854 : } // y
855 : #endif
856 : #if (MD_DIM > 2)
857 : } // z
858 : #endif
859 : // end iteration();
860 : a.endCellIteration();
861 136 : }
862 :
863 : template <class A> void simplemd::MoleculeContainer::iterateCellPairsParallel(A& a) {
864 : const tarch::la::Vector<MD_DIM, unsigned int> pairIterationStart(0);
865 : const tarch::la::Vector<MD_DIM, unsigned int> pairIterationLength(getLocalNumberOfCells() + getLocalIndexOfFirstCell());
866 : iterateCellPairsParallel(a, pairIterationStart, pairIterationLength);
867 : }
868 :
869 : template <class A>
870 140 : void simplemd::MoleculeContainer::iterateCellPairs(A& a, const tarch::la::Vector<MD_DIM, unsigned int>& lowerLeftFrontCell,
871 : const tarch::la::Vector<MD_DIM, unsigned int>& cellRange) {
872 : if constexpr (A::IsParallel) {
873 136 : iterateCellPairsParallel(a, lowerLeftFrontCell, cellRange);
874 : } else {
875 4 : iterateCellPairsSerial(a, lowerLeftFrontCell, cellRange);
876 : }
877 140 : }
878 :
879 140 : template <class A> void simplemd::MoleculeContainer::iterateCellPairs(A& a) {
880 140 : const tarch::la::Vector<MD_DIM, unsigned int> pairIterationStart(0);
881 140 : const tarch::la::Vector<MD_DIM, unsigned int> pairIterationLength(getLocalNumberOfCells() + getLocalIndexOfFirstCell());
882 140 : iterateCellPairs(a, pairIterationStart, pairIterationLength);
883 140 : }
884 :
885 : template <class A>
886 8 : void simplemd::MoleculeContainer::iterateCells(A& a, const tarch::la::Vector<MD_DIM, unsigned int>& lowerLeftFrontCell,
887 : const tarch::la::Vector<MD_DIM, unsigned int>& cellRange) {
888 : if constexpr (A::IsParallel) {
889 4 : iterateCellsParallel(a, lowerLeftFrontCell, cellRange);
890 : } else {
891 4 : iterateCellsSerial(a, lowerLeftFrontCell, cellRange);
892 : }
893 8 : }
894 :
895 8 : template <class A> void simplemd::MoleculeContainer::iterateCells(A& a) { iterateCells(a, _ghostCellLayerThickness, _numLocalCellsNoGhost); }
896 :
897 : #endif // _MOLECULARDYNAMICS_MOLECULARCONTAINER_H_
|
