Commit 15bb09e1 authored by Sebastian Eibl's avatar Sebastian Eibl
Browse files

introduced SparseLinkedCells

parent 7d6a8e5a
......@@ -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;
......
......@@ -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/")
......
# -*- coding: utf-8 -*-
import numpy as np
from ..utility import generateFile
class SparseLinkedCells:
def generate(self, path):
generateFile(path, 'data/SparseLinkedCells.templ.h')
......@@ -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']
# -*- 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)
......@@ -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',
......
......@@ -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);
......
//======================================================================================================================
//
// 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
......@@ -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)));
}
}
......
//======================================================================================================================
//
// 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::