apps/showcases/CMakeLists.txt

 add_subdirectory( BidisperseFluidizedBed )
+add_subdirectory( ChargedParticles )
 add_subdirectory( CombinedResolvedUnresolved )
 add_subdirectory( FluidizedBed )
 add_subdirectory( FreeSurface )

apps/showcases/ChargedParticles/CMakeLists.txt

+waLBerla_link_files_to_builddir( "*.prm" )
+
+waLBerla_add_executable ( NAME ChargedParticles FILES ChargedParticles.cpp
+        DEPENDS blockforest boundary core pde domain_decomposition field lbm lbm_mesapd_coupling mesa_pd timeloop vtk )
+
+waLBerla_add_executable ( NAME TestIterativeSolvers FILES TestIterativeSolvers.cpp
+        DEPENDS blockforest boundary core pde domain_decomposition field timeloop vtk )

apps/showcases/ChargedParticles/ChargeForce.h

+//
+// Created by RichardAngersbach on 16.01.2023.
+//
+
+#ifndef WALBERLA_CHARGEFORCE_H
+#define WALBERLA_CHARGEFORCE_H
+
+namespace walberla
+{
+using Stencil_T = stencil::D3Q7;
+typedef GhostLayerField< real_t, 1 > ScalarField_T;
+typedef GhostLayerField< real_t, 3 > VectorField_T;
+
+class ChargeForceUpdate
+{
+ public:
+   ChargeForceUpdate(const BlockDataID& potential, const BlockDataID& chargeForce, const std::shared_ptr< StructuredBlockForest >& blocks)
+      : potential_(potential), chargeForce_(chargeForce), blocks_(blocks) {}
+
+   void operator()() {
+      // get charge force with FD gradient from electric potential
+      for (auto block = blocks_->begin(); block!=blocks_->end(); ++block) {
+         VectorField_T* chargeForce = block->getData< VectorField_T >(chargeForce_);
+         ScalarField_T* potential = block->getData< ScalarField_T >(potential_);
+
+         WALBERLA_FOR_ALL_CELLS_XYZ(potential,
+            chargeForce->get(x, y, z, 0) = (real_c(1) / (real_c(2) * (blocks_->dx()))) * (potential->get(x + 1, y, z) - potential->get(x - 1, y, z));
+            chargeForce->get(x, y, z, 1) = (real_c(1) / (real_c(2) * (blocks_->dy()))) * (potential->get(x, y + 1, z) - potential->get(x, y - 1, z));
+            chargeForce->get(x, y, z, 2) = (real_c(1) / (real_c(2) * (blocks_->dz()))) * (potential->get(x, y, z + 1) - potential->get(x, y, z - 1));
+         )
+      }
+   }
+
+ private:
+   BlockDataID potential_;
+   BlockDataID chargeForce_;
+
+   std::shared_ptr< StructuredBlockForest > blocks_;
+};
+
+} // namespace walberla
+
+#endif // WALBERLA_CHARGEFORCE_H

apps/showcases/ChargedParticles/DirichletDomainBoundary.h

+#ifndef WALBERLA_DIRICHLETDOMAINBOUNDARY_H
+#define WALBERLA_DIRICHLETDOMAINBOUNDARY_H
+
+namespace walberla
+{
+template< typename PdeField >
+class DirichletDomainBoundary
+{
+ public:
+   DirichletDomainBoundary(StructuredBlockStorage& blocks, const BlockDataID& fieldId)
+      : blocks_(blocks), fieldId_(fieldId)
+   {
+      for (uint_t i = 0; i != stencil::D3Q6::Size; ++i)
+      {
+         includeBoundary_[i] = true;
+         values_[i] = real_t(0);
+      }
+   }
+
+   void includeBoundary(const stencil::Direction& direction) { includeBoundary_[stencil::D3Q6::idx[direction]] = true; }
+   void excludeBoundary(const stencil::Direction& direction)
+   {
+      includeBoundary_[stencil::D3Q6::idx[direction]] = false;
+   }
+
+   void setValue( const real_t value ) { for( uint_t i = 0; i != stencil::D3Q6::Size; ++i ) values_[i] = value; }
+   void setValue( const stencil::Direction & direction, const real_t value ) { values_[stencil::D3Q6::idx[direction]] = value; }
+
+   void operator()();
+
+ protected:
+   // first-order dirichlet boundary conditions with scalar value
+   void apply(PdeField* p, const CellInterval& interval, const cell_idx_t cx, const cell_idx_t cy, const cell_idx_t cz,
+              const real_t value) const;
+
+   StructuredBlockStorage& blocks_;
+   BlockDataID fieldId_;
+
+   bool includeBoundary_[stencil::D3Q6::Size];
+
+   real_t values_[ stencil::D3Q6::Size ];
+
+}; // class DirichletDomainBoundary
+
+template< typename PdeField >
+void DirichletDomainBoundary< PdeField >::operator()()
+{
+   for (auto block = blocks_.begin(); block != blocks_.end(); ++block)
+   {
+      PdeField* p = block->template getData< PdeField >(fieldId_);
+
+      if (includeBoundary_[stencil::D3Q6::idx[stencil::W]] && blocks_.atDomainXMinBorder(*block))
+      {
+         apply(p,
+               CellInterval(
+                  cell_idx_t(-1), cell_idx_t(0)                         , cell_idx_t(0),
+                  cell_idx_t(-1), cell_idx_c(p->ySize()) - cell_idx_t(1), cell_idx_c(p->zSize()) - cell_idx_t(1)),
+               cell_idx_t(1), cell_idx_t(0), cell_idx_t(0), values_[ stencil::D3Q6::idx[ stencil::W ] ]);
+      }
+      if (includeBoundary_[stencil::D3Q6::idx[stencil::E]] && blocks_.atDomainXMaxBorder(*block))
+      {
+         apply(p,
+               CellInterval(
+                  cell_idx_c(p->xSize()), cell_idx_t(0)                         , cell_idx_t(0),
+                  cell_idx_c(p->xSize()), cell_idx_c(p->ySize()) - cell_idx_t(1), cell_idx_c(p->zSize()) - cell_idx_t(1)),
+               cell_idx_t(-1), cell_idx_t(0), cell_idx_t(0), values_[ stencil::D3Q6::idx[ stencil::E ] ]);
+      }
+
+      if (includeBoundary_[stencil::D3Q6::idx[stencil::S]] && blocks_.atDomainYMinBorder(*block))
+      {
+         apply(p,
+               CellInterval(
+                  cell_idx_t(0)                         , cell_idx_t(-1), cell_idx_t(0),
+                  cell_idx_c(p->xSize()) - cell_idx_t(1), cell_idx_t(-1), cell_idx_c(p->zSize()) - cell_idx_t(1)),
+               cell_idx_t(0), cell_idx_t(1), cell_idx_t(0), values_[ stencil::D3Q6::idx[ stencil::S ] ]);
+      }
+      if (includeBoundary_[stencil::D3Q6::idx[stencil::N]] && blocks_.atDomainYMaxBorder(*block))
+      {
+         apply(p,
+               CellInterval(
+                  cell_idx_t(0)                         , cell_idx_c(p->ySize()), cell_idx_t(0),
+                  cell_idx_c(p->xSize()) - cell_idx_t(1), cell_idx_c(p->ySize()), cell_idx_c(p->zSize()) - cell_idx_t(1)),
+               cell_idx_t(0), cell_idx_t(-1), cell_idx_t(0), values_[ stencil::D3Q6::idx[ stencil::N ] ]);
+      }
+
+      if (includeBoundary_[stencil::D3Q6::idx[stencil::B]] && blocks_.atDomainZMinBorder(*block))
+      {
+         apply(p,
+               CellInterval(
+                  cell_idx_t(0)                         , cell_idx_t(0)                         , cell_idx_t(-1),
+                  cell_idx_c(p->xSize()) - cell_idx_t(1), cell_idx_c(p->ySize()) - cell_idx_t(1), cell_idx_t(-1)),
+               cell_idx_t(0), cell_idx_t(0), cell_idx_t(1), values_[ stencil::D3Q6::idx[ stencil::B ] ]);
+      }
+      if (includeBoundary_[stencil::D3Q6::idx[stencil::T]] && blocks_.atDomainZMaxBorder(*block))
+      {
+         apply(p,
+               CellInterval(
+                  cell_idx_t(0)                         , cell_idx_t(0)                         , cell_idx_c(p->zSize()),
+                  cell_idx_c(p->xSize()) - cell_idx_t(1), cell_idx_c(p->ySize()) - cell_idx_t(1), cell_idx_c(p->zSize())),
+               cell_idx_t(0), cell_idx_t(0), cell_idx_t(-1), values_[ stencil::D3Q6::idx[ stencil::T ] ]);
+      }
+   }
+}
+
+template< typename PdeField >
+void DirichletDomainBoundary< PdeField >::apply(PdeField* p, const CellInterval& interval, const cell_idx_t cx,
+                                                const cell_idx_t cy, const cell_idx_t cz, const real_t value) const
+{
+   WALBERLA_FOR_ALL_CELLS_IN_INTERVAL_XYZ(interval,
+                                          p->get(x, y, z) = real_c(2) * value - p->get(x + cx, y + cy, z + cz);)
+}
+
+template< typename PdeField >
+class DirichletFunctionDomainBoundary
+{
+ public:
+
+   using ApplyFunction = std::function< void (IBlock* block, PdeField* p, const CellInterval& interval, const cell_idx_t cx, const cell_idx_t cy, const cell_idx_t cz) >;
+
+   DirichletFunctionDomainBoundary(StructuredBlockStorage& blocks, const BlockDataID& fieldId)
+      : blocks_(blocks), fieldId_(fieldId)
+   {
+      for (uint_t i = 0; i != stencil::D3Q6::Size; ++i)
+      {
+         includeBoundary_[i] = true;
+
+         applyFunctions_[i] = {};
+      }
+   }
+
+   void includeBoundary(const stencil::Direction& direction) { includeBoundary_[stencil::D3Q6::idx[direction]] = true; }
+   void excludeBoundary(const stencil::Direction& direction)
+   {
+      includeBoundary_[stencil::D3Q6::idx[direction]] = false;
+   }
+
+   void setFunction( const ApplyFunction func ) { for( uint_t i = 0; i != stencil::D3Q6::Size; ++i ) applyFunctions_[i] = func; }
+   void setFunction( const stencil::Direction & direction, const ApplyFunction func ) { applyFunctions_[stencil::D3Q6::idx[direction]] = func; }
+
+   void operator()();
+
+ protected:
+   StructuredBlockStorage& blocks_;
+   BlockDataID fieldId_;
+
+   bool includeBoundary_[stencil::D3Q6::Size];
+
+   // user-defined apply function
+   ApplyFunction applyFunctions_[ stencil::D3Q6::Size ];
+
+}; // class DirichletFunctionDomainBoundary
+
+template< typename PdeField >
+void DirichletFunctionDomainBoundary< PdeField >::operator()()
+{
+   for (auto blockIt = blocks_.begin(); blockIt != blocks_.end(); ++blockIt)
+   {
+      auto * block = static_cast<blockforest::Block*> (&(*blockIt));
+
+      PdeField* p = block->template getData< PdeField >(fieldId_);
+
+      if (applyFunctions_[stencil::D3Q6::idx[stencil::W]] && includeBoundary_[stencil::D3Q6::idx[stencil::W]] && blocks_.atDomainXMinBorder(*block))
+      {
+         applyFunctions_[stencil::D3Q6::idx[stencil::W]](block, p,
+               CellInterval(
+                  cell_idx_t(-1), cell_idx_t(0)                         , cell_idx_t(0),
+                  cell_idx_t(-1), cell_idx_c(p->ySize()) - cell_idx_t(1), cell_idx_c(p->zSize()) - cell_idx_t(1)),
+               cell_idx_t(1), cell_idx_t(0), cell_idx_t(0));
+      }
+      if (applyFunctions_[stencil::D3Q6::idx[stencil::E]] && includeBoundary_[stencil::D3Q6::idx[stencil::E]] && blocks_.atDomainXMaxBorder(*block))
+      {
+         applyFunctions_[stencil::D3Q6::idx[stencil::E]](block, p,
+               CellInterval(
+                  cell_idx_c(p->xSize()), cell_idx_t(0)                         , cell_idx_t(0),
+                  cell_idx_c(p->xSize()), cell_idx_c(p->ySize()) - cell_idx_t(1), cell_idx_c(p->zSize()) - cell_idx_t(1)),
+               cell_idx_t(-1), cell_idx_t(0), cell_idx_t(0));
+      }
+
+      if (applyFunctions_[stencil::D3Q6::idx[stencil::S]] && includeBoundary_[stencil::D3Q6::idx[stencil::S]] && blocks_.atDomainYMinBorder(*block))
+      {
+         applyFunctions_[stencil::D3Q6::idx[stencil::S]](block, p,
+               CellInterval(
+                  cell_idx_t(0)                         , cell_idx_t(-1), cell_idx_t(0),
+                  cell_idx_c(p->xSize()) - cell_idx_t(1), cell_idx_t(-1), cell_idx_c(p->zSize()) - cell_idx_t(1)),
+               cell_idx_t(0), cell_idx_t(1), cell_idx_t(0));
+      }
+      if (applyFunctions_[stencil::D3Q6::idx[stencil::N]] && includeBoundary_[stencil::D3Q6::idx[stencil::N]] && blocks_.atDomainYMaxBorder(*block))
+      {
+         applyFunctions_[stencil::D3Q6::idx[stencil::N]](block, p,
+               CellInterval(
+                  cell_idx_t(0)                         , cell_idx_c(p->ySize()), cell_idx_t(0),
+                  cell_idx_c(p->xSize()) - cell_idx_t(1), cell_idx_c(p->ySize()), cell_idx_c(p->zSize()) - cell_idx_t(1)),
+               cell_idx_t(0), cell_idx_t(-1), cell_idx_t(0));
+      }
+
+      if (applyFunctions_[stencil::D3Q6::idx[stencil::B]] && includeBoundary_[stencil::D3Q6::idx[stencil::B]] && blocks_.atDomainZMinBorder(*block))
+      {
+         applyFunctions_[stencil::D3Q6::idx[stencil::B]](block, p,
+               CellInterval(
+                  cell_idx_t(0)                         , cell_idx_t(0)                         , cell_idx_t(-1),
+                  cell_idx_c(p->xSize()) - cell_idx_t(1), cell_idx_c(p->ySize()) - cell_idx_t(1), cell_idx_t(-1)),
+               cell_idx_t(0), cell_idx_t(0), cell_idx_t(1));
+      }
+      if (applyFunctions_[stencil::D3Q6::idx[stencil::T]] && includeBoundary_[stencil::D3Q6::idx[stencil::T]] && blocks_.atDomainZMaxBorder(*block))
+      {
+         applyFunctions_[stencil::D3Q6::idx[stencil::T]](block, p,
+               CellInterval(
+                  cell_idx_t(0)                         , cell_idx_t(0)                         , cell_idx_c(p->zSize()),
+                  cell_idx_c(p->xSize()) - cell_idx_t(1), cell_idx_c(p->ySize()) - cell_idx_t(1), cell_idx_c(p->zSize())),
+               cell_idx_t(0), cell_idx_t(0), cell_idx_t(-1));
+      }
+   }
+}
+
+} // namespace walberla
+
+#endif // WALBERLA_DIRICHLETDOMAINBOUNDARY_H

apps/showcases/ChargedParticles/PoissonSolver.h

+#ifndef WALBERLA_POISSONSOLVER_H
+#define WALBERLA_POISSONSOLVER_H
+
+#include "pde/all.h"
+#include "DirichletDomainBoundary.h"
+
+enum Enum { WALBERLA_JACOBI, WALBERLA_SOR, DAMPED_JACOBI };
+
+namespace walberla
+{
+
+// typedefs
+using Stencil_T = stencil::D3Q7;
+typedef GhostLayerField< real_t, 1 > ScalarField_T;
+
+template< Enum solver, bool useDirichlet >
+class PoissonSolver
+{
+ public:
+
+   void dampedJacobiSweep(IBlock* const block) {
+      ScalarField_T* srcField = block->getData< ScalarField_T >(src_);
+      ScalarField_T* dstField = block->getData< ScalarField_T >(dst_);
+      ScalarField_T* rhsField = block->getData< ScalarField_T >(rhs_);
+
+      const real_t omega          = real_c(0.6);
+      const real_t invLaplaceDiag = real_t(1) / laplaceWeights_[Stencil_T::idx[stencil::C]];
+
+      WALBERLA_ASSERT_GREATER_EQUAL(srcField->nrOfGhostLayers(), 1);
+
+      WALBERLA_FOR_ALL_CELLS_XYZ(srcField,
+                                 real_t stencilTimesSrc = real_c(0);
+                                 for (auto dir = Stencil_T::begin(); dir != Stencil_T::end(); ++dir)
+                                    stencilTimesSrc += laplaceWeights_[dir.toIdx()] * srcField->getNeighbor(x, y, z, *dir);
+
+                                 dstField->get(x, y, z) = srcField->get(x, y, z) + omega * invLaplaceDiag * (rhsField->get(x, y, z) - stencilTimesSrc);
+      )
+
+      srcField->swapDataPointers(dstField);
+   }
+
+   PoissonSolver(const BlockDataID& src, const BlockDataID& dst, const BlockDataID& rhs,
+                 const std::shared_ptr< StructuredBlockForest >& blocks,
+                 uint_t iterations = uint_t(1000),
+                 real_t residualNormThreshold = real_c(1e-4),
+                 uint_t residualCheckFrequency = uint_t(100),
+                 const std::function< void () >& boundaryHandling = {})
+      : src_(src), dst_(dst), rhs_(rhs), blocks_(blocks), boundaryHandling_(boundaryHandling) {
+
+      // stencil weights
+
+      laplaceWeights_ = std::vector < real_t > (Stencil_T::Size, real_c(0));
+      laplaceWeights_[Stencil_T::idx[stencil::C]] = real_t( 2) / (blocks_->dx() * blocks_->dx()) +
+                                             real_t( 2) / (blocks_->dy() * blocks_->dy()) +
+                                             real_t( 2) / (blocks_->dz() * blocks_->dz());
+      laplaceWeights_[Stencil_T::idx[stencil::T]] = real_t(-1) / (blocks_->dz() * blocks_->dz());
+      laplaceWeights_[Stencil_T::idx[stencil::B]] = real_t(-1) / (blocks_->dz() * blocks_->dz());
+      laplaceWeights_[Stencil_T::idx[stencil::N]] = real_t(-1) / (blocks_->dy() * blocks_->dy());
+      laplaceWeights_[Stencil_T::idx[stencil::S]] = real_t(-1) / (blocks_->dy() * blocks_->dy());
+      laplaceWeights_[Stencil_T::idx[stencil::E]] = real_t(-1) / (blocks_->dx() * blocks_->dx());
+      laplaceWeights_[Stencil_T::idx[stencil::W]] = real_t(-1) / (blocks_->dx() * blocks_->dx());
+
+      // communication
+
+      commScheme_ = make_shared< blockforest::communication::UniformBufferedScheme< Stencil_T > >(blocks_);
+      commScheme_->addPackInfo(make_shared< field::communication::PackInfo< ScalarField_T > >(src_));
+      commScheme_->addPackInfo(make_shared< field::communication::PackInfo< ScalarField_T > >(rhs_));
+
+      // boundary handling
+
+      if (!boundaryHandling_) {
+         if constexpr (useDirichlet) {
+            // dirichlet BCs
+            boundaryHandling_ = DirichletDomainBoundary< ScalarField_T >(*blocks_, src_);
+         } else {
+            // neumann BCs
+            boundaryHandling_ = pde::NeumannDomainBoundary< ScalarField_T >(*blocks_, src_);
+         }
+      }
+
+      // res norm
+
+      residualNorm_ = make_shared< pde::ResidualNorm< Stencil_T > >(blocks_->getBlockStorage(), src_, rhs_, laplaceWeights_);
+
+      // jacobi
+
+      jacobiFixedSweep_ = make_shared < pde::JacobiFixedStencil< Stencil_T > >(src_, dst_, rhs_, laplaceWeights_);
+
+      // use custom impl with damping or jacobi from waLBerla
+      std::function< void ( IBlock * ) > jacSweep = {};
+      if (solver == DAMPED_JACOBI) {
+         jacSweep = [this](IBlock* block) { dampedJacobiSweep(block); };
+      } else {
+         jacSweep = *jacobiFixedSweep_;
+      }
+
+      jacobiIteration_ = std::make_unique< pde::JacobiIteration >(
+         blocks_->getBlockStorage(), iterations,
+         *commScheme_,
+         jacSweep,
+         *residualNorm_, residualNormThreshold, residualCheckFrequency);
+
+      jacobiIteration_->addBoundaryHandling(boundaryHandling_);
+
+      // SOR
+
+      real_t omega = real_t(2) / real_t(3);
+
+      sorFixedSweep_ = make_shared< pde::SORFixedStencil< Stencil_T > >(blocks, src_, rhs_, laplaceWeights_, omega);
+
+      sorIteration_ = std::make_unique< pde::RBGSIteration >(
+         blocks_->getBlockStorage(), iterations,
+         *commScheme_,
+         sorFixedSweep_->getRedSweep(), sorFixedSweep_->getBlackSweep(),
+         *residualNorm_, residualNormThreshold, residualCheckFrequency);
+
+      sorIteration_->addBoundaryHandling(boundaryHandling);
+   }
+
+   // get approximate solution of electric potential
+   void operator()() {
+      if constexpr (solver != WALBERLA_SOR) {
+         (*jacobiIteration_)();
+      } else {
+         (*sorIteration_)();
+      }
+   }
+
+ private:
+   BlockDataID src_;
+   BlockDataID dst_;
+   BlockDataID rhs_;
+
+   std::vector< real_t > laplaceWeights_;
+   std::shared_ptr< StructuredBlockForest > blocks_;
+   std::shared_ptr< blockforest::communication::UniformBufferedScheme< Stencil_T > > commScheme_;
+
+   std::function< void () > boundaryHandling_;
+
+   std::shared_ptr < pde::ResidualNorm< Stencil_T > > residualNorm_;
+
+   std::shared_ptr< pde::JacobiFixedStencil< Stencil_T > > jacobiFixedSweep_;
+   std::unique_ptr< pde::JacobiIteration > jacobiIteration_;
+
+   std::shared_ptr< pde::SORFixedStencil< Stencil_T  > > sorFixedSweep_;
+   std::unique_ptr< pde::RBGSIteration > sorIteration_;
+};
+
+} // namespace walberla
+
+#endif // WALBERLA_POISSONSOLVER_H

apps/showcases/ChargedParticles/TestIterativeSolvers.cpp

+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "boundary/all.h"
+
+#include "core/Environment.h"
+#include "core/grid_generator/SCIterator.h"
+#include "core/logging/all.h"
+#include "core/waLBerlaBuildInfo.h"
+
+#include "field/AddToStorage.h"
+#include "field/vtk/all.h"
+
+#include <core/math/all.h>
+
+#include "PoissonSolver.h"
+
+namespace walberla {
+
+enum Testcase { TEST_DIRICHLET_1, TEST_DIRICHLET_2, TEST_NEUMANN };
+
+using ScalarField_T = GhostLayerField< real_t, 1 >;
+
+template < typename PdeField, Testcase testcase >
+void applyDirichletFunction(const shared_ptr< StructuredBlockStorage > & blocks, const math::AABB & domainAABB, const stencil::Direction& direction,
+                            IBlock* block, PdeField* p, const CellInterval& interval, const cell_idx_t cx, const cell_idx_t cy, const cell_idx_t cz) {
+
+   WALBERLA_FOR_ALL_CELLS_IN_INTERVAL_XYZ(interval,
+      real_t boundaryCoord_x = 0.;
+      real_t boundaryCoord_y = 0.;
+      real_t boundaryCoord_z = 0.;
+
+      const auto cellAABB = blocks->getBlockLocalCellAABB(*block, Cell(x, y, z));
+      auto cellCenter     = cellAABB.center();
+
+      // snap cell position to actual domain position
+      switch (direction) {
+         case stencil::W:
+            boundaryCoord_x = domainAABB.xMin();
+            boundaryCoord_y = cellCenter[1];
+            boundaryCoord_z = cellCenter[2];
+            break;
+         case stencil::E:
+            boundaryCoord_x = domainAABB.xMax();
+            boundaryCoord_y = cellCenter[1];
+            boundaryCoord_z = cellCenter[2];
+            break;
+         case stencil::S:
+            boundaryCoord_x = cellCenter[0];
+            boundaryCoord_y = domainAABB.yMin();
+            boundaryCoord_z = cellCenter[2];
+            break;
+         case stencil::N:
+            boundaryCoord_x = cellCenter[0];
+            boundaryCoord_y = domainAABB.yMax();
+            boundaryCoord_z = cellCenter[2];
+            break;
+         case stencil::B:
+            boundaryCoord_x = cellCenter[0];
+            boundaryCoord_y = cellCenter[1];
+            boundaryCoord_z = domainAABB.zMin();
+            break;
+         case stencil::T:
+            boundaryCoord_x = cellCenter[0];
+            boundaryCoord_y = cellCenter[1];
+            boundaryCoord_z = domainAABB.zMax();
+            break;
+         default:
+            WALBERLA_ABORT("Unknown direction");
+      }
+
+      // use positions normalized to unit cube
+      boundaryCoord_x /= domainAABB.size(0);
+      boundaryCoord_y /= domainAABB.size(1);
+      boundaryCoord_z /= domainAABB.size(2);
+
+      auto funcVal = real_c(0);
+      switch (testcase) {
+         case TEST_DIRICHLET_1:
+            funcVal = (boundaryCoord_x * boundaryCoord_x) - real_c(0.5) * (boundaryCoord_y * boundaryCoord_y) - real_c(0.5) * (boundaryCoord_z * boundaryCoord_z);
+            break;
+         case TEST_DIRICHLET_2:
+            funcVal = real_c( sin ( math::pi * boundaryCoord_x ) ) *
+                      real_c( sin ( math::pi * boundaryCoord_y ) ) *
+                      real_c( math::root_two * math::pi * boundaryCoord_z );
+            break;
+         default:
+            WALBERLA_ABORT("Unknown testcase");
+      }
+      p->get(x, y, z) = real_c(2) * funcVal - p->get(x + cx, y + cy, z + cz);
+   )
+}
+
+void resetSolution(const shared_ptr< StructuredBlockStorage > & blocks, BlockDataID & solution, BlockDataID & solutionCpy) {
+   for (auto block = blocks->begin(); block != blocks->end(); ++block) {
+      ScalarField_T* solutionField    = block->getData< ScalarField_T >(solution);
+      ScalarField_T* solutionFieldCpy = block->getData< ScalarField_T >(solutionCpy);
+
+      // reset fields
+      solutionField->set(real_c(0));
+      solutionFieldCpy->set(real_c(0));
+   }
+}
+
+void resetRHS(const shared_ptr< StructuredBlockStorage > & blocks, BlockDataID & rhs) {
+   for (auto block = blocks->begin(); block != blocks->end(); ++block) {
+      ScalarField_T* rhsField = block->getData< ScalarField_T >(rhs);
+
+      // reset field
+      rhsField->set(real_c(0));
+   }
+}
+
+// solve two different scenarios (dirichlet scenario and neumann scenario) with different analytical solutions and setups
+template < Testcase testcase >
+void solve(const shared_ptr< StructuredBlockForest > & blocks,
+           const math::AABB & domainAABB, BlockDataID & solution, BlockDataID & solutionCpy, BlockDataID & rhs) {
+
+   const bool useDirichlet = testcase == TEST_DIRICHLET_1 || testcase == TEST_DIRICHLET_2;
+
+   // set boundary handling depending on scenario
+   std::function< void () > boundaryHandling = {};
+
+   if constexpr (useDirichlet) {
+      // set dirichlet function per domain face
+      auto dirichletFunction = DirichletFunctionDomainBoundary< ScalarField_T >(*blocks, solution);
+
+#define GET_BOUNDARY_LAMBDA(dir) \
+         [&blocks, &domainAABB](IBlock* block, ScalarField_T* p, const CellInterval& interval, const cell_idx_t cx, const cell_idx_t cy, const cell_idx_t cz) { \
+            applyDirichletFunction< ScalarField_T, testcase >(blocks, domainAABB, dir, block, p, interval, cx, cy, cz); \
+         }
+
+      dirichletFunction.setFunction(stencil::W, GET_BOUNDARY_LAMBDA(stencil::W));
+      dirichletFunction.setFunction(stencil::E, GET_BOUNDARY_LAMBDA(stencil::E));
+      dirichletFunction.setFunction(stencil::S, GET_BOUNDARY_LAMBDA(stencil::S));
+      dirichletFunction.setFunction(stencil::N, GET_BOUNDARY_LAMBDA(stencil::N));
+      dirichletFunction.setFunction(stencil::B, GET_BOUNDARY_LAMBDA(stencil::B));
+      dirichletFunction.setFunction(stencil::T, GET_BOUNDARY_LAMBDA(stencil::T));
+
+      boundaryHandling = dirichletFunction;
+   }
+
+   // solvers: Jacobi and SOR
+
+   auto numIter = uint_c(50000);
+   auto resThres = real_c(1e-10);
+   auto resCheckFreq = uint_c(1000);
+
+   auto poissonSolverJacobi = PoissonSolver< WALBERLA_JACOBI, useDirichlet > (solution, solutionCpy, rhs, blocks, numIter, resThres, resCheckFreq, boundaryHandling);
+   auto poissonSolverDampedJac = PoissonSolver< DAMPED_JACOBI, useDirichlet > (solution, solutionCpy, rhs, blocks, numIter, resThres, resCheckFreq, boundaryHandling);
+   auto poissonSolverSOR = PoissonSolver< WALBERLA_SOR, useDirichlet > (solution, solutionCpy, rhs, blocks, numIter, resThres, resCheckFreq, boundaryHandling);
+
+   // calc error depending on scenario
+
+   auto computeMaxError = [&blocks, &solution, &domainAABB]() {
+      real_t error = real_c(0);
+
+      for (auto block = blocks->begin(); block != blocks->end(); ++block) {
+         ScalarField_T* solutionField = block->getData< ScalarField_T >(solution);
+
+         real_t blockResult = real_t(0);
+         WALBERLA_FOR_ALL_CELLS_XYZ_OMP(solutionField, omp parallel for schedule(static) reduction(max: blockResult),
+                                    const auto cellAABB = blocks->getBlockLocalCellAABB(*block, Cell(x,y,z));
+                                    auto cellCenter = cellAABB.center();
+
+                                    // use positions normalized to unit cube
+                                    real_t scaleX = real_c(1) / domainAABB.size(0);
+                                    real_t scaleY = real_c(1) / domainAABB.size(1);
+                                    real_t scaleZ = real_c(1) / domainAABB.size(2);
+
+                                    real_t posX = cellCenter[0] * scaleX;
+                                    real_t posY = cellCenter[1] * scaleY;
+                                    real_t posZ = cellCenter[2] * scaleZ;
+
+                                    real_t analyticalSol;
+
+                                    // analytical solution of problem with neumann/dirichlet boundaries
+                                    switch (testcase) {
+                                       case TEST_DIRICHLET_1:
+                                          analyticalSol =                 (posX * posX)
+                                                          - real_c(0.5) * (posY * posY)
+                                                          - real_c(0.5) * (posZ * posZ);
+                                          break;
+                                       case TEST_DIRICHLET_2:
+                                          analyticalSol = real_c( sin ( math::pi * posX ) ) *
+                                                          real_c( sin ( math::pi * posY ) ) *
+                                                          real_c( sinh ( math::root_two * math::pi * posZ ) );
+                                          break;
+                                       case TEST_NEUMANN:
+                                          analyticalSol = real_c( cos ( real_c(2) * math::pi * posX ) ) *
+                                                          real_c( cos ( real_c(2) * math::pi * posY ) ) *
+                                                          real_c( cos ( real_c(2) * math::pi * posZ ) );
+                                          break;
+                                       default:
+                                          WALBERLA_ABORT("Unknown testcase");
+                                    }
+
+                                    real_t currErr = real_c(fabs(solutionField->get(x, y, z) - analyticalSol));
+                                    blockResult = std::max(blockResult, currErr);
+         )
+         error = std::max(error, blockResult);
+      }
+      mpi::allReduceInplace( error, mpi::MAX );
+
+      return error;
+   };
+
+   // init rhs depending on scenario
+
+   for (auto block = blocks->begin(); block != blocks->end(); ++block) {
+      ScalarField_T* rhsField = block->getData< ScalarField_T >(rhs);
+
+      WALBERLA_FOR_ALL_CELLS_XYZ(
+         rhsField,
+
+         const auto cellAABB = blocks->getBlockLocalCellAABB(*block, Cell(x, y, z));
+         auto cellCenter = cellAABB.center();
+
+         // use positions normalized to unit cube
+         real_t scaleX = real_c(1) / domainAABB.size(0);
+         real_t scaleY = real_c(1) / domainAABB.size(1);
+         real_t scaleZ = real_c(1) / domainAABB.size(2);
+
+         real_t posX = cellCenter[0] * scaleX;
+         real_t posY = cellCenter[1] * scaleY;
+         real_t posZ = cellCenter[2] * scaleZ;
+
+         switch (testcase) {
+            case TEST_DIRICHLET_1:
+               rhsField->get(x, y, z) = real_c(0);
+               break;
+            case TEST_DIRICHLET_2:
+               rhsField->get(x, y, z) = real_c( -( math::pi * math::pi ) * ( -( scaleX * scaleX ) - ( scaleY * scaleY ) + real_c(2) * ( scaleZ * scaleZ ) ) ) *
+                                        real_c( sin ( math::pi * posX ) ) *
+                                        real_c( sin ( math::pi * posY ) ) *
+                                        real_c( sinh ( math::root_two * math::pi * posZ ) );
+               break;
+            case TEST_NEUMANN:
+               rhsField->get(x, y, z) = real_c(4) * math::pi * math::pi *
+                                        real_c( (pow(scaleX, 2) + pow(scaleY, 2) + pow(scaleZ, 2)) ) *
+                                        real_c( cos(real_c(2) * math::pi * posX) ) *
+                                        real_c( cos(real_c(2) * math::pi * posY) ) *
+                                        real_c( cos(real_c(2) * math::pi * posZ) );
+               break;
+            default:
+               WALBERLA_ABORT("Unknown testcase");
+         }
+      )
+   }
+
+   auto initError = computeMaxError();
+   WALBERLA_LOG_INFO_ON_ROOT("Initial error is: " << initError);
+
+   // solve with jacobi
+   WALBERLA_LOG_INFO_ON_ROOT("-- Solve using Jacobi --");
+   poissonSolverJacobi();
+   auto errJac = computeMaxError();
+   WALBERLA_LOG_INFO_ON_ROOT("Error after Jacobi solver is: " << errJac);
+
+   // solve with damped jacobi
+   WALBERLA_LOG_INFO_ON_ROOT("-- Solve using (damped) Jacobi --");
+   resetSolution(blocks, solution, solutionCpy); // reset solutions and solve anew
+   poissonSolverDampedJac();
+   auto errDampedJac = computeMaxError();
+   WALBERLA_LOG_INFO_ON_ROOT("Error after (damped) Jacobi solver is: " << errDampedJac);
+
+   // solve with SOR
+   WALBERLA_LOG_INFO_ON_ROOT("-- Solve using SOR --");
+   resetSolution(blocks, solution, solutionCpy); // reset solutions and solve anew
+   poissonSolverSOR();
+   auto errSOR = computeMaxError();
+   WALBERLA_LOG_INFO_ON_ROOT("Error after SOR solver is: " << errSOR);
+}
+
+// solve two different charged particle scenarios (dirichlet scenario and neumann scenario) with different setups
+template < bool useDirichlet >
+void solveChargedParticles(const shared_ptr< StructuredBlockForest > & blocks,
+                           const math::AABB & domainAABB, BlockDataID & solution, BlockDataID & solutionCpy, BlockDataID & rhs) {
+
+   // solvers: Jacobi and SOR
+
+   auto numIter = uint_c(20000);
+   auto resThres = real_c(1e-5);
+   auto resCheckFreq = uint_c(1000);
+
+   auto poissonSolverJacobi = PoissonSolver< DAMPED_JACOBI, useDirichlet > (solution, solutionCpy, rhs, blocks, numIter, resThres, resCheckFreq);
+   auto poissonSolverSOR = PoissonSolver< WALBERLA_SOR, useDirichlet > (solution, solutionCpy, rhs, blocks, numIter, resThres, resCheckFreq);
+
+   // init rhs with two charged particles
+
+   for (auto block = blocks->begin(); block != blocks->end(); ++block) {
+      ScalarField_T* rhsField = block->getData< ScalarField_T >(rhs);
+
+      WALBERLA_FOR_ALL_CELLS_XYZ(
+         rhsField,
+
+         const auto cellAABB = blocks->getBlockLocalCellAABB(*block, Cell(x, y, z));
+         auto cellCenter = cellAABB.center();
+
+         const real_t x0 = domainAABB.xMin() + real_c(0.45) * domainAABB.size(0);
+         const real_t y0 = domainAABB.yMin() + real_c(0.45) * domainAABB.size(1);
+         const real_t z0 = domainAABB.zMin() + real_c(0.45) * domainAABB.size(2);
+         const real_t r0 = real_c(0.08) * domainAABB.size(0);
+         const real_t s0 = real_c(1);
+
+         const real_t x1 = domainAABB.xMin() + real_c(0.65) * domainAABB.size(0);
+         const real_t y1 = domainAABB.yMin() + real_c(0.65) * domainAABB.size(1);
+         const real_t z1 = domainAABB.zMin() + real_c(0.65) * domainAABB.size(2);
+         const real_t r1 = real_c(0.08) * domainAABB.size(0);
+         const real_t s1 = real_c(1);
+
+         real_t posX = cellCenter[0];
+         real_t posY = cellCenter[1];
+         real_t posZ = cellCenter[2];
+
+         if ( ( real_c( pow( posX - x0, 2 ) ) + real_c( pow( posY - y0, 2 ) ) + real_c( pow( posZ - z0, 2 ) ) ) < real_c( pow( r0, 2 ) ) ) {
+            auto relDistPos0 = real_c( sqrt( real_c( pow( ( posX - x0 ) / r0, 2 ) ) + real_c( pow( ( posY - y0 ) / r0, 2 ) ) + real_c( pow( ( posZ - z0 ) / r0, 2 ) ) ) );
+            rhsField->get(x, y, z) = -s0 * ( real_c(1) - relDistPos0 );
+         } else if ( ( real_c ( pow( posX - x1, 2 ) ) + real_c ( pow( posY - y1, 2 ) ) + real_c ( pow( posZ - z1, 2 ) ) ) < real_c ( pow( r1, 2 ) ) ) {
+            auto relDistPos1 = real_c( sqrt( real_c( pow( ( posX - x1 ) / r1, 2 ) ) + real_c( pow( ( posY - y1 ) / r1, 2 ) ) + real_c( pow( ( posZ - z1 ) / r1, 2 ) ) ) );
+            rhsField->get(x, y, z) = -s1 * ( real_c(1) - relDistPos1);
+         } else {
+            rhsField->get(x, y, z) = real_c(0);
+         }
+      )
+   }
+
+   // solve with jacobi
+   poissonSolverJacobi();
+
+   // solve with SOR
+   resetSolution(blocks, solution, solutionCpy); // reset solutions and solve anew
+   poissonSolverSOR();
+}
+
+int main(int argc, char** argv)
+{
+   Environment env(argc, argv);
+
+   WALBERLA_LOG_INFO_ON_ROOT("waLBerla revision: " << std::string(WALBERLA_GIT_SHA1).substr(0, 8));
+
+   logging::Logging::instance()->setLogLevel(logging::Logging::LogLevel::DETAIL);
+   logging::Logging::instance()->includeLoggingToFile( "log" );
+
+   ///////////////////////////
+   // BLOCK STRUCTURE SETUP //
+   ///////////////////////////
+
+   auto domainAABB = math::AABB(0, 0, 0, 125, 50, 250);
+   WALBERLA_LOG_INFO_ON_ROOT("Domain sizes are: x = " << domainAABB.size(0) << ", y = " << domainAABB.size(1) << ", z = " << domainAABB.size(2));
+
+   shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid(
+      domainAABB,
+      uint_c( 1), uint_c( 1), uint_c( 1), // number of blocks in x,y,z direction
+      uint_c( 125), uint_c( 50), uint_c( 250), // how many cells per block (x,y,z)
+      true,                               // max blocks per process
+      false, false, false,                // periodicity
+      false);
+
+   BlockDataID rhs = field::addToStorage< ScalarField_T >(blocks, "rhs", 0, field::fzyx, 1);
+   BlockDataID solution = field::addToStorage< ScalarField_T >(blocks, "solution", 0, field::fzyx, 1);
+   BlockDataID solutionCpy = field::addCloneToStorage< ScalarField_T >(blocks, solution, "solution (copy)");
+
+   // first solve neumann problem...
+   WALBERLA_LOG_INFO_ON_ROOT("Run analytical test cases...")
+   WALBERLA_LOG_INFO_ON_ROOT("- Solving analytical Neumann problem with Jacobi and SOR... -")
+   solve< TEST_NEUMANN > (blocks, domainAABB, solution, solutionCpy, rhs);
+
+   // ... then solve dirichlet problem
+   resetRHS(blocks, rhs); // reset fields and solve anew
+   resetSolution(blocks, solution, solutionCpy);
+   WALBERLA_LOG_INFO_ON_ROOT("- Solving analytical Dirichlet problem (1) with Jacobi and SOR... -")
+   solve< TEST_DIRICHLET_1 > (blocks, domainAABB, solution, solutionCpy, rhs);
+
+   resetRHS(blocks, rhs); // reset fields and solve anew
+   resetSolution(blocks, solution, solutionCpy);
+   WALBERLA_LOG_INFO_ON_ROOT("- Solving analytical Dirichlet problem (2) with Jacobi and SOR... -")
+   solve< TEST_DIRICHLET_2 > (blocks, domainAABB, solution, solutionCpy, rhs);
+
+   // ... charged particle test
+
+   WALBERLA_LOG_INFO_ON_ROOT("Run charged particle test cases...")
+   resetRHS(blocks, rhs); // reset fields and solve anew
+   resetSolution(blocks, solution, solutionCpy);
+
+   // neumann
+   WALBERLA_LOG_INFO_ON_ROOT("- Run charged particles with Neumann boundaries... -")
+   solveChargedParticles < false > (blocks, domainAABB, solution, solutionCpy, rhs);
+
+   // dirichlet
+   WALBERLA_LOG_INFO_ON_ROOT("- Run charged particles with Dirichlet (val = 0) boundaries... -")
+   resetSolution(blocks, solution, solutionCpy); // reset solutions and solve anew
+   solveChargedParticles < true > (blocks, domainAABB, solution, solutionCpy, rhs);
+
+   return EXIT_SUCCESS;
+}
+
+} // namespace walberla
+
+int main(int argc, char** argv) { walberla::main(argc, argv); }
+

apps/showcases/ChargedParticles/Utility.h

+//======================================================================================================================
+//
+//  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
+//! \author Samuel Kemmler <samuel.kemmler@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "lbm_mesapd_coupling/DataTypes.h"
+
+namespace walberla
+{
+namespace charged_particles
+{
+
+template< typename BlockStorage_T >
+void computeChargeDensity(const shared_ptr< BlockStorage_T >& blocks,
+                          const BlockDataID& particleAndVolumeFractionFieldID, const BlockDataID& chargeDensityFieldID)
+{
+   // TODO: compute physically correct charge density here using the particle charges, have a look at src/lbm_mesapd_coupling/partially_saturated_cells_method/ParticleAndVolumeFractionMapping.h for how to iterate over particles and cells
+   for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
+   {
+      lbm_mesapd_coupling::psm::ParticleAndVolumeFractionField_T* particleAndVolumeFractionField =
+         blockIt->template getData< lbm_mesapd_coupling::psm::ParticleAndVolumeFractionField_T >(
+            particleAndVolumeFractionFieldID);
+      GhostLayerField< real_t, 1 >* chargeDensityField =
+         blockIt->template getData< GhostLayerField< real_t, 1 > >(chargeDensityFieldID);
+
+      WALBERLA_FOR_ALL_CELLS_XYZ(particleAndVolumeFractionField, chargeDensityField->get(x, y, z) = 0.0;
+                                 for (auto& e
+                                      : particleAndVolumeFractionField->get(x, y, z))
+                                    // TODO: try out with former implementation:
+                                    // chargeDensityField->get(x, y, z) -= e.second
+                                    if (e.second > real_c(0))
+                                       chargeDensityField->get(x, y, z) = real_c(-1);) // rhs depends on the negative charge density
+   }
+}
+
+} // namespace charged_particles
+} // namespace walberla

apps/showcases/ChargedParticles/input.prm

+PhysicalSetup // all to be specified in SI units!
+{
+    xSize 0.01; // = width
+    ySize 0.005; // = depth
+    zSize 0.02; // = height
+
+    periodicInX false;
+    periodicInY false;
+
+    runtime 10.0;
+
+    uInflow 0.01;
+    gravitationalAcceleration 9.81;
+
+    kinematicViscosityFluid 1e-5;
+    densityFluid 1000.;
+
+    particleDiameter 0.002;
+    densityParticle 1100.;
+    dynamicFrictionCoefficient 0.15;
+    coefficientOfRestitution 0.6;
+
+    particleGenerationSpacing 0.00401;
+}
+
+NumericalSetup
+{
+    dx 0.0002; // in m
+    uInflow 0.01; // in LBM units, should be smaller than 0.1, this then determines dt
+
+    // product of number of blocks should be equal to number of used processes
+    numXBlocks 2;
+    numYBlocks 1;
+    numZBlocks 2;
+
+    useLubricationForces true;
+    numberOfParticleSubCycles 10;
+}
+
+Output
+{
+    infoSpacing 0.01; // in s
+
+    vtkSpacingParticles 0.1; // in s
+    vtkSpacingFluid 0.1; // in s
+    vtkFolder vtk_out;
+}

apps/showcases/FluidizedBed/CMakeLists.txt

 waLBerla_link_files_to_builddir( "*.prm" )
 
-waLBerla_add_executable ( NAME FluidizedBed FILES FluidizedBed.cpp
-                          DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+waLBerla_add_executable ( NAME FluidizedBedMEM FILES FluidizedBedMEM.cpp
+                          DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd timeloop vtk )
+
+waLBerla_add_executable ( NAME FluidizedBedPSM FILES FluidizedBedPSM.cpp
+        DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd timeloop vtk )

src/lbm_mesapd_coupling/partially_saturated_cells_method/ParticleAndVolumeFractionMapping.h

@@ -77,6 +77,22 @@ class ParticleAndVolumeFractionMapping
       {
          if (mappingParticleSelector_(idx, *ac_)) { update(idx); }
       }
+
+      // normalize the sum of all overlap fractions of a cell to 1
+      for (auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt)
+      {
+         ParticleAndVolumeFractionField_T* particleAndVolumeFractionField =
+            blockIt->getData< ParticleAndVolumeFractionField_T >(particleAndVolumeFractionFieldID_);
+
+         WALBERLA_FOR_ALL_CELLS_XYZ(
+            particleAndVolumeFractionField, real_t fractionSum = 0.0;
+            for (auto& e
+                 : particleAndVolumeFractionField->get(x, y, z)) fractionSum += e.second;
+            if (fractionSum > 1.0) {
+               for (auto& e : particleAndVolumeFractionField->get(x, y, z))
+                  e.second /= fractionSum;
+            })
+      }
    }
 
  private: