From 15bb09e17f07898f432355a065c59fee26bf6152 Mon Sep 17 00:00:00 2001
From: Sebastian Eibl <sebastian.eibl@fau.de>
Date: Mon, 2 Sep 2019 15:16:31 +0200
Subject: [PATCH] introduced SparseLinkedCells
---
apps/benchmarks/GranularGas/GranularGas.cfg | 7 +-
python/mesa_pd.py | 3 +
python/mesa_pd/data/SparseLinkedCells.py | 8 +
python/mesa_pd/data/__init__.py | 4 +-
.../InsertParticleIntoSparseLinkedCells.py | 19 +
python/mesa_pd/kernel/__init__.py | 2 +
.../templates/data/LinkedCells.templ.h | 47 +-
.../templates/data/SparseLinkedCells.templ.h | 302 +++++++++++++
.../InsertParticleIntoLinkedCells.templ.h | 4 +-
...nsertParticleIntoSparseLinkedCells.templ.h | 98 +++++
src/mesa_pd/data/LinkedCells.h | 55 ++-
src/mesa_pd/data/SparseLinkedCells.h | 400 ++++++++++++++++++
.../kernel/InsertParticleIntoLinkedCells.h | 4 +-
.../InsertParticleIntoSparseLinkedCells.h | 95 +++++
tests/mesa_pd/CMakeLists.txt | 6 +
tests/mesa_pd/data/LinkedCells.cpp | 65 +++
tests/mesa_pd/data/SparseLinkedCells.cpp | 65 +++
tests/mesa_pd/kernel/GenerateLinkedCells.cpp | 4 +-
18 files changed, 1155 insertions(+), 33 deletions(-)
create mode 100644 python/mesa_pd/data/SparseLinkedCells.py
create mode 100644 python/mesa_pd/kernel/InsertParticleIntoSparseLinkedCells.py
create mode 100644 python/mesa_pd/templates/data/SparseLinkedCells.templ.h
create mode 100644 python/mesa_pd/templates/kernel/InsertParticleIntoSparseLinkedCells.templ.h
create mode 100644 src/mesa_pd/data/SparseLinkedCells.h
create mode 100644 src/mesa_pd/kernel/InsertParticleIntoSparseLinkedCells.h
create mode 100644 tests/mesa_pd/data/LinkedCells.cpp
create mode 100644 tests/mesa_pd/data/SparseLinkedCells.cpp
diff --git a/apps/benchmarks/GranularGas/GranularGas.cfg b/apps/benchmarks/GranularGas/GranularGas.cfg
index 0733512f4..2c0c47552 100644
--- a/apps/benchmarks/GranularGas/GranularGas.cfg
+++ b/apps/benchmarks/GranularGas/GranularGas.cfg
@@ -2,11 +2,14 @@ GranularGas
{
simulationCorner < 0, 0, 0 >;
simulationDomain < 80, 80, 80 >;
- blocks < 2,2,2 >;
+ blocks < 1,1,1 >;
isPeriodic < 0, 0, 0 >;
- initialRefinementLevel 1;
+ initialRefinementLevel 3;
sorting linear;
+ LBAlgorithm Morton;
+ baseWeight 1;
+
normal <1,1,1>;
radius 0.6;
spacing 1.0;
diff --git a/python/mesa_pd.py b/python/mesa_pd.py
index c8f000e51..b010eb223 100755
--- a/python/mesa_pd.py
+++ b/python/mesa_pd.py
@@ -32,6 +32,7 @@ if __name__ == '__main__':
ps = data.ParticleStorage()
ch = data.ContactHistory()
lc = data.LinkedCells()
+ slc = data.SparseLinkedCells()
ss = data.ShapeStorage(ps, shapes)
cs = data.ContactStorage()
@@ -95,6 +96,7 @@ if __name__ == '__main__':
kernels.append( kernel.ForceLJ() )
kernels.append( kernel.HeatConduction() )
kernels.append( kernel.InsertParticleIntoLinkedCells() )
+ kernels.append( kernel.InsertParticleIntoSparseLinkedCells() )
kernels.append( kernel.LinearSpringDashpot() )
kernels.append( kernel.NonLinearSpringDashpot() )
kernels.append( kernel.SingleCast(shapes) )
@@ -121,6 +123,7 @@ if __name__ == '__main__':
ps.generate(args.path + "/src/mesa_pd/")
ch.generate(args.path + "/src/mesa_pd/")
lc.generate(args.path + "/src/mesa_pd/")
+ slc.generate(args.path + "/src/mesa_pd/")
ss.generate(args.path + "/src/mesa_pd/")
cs.generate(args.path + "/src/mesa_pd/")
diff --git a/python/mesa_pd/data/SparseLinkedCells.py b/python/mesa_pd/data/SparseLinkedCells.py
new file mode 100644
index 000000000..516acfd82
--- /dev/null
+++ b/python/mesa_pd/data/SparseLinkedCells.py
@@ -0,0 +1,8 @@
+# -*- coding: utf-8 -*-
+
+import numpy as np
+from ..utility import generateFile
+
+class SparseLinkedCells:
+ def generate(self, path):
+ generateFile(path, 'data/SparseLinkedCells.templ.h')
diff --git a/python/mesa_pd/data/__init__.py b/python/mesa_pd/data/__init__.py
index a8b066957..d09dc6e57 100644
--- a/python/mesa_pd/data/__init__.py
+++ b/python/mesa_pd/data/__init__.py
@@ -5,9 +5,11 @@ from .ContactStorage import ContactStorage
from .LinkedCells import LinkedCells
from .ParticleStorage import ParticleStorage
from .ShapeStorage import ShapeStorage
+from .SparseLinkedCells import SparseLinkedCells
__all__ = ['ContactHistory',
'ContactStorage',
'GeometryStorage',
'LinkedCells',
- 'ParticleStorage']
+ 'ParticleStorage',
+ 'SparseLinkedCells']
diff --git a/python/mesa_pd/kernel/InsertParticleIntoSparseLinkedCells.py b/python/mesa_pd/kernel/InsertParticleIntoSparseLinkedCells.py
new file mode 100644
index 000000000..9bb841b01
--- /dev/null
+++ b/python/mesa_pd/kernel/InsertParticleIntoSparseLinkedCells.py
@@ -0,0 +1,19 @@
+# -*- coding: utf-8 -*-
+
+from mesa_pd.accessor import Accessor
+from mesa_pd.utility import generateFile
+
+class InsertParticleIntoSparseLinkedCells:
+ def __init__(self):
+ self.accessor = Accessor()
+ self.accessor.require("position", "walberla::mesa_pd::Vec3", access="g")
+ self.accessor.require("flags", "walberla::mesa_pd::data::particle_flags::FlagT", access="g")
+ self.accessor.require("nextParticle", "size_t", access="gs" )
+
+ def getRequirements(self):
+ return self.accessor
+
+ def generate(self, path):
+ context = dict()
+ context["interface"] = self.accessor.properties
+ generateFile(path, 'kernel/InsertParticleIntoSparseLinkedCells.templ.h', context)
diff --git a/python/mesa_pd/kernel/__init__.py b/python/mesa_pd/kernel/__init__.py
index 9ce98cad5..758512e84 100644
--- a/python/mesa_pd/kernel/__init__.py
+++ b/python/mesa_pd/kernel/__init__.py
@@ -6,6 +6,7 @@ from .ExplicitEulerWithShape import ExplicitEulerWithShape
from .ForceLJ import ForceLJ
from .HeatConduction import HeatConduction
from .InsertParticleIntoLinkedCells import InsertParticleIntoLinkedCells
+from .InsertParticleIntoSparseLinkedCells import InsertParticleIntoSparseLinkedCells
from .LinearSpringDashpot import LinearSpringDashpot
from .NonLinearSpringDashpot import NonLinearSpringDashpot
from .SingleCast import SingleCast
@@ -21,6 +22,7 @@ __all__ = ['DoubleCast',
'ForceLJ',
'HeatConduction',
'InsertParticleIntoLinkedCells',
+ 'InsertParticleIntoSparseLinkedCells',
'LinearSpringDashpot',
'NonLinearSpringDashpot',
'SingleCast',
diff --git a/python/mesa_pd/templates/data/LinkedCells.templ.h b/python/mesa_pd/templates/data/LinkedCells.templ.h
index c5107a6ef..3cd762e5e 100644
--- a/python/mesa_pd/templates/data/LinkedCells.templ.h
+++ b/python/mesa_pd/templates/data/LinkedCells.templ.h
@@ -93,7 +93,10 @@ struct LinkedCells
};
inline
-math::AABB getCellAABB(const LinkedCells& ll, const int hash0, const int hash1, const int hash2)
+math::AABB getCellAABB(const LinkedCells& ll,
+ const int64_t hash0,
+ const int64_t hash1,
+ const int64_t hash2)
{
{%- for dim in range(3) %}
WALBERLA_ASSERT_GREATER_EQUAL(hash{{dim}}, 0);
@@ -110,13 +113,35 @@ math::AABB getCellAABB(const LinkedCells& ll, const int hash0, const int hash1,
}
inline
-int getCellIdx(const LinkedCells& ll, const int hash0, const int hash1, const int hash2)
+uint_t getCellIdx(const LinkedCells& ll,
+ const int64_t hash0,
+ const int64_t hash1,
+ const int64_t hash2)
{
{%- for dim in range(3) %}
WALBERLA_ASSERT_GREATER_EQUAL(hash{{dim}}, 0);
WALBERLA_ASSERT_LESS(hash{{dim}}, ll.numCellsPerDim_[{{dim}}]);
{%- endfor %}
- return hash2 * ll.numCellsPerDim_[1] * ll.numCellsPerDim_[0] + hash1 * ll.numCellsPerDim_[0] + hash0;
+ return uint_c(hash2 * ll.numCellsPerDim_[1] * ll.numCellsPerDim_[0] + hash1 * ll.numCellsPerDim_[0] + hash0);
+}
+
+inline
+void getCellCoordinates(const LinkedCells& ll,
+ const uint64_t idx,
+ int64_t& hash0,
+ int64_t& hash1,
+ int64_t& hash2)
+{
+ hash2 = int64_c(idx) / (ll.numCellsPerDim_[1] * ll.numCellsPerDim_[0]);
+ hash1 = (int64_c(idx) - (hash2 * ll.numCellsPerDim_[1] * ll.numCellsPerDim_[0])) / (ll.numCellsPerDim_[0]);
+ hash0 = int64_c(idx) - hash2 * ll.numCellsPerDim_[1] * ll.numCellsPerDim_[0] - hash1 * ll.numCellsPerDim_[0];
+
+ WALBERLA_ASSERT_GREATER_EQUAL(hash0, 0);
+ WALBERLA_ASSERT_LESS(hash0, ll.numCellsPerDim_[0]);
+ WALBERLA_ASSERT_GREATER_EQUAL(hash1, 0);
+ WALBERLA_ASSERT_LESS(hash1, ll.numCellsPerDim_[1]);
+ WALBERLA_ASSERT_GREATER_EQUAL(hash2, 0);
+ WALBERLA_ASSERT_LESS(hash2, ll.numCellsPerDim_[2]);
}
inline
@@ -143,15 +168,17 @@ LinkedCells::LinkedCells(const math::AABB& domain, const Vec3& cellDiameter)
WALBERLA_CHECK_GREATER_EQUAL(numCellsPerDim_[{{dim}}], 0);
{%- endfor %}
+
+ std::fill(cells_.begin(), cells_.end(), -1);
}
void LinkedCells::clear()
{
const uint64_t cellsSize = cells_.size();
//clear existing linked cells
- #ifdef _OPENMP
- #pragma omp parallel for schedule(static)
- #endif
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static)
+#endif
for (int64_t i = 0; i < int64_c(cellsSize); ++i)
cells_[uint64_c(i)] = -1;
infiniteParticles_ = -1;
@@ -172,8 +199,8 @@ inline void LinkedCells::forEachParticlePair{%- if half %}Half{%- endif %}(const
{
for (int x = 0; x < numCellsPerDim_[0]; ++x)
{
- const int cell_idx = getCellIdx(*this, x, y, z); ///< current cell index
- int p_idx = cells_[uint_c(cell_idx)]; ///< current particle index
+ const uint_t cell_idx = getCellIdx(*this, x, y, z); ///< current cell index
+ int p_idx = cells_[cell_idx]; ///< current particle index
int np_idx = -1; ///< particle to be checked against
while (p_idx != -1)
@@ -212,11 +239,11 @@ inline void LinkedCells::forEachParticlePair{%- if half %}Half{%- endif %}(const
if (ny >= numCellsPerDim_[1]) continue;
if (nz >= numCellsPerDim_[2]) continue;
- const int ncell_idx = getCellIdx(*this, nx, ny, nz); ///< neighbor cell index
+ const uint_t ncell_idx = getCellIdx(*this, nx, ny, nz); ///< neighbor cell index
WALBERLA_ASSERT_GREATER_EQUAL(p_idx, 0);
WALBERLA_ASSERT_LESS(p_idx, acForLC.size());
- np_idx = cells_[uint_c(ncell_idx)]; ///< neighbor particle index
+ np_idx = cells_[ncell_idx]; ///< neighbor particle index
while (np_idx != -1)
{
WALBERLA_ASSERT_GREATER_EQUAL(np_idx, 0);
diff --git a/python/mesa_pd/templates/data/SparseLinkedCells.templ.h b/python/mesa_pd/templates/data/SparseLinkedCells.templ.h
new file mode 100644
index 000000000..139dbb0a4
--- /dev/null
+++ b/python/mesa_pd/templates/data/SparseLinkedCells.templ.h
@@ -0,0 +1,302 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file SparseLinkedCells.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/IAccessor.h>
+#include <mesa_pd/data/ParticleStorage.h>
+
+#include <core/Abort.h>
+#include <core/debug/Debug.h>
+#include <core/math/AABB.h>
+#include <stencil/D3Q27.h>
+
+#include <atomic>
+#include <cmath>
+#include <vector>
+
+namespace walberla {
+namespace mesa_pd {
+namespace data {
+
+struct SparseLinkedCells
+{
+ SparseLinkedCells(const math::AABB& domain, const real_t cellDiameter)
+ : SparseLinkedCells(domain, Vec3(cellDiameter,cellDiameter,cellDiameter))
+ {}
+ SparseLinkedCells(const math::AABB& domain, const Vec3& cellDiameter);
+
+ void clear();
+
+ /**
+ * Calls the provided functor \p func for all particle pairs.
+ *
+ * Additional arguments can be provided. No pairs with twice the same particle.
+ * Call syntax for the provided functor
+ * \code
+ * func( *this, i, j, std::forward<Args>(args)... );
+ * \endcode
+ * \param openmp enables/disables OpenMP parallelization of the kernel call
+ */
+ template <typename Selector, typename Accessor, typename Func, typename... Args>
+ void forEachParticlePair(const bool openmp,
+ const Selector& selector,
+ Accessor& acForLC,
+ Func&& func,
+ Args&&... args) const;
+ /**
+ * Calls the provided functor \p func for all particle pairs.
+ *
+ * Additional arguments can be provided. No pairs with twice the same particle are generated.
+ * No pair is called twice!
+ * Call syntax for the provided functor
+ * \code
+ * func( *this, i, j, std::forward<Args>(args)... );
+ * \endcode
+ * \param openmp enables/disables OpenMP parallelization of the kernel call
+ */
+ template <typename Selector, typename Accessor, typename Func, typename... Args>
+ void forEachParticlePairHalf(const bool openmp,
+ const Selector& selector,
+ Accessor& acForLC,
+ Func&& func,
+ Args&&... args) const;
+
+ math::AABB domain_ {}; ///< local domain covered by this data structure
+ Vector3<int> numCellsPerDim_ {}; ///< number of linked cells per dimension
+ Vec3 cellDiameter_ {};
+ Vec3 invCellDiameter_ {};
+ std::atomic<int> infiniteParticles_ {}; ///< data structure for particles to large for the cells
+ std::vector< std::atomic<int> > cells_ {}; ///< actual cell data structure
+ std::vector<size_t> nonEmptyCells_ {}; ///< list of cells containing particles
+};
+
+inline
+math::AABB getCellAABB(const SparseLinkedCells& ll,
+ const int64_t hash0,
+ const int64_t hash1,
+ const int64_t hash2)
+{
+ {%- for dim in range(3) %}
+ WALBERLA_ASSERT_GREATER_EQUAL(hash{{dim}}, 0);
+ WALBERLA_ASSERT_LESS(hash{{dim}}, ll.numCellsPerDim_[{{dim}}]);
+ {%- endfor %}
+ const auto& minCorner = ll.domain_.minCorner();
+ const real_t xMin = ll.cellDiameter_[0] * real_c(hash0) + minCorner[0];
+ const real_t yMin = ll.cellDiameter_[1] * real_c(hash1) + minCorner[1];
+ const real_t zMin = ll.cellDiameter_[2] * real_c(hash2) + minCorner[2];
+ const real_t xMax = ll.cellDiameter_[0] * real_c(hash0 + 1) + minCorner[0];
+ const real_t yMax = ll.cellDiameter_[1] * real_c(hash1 + 1) + minCorner[1];
+ const real_t zMax = ll.cellDiameter_[2] * real_c(hash2 + 1) + minCorner[2];
+ return math::AABB(xMin, yMin, zMin, xMax, yMax, zMax);
+}
+
+inline
+uint_t getCellIdx(const SparseLinkedCells& ll,
+ const int64_t hash0,
+ const int64_t hash1,
+ const int64_t hash2)
+{
+ {%- for dim in range(3) %}
+ WALBERLA_ASSERT_GREATER_EQUAL(hash{{dim}}, 0);
+ WALBERLA_ASSERT_LESS(hash{{dim}}, ll.numCellsPerDim_[{{dim}}]);
+ {%- endfor %}
+ return uint_c(hash2 * ll.numCellsPerDim_[1] * ll.numCellsPerDim_[0] + hash1 * ll.numCellsPerDim_[0] + hash0);
+}
+
+inline
+void getCellCoordinates(const SparseLinkedCells& ll,
+ const uint64_t idx,
+ int64_t& hash0,
+ int64_t& hash1,
+ int64_t& hash2)
+{
+ hash2 = int64_c(idx) / (ll.numCellsPerDim_[1] * ll.numCellsPerDim_[0]);
+ hash1 = (int64_c(idx) - (hash2 * ll.numCellsPerDim_[1] * ll.numCellsPerDim_[0])) / (ll.numCellsPerDim_[0]);
+ hash0 = int64_c(idx) - hash2 * ll.numCellsPerDim_[1] * ll.numCellsPerDim_[0] - hash1 * ll.numCellsPerDim_[0];
+
+ WALBERLA_ASSERT_GREATER_EQUAL(hash0, 0);
+ WALBERLA_ASSERT_LESS(hash0, ll.numCellsPerDim_[0]);
+ WALBERLA_ASSERT_GREATER_EQUAL(hash1, 0);
+ WALBERLA_ASSERT_LESS(hash1, ll.numCellsPerDim_[1]);
+ WALBERLA_ASSERT_GREATER_EQUAL(hash2, 0);
+ WALBERLA_ASSERT_LESS(hash2, ll.numCellsPerDim_[2]);
+}
+
+inline
+SparseLinkedCells::SparseLinkedCells(const math::AABB& domain, const Vec3& cellDiameter)
+ : domain_(domain)
+ , numCellsPerDim_( static_cast<int>(std::ceil( domain.sizes()[0] / cellDiameter[0])),
+ static_cast<int>(std::ceil( domain.sizes()[1] / cellDiameter[1])),
+ static_cast<int>(std::ceil( domain.sizes()[2] / cellDiameter[2])) )
+ , cellDiameter_( domain.sizes()[0] / real_c(numCellsPerDim_[0]),
+ domain.sizes()[1] / real_c(numCellsPerDim_[1]),
+ domain.sizes()[2] / real_c(numCellsPerDim_[2]) )
+ , invCellDiameter_( real_t(1) / cellDiameter_[0], real_t(1) / cellDiameter_[1], real_t(1) / cellDiameter_[2] )
+ , cells_(uint_c(numCellsPerDim_[0]*numCellsPerDim_[1]*numCellsPerDim_[2]))
+ , nonEmptyCells_(uint_c(numCellsPerDim_[0]*numCellsPerDim_[1]*numCellsPerDim_[2]))
+{
+ //precondition
+ {%- for dim in range(3) %}
+ WALBERLA_CHECK_GREATER_EQUAL(cellDiameter[{{dim}}], real_t(0));
+ {%- endfor %}
+
+ //postcondition
+ {%- for dim in range(3) %}
+ WALBERLA_CHECK_GREATER_EQUAL(cellDiameter_[{{dim}}], real_t(0));
+ WALBERLA_CHECK_LESS_EQUAL(cellDiameter_[{{dim}}], cellDiameter[{{dim}}]);
+
+ WALBERLA_CHECK_GREATER_EQUAL(numCellsPerDim_[{{dim}}], 0);
+ {%- endfor %}
+
+ nonEmptyCells_.clear();
+ std::fill(cells_.begin(), cells_.end(), -1);
+}
+
+void SparseLinkedCells::clear()
+{
+ for (const auto v : nonEmptyCells_)
+ {
+ WALBERLA_ASSERT_LESS(v, cells_.size());
+ cells_[v] = -1;
+ }
+ nonEmptyCells_.clear();
+ infiniteParticles_ = -1;
+}
+
+{%- for half in [False, True] %}
+template <typename Selector, typename Accessor, typename Func, typename... Args>
+inline void SparseLinkedCells::forEachParticlePair{%- if half %}Half{%- endif %}(const bool openmp, const Selector& selector, Accessor& acForLC, Func&& func, Args&&... args) const
+{
+ static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
+ WALBERLA_UNUSED(openmp);
+
+ for (const auto cellIdx : nonEmptyCells_)
+ {
+ int64_t x = 0;
+ int64_t y = 0;
+ int64_t z = 0;
+ getCellCoordinates(*this, cellIdx, x, y, z);
+ int p_idx = cells_[cellIdx]; ///< current particle index
+ int np_idx = -1; ///< particle to be checked against
+
+ while (p_idx != -1)
+ {
+ WALBERLA_ASSERT_GREATER_EQUAL(p_idx, 0);
+ WALBERLA_ASSERT_LESS(p_idx, acForLC.size());
+
+ // check particles in own cell
+ np_idx = acForLC.getNextParticle(uint_c(p_idx)); ///< neighbor particle index
+ while (np_idx != -1)
+ {
+ WALBERLA_ASSERT_GREATER_EQUAL(np_idx, 0);
+ WALBERLA_ASSERT_LESS(np_idx, acForLC.size());
+
+ if (selector(uint_c(p_idx), uint_c(np_idx), acForLC))
+ {
+ func(uint_c(p_idx), uint_c(np_idx), std::forward<Args>(args)...);
+ {%- if not half %}
+ func(uint_c(np_idx), uint_c(p_idx), std::forward<Args>(args)...);
+ {%- endif %}
+ }
+
+ // go to next particle
+ np_idx = acForLC.getNextParticle(uint_c(np_idx));
+ }
+
+ // check particles in infiniteParticles list
+ np_idx = infiniteParticles_; ///< neighbor particle index
+ while (np_idx != -1)
+ {
+ WALBERLA_ASSERT_GREATER_EQUAL(np_idx, 0);
+ WALBERLA_ASSERT_LESS(np_idx, acForLC.size());
+
+ if (selector(uint_c(p_idx), uint_c(np_idx), acForLC))
+ {
+ func(uint_c(p_idx), uint_c(np_idx), std::forward<Args>(args)...);
+ {%- if not half %}
+ func(uint_c(np_idx), uint_c(p_idx), std::forward<Args>(args)...);
+ {%- endif %}
+ }
+
+ // go to next particle
+ np_idx = acForLC.getNextParticle(uint_c(np_idx));
+ }
+
+ // go to next particle
+ p_idx = acForLC.getNextParticle(uint_c(p_idx));
+ }
+
+ // check particles in neighboring cells (only positive ones)
+ for (auto dir : stencil::D3Q27::dir_pos)
+ {
+ const int64_t nx = x + int64_c(stencil::cx[dir]);
+ const int64_t ny = y + int64_c(stencil::cy[dir]);
+ const int64_t nz = z + int64_c(stencil::cz[dir]);
+ if (nx < 0) continue;
+ if (ny < 0) continue;
+ if (nz < 0) continue;
+ if (nx >= numCellsPerDim_[0]) continue;
+ if (ny >= numCellsPerDim_[1]) continue;
+ if (nz >= numCellsPerDim_[2]) continue;
+
+ const uint64_t ncell_idx = getCellIdx(*this, nx, ny, nz); ///< neighbor cell index
+
+ p_idx = cells_[cellIdx]; ///< current particle index
+ WALBERLA_ASSERT_GREATER_EQUAL(p_idx, 0);
+ WALBERLA_ASSERT_LESS(p_idx, acForLC.size());
+ while (p_idx != -1)
+ {
+ np_idx = cells_[ncell_idx]; ///< neighbor particle index
+ while (np_idx != -1)
+ {
+ WALBERLA_ASSERT_GREATER_EQUAL(np_idx, 0);
+ WALBERLA_ASSERT_LESS(np_idx, acForLC.size());
+
+ if (selector(uint_c(p_idx), uint_c(np_idx), acForLC))
+ {
+ func(uint_c(p_idx), uint_c(np_idx), std::forward<Args>(args)...);
+ {%- if not half %}
+ func(uint_c(np_idx), uint_c(p_idx), std::forward<Args>(args)...);
+ {%- endif %}
+ }
+
+ // go to next particle
+ np_idx = acForLC.getNextParticle(uint_c(np_idx));
+ }
+
+ // go to next particle
+ p_idx = acForLC.getNextParticle(uint_c(p_idx));
+ }
+ }
+ }
+}
+{%- endfor %}
+
+} //namespace data
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/python/mesa_pd/templates/kernel/InsertParticleIntoLinkedCells.templ.h b/python/mesa_pd/templates/kernel/InsertParticleIntoLinkedCells.templ.h
index ac2bcc557..d80b9cfd8 100644
--- a/python/mesa_pd/templates/kernel/InsertParticleIntoLinkedCells.templ.h
+++ b/python/mesa_pd/templates/kernel/InsertParticleIntoLinkedCells.templ.h
@@ -84,8 +84,8 @@ inline void InsertParticleIntoLinkedCells::operator()(const size_t p_idx, Access
if (hash{{dim}} < 0) hash{{dim}} = 0;
if (hash{{dim}} >= lc.numCellsPerDim_[{{dim}}]) hash{{dim}} = lc.numCellsPerDim_[{{dim}}] - 1;
{%- endfor %}
- int cell_idx = getCellIdx(lc, hash0, hash1, hash2);
- ac.setNextParticle(p_idx, lc.cells_[uint_c(cell_idx)].exchange(int_c(p_idx)));
+ uint_t cell_idx = getCellIdx(lc, hash0, hash1, hash2);
+ ac.setNextParticle(p_idx, lc.cells_[cell_idx].exchange(int_c(p_idx)));
}
}
diff --git a/python/mesa_pd/templates/kernel/InsertParticleIntoSparseLinkedCells.templ.h b/python/mesa_pd/templates/kernel/InsertParticleIntoSparseLinkedCells.templ.h
new file mode 100644
index 000000000..ca9f6ba01
--- /dev/null
+++ b/python/mesa_pd/templates/kernel/InsertParticleIntoSparseLinkedCells.templ.h
@@ -0,0 +1,98 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file InsertParticleIntoSparseLinkedCells.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/IAccessor.h>
+#include <mesa_pd/data/SparseLinkedCells.h>
+
+#include <vector>
+
+namespace walberla {
+namespace mesa_pd {
+namespace kernel {
+
+/**
+ * Inserts a particle into the data::SparseLinkedCells data structure
+ *
+ * \attention Make sure to data::SparseLinkedCells::clear() the data structure before
+ * reinserting new particles.
+ *
+ * This kernel requires the following particle accessor interface
+ * \code
+ {%- for prop in interface %}
+ {%- if 'g' in prop.access %}
+ * const {{prop.type}}& get{{prop.name | capFirst}}(const size_t p_idx) const;
+ {%- endif %}
+ {%- if 's' in prop.access %}
+ * void set{{prop.name | capFirst}}(const size_t p_idx, const {{prop.type}}& v);
+ {%- endif %}
+ {%- if 'r' in prop.access %}
+ * {{prop.type}}& get{{prop.name | capFirst}}Ref(const size_t p_idx);
+ {%- endif %}
+ *
+ {%- endfor %}
+ * \endcode
+ * \ingroup mesa_pd_kernel
+ */
+class InsertParticleIntoSparseLinkedCells
+{
+public:
+ template <typename Accessor>
+ void operator()(const size_t p_idx, Accessor& ac, data::SparseLinkedCells& lc) const;
+};
+
+template <typename Accessor>
+inline void InsertParticleIntoSparseLinkedCells::operator()(const size_t p_idx, Accessor& ac, data::SparseLinkedCells& lc) const
+{
+ static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
+
+ const auto& minCorner = lc.domain_.minCorner();
+ if (data::particle_flags::isSet(ac.getFlags(p_idx), data::particle_flags::INFINITE))
+ {
+ ac.setNextParticle(p_idx, lc.infiniteParticles_.exchange(int_c(p_idx)));
+ } else
+ {
+ {%- for dim in range(3) %}
+ int hash{{dim}} = static_cast<int>(std::floor((ac.getPosition(p_idx)[{{dim}}] - minCorner[{{dim}}]) * lc.invCellDiameter_[{{dim}}]));
+ {%- endfor %}
+ {%- for dim in range(3) %}
+ if (hash{{dim}} < 0) hash{{dim}} = 0;
+ if (hash{{dim}} >= lc.numCellsPerDim_[{{dim}}]) hash{{dim}} = lc.numCellsPerDim_[{{dim}}] - 1;
+ {%- endfor %}
+ uint64_t cell_idx = getCellIdx(lc, hash0, hash1, hash2);
+ ac.setNextParticle(p_idx, lc.cells_[cell_idx].exchange(int_c(p_idx)));
+ if (ac.getNextParticle(p_idx) == -1)
+ {
+ lc.nonEmptyCells_.emplace_back(cell_idx);
+ }
+ }
+}
+
+} //namespace kernel
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/src/mesa_pd/data/LinkedCells.h b/src/mesa_pd/data/LinkedCells.h
index cdcec8172..ead0d25d0 100644
--- a/src/mesa_pd/data/LinkedCells.h
+++ b/src/mesa_pd/data/LinkedCells.h
@@ -93,7 +93,10 @@ struct LinkedCells
};
inline
-math::AABB getCellAABB(const LinkedCells& ll, const int hash0, const int hash1, const int hash2)
+math::AABB getCellAABB(const LinkedCells& ll,
+ const int64_t hash0,
+ const int64_t hash1,
+ const int64_t hash2)
{
WALBERLA_ASSERT_GREATER_EQUAL(hash0, 0);
WALBERLA_ASSERT_LESS(hash0, ll.numCellsPerDim_[0]);
@@ -112,7 +115,10 @@ math::AABB getCellAABB(const LinkedCells& ll, const int hash0, const int hash1,
}
inline
-int getCellIdx(const LinkedCells& ll, const int hash0, const int hash1, const int hash2)
+uint_t getCellIdx(const LinkedCells& ll,
+ const int64_t hash0,
+ const int64_t hash1,
+ const int64_t hash2)
{
WALBERLA_ASSERT_GREATER_EQUAL(hash0, 0);
WALBERLA_ASSERT_LESS(hash0, ll.numCellsPerDim_[0]);
@@ -120,7 +126,26 @@ int getCellIdx(const LinkedCells& ll, const int hash0, const int hash1, const in
WALBERLA_ASSERT_LESS(hash1, ll.numCellsPerDim_[1]);
WALBERLA_ASSERT_GREATER_EQUAL(hash2, 0);
WALBERLA_ASSERT_LESS(hash2, ll.numCellsPerDim_[2]);
- return hash2 * ll.numCellsPerDim_[1] * ll.numCellsPerDim_[0] + hash1 * ll.numCellsPerDim_[0] + hash0;
+ return uint_c(hash2 * ll.numCellsPerDim_[1] * ll.numCellsPerDim_[0] + hash1 * ll.numCellsPerDim_[0] + hash0);
+}
+
+inline
+void getCellCoordinates(const LinkedCells& ll,
+ const uint64_t idx,
+ int64_t& hash0,
+ int64_t& hash1,
+ int64_t& hash2)
+{
+ hash2 = int64_c(idx) / (ll.numCellsPerDim_[1] * ll.numCellsPerDim_[0]);
+ hash1 = (int64_c(idx) - (hash2 * ll.numCellsPerDim_[1] * ll.numCellsPerDim_[0])) / (ll.numCellsPerDim_[0]);
+ hash0 = int64_c(idx) - hash2 * ll.numCellsPerDim_[1] * ll.numCellsPerDim_[0] - hash1 * ll.numCellsPerDim_[0];
+
+ WALBERLA_ASSERT_GREATER_EQUAL(hash0, 0);
+ WALBERLA_ASSERT_LESS(hash0, ll.numCellsPerDim_[0]);
+ WALBERLA_ASSERT_GREATER_EQUAL(hash1, 0);
+ WALBERLA_ASSERT_LESS(hash1, ll.numCellsPerDim_[1]);
+ WALBERLA_ASSERT_GREATER_EQUAL(hash2, 0);
+ WALBERLA_ASSERT_LESS(hash2, ll.numCellsPerDim_[2]);
}
inline
@@ -153,15 +178,17 @@ LinkedCells::LinkedCells(const math::AABB& domain, const Vec3& cellDiameter)
WALBERLA_CHECK_LESS_EQUAL(cellDiameter_[2], cellDiameter[2]);
WALBERLA_CHECK_GREATER_EQUAL(numCellsPerDim_[2], 0);
+
+ std::fill(cells_.begin(), cells_.end(), -1);
}
void LinkedCells::clear()
{
const uint64_t cellsSize = cells_.size();
//clear existing linked cells
- #ifdef _OPENMP
- #pragma omp parallel for schedule(static)
- #endif
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static)
+#endif
for (int64_t i = 0; i < int64_c(cellsSize); ++i)
cells_[uint64_c(i)] = -1;
infiniteParticles_ = -1;
@@ -180,8 +207,8 @@ inline void LinkedCells::forEachParticlePair(const bool openmp, const Selector&
{
for (int x = 0; x < numCellsPerDim_[0]; ++x)
{
- const int cell_idx = getCellIdx(*this, x, y, z); ///< current cell index
- int p_idx = cells_[uint_c(cell_idx)]; ///< current particle index
+ const uint_t cell_idx = getCellIdx(*this, x, y, z); ///< current cell index
+ int p_idx = cells_[cell_idx]; ///< current particle index
int np_idx = -1; ///< particle to be checked against
while (p_idx != -1)
@@ -218,11 +245,11 @@ inline void LinkedCells::forEachParticlePair(const bool openmp, const Selector&
if (ny >= numCellsPerDim_[1]) continue;
if (nz >= numCellsPerDim_[2]) continue;
- const int ncell_idx = getCellIdx(*this, nx, ny, nz); ///< neighbor cell index
+ const uint_t ncell_idx = getCellIdx(*this, nx, ny, nz); ///< neighbor cell index
WALBERLA_ASSERT_GREATER_EQUAL(p_idx, 0);
WALBERLA_ASSERT_LESS(p_idx, acForLC.size());
- np_idx = cells_[uint_c(ncell_idx)]; ///< neighbor particle index
+ np_idx = cells_[ncell_idx]; ///< neighbor particle index
while (np_idx != -1)
{
WALBERLA_ASSERT_GREATER_EQUAL(np_idx, 0);
@@ -275,8 +302,8 @@ inline void LinkedCells::forEachParticlePairHalf(const bool openmp, const Select
{
for (int x = 0; x < numCellsPerDim_[0]; ++x)
{
- const int cell_idx = getCellIdx(*this, x, y, z); ///< current cell index
- int p_idx = cells_[uint_c(cell_idx)]; ///< current particle index
+ const uint_t cell_idx = getCellIdx(*this, x, y, z); ///< current cell index
+ int p_idx = cells_[cell_idx]; ///< current particle index
int np_idx = -1; ///< particle to be checked against
while (p_idx != -1)
@@ -312,11 +339,11 @@ inline void LinkedCells::forEachParticlePairHalf(const bool openmp, const Select
if (ny >= numCellsPerDim_[1]) continue;
if (nz >= numCellsPerDim_[2]) continue;
- const int ncell_idx = getCellIdx(*this, nx, ny, nz); ///< neighbor cell index
+ const uint_t ncell_idx = getCellIdx(*this, nx, ny, nz); ///< neighbor cell index
WALBERLA_ASSERT_GREATER_EQUAL(p_idx, 0);
WALBERLA_ASSERT_LESS(p_idx, acForLC.size());
- np_idx = cells_[uint_c(ncell_idx)]; ///< neighbor particle index
+ np_idx = cells_[ncell_idx]; ///< neighbor particle index
while (np_idx != -1)
{
WALBERLA_ASSERT_GREATER_EQUAL(np_idx, 0);
diff --git a/src/mesa_pd/data/SparseLinkedCells.h b/src/mesa_pd/data/SparseLinkedCells.h
new file mode 100644
index 000000000..7185c2988
--- /dev/null
+++ b/src/mesa_pd/data/SparseLinkedCells.h
@@ -0,0 +1,400 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file SparseLinkedCells.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/IAccessor.h>
+#include <mesa_pd/data/ParticleStorage.h>
+
+#include <core/Abort.h>
+#include <core/debug/Debug.h>
+#include <core/math/AABB.h>
+#include <stencil/D3Q27.h>
+
+#include <atomic>
+#include <cmath>
+#include <vector>
+
+namespace walberla {
+namespace mesa_pd {
+namespace data {
+
+struct SparseLinkedCells
+{
+ SparseLinkedCells(const math::AABB& domain, const real_t cellDiameter)
+ : SparseLinkedCells(domain, Vec3(cellDiameter,cellDiameter,cellDiameter))
+ {}
+ SparseLinkedCells(const math::AABB& domain, const Vec3& cellDiameter);
+
+ void clear();
+
+ /**
+ * Calls the provided functor \p func for all particle pairs.
+ *
+ * Additional arguments can be provided. No pairs with twice the same particle.
+ * Call syntax for the provided functor
+ * \code
+ * func( *this, i, j, std::forward<Args>(args)... );
+ * \endcode
+ * \param openmp enables/disables OpenMP parallelization of the kernel call
+ */
+ template <typename Selector, typename Accessor, typename Func, typename... Args>
+ void forEachParticlePair(const bool openmp,
+ const Selector& selector,
+ Accessor& acForLC,
+ Func&& func,
+ Args&&... args) const;
+ /**
+ * Calls the provided functor \p func for all particle pairs.
+ *
+ * Additional arguments can be provided. No pairs with twice the same particle are generated.
+ * No pair is called twice!
+ * Call syntax for the provided functor
+ * \code
+ * func( *this, i, j, std::forward<Args>(args)... );
+ * \endcode
+ * \param openmp enables/disables OpenMP parallelization of the kernel call
+ */
+ template <typename Selector, typename Accessor, typename Func, typename... Args>
+ void forEachParticlePairHalf(const bool openmp,
+ const Selector& selector,
+ Accessor& acForLC,
+ Func&& func,
+ Args&&... args) const;
+
+ math::AABB domain_ {}; ///< local domain covered by this data structure
+ Vector3<int> numCellsPerDim_ {}; ///< number of linked cells per dimension
+ Vec3 cellDiameter_ {};
+ Vec3 invCellDiameter_ {};
+ std::atomic<int> infiniteParticles_ {}; ///< data structure for particles to large for the cells
+ std::vector< std::atomic<int> > cells_ {}; ///< actual cell data structure
+ std::vector<size_t> nonEmptyCells_ {}; ///< list of cells containing particles
+};
+
+inline
+math::AABB getCellAABB(const SparseLinkedCells& ll,
+ const int64_t hash0,
+ const int64_t hash1,
+ const int64_t hash2)
+{
+ WALBERLA_ASSERT_GREATER_EQUAL(hash0, 0);
+ WALBERLA_ASSERT_LESS(hash0, ll.numCellsPerDim_[0]);
+ WALBERLA_ASSERT_GREATER_EQUAL(hash1, 0);
+ WALBERLA_ASSERT_LESS(hash1, ll.numCellsPerDim_[1]);
+ WALBERLA_ASSERT_GREATER_EQUAL(hash2, 0);
+ WALBERLA_ASSERT_LESS(hash2, ll.numCellsPerDim_[2]);
+ const auto& minCorner = ll.domain_.minCorner();
+ const real_t xMin = ll.cellDiameter_[0] * real_c(hash0) + minCorner[0];
+ const real_t yMin = ll.cellDiameter_[1] * real_c(hash1) + minCorner[1];
+ const real_t zMin = ll.cellDiameter_[2] * real_c(hash2) + minCorner[2];
+ const real_t xMax = ll.cellDiameter_[0] * real_c(hash0 + 1) + minCorner[0];
+ const real_t yMax = ll.cellDiameter_[1] * real_c(hash1 + 1) + minCorner[1];
+ const real_t zMax = ll.cellDiameter_[2] * real_c(hash2 + 1) + minCorner[2];
+ return math::AABB(xMin, yMin, zMin, xMax, yMax, zMax);
+}
+
+inline
+uint_t getCellIdx(const SparseLinkedCells& ll,
+ const int64_t hash0,
+ const int64_t hash1,
+ const int64_t hash2)
+{
+ WALBERLA_ASSERT_GREATER_EQUAL(hash0, 0);
+ WALBERLA_ASSERT_LESS(hash0, ll.numCellsPerDim_[0]);
+ WALBERLA_ASSERT_GREATER_EQUAL(hash1, 0);
+ WALBERLA_ASSERT_LESS(hash1, ll.numCellsPerDim_[1]);
+ WALBERLA_ASSERT_GREATER_EQUAL(hash2, 0);
+ WALBERLA_ASSERT_LESS(hash2, ll.numCellsPerDim_[2]);
+ return uint_c(hash2 * ll.numCellsPerDim_[1] * ll.numCellsPerDim_[0] + hash1 * ll.numCellsPerDim_[0] + hash0);
+}
+
+inline
+void getCellCoordinates(const SparseLinkedCells& ll,
+ const uint64_t idx,
+ int64_t& hash0,
+ int64_t& hash1,
+ int64_t& hash2)
+{
+ hash2 = int64_c(idx) / (ll.numCellsPerDim_[1] * ll.numCellsPerDim_[0]);
+ hash1 = (int64_c(idx) - (hash2 * ll.numCellsPerDim_[1] * ll.numCellsPerDim_[0])) / (ll.numCellsPerDim_[0]);
+ hash0 = int64_c(idx) - hash2 * ll.numCellsPerDim_[1] * ll.numCellsPerDim_[0] - hash1 * ll.numCellsPerDim_[0];
+
+ WALBERLA_ASSERT_GREATER_EQUAL(hash0, 0);
+ WALBERLA_ASSERT_LESS(hash0, ll.numCellsPerDim_[0]);
+ WALBERLA_ASSERT_GREATER_EQUAL(hash1, 0);
+ WALBERLA_ASSERT_LESS(hash1, ll.numCellsPerDim_[1]);
+ WALBERLA_ASSERT_GREATER_EQUAL(hash2, 0);
+ WALBERLA_ASSERT_LESS(hash2, ll.numCellsPerDim_[2]);
+}
+
+inline
+SparseLinkedCells::SparseLinkedCells(const math::AABB& domain, const Vec3& cellDiameter)
+ : domain_(domain)
+ , numCellsPerDim_( static_cast<int>(std::ceil( domain.sizes()[0] / cellDiameter[0])),
+ static_cast<int>(std::ceil( domain.sizes()[1] / cellDiameter[1])),
+ static_cast<int>(std::ceil( domain.sizes()[2] / cellDiameter[2])) )
+ , cellDiameter_( domain.sizes()[0] / real_c(numCellsPerDim_[0]),
+ domain.sizes()[1] / real_c(numCellsPerDim_[1]),
+ domain.sizes()[2] / real_c(numCellsPerDim_[2]) )
+ , invCellDiameter_( real_t(1) / cellDiameter_[0], real_t(1) / cellDiameter_[1], real_t(1) / cellDiameter_[2] )
+ , cells_(uint_c(numCellsPerDim_[0]*numCellsPerDim_[1]*numCellsPerDim_[2]))
+ , nonEmptyCells_(uint_c(numCellsPerDim_[0]*numCellsPerDim_[1]*numCellsPerDim_[2]))
+{
+ //precondition
+ WALBERLA_CHECK_GREATER_EQUAL(cellDiameter[0], real_t(0));
+ WALBERLA_CHECK_GREATER_EQUAL(cellDiameter[1], real_t(0));
+ WALBERLA_CHECK_GREATER_EQUAL(cellDiameter[2], real_t(0));
+
+ //postcondition
+ WALBERLA_CHECK_GREATER_EQUAL(cellDiameter_[0], real_t(0));
+ WALBERLA_CHECK_LESS_EQUAL(cellDiameter_[0], cellDiameter[0]);
+
+ WALBERLA_CHECK_GREATER_EQUAL(numCellsPerDim_[0], 0);
+ WALBERLA_CHECK_GREATER_EQUAL(cellDiameter_[1], real_t(0));
+ WALBERLA_CHECK_LESS_EQUAL(cellDiameter_[1], cellDiameter[1]);
+
+ WALBERLA_CHECK_GREATER_EQUAL(numCellsPerDim_[1], 0);
+ WALBERLA_CHECK_GREATER_EQUAL(cellDiameter_[2], real_t(0));
+ WALBERLA_CHECK_LESS_EQUAL(cellDiameter_[2], cellDiameter[2]);
+
+ WALBERLA_CHECK_GREATER_EQUAL(numCellsPerDim_[2], 0);
+
+ nonEmptyCells_.clear();
+ std::fill(cells_.begin(), cells_.end(), -1);
+}
+
+void SparseLinkedCells::clear()
+{
+ for (const auto v : nonEmptyCells_)
+ {
+ WALBERLA_ASSERT_LESS(v, cells_.size());
+ cells_[v] = -1;
+ }
+ nonEmptyCells_.clear();
+ infiniteParticles_ = -1;
+}
+template <typename Selector, typename Accessor, typename Func, typename... Args>
+inline void SparseLinkedCells::forEachParticlePair(const bool openmp, const Selector& selector, Accessor& acForLC, Func&& func, Args&&... args) const
+{
+ static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
+ WALBERLA_UNUSED(openmp);
+
+ for (const auto cellIdx : nonEmptyCells_)
+ {
+ int64_t x = 0;
+ int64_t y = 0;
+ int64_t z = 0;
+ getCellCoordinates(*this, cellIdx, x, y, z);
+ int p_idx = cells_[cellIdx]; ///< current particle index
+ int np_idx = -1; ///< particle to be checked against
+
+ while (p_idx != -1)
+ {
+ WALBERLA_ASSERT_GREATER_EQUAL(p_idx, 0);
+ WALBERLA_ASSERT_LESS(p_idx, acForLC.size());
+
+ // check particles in own cell
+ np_idx = acForLC.getNextParticle(uint_c(p_idx)); ///< neighbor particle index
+ while (np_idx != -1)
+ {
+ WALBERLA_ASSERT_GREATER_EQUAL(np_idx, 0);
+ WALBERLA_ASSERT_LESS(np_idx, acForLC.size());
+
+ if (selector(uint_c(p_idx), uint_c(np_idx), acForLC))
+ {
+ func(uint_c(p_idx), uint_c(np_idx), std::forward<Args>(args)...);
+ func(uint_c(np_idx), uint_c(p_idx), std::forward<Args>(args)...);
+ }
+
+ // go to next particle
+ np_idx = acForLC.getNextParticle(uint_c(np_idx));
+ }
+
+ // check particles in infiniteParticles list
+ np_idx = infiniteParticles_; ///< neighbor particle index
+ while (np_idx != -1)
+ {
+ WALBERLA_ASSERT_GREATER_EQUAL(np_idx, 0);
+ WALBERLA_ASSERT_LESS(np_idx, acForLC.size());
+
+ if (selector(uint_c(p_idx), uint_c(np_idx), acForLC))
+ {
+ func(uint_c(p_idx), uint_c(np_idx), std::forward<Args>(args)...);
+ func(uint_c(np_idx), uint_c(p_idx), std::forward<Args>(args)...);
+ }
+
+ // go to next particle
+ np_idx = acForLC.getNextParticle(uint_c(np_idx));
+ }
+
+ // go to next particle
+ p_idx = acForLC.getNextParticle(uint_c(p_idx));
+ }
+
+ // check particles in neighboring cells (only positive ones)
+ for (auto dir : stencil::D3Q27::dir_pos)
+ {
+ const int64_t nx = x + int64_c(stencil::cx[dir]);
+ const int64_t ny = y + int64_c(stencil::cy[dir]);
+ const int64_t nz = z + int64_c(stencil::cz[dir]);
+ if (nx < 0) continue;
+ if (ny < 0) continue;
+ if (nz < 0) continue;
+ if (nx >= numCellsPerDim_[0]) continue;
+ if (ny >= numCellsPerDim_[1]) continue;
+ if (nz >= numCellsPerDim_[2]) continue;
+
+ const uint64_t ncell_idx = getCellIdx(*this, nx, ny, nz); ///< neighbor cell index
+
+ p_idx = cells_[cellIdx]; ///< current particle index
+ WALBERLA_ASSERT_GREATER_EQUAL(p_idx, 0);
+ WALBERLA_ASSERT_LESS(p_idx, acForLC.size());
+ while (p_idx != -1)
+ {
+ np_idx = cells_[ncell_idx]; ///< neighbor particle index
+ while (np_idx != -1)
+ {
+ WALBERLA_ASSERT_GREATER_EQUAL(np_idx, 0);
+ WALBERLA_ASSERT_LESS(np_idx, acForLC.size());
+
+ if (selector(uint_c(p_idx), uint_c(np_idx), acForLC))
+ {
+ func(uint_c(p_idx), uint_c(np_idx), std::forward<Args>(args)...);
+ func(uint_c(np_idx), uint_c(p_idx), std::forward<Args>(args)...);
+ }
+
+ // go to next particle
+ np_idx = acForLC.getNextParticle(uint_c(np_idx));
+ }
+
+ // go to next particle
+ p_idx = acForLC.getNextParticle(uint_c(p_idx));
+ }
+ }
+ }
+}
+template <typename Selector, typename Accessor, typename Func, typename... Args>
+inline void SparseLinkedCells::forEachParticlePairHalf(const bool openmp, const Selector& selector, Accessor& acForLC, Func&& func, Args&&... args) const
+{
+ static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
+ WALBERLA_UNUSED(openmp);
+
+ for (const auto cellIdx : nonEmptyCells_)
+ {
+ int64_t x = 0;
+ int64_t y = 0;
+ int64_t z = 0;
+ getCellCoordinates(*this, cellIdx, x, y, z);
+ int p_idx = cells_[cellIdx]; ///< current particle index
+ int np_idx = -1; ///< particle to be checked against
+
+ while (p_idx != -1)
+ {
+ WALBERLA_ASSERT_GREATER_EQUAL(p_idx, 0);
+ WALBERLA_ASSERT_LESS(p_idx, acForLC.size());
+
+ // check particles in own cell
+ np_idx = acForLC.getNextParticle(uint_c(p_idx)); ///< neighbor particle index
+ while (np_idx != -1)
+ {
+ WALBERLA_ASSERT_GREATER_EQUAL(np_idx, 0);
+ WALBERLA_ASSERT_LESS(np_idx, acForLC.size());
+
+ if (selector(uint_c(p_idx), uint_c(np_idx), acForLC))
+ {
+ func(uint_c(p_idx), uint_c(np_idx), std::forward<Args>(args)...);
+ }
+
+ // go to next particle
+ np_idx = acForLC.getNextParticle(uint_c(np_idx));
+ }
+
+ // check particles in infiniteParticles list
+ np_idx = infiniteParticles_; ///< neighbor particle index
+ while (np_idx != -1)
+ {
+ WALBERLA_ASSERT_GREATER_EQUAL(np_idx, 0);
+ WALBERLA_ASSERT_LESS(np_idx, acForLC.size());
+
+ if (selector(uint_c(p_idx), uint_c(np_idx), acForLC))
+ {
+ func(uint_c(p_idx), uint_c(np_idx), std::forward<Args>(args)...);
+ }
+
+ // go to next particle
+ np_idx = acForLC.getNextParticle(uint_c(np_idx));
+ }
+
+ // go to next particle
+ p_idx = acForLC.getNextParticle(uint_c(p_idx));
+ }
+
+ // check particles in neighboring cells (only positive ones)
+ for (auto dir : stencil::D3Q27::dir_pos)
+ {
+ const int64_t nx = x + int64_c(stencil::cx[dir]);
+ const int64_t ny = y + int64_c(stencil::cy[dir]);
+ const int64_t nz = z + int64_c(stencil::cz[dir]);
+ if (nx < 0) continue;
+ if (ny < 0) continue;
+ if (nz < 0) continue;
+ if (nx >= numCellsPerDim_[0]) continue;
+ if (ny >= numCellsPerDim_[1]) continue;
+ if (nz >= numCellsPerDim_[2]) continue;
+
+ const uint64_t ncell_idx = getCellIdx(*this, nx, ny, nz); ///< neighbor cell index
+
+ p_idx = cells_[cellIdx]; ///< current particle index
+ WALBERLA_ASSERT_GREATER_EQUAL(p_idx, 0);
+ WALBERLA_ASSERT_LESS(p_idx, acForLC.size());
+ while (p_idx != -1)
+ {
+ np_idx = cells_[ncell_idx]; ///< neighbor particle index
+ while (np_idx != -1)
+ {
+ WALBERLA_ASSERT_GREATER_EQUAL(np_idx, 0);
+ WALBERLA_ASSERT_LESS(np_idx, acForLC.size());
+
+ if (selector(uint_c(p_idx), uint_c(np_idx), acForLC))
+ {
+ func(uint_c(p_idx), uint_c(np_idx), std::forward<Args>(args)...);
+ }
+
+ // go to next particle
+ np_idx = acForLC.getNextParticle(uint_c(np_idx));
+ }
+
+ // go to next particle
+ p_idx = acForLC.getNextParticle(uint_c(p_idx));
+ }
+ }
+ }
+}
+
+} //namespace data
+} //namespace mesa_pd
+} //namespace walberla
\ No newline at end of file
diff --git a/src/mesa_pd/kernel/InsertParticleIntoLinkedCells.h b/src/mesa_pd/kernel/InsertParticleIntoLinkedCells.h
index 8ec5b660a..6f693bce6 100644
--- a/src/mesa_pd/kernel/InsertParticleIntoLinkedCells.h
+++ b/src/mesa_pd/kernel/InsertParticleIntoLinkedCells.h
@@ -81,8 +81,8 @@ inline void InsertParticleIntoLinkedCells::operator()(const size_t p_idx, Access
if (hash1 >= lc.numCellsPerDim_[1]) hash1 = lc.numCellsPerDim_[1] - 1;
if (hash2 < 0) hash2 = 0;
if (hash2 >= lc.numCellsPerDim_[2]) hash2 = lc.numCellsPerDim_[2] - 1;
- int cell_idx = getCellIdx(lc, hash0, hash1, hash2);
- ac.setNextParticle(p_idx, lc.cells_[uint_c(cell_idx)].exchange(int_c(p_idx)));
+ uint_t cell_idx = getCellIdx(lc, hash0, hash1, hash2);
+ ac.setNextParticle(p_idx, lc.cells_[cell_idx].exchange(int_c(p_idx)));
}
}
diff --git a/src/mesa_pd/kernel/InsertParticleIntoSparseLinkedCells.h b/src/mesa_pd/kernel/InsertParticleIntoSparseLinkedCells.h
new file mode 100644
index 000000000..a473a6fbb
--- /dev/null
+++ b/src/mesa_pd/kernel/InsertParticleIntoSparseLinkedCells.h
@@ -0,0 +1,95 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file InsertParticleIntoSparseLinkedCells.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/IAccessor.h>
+#include <mesa_pd/data/SparseLinkedCells.h>
+
+#include <vector>
+
+namespace walberla {
+namespace mesa_pd {
+namespace kernel {
+
+/**
+ * Inserts a particle into the data::SparseLinkedCells data structure
+ *
+ * \attention Make sure to data::SparseLinkedCells::clear() the data structure before
+ * reinserting new particles.
+ *
+ * This kernel requires the following particle accessor interface
+ * \code
+ * const walberla::mesa_pd::Vec3& getPosition(const size_t p_idx) const;
+ *
+ * const walberla::mesa_pd::data::particle_flags::FlagT& getFlags(const size_t p_idx) const;
+ *
+ * const size_t& getNextParticle(const size_t p_idx) const;
+ * void setNextParticle(const size_t p_idx, const size_t& v);
+ *
+ * \endcode
+ * \ingroup mesa_pd_kernel
+ */
+class InsertParticleIntoSparseLinkedCells
+{
+public:
+ template <typename Accessor>
+ void operator()(const size_t p_idx, Accessor& ac, data::SparseLinkedCells& lc) const;
+};
+
+template <typename Accessor>
+inline void InsertParticleIntoSparseLinkedCells::operator()(const size_t p_idx, Accessor& ac, data::SparseLinkedCells& lc) const
+{
+ static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
+
+ const auto& minCorner = lc.domain_.minCorner();
+ if (data::particle_flags::isSet(ac.getFlags(p_idx), data::particle_flags::INFINITE))
+ {
+ ac.setNextParticle(p_idx, lc.infiniteParticles_.exchange(int_c(p_idx)));
+ } else
+ {
+ int hash0 = static_cast<int>(std::floor((ac.getPosition(p_idx)[0] - minCorner[0]) * lc.invCellDiameter_[0]));
+ int hash1 = static_cast<int>(std::floor((ac.getPosition(p_idx)[1] - minCorner[1]) * lc.invCellDiameter_[1]));
+ int hash2 = static_cast<int>(std::floor((ac.getPosition(p_idx)[2] - minCorner[2]) * lc.invCellDiameter_[2]));
+ if (hash0 < 0) hash0 = 0;
+ if (hash0 >= lc.numCellsPerDim_[0]) hash0 = lc.numCellsPerDim_[0] - 1;
+ if (hash1 < 0) hash1 = 0;
+ if (hash1 >= lc.numCellsPerDim_[1]) hash1 = lc.numCellsPerDim_[1] - 1;
+ if (hash2 < 0) hash2 = 0;
+ if (hash2 >= lc.numCellsPerDim_[2]) hash2 = lc.numCellsPerDim_[2] - 1;
+ uint64_t cell_idx = getCellIdx(lc, hash0, hash1, hash2);
+ ac.setNextParticle(p_idx, lc.cells_[cell_idx].exchange(int_c(p_idx)));
+ if (ac.getNextParticle(p_idx) == -1)
+ {
+ lc.nonEmptyCells_.emplace_back(cell_idx);
+ }
+ }
+}
+
+} //namespace kernel
+} //namespace mesa_pd
+} //namespace walberla
\ No newline at end of file
diff --git a/tests/mesa_pd/CMakeLists.txt b/tests/mesa_pd/CMakeLists.txt
index e8945e940..cdaf984c4 100644
--- a/tests/mesa_pd/CMakeLists.txt
+++ b/tests/mesa_pd/CMakeLists.txt
@@ -40,9 +40,15 @@ waLBerla_execute_test( NAME MESA_PD_ContactDetection PROCESSES 8 )
waLBerla_compile_test( NAME MESA_PD_Data_Flags FILES data/Flags.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_Data_Flags )
+waLBerla_compile_test( NAME MESA_PD_Data_LinkedCells FILES data/LinkedCells.cpp DEPENDS core )
+waLBerla_execute_test( NAME MESA_PD_Data_LinkedCells )
+
waLBerla_compile_test( NAME MESA_PD_Data_ParticleStorage FILES data/ParticleStorage.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_Data_ParticleStorage )
+waLBerla_compile_test( NAME MESA_PD_Data_SparseLinkedCells FILES data/SparseLinkedCells.cpp DEPENDS core )
+waLBerla_execute_test( NAME MESA_PD_Data_SparseLinkedCells )
+
waLBerla_compile_test( NAME MESA_PD_Domain_BlockForestDomain FILES domain/BlockForestDomain.cpp DEPENDS blockforest core )
waLBerla_execute_test( NAME MESA_PD_Domain_BlockForestDomain )
diff --git a/tests/mesa_pd/data/LinkedCells.cpp b/tests/mesa_pd/data/LinkedCells.cpp
new file mode 100644
index 000000000..0653c116d
--- /dev/null
+++ b/tests/mesa_pd/data/LinkedCells.cpp
@@ -0,0 +1,65 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file LinkedCells.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include <mesa_pd/data/LinkedCells.h>
+
+#include <core/Environment.h>
+#include <core/logging/Logging.h>
+
+#include <algorithm>
+#include <iostream>
+
+namespace walberla {
+namespace mesa_pd {
+
+int main( int argc, char ** argv )
+{
+ Environment env(argc, argv);
+ WALBERLA_UNUSED(env);
+ mpi::MPIManager::instance()->useWorldComm();
+
+ data::LinkedCells lc(math::AABB(real_t(0),real_t(0),real_t(0),
+ real_t(4.5),real_t(5.5),real_t(6.5)),
+ real_t(1));
+
+ WALBERLA_CHECK_EQUAL(lc.numCellsPerDim_[0], 5);
+ WALBERLA_CHECK_EQUAL(lc.numCellsPerDim_[1], 6);
+ WALBERLA_CHECK_EQUAL(lc.numCellsPerDim_[2], 7);
+
+ auto cellIdx = data::getCellIdx(lc, 2, 1, 3);
+ WALBERLA_CHECK_EQUAL(cellIdx, 90 + 5 + 2);
+ int64_t x = 0;
+ int64_t y = 0;
+ int64_t z = 0;
+ data::getCellCoordinates(lc, cellIdx, x, y, z);
+ WALBERLA_CHECK_EQUAL(x, 2);
+ WALBERLA_CHECK_EQUAL(y, 1);
+ WALBERLA_CHECK_EQUAL(z, 3);
+
+ return EXIT_SUCCESS;
+}
+
+} //namespace mesa_pd
+} //namespace walberla
+
+int main( int argc, char ** argv )
+{
+ return walberla::mesa_pd::main(argc, argv);
+}
diff --git a/tests/mesa_pd/data/SparseLinkedCells.cpp b/tests/mesa_pd/data/SparseLinkedCells.cpp
new file mode 100644
index 000000000..090e739e2
--- /dev/null
+++ b/tests/mesa_pd/data/SparseLinkedCells.cpp
@@ -0,0 +1,65 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file SparseLinkedCells.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include <mesa_pd/data/SparseLinkedCells.h>
+
+#include <core/Environment.h>
+#include <core/logging/Logging.h>
+
+#include <algorithm>
+#include <iostream>
+
+namespace walberla {
+namespace mesa_pd {
+
+int main( int argc, char ** argv )
+{
+ Environment env(argc, argv);
+ WALBERLA_UNUSED(env);
+ mpi::MPIManager::instance()->useWorldComm();
+
+ data::SparseLinkedCells lc(math::AABB(real_t(0),real_t(0),real_t(0),
+ real_t(4.5),real_t(5.5),real_t(6.5)),
+ real_t(1));
+
+ WALBERLA_CHECK_EQUAL(lc.numCellsPerDim_[0], 5);
+ WALBERLA_CHECK_EQUAL(lc.numCellsPerDim_[1], 6);
+ WALBERLA_CHECK_EQUAL(lc.numCellsPerDim_[2], 7);
+
+ auto cellIdx = data::getCellIdx(lc, 2, 1, 3);
+ WALBERLA_CHECK_EQUAL(cellIdx, 90 + 5 + 2);
+ int64_t x = 0;
+ int64_t y = 0;
+ int64_t z = 0;
+ data::getCellCoordinates(lc, cellIdx, x, y, z);
+ WALBERLA_CHECK_EQUAL(x, 2);
+ WALBERLA_CHECK_EQUAL(y, 1);
+ WALBERLA_CHECK_EQUAL(z, 3);
+
+ return EXIT_SUCCESS;
+}
+
+} //namespace mesa_pd
+} //namespace walberla
+
+int main( int argc, char ** argv )
+{
+ return walberla::mesa_pd::main(argc, argv);
+}
diff --git a/tests/mesa_pd/kernel/GenerateLinkedCells.cpp b/tests/mesa_pd/kernel/GenerateLinkedCells.cpp
index ef4724b5e..c54473cb0 100644
--- a/tests/mesa_pd/kernel/GenerateLinkedCells.cpp
+++ b/tests/mesa_pd/kernel/GenerateLinkedCells.cpp
@@ -80,9 +80,9 @@ int main( int argc, char ** argv )
for (int y = 0; y < linkedCells.numCellsPerDim_[1]; ++y)
for (int z = 0; z < linkedCells.numCellsPerDim_[2]; ++z)
{
- const int cell_idx = getCellIdx(linkedCells, x, y, z);
+ const uint_t cell_idx = getCellIdx(linkedCells, x, y, z);
auto aabb = getCellAABB(linkedCells, x, y, z);
- int p_idx = linkedCells.cells_[uint_c(cell_idx)];
+ int p_idx = linkedCells.cells_[cell_idx];
while (p_idx != -1)
{
++particleCounter;
--
GitLab