From c5d065bb2debc4921fcee464fbf709b814aad7cb Mon Sep 17 00:00:00 2001
From: Samuel Kemmler <samuel.kemmler@fau.de>
Date: Wed, 1 Jun 2022 14:07:04 +0200
Subject: [PATCH] Add the PSM to lbm_mesapd_coupling
---
.clang-format | 28 +-
src/geometry/bodies/DynamicBody.h | 8 +-
src/lbm_mesapd_coupling/CMakeLists.txt | 3 +-
src/lbm_mesapd_coupling/DataTypes.h | 29 +-
.../overlapping/OverlapFraction.h | 117 +++
.../overlapping/shapes/BoxWithOverlap.h | 62 ++
.../shapes/CylindricalBoundaryWithOverlap.h | 64 ++
.../overlapping/shapes/EllipsoidWithOverlap.h | 63 ++
.../overlapping/shapes/HalfSpaceWithOverlap.h | 62 ++
.../overlapping/shapes/SphereWithOverlap.h | 61 ++
.../CMakeLists.txt | 6 +
.../PSMSweep.h | 686 ++++++++++++++
.../PSMUtility.h | 159 ++++
.../ParticleAndVolumeFractionMapping.h | 124 +++
tests/lbm_mesapd_coupling/CMakeLists.txt | 25 +
.../DragForceSphere.cpp | 503 ++++++++++
.../ParticleAndVolumeFractionMapping.cpp | 212 +++++
.../SettlingSphere.cpp | 877 ++++++++++++++++++
.../TorqueSphere.cpp | 561 +++++++++++
.../Utility.h | 85 ++
20 files changed, 3710 insertions(+), 25 deletions(-)
create mode 100644 src/lbm_mesapd_coupling/overlapping/OverlapFraction.h
create mode 100644 src/lbm_mesapd_coupling/overlapping/shapes/BoxWithOverlap.h
create mode 100644 src/lbm_mesapd_coupling/overlapping/shapes/CylindricalBoundaryWithOverlap.h
create mode 100644 src/lbm_mesapd_coupling/overlapping/shapes/EllipsoidWithOverlap.h
create mode 100644 src/lbm_mesapd_coupling/overlapping/shapes/HalfSpaceWithOverlap.h
create mode 100644 src/lbm_mesapd_coupling/overlapping/shapes/SphereWithOverlap.h
create mode 100644 src/lbm_mesapd_coupling/partially_saturated_cells_method/CMakeLists.txt
create mode 100644 src/lbm_mesapd_coupling/partially_saturated_cells_method/PSMSweep.h
create mode 100644 src/lbm_mesapd_coupling/partially_saturated_cells_method/PSMUtility.h
create mode 100644 src/lbm_mesapd_coupling/partially_saturated_cells_method/ParticleAndVolumeFractionMapping.h
create mode 100644 tests/lbm_mesapd_coupling/partially_saturated_cells_method/DragForceSphere.cpp
create mode 100644 tests/lbm_mesapd_coupling/partially_saturated_cells_method/ParticleAndVolumeFractionMapping.cpp
create mode 100644 tests/lbm_mesapd_coupling/partially_saturated_cells_method/SettlingSphere.cpp
create mode 100644 tests/lbm_mesapd_coupling/partially_saturated_cells_method/TorqueSphere.cpp
create mode 100644 tests/lbm_mesapd_coupling/partially_saturated_cells_method/Utility.h
diff --git a/.clang-format b/.clang-format
index 3750f7b5b..5f0ff6558 100644
--- a/.clang-format
+++ b/.clang-format
@@ -79,30 +79,34 @@ IncludeCategories:
Priority: 12
- Regex: '^"lbm/'
Priority: 13
- - Regex: '^"mesh/'
+ - Regex: '^"lbm_mesapd_coupling/'
Priority: 14
- - Regex: '^"pde/'
+ - Regex: '^"mesh/'
Priority: 15
- - Regex: '^"pe/'
+ - Regex: '^"mesa_pd/'
Priority: 16
- - Regex: '^"pe_coupling/'
+ - Regex: '^"pde/'
Priority: 17
- - Regex: '^"postprocessing/'
+ - Regex: '^"pe/'
Priority: 18
- - Regex: '^"python_coupling/'
+ - Regex: '^"pe_coupling/'
Priority: 19
- - Regex: '^"simd/'
+ - Regex: '^"postprocessing/'
Priority: 20
- - Regex: '^"stencil/'
+ - Regex: '^"python_coupling/'
Priority: 21
- - Regex: '^"timeloop/'
+ - Regex: '^"simd/'
Priority: 22
- - Regex: '^"vtk/'
+ - Regex: '^"stencil/'
Priority: 23
- - Regex: '^<boost/'
+ - Regex: '^"timeloop/'
Priority: 24
- - Regex: '^<'
+ - Regex: '^"vtk/'
Priority: 25
+ - Regex: '^<boost/'
+ Priority: 26
+ - Regex: '^<'
+ Priority: 27
IndentCaseLabels: false
IndentPPDirectives: AfterHash
IndentWidth: 3
diff --git a/src/geometry/bodies/DynamicBody.h b/src/geometry/bodies/DynamicBody.h
index 1ebd1670e..bb656cd69 100644
--- a/src/geometry/bodies/DynamicBody.h
+++ b/src/geometry/bodies/DynamicBody.h
@@ -32,8 +32,12 @@ class AbstractBody {
public:
virtual ~AbstractBody() = default;
virtual bool contains (const Vector3<real_t> & point ) const = 0;
- virtual FastOverlapResult fastOverlapCheck ( const Vector3<real_t> & cellMidpoint, const Vector3<real_t> & dx ) const = 0;
- virtual FastOverlapResult fastOverlapCheck ( const AABB & box ) const = 0;
+ virtual FastOverlapResult fastOverlapCheck(const Vector3< real_t >& /*cellMidpoint*/,
+ const Vector3< real_t >& /*dx*/) const
+ {
+ return geometry::DONT_KNOW;
+ }
+ virtual FastOverlapResult fastOverlapCheck(const AABB& /*box*/) const { return geometry::DONT_KNOW; }
};
// this is an adaptor from static to dynamic polymorphism
diff --git a/src/lbm_mesapd_coupling/CMakeLists.txt b/src/lbm_mesapd_coupling/CMakeLists.txt
index 76f78cca8..bde22da3f 100644
--- a/src/lbm_mesapd_coupling/CMakeLists.txt
+++ b/src/lbm_mesapd_coupling/CMakeLists.txt
@@ -14,5 +14,6 @@ target_sources( lbm_mesapd_coupling
add_subdirectory( amr )
add_subdirectory( momentum_exchange_method )
+add_subdirectory( partially_saturated_cells_method )
add_subdirectory( utility )
-add_subdirectory( mapping )
\ No newline at end of file
+add_subdirectory( mapping )
diff --git a/src/lbm_mesapd_coupling/DataTypes.h b/src/lbm_mesapd_coupling/DataTypes.h
index a71680613..f543a1496 100644
--- a/src/lbm_mesapd_coupling/DataTypes.h
+++ b/src/lbm_mesapd_coupling/DataTypes.h
@@ -1,15 +1,15 @@
//======================================================================================================================
//
-// This file is part of waLBerla. waLBerla is free software: you can
+// 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
+// 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
+//
+// 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/>.
//
@@ -21,17 +21,26 @@
#pragma once
-#include "field/GhostLayerField.h"
#include "core/DataTypes.h"
-namespace walberla {
-namespace lbm_mesapd_coupling {
+#include "field/GhostLayerField.h"
+
+namespace walberla
+{
+namespace lbm_mesapd_coupling
+{
/*
* Typedefs specific to the lbm - mesa_pd coupling
*/
using ParticleField_T = walberla::GhostLayerField< walberla::id_t, 1 >;
+namespace psm
+{
+// store the particle uid together with the overlap fraction
+using ParticleAndVolumeFraction_T = std::pair< id_t, real_t >;
+using ParticleAndVolumeFractionField_T = GhostLayerField< std::vector< ParticleAndVolumeFraction_T >, 1 >;
+} // namespace psm
} // namespace lbm_mesapd_coupling
} // namespace walberla
diff --git a/src/lbm_mesapd_coupling/overlapping/OverlapFraction.h b/src/lbm_mesapd_coupling/overlapping/OverlapFraction.h
new file mode 100644
index 000000000..3dbf68e03
--- /dev/null
+++ b/src/lbm_mesapd_coupling/overlapping/OverlapFraction.h
@@ -0,0 +1,117 @@
+//======================================================================================================================
+//
+// 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 OverlapFraction.h
+//! \ingroup lbm_mesapd_coupling
+//! \author Samuel Kemmler <samuel.kemmler@fau.de>
+//! \brief Functor that provides overlap fraction computations for different MESA-PD shapes (used for SingleCast)
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/DataTypes.h"
+
+#include "geometry/bodies/BodyOverlapFunctions.h"
+
+#include "mesa_pd/common/Contains.h"
+#include "mesa_pd/data/DataTypes.h"
+#include "mesa_pd/data/shape/Sphere.h"
+
+#include "shapes/BoxWithOverlap.h"
+#include "shapes/CylindricalBoundaryWithOverlap.h"
+#include "shapes/EllipsoidWithOverlap.h"
+#include "shapes/HalfSpaceWithOverlap.h"
+#include "shapes/SphereWithOverlap.h"
+
+namespace walberla
+{
+namespace lbm_mesapd_coupling
+{
+namespace psm
+{
+
+struct OverlapFractionFunctor
+{
+ template< typename ParticleAccessor_T, typename Shape_T >
+ real_t operator()(const size_t /*particleIdx*/, const Shape_T& /*shape*/,
+ const shared_ptr< ParticleAccessor_T >& /*ac*/, const Vector3< real_t >& /*point*/,
+ const Vector3< real_t >& /*dxVec*/, uint_t /*superSamplingDepth*/)
+ {
+ WALBERLA_ABORT("OverlapFraction not implemented!");
+ }
+
+ template< typename ParticleAccessor_T >
+ real_t operator()(const size_t particleIdx, const mesa_pd::data::Sphere& sphere,
+ const shared_ptr< ParticleAccessor_T >& ac, const Vector3< real_t >& point,
+ const Vector3< real_t >& dxVec, uint_t superSamplingDepth)
+ {
+ WALBERLA_STATIC_ASSERT((std::is_base_of< mesa_pd::data::IAccessor, ParticleAccessor_T >::value));
+
+ SphereWithOverlap sphereWithOverlap(particleIdx, ac, sphere);
+ return geometry::overlapFraction< geometry::AbstractBody >(sphereWithOverlap, point, dxVec, superSamplingDepth);
+ }
+
+ template< typename ParticleAccessor_T >
+ real_t operator()(const size_t particleIdx, const mesa_pd::data::HalfSpace& halfSphere,
+ const shared_ptr< ParticleAccessor_T >& ac, const Vector3< real_t >& point, real_t dxVec,
+ uint_t superSamplingDepth)
+ {
+ WALBERLA_STATIC_ASSERT((std::is_base_of< mesa_pd::data::IAccessor, ParticleAccessor_T >::value));
+
+ HalfSpaceWithOverlap halfSpaceWithOverlap(particleIdx, ac, halfSphere);
+ return geometry::overlapFraction< geometry::AbstractBody >(halfSpaceWithOverlap, point, dxVec,
+ superSamplingDepth);
+ }
+
+ template< typename ParticleAccessor_T >
+ real_t operator()(const size_t particleIdx, const mesa_pd::data::CylindricalBoundary& cylindricalBoundary,
+ const shared_ptr< ParticleAccessor_T >& ac, const Vector3< real_t >& point, real_t dxVec,
+ uint_t superSamplingDepth)
+ {
+ WALBERLA_STATIC_ASSERT((std::is_base_of< mesa_pd::data::IAccessor, ParticleAccessor_T >::value));
+
+ CylindricalBoundaryWithOverlap cylindricalBoundaryWithOverlap(particleIdx, ac, cylindricalBoundary);
+ return geometry::overlapFraction< geometry::AbstractBody >(cylindricalBoundaryWithOverlap, point, dxVec,
+ superSamplingDepth);
+ }
+
+ template< typename ParticleAccessor_T >
+ real_t operator()(const size_t particleIdx, const mesa_pd::data::Box& box,
+ const shared_ptr< ParticleAccessor_T >& ac, const Vector3< real_t >& point, real_t dxVec,
+ uint_t superSamplingDepth)
+ {
+ WALBERLA_STATIC_ASSERT((std::is_base_of< mesa_pd::data::IAccessor, ParticleAccessor_T >::value));
+
+ BoxWithOverlap boxWithOverlap(particleIdx, ac, box);
+ return geometry::overlapFraction< geometry::AbstractBody >(boxWithOverlap, point, dxVec, superSamplingDepth);
+ }
+
+ template< typename ParticleAccessor_T >
+ real_t operator()(const size_t particleIdx, const mesa_pd::data::Ellipsoid& ellipsoid,
+ const shared_ptr< ParticleAccessor_T >& ac, const Vector3< real_t >& point, real_t dxVec,
+ uint_t superSamplingDepth)
+ {
+ WALBERLA_STATIC_ASSERT((std::is_base_of< mesa_pd::data::IAccessor, ParticleAccessor_T >::value));
+
+ EllipsoidWithOverlap ellipsoidWithOverlap(particleIdx, ac, ellipsoid);
+ return geometry::overlapFraction< geometry::AbstractBody >(ellipsoidWithOverlap, point, dxVec,
+ superSamplingDepth);
+ }
+};
+
+} // namespace psm
+} // namespace lbm_mesapd_coupling
+} // namespace walberla
diff --git a/src/lbm_mesapd_coupling/overlapping/shapes/BoxWithOverlap.h b/src/lbm_mesapd_coupling/overlapping/shapes/BoxWithOverlap.h
new file mode 100644
index 000000000..6bf00ca38
--- /dev/null
+++ b/src/lbm_mesapd_coupling/overlapping/shapes/BoxWithOverlap.h
@@ -0,0 +1,62 @@
+//======================================================================================================================
+//
+// 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 BoxWithOverlap.h
+//! \ingroup lbm_mesapd_coupling
+//! \author Samuel Kemmler <samuel.kemmler@fau.de>
+//! \brief Wrapper class that provides a "contains" function for MESA-PD boxes
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/math/Vector3.h"
+
+#include "geometry/bodies/BodyOverlapFunctions.h"
+#include "geometry/bodies/DynamicBody.h"
+
+#include "mesa_pd/common/Contains.h"
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/shape/Box.h"
+
+namespace walberla
+{
+namespace lbm_mesapd_coupling
+{
+namespace psm
+{
+
+template< typename ParticleAccessor_T >
+class BoxWithOverlap : public geometry::AbstractBody
+{
+ public:
+ BoxWithOverlap(const size_t idx, const shared_ptr< ParticleAccessor_T >& ac, const mesa_pd::data::Box& box)
+ : idx_(idx), ac_(ac), box_(box)
+ {}
+
+ bool contains(const Vector3< real_t >& point) const
+ {
+ return mesa_pd::isPointInsideBoxBF(mesa_pd::transformPositionFromWFtoBF(idx_, ac_, point), box_.getEdgeLength());
+ }
+
+ private:
+ size_t idx_;
+ shared_ptr< ParticleAccessor_T > ac_;
+ mesa_pd::data::Box box_;
+};
+
+} // namespace psm
+} // namespace lbm_mesapd_coupling
+} // namespace walberla
diff --git a/src/lbm_mesapd_coupling/overlapping/shapes/CylindricalBoundaryWithOverlap.h b/src/lbm_mesapd_coupling/overlapping/shapes/CylindricalBoundaryWithOverlap.h
new file mode 100644
index 000000000..56a14bfec
--- /dev/null
+++ b/src/lbm_mesapd_coupling/overlapping/shapes/CylindricalBoundaryWithOverlap.h
@@ -0,0 +1,64 @@
+//======================================================================================================================
+//
+// 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 CylindricalBoundaryWithOverlap.h
+//! \ingroup lbm_mesapd_coupling
+//! \author Samuel Kemmler <samuel.kemmler@fau.de>
+//! \brief Wrapper class that provides a "contains" function for MESA-PD cylindrical boundaries
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/math/Vector3.h"
+
+#include "geometry/bodies/BodyOverlapFunctions.h"
+#include "geometry/bodies/DynamicBody.h"
+
+#include "mesa_pd/common/Contains.h"
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/shape/CylindricalBoundary.h"
+
+namespace walberla
+{
+namespace lbm_mesapd_coupling
+{
+namespace psm
+{
+
+template< typename ParticleAccessor_T >
+class CylindricalBoundaryWithOverlap : public geometry::AbstractBody
+{
+ public:
+ CylindricalBoundaryWithOverlap(const size_t idx, const shared_ptr< ParticleAccessor_T >& ac,
+ const mesa_pd::data::CylindricalBoundary& cylindricalBoundary)
+ : idx_(idx), ac_(ac), cylindricalBoundary_(cylindricalBoundary)
+ {}
+
+ bool contains(const Vector3< real_t >& point) const
+ {
+ return mesa_pd::isPointInsideCylindricalBoundary(point, ac_->getPosition(idx_), cylindricalBoundary_.getRadius(),
+ cylindricalBoundary_.getAxis());
+ }
+
+ private:
+ size_t idx_;
+ shared_ptr< ParticleAccessor_T > ac_;
+ mesa_pd::data::CylindricalBoundary cylindricalBoundary_;
+};
+
+} // namespace psm
+} // namespace lbm_mesapd_coupling
+} // namespace walberla
diff --git a/src/lbm_mesapd_coupling/overlapping/shapes/EllipsoidWithOverlap.h b/src/lbm_mesapd_coupling/overlapping/shapes/EllipsoidWithOverlap.h
new file mode 100644
index 000000000..3c8ffe185
--- /dev/null
+++ b/src/lbm_mesapd_coupling/overlapping/shapes/EllipsoidWithOverlap.h
@@ -0,0 +1,63 @@
+//======================================================================================================================
+//
+// 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 EllipsoidWithOverlap.h
+//! \ingroup lbm_mesapd_coupling
+//! \author Samuel Kemmler <samuel.kemmler@fau.de>
+//! \brief Wrapper class that provides a "contains" function for MESA-PD ellipsoids
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/math/Vector3.h"
+
+#include "geometry/bodies/BodyOverlapFunctions.h"
+
+#include "mesa_pd/common/Contains.h"
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/shape/Ellipsoid.h"
+
+namespace walberla
+{
+namespace lbm_mesapd_coupling
+{
+namespace psm
+{
+
+template< typename ParticleAccessor_T >
+class EllipsoidWithOverlap : public geometry::AbstractBody
+{
+ public:
+ EllipsoidWithOverlap(const size_t idx, const shared_ptr< ParticleAccessor_T >& ac,
+ const mesa_pd::data::Ellipsoid& ellipsoid)
+ : idx_(idx), ac_(ac), ellipsoid_(ellipsoid)
+ {}
+
+ bool contains(const Vector3< real_t >& point) const
+ {
+ return mesa_pd::isPointInsideEllipsoidBF(mesa_pd::transformPositionFromWFtoBF(idx_, ac_, point),
+ ellipsoid_.getSemiAxes());
+ }
+
+ private:
+ size_t idx_;
+ shared_ptr< ParticleAccessor_T > ac_;
+ mesa_pd::data::Ellipsoid ellipsoid_;
+};
+
+} // namespace psm
+} // namespace lbm_mesapd_coupling
+} // namespace walberla
diff --git a/src/lbm_mesapd_coupling/overlapping/shapes/HalfSpaceWithOverlap.h b/src/lbm_mesapd_coupling/overlapping/shapes/HalfSpaceWithOverlap.h
new file mode 100644
index 000000000..9f50d6d37
--- /dev/null
+++ b/src/lbm_mesapd_coupling/overlapping/shapes/HalfSpaceWithOverlap.h
@@ -0,0 +1,62 @@
+//======================================================================================================================
+//
+// 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 HalfSpaceWithOverlap.h
+//! \ingroup lbm_mesapd_coupling
+//! \author Samuel Kemmler <samuel.kemmler@fau.de>
+//! \brief Wrapper class that provides a "contains" function for MESA-PD half spaces
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/math/Vector3.h"
+
+#include "geometry/bodies/BodyOverlapFunctions.h"
+
+#include "mesa_pd/common/Contains.h"
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/shape/HalfSpace.h"
+
+namespace walberla
+{
+namespace lbm_mesapd_coupling
+{
+namespace psm
+{
+
+template< typename ParticleAccessor_T >
+class HalfSpaceWithOverlap : public geometry::AbstractBody
+{
+ public:
+ HalfSpaceWithOverlap(const size_t idx, const shared_ptr< ParticleAccessor_T >& ac,
+ const mesa_pd::data::HalfSpace& halfSpace)
+ : idx_(idx), ac_(ac), halfSpace_(halfSpace)
+ {}
+
+ bool contains(const Vector3< real_t >& point) const
+ {
+ return mesa_pd::isPointInsideHalfSpace(point, ac_->getPosition(idx_), halfSpace_.getNormal());
+ }
+
+ private:
+ size_t idx_;
+ shared_ptr< ParticleAccessor_T > ac_;
+ mesa_pd::data::HalfSpace halfSpace_;
+};
+
+} // namespace psm
+} // namespace lbm_mesapd_coupling
+} // namespace walberla
diff --git a/src/lbm_mesapd_coupling/overlapping/shapes/SphereWithOverlap.h b/src/lbm_mesapd_coupling/overlapping/shapes/SphereWithOverlap.h
new file mode 100644
index 000000000..1db48699c
--- /dev/null
+++ b/src/lbm_mesapd_coupling/overlapping/shapes/SphereWithOverlap.h
@@ -0,0 +1,61 @@
+//======================================================================================================================
+//
+// 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 SphereWithOverlap.h
+//! \ingroup lbm_mesapd_coupling
+//! \author Samuel Kemmler <samuel.kemmler@fau.de>
+//! \brief Wrapper class that provides a "contains" function for MESA-PD spheres
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/math/Vector3.h"
+
+#include "geometry/bodies/BodyOverlapFunctions.h"
+
+#include "mesa_pd/common/Contains.h"
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/shape/Sphere.h"
+
+namespace walberla
+{
+namespace lbm_mesapd_coupling
+{
+namespace psm
+{
+
+template< typename ParticleAccessor_T >
+class SphereWithOverlap : public geometry::AbstractBody
+{
+ public:
+ SphereWithOverlap(const size_t idx, const shared_ptr< ParticleAccessor_T >& ac, const mesa_pd::data::Sphere& sphere)
+ : idx_(idx), ac_(ac), sphere_(sphere)
+ {}
+
+ bool contains(const Vector3< real_t >& point) const
+ {
+ return mesa_pd::isPointInsideSphere(point, ac_->getPosition(idx_), sphere_.getRadius());
+ }
+
+ private:
+ size_t idx_;
+ shared_ptr< ParticleAccessor_T > ac_;
+ mesa_pd::data::Sphere sphere_;
+};
+
+} // namespace psm
+} // namespace lbm_mesapd_coupling
+} // namespace walberla
diff --git a/src/lbm_mesapd_coupling/partially_saturated_cells_method/CMakeLists.txt b/src/lbm_mesapd_coupling/partially_saturated_cells_method/CMakeLists.txt
new file mode 100644
index 000000000..65b1c468c
--- /dev/null
+++ b/src/lbm_mesapd_coupling/partially_saturated_cells_method/CMakeLists.txt
@@ -0,0 +1,6 @@
+target_sources( lbm_mesapd_coupling
+ PRIVATE
+ PSMSweep.h
+ ParticleAndVolumeFractionMapping.h
+ PSMUtility.h
+ )
diff --git a/src/lbm_mesapd_coupling/partially_saturated_cells_method/PSMSweep.h b/src/lbm_mesapd_coupling/partially_saturated_cells_method/PSMSweep.h
new file mode 100644
index 000000000..8ad996ba8
--- /dev/null
+++ b/src/lbm_mesapd_coupling/partially_saturated_cells_method/PSMSweep.h
@@ -0,0 +1,686 @@
+//======================================================================================================================
+//
+// 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 PSMSweep.h
+//! \ingroup lbm_mesapd_coupling
+//! \author Samuel Kemmler <samuel.kemmler@fau.de>
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \brief Modification of pe_coupling/partially_saturated_cells_method/PSMSweep.h
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "domain_decomposition/StructuredBlockStorage.h"
+
+#include "field/GhostLayerField.h"
+
+#include "lbm/lattice_model/all.h"
+#include "lbm/sweeps/StreamPull.h"
+#include "lbm/sweeps/SweepBase.h"
+
+#include "lbm_mesapd_coupling/DataTypes.h"
+#include "lbm_mesapd_coupling/utility/ParticleFunctions.h"
+
+#include "mesa_pd/common/ParticleFunctions.h"
+
+namespace walberla
+{
+namespace lbm_mesapd_coupling
+{
+namespace psm
+{
+
+// functions to calculate the cell coverage weighting B
+template< int Weighting_T >
+real_t calculateWeighting(const real_t& /*epsilon*/, const real_t& /*tau*/)
+{
+ WALBERLA_STATIC_ASSERT(Weighting_T == 1 || Weighting_T == 2);
+ exit(EXIT_FAILURE);
+}
+template<>
+inline real_t calculateWeighting< 1 >(const real_t& epsilon, const real_t& /*tau*/)
+{
+ return epsilon;
+}
+template<>
+inline real_t calculateWeighting< 2 >(const real_t& epsilon, const real_t& tau)
+{
+ return epsilon * (tau - real_c(0.5)) / ((real_c(1) - epsilon) + (tau - real_c(0.5)));
+}
+
+/*!\brief LBM sweep for the partially saturated cells method
+ *
+ * Literature:
+ * Original idea: Noble et al. - A lattice-Boltzmann method for partially saturated computational cells, 1998
+ * Overview paper: Owen et al. - An efficient framework for fluid-structure interaction using the lattice Boltzmann
+ * method and immersed moving boundaries, 2010
+ *
+ * Based on the overview paper, at least three different versions for the calculation of the solid collision term
+ * (omega_S) and two different versions for the cell coverage weighting function (B) are available. They can be selected
+ * via the template parameters SolidCollision_T ( = {1,2,3} ) and Weighting_T ( = {1,2} ). SolidCollision_T = 1: Eq. 28
+ * Weighting_T = 1: B = solid volume fraction SolidCollision_T = 2: Eq. 29 Weighting_T = 2: Eq. 30
+ * SolidCollision_T = 3: Eq. 33
+ *
+ * For the calculation of the force acting on the particle, an additional minus sign has to be added compared to Eqs. 31
+ * and 32 in the paper.
+ *
+ */
+template< typename LatticeModel_T, typename Filter_T, typename DensityVelocityIn_T, typename DensityVelocityOut_T,
+ int SolidCollision_T, int Weighting_T, typename ParticleAccessor_T >
+class PSMSweep : public lbm::SweepBase< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut_T >
+{
+ public:
+ using PdfField_T =
+ typename lbm::SweepBase< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut_T >::PdfField_T;
+ using Stencil_T = typename LatticeModel_T::Stencil;
+
+ PSMSweep(const BlockDataID& pdfFieldID, const BlockDataID& particleAndVolumeFractionFieldID,
+ const shared_ptr< StructuredBlockStorage >& blockStorage, const shared_ptr< ParticleAccessor_T >& ac,
+ const Filter_T& _filter = walberla::field::DefaultEvaluationFilter(),
+ const DensityVelocityIn_T& _densityVelocityIn = lbm::DefaultDensityEquilibriumVelocityCalculation(),
+ const DensityVelocityOut_T& _densityVelocityOut = lbm::DefaultDensityVelocityCallback())
+ : lbm::SweepBase< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut_T >(
+ pdfFieldID, _filter, _densityVelocityIn, _densityVelocityOut),
+ particleAndVolumeFractionFieldID_(particleAndVolumeFractionFieldID), blockStorage_(blockStorage), accessor_(ac)
+ {}
+
+ PSMSweep(const BlockDataID& src, const BlockDataID& dst, const BlockDataID& particleAndVolumeFractionFieldID,
+ const shared_ptr< StructuredBlockStorage >& blockStorage, const shared_ptr< ParticleAccessor_T >& ac,
+ const Filter_T& _filter = walberla::field::DefaultEvaluationFilter(),
+ const DensityVelocityIn_T& _densityVelocityIn = lbm::DefaultDensityEquilibriumVelocityCalculation(),
+ const DensityVelocityOut_T& _densityVelocityOut = lbm::DefaultDensityVelocityCallback())
+ : lbm::SweepBase< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut_T >(
+ src, dst, _filter, _densityVelocityIn, _densityVelocityOut),
+ particleAndVolumeFractionFieldID_(particleAndVolumeFractionFieldID), blockStorage_(blockStorage), accessor_(ac)
+ {}
+
+ void operator()(IBlock* const block, const uint_t numberOfGhostLayersToInclude = uint_t(0))
+ {
+ streamCollide(block, numberOfGhostLayersToInclude);
+ }
+
+ void streamCollide(IBlock* const block, const uint_t numberOfGhostLayersToInclude = uint_t(0));
+
+ void stream(IBlock* const block, const uint_t numberOfGhostLayersToInclude = uint_t(0));
+ void collide(IBlock* const block, const uint_t numberOfGhostLayersToInclude = uint_t(0));
+
+ inline ParticleAndVolumeFractionField_T* getParticleAndVolumeFractionField(IBlock* const block) const
+ {
+ WALBERLA_ASSERT_NOT_NULLPTR(block);
+ ParticleAndVolumeFractionField_T* particleAndVolumeFractionField =
+ block->getData< ParticleAndVolumeFractionField_T >(particleAndVolumeFractionFieldID_);
+ WALBERLA_ASSERT_NOT_NULLPTR(particleAndVolumeFractionField);
+ return particleAndVolumeFractionField;
+ }
+
+ inline void getFields(IBlock* const block, PdfField_T*& src, PdfField_T*& dst,
+ ParticleAndVolumeFractionField_T*& particleAndVolumeFractionField)
+ {
+ WALBERLA_ASSERT_NOT_NULLPTR(block);
+
+ src = this->getSrcField(block);
+ dst = this->getDstField(block, src);
+ particleAndVolumeFractionField = getParticleAndVolumeFractionField(block);
+
+ WALBERLA_ASSERT_NOT_NULLPTR(src);
+ WALBERLA_ASSERT_NOT_NULLPTR(dst);
+ WALBERLA_ASSERT_NOT_NULLPTR(particleAndVolumeFractionField);
+
+ WALBERLA_ASSERT_EQUAL(src->xyzSize(), dst->xyzSize());
+ WALBERLA_ASSERT_EQUAL(src->xyzSize(), particleAndVolumeFractionField->xyzSize());
+ }
+
+ private:
+ const BlockDataID particleAndVolumeFractionFieldID_;
+ shared_ptr< StructuredBlockStorage > blockStorage_;
+ const shared_ptr< ParticleAccessor_T > accessor_;
+};
+
+template< typename LatticeModel_T, typename Filter_T, typename DensityVelocityIn_T, typename DensityVelocityOut_T,
+ int SolidCollision_T, int Weighting_T, typename ParticleAccessor_T >
+void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut_T, SolidCollision_T, Weighting_T,
+ ParticleAccessor_T >::streamCollide(IBlock* const block, const uint_t numberOfGhostLayersToInclude)
+{
+ PdfField_T* src(nullptr);
+ PdfField_T* dst(nullptr);
+ ParticleAndVolumeFractionField_T* particleAndVolumeFractionField(nullptr);
+
+ getFields(block, src, dst, particleAndVolumeFractionField);
+
+ const real_t dxCurrentLevel = blockStorage_->dx(blockStorage_->getLevel(*block));
+ const real_t lengthScalingFactor = dxCurrentLevel;
+ const real_t forceScalingFactor = lengthScalingFactor * lengthScalingFactor;
+
+ WALBERLA_ASSERT_GREATER(src->nrOfGhostLayers(), numberOfGhostLayersToInclude);
+ WALBERLA_ASSERT_GREATER_EQUAL(dst->nrOfGhostLayers(), numberOfGhostLayersToInclude);
+
+ const auto& lm = src->latticeModel();
+ dst->resetLatticeModel(
+ lm); /* required so that member functions for getting density and equilibrium velocity can be called for dst! */
+
+ this->filter(*block);
+ this->densityVelocityIn(*block);
+ this->densityVelocityOut(*block);
+
+ const auto& collisionModel = lm.collisionModel();
+
+ const real_t omega = collisionModel.omega();
+ const real_t tau = real_c(1) / omega;
+
+ const cell_idx_t xSize = cell_idx_c(src->xSize());
+ const cell_idx_t ySize = cell_idx_c(src->ySize());
+ const cell_idx_t zSize = cell_idx_c(src->zSize());
+ const cell_idx_t gl = cell_idx_c(numberOfGhostLayersToInclude);
+ for (cell_idx_t z = -gl; z < (zSize + gl); ++z)
+ {
+ for (cell_idx_t y = -gl; y < (ySize + gl); ++y)
+ {
+ for (cell_idx_t x = -gl; x < (xSize + gl); ++x)
+ {
+ if (this->filter(x, y, z))
+ {
+ using namespace stencil;
+
+ real_t pdfs[Stencil_T::Size];
+
+ // stream pull & temporal storage of PDFs
+ for (auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d)
+ {
+ dst->get(x, y, z, d.toIdx()) = src->get(x - d.cx(), y - d.cy(), z - d.cz(), d.toIdx());
+ pdfs[d.toIdx()] = dst->get(x, y, z, d.toIdx());
+ }
+
+ Vector3< real_t > velocity;
+ real_t rho = this->densityVelocityIn(velocity, dst, x, y, z);
+
+ this->densityVelocityOut(x, y, z, lm, velocity, rho);
+
+ // equilibrium distributions
+ auto pdfs_equ = lbm::EquilibriumDistribution< LatticeModel_T >::get(velocity, rho);
+
+ // possible external forcing on fluid
+ const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >(
+ x, y, z, velocity, rho, omega, collisionModel.omega_bulk(), collisionModel.omega_odd());
+
+ // check if particle is present
+ if (particleAndVolumeFractionField->get(x, y, z).size() != size_t(0))
+ {
+ // total coverage ratio in the cell
+ real_t Bn = real_t(0);
+
+ // averaged solid collision operator for all intersecting bodies s
+ // = \sum_s B_s * \Omega_s_i
+ std::vector< real_t > omega_n(Stencil_T::Size, real_t(0));
+
+ // get center of cell
+ Vector3< real_t > cellCenter = blockStorage_->getBlockLocalCellCenter(*block, Cell(x, y, z));
+
+ for (auto particleFracIt = particleAndVolumeFractionField->get(x, y, z).begin();
+ particleFracIt != particleAndVolumeFractionField->get(x, y, z).end(); ++particleFracIt)
+ {
+ real_t omega_s(real_c(0));
+ Vector3< real_t > forceOnParticle(real_c(0));
+
+ real_t Bs = calculateWeighting< Weighting_T >((*particleFracIt).second, tau);
+ Bn += Bs;
+
+ // particle velocity at cell center
+ const size_t idx = accessor_->uidToIdx(particleFracIt->first);
+ WALBERLA_ASSERT_UNEQUAL(idx, accessor_->getInvalidIdx(), "Index of particle is invalid!");
+ const auto particleVelocity = mesa_pd::getVelocityAtWFPoint(idx, *accessor_, cellCenter);
+
+ // equilibrium distributions with solid velocity
+ auto pdfs_equ_solid = lbm::EquilibriumDistribution< LatticeModel_T >::get(particleVelocity, rho);
+
+ for (auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d)
+ {
+ // Different solid collision operators available
+ if (SolidCollision_T == 1)
+ {
+ omega_s =
+ pdfs[d.toInvIdx()] - pdfs_equ[d.toInvIdx()] + pdfs_equ_solid[d.toIdx()] - pdfs[d.toIdx()];
+ }
+ else if (SolidCollision_T == 2)
+ {
+ omega_s = pdfs_equ_solid[d.toIdx()] - pdfs[d.toIdx()] +
+ (real_c(1) - omega) * (pdfs[d.toIdx()] - pdfs_equ[d.toIdx()]);
+ }
+ else if (SolidCollision_T == 3)
+ {
+ omega_s = pdfs[d.toInvIdx()] - pdfs_equ_solid[d.toInvIdx()] + pdfs_equ_solid[d.toIdx()] -
+ pdfs[d.toIdx()];
+ }
+ else
+ {
+ WALBERLA_STATIC_ASSERT(SolidCollision_T >= 1 && SolidCollision_T <= 3);
+ exit(EXIT_FAILURE);
+ }
+ real_t BsOmegaS = Bs * omega_s;
+
+ omega_n[d.toIdx()] += BsOmegaS;
+
+ forceOnParticle[0] -= BsOmegaS * real_c(d.cx());
+ forceOnParticle[1] -= BsOmegaS * real_c(d.cy());
+ forceOnParticle[2] -= BsOmegaS * real_c(d.cz());
+ }
+
+ // scale force when using refinement with (dx)^3 / dt
+ forceOnParticle *= forceScalingFactor;
+
+ // only if cell inside inner domain
+ if (dst->isInInnerPart(Cell(x, y, z)))
+ {
+ // apply force (and automatically torque) on particle
+ lbm_mesapd_coupling::addHydrodynamicForceAtWFPosAtomic(idx, *accessor_, forceOnParticle,
+ cellCenter);
+ }
+ }
+
+ // collide step
+ for (auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d)
+ {
+ // external forcing
+ const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >(
+ x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[d.toIdx()], real_c(d.cx()),
+ real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk(), collisionModel.omega_odd());
+
+ dst->get(x, y, z, d.toIdx()) =
+ pdfs[d.toIdx()] - omega * (real_c(1) - Bn) * (pdfs[d.toIdx()] - pdfs_equ[d.toIdx()]) // SRT
+ + omega_n[d.toIdx()] + (real_c(1) - Bn) * forceTerm;
+ }
+ }
+ else
+ {
+ // SRT collide step
+ for (auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d)
+ {
+ // external forcing
+ const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >(
+ x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[d.toIdx()], real_c(d.cx()),
+ real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk(), collisionModel.omega_odd());
+
+ dst->get(x, y, z, d.toIdx()) =
+ pdfs[d.toIdx()] - omega * (pdfs[d.toIdx()] - pdfs_equ[d.toIdx()]) + forceTerm;
+ }
+ }
+ }
+ }
+ }
+ }
+ src->swapDataPointers(dst);
+}
+
+template< typename LatticeModel_T, typename Filter_T, typename DensityVelocityIn_T, typename DensityVelocityOut_T,
+ int SolidCollision_T, int Weighting_T, typename ParticleAccessor_T >
+void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut_T, SolidCollision_T, Weighting_T,
+ ParticleAccessor_T >::stream(IBlock* const block, const uint_t numberOfGhostLayersToInclude)
+{
+ PdfField_T* src(nullptr);
+ PdfField_T* dst(nullptr);
+ lbm::SweepBase< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut_T >::getFields(block, src, dst);
+ lbm::StreamPull< LatticeModel_T >::execute(src, dst, block, this->filter_, numberOfGhostLayersToInclude);
+}
+
+template< typename LatticeModel_T, typename Filter_T, typename DensityVelocityIn_T, typename DensityVelocityOut_T,
+ int SolidCollision_T, int Weighting_T, typename ParticleAccessor_T >
+void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut_T, SolidCollision_T, Weighting_T,
+ ParticleAccessor_T >::collide(IBlock* const block, const uint_t numberOfGhostLayersToInclude)
+{
+ PdfField_T* src = this->getSrcField(block);
+ ParticleAndVolumeFractionField_T* particleAndVolumeFractionField = getParticleAndVolumeFractionField(block);
+
+ WALBERLA_ASSERT_GREATER(src->nrOfGhostLayers(), numberOfGhostLayersToInclude);
+
+ const real_t dxCurrentLevel = blockStorage_->dx(blockStorage_->getLevel(*block));
+ const real_t lengthScalingFactor = dxCurrentLevel;
+ const real_t forceScalingFactor = lengthScalingFactor * lengthScalingFactor;
+
+ this->filter(*block);
+ this->densityVelocityIn(*block);
+ this->densityVelocityOut(*block);
+
+ const auto& lm = src->latticeModel();
+ const auto& collisionModel = lm.collisionModel();
+
+ const real_t omega = collisionModel.omega();
+ const real_t tau = real_c(1) / omega;
+
+ const cell_idx_t xSize = cell_idx_c(src->xSize());
+ const cell_idx_t ySize = cell_idx_c(src->ySize());
+ const cell_idx_t zSize = cell_idx_c(src->zSize());
+ const cell_idx_t gl = cell_idx_c(numberOfGhostLayersToInclude);
+ for (cell_idx_t z = -gl; z < (zSize + gl); ++z)
+ {
+ for (cell_idx_t y = -gl; y < (ySize + gl); ++y)
+ {
+ for (cell_idx_t x = -gl; x < (xSize + gl); ++x)
+ {
+ if (this->filter(x, y, z))
+ {
+ using namespace stencil;
+
+ real_t pdfs[Stencil_T::Size];
+
+ // temporal storage of PDFs
+ for (auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d)
+ {
+ pdfs[d.toIdx()] = src->get(x, y, z, d.toIdx());
+ }
+
+ Vector3< real_t > velocity;
+ real_t rho = this->densityVelocityIn(velocity, src, x, y, z);
+
+ this->densityVelocityOut(x, y, z, lm, velocity, rho);
+
+ // equilibrium distributions
+ auto pdfs_equ = lbm::EquilibriumDistribution< LatticeModel_T >::get(velocity, rho);
+
+ // possible external forcing on fluid
+ const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >(
+ x, y, z, velocity, rho, omega, collisionModel.omega_bulk(), collisionModel.omega_odd());
+
+ // check if a particle is present
+ if (particleAndVolumeFractionField->get(x, y, z).size() != size_t(0))
+ {
+ // total coverage ratio in the cell
+ real_t Bn = real_t(0);
+
+ // averaged solid collision operator for all intersecting bodies s
+ // = \sum_s B_s * \Omega_s_i
+ std::vector< real_t > omega_n(Stencil_T::Size, real_t(0));
+
+ // get center of cell
+ Vector3< real_t > cellCenter = blockStorage_->getBlockLocalCellCenter(*block, Cell(x, y, z));
+
+ for (auto particleFracIt = particleAndVolumeFractionField->get(x, y, z).begin();
+ particleFracIt != particleAndVolumeFractionField->get(x, y, z).end(); ++particleFracIt)
+ {
+ real_t omega_s(real_c(0));
+ Vector3< real_t > forceOnParticle(real_c(0));
+
+ real_t Bs = calculateWeighting< Weighting_T >((*particleFracIt).second, tau);
+ Bn += Bs;
+
+ // particle velocity at cell center
+ const size_t idx = accessor_->uidToIdx(particleFracIt->first);
+ WALBERLA_ASSERT_UNEQUAL(idx, accessor_->getInvalidIdx(), "Index of particle is invalid!");
+ const auto particleVelocity = mesa_pd::getVelocityAtWFPoint(idx, *accessor_, cellCenter);
+
+ // equilibrium distributions with solid velocity
+ auto pdfs_equ_solid = lbm::EquilibriumDistribution< LatticeModel_T >::get(particleVelocity, rho);
+
+ for (auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d)
+ {
+ // Different solid collision operators available
+ if (SolidCollision_T == 1)
+ {
+ omega_s =
+ pdfs[d.toInvIdx()] - pdfs_equ[d.toInvIdx()] + pdfs_equ_solid[d.toIdx()] - pdfs[d.toIdx()];
+ }
+ else if (SolidCollision_T == 2)
+ {
+ omega_s = pdfs_equ_solid[d.toIdx()] - pdfs[d.toIdx()] +
+ (real_c(1) - omega) * (pdfs[d.toIdx()] - pdfs_equ[d.toIdx()]);
+ }
+ else if (SolidCollision_T == 3)
+ {
+ omega_s = pdfs[d.toInvIdx()] - pdfs_equ_solid[d.toInvIdx()] + pdfs_equ_solid[d.toIdx()] -
+ pdfs[d.toIdx()];
+ }
+ else
+ {
+ WALBERLA_STATIC_ASSERT(SolidCollision_T >= 1 && SolidCollision_T <= 3);
+ exit(EXIT_FAILURE);
+ }
+ real_t BsOmegaS = Bs * omega_s;
+
+ omega_n[d.toIdx()] += BsOmegaS;
+
+ forceOnParticle[0] -= BsOmegaS * real_c(d.cx());
+ forceOnParticle[1] -= BsOmegaS * real_c(d.cy());
+ forceOnParticle[2] -= BsOmegaS * real_c(d.cz());
+ }
+
+ // scale force when using refinement with (dx)^3 / dt
+ forceOnParticle *= forceScalingFactor;
+
+ // only if cell inside inner domain
+ if (src->isInInnerPart(Cell(x, y, z)))
+ {
+ // apply force (and automatically torque) on particle
+ lbm_mesapd_coupling::addHydrodynamicForceAtWFPosAtomic(idx, *accessor_, forceOnParticle,
+ cellCenter);
+ }
+ }
+
+ // collide step
+ for (auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d)
+ {
+ // external forcing
+ const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >(
+ x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[d.toIdx()], real_c(d.cx()),
+ real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk(), collisionModel.omega_odd());
+
+ src->get(x, y, z, d.toIdx()) =
+ pdfs[d.toIdx()] - omega * (real_c(1) - Bn) * (pdfs[d.toIdx()] - pdfs_equ[d.toIdx()]) // SRT
+ + omega_n[d.toIdx()] + (real_c(1) - Bn) * forceTerm;
+ }
+ }
+ else
+ {
+ // SRT collide step
+ for (auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d)
+ {
+ // external forcing
+ const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >(
+ x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[d.toIdx()], real_c(d.cx()),
+ real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk(), collisionModel.omega_odd());
+
+ src->get(x, y, z, d.toIdx()) =
+ pdfs[d.toIdx()] - omega * (pdfs[d.toIdx()] - pdfs_equ[d.toIdx()]) + forceTerm;
+ }
+ }
+ }
+ }
+ }
+ }
+}
+
+/////////////////////////////
+// makePSMSweep FUNCTIONS //
+////////////////////////////
+
+template< typename LatticeModel_T, typename Filter_T, typename DensityVelocityIn_T, typename DensityVelocityOut_T,
+ int SolidCollision_T, int Weighting_T, typename ParticleAccessor_T >
+shared_ptr< PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut_T, SolidCollision_T,
+ Weighting_T, ParticleAccessor_T > >
+ makePSMSweep(const BlockDataID& pdfFieldID, const BlockDataID& particleAndVolumeFractionFieldID,
+ const shared_ptr< StructuredBlockStorage >& blockStorage, const shared_ptr< ParticleAccessor_T >& ac,
+ const Filter_T& filter, const DensityVelocityIn_T& densityVelocityIn,
+ const DensityVelocityOut_T& densityVelocityOut)
+{
+ using PSMS_T = PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut_T, SolidCollision_T,
+ Weighting_T, ParticleAccessor_T >;
+ return shared_ptr< PSMS_T >(new PSMS_T(pdfFieldID, particleAndVolumeFractionFieldID, blockStorage, ac, filter,
+ densityVelocityIn, densityVelocityOut));
+}
+
+template< typename LatticeModel_T, typename Filter_T, typename DensityVelocityIn_T, typename DensityVelocityOut_T,
+ int SolidCollision_T, int Weighting_T, typename ParticleAccessor_T >
+shared_ptr< PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut_T, SolidCollision_T,
+ Weighting_T, ParticleAccessor_T > >
+ makePSMSweep(const BlockDataID& srcID, const BlockDataID& dstID, const BlockDataID& particleAndVolumeFractionFieldID,
+ const shared_ptr< StructuredBlockStorage >& blockStorage, const shared_ptr< ParticleAccessor_T >& ac,
+ const Filter_T& filter, const DensityVelocityIn_T& densityVelocityIn,
+ const DensityVelocityOut_T& densityVelocityOut)
+{
+ using PSMS_T = PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut_T, SolidCollision_T,
+ Weighting_T, ParticleAccessor_T >;
+ return shared_ptr< PSMS_T >(new PSMS_T(srcID, dstID, particleAndVolumeFractionFieldID, blockStorage, ac, filter,
+ densityVelocityIn, densityVelocityOut));
+}
+
+// only block data IDs of PDF data
+
+template< typename LatticeModel_T, int SolidCollision_T, int Weighting_T, typename ParticleAccessor_T >
+shared_ptr< PSMSweep< LatticeModel_T, walberla::field::DefaultEvaluationFilter,
+ lbm::DefaultDensityEquilibriumVelocityCalculation, lbm::DefaultDensityVelocityCallback,
+ SolidCollision_T, Weighting_T, ParticleAccessor_T > >
+ makePSMSweep(const BlockDataID& pdfFieldID, const BlockDataID& particleAndVolumeFractionFieldID,
+ const shared_ptr< StructuredBlockStorage >& blockStorage, const shared_ptr< ParticleAccessor_T >& ac)
+{
+ return makePSMSweep< LatticeModel_T, walberla::field::DefaultEvaluationFilter,
+ lbm::DefaultDensityEquilibriumVelocityCalculation, lbm::DefaultDensityVelocityCallback,
+ SolidCollision_T, Weighting_T >(
+ pdfFieldID, particleAndVolumeFractionFieldID, blockStorage, ac, walberla::field::DefaultEvaluationFilter(),
+ lbm::DefaultDensityEquilibriumVelocityCalculation(), lbm::DefaultDensityVelocityCallback());
+}
+
+template< typename LatticeModel_T, int SolidCollision_T, int Weighting_T, typename ParticleAccessor_T >
+shared_ptr< PSMSweep< LatticeModel_T, walberla::field::DefaultEvaluationFilter,
+ lbm::DefaultDensityEquilibriumVelocityCalculation, lbm::DefaultDensityVelocityCallback,
+ SolidCollision_T, Weighting_T, ParticleAccessor_T > >
+ makePSMSweep(const BlockDataID& srcID, const BlockDataID& dstID, const BlockDataID& particleAndVolumeFractionFieldID,
+ const shared_ptr< StructuredBlockStorage >& blockStorage, const shared_ptr< ParticleAccessor_T >& ac)
+{
+ return makePSMSweep< LatticeModel_T, walberla::field::DefaultEvaluationFilter,
+ lbm::DefaultDensityEquilibriumVelocityCalculation, lbm::DefaultDensityVelocityCallback,
+ SolidCollision_T, Weighting_T >(
+ srcID, dstID, particleAndVolumeFractionFieldID, blockStorage, ac, walberla::field::DefaultEvaluationFilter(),
+ lbm::DefaultDensityEquilibriumVelocityCalculation(), lbm::DefaultDensityVelocityCallback());
+}
+
+// block data IDs of PDF data + flag field as filter
+
+template< typename LatticeModel_T, typename FlagField_T, int SolidCollision_T, int Weighting_T,
+ typename ParticleAccessor_T >
+shared_ptr< PSMSweep< LatticeModel_T, walberla::field::FlagFieldEvaluationFilter< FlagField_T >,
+ lbm::DefaultDensityEquilibriumVelocityCalculation, lbm::DefaultDensityVelocityCallback,
+ SolidCollision_T, Weighting_T, ParticleAccessor_T > >
+ makePSMSweep(const BlockDataID& pdfFieldID, const BlockDataID& particleAndVolumeFractionFieldID,
+ const shared_ptr< StructuredBlockStorage >& blockStorage, const shared_ptr< ParticleAccessor_T >& ac,
+ const ConstBlockDataID& flagFieldID, const Set< FlagUID >& cellsToEvaluate)
+{
+ return makePSMSweep< LatticeModel_T, walberla::field::FlagFieldEvaluationFilter< FlagField_T >,
+ lbm::DefaultDensityEquilibriumVelocityCalculation, lbm::DefaultDensityVelocityCallback,
+ SolidCollision_T, Weighting_T >(
+ pdfFieldID, particleAndVolumeFractionFieldID, blockStorage, ac,
+ walberla::field::FlagFieldEvaluationFilter< FlagField_T >(flagFieldID, cellsToEvaluate),
+ lbm::DefaultDensityEquilibriumVelocityCalculation(), lbm::DefaultDensityVelocityCallback());
+}
+
+template< typename LatticeModel_T, typename FlagField_T, int SolidCollision_T, int Weighting_T,
+ typename ParticleAccessor_T >
+shared_ptr< PSMSweep< LatticeModel_T, walberla::field::FlagFieldEvaluationFilter< FlagField_T >,
+ lbm::DefaultDensityEquilibriumVelocityCalculation, lbm::DefaultDensityVelocityCallback,
+ SolidCollision_T, Weighting_T, ParticleAccessor_T > >
+ makePSMSweep(const BlockDataID& srcID, const BlockDataID& dstID, const BlockDataID& particleAndVolumeFractionFieldID,
+ const shared_ptr< StructuredBlockStorage >& blockStorage, const shared_ptr< ParticleAccessor_T >& ac,
+ const ConstBlockDataID& flagFieldID, const Set< FlagUID >& cellsToEvaluate)
+{
+ return makePSMSweep< LatticeModel_T, walberla::field::FlagFieldEvaluationFilter< FlagField_T >,
+ lbm::DefaultDensityEquilibriumVelocityCalculation, lbm::DefaultDensityVelocityCallback,
+ SolidCollision_T, Weighting_T >(
+ srcID, dstID, particleAndVolumeFractionFieldID, blockStorage, ac,
+ walberla::field::FlagFieldEvaluationFilter< FlagField_T >(flagFieldID, cellsToEvaluate),
+ lbm::DefaultDensityEquilibriumVelocityCalculation(), lbm::DefaultDensityVelocityCallback());
+}
+
+// block data IDs of PDF data + flag field as filter + block data ID of velocity field (out)
+
+template< typename LatticeModel_T, typename FlagField_T, typename VelocityField_T, int SolidCollision_T,
+ int Weighting_T, typename ParticleAccessor_T >
+shared_ptr< PSMSweep< LatticeModel_T, walberla::field::FlagFieldEvaluationFilter< FlagField_T >,
+ lbm::DefaultDensityEquilibriumVelocityCalculation, lbm::VelocityCallback< VelocityField_T >,
+ SolidCollision_T, Weighting_T, ParticleAccessor_T > >
+ makePSMSweep(const BlockDataID& pdfFieldID, const BlockDataID& particleAndVolumeFractionFieldID,
+ const shared_ptr< StructuredBlockStorage >& blockStorage, const shared_ptr< ParticleAccessor_T >& ac,
+ const ConstBlockDataID& flagFieldID, const Set< FlagUID >& cellsToEvaluate,
+ const BlockDataID& velocityFieldID)
+{
+ return makePSMSweep< LatticeModel_T, walberla::field::FlagFieldEvaluationFilter< FlagField_T >,
+ lbm::DefaultDensityEquilibriumVelocityCalculation, lbm::VelocityCallback< VelocityField_T >,
+ SolidCollision_T, Weighting_T >(
+ pdfFieldID, particleAndVolumeFractionFieldID, blockStorage, ac,
+ walberla::field::FlagFieldEvaluationFilter< FlagField_T >(flagFieldID, cellsToEvaluate),
+ lbm::DefaultDensityEquilibriumVelocityCalculation(), lbm::VelocityCallback< VelocityField_T >(velocityFieldID));
+}
+
+template< typename LatticeModel_T, typename FlagField_T, typename VelocityField_T, int SolidCollision_T,
+ int Weighting_T, typename ParticleAccessor_T >
+shared_ptr< PSMSweep< LatticeModel_T, walberla::field::FlagFieldEvaluationFilter< FlagField_T >,
+ lbm::DefaultDensityEquilibriumVelocityCalculation, lbm::VelocityCallback< VelocityField_T >,
+ SolidCollision_T, Weighting_T, ParticleAccessor_T > >
+ makePSMSweep(const BlockDataID& srcID, const BlockDataID& dstID, const BlockDataID& particleAndVolumeFractionFieldID,
+ const shared_ptr< StructuredBlockStorage >& blockStorage, const shared_ptr< ParticleAccessor_T >& ac,
+ const ConstBlockDataID& flagFieldID, const Set< FlagUID >& cellsToEvaluate,
+ const BlockDataID& velocityFieldID)
+{
+ return makePSMSweep< LatticeModel_T, walberla::field::FlagFieldEvaluationFilter< FlagField_T >,
+ lbm::DefaultDensityEquilibriumVelocityCalculation, lbm::VelocityCallback< VelocityField_T >,
+ SolidCollision_T, Weighting_T >(
+ srcID, dstID, particleAndVolumeFractionFieldID, blockStorage, ac,
+ walberla::field::FlagFieldEvaluationFilter< FlagField_T >(flagFieldID, cellsToEvaluate),
+ lbm::DefaultDensityEquilibriumVelocityCalculation(), lbm::VelocityCallback< VelocityField_T >(velocityFieldID));
+}
+
+// block data IDs of PDF data + flag field as filter + block data IDs of velocity and density field (out)
+
+template< typename LatticeModel_T, typename FlagField_T, typename VelocityField_T, typename DensityField_T,
+ int SolidCollision_T, int Weighting_T, typename ParticleAccessor_T >
+shared_ptr< PSMSweep< LatticeModel_T, walberla::field::FlagFieldEvaluationFilter< FlagField_T >,
+ lbm::DefaultDensityEquilibriumVelocityCalculation,
+ lbm::DensityVelocityCallback< VelocityField_T, DensityField_T >, SolidCollision_T, Weighting_T,
+ ParticleAccessor_T > >
+ makePSMSweep(const BlockDataID& pdfFieldID, const BlockDataID& particleAndVolumeFractionFieldID,
+ const shared_ptr< StructuredBlockStorage >& blockStorage, const shared_ptr< ParticleAccessor_T >& ac,
+ const ConstBlockDataID& flagFieldID, const Set< FlagUID >& cellsToEvaluate,
+ const BlockDataID& velocityFieldID, const BlockDataID& densityFieldID)
+{
+ return makePSMSweep< LatticeModel_T, walberla::field::FlagFieldEvaluationFilter< FlagField_T >,
+ lbm::DefaultDensityEquilibriumVelocityCalculation,
+ lbm::DensityVelocityCallback< VelocityField_T, DensityField_T >, SolidCollision_T,
+ Weighting_T >(
+ pdfFieldID, particleAndVolumeFractionFieldID, blockStorage, ac,
+ walberla::field::FlagFieldEvaluationFilter< FlagField_T >(flagFieldID, cellsToEvaluate),
+ lbm::DefaultDensityEquilibriumVelocityCalculation(),
+ lbm::DensityVelocityCallback< VelocityField_T, DensityField_T >(velocityFieldID, densityFieldID));
+}
+
+template< typename LatticeModel_T, typename FlagField_T, typename VelocityField_T, typename DensityField_T,
+ int SolidCollision_T, int Weighting_T, typename ParticleAccessor_T >
+shared_ptr< PSMSweep< LatticeModel_T, walberla::field::FlagFieldEvaluationFilter< FlagField_T >,
+ lbm::DefaultDensityEquilibriumVelocityCalculation,
+ lbm::DensityVelocityCallback< VelocityField_T, DensityField_T >, SolidCollision_T, Weighting_T,
+ ParticleAccessor_T > >
+ makePSMSweep(const BlockDataID& srcID, const BlockDataID& dstID, const BlockDataID& particleAndVolumeFractionFieldID,
+ const shared_ptr< StructuredBlockStorage >& blockStorage, const shared_ptr< ParticleAccessor_T >& ac,
+ const ConstBlockDataID& flagFieldID, const Set< FlagUID >& cellsToEvaluate,
+ const BlockDataID& velocityFieldID, const BlockDataID& densityFieldID)
+{
+ return makePSMSweep< LatticeModel_T, walberla::field::FlagFieldEvaluationFilter< FlagField_T >,
+ lbm::DefaultDensityEquilibriumVelocityCalculation,
+ lbm::DensityVelocityCallback< VelocityField_T, DensityField_T >, SolidCollision_T,
+ Weighting_T >(
+ srcID, dstID, particleAndVolumeFractionFieldID, blockStorage, ac,
+ walberla::field::FlagFieldEvaluationFilter< FlagField_T >(flagFieldID, cellsToEvaluate),
+ lbm::DefaultDensityEquilibriumVelocityCalculation(),
+ lbm::DensityVelocityCallback< VelocityField_T, DensityField_T >(velocityFieldID, densityFieldID));
+}
+
+} // namespace psm
+} // namespace lbm_mesapd_coupling
+} // namespace walberla
diff --git a/src/lbm_mesapd_coupling/partially_saturated_cells_method/PSMUtility.h b/src/lbm_mesapd_coupling/partially_saturated_cells_method/PSMUtility.h
new file mode 100644
index 000000000..854d27f99
--- /dev/null
+++ b/src/lbm_mesapd_coupling/partially_saturated_cells_method/PSMUtility.h
@@ -0,0 +1,159 @@
+//======================================================================================================================
+//
+// 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 PSMUtility.h
+//! \ingroup lbm_mesapd_coupling
+//! \author Samuel Kemmler <samuel.kemmler@fau.de>
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \brief Modification of pe_coupling/partially_saturated_cells_method/PSMUtility.h
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "domain_decomposition/StructuredBlockStorage.h"
+
+#include "field/Field.h"
+#include "field/GhostLayerField.h"
+
+#include "lbm/field/PdfField.h"
+
+namespace walberla
+{
+namespace lbm_mesapd_coupling
+{
+namespace psm
+{
+
+/*!\brief Evaluates the macroscopic velocity for a given cell when using the PSM method.
+ *
+ * returns the cell local macroscopic velocity for the PSM method
+ * cell local velocity = \f$ (1-\sum_s B_s) * u + \sum_s B_s * v_s \f$
+ * u = fluid velocity
+ * v_s = velocity of object s at cell center
+ *
+ * lbm::PdfField< LatticeModel_T > is typically typedef'ed as PdfField_T;
+ * GhostLayerField< std::vector< std::pair< walberla::id_t, real_t > >, 1 > is typically typedef'ed as
+ * ParticleAndVolumeFractionField_T;
+ *
+ * Weighting_T is like in the PSMSweep.
+ */
+template< typename LatticeModel_T, int Weighting_T, typename ParticleAccessor_T >
+Vector3< real_t > getPSMMacroscopicVelocity(
+ const IBlock& block, lbm::PdfField< LatticeModel_T >* pdfField,
+ GhostLayerField< std::vector< std::pair< walberla::size_t, real_t > >, 1 >* particleAndVolumeFractionField,
+ StructuredBlockStorage& blockStorage, const Cell& cell, const ParticleAccessor_T& ac)
+{
+ static_assert(LatticeModel_T::compressible == false, "Only works with incompressible models!");
+ WALBERLA_ASSERT_NOT_NULLPTR(pdfField);
+ WALBERLA_ASSERT_NOT_NULLPTR(particleAndVolumeFractionField);
+
+ const real_t ralaxationTime =
+ real_c(1) / pdfField->latticeModel().collisionModel().omega(cell.x(), cell.y(), cell.z());
+
+ Vector3< real_t > velocity(real_t(0));
+ real_t totalSolidWeightingInCell = real_t(0);
+
+ for (auto particleFracIt = particleAndVolumeFractionField->get(cell).begin();
+ particleFracIt != particleAndVolumeFractionField->get(cell).end(); ++particleFracIt)
+ {
+ const Vector3< real_t > coordsCellCenter = blockStorage.getBlockLocalCellCenter(block, cell);
+
+ const real_t eps = (*particleFracIt).second;
+
+ const real_t Bs = calculateWeighting< Weighting_T >(eps, ralaxationTime);
+
+ const size_t idx = ac.uidToIdx(particleFracIt->first);
+ WALBERLA_ASSERT_UNEQUAL(idx, ac.getInvalidIdx(), "Index of particle is invalid!");
+ velocity += Bs * mesa_pd::getVelocityAtWFPoint(idx, ac, coordsCellCenter);
+
+ totalSolidWeightingInCell += Bs;
+ }
+
+ velocity += pdfField->getVelocity(cell) * (real_c(1) - totalSolidWeightingInCell);
+
+ return velocity;
+}
+
+/*!\brief Initializes the PDF field inside the bodies according to the velocities of the bodies.
+ *
+ * As the Partially Saturated Cells method relies on executing the LBM sweep also inside the bodies, it is good practice
+ * (and for some PSM variants also required) to initialize the PDF field ( i.e. the velocity ) in agreement with
+ * possible initial velocities of the bodies. This is also the case in the presence of external forces acting on the
+ * fluid, as these will often shift the macroscopic velocities during the initialization of the PDF field.
+ *
+ * Note, that the ParticleAndVolumeFractionMapping for PSM has to be called first to have a valid field.
+ *
+ * Only the velocity of cells intersecting with bodies is set, pure fluid cells remain unchanged.
+ */
+template< typename LatticeModel_T, int Weighting_T, typename ParticleAccessor_T >
+void initializeDomainForPSM(StructuredBlockStorage& blockStorage, const BlockDataID& pdfFieldID,
+ const BlockDataID& particleAndVolumeFractionFieldID, const ParticleAccessor_T& ac)
+{
+ using PdfField_T = lbm::PdfField< LatticeModel_T >;
+
+ // iterate all blocks with an iterator 'block'
+ for (auto blockIt = blockStorage.begin(); blockIt != blockStorage.end(); ++blockIt)
+ {
+ // get the field data out of the block
+ PdfField_T* pdfField = blockIt->getData< PdfField_T >(pdfFieldID);
+ ParticleAndVolumeFractionField_T* particleAndVolumeFractionField =
+ blockIt->getData< ParticleAndVolumeFractionField_T >(particleAndVolumeFractionFieldID);
+
+ auto xyzFieldSize = pdfField->xyzSize();
+ for (auto fieldIt = xyzFieldSize.begin(); fieldIt != xyzFieldSize.end(); ++fieldIt)
+ {
+ Cell cell(*fieldIt);
+
+ const real_t ralaxationTime =
+ real_c(1) / pdfField->latticeModel().collisionModel().omega(cell.x(), cell.y(), cell.z());
+
+ Vector3< real_t > weightedAverageParticleVelocityInCell(real_t(0));
+ real_t totalSolidWeightingInCell = real_t(0);
+
+ for (auto particleFracIt = particleAndVolumeFractionField->get(cell).begin();
+ particleFracIt != particleAndVolumeFractionField->get(cell).end(); ++particleFracIt)
+ {
+ const Vector3< real_t > coordsCellCenter = blockStorage.getBlockLocalCellCenter(*blockIt, cell);
+
+ const real_t eps = (*particleFracIt).second;
+
+ const real_t Bs = calculateWeighting< Weighting_T >(eps, ralaxationTime);
+
+ const size_t idx = ac.uidToIdx(particleFracIt->first);
+ WALBERLA_ASSERT_UNEQUAL(idx, ac.getInvalidIdx(), "Index of particle is invalid!");
+ weightedAverageParticleVelocityInCell += Bs * mesa_pd::getVelocityAtWFPoint(idx, ac, coordsCellCenter);
+
+ totalSolidWeightingInCell += Bs;
+ }
+
+ if (totalSolidWeightingInCell > real_t(0))
+ {
+ Vector3< real_t > fluidVelocityInCell(real_t(0));
+ const real_t rho = pdfField->getDensityAndVelocity(fluidVelocityInCell, cell);
+
+ // set the PDFs to equilibrium with the density rho and the average velocity of all intersecting bodies
+ pdfField->setToEquilibrium(cell,
+ fluidVelocityInCell * (real_c(1) - totalSolidWeightingInCell) +
+ weightedAverageParticleVelocityInCell,
+ rho);
+ }
+ }
+ }
+}
+
+} // namespace psm
+} // namespace lbm_mesapd_coupling
+} // namespace walberla
diff --git a/src/lbm_mesapd_coupling/partially_saturated_cells_method/ParticleAndVolumeFractionMapping.h b/src/lbm_mesapd_coupling/partially_saturated_cells_method/ParticleAndVolumeFractionMapping.h
new file mode 100644
index 000000000..021d3d176
--- /dev/null
+++ b/src/lbm_mesapd_coupling/partially_saturated_cells_method/ParticleAndVolumeFractionMapping.h
@@ -0,0 +1,124 @@
+//======================================================================================================================
+//
+// 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 ParticleAndVolumeFractionMapping.h
+//! \ingroup lbm_mesapd_coupling
+//! \author Samuel Kemmler <samuel.kemmler@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "domain_decomposition/StructuredBlockStorage.h"
+
+#include "field/GhostLayerField.h"
+
+#include "lbm_mesapd_coupling/DataTypes.h"
+#include "lbm_mesapd_coupling/mapping/ParticleBoundingBox.h"
+#include "lbm_mesapd_coupling/overlapping/OverlapFraction.h"
+#include "lbm_mesapd_coupling/utility/ParticleSelector.h"
+
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/kernel/SingleCast.h"
+
+#include <functional>
+#include <mesa_pd/data/ParticleStorage.h>
+
+namespace walberla
+{
+namespace lbm_mesapd_coupling
+{
+namespace psm
+{
+
+template< typename ParticleAccessor_T, typename ParticleSelector_T >
+class ParticleAndVolumeFractionMapping
+{
+ public:
+ ParticleAndVolumeFractionMapping(const shared_ptr< StructuredBlockStorage >& blockStorage,
+ const shared_ptr< ParticleAccessor_T >& ac,
+ const ParticleSelector_T& mappingParticleSelector,
+ const BlockDataID& particleAndVolumeFractionFieldID,
+ const uint_t superSamplingDepth = uint_t(4))
+ : blockStorage_(blockStorage), ac_(ac), mappingParticleSelector_(mappingParticleSelector),
+ particleAndVolumeFractionFieldID_(particleAndVolumeFractionFieldID), superSamplingDepth_(superSamplingDepth)
+ {
+ static_assert(std::is_base_of< mesa_pd::data::IAccessor, ParticleAccessor_T >::value,
+ "Provide a valid accessor as template");
+ }
+
+ void operator()()
+ {
+ // clear the field
+ for (auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt)
+ {
+ ParticleAndVolumeFractionField_T* particleAndVolumeFractionField =
+ blockIt->getData< ParticleAndVolumeFractionField_T >(particleAndVolumeFractionFieldID_);
+
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(particleAndVolumeFractionField,
+ particleAndVolumeFractionField->nrOfGhostLayers(),
+ (particleAndVolumeFractionField->get(x, y, z)).clear();)
+ }
+
+ for (size_t idx = 0; idx < ac_->size(); ++idx)
+ {
+ if (mappingParticleSelector_(idx, *ac_)) { update(idx); }
+ }
+ }
+
+ private:
+ void update(const size_t idx)
+ {
+ for (auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt)
+ {
+ ParticleAndVolumeFractionField_T* particleAndVolumeFractionField =
+ blockIt->getData< ParticleAndVolumeFractionField_T >(particleAndVolumeFractionFieldID_);
+
+ CellInterval cellBB =
+ getParticleCellBB(idx, *ac_, *blockIt, *blockStorage_, particleAndVolumeFractionField->nrOfGhostLayers());
+
+ uint_t level = blockStorage_->getLevel(*blockIt);
+ Vector3< real_t > dxVec(blockStorage_->dx(level), blockStorage_->dy(level), blockStorage_->dz(level));
+
+ // compute the overlap fraction for each cell that has to be considered and save the information
+ for (auto cellIt = cellBB.begin(); cellIt != cellBB.end(); ++cellIt)
+ {
+ Cell cell(*cellIt);
+
+ Vector3< real_t > cellCenter;
+ cellCenter = blockStorage_->getBlockLocalCellCenter(*blockIt, cell);
+
+ real_t fraction = singleCast_(idx, *ac_, overlapFractionFctr_, ac_, cellCenter, dxVec, superSamplingDepth_);
+
+ id_t particleUid = ac_->getUid(idx);
+ if (fraction > real_t(0)) { particleAndVolumeFractionField->get(cell).emplace_back(particleUid, fraction); }
+ }
+ }
+ }
+
+ shared_ptr< StructuredBlockStorage > blockStorage_;
+ const shared_ptr< ParticleAccessor_T > ac_;
+ ParticleSelector_T mappingParticleSelector_;
+ const BlockDataID particleAndVolumeFractionFieldID_;
+ const uint_t superSamplingDepth_;
+
+ mesa_pd::kernel::SingleCast singleCast_;
+ OverlapFractionFunctor overlapFractionFctr_;
+};
+
+} // namespace psm
+} // namespace lbm_mesapd_coupling
+} // namespace walberla
diff --git a/tests/lbm_mesapd_coupling/CMakeLists.txt b/tests/lbm_mesapd_coupling/CMakeLists.txt
index eeb208530..94438a2d9 100644
--- a/tests/lbm_mesapd_coupling/CMakeLists.txt
+++ b/tests/lbm_mesapd_coupling/CMakeLists.txt
@@ -47,3 +47,28 @@ waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_UTIL_HydForceMultBlocks_VVAvg CO
waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_UTIL_HydForceMultBlocks_EulerNoAvg COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_UTIL_HydForceMultBlocks> --noForceAveraging PROCESSES 4 )
waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_UTIL_HydForceMultBlocks_VVNoAvg COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_UTIL_HydForceMultBlocks> --noForceAveraging --useVV PROCESSES 4 )
+waLBerla_compile_test( NAME LBM_MESAPD_COUPLING_PSM_BodyAndVolumeFractionMapping FILES partially_saturated_cells_method/ParticleAndVolumeFractionMapping.cpp DEPENDS blockforest boundary core field lbm_mesapd_coupling stencil mesa_pd )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_PSM_BodyAndVolumeFractionMapping PROCESSES 27 )
+
+waLBerla_compile_test( NAME LBM_MESAPD_COUPLING_PSM_DragForceSphere FILES partially_saturated_cells_method/DragForceSphere.cpp DEPENDS blockforest boundary core field lbm lbm_mesapd_coupling timeloop mesa_pd )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_PSM_DragForceSphereFuncTest COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_PSM_DragForceSphere> --funcTest PROCESSES 8 )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_PSM_DragForceSphere COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_PSM_DragForceSphere> PROCESSES 8 LABELS longrun CONFIGURATIONS Release RelWithDbgInfo)
+
+waLBerla_compile_test( NAME LBM_MESAPD_COUPLING_PSM_SettlingSphere FILES partially_saturated_cells_method/SettlingSphere.cpp DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling timeloop mesa_pd )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_PSM_SettlingSphereFuncTest COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_PSM_SettlingSphere> --funcTest PROCESSES 4 )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_PSM_SettlingSphere COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_PSM_SettlingSphere> --resolution 70 PROCESSES 4 LABELS longrun CONFIGURATIONS Release RelWithDbgInfo )
+
+waLBerla_compile_test( NAME LBM_MESAPD_COUPLING_PSM_TorqueSphere FILES partially_saturated_cells_method/TorqueSphere.cpp DEPENDS blockforest boundary core domain_decomposition field lbm stencil timeloop )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_PSM_TorqueSphereSC1W1FuncTest COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_PSM_TorqueSphere> --funcTest --SC1W1 PROCESSES 1 )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_PSM_TorqueSphereSC1W1SingleTest COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_PSM_TorqueSphere> --SC1W1 PROCESSES 1 LABELS longrun CONFIGURATIONS Release RelWithDbgInfo )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_PSM_TorqueSphereSC1W1ParallelTest COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_PSM_TorqueSphere> --SC1W1 PROCESSES 8 LABELS verylongrun CONFIGURATIONS Release RelWithDbgInfo )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_PSM_TorqueSphereSC2W1FuncTest COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_PSM_TorqueSphere> --funcTest --SC2W1 PROCESSES 1 )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_PSM_TorqueSphereSC2W1SingleTest COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_PSM_TorqueSphere> --SC2W1 PROCESSES 1 LABELS longrun CONFIGURATIONS Release RelWithDbgInfo )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_PSM_TorqueSphereSC3W1FuncTest COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_PSM_TorqueSphere> --funcTest --SC3W1 PROCESSES 1 )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_PSM_TorqueSphereSC3W1SingleTest COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_PSM_TorqueSphere> --SC3W1 PROCESSES 1 LABELS longrun CONFIGURATIONS Release RelWithDbgInfo )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_PSM_TorqueSphereSC1W2FuncTest COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_PSM_TorqueSphere> --funcTest --SC1W2 PROCESSES 1 )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_PSM_TorqueSphereSC1W2SingleTest COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_PSM_TorqueSphere> --SC1W2 PROCESSES 1 LABELS longrun CONFIGURATIONS Release RelWithDbgInfo )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_PSM_TorqueSphereSC2W2FuncTest COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_PSM_TorqueSphere> --funcTest --SC2W2 PROCESSES 1 )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_PSM_TorqueSphereSC2W2SingleTest COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_PSM_TorqueSphere> --SC2W2 PROCESSES 1 LABELS longrun CONFIGURATIONS Release RelWithDbgInfo )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_PSM_TorqueSphereSC3W2FuncTest COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_PSM_TorqueSphere> --funcTest --SC3W2 PROCESSES 1 )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_PSM_TorqueSphereSC3W2SingleTest COMMAND $<TARGET_FILE:LBM_MESAPD_COUPLING_PSM_TorqueSphere> --SC3W2 PROCESSES 1 LABELS longrun CONFIGURATIONS Release RelWithDbgInfo )
diff --git a/tests/lbm_mesapd_coupling/partially_saturated_cells_method/DragForceSphere.cpp b/tests/lbm_mesapd_coupling/partially_saturated_cells_method/DragForceSphere.cpp
new file mode 100644
index 000000000..fc7f56dd3
--- /dev/null
+++ b/tests/lbm_mesapd_coupling/partially_saturated_cells_method/DragForceSphere.cpp
@@ -0,0 +1,503 @@
+//======================================================================================================================
+//
+// 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 DragForceSphere.cpp
+//! \ingroup lbm_mesapd_coupling
+//! \author Samuel Kemmler <samuel.kemmler@fau.de>
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \brief Modification of momentum_exchange_method/DragForceSphere.cpp
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "boundary/all.h"
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/SharedFunctor.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/logging/Logging.h"
+#include "core/mpi/MPIManager.h"
+#include "core/mpi/Reduce.h"
+#include "core/timing/RemainingTimeLogger.h"
+
+#include "field/AddToStorage.h"
+
+#include "lbm/communication/PdfFieldPackInfo.h"
+#include "lbm/field/AddToStorage.h"
+#include "lbm/field/PdfField.h"
+#include "lbm/lattice_model/D3Q19.h"
+#include "lbm/lattice_model/ForceModel.h"
+#include "lbm/sweeps/SweepWrappers.h"
+
+#include "lbm_mesapd_coupling/DataTypes.h"
+#include "lbm_mesapd_coupling/momentum_exchange_method/MovingParticleMapping.h"
+#include "lbm_mesapd_coupling/momentum_exchange_method/boundary/SimpleBB.h"
+#include "lbm_mesapd_coupling/partially_saturated_cells_method/PSMSweep.h"
+#include "lbm_mesapd_coupling/partially_saturated_cells_method/PSMUtility.h"
+#include "lbm_mesapd_coupling/partially_saturated_cells_method/ParticleAndVolumeFractionMapping.h"
+#include "lbm_mesapd_coupling/utility/ParticleSelector.h"
+#include "lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
+
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/ShapeStorage.h"
+#include "mesa_pd/domain/BlockForestDomain.h"
+#include "mesa_pd/kernel/ParticleSelector.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+#include <iostream>
+#include <vector>
+
+#include "Utility.h"
+
+namespace drag_force_sphere_psm
+{
+
+///////////
+// USING //
+///////////
+
+using namespace walberla;
+using walberla::uint_t;
+
+using ForceModel_T = lbm::force_model::LuoConstant;
+using LatticeModel_T = lbm::D3Q19< lbm::collision_model::TRT, false, ForceModel_T >;
+
+using Stencil_T = LatticeModel_T::Stencil;
+using PdfField_T = lbm::PdfField< LatticeModel_T >;
+
+////////////////
+// PARAMETERS //
+////////////////
+
+struct Setup
+{
+ uint_t checkFrequency;
+ real_t visc;
+ real_t tau;
+ real_t radius;
+ uint_t length;
+ real_t chi;
+ real_t extForce;
+ real_t analyticalDrag;
+};
+
+template< typename ParticleAccessor_T >
+class DragForceEvaluator
+{
+ public:
+ DragForceEvaluator(SweepTimeloop* timeloop, Setup* setup, const shared_ptr< StructuredBlockStorage >& blocks,
+ const BlockDataID& pdfFieldID, const shared_ptr< ParticleAccessor_T >& ac,
+ walberla::id_t sphereID)
+ : timeloop_(timeloop), setup_(setup), blocks_(blocks), pdfFieldID_(pdfFieldID), ac_(ac), sphereID_(sphereID),
+ normalizedDragOld_(0.0), normalizedDragNew_(0.0)
+ {
+ // calculate the analytical drag force value based on the series expansion of chi
+ // see also Sangani - Slow flow through a periodic array of spheres, IJMF 1982. Eq. 60 and Table 1
+ real_t analyticalDrag = real_c(0);
+ real_t tempChiPowS = real_c(1);
+
+ // coefficients to calculate the drag in a series expansion
+ real_t dragCoefficients[31] = { real_c(1.000000), real_c(1.418649), real_c(2.012564), real_c(2.331523),
+ real_c(2.564809), real_c(2.584787), real_c(2.873609), real_c(3.340163),
+ real_c(3.536763), real_c(3.504092), real_c(3.253622), real_c(2.689757),
+ real_c(2.037769), real_c(1.809341), real_c(1.877347), real_c(1.534685),
+ real_c(0.9034708), real_c(0.2857896), real_c(-0.5512626), real_c(-1.278724),
+ real_c(1.013350), real_c(5.492491), real_c(4.615388), real_c(-0.5736023),
+ real_c(-2.865924), real_c(-4.709215), real_c(-6.870076), real_c(0.1455304),
+ real_c(12.51891), real_c(9.742811), real_c(-5.566269) };
+
+ for (uint_t s = 0; s <= uint_t(30); ++s)
+ {
+ analyticalDrag += dragCoefficients[s] * tempChiPowS;
+ tempChiPowS *= setup->chi;
+ }
+ setup_->analyticalDrag = analyticalDrag;
+ }
+
+ // evaluate the acting drag force
+ void operator()()
+ {
+ const uint_t timestep(timeloop_->getCurrentTimeStep() + 1);
+
+ if (timestep % setup_->checkFrequency != 0) return;
+
+ // get force in x-direction acting on the sphere
+ real_t forceX = computeDragForce();
+ // get average volumetric flowrate in the domain
+ real_t uBar = computeAverageVel();
+
+ // f_total = f_drag + f_buoyancy
+ real_t totalForce =
+ forceX + real_c(4.0 / 3.0) * math::pi * setup_->radius * setup_->radius * setup_->radius * setup_->extForce;
+
+ real_t normalizedDragForce = totalForce / real_c(6.0 * math::pi * setup_->visc * setup_->radius * uBar);
+
+ // update drag force values
+ normalizedDragOld_ = normalizedDragNew_;
+ normalizedDragNew_ = normalizedDragForce;
+ }
+
+ // return the relative temporal change in the normalized drag
+ real_t getDragForceDiff() const { return std::fabs((normalizedDragNew_ - normalizedDragOld_) / normalizedDragNew_); }
+
+ // return the drag force
+ real_t getDragForce() const { return normalizedDragNew_; }
+
+ void logResultToFile(const std::string& filename) const
+ {
+ // write to file if desired
+ // format: length tau viscosity simulatedDrag analyticalDrag\n
+ WALBERLA_ROOT_SECTION()
+ {
+ std::ofstream file;
+ file.open(filename.c_str(), std::ofstream::app);
+ file.precision(8);
+ file << setup_->length << " " << setup_->tau << " " << setup_->visc << " " << normalizedDragNew_ << " "
+ << setup_->analyticalDrag << "\n";
+ file.close();
+ }
+ }
+
+ private:
+ // obtain the drag force acting on the sphere by summing up all the process local parts of fX
+ real_t computeDragForce()
+ {
+ size_t idx = ac_->uidToIdx(sphereID_);
+ WALBERLA_ASSERT_UNEQUAL(idx, ac_->getInvalidIdx(), "Index of particle is invalid!");
+ real_t force = real_t(0);
+ if (idx != ac_->getInvalidIdx()) { force = ac_->getHydrodynamicForce(idx)[0]; }
+
+ WALBERLA_MPI_SECTION() { mpi::allReduceInplace(force, mpi::SUM); }
+
+ return force;
+ }
+
+ // calculate the average velocity in forcing direction (here: x) inside the domain (assuming dx=1)
+ real_t computeAverageVel()
+ {
+ auto velocity_sum = real_t(0);
+ // iterate all blocks stored locally on this process
+ for (auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt)
+ {
+ // retrieve the pdf field and the flag field from the block
+ PdfField_T* pdfField = blockIt->getData< PdfField_T >(pdfFieldID_);
+
+ // get the flag that marks a cell as being fluid
+
+ auto xyzField = pdfField->xyzSize();
+ for (auto cell : xyzField)
+ {
+ velocity_sum += pdfField->getVelocity(cell)[0];
+ }
+ }
+
+ WALBERLA_MPI_SECTION() { mpi::allReduceInplace(velocity_sum, mpi::SUM); }
+
+ return velocity_sum / real_c(setup_->length * setup_->length * setup_->length);
+ }
+
+ SweepTimeloop* timeloop_;
+
+ Setup* setup_;
+
+ shared_ptr< StructuredBlockStorage > blocks_;
+ const BlockDataID pdfFieldID_;
+
+ shared_ptr< ParticleAccessor_T > ac_;
+ const walberla::id_t sphereID_;
+
+ // drag coefficient
+ real_t normalizedDragOld_;
+ real_t normalizedDragNew_;
+};
+
+//////////
+// MAIN //
+//////////
+
+//*******************************************************************************************************************
+/*!\brief Testcase that checks the drag force acting on a fixed sphere in the center of a cubic domain in Stokes flow
+ *
+ * The drag force for this problem (often denoted as Simple Cubic setup) is given by a semi-analytical series expansion.
+ * The cubic domain is periodic in all directions, making it a physically infinite periodic array of spheres.
+ \verbatim
+ _______________
+ ->| |->
+ ->| ___ |->
+ W ->| / \ |-> E
+ E ->| | x | |-> A
+ S ->| \___/ |-> S
+ T ->| |-> T
+ ->|_______________|->
+
+ \endverbatim
+ *
+ * The collision model used for the LBM is TRT with a relaxation parameter tau=1.5 and the magic parameter 3/16.
+ * The Stokes approximation of the equilibrium PDFs is used.
+ * The flow is driven by a constant particle force of 1e-5.
+ * The domain is 32x32x32, and the sphere has a diameter of 16 cells ( chi * domainlength )
+ * The simulation is run until the relative change in the dragforce between 100 time steps is less than 1e-5.
+ * The RPD is not used since the sphere is kept fixed and the force is explicitly reset after each time step.
+ * To avoid periodicity constrain problems, the sphere is set as global.
+ *
+ */
+//*******************************************************************************************************************
+
+int main(int argc, char** argv)
+{
+ debug::enterTestMode();
+
+ mpi::Environment env(argc, argv);
+
+ auto processes = MPIManager::instance()->numProcesses();
+
+ if (processes != 1 && processes != 2 && processes != 4 && processes != 8)
+ {
+ std::cerr << "Number of processes must be equal to either 1, 2, 4, or 8!" << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ ///////////////////
+ // Customization //
+ ///////////////////
+
+ bool shortrun = false;
+ bool funcTest = false;
+ bool logging = false;
+ real_t tau = real_c(1.5);
+ uint_t length = uint_c(32);
+
+ for (int i = 1; i < argc; ++i)
+ {
+ if (std::strcmp(argv[i], "--shortrun") == 0)
+ {
+ shortrun = true;
+ continue;
+ }
+ if (std::strcmp(argv[i], "--funcTest") == 0)
+ {
+ funcTest = true;
+ continue;
+ }
+ if (std::strcmp(argv[i], "--logging") == 0)
+ {
+ logging = true;
+ continue;
+ }
+ if (std::strcmp(argv[i], "--tau") == 0)
+ {
+ tau = real_c(std::atof(argv[++i]));
+ continue;
+ }
+ if (std::strcmp(argv[i], "--length") == 0)
+ {
+ length = uint_c(std::atof(argv[++i]));
+ continue;
+ }
+ WALBERLA_ABORT("Unrecognized command line argument found: " << argv[i]);
+ }
+
+ ///////////////////////////
+ // SIMULATION PROPERTIES //
+ ///////////////////////////
+
+ Setup setup;
+
+ setup.length = length; // length of the cubic domain in lattice cells
+ setup.chi = real_c(0.5); // porosity parameter: diameter / length
+ setup.tau = tau; // relaxation time
+ setup.extForce = real_c(1e-7); // constant particle force in lattice units
+ setup.checkFrequency = uint_t(100); // evaluate the drag force only every checkFrequency time steps
+ setup.radius = real_c(0.5) * setup.chi * real_c(setup.length); // sphere radius
+ setup.visc = (setup.tau - real_c(0.5)) / real_c(3); // viscosity in lattice units
+ const real_t omega = real_c(1) / setup.tau; // relaxation rate
+ const real_t dx = real_c(1); // lattice dx
+ const real_t convergenceLimit = real_c(1e-7); // tolerance for relative change in drag force
+ const uint_t timesteps =
+ funcTest ? 1 : (shortrun ? uint_c(150) : uint_c(50000)); // maximum number of time steps for the whole simulation
+
+ ///////////////////////////
+ // BLOCK STRUCTURE SETUP //
+ ///////////////////////////
+
+ const uint_t XBlocks = (processes >= 2) ? uint_t(2) : uint_t(1);
+ const uint_t YBlocks = (processes >= 4) ? uint_t(2) : uint_t(1);
+ const uint_t ZBlocks = (processes == 8) ? uint_t(2) : uint_t(1);
+ const uint_t XCells = setup.length / XBlocks;
+ const uint_t YCells = setup.length / YBlocks;
+ const uint_t ZCells = setup.length / ZBlocks;
+
+ // create fully periodic domain
+ auto blocks = blockforest::createUniformBlockGrid(XBlocks, YBlocks, ZBlocks, XCells, YCells, ZCells, dx, true, true,
+ true, true);
+
+ /////////
+ // RPD //
+ /////////
+
+ mesa_pd::domain::BlockForestDomain domain(blocks->getBlockForestPointer());
+
+ // init data structures
+ auto ps = std::make_shared< mesa_pd::data::ParticleStorage >(1);
+ auto ss = std::make_shared< mesa_pd::data::ShapeStorage >();
+ using ParticleAccessor_T = mesa_pd::data::ParticleAccessorWithShape;
+ auto accessor = make_shared< ParticleAccessor_T >(ps, ss);
+ auto sphereShape = ss->create< mesa_pd::data::Sphere >(setup.radius);
+
+ //////////////////
+ // RPD COUPLING //
+ //////////////////
+
+ // connect to pe
+ const real_t overlap = real_t(1.5) * dx;
+
+ if (setup.radius > real_c(setup.length) * real_t(0.5) - overlap)
+ {
+ std::cerr << "Periodic sphere is too large and would lead to incorrect mapping!" << std::endl;
+ // solution: create the periodic copies explicitly
+ return EXIT_FAILURE;
+ }
+
+ // create the sphere in the middle of the domain
+ // it is global and thus present on all processes
+ Vector3< real_t > position(real_c(setup.length) * real_c(0.5));
+ walberla::id_t sphereID;
+ {
+ mesa_pd::data::Particle&& p = *ps->create(true);
+ p.setPosition(position);
+ p.setInteractionRadius(setup.radius);
+ p.setOwner(mpi::MPIManager::instance()->rank());
+ p.setShapeID(sphereShape);
+ sphereID = p.getUid();
+ }
+
+ ///////////////////////
+ // ADD DATA TO BLOCKS //
+ ////////////////////////
+
+ // create the lattice model
+ LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::TRT::constructWithMagicNumber(omega),
+ ForceModel_T(Vector3< real_t >(setup.extForce, 0, 0)));
+
+ // add PDF field
+ BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >(
+ blocks, "pdf field (zyxf)", latticeModel, Vector3< real_t >(real_t(0)), real_t(1), uint_t(1), field::zyxf);
+
+ ///////////////
+ // TIME LOOP //
+ ///////////////
+
+ // create the timeloop
+ SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);
+
+ // setup of the LBM communication for synchronizing the pdf field between neighboring blocks
+ blockforest::communication::UniformBufferedScheme< Stencil_T > optimizedPDFCommunicationScheme(blocks);
+ optimizedPDFCommunicationScheme.addPackInfo(
+ make_shared< lbm::PdfFieldPackInfo< LatticeModel_T > >(pdfFieldID)); // optimized sync
+
+ // initially map particles into the LBM simulation
+ BlockDataID particleAndVolumeFractionFieldID =
+ field::addToStorage< lbm_mesapd_coupling::psm::ParticleAndVolumeFractionField_T >(
+ blocks, "particle and volume fraction field",
+ std::vector< lbm_mesapd_coupling::psm::ParticleAndVolumeFraction_T >(), field::zyxf, 0);
+ lbm_mesapd_coupling::psm::ParticleAndVolumeFractionMapping particleMapping(
+ blocks, accessor, lbm_mesapd_coupling::GlobalParticlesSelector(), particleAndVolumeFractionFieldID, 4);
+ particleMapping();
+
+ lbm_mesapd_coupling::psm::initializeDomainForPSM< LatticeModel_T, 1 >(*blocks, pdfFieldID,
+ particleAndVolumeFractionFieldID, *accessor);
+
+ walberla::lbm_mesapd_coupling::FractionFieldSum fractionFieldSum(blocks, particleAndVolumeFractionFieldID);
+ // check that the sum of all fractions is roughly the volume of the particle
+ WALBERLA_CHECK_LESS(
+ std::fabs(4.0 / 3.0 * math::pi * setup.radius * setup.radius * setup.radius - fractionFieldSum()), real_c(1.0));
+
+ // since external forcing is applied, the evaluation of the velocity has to be carried out directly after the
+ // streaming step however, the default sweep is a stream - collide step, i.e. after the sweep, the velocity
+ // evaluation is not correct solution: split the sweep explicitly into collide and stream
+ auto sweep = lbm_mesapd_coupling::psm::makePSMSweep< LatticeModel_T, 1, 1 >(
+ pdfFieldID, particleAndVolumeFractionFieldID, blocks, accessor);
+
+ // collision sweep
+ timeloop.add() << Sweep(lbm::makeCollideSweep(sweep), "cell-wise LB sweep (collide)");
+
+ // add LBM communication function and streaming & force evaluation
+ using DragForceEval_T = DragForceEvaluator< ParticleAccessor_T >;
+ auto forceEval = make_shared< DragForceEval_T >(&timeloop, &setup, blocks, pdfFieldID, accessor, sphereID);
+ timeloop.add() << BeforeFunction(optimizedPDFCommunicationScheme, "LBM Communication")
+ << Sweep(lbm::makeStreamSweep(sweep), "cell-wise LB sweep (stream)")
+ << AfterFunction(SharedFunctor< DragForceEval_T >(forceEval), "drag force evaluation");
+
+ // resetting force
+ timeloop.addFuncAfterTimeStep(
+ [ps, accessor]() {
+ ps->forEachParticle(false, mesa_pd::kernel::SelectAll(), *accessor,
+ lbm_mesapd_coupling::ResetHydrodynamicForceTorqueKernel(), *accessor);
+ },
+ "reset force on sphere");
+
+ timeloop.addFuncAfterTimeStep(RemainingTimeLogger(timeloop.getNrOfTimeSteps()), "Remaining Time Logger");
+
+ ////////////////////////
+ // EXECUTE SIMULATION //
+ ////////////////////////
+
+ WcTimingPool timeloopTiming;
+
+ // time loop
+ for (uint_t i = 0; i < timesteps; ++i)
+ {
+ // perform a single simulation step
+ timeloop.singleStep(timeloopTiming);
+
+ // check if the relative change in the normalized drag force is below the specified convergence criterion
+ if (i > setup.checkFrequency && forceEval->getDragForceDiff() < convergenceLimit)
+ {
+ // if simulation has converged, terminate simulation
+ break;
+ }
+ }
+
+ timeloopTiming.logResultOnRoot();
+
+ if (!funcTest && !shortrun)
+ {
+ // check the result
+ real_t relErr = std::fabs((setup.analyticalDrag - forceEval->getDragForce()) / setup.analyticalDrag);
+ if (logging)
+ {
+ WALBERLA_ROOT_SECTION()
+ {
+ std::cout << "Analytical drag: " << setup.analyticalDrag << "\n"
+ << "Simulated drag: " << forceEval->getDragForce() << "\n"
+ << "Relative error: " << relErr << "\n";
+ }
+ forceEval->logResultToFile("log_DragForceSphere.txt");
+ }
+ // the relative error has to be below 10%
+ WALBERLA_CHECK_LESS(relErr, real_c(0.1));
+ }
+
+ return 0;
+}
+
+} // namespace drag_force_sphere_psm
+
+int main(int argc, char** argv) { drag_force_sphere_psm::main(argc, argv); }
diff --git a/tests/lbm_mesapd_coupling/partially_saturated_cells_method/ParticleAndVolumeFractionMapping.cpp b/tests/lbm_mesapd_coupling/partially_saturated_cells_method/ParticleAndVolumeFractionMapping.cpp
new file mode 100644
index 000000000..eb03d3537
--- /dev/null
+++ b/tests/lbm_mesapd_coupling/partially_saturated_cells_method/ParticleAndVolumeFractionMapping.cpp
@@ -0,0 +1,212 @@
+//======================================================================================================================
+//
+// 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 ParticleAndVolumeFractionMapping.cpp
+//! \ingroup lbm_mesapd_coupling
+//! \author Samuel Kemmler <samuel.kemmler@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/debug/TestSubsystem.h"
+
+#include "field/AddToStorage.h"
+
+#include "lbm_mesapd_coupling/partially_saturated_cells_method/ParticleAndVolumeFractionMapping.h"
+#include "lbm_mesapd_coupling/utility/ParticleSelector.h"
+
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/ShapeStorage.h"
+#include "mesa_pd/domain/BlockForestDomain.h"
+#include "mesa_pd/kernel/SemiImplicitEuler.h"
+#include "mesa_pd/mpi/SyncNextNeighbors.h"
+
+#include <functional>
+#include <memory>
+
+#include "Utility.h"
+
+namespace particle_volume_fraction_check
+{
+
+///////////
+// USING //
+///////////
+
+using namespace walberla;
+using walberla::uint_t;
+using namespace lbm_mesapd_coupling;
+
+////////////////
+// Parameters //
+////////////////
+
+struct Setup
+{
+ // domain size (in lattice cells) in x, y and z direction
+ uint_t xlength;
+ uint_t ylength;
+ uint_t zlength;
+
+ // number of block in x, y and z, direction
+ Vector3< uint_t > nBlocks;
+
+ // cells per block in x, y and z direction
+ Vector3< uint_t > cellsPerBlock;
+
+ real_t sphereDiam;
+
+ uint_t timesteps;
+};
+
+//////////
+// MAIN //
+//////////
+
+//*******************************************************************************************************************
+/*!\brief Testcase that checks if ParticleAndVolumeFractionMapping.h works as intended
+ *
+ * A sphere particle is placed inside the domain and is moving with a constant velocity. The overlap fraction is
+ * computed for all cells in each time step. If the mapping is correct, the sum over all fractions should be roughly
+ * equivalent to the volume of the sphere.
+ *
+ */
+//*******************************************************************************************************************
+
+int main(int argc, char** argv)
+{
+ debug::enterTestMode();
+
+ mpi::Environment env(argc, argv);
+
+ auto processes = MPIManager::instance()->numProcesses();
+
+ if (processes != 27)
+ {
+ std::cerr << "Number of processes must be 27!" << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ ///////////////////////////
+ // SIMULATION PROPERTIES //
+ ///////////////////////////
+
+ Setup setup;
+
+ setup.sphereDiam = real_c(12);
+ setup.zlength = uint_c(4 * setup.sphereDiam);
+ setup.xlength = setup.zlength;
+ setup.ylength = setup.zlength;
+
+ const real_t sphereRadius = real_c(0.5) * setup.sphereDiam;
+ const real_t dx = real_c(1);
+
+ setup.timesteps = 100;
+
+ ///////////////////////////
+ // BLOCK STRUCTURE SETUP //
+ ///////////////////////////
+
+ setup.nBlocks[0] = uint_c(3);
+ setup.nBlocks[1] = uint_c(3);
+ setup.nBlocks[2] = uint_c(3);
+ setup.cellsPerBlock[0] = setup.xlength / setup.nBlocks[0];
+ setup.cellsPerBlock[1] = setup.ylength / setup.nBlocks[1];
+ setup.cellsPerBlock[2] = setup.zlength / setup.nBlocks[2];
+
+ auto blocks =
+ blockforest::createUniformBlockGrid(setup.nBlocks[0], setup.nBlocks[1], setup.nBlocks[2], setup.cellsPerBlock[0],
+ setup.cellsPerBlock[1], setup.cellsPerBlock[2], dx, true, true, true, true);
+
+ ////////////
+ // MesaPD //
+ ////////////
+
+ auto mesapdDomain = std::make_shared< mesa_pd::domain::BlockForestDomain >(blocks->getBlockForestPointer());
+ auto ps = std::make_shared< mesa_pd::data::ParticleStorage >(1);
+ auto ss = std::make_shared< mesa_pd::data::ShapeStorage >();
+ using ParticleAccessor_T = mesa_pd::data::ParticleAccessorWithShape;
+ auto accessor = walberla::make_shared< ParticleAccessor_T >(ps, ss);
+
+ // set up synchronization
+ std::function< void(void) > syncCall = [&]() {
+ const real_t overlap = real_t(1.5) * dx;
+ mesa_pd::mpi::SyncNextNeighbors syncNextNeighborFunc;
+ syncNextNeighborFunc(*ps, *mesapdDomain, overlap);
+ };
+
+ // add the sphere in the center of the domain
+ Vector3< real_t > position(real_c(setup.xlength) * real_c(0.5), real_c(setup.ylength) * real_c(0.5),
+ real_c(setup.zlength) * real_c(0.5));
+ Vector3< real_t > velocity(real_c(0.1), real_c(0.1), real_c(0.1));
+ auto sphereShape = ss->create< mesa_pd::data::Sphere >(sphereRadius);
+
+ if (mesapdDomain->isContainedInProcessSubdomain(uint_c(walberla::mpi::MPIManager::instance()->rank()), position))
+ {
+ auto sphereParticle = ps->create();
+
+ sphereParticle->setShapeID(sphereShape);
+ sphereParticle->setType(0);
+ sphereParticle->setPosition(position);
+ sphereParticle->setLinearVelocity(velocity);
+ sphereParticle->setOwner(walberla::MPIManager::instance()->rank());
+ sphereParticle->setInteractionRadius(sphereRadius);
+ }
+
+ syncCall();
+
+ ///////////////////////
+ // ADD DATA TO BLOCKS //
+ ////////////////////////
+
+ // add particle and volume fraction field (needed for the PSM)
+ BlockDataID particleAndVolumeFractionFieldID =
+ field::addToStorage< lbm_mesapd_coupling::psm::ParticleAndVolumeFractionField_T >(
+ blocks, "particle and volume fraction field",
+ std::vector< lbm_mesapd_coupling::psm::ParticleAndVolumeFraction_T >(), field::zyxf, 0);
+
+ // calculate fraction
+ lbm_mesapd_coupling::psm::ParticleAndVolumeFractionMapping particleMapping(
+ blocks, accessor, lbm_mesapd_coupling::RegularParticlesSelector(), particleAndVolumeFractionFieldID, 4);
+ particleMapping();
+
+ FractionFieldSum fractionFieldSum(blocks, particleAndVolumeFractionFieldID);
+ auto selector = mesa_pd::kernel::SelectMaster();
+ mesa_pd::kernel::SemiImplicitEuler particleIntegration(1.0);
+
+ for (uint_t i = 0; i < setup.timesteps; ++i)
+ {
+ // check that the sum over all fractions is roughly the volume of the sphere
+ real_t sum = fractionFieldSum();
+ WALBERLA_CHECK_LESS(std::fabs(4.0 / 3.0 * math::pi * sphereRadius * sphereRadius * sphereRadius - sum),
+ real_c(1.0));
+
+ // update position
+ ps->forEachParticle(false, selector, *accessor, particleIntegration, *accessor);
+ syncCall();
+ // update fraction mapping
+ particleMapping();
+ }
+
+ return EXIT_SUCCESS;
+}
+
+} // namespace particle_volume_fraction_check
+
+int main(int argc, char** argv) { particle_volume_fraction_check::main(argc, argv); }
diff --git a/tests/lbm_mesapd_coupling/partially_saturated_cells_method/SettlingSphere.cpp b/tests/lbm_mesapd_coupling/partially_saturated_cells_method/SettlingSphere.cpp
new file mode 100644
index 000000000..7fb9d96b4
--- /dev/null
+++ b/tests/lbm_mesapd_coupling/partially_saturated_cells_method/SettlingSphere.cpp
@@ -0,0 +1,877 @@
+//======================================================================================================================
+//
+// 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 SettlingSphere.cpp
+//! \ingroup lbm_mesapd_coupling
+//! \author Samuel Kemmler <samuel.kemmler@fau.de>
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \brief Modification of momentum_exchange_method/SettlingSphere.cpp
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "boundary/all.h"
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/debug/Debug.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/logging/all.h"
+#include "core/math/all.h"
+#include "core/timing/RemainingTimeLogger.h"
+
+#include "domain_decomposition/SharedSweep.h"
+
+#include "field/AddToStorage.h"
+#include "field/communication/PackInfo.h"
+#include "field/vtk/all.h"
+
+#include "lbm/boundary/all.h"
+#include "lbm/communication/PdfFieldPackInfo.h"
+#include "lbm/field/AddToStorage.h"
+#include "lbm/field/PdfField.h"
+#include "lbm/lattice_model/D3Q19.h"
+#include "lbm/sweeps/SweepWrappers.h"
+#include "lbm/vtk/all.h"
+
+#include "lbm_mesapd_coupling/DataTypes.h"
+#include "lbm_mesapd_coupling/mapping/ParticleMapping.h"
+#include "lbm_mesapd_coupling/partially_saturated_cells_method/PSMSweep.h"
+#include "lbm_mesapd_coupling/partially_saturated_cells_method/PSMUtility.h"
+#include "lbm_mesapd_coupling/partially_saturated_cells_method/ParticleAndVolumeFractionMapping.h"
+#include "lbm_mesapd_coupling/utility/AddForceOnParticlesKernel.h"
+#include "lbm_mesapd_coupling/utility/AddHydrodynamicInteractionKernel.h"
+#include "lbm_mesapd_coupling/utility/AverageHydrodynamicForceTorqueKernel.h"
+#include "lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h"
+#include "lbm_mesapd_coupling/utility/LubricationCorrectionKernel.h"
+#include "lbm_mesapd_coupling/utility/ParticleSelector.h"
+#include "lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
+
+#include "mesa_pd/collision_detection/AnalyticContactDetection.h"
+#include "mesa_pd/data/DataTypes.h"
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/ShapeStorage.h"
+#include "mesa_pd/data/shape/HalfSpace.h"
+#include "mesa_pd/data/shape/Sphere.h"
+#include "mesa_pd/domain/BlockForestDomain.h"
+#include "mesa_pd/kernel/DoubleCast.h"
+#include "mesa_pd/kernel/ExplicitEuler.h"
+#include "mesa_pd/kernel/ParticleSelector.h"
+#include "mesa_pd/kernel/SpringDashpot.h"
+#include "mesa_pd/kernel/VelocityVerlet.h"
+#include "mesa_pd/mpi/ContactFilter.h"
+#include "mesa_pd/mpi/ReduceProperty.h"
+#include "mesa_pd/mpi/SyncNextNeighbors.h"
+#include "mesa_pd/mpi/notifications/ForceTorqueNotification.h"
+#include "mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h"
+#include "mesa_pd/vtk/ParticleVtkOutput.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+#include "vtk/all.h"
+
+#include <functional>
+
+#include "Utility.h"
+
+namespace settling_sphere
+{
+
+///////////
+// USING //
+///////////
+
+using namespace walberla;
+using walberla::uint_t;
+
+using LatticeModel_T = lbm::D3Q19< lbm::collision_model::TRT >;
+
+using Stencil_T = LatticeModel_T::Stencil;
+using PdfField_T = lbm::PdfField< LatticeModel_T >;
+
+using flag_t = walberla::uint8_t;
+using FlagField_T = FlagField< flag_t >;
+
+const uint_t FieldGhostLayers = 1;
+
+///////////
+// FLAGS //
+///////////
+
+const FlagUID Fluid_Flag("fluid");
+const FlagUID NoSlip_Flag("no slip");
+
+/////////////////////////////////////
+// BOUNDARY HANDLING CUSTOMIZATION //
+/////////////////////////////////////
+template< typename ParticleAccessor_T >
+class MyBoundaryHandling
+{
+ public:
+ using NoSlip_T = lbm::NoSlip< LatticeModel_T, flag_t >;
+ using Type = BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T >;
+
+ MyBoundaryHandling(const BlockDataID& flagFieldID, const BlockDataID& pdfFieldID,
+ const shared_ptr< ParticleAccessor_T >& ac)
+ : flagFieldID_(flagFieldID), pdfFieldID_(pdfFieldID), ac_(ac)
+ {}
+
+ Type* operator()(IBlock* const block, const StructuredBlockStorage* const /*storage*/) const
+ {
+ WALBERLA_ASSERT_NOT_NULLPTR(block);
+
+ auto* flagField = block->getData< FlagField_T >(flagFieldID_);
+ auto* pdfField = block->getData< PdfField_T >(pdfFieldID_);
+
+ const auto fluid =
+ flagField->flagExists(Fluid_Flag) ? flagField->getFlag(Fluid_Flag) : flagField->registerFlag(Fluid_Flag);
+
+ Type* handling =
+ new Type("moving obstacle boundary handling", flagField, fluid, NoSlip_T("NoSlip", NoSlip_Flag, pdfField));
+
+ handling->fillWithDomain(FieldGhostLayers);
+
+ return handling;
+ }
+
+ private:
+ const BlockDataID flagFieldID_;
+ const BlockDataID pdfFieldID_;
+
+ shared_ptr< ParticleAccessor_T > ac_;
+};
+//*******************************************************************************************************************
+
+//*******************************************************************************************************************
+/*!\brief Evaluating the position and velocity of the sphere
+ *
+ */
+//*******************************************************************************************************************
+template< typename ParticleAccessor_T >
+class SpherePropertyLogger
+{
+ public:
+ SpherePropertyLogger(const shared_ptr< ParticleAccessor_T >& ac, walberla::id_t sphereUid,
+ const std::string& fileName, bool fileIO, real_t dx_SI, real_t dt_SI, real_t diameter,
+ real_t gravitationalForceMag)
+ : ac_(ac), sphereUid_(sphereUid), fileName_(fileName), fileIO_(fileIO), dx_SI_(dx_SI), dt_SI_(dt_SI),
+ diameter_(diameter), gravitationalForceMag_(gravitationalForceMag), position_(real_t(0)),
+ maxVelocity_(real_t(0))
+ {
+ if (fileIO_)
+ {
+ WALBERLA_ROOT_SECTION()
+ {
+ std::ofstream file;
+ file.open(fileName_.c_str());
+ file << "#\t t\t posX\t posY\t gapZ\t velX\t velY\t velZ\n";
+ file.close();
+ }
+ }
+ }
+
+ void operator()(const uint_t timestep)
+ {
+ Vector3< real_t > pos(real_t(0));
+ Vector3< real_t > transVel(real_t(0));
+ Vector3< real_t > hydForce(real_t(0));
+
+ size_t idx = ac_->uidToIdx(sphereUid_);
+ if (idx != ac_->getInvalidIdx())
+ {
+ if (!mesa_pd::data::particle_flags::isSet(ac_->getFlags(idx), mesa_pd::data::particle_flags::GHOST))
+ {
+ pos = ac_->getPosition(idx);
+ transVel = ac_->getLinearVelocity(idx);
+ hydForce = ac_->getHydrodynamicForce(idx);
+ }
+ }
+
+ WALBERLA_MPI_SECTION()
+ {
+ mpi::allReduceInplace(pos[0], mpi::SUM);
+ mpi::allReduceInplace(pos[1], mpi::SUM);
+ mpi::allReduceInplace(pos[2], mpi::SUM);
+
+ mpi::allReduceInplace(transVel[0], mpi::SUM);
+ mpi::allReduceInplace(transVel[1], mpi::SUM);
+ mpi::allReduceInplace(transVel[2], mpi::SUM);
+
+ mpi::allReduceInplace(hydForce[0], mpi::SUM);
+ mpi::allReduceInplace(hydForce[1], mpi::SUM);
+ mpi::allReduceInplace(hydForce[2], mpi::SUM);
+ }
+
+ position_ = pos[2];
+ maxVelocity_ = std::max(maxVelocity_, -transVel[2]);
+
+ if (fileIO_) writeToFile(timestep, pos, transVel, hydForce);
+ }
+
+ real_t getPosition() const { return position_; }
+
+ real_t getMaxVelocity() const { return maxVelocity_; }
+
+ private:
+ void writeToFile(const uint_t timestep, const Vector3< real_t >& position, const Vector3< real_t >& velocity,
+ const Vector3< real_t >& hydForce)
+ {
+ WALBERLA_ROOT_SECTION()
+ {
+ std::ofstream file;
+ file.open(fileName_.c_str(), std::ofstream::app);
+
+ auto scaledPosition = position / diameter_;
+ auto velocity_SI = velocity * dx_SI_ / dt_SI_;
+ auto normalizedHydForce = hydForce / gravitationalForceMag_;
+
+ file << timestep << "\t" << real_c(timestep) * dt_SI_ << "\t"
+ << "\t" << scaledPosition[0] << "\t" << scaledPosition[1] << "\t" << scaledPosition[2] - real_t(0.5)
+ << "\t" << velocity_SI[0] << "\t" << velocity_SI[1] << "\t" << velocity_SI[2] << "\t"
+ << normalizedHydForce[0] << "\t" << normalizedHydForce[1] << "\t" << normalizedHydForce[2] << "\n";
+ file.close();
+ }
+ }
+
+ shared_ptr< ParticleAccessor_T > ac_;
+ const walberla::id_t sphereUid_;
+ std::string fileName_;
+ bool fileIO_;
+ real_t dx_SI_, dt_SI_, diameter_, gravitationalForceMag_;
+
+ real_t position_;
+ real_t maxVelocity_;
+};
+
+void createPlaneSetup(const shared_ptr< mesa_pd::data::ParticleStorage >& ps,
+ const shared_ptr< mesa_pd::data::ShapeStorage >& ss, const math::AABB& simulationDomain)
+{
+ // create bounding planes
+ mesa_pd::data::Particle p0 = *ps->create(true);
+ p0.setPosition(simulationDomain.minCorner());
+ p0.setInteractionRadius(std::numeric_limits< real_t >::infinity());
+ p0.setShapeID(ss->create< mesa_pd::data::HalfSpace >(Vector3< real_t >(0, 0, 1)));
+ p0.setOwner(mpi::MPIManager::instance()->rank());
+ p0.setType(0);
+ mesa_pd::data::particle_flags::set(p0.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+ mesa_pd::data::particle_flags::set(p0.getFlagsRef(), mesa_pd::data::particle_flags::FIXED);
+
+ mesa_pd::data::Particle p1 = *ps->create(true);
+ p1.setPosition(simulationDomain.maxCorner());
+ p1.setInteractionRadius(std::numeric_limits< real_t >::infinity());
+ p1.setShapeID(ss->create< mesa_pd::data::HalfSpace >(Vector3< real_t >(0, 0, -1)));
+ p1.setOwner(mpi::MPIManager::instance()->rank());
+ p1.setType(0);
+ mesa_pd::data::particle_flags::set(p1.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+ mesa_pd::data::particle_flags::set(p1.getFlagsRef(), mesa_pd::data::particle_flags::FIXED);
+
+ mesa_pd::data::Particle p2 = *ps->create(true);
+ p2.setPosition(simulationDomain.minCorner());
+ p2.setInteractionRadius(std::numeric_limits< real_t >::infinity());
+ p2.setShapeID(ss->create< mesa_pd::data::HalfSpace >(Vector3< real_t >(1, 0, 0)));
+ p2.setOwner(mpi::MPIManager::instance()->rank());
+ p2.setType(0);
+ mesa_pd::data::particle_flags::set(p2.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+ mesa_pd::data::particle_flags::set(p2.getFlagsRef(), mesa_pd::data::particle_flags::FIXED);
+
+ mesa_pd::data::Particle p3 = *ps->create(true);
+ p3.setPosition(simulationDomain.maxCorner());
+ p3.setInteractionRadius(std::numeric_limits< real_t >::infinity());
+ p3.setShapeID(ss->create< mesa_pd::data::HalfSpace >(Vector3< real_t >(-1, 0, 0)));
+ p3.setOwner(mpi::MPIManager::instance()->rank());
+ p3.setType(0);
+ mesa_pd::data::particle_flags::set(p3.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+ mesa_pd::data::particle_flags::set(p3.getFlagsRef(), mesa_pd::data::particle_flags::FIXED);
+
+ mesa_pd::data::Particle p4 = *ps->create(true);
+ p4.setPosition(simulationDomain.minCorner());
+ p4.setInteractionRadius(std::numeric_limits< real_t >::infinity());
+ p4.setShapeID(ss->create< mesa_pd::data::HalfSpace >(Vector3< real_t >(0, 1, 0)));
+ p4.setOwner(mpi::MPIManager::instance()->rank());
+ p4.setType(0);
+ mesa_pd::data::particle_flags::set(p4.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+ mesa_pd::data::particle_flags::set(p4.getFlagsRef(), mesa_pd::data::particle_flags::FIXED);
+
+ mesa_pd::data::Particle p5 = *ps->create(true);
+ p5.setPosition(simulationDomain.maxCorner());
+ p5.setInteractionRadius(std::numeric_limits< real_t >::infinity());
+ p5.setShapeID(ss->create< mesa_pd::data::HalfSpace >(Vector3< real_t >(0, -1, 0)));
+ p5.setOwner(mpi::MPIManager::instance()->rank());
+ p5.setType(0);
+ mesa_pd::data::particle_flags::set(p5.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+ mesa_pd::data::particle_flags::set(p5.getFlagsRef(), mesa_pd::data::particle_flags::FIXED);
+}
+
+//////////
+// MAIN //
+//////////
+
+//*******************************************************************************************************************
+/*!\brief Testcase that simulates the settling of a sphere inside a rectangular column filled with viscous fluid
+ *
+ * see: ten Cate, Nieuwstad, Derksen, Van den Akker - "Particle imaging velocimetry experiments and lattice-Boltzmann
+ * simulations on a single sphere settling under gravity" (2002), Physics of Fluids, doi: 10.1063/1.1512918
+ */
+//*******************************************************************************************************************
+
+int main(int argc, char** argv)
+{
+ debug::enterTestMode();
+
+ mpi::Environment env(argc, argv);
+
+ ///////////////////
+ // Customization //
+ ///////////////////
+
+ // simulation control
+ bool shortrun = false;
+ bool funcTest = false;
+ bool fileIO = false;
+ uint_t vtkIOFreq = 0;
+ std::string baseFolder = "vtk_out_SettlingSphere";
+
+ // physical setup
+ uint_t fluidType = 1;
+
+ // numerical parameters
+ uint_t numberOfCellsInHorizontalDirection = uint_t(135);
+ bool averageForceTorqueOverTwoTimeSteps = true;
+ uint_t numRPDSubCycles = uint_t(1);
+ bool useVelocityVerlet = false;
+
+ for (int i = 1; i < argc; ++i)
+ {
+ if (std::strcmp(argv[i], "--shortrun") == 0)
+ {
+ shortrun = true;
+ continue;
+ }
+ if (std::strcmp(argv[i], "--funcTest") == 0)
+ {
+ funcTest = true;
+ continue;
+ }
+ if (std::strcmp(argv[i], "--fileIO") == 0)
+ {
+ fileIO = true;
+ continue;
+ }
+ if (std::strcmp(argv[i], "--vtkIOFreq") == 0)
+ {
+ vtkIOFreq = uint_c(std::atof(argv[++i]));
+ continue;
+ }
+ if (std::strcmp(argv[i], "--fluidType") == 0)
+ {
+ fluidType = uint_c(std::atof(argv[++i]));
+ continue;
+ }
+ if (std::strcmp(argv[i], "--numRPDSubCycles") == 0)
+ {
+ numRPDSubCycles = uint_c(std::atof(argv[++i]));
+ continue;
+ }
+ if (std::strcmp(argv[i], "--resolution") == 0)
+ {
+ numberOfCellsInHorizontalDirection = uint_c(std::atof(argv[++i]));
+ continue;
+ }
+ if (std::strcmp(argv[i], "--noForceAveraging") == 0)
+ {
+ averageForceTorqueOverTwoTimeSteps = false;
+ continue;
+ }
+ if (std::strcmp(argv[i], "--baseFolder") == 0)
+ {
+ baseFolder = argv[++i];
+ continue;
+ }
+ if (std::strcmp(argv[i], "--useVV") == 0)
+ {
+ useVelocityVerlet = true;
+ continue;
+ }
+ WALBERLA_ABORT("Unrecognized command line argument found: " << argv[i]);
+ }
+
+ if (funcTest) { walberla::logging::Logging::instance()->setLogLevel(logging::Logging::LogLevel::WARNING); }
+
+ if (fileIO)
+ {
+ // create base directory if it does not yet exist
+ filesystem::path tpath(baseFolder);
+ if (!filesystem::exists(tpath)) filesystem::create_directory(tpath);
+ }
+
+ //////////////////////////////////////
+ // SIMULATION PROPERTIES in SI units//
+ //////////////////////////////////////
+
+ // values are mainly taken from the reference paper
+ const real_t diameter_SI = real_t(15e-3);
+ const real_t densitySphere_SI = real_t(1120);
+
+ real_t densityFluid_SI, dynamicViscosityFluid_SI;
+ real_t expectedSettlingVelocity_SI;
+ switch (fluidType)
+ {
+ case 1:
+ // Re_p around 1.5
+ densityFluid_SI = real_t(970);
+ dynamicViscosityFluid_SI = real_t(373e-3);
+ expectedSettlingVelocity_SI = real_t(0.035986);
+ break;
+ case 2:
+ // Re_p around 4.1
+ densityFluid_SI = real_t(965);
+ dynamicViscosityFluid_SI = real_t(212e-3);
+ expectedSettlingVelocity_SI = real_t(0.05718);
+ break;
+ case 3:
+ // Re_p around 11.6
+ densityFluid_SI = real_t(962);
+ dynamicViscosityFluid_SI = real_t(113e-3);
+ expectedSettlingVelocity_SI = real_t(0.087269);
+ break;
+ case 4:
+ // Re_p around 31.9
+ densityFluid_SI = real_t(960);
+ dynamicViscosityFluid_SI = real_t(58e-3);
+ expectedSettlingVelocity_SI = real_t(0.12224);
+ break;
+ default:
+ WALBERLA_ABORT("Only four different fluids are supported! Choose type between 1 and 4.");
+ }
+ const real_t kinematicViscosityFluid_SI = dynamicViscosityFluid_SI / densityFluid_SI;
+
+ const real_t gravitationalAcceleration_SI = real_t(9.81);
+ Vector3< real_t > domainSize_SI(real_t(100e-3), real_t(100e-3), real_t(160e-3));
+ // shift starting gap a bit upwards to match the reported (plotted) values
+ const real_t startingGapSize_SI = real_t(120e-3) + real_t(0.25) * diameter_SI;
+
+ WALBERLA_LOG_INFO_ON_ROOT("Setup (in SI units):");
+ WALBERLA_LOG_INFO_ON_ROOT(" - fluid type = " << fluidType);
+ WALBERLA_LOG_INFO_ON_ROOT(" - domain size = " << domainSize_SI);
+ WALBERLA_LOG_INFO_ON_ROOT(" - sphere: diameter = " << diameter_SI << ", density = " << densitySphere_SI
+ << ", starting gap size = " << startingGapSize_SI);
+ WALBERLA_LOG_INFO_ON_ROOT(" - fluid: density = " << densityFluid_SI << ", dyn. visc = " << dynamicViscosityFluid_SI
+ << ", kin. visc = " << kinematicViscosityFluid_SI);
+ WALBERLA_LOG_INFO_ON_ROOT(" - expected settling velocity = "
+ << expectedSettlingVelocity_SI << " --> Re_p = "
+ << expectedSettlingVelocity_SI * diameter_SI / kinematicViscosityFluid_SI);
+
+ //////////////////////////
+ // NUMERICAL PARAMETERS //
+ //////////////////////////
+
+ const real_t dx_SI = domainSize_SI[0] / real_c(numberOfCellsInHorizontalDirection);
+ const Vector3< uint_t > domainSize(uint_c(floor(domainSize_SI[0] / dx_SI + real_t(0.5))),
+ uint_c(floor(domainSize_SI[1] / dx_SI + real_t(0.5))),
+ uint_c(floor(domainSize_SI[2] / dx_SI + real_t(0.5))));
+ const real_t diameter = diameter_SI / dx_SI;
+ const real_t sphereVolume = math::pi / real_t(6) * diameter * diameter * diameter;
+
+ const real_t expectedSettlingVelocity = real_t(0.01);
+ const real_t dt_SI = expectedSettlingVelocity / expectedSettlingVelocity_SI * dx_SI;
+
+ const real_t viscosity = kinematicViscosityFluid_SI * dt_SI / (dx_SI * dx_SI);
+ const real_t relaxationTime = real_t(1) / lbm::collision_model::omegaFromViscosity(viscosity);
+
+ const real_t gravitationalAcceleration = gravitationalAcceleration_SI * dt_SI * dt_SI / dx_SI;
+
+ const real_t densityFluid = real_t(1);
+ const real_t densitySphere = densityFluid * densitySphere_SI / densityFluid_SI;
+
+ const real_t dx = real_t(1);
+
+ const uint_t timesteps = funcTest ? 1 : (shortrun ? uint_t(200) : uint_t(250000));
+
+ WALBERLA_LOG_INFO_ON_ROOT(" - dx_SI = " << dx_SI << ", dt_SI = " << dt_SI);
+ WALBERLA_LOG_INFO_ON_ROOT("Setup (in simulation, i.e. lattice, units):");
+ WALBERLA_LOG_INFO_ON_ROOT(" - domain size = " << domainSize);
+ WALBERLA_LOG_INFO_ON_ROOT(" - sphere: diameter = " << diameter << ", density = " << densitySphere);
+ WALBERLA_LOG_INFO_ON_ROOT(" - fluid: density = " << densityFluid << ", relaxation time (tau) = " << relaxationTime
+ << ", kin. visc = " << viscosity);
+ WALBERLA_LOG_INFO_ON_ROOT(" - gravitational acceleration = " << gravitationalAcceleration);
+ WALBERLA_LOG_INFO_ON_ROOT(" - expected settling velocity = " << expectedSettlingVelocity << " --> Re_p = "
+ << expectedSettlingVelocity * diameter / viscosity);
+ WALBERLA_LOG_INFO_ON_ROOT(" - integrator = " << (useVelocityVerlet ? "Velocity Verlet" : "Explicit Euler"));
+
+ if (vtkIOFreq > 0)
+ {
+ WALBERLA_LOG_INFO_ON_ROOT(" - writing vtk files to folder \"" << baseFolder << "\" with frequency " << vtkIOFreq);
+ }
+
+ ///////////////////////////
+ // BLOCK STRUCTURE SETUP //
+ ///////////////////////////
+
+ Vector3< uint_t > numberOfBlocksPerDirection(uint_t(1), uint_t(1), uint_t(4));
+ Vector3< uint_t > cellsPerBlockPerDirection(domainSize[0] / numberOfBlocksPerDirection[0],
+ domainSize[1] / numberOfBlocksPerDirection[1],
+ domainSize[2] / numberOfBlocksPerDirection[2]);
+
+ for (uint_t i = 0; i < 3; ++i)
+ {
+ WALBERLA_CHECK_EQUAL(cellsPerBlockPerDirection[i] * numberOfBlocksPerDirection[i], domainSize[i],
+ "Unmatching domain decomposition in direction " << i << "!");
+ }
+
+ auto blocks = blockforest::createUniformBlockGrid(numberOfBlocksPerDirection[0], numberOfBlocksPerDirection[1],
+ numberOfBlocksPerDirection[2], cellsPerBlockPerDirection[0],
+ cellsPerBlockPerDirection[1], cellsPerBlockPerDirection[2], dx, 0,
+ false, false, false, false, false, // periodicity
+ false);
+
+ WALBERLA_LOG_INFO_ON_ROOT("Domain decomposition:");
+ WALBERLA_LOG_INFO_ON_ROOT(" - blocks per direction = " << numberOfBlocksPerDirection);
+ WALBERLA_LOG_INFO_ON_ROOT(" - cells per block = " << cellsPerBlockPerDirection);
+
+ // write domain decomposition to file
+ if (vtkIOFreq > 0) { vtk::writeDomainDecomposition(blocks, "initial_domain_decomposition", baseFolder); }
+
+ //////////////////
+ // RPD COUPLING //
+ //////////////////
+
+ auto rpdDomain = std::make_shared< mesa_pd::domain::BlockForestDomain >(blocks->getBlockForestPointer());
+
+ // init data structures
+ auto ps = walberla::make_shared< mesa_pd::data::ParticleStorage >(1);
+ auto ss = walberla::make_shared< mesa_pd::data::ShapeStorage >();
+ using ParticleAccessor_T = mesa_pd::data::ParticleAccessorWithShape;
+ auto accessor = walberla::make_shared< ParticleAccessor_T >(ps, ss);
+
+ // bounding planes
+ createPlaneSetup(ps, ss, blocks->getDomain());
+
+ // create sphere and store Uid
+ Vector3< real_t > initialPosition(real_t(0.5) * real_c(domainSize[0]), real_t(0.5) * real_c(domainSize[1]),
+ startingGapSize_SI / dx_SI + real_t(0.5) * diameter);
+ auto sphereShape = ss->create< mesa_pd::data::Sphere >(diameter * real_t(0.5));
+ ss->shapes[sphereShape]->updateMassAndInertia(densitySphere);
+
+ walberla::id_t sphereUid = 0;
+ if (rpdDomain->isContainedInProcessSubdomain(uint_c(mpi::MPIManager::instance()->rank()), initialPosition))
+ {
+ mesa_pd::data::Particle&& p = *ps->create();
+ p.setPosition(initialPosition);
+ p.setInteractionRadius(diameter * real_t(0.5));
+ p.setOwner(mpi::MPIManager::instance()->rank());
+ p.setShapeID(sphereShape);
+ sphereUid = p.getUid();
+ }
+ mpi::allReduceInplace(sphereUid, mpi::SUM);
+
+ ///////////////////////
+ // ADD DATA TO BLOCKS //
+ ////////////////////////
+
+ // create the lattice model
+ LatticeModel_T latticeModel =
+ LatticeModel_T(lbm::collision_model::TRT::constructWithMagicNumber(real_t(1) / relaxationTime));
+
+ // add PDF field
+ BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >(
+ blocks, "pdf field (zyxf)", latticeModel, Vector3< real_t >(real_t(0)), real_t(1), uint_t(1), field::zyxf);
+ // add flag field
+ BlockDataID flagFieldID = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
+
+ // add boundary handling
+ using BoundaryHandling_T = MyBoundaryHandling< ParticleAccessor_T >::Type;
+ BlockDataID boundaryHandlingID = blocks->addStructuredBlockData< BoundaryHandling_T >(
+ MyBoundaryHandling< ParticleAccessor_T >(flagFieldID, pdfFieldID, accessor), "boundary handling");
+
+ // set up RPD functionality
+ std::function< void(void) > syncCall = [ps, rpdDomain]() {
+ const real_t overlap = real_t(1.5);
+ mesa_pd::mpi::SyncNextNeighbors syncNextNeighborFunc;
+ syncNextNeighborFunc(*ps, *rpdDomain, overlap);
+ };
+
+ syncCall();
+
+ mesa_pd::kernel::ExplicitEuler explEulerIntegrator(real_t(1) / real_t(numRPDSubCycles));
+ mesa_pd::kernel::VelocityVerletPreForceUpdate vvIntegratorPreForce(real_t(1) / real_t(numRPDSubCycles));
+ mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(real_t(1) / real_t(numRPDSubCycles));
+
+ mesa_pd::kernel::SpringDashpot collisionResponse(1);
+ mesa_pd::mpi::ReduceProperty reduceProperty;
+
+ // set up coupling functionality
+ lbm_mesapd_coupling::RegularParticlesSelector sphereSelector;
+ Vector3< real_t > gravitationalForce(real_t(0), real_t(0),
+ -(densitySphere - densityFluid) * gravitationalAcceleration * sphereVolume);
+ lbm_mesapd_coupling::AddForceOnParticlesKernel addGravitationalForce(gravitationalForce);
+ lbm_mesapd_coupling::AddHydrodynamicInteractionKernel addHydrodynamicInteraction;
+ lbm_mesapd_coupling::ResetHydrodynamicForceTorqueKernel resetHydrodynamicForceTorque;
+ lbm_mesapd_coupling::AverageHydrodynamicForceTorqueKernel averageHydrodynamicForceTorque;
+ lbm_mesapd_coupling::LubricationCorrectionKernel lubricationCorrectionKernel(
+ viscosity, [](real_t r) { return real_t(0.0016) * r; });
+ lbm_mesapd_coupling::ParticleMappingKernel< BoundaryHandling_T > particleMappingKernel(blocks, boundaryHandlingID);
+
+ ///////////////
+ // TIME LOOP //
+ ///////////////
+
+ // map planes into the LBM simulation -> act as no-slip boundaries
+ ps->forEachParticle(false, lbm_mesapd_coupling::GlobalParticlesSelector(), *accessor, particleMappingKernel,
+ *accessor, NoSlip_Flag);
+
+ // map particles into the LBM simulation
+ BlockDataID particleAndVolumeFractionFieldID =
+ field::addToStorage< lbm_mesapd_coupling::psm::ParticleAndVolumeFractionField_T >(
+ blocks, "particle and volume fraction field",
+ std::vector< lbm_mesapd_coupling::psm::ParticleAndVolumeFraction_T >(), field::zyxf, 0);
+ lbm_mesapd_coupling::psm::ParticleAndVolumeFractionMapping particleMapping(blocks, accessor, sphereSelector,
+ particleAndVolumeFractionFieldID, 4);
+ particleMapping();
+
+ lbm_mesapd_coupling::psm::initializeDomainForPSM< LatticeModel_T, 1 >(*blocks, pdfFieldID,
+ particleAndVolumeFractionFieldID, *accessor);
+
+ walberla::lbm_mesapd_coupling::FractionFieldSum fractionFieldSum(blocks, particleAndVolumeFractionFieldID);
+ // check that the sum of all fractions is roughly the volume of the particle
+ real_t particleRadius = diameter * real_t(0.5);
+ WALBERLA_CHECK_LESS(
+ std::fabs(4.0 / 3.0 * math::pi * particleRadius * particleRadius * particleRadius - fractionFieldSum()),
+ real_c(1.0));
+
+ // setup of the LBM communication for synchronizing the pdf field between neighboring blocks
+ std::function< void() > commFunction;
+ blockforest::communication::UniformBufferedScheme< Stencil_T > scheme(blocks);
+ scheme.addPackInfo(make_shared< lbm::PdfFieldPackInfo< LatticeModel_T > >(pdfFieldID));
+ commFunction = scheme;
+
+ // create the timeloop
+ SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);
+
+ auto bhSweep = BoundaryHandling_T::getBlockSweep(boundaryHandlingID);
+ auto lbmSweep = lbm_mesapd_coupling::psm::makePSMSweep< LatticeModel_T, FlagField_T, 1, 1 >(
+ pdfFieldID, particleAndVolumeFractionFieldID, blocks, accessor, flagFieldID, Fluid_Flag);
+
+ timeloop.addFuncBeforeTimeStep(RemainingTimeLogger(timeloop.getNrOfTimeSteps()), "Remaining Time Logger");
+
+ // vtk output
+ if (vtkIOFreq != uint_t(0))
+ {
+ // spheres
+ auto particleVtkOutput = make_shared< mesa_pd::vtk::ParticleVtkOutput >(ps);
+ particleVtkOutput->addOutput< mesa_pd::data::SelectParticleOwner >("owner");
+ particleVtkOutput->addOutput< mesa_pd::data::SelectParticleLinearVelocity >("velocity");
+ auto particleVtkWriter =
+ vtk::createVTKOutput_PointData(particleVtkOutput, "Particles", vtkIOFreq, baseFolder, "simulation_step");
+ timeloop.addFuncBeforeTimeStep(vtk::writeFiles(particleVtkWriter), "VTK (sphere data)");
+
+ // flag field (written only once in the first time step, ghost layers are also written)
+ // auto flagFieldVTK = vtk::createVTKOutput_BlockData( blocks, "flag_field", timesteps, FieldGhostLayers, false,
+ // baseFolder ); flagFieldVTK->addCellDataWriter( make_shared< field::VTKWriter< FlagField_T > >( flagFieldID,
+ // "FlagField" ) ); timeloop.addFuncBeforeTimeStep( vtk::writeFiles( flagFieldVTK ), "VTK (flag field data)" );
+
+ // pdf field
+ auto pdfFieldVTK = vtk::createVTKOutput_BlockData(blocks, "fluid_field", vtkIOFreq, 0, false, baseFolder);
+
+ blockforest::communication::UniformBufferedScheme< stencil::D3Q27 > pdfGhostLayerSync(blocks);
+ pdfGhostLayerSync.addPackInfo(make_shared< field::communication::PackInfo< PdfField_T > >(pdfFieldID));
+ pdfFieldVTK->addBeforeFunction(pdfGhostLayerSync);
+
+ field::FlagFieldCellFilter< FlagField_T > fluidFilter(flagFieldID);
+ fluidFilter.addFlag(Fluid_Flag);
+ pdfFieldVTK->addCellInclusionFilter(fluidFilter);
+
+ pdfFieldVTK->addCellDataWriter(
+ make_shared< lbm::VelocityVTKWriter< LatticeModel_T, float > >(pdfFieldID, "VelocityFromPDF"));
+ pdfFieldVTK->addCellDataWriter(
+ make_shared< lbm::DensityVTKWriter< LatticeModel_T, float > >(pdfFieldID, "DensityFromPDF"));
+
+ timeloop.addFuncBeforeTimeStep(vtk::writeFiles(pdfFieldVTK), "VTK (fluid field data)");
+ }
+
+ bool useStreamCollide = true;
+
+ if (useStreamCollide)
+ {
+ // add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip
+ // treatment)
+ timeloop.add() << BeforeFunction(commFunction, "LBM Communication") << Sweep(bhSweep, "Boundary Handling");
+
+ // stream + collide LBM step
+ timeloop.add() << Sweep(makeSharedSweep(lbmSweep), "cell-wise LB sweep");
+ }
+ else
+ {
+ // collide LBM step
+ timeloop.add() << Sweep(makeCollideSweep(lbmSweep), "cell-wise collide LB sweep");
+
+ // add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip
+ // treatment)
+ timeloop.add() << BeforeFunction(commFunction, "LBM Communication")
+ << Sweep(BoundaryHandling_T::getBlockSweep(boundaryHandlingID), "Boundary Handling");
+
+ // stream LBM step
+ timeloop.add() << Sweep(makeStreamSweep(lbmSweep), "cell-wise stream LB sweep");
+ }
+
+ // evaluation functionality
+ std::string loggingFileName(baseFolder + "/LoggingSettlingSphere_");
+ loggingFileName += std::to_string(fluidType);
+ loggingFileName += ".txt";
+ if (fileIO) { WALBERLA_LOG_INFO_ON_ROOT(" - writing logging output to file \"" << loggingFileName << "\""); }
+ SpherePropertyLogger< ParticleAccessor_T > logger(accessor, sphereUid, loggingFileName, fileIO, dx_SI, dt_SI,
+ diameter, -gravitationalForce[2]);
+
+ ////////////////////////
+ // EXECUTE SIMULATION //
+ ////////////////////////
+
+ WcTimingPool timeloopTiming;
+
+ real_t terminationPosition = real_t(0.51) * diameter; // right before sphere touches the bottom wall
+
+ const bool useOpenMP = false;
+
+ // time loop
+ for (uint_t i = 0; i < timesteps; ++i)
+ {
+ // perform a single simulation step -> this contains LBM and setting of the hydrodynamic interactions
+ timeloop.singleStep(timeloopTiming);
+
+ timeloopTiming["RPD"].start();
+
+ // -> full hydrodynamic force/torque info is available on local particle
+ reduceProperty.operator()< mesa_pd::HydrodynamicForceTorqueNotification >(*ps);
+
+ if (averageForceTorqueOverTwoTimeSteps)
+ {
+ if (i == 0)
+ {
+ lbm_mesapd_coupling::InitializeHydrodynamicForceTorqueForAveragingKernel
+ initializeHydrodynamicForceTorqueForAveragingKernel;
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor,
+ initializeHydrodynamicForceTorqueForAveragingKernel, *accessor);
+ }
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, averageHydrodynamicForceTorque,
+ *accessor);
+ }
+
+ for (auto subCycle = uint_t(0); subCycle < numRPDSubCycles; ++subCycle)
+ {
+ if (useVelocityVerlet)
+ {
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, vvIntegratorPreForce, *accessor);
+ syncCall();
+ }
+
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addHydrodynamicInteraction,
+ *accessor);
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addGravitationalForce, *accessor);
+
+ // lubrication correction
+ ps->forEachParticlePairHalf(
+ useOpenMP, mesa_pd::kernel::ExcludeInfiniteInfinite(), *accessor,
+ [&lubricationCorrectionKernel, rpdDomain](const size_t idx1, const size_t idx2, auto& ac) {
+ // TODO change this to storing copy, not reference
+ mesa_pd::collision_detection::AnalyticContactDetection acd;
+ acd.getContactThreshold() = lubricationCorrectionKernel.getNormalCutOffDistance();
+ mesa_pd::kernel::DoubleCast double_cast;
+ mesa_pd::mpi::ContactFilter contact_filter;
+ if (double_cast(idx1, idx2, ac, acd, ac))
+ {
+ if (contact_filter(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), *rpdDomain))
+ {
+ double_cast(idx1, idx2, ac, lubricationCorrectionKernel, ac, acd.getContactNormal(),
+ acd.getPenetrationDepth());
+ }
+ }
+ },
+ *accessor);
+
+ // one could add linked cells here
+
+ // collision response
+ ps->forEachParticlePairHalf(
+ useOpenMP, mesa_pd::kernel::ExcludeInfiniteInfinite(), *accessor,
+ [collisionResponse, rpdDomain](const size_t idx1, const size_t idx2, auto& ac) {
+ mesa_pd::collision_detection::AnalyticContactDetection acd;
+ mesa_pd::kernel::DoubleCast double_cast;
+ mesa_pd::mpi::ContactFilter contact_filter;
+ if (double_cast(idx1, idx2, ac, acd, ac))
+ {
+ if (contact_filter(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), *rpdDomain))
+ {
+ collisionResponse(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), acd.getContactNormal(),
+ acd.getPenetrationDepth());
+ }
+ }
+ },
+ *accessor);
+
+ reduceProperty.operator()< mesa_pd::ForceTorqueNotification >(*ps);
+
+ if (useVelocityVerlet)
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, vvIntegratorPostForce, *accessor);
+ else
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, explEulerIntegrator, *accessor);
+
+ syncCall();
+ particleMapping();
+ WALBERLA_CHECK_LESS(
+ std::fabs(4.0 / 3.0 * math::pi * particleRadius * particleRadius * particleRadius - fractionFieldSum()),
+ real_c(1.0));
+ }
+
+ timeloopTiming["RPD"].end();
+
+ // evaluation
+ timeloopTiming["Logging"].start();
+ logger(i);
+ timeloopTiming["Logging"].end();
+
+ // reset after logging here
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, resetHydrodynamicForceTorque, *accessor);
+
+ // check for termination
+ if (logger.getPosition() < terminationPosition)
+ {
+ WALBERLA_LOG_INFO_ON_ROOT("Sphere reached terminal position " << logger.getPosition() << " after " << i
+ << " timesteps!");
+ break;
+ }
+ }
+
+ timeloopTiming.logResultOnRoot();
+
+ // check the result
+ if (!funcTest && !shortrun)
+ {
+ real_t relErr = std::fabs(expectedSettlingVelocity - logger.getMaxVelocity()) / expectedSettlingVelocity;
+ WALBERLA_LOG_INFO_ON_ROOT("Expected maximum settling velocity: " << expectedSettlingVelocity);
+ WALBERLA_LOG_INFO_ON_ROOT("Simulated maximum settling velocity: " << logger.getMaxVelocity());
+ WALBERLA_LOG_INFO_ON_ROOT("Relative error: " << relErr);
+
+ // the relative error has to be below 10%
+ WALBERLA_CHECK_LESS(relErr, real_t(0.1));
+ }
+
+ return EXIT_SUCCESS;
+}
+
+} // namespace settling_sphere
+
+int main(int argc, char** argv) { settling_sphere::main(argc, argv); }
diff --git a/tests/lbm_mesapd_coupling/partially_saturated_cells_method/TorqueSphere.cpp b/tests/lbm_mesapd_coupling/partially_saturated_cells_method/TorqueSphere.cpp
new file mode 100644
index 000000000..9b965c164
--- /dev/null
+++ b/tests/lbm_mesapd_coupling/partially_saturated_cells_method/TorqueSphere.cpp
@@ -0,0 +1,561 @@
+//======================================================================================================================
+//
+// 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 TorqueSphere.cpp
+//! \ingroup lbm_mesapd_coupling
+//! \author Samuel Kemmler <samuel.kemmler@fau.de>
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \brief Modification of pe_coupling/partially_saturated_cells_method/TorqueSpherePSM.cpp
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "boundary/all.h"
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/SharedFunctor.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/logging/Logging.h"
+#include "core/mpi/MPIManager.h"
+#include "core/mpi/Reduce.h"
+#include "core/timing/RemainingTimeLogger.h"
+
+#include "domain_decomposition/SharedSweep.h"
+
+#include "field/AddToStorage.h"
+#include "field/communication/PackInfo.h"
+
+#include "lbm/field/AddToStorage.h"
+#include "lbm/field/PdfField.h"
+#include "lbm/lattice_model/D3Q19.h"
+#include "lbm/lattice_model/ForceModel.h"
+
+#include "lbm_mesapd_coupling/DataTypes.h"
+#include "lbm_mesapd_coupling/partially_saturated_cells_method/PSMSweep.h"
+#include "lbm_mesapd_coupling/partially_saturated_cells_method/PSMUtility.h"
+#include "lbm_mesapd_coupling/partially_saturated_cells_method/ParticleAndVolumeFractionMapping.h"
+#include "lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
+
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/ShapeStorage.h"
+#include "mesa_pd/domain/BlockForestDomain.h"
+#include "mesa_pd/mpi/SyncNextNeighbors.h"
+
+#include "stencil/D3Q27.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+#include <iostream>
+#include <vector>
+
+namespace torque_sphere_psm
+{
+
+///////////
+// USING //
+///////////
+
+using namespace walberla;
+using walberla::uint_t;
+
+// PDF field, flag field & particle field
+typedef lbm::D3Q19< lbm::collision_model::SRT, false > LatticeModel_T;
+
+using Stencil_T = LatticeModel_T::Stencil;
+using PdfField_T = lbm::PdfField< LatticeModel_T >;
+
+using flag_t = walberla::uint8_t;
+using FlagField_T = FlagField< flag_t >;
+
+///////////
+// FLAGS //
+///////////
+
+const FlagUID Fluid_Flag("fluid");
+
+////////////////
+// PARAMETERS //
+////////////////
+
+struct Setup
+{
+ uint_t checkFrequency;
+ real_t visc;
+ real_t tau;
+ real_t radius;
+ uint_t length;
+ real_t phi;
+ real_t angularVel;
+ real_t analyticalTorque;
+};
+
+//*******************************************************************************************************************
+/*!\brief Evaluating the torque on a sphere, rotating with a constant angular velocity
+ */
+//*******************************************************************************************************************
+template< typename ParticleAccessor_T >
+class TorqueEval
+{
+ public:
+ TorqueEval(SweepTimeloop* timeloop, Setup* setup, const shared_ptr< StructuredBlockStorage >& blocks,
+ const shared_ptr< ParticleAccessor_T >& ac, bool fileIO)
+ : timeloop_(timeloop), setup_(setup), blocks_(blocks), ac_(ac), fileIO_(fileIO), torqueOld_(0.0), torqueNew_(0.0)
+ {
+ // calculate the (semi)analytical torque value
+ // see also Hofmann et al. - Hydrodynamic interactions in colloidal crystals:(II). Application to dense cubic and
+ // tetragonal arrays (1999), Eqs. 5.1 and 5.5
+ const real_t S = real_c(1.95708);
+ setup_->analyticalTorque =
+ -setup_->visc *
+ (real_c(6) * setup_->phi / (real_c(1) - setup_->phi - S * std::pow(setup_->phi, real_c(10. / 3.)))) *
+ setup_->angularVel * real_c(setup_->length * setup_->length * setup_->length);
+
+ if (fileIO_)
+ {
+ std::ofstream file;
+ filename_ = "TorqueSpherePSM.txt";
+ WALBERLA_ROOT_SECTION()
+ {
+ file.open(filename_.c_str());
+ file << "#\t torqueSim\t torqueAnaly\n";
+ file.close();
+ }
+ }
+ }
+
+ // evaluate the acting torque
+ void operator()()
+ {
+ const uint_t timestep(timeloop_->getCurrentTimeStep() + 1);
+
+ if (timestep % setup_->checkFrequency != 0) return;
+
+ // update torque values
+ torqueOld_ = torqueNew_;
+ torqueNew_ = calculateTorque();
+
+ // write to file if desired
+ WALBERLA_ROOT_SECTION()
+ {
+ if (fileIO_)
+ {
+ std::ofstream file;
+ file.open(filename_.c_str(), std::ofstream::app);
+ file.setf(std::ios::unitbuf);
+ file.precision(15);
+ file << timestep << " " << torqueNew_ << " " << setup_->analyticalTorque << "\n";
+ file.close();
+ }
+ }
+ }
+
+ // obtain the torque acting on the sphere by summing up all the process local parts
+ real_t calculateTorque()
+ {
+ real_t torque = real_c(0);
+ for (auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt)
+ {
+ for (size_t idx = 0; idx < ac_->size(); ++idx)
+ {
+ torque += ac_->getHydrodynamicTorque(idx)[1];
+ }
+ }
+
+ WALBERLA_MPI_SECTION() { mpi::allReduceInplace(torque, mpi::SUM); }
+ return torque;
+ }
+
+ // return the relative temporal change in the torque
+ real_t getTorqueDiff() const { return std::fabs((torqueNew_ - torqueOld_) / torqueNew_); }
+
+ // return the torque
+ real_t getTorque() const { return torqueNew_; }
+
+ private:
+ SweepTimeloop* timeloop_;
+
+ Setup* setup_;
+
+ shared_ptr< StructuredBlockStorage > blocks_;
+ shared_ptr< ParticleAccessor_T > ac_;
+
+ bool fileIO_;
+ std::string filename_;
+
+ real_t torqueOld_;
+ real_t torqueNew_;
+};
+
+//////////
+// MAIN //
+//////////
+
+//*******************************************************************************************************************
+/*!\brief Testcase that checks the torque acting on a constantly rotating sphere in the center of a cubic domain
+ *
+ * The torque for this problem (often denoted as Simple Cubic setup) is given by a semi-analytical formula.
+ * The cubic domain is periodic in all directions, making it a physically infinite periodic array of spheres.
+ \verbatim
+ _______________
+ | <- |
+ | ___ |
+ | / \ |
+ | | x | |
+ | \___/ |
+ | -> |
+ |_______________|
+
+ \endverbatim
+ *
+ * The collision model used for the LBM is TRT with a relaxation parameter tau=1.5 and the magic parameter 3/16.
+ * The Stokes approximation of the equilibrium PDFs is used.
+ * The sphere rotates with a angular velocity of 1e-5.
+ * The domain is 32x32x32, and the sphere has a diameter of around 27 cells ( chi * domainlength )
+ * The simulation is run until the relative change in the torque between 100 time steps is less than 1e-5.
+ * The pe is not used since the angular velocity is kept constant and the force is explicitly reset after each time
+ step.
+ *
+ */
+//*******************************************************************************************************************
+
+int main(int argc, char** argv)
+{
+ debug::enterTestMode();
+
+ mpi::Environment env(argc, argv);
+
+ auto processes = MPIManager::instance()->numProcesses();
+
+ if (processes != 1 && processes != 2 && processes != 4 && processes != 8)
+ {
+ std::cerr << "Number of processes must be equal to either 1, 2, 4, or 8!" << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ ///////////////////
+ // Customization //
+ ///////////////////
+
+ bool shortrun = false;
+ bool funcTest = false;
+ bool fileIO = false;
+ bool SC1W1 = false;
+ bool SC2W1 = false;
+ bool SC3W1 = false;
+ bool SC1W2 = false;
+ bool SC2W2 = false;
+ bool SC3W2 = false;
+
+ for (int i = 1; i < argc; ++i)
+ {
+ if (std::strcmp(argv[i], "--shortrun") == 0)
+ {
+ shortrun = true;
+ continue;
+ }
+ if (std::strcmp(argv[i], "--funcTest") == 0)
+ {
+ funcTest = true;
+ continue;
+ }
+ if (std::strcmp(argv[i], "--fileIO") == 0)
+ {
+ fileIO = true;
+ continue;
+ }
+ if (std::strcmp(argv[i], "--SC1W1") == 0)
+ {
+ SC1W1 = true;
+ continue;
+ }
+ if (std::strcmp(argv[i], "--SC2W1") == 0)
+ {
+ SC2W1 = true;
+ continue;
+ }
+ if (std::strcmp(argv[i], "--SC3W1") == 0)
+ {
+ SC3W1 = true;
+ continue;
+ }
+ if (std::strcmp(argv[i], "--SC1W2") == 0)
+ {
+ SC1W2 = true;
+ continue;
+ }
+ if (std::strcmp(argv[i], "--SC2W2") == 0)
+ {
+ SC2W2 = true;
+ continue;
+ }
+ if (std::strcmp(argv[i], "--SC3W2") == 0)
+ {
+ SC3W2 = true;
+ continue;
+ }
+ WALBERLA_ABORT("Unrecognized command line argument found: " << argv[i]);
+ }
+
+ if (!SC1W1 && !SC2W1 && !SC3W1 && !SC1W2 && !SC2W2 && !SC3W2)
+ {
+ std::cerr << "Specify the model (--SC_W_) you want to use for the partially saturated cells method!" << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ ///////////////////////////
+ // SIMULATION PROPERTIES //
+ ///////////////////////////
+
+ Setup setup;
+
+ setup.length = uint_t(32); // length of the cubic domain in lattice cells
+ const real_t chi = real_c(0.85); // porosity parameter: diameter / length
+ setup.tau = real_c(1.5); // relaxation time
+ setup.angularVel = real_c(1e-5); // angular velocity of the sphere
+ setup.checkFrequency = uint_t(100); // evaluate the torque only every checkFrequency time steps
+ setup.radius = real_c(0.5) * chi * real_c(setup.length); // sphere radius
+ setup.visc = (setup.tau - real_c(0.5)) / real_c(3); // viscosity in lattice units
+ setup.phi = real_c(4.0 / 3.0) * math::pi * setup.radius * setup.radius * setup.radius /
+ (real_c(setup.length * setup.length * setup.length)); // solid volume fraction
+ const real_t omega = real_c(1) / setup.tau; // relaxation rate
+ const real_t dx = real_c(1); // lattice dx
+ const real_t convergenceLimit = real_c(1e-5); // tolerance for relative change in torque
+ const uint_t timesteps =
+ funcTest ? 1 : (shortrun ? uint_c(150) : uint_c(5000)); // maximum number of time steps for the whole simulation
+
+ ///////////////////////////
+ // BLOCK STRUCTURE SETUP //
+ ///////////////////////////
+
+ const uint_t XBlocks = (processes >= 2) ? uint_t(2) : uint_t(1);
+ const uint_t YBlocks = (processes >= 4) ? uint_t(2) : uint_t(1);
+ const uint_t ZBlocks = (processes == 8) ? uint_t(2) : uint_t(1);
+ const uint_t XCells = setup.length / XBlocks;
+ const uint_t YCells = setup.length / YBlocks;
+ const uint_t ZCells = setup.length / ZBlocks;
+
+ // create fully periodic domain
+ auto blocks = blockforest::createUniformBlockGrid(XBlocks, YBlocks, ZBlocks, XCells, YCells, ZCells, dx, true, true,
+ true, true);
+
+ ////////
+ // PE //
+ ////////
+
+ auto mesapdDomain = std::make_shared< mesa_pd::domain::BlockForestDomain >(blocks->getBlockForestPointer());
+ auto ps = std::make_shared< mesa_pd::data::ParticleStorage >(1);
+ auto ss = std::make_shared< mesa_pd::data::ShapeStorage >();
+ using ParticleAccessor_T = mesa_pd::data::ParticleAccessorWithShape;
+ auto accessor = walberla::make_shared< ParticleAccessor_T >(ps, ss);
+
+ /////////////////
+ // PE COUPLING //
+ /////////////////
+
+ // connect to pe
+ const real_t overlap = real_c(1.5) * dx;
+
+ if (setup.radius > real_c(setup.length) * real_c(0.5) - overlap)
+ {
+ std::cerr << "Periodic sphere is too large!" << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ // create the sphere in the middle of the domain
+ Vector3< real_t > position(real_c(setup.length) * real_c(0.5));
+ auto sphereShape = ss->create< mesa_pd::data::Sphere >(setup.radius);
+
+ if (mesapdDomain->isContainedInProcessSubdomain(uint_c(walberla::mpi::MPIManager::instance()->rank()), position))
+ {
+ auto sphereParticle = ps->create();
+ sphereParticle->setShapeID(sphereShape);
+ sphereParticle->setType(0);
+ sphereParticle->setPosition(position);
+ sphereParticle->setAngularVelocity(Vector3(real_c(0), setup.angularVel, real_c(0)));
+ sphereParticle->setOwner(walberla::MPIManager::instance()->rank());
+ sphereParticle->setInteractionRadius(setup.radius);
+ }
+
+ // synchronize often enough for large bodies
+ std::function< void(void) > syncCall = [&]() {
+ // const real_t overlap = real_t(1.5) * dx;
+ mesa_pd::mpi::SyncNextNeighbors syncNextNeighborFunc;
+ syncNextNeighborFunc(*ps, *mesapdDomain, overlap);
+ };
+
+ syncCall();
+
+ ///////////////////////
+ // ADD DATA TO BLOCKS //
+ ////////////////////////
+
+ // create the lattice model
+ LatticeModel_T latticeModel = LatticeModel_T(omega);
+
+ // add PDF field ( uInit = <0,0,0>, rhoInit = 1 )
+ BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >(
+ blocks, "pdf field (zyxf)", latticeModel, Vector3< real_t >(real_c(0), real_c(0), real_c(0)), real_c(1),
+ uint_t(1), field::zyxf);
+
+ // add particle and volume fraction field
+ BlockDataID particleAndVolumeFractionFieldID =
+ field::addToStorage< lbm_mesapd_coupling::psm::ParticleAndVolumeFractionField_T >(
+ blocks, "particle and volume fraction field",
+ std::vector< lbm_mesapd_coupling::psm::ParticleAndVolumeFraction_T >(), field::zyxf, 0);
+ // map bodies and calculate solid volume fraction initially
+ lbm_mesapd_coupling::psm::ParticleAndVolumeFractionMapping particleMapping(
+ blocks, accessor, lbm_mesapd_coupling::RegularParticlesSelector(), particleAndVolumeFractionFieldID, 4);
+ particleMapping();
+
+ // initialize the PDF field for PSM
+ if (SC1W1 || SC2W1 || SC3W1)
+ {
+ lbm_mesapd_coupling::psm::initializeDomainForPSM< LatticeModel_T, 1 >(
+ *blocks, pdfFieldID, particleAndVolumeFractionFieldID, *accessor);
+ }
+ else
+ {
+ lbm_mesapd_coupling::psm::initializeDomainForPSM< LatticeModel_T, 2 >(
+ *blocks, pdfFieldID, particleAndVolumeFractionFieldID, *accessor);
+ }
+
+ ///////////////
+ // TIME LOOP //
+ ///////////////
+
+ // create the timeloop
+ SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);
+
+ // setup of the LBM communication for synchronizing the pdf field between neighboring blocks
+ std::function< void() > commFunction;
+
+ blockforest::communication::UniformBufferedScheme< stencil::D3Q27 > scheme(blocks);
+ scheme.addPackInfo(make_shared< field::communication::PackInfo< PdfField_T > >(pdfFieldID));
+ commFunction = scheme;
+
+ using TorqueEval_T = TorqueEval< ParticleAccessor_T >;
+ shared_ptr< TorqueEval_T > torqueEval = make_shared< TorqueEval_T >(&timeloop, &setup, blocks, accessor, fileIO);
+
+ if (SC1W1)
+ {
+ auto sweep = lbm_mesapd_coupling::psm::makePSMSweep< LatticeModel_T, 1, 1 >(
+ pdfFieldID, particleAndVolumeFractionFieldID, blocks, accessor);
+ // communication, streaming and force evaluation
+ timeloop.add() << BeforeFunction(commFunction, "LBM Communication")
+ << Sweep(makeSharedSweep(sweep), "cell-wise LB sweep")
+ << AfterFunction(SharedFunctor< TorqueEval_T >(torqueEval), "torque evaluation");
+ }
+ else if (SC2W1)
+ {
+ auto sweep = lbm_mesapd_coupling::psm::makePSMSweep< LatticeModel_T, 2, 1 >(
+ pdfFieldID, particleAndVolumeFractionFieldID, blocks, accessor);
+ // communication, streaming and force evaluation
+ timeloop.add() << BeforeFunction(commFunction, "LBM Communication")
+ << Sweep(makeSharedSweep(sweep), "cell-wise LB sweep")
+ << AfterFunction(SharedFunctor< TorqueEval_T >(torqueEval), "torque evaluation");
+ }
+ else if (SC3W1)
+ {
+ auto sweep = lbm_mesapd_coupling::psm::makePSMSweep< LatticeModel_T, 3, 1 >(
+ pdfFieldID, particleAndVolumeFractionFieldID, blocks, accessor);
+ // communication, streaming and force evaluation
+ timeloop.add() << BeforeFunction(commFunction, "LBM Communication")
+ << Sweep(makeSharedSweep(sweep), "cell-wise LB sweep")
+ << AfterFunction(SharedFunctor< TorqueEval_T >(torqueEval), "torque evaluation");
+ }
+ else if (SC1W2)
+ {
+ auto sweep = lbm_mesapd_coupling::psm::makePSMSweep< LatticeModel_T, 1, 2 >(
+ pdfFieldID, particleAndVolumeFractionFieldID, blocks, accessor);
+ // communication, streaming and force evaluation
+ timeloop.add() << BeforeFunction(commFunction, "LBM Communication")
+ << Sweep(makeSharedSweep(sweep), "cell-wise LB sweep")
+ << AfterFunction(SharedFunctor< TorqueEval_T >(torqueEval), "torque evaluation");
+ }
+ else if (SC2W2)
+ {
+ auto sweep = lbm_mesapd_coupling::psm::makePSMSweep< LatticeModel_T, 2, 2 >(
+ pdfFieldID, particleAndVolumeFractionFieldID, blocks, accessor);
+ // communication, streaming and force evaluation
+ timeloop.add() << BeforeFunction(commFunction, "LBM Communication")
+ << Sweep(makeSharedSweep(sweep), "cell-wise LB sweep")
+ << AfterFunction(SharedFunctor< TorqueEval_T >(torqueEval), "torque evaluation");
+ }
+ else if (SC3W2)
+ {
+ auto sweep = lbm_mesapd_coupling::psm::makePSMSweep< LatticeModel_T, 3, 2 >(
+ pdfFieldID, particleAndVolumeFractionFieldID, blocks, accessor);
+ // communication, streaming and force evaluation
+ timeloop.add() << BeforeFunction(commFunction, "LBM Communication")
+ << Sweep(makeSharedSweep(sweep), "cell-wise LB sweep")
+ << AfterFunction(SharedFunctor< TorqueEval_T >(torqueEval), "torque evaluation");
+ }
+
+ lbm_mesapd_coupling::ResetHydrodynamicForceTorqueKernel resetHydrodynamicForceTorque;
+
+ timeloop.addFuncAfterTimeStep(RemainingTimeLogger(timeloop.getNrOfTimeSteps()), "Remaining Time Logger");
+
+ ////////////////////////
+ // EXECUTE SIMULATION //
+ ////////////////////////
+
+ WcTimingPool timeloopTiming;
+
+ // time loop
+ for (uint_t i = 0; i < timesteps; ++i)
+ {
+ // perform a single simulation step
+ timeloop.singleStep(timeloopTiming);
+
+ // resetting force
+ ps->forEachParticle(false, mesa_pd::kernel::SelectAll(), *accessor, resetHydrodynamicForceTorque, *accessor);
+
+ // check if the relative change in the torque is below the specified convergence criterion
+ if (i > setup.checkFrequency && torqueEval->getTorqueDiff() < convergenceLimit)
+ {
+ // if simulation has converged, terminate simulation
+ break;
+ }
+ }
+
+ timeloopTiming.logResultOnRoot();
+
+ // check the result
+ if (!funcTest && !shortrun)
+ {
+ real_t relErr = std::fabs((setup.analyticalTorque - torqueEval->getTorque()) / setup.analyticalTorque);
+ if (fileIO)
+ {
+ WALBERLA_ROOT_SECTION()
+ {
+ std::cout << "Analytical torque: " << setup.analyticalTorque << "\n"
+ << "Simulated torque: " << torqueEval->getTorque() << "\n"
+ << "Relative error: " << relErr << "\n";
+ }
+ }
+ // the relative error has to be below 10% (25% for SC2)
+ WALBERLA_CHECK_LESS(relErr, (SC2W1 || SC2W2) ? real_c(0.25) : real_c(0.1));
+ }
+
+ return 0;
+}
+
+} // namespace torque_sphere_psm
+
+int main(int argc, char** argv) { torque_sphere_psm::main(argc, argv); }
diff --git a/tests/lbm_mesapd_coupling/partially_saturated_cells_method/Utility.h b/tests/lbm_mesapd_coupling/partially_saturated_cells_method/Utility.h
new file mode 100644
index 000000000..f24a4332a
--- /dev/null
+++ b/tests/lbm_mesapd_coupling/partially_saturated_cells_method/Utility.h
@@ -0,0 +1,85 @@
+//======================================================================================================================
+//
+// 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 Utility.h
+//! \ingroup lbm_mesapd_coupling
+//! \author Samuel Kemmler <samuel.kemmler@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/ShapeStorage.h"
+
+namespace walberla
+{
+namespace lbm_mesapd_coupling
+{
+
+//*******************************************************************************************************************
+/*!\brief Calculating the sum over all fraction values. This can be used as a sanity check since it has to be roughly
+ * equal to the volume of all particles.
+ *
+ */
+//*******************************************************************************************************************
+class FractionFieldSum
+{
+ public:
+ FractionFieldSum(const shared_ptr< StructuredBlockStorage >& blockStorage,
+ BlockDataID particleAndVolumeFractionFieldID)
+ : blockStorage_(blockStorage), particleAndVolumeFractionFieldID_(particleAndVolumeFractionFieldID)
+ {}
+
+ real_t operator()()
+ {
+ real_t sum = 0.0;
+
+ for (auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt)
+ {
+ psm::ParticleAndVolumeFractionField_T* particleAndVolumeFractionField =
+ blockIt->getData< psm::ParticleAndVolumeFractionField_T >(particleAndVolumeFractionFieldID_);
+
+ const cell_idx_t xSize = cell_idx_c(particleAndVolumeFractionField->xSize());
+ const cell_idx_t ySize = cell_idx_c(particleAndVolumeFractionField->ySize());
+ const cell_idx_t zSize = cell_idx_c(particleAndVolumeFractionField->zSize());
+
+ for (cell_idx_t z = 0; z < zSize; ++z)
+ {
+ for (cell_idx_t y = 0; y < ySize; ++y)
+ {
+ for (cell_idx_t x = 0; x < xSize; ++x)
+ {
+ for (auto particleAndVolumeFraction : particleAndVolumeFractionField->get(x, y, z))
+ {
+ sum += particleAndVolumeFraction.second;
+ }
+ }
+ }
+ }
+ }
+
+ WALBERLA_MPI_SECTION() { mpi::allReduceInplace(sum, mpi::SUM); }
+
+ return sum;
+ }
+
+ private:
+ shared_ptr< StructuredBlockStorage > blockStorage_;
+ BlockDataID particleAndVolumeFractionFieldID_;
+};
+
+} // namespace lbm_mesapd_coupling
+} // namespace walberla
--
GitLab