diff --git a/apps/showcases/ChargedParticles/ChargedParticles.cpp b/apps/showcases/ChargedParticles/ChargedParticles.cpp
index 00b09979efbb3d0b11649c656cfd9bd9337cde87..8bdf6d8770dab0314c78f8c978b00c665443a50a 100644
--- a/apps/showcases/ChargedParticles/ChargedParticles.cpp
+++ b/apps/showcases/ChargedParticles/ChargedParticles.cpp
@@ -86,14 +86,15 @@
 
 #include "vtk/all.h"
 
+#include "./poisson_solver/CustomBoundary.h"
+#include "./poisson_solver/PotentialValidationCustomBoundary.h"
 #include "AddElectrostaticInteractionKernel.h"
+#include "ChargeDensity.h"
 #include "ChargeForce.h"
-#include "poisson_solver/PoissonSolver.h"
+#include "PostProcessingUtilities.h"
 #include "ResetElectrostaticForceKernel.h"
-#include "ChargeDensity.h"
-#include "postProcessingUtilities.h"
-#include "./poisson_solver/PotentialValidationCustomBoundary.h"
-#include "./poisson_solver/CustomBoundary.h"
+#include "poisson_solver/PoissonSolver.h"
+#include "ErrorNorms.h"
 
 namespace charged_particles
 {
@@ -278,32 +279,35 @@ void createPlaneSetup(const shared_ptr< mesa_pd::data::ParticleStorage >& ps,
 
 struct ParticleInfo
 {
-   real_t averageVelocity = 0_r;
-   real_t maximumVelocity = 0_r;
-   uint_t numParticles    = 0;
-   real_t maximumHeight   = 0_r;
-   real_t particleVolume  = 0_r;
-   real_t heightOfMass    = 0_r;
+   real_t averageVelocity           = 0_r;
+   real_t ensembledAverageVelocityZ = 0_r;
+   real_t maximumVelocity           = 0_r;
+   uint_t numParticles              = 0;
+   real_t maximumHeight             = 0_r;
+   real_t particleVolume            = 0_r;
+   real_t heightOfMass              = 0_r;
 
    void allReduce()
    {
       walberla::mpi::allReduceInplace(numParticles, walberla::mpi::SUM);
       walberla::mpi::allReduceInplace(averageVelocity, walberla::mpi::SUM);
+      walberla::mpi::allReduceInplace(ensembledAverageVelocityZ, walberla::mpi::SUM);
       walberla::mpi::allReduceInplace(maximumVelocity, walberla::mpi::MAX);
       walberla::mpi::allReduceInplace(maximumHeight, walberla::mpi::MAX);
       walberla::mpi::allReduceInplace(particleVolume, walberla::mpi::SUM);
       walberla::mpi::allReduceInplace(heightOfMass, walberla::mpi::SUM);
 
       averageVelocity /= real_c(numParticles);
+      ensembledAverageVelocityZ /= real_c(numParticles);
       heightOfMass /= particleVolume;
    }
 };
 
 std::ostream& operator<<(std::ostream& os, ParticleInfo const& m)
 {
-   return os << "Particle Info: uAvg = " << m.averageVelocity << ", uMax = " << m.maximumVelocity
-             << ", numParticles = " << m.numParticles << ", zMax = " << m.maximumHeight << ", Vp = " << m.particleVolume
-             << ", zMass = " << m.heightOfMass;
+   return os << "Ensembled Avg Z = " << m.ensembledAverageVelocityZ << ", Particle Info: uAvg = " << m.averageVelocity
+             << ", uMax = " << m.maximumVelocity << ", numParticles = " << m.numParticles
+             << ", zMax = " << m.maximumHeight << ", Vp = " << m.particleVolume << ", zMass = " << m.heightOfMass;
 }
 
 template< typename Accessor_T >
@@ -319,9 +323,11 @@ ParticleInfo evaluateParticleInfo(const Accessor_T& ac)
 
       ++info.numParticles;
       real_t velMagnitude   = ac.getLinearVelocity(i).length();
+      real_t velocityZ      = ac.getLinearVelocity(i)[2];
       real_t particleVolume = ac.getShape(i)->getVolume();
       real_t height         = ac.getPosition(i)[2];
       info.averageVelocity += velMagnitude;
+      info.ensembledAverageVelocityZ += velocityZ;
       info.maximumVelocity = std::max(info.maximumVelocity, velMagnitude);
       info.maximumHeight   = std::max(info.maximumHeight, height);
       info.particleVolume += particleVolume;
@@ -444,6 +450,8 @@ int main(int argc, char** argv)
    const uint_t numXBlocks                = numericalSetup.getParameter< uint_t >("numXBlocks");
    const uint_t numYBlocks                = numericalSetup.getParameter< uint_t >("numYBlocks");
    const uint_t numZBlocks                = numericalSetup.getParameter< uint_t >("numZBlocks");
+   const real_t absResThreshold           = numericalSetup.getParameter< real_t >("absResThreshold");
+   const uint_t jacobiIterations          = numericalSetup.getParameter< uint_t >("jacobiIterations");
    const bool useLubricationForces        = numericalSetup.getParameter< bool >("useLubricationForces");
    const uint_t numberOfParticleSubCycles = numericalSetup.getParameter< uint_t >("numberOfParticleSubCycles");
 
@@ -462,7 +470,7 @@ int main(int argc, char** argv)
    if (on == 0) { useIntegrator = false; }
 
    Config::BlockHandle ParticleGenerationSchemes = cfgFile->getBlock("ParticleGenerationScheme");
-   const bool extendSimulationDomain        = ParticleGenerationSchemes.getParameter< bool >("extendSimulationDomain");
+   const bool extendSimulationDomain = ParticleGenerationSchemes.getParameter< bool >("extendSimulationDomain");
 
    Config::BlockHandle boundaryTypes = cfgFile->getBlock("BoundaryTypes");
    std::string boundaryTypeNorth     = boundaryTypes.getParameter< std::string >("North_type");
@@ -514,30 +522,32 @@ int main(int argc, char** argv)
 
    Vector3< uint_t > domainSize(uint_c((xSize_SI / dx_SI)), uint_c((ySize_SI / dx_SI)), uint_c((zSize_SI / dx_SI)));
 
-   WALBERLA_LOG_INFO_ON_ROOT("domain x size:" << " " << std::ceil(xSize_SI / dx_SI));
-   WALBERLA_LOG_INFO_ON_ROOT("domain x size without ceil:" << " " << (xSize_SI / dx_SI));
-   WALBERLA_LOG_INFO_ON_ROOT("ceil of 800 :" << " " << std::ceil(800));
-
-   WALBERLA_CHECK_FLOAT_EQUAL(real_t(domainSize[0]) * dx_SI, xSize_SI, "domain size in x is not divisible by given dx");
-   WALBERLA_CHECK_FLOAT_EQUAL(real_t(domainSize[1]) * dx_SI, ySize_SI, "domain size in y is not divisible by given dx");
-   WALBERLA_CHECK_FLOAT_EQUAL(real_t(domainSize[2]) * dx_SI, zSize_SI, "domain size in z is not divisible by given dx");
-
    Vector3< uint_t > cellsPerBlockPerDirection(domainSize[0] / numXBlocks, domainSize[1] / numYBlocks,
                                                domainSize[2] / numZBlocks);
 
-   WALBERLA_CHECK_EQUAL(domainSize[0], cellsPerBlockPerDirection[0] * numXBlocks,
-                        "number of cells in x of " << domainSize[0]
-                                                   << " is not divisible by given number of blocks in x direction");
-   WALBERLA_CHECK_EQUAL(domainSize[1], cellsPerBlockPerDirection[1] * numYBlocks,
-                        "number of cells in y of " << domainSize[1]
-                                                   << " is not divisible by given number of blocks in y direction");
-   WALBERLA_CHECK_EQUAL(domainSize[2], cellsPerBlockPerDirection[2] * numZBlocks,
-                        "number of cells in z of " << domainSize[2]
-                                                   << " is not divisible by given number of blocks in z direction");
-
-   WALBERLA_CHECK_GREATER(
-      particleDiameter_SI / dx_SI, 5_r,
-      "Your numerical resolution is below 5 cells per diameter and thus too small for such simulations!");
+   if (simulationName != "pivVelocityValidation")
+   {
+      WALBERLA_CHECK_FLOAT_EQUAL(real_t(domainSize[0]) * dx_SI, xSize_SI,
+                                 "domain size in x is not divisible by given dx");
+      WALBERLA_CHECK_FLOAT_EQUAL(real_t(domainSize[1]) * dx_SI, ySize_SI,
+                                 "domain size in y is not divisible by given dx");
+      WALBERLA_CHECK_FLOAT_EQUAL(real_t(domainSize[2]) * dx_SI, zSize_SI,
+                                 "domain size in z is not divisible by given dx");
+
+      WALBERLA_CHECK_EQUAL(domainSize[0], cellsPerBlockPerDirection[0] * numXBlocks,
+                           "number of cells in x of " << domainSize[0]
+                                                      << " is not divisible by given number of blocks in x direction");
+      WALBERLA_CHECK_EQUAL(domainSize[1], cellsPerBlockPerDirection[1] * numYBlocks,
+                           "number of cells in y of " << domainSize[1]
+                                                      << " is not divisible by given number of blocks in y direction");
+      WALBERLA_CHECK_EQUAL(domainSize[2], cellsPerBlockPerDirection[2] * numZBlocks,
+                           "number of cells in z of " << domainSize[2]
+                                                      << " is not divisible by given number of blocks in z direction");
+
+      WALBERLA_CHECK_GREATER(
+         particleDiameter_SI / dx_SI, 5_r,
+         "Your numerical resolution is below 5 cells per diameter and thus too small for such simulations!");
+   }
 
    const real_t densityRatio           = densityParticle_SI / densityFluid_SI;
    const real_t ReynoldsNumberParticle = uInflow_SI * particleDiameter_SI / kinematicViscosityFluid_SI;
@@ -567,28 +577,31 @@ int main(int argc, char** argv)
    // relative_permitivity = 78 and vacuum_permitivity SI units ->  A2⋅s4⋅kg−1⋅m−3
 
    const real_t vacuum_permitivity_SI = real_c(78.5 * 8.8541878128 * pow(10, -12));
-   const real_t Ampere_unit           = ((densityFluid_SI) * real_c(pow(dx_SI, 5))) / (real_c(pow(dt_SI, 3)) * V);
-   real_t vacuum_permitivity =
-      ((densityFluid_SI * real_c(pow(dx_SI, 3))) * real_c(pow(dx_SI, 3))) / (real_c(pow(dt_SI, 4)) * Ampere_unit * Ampere_unit);
+   const real_t Ampere_unit           = ((densityFluid_SI) *real_c(pow(dx_SI, 5))) / (real_c(pow(dt_SI, 3)) * V);
+   real_t vacuum_permitivity          = ((densityFluid_SI * real_c(pow(dx_SI, 3))) * real_c(pow(dx_SI, 3))) /
+                               (real_c(pow(dt_SI, 4)) * Ampere_unit * Ampere_unit);
 
    vacuum_permitivity = vacuum_permitivity * (vacuum_permitivity_SI);
 
    // now for the charge: read from the parameter file
    const real_t elementaryCharge =
-           real_c(1.60217663 * pow(10, -19)); // 1 elementary charge = 1.60217663 * pow(10, -19) coloumbs
+      real_c(1.60217663 * pow(10, -19)); // 1 elementary charge = 1.60217663 * pow(10, -19) coloumbs
    maxCharge_SI           = maxCharge_SI * elementaryCharge;
    minCharge_SI           = minCharge_SI * elementaryCharge;
    const real_t maxCharge = ((maxCharge_SI) * (V * dt_SI * dt_SI)) / (densityFluid_SI * real_c(pow(dx_SI, 5)));
    const real_t minCharge = ((minCharge_SI) * (V * dt_SI * dt_SI)) / (densityFluid_SI * real_c(pow(dx_SI, 5)));
 
-   WALBERLA_LOG_INFO_ON_ROOT("permitivity in LBM units" << " " << vacuum_permitivity);
-   WALBERLA_LOG_INFO_ON_ROOT("Max charge in LBM units" << " " << maxCharge);
-   WALBERLA_LOG_INFO_ON_ROOT("Min charge in LBM units" << " " << minCharge);
+   WALBERLA_LOG_INFO_ON_ROOT("permitivity in LBM units"
+                             << " " << vacuum_permitivity);
+   WALBERLA_LOG_INFO_ON_ROOT("Max charge in LBM units"
+                             << " " << maxCharge);
+   WALBERLA_LOG_INFO_ON_ROOT("Min charge in LBM units"
+                             << " " << minCharge);
 
    // this is just for verification of the SI to LBM conversions is correctly done or not
-   const real_t q_unit = (densityFluid_SI * real_c(pow(dx_SI, 5))) / (V * dt_SI * dt_SI);
-   const real_t epsilon_unit =
-      (real_c(pow(dt_SI, 4)) * Ampere_unit * Ampere_unit) / ((densityFluid_SI * real_c(pow(dx_SI, 3))) * real_c(pow(dx_SI, 3)));
+   const real_t q_unit       = (densityFluid_SI * real_c(pow(dx_SI, 5))) / (V * dt_SI * dt_SI);
+   const real_t epsilon_unit = (real_c(pow(dt_SI, 4)) * Ampere_unit * Ampere_unit) /
+                               ((densityFluid_SI * real_c(pow(dx_SI, 3))) * real_c(pow(dx_SI, 3)));
    const real_t potential_unit = q_unit / (epsilon_unit * dx_SI);
 
    WALBERLA_UNUSED(potential_unit);
@@ -597,8 +610,8 @@ int main(int argc, char** argv)
    // std::cout << "domain size" << " " << domainSize[0] << std::endl;
    // std::cout << "potential at boundary in  Lattice units"<< " "
    // <<(1/(4*3.14*vacuum_permitivity))*(maxCharge/(domainSize[0]/2)) << std::endl; std::cout << "potential at boundary
-   // in  SI units"<< " " << (1/(4*3.14*vacuum_permitivity*epsilon_unit))*((maxCharge*q_unit)/(domainSize[0]/2*dx_SI)) <<
-   // std::endl; std::cout << "potential at boundary in  SI units"<< " " <<
+   // in  SI units"<< " " << (1/(4*3.14*vacuum_permitivity*epsilon_unit))*((maxCharge*q_unit)/(domainSize[0]/2*dx_SI))
+   // << std::endl; std::cout << "potential at boundary in  SI units"<< " " <<
    // (1/(4*3.14*vacuum_permitivity_SI))*((maxCharge_SI)/(domainSize[0]/2*dx_SI)) << std::endl;
 
    // verification ends here
@@ -608,7 +621,6 @@ int main(int argc, char** argv)
    const uint_t vtkSpacingParticles = uint_c(std::ceil(vtkSpacingParticles_SI / dt_SI));
    const uint_t vtkSpacingFluid     = uint_c(std::ceil(vtkSpacingFluid_SI / dt_SI));
 
-
    const real_t poissonsRatio         = real_t(0.22);
    const real_t kappa                 = real_t(2) * (real_t(1) - poissonsRatio) / (real_t(2) - poissonsRatio);
    const real_t particleCollisionTime = 4_r * diameter;
@@ -668,8 +680,13 @@ int main(int argc, char** argv)
 
    auto poissonSolver = PoissonSolver< DAMPED_JACOBI >(/* src */ potentialFieldID, /* dst */ potentialFieldCopyID,
                                                        /* rhs */ chargeDensityFieldID,
-                                                       blocks, boundaryHandling,
-                                                       uint_c(1000), real_c(1e-16), uint_c(1000));
+                                                       /* blockforest */blocks,
+                                                       /* boundary conditions */ boundaryHandling, boundaryConditions,
+                                                       /* iterations */ uint_c(jacobiIterations),
+                                                       /* use abs (true) or rel (false) threshold */ true,
+                                                       /* abs res threshold */ real_c(absResThreshold),
+                                                       /* rel res threshold */ real_c(1e-6),
+                                                       /* res check freq */ uint_c(1000));
    //////////////////
    // RPD COUPLING //
    //////////////////
@@ -796,6 +813,7 @@ int main(int argc, char** argv)
       else if (simulationCase == StokesFlow)
       {
          WALBERLA_LOG_INFO_ON_ROOT("Stokes Flow Validation Test Case Simulation is Chosen");
+         WALBERLA_LOG_INFO_ON_ROOT("Gravitational Force is turned off for this simulation");
          particleLocation = Vector3< real_t >(simulationDomain.center()[0], simulationDomain.center()[1],
                                               2 * simulationDomain.center()[2] - 100);
          inflowVec        = Vector3< real_t >(0_r, 0_r, 0);
@@ -804,6 +822,7 @@ int main(int argc, char** argv)
       else if (simulationCase == moderateReynoldsTerminalVelocity)
       {
          WALBERLA_LOG_INFO_ON_ROOT("Moderate Reynolds Number Velocity Validation Test Case Simulation is Chosen");
+         WALBERLA_LOG_INFO_ON_ROOT("Gravitational Force is turned off for this simulation");
          const real_t startingGapSize_SI = real_t(120e-3) + real_t(0.25) * particleDiameter_SI;
          WALBERLA_LOG_INFO_ON_ROOT("starting gap si"
                                    << " " << startingGapSize_SI);
@@ -830,7 +849,7 @@ int main(int argc, char** argv)
 
    case Showcase: {
       WALBERLA_LOG_INFO_ON_ROOT("Charged Particle Showcase Simulation is Chosen");
-      generationDomain = simulationDomain.createFromMinMaxCorner(0, 0, 0, 200, 200, 800);
+      generationDomain = math::GenericAABB< real_t >::createFromMinMaxCorner(0, 0, 0, 200, 200, 800);
       inflowVec        = Vector3< real_t >(0_r, 0_r, uInflow);
 
       uint_t particleCount = 0;
@@ -858,7 +877,7 @@ int main(int argc, char** argv)
 
       WALBERLA_LOG_INFO_ON_ROOT("Initiating User Defined Simulation");
       if (extendSimulationDomain) { generationDomain = simulationDomain.getExtended(-0.5_r * Spacing); }
-      else { generationDomain = simulationDomain.createFromMinMaxCorner(0, 0, 0, 200, 200, 800); }
+      else { generationDomain = math::GenericAABB< real_t >::createFromMinMaxCorner(0, 0, 0, 200, 200, 800); }
       inflowVec = Vector3< real_t >(0_r, 0_r, uInflow);
 
       for (auto pt : grid_generator::SCGrid(generationDomain, generationDomain.center(), Spacing))
@@ -955,11 +974,12 @@ int main(int argc, char** argv)
          blocks, "particle and volume fraction field",
          std::vector< lbm_mesapd_coupling::psm::ParticleAndVolumeFraction_T >(), field::fzyx, 0);
 
-   auto chargeForceUpdate = ChargeForceUpdate(blocks, potentialFieldID, electrostaticForceFieldID,
-                                  particleAndVolumeFractionFieldID, chargeDensityFieldID, accessor, vacuum_permitivity);
+   auto chargeForceUpdate =
+      ChargeForceUpdate(blocks, potentialFieldID, electrostaticForceFieldID, particleAndVolumeFractionFieldID,
+                        chargeDensityFieldID, accessor, vacuum_permitivity);
 
-   auto chargeDensityUpdate = ChargeDensityUpdate(blocks, particleAndVolumeFractionFieldID, chargeDensityFieldID,
-                                                  accessor, vacuum_permitivity);
+   auto chargeDensityUpdate =
+      ChargeDensityUpdate(blocks, particleAndVolumeFractionFieldID, chargeDensityFieldID, accessor, vacuum_permitivity);
 
    lbm_mesapd_coupling::psm::ParticleAndVolumeFractionMapping particleMapping(
       blocks, accessor, lbm_mesapd_coupling::RegularParticlesSelector(), particleAndVolumeFractionFieldID, 3);
@@ -968,6 +988,10 @@ int main(int argc, char** argv)
    lbm_mesapd_coupling::psm::initializeDomainForPSM< LatticeModel_T, 1 >(*blocks, pdfFieldID,
                                                                          particleAndVolumeFractionFieldID, *accessor);
 
+   // after particle-volume field setup: init RHS of PDE and compute init residual
+   chargeDensityUpdate();
+   poissonSolver.computeInitialResidual();
+
    // setup of the LBM communication for synchronizing the pdf field between neighboring blocks
    blockforest::communication::UniformBufferedScheme< Stencil_T > optimizedPDFCommunicationScheme(blocks);
    optimizedPDFCommunicationScheme.addPackInfo(
@@ -1040,12 +1064,11 @@ int main(int argc, char** argv)
                   particleAndVolumeFractionFieldID);
             ScalarField_T* BField = blockIt->getData< ScalarField_T >(BFieldID);
 
-             WALBERLA_FOR_ALL_CELLS_XYZ(particleAndVolumeFractionField,
-                                        BField->get(x, y, z) = 0.0;
+            WALBERLA_FOR_ALL_CELLS_XYZ(particleAndVolumeFractionField, BField->get(x, y, z) = 0.0;
 
-                                        for (auto &e: particleAndVolumeFractionField->get(x, y, z))
-                                            BField->get(x, y, z) += e.second;
-             )
+                                       for (auto& e
+                                            : particleAndVolumeFractionField->get(x, y, z)) BField->get(x, y, z) +=
+                                       e.second;)
          }
       });
 
@@ -1091,7 +1114,9 @@ int main(int argc, char** argv)
    }
 
    if (vtkSpacingFluid != uint_t(0) || vtkSpacingParticles != uint_t(0))
-   { vtk::writeDomainDecomposition(blocks, "domain_decomposition", vtkFolder); }
+   {
+      vtk::writeDomainDecomposition(blocks, "domain_decomposition", vtkFolder);
+   }
 
    // add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip
    // treatment)
@@ -1110,12 +1135,11 @@ int main(int argc, char** argv)
    ////////////////////////
 
    WcTimingPool timeloopTiming;
-   const bool useOpenMP = false;
-   uint_t numBins = uint_c(simulationDomain.zSize()/(diameter/2));
-   std::vector< real_t >  hydroForceGlobal(numBins, real_t(0));
-   std::vector< real_t >  collisionForceGlobal(numBins, real_t(0));
-   std::vector< real_t >  binCount(numBins, real_t(0));
-
+   const bool useOpenMP  = false;
+   uint_t collisionCount = 0;
+   uint_t tdiff          = 0;
+   Vector3< real_t > collisionStress(real_t(0));
+   Vector3< real_t > hydrodynamicStress(real_t(0));
 
    // time loop
    for (uint_t timeStep = 0; timeStep < numTimeSteps; ++timeStep)
@@ -1138,8 +1162,8 @@ int main(int argc, char** argv)
 
       // Compute electrostatic force field from electric potential (using finite differences)
       chargeForceUpdate();
-
       reduceProperty.operator()< mesa_pd::ElectrostaticForceNotification >(*ps);
+      //WriteElectrostaticForces< ParticleAccessor_T >(accessor, timeStep);
 
       if (timeStep == 0)
       {
@@ -1151,6 +1175,8 @@ int main(int argc, char** argv)
       ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, averageHydrodynamicForceTorque,
                           *accessor);
 
+      if (timeStep >= uint_t(0.75 * real_c(numTimeSteps))) { tdiff += 1; }
+
       for (auto subCycle = uint_t(0); subCycle < numberOfParticleSubCycles; ++subCycle)
       {
          timeloopTiming["RPD"].start();
@@ -1161,6 +1187,42 @@ int main(int argc, char** argv)
          }
          syncCall();
 
+         // collision response
+         ps->forEachParticlePairHalf(
+            useOpenMP, mesa_pd::kernel::ExcludeInfiniteInfinite(), *accessor,
+            [&collisionResponse, &rpdDomain, timeStepSizeRPD, coefficientOfRestitution, particleCollisionTime, kappa,
+             timeStep, numTimeSteps, &collisionCount](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))
+                  {
+                     if (timeStep >= uint_t(0.75 * real_c(numTimeSteps))) { collisionCount += 2; }
+                     auto meff = real_t(1) / (ac.getInvMass(idx1) + ac.getInvMass(idx2));
+                     collisionResponse.setStiffnessAndDamping(0, 0, coefficientOfRestitution, particleCollisionTime,
+                                                              kappa, meff);
+                     collisionResponse(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), acd.getContactNormal(),
+                                       acd.getPenetrationDepth(), timeStepSizeRPD);
+                  }
+               }
+            },
+            *accessor);
+
+         // computing the collision stresses
+         if (timeStep >= uint_t(0.75 * real_c(numTimeSteps)) && subCycle == numberOfParticleSubCycles - 1)
+         {
+            ps->forEachParticle(
+               useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor,
+               [&collisionStress, &diameter](const size_t idx1, auto& ac) {
+                  collisionStress[0] += real_c(fabs(ac.getForce(idx1)[0])) / ((math::pi) * diameter * diameter);
+                  collisionStress[1] += real_c(fabs(ac.getForce(idx1)[1])) / ((math::pi) * diameter * diameter);
+                  collisionStress[2] += real_c(fabs(ac.getForce(idx1)[2])) / ((math::pi) * diameter * diameter);
+               },
+               *accessor);
+         }
+
          if (useLubricationForces)
          {
             // lubrication correction
@@ -1183,28 +1245,6 @@ int main(int argc, char** argv)
                *accessor);
          }
 
-         // collision response
-         ps->forEachParticlePairHalf(
-            useOpenMP, mesa_pd::kernel::ExcludeInfiniteInfinite(), *accessor,
-            [&collisionResponse, &rpdDomain, timeStepSizeRPD, coefficientOfRestitution, particleCollisionTime,
-             kappa](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))
-                  {
-                     auto meff = real_t(1) / (ac.getInvMass(idx1) + ac.getInvMass(idx2));
-                     collisionResponse.setStiffnessAndDamping(0, 0, coefficientOfRestitution, particleCollisionTime,
-                                                              kappa, meff);
-                     collisionResponse(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), acd.getContactNormal(),
-                                       acd.getPenetrationDepth(), timeStepSizeRPD);
-                  }
-               }
-            },
-            *accessor);
-
          reduceAndSwapContactHistory(*ps);
 
          // add hydrodynamic force
@@ -1217,9 +1257,9 @@ int main(int argc, char** argv)
          // add electrostatic force
          AddElectrostaticInteractionKernel addElectrostaticInteraction;
          ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addElectrostaticInteraction,
-                           *accessor);
+                             *accessor);
 
-         if (withoutGravity == false)
+         if (!withoutGravity)
          {
             ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addGravitationalForce, *accessor);
          }
@@ -1233,13 +1273,31 @@ int main(int argc, char** argv)
          syncCall();
 
          timeloopTiming["RPD"].end();
+      } // end sub cycle
+
+      // computing the hydrodynamic stresses
+      if (timeStep >= uint_t(0.75 * real_c(numTimeSteps)))
+      {
+         ps->forEachParticle(
+            useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor,
+            [&hydrodynamicStress, diameter, gravitationalForce](const size_t idx1, auto& ac) {
+               hydrodynamicStress[0] +=
+                  real_c(fabs(ac.getHydrodynamicForce(idx1)[0] + gravitationalForce[0])) / ((math::pi) * diameter * diameter);
+               hydrodynamicStress[1] +=
+                  real_c(fabs(ac.getHydrodynamicForce(idx1)[1] + gravitationalForce[1])) / ((math::pi) * diameter * diameter);
+               hydrodynamicStress[2] +=
+                  real_c(fabs(ac.getHydrodynamicForce(idx1)[2] + gravitationalForce[2])) / ((math::pi) * diameter * diameter);
+            },
+            *accessor);
       }
 
       ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, resetHydrodynamicForceTorque, *accessor);
 
       // TODO: write and add resetElectrostaticForce, see above   ---> completed
       auto particleInfo = evaluateParticleInfo(*accessor);
-      auto fluidInfo = evaluateFluidInfo< BoundaryHandling_T >(blocks, pdfFieldID, boundaryHandlingID);
+      auto fluidInfo    = evaluateFluidInfo< BoundaryHandling_T >(blocks, pdfFieldID, boundaryHandlingID);
+
+      if (timeStep % 1000 == 0) { WriteEnsembledVelocityToFile(timeStep, particleInfo); }
 
       if (timeStep % infoSpacing == 0)
       {
@@ -1260,38 +1318,67 @@ int main(int argc, char** argv)
          }
       }
 
+      if (simulationName == "pivVelocityValidation")
+      {
+         WALBERLA_LOG_INFO_ON_ROOT("sphere position info:"
+                                   << " " << particleInfo.heightOfMass << " "
+                                   << "velocity info is:"
+                                   << " " << particleInfo.averageVelocity << " "
+                                   << "numParticles" << particleInfo.numParticles);
+      }
+
       // Solid Volume Fractions printing //
-      auto solidVolFrac = computeSolidVolumeFraction< StructuredBlockStorage, FlagField_T >(
+      auto solidVolFrac = ComputeSolidVolumeFraction< StructuredBlockStorage, FlagField_T >(
          blocks, flagFieldID, particleAndVolumeFractionFieldID, simulationDomain);
 
-      if (timeStep < 10)
+      if (timeStep > 120000)
       {
-         if (timeStep % 1 == 0)
+         if (timeStep % 1000 == 0)
          {
             solidVolFrac(timeStep);
-            computeParticleProperties< ParticleAccessor_T >(accessor, timeStep);
+            ComputeParticleProperties< ParticleAccessor_T >(accessor, timeStep);
          }
       }
+   } // end time loop
 
-      computeParticleStresses< ParticleAccessor_T >(accessor, hydroForceGlobal, collisionForceGlobal, binCount,
-                                                    gravitationalForce);
-   }
-   walberla::mpi::allReduceInplace(hydroForceGlobal, walberla::mpi::SUM);
-   walberla::mpi::allReduceInplace(binCount, walberla::mpi::SUM);
-   walberla::mpi::allReduceInplace(collisionForceGlobal, walberla::mpi::SUM);
+   walberla::mpi::allReduceInplace(collisionCount, walberla::mpi::SUM);
 
-   for (uint_t i = 0; i < numBins; i++)
+   // collision frequency time and ensembled average:
+   uint_t numParticles = 2000;
+   if (simulationCase == Showcase)
    {
-      if (binCount[i] > 0)
-      {
-         hydroForceGlobal[i]     = (hydroForceGlobal[i] * diameter * diameter) / (densityFluid * viscosity*viscosity);
-         collisionForceGlobal[i] = (collisionForceGlobal[i] * diameter * diameter) / (densityFluid * viscosity*viscosity);
-         hydroForceGlobal[i] /= binCount[i];
-         collisionForceGlobal[i] /= binCount[i];
-      }
+      WALBERLA_LOG_INFO_ON_ROOT("Coliision count is:  " << collisionCount);
+      WALBERLA_LOG_INFO_ON_ROOT("tdiff is:  " << tdiff);
+      WALBERLA_LOG_INFO_ON_ROOT("num particles is:  " << numParticles);
+      WALBERLA_LOG_INFO_ON_ROOT("subcycles is:  " << numberOfParticleSubCycles);
+      real_t collisionfreq = real_c(collisionCount) / real_c((numParticles * tdiff * numberOfParticleSubCycles));
+      WALBERLA_LOG_INFO_ON_ROOT("Coliision Frequency is:  " << collisionfreq);
+
+      // collision stress reduction
+      walberla::mpi::allReduceInplace(collisionStress[0], walberla::mpi::SUM);
+      walberla::mpi::allReduceInplace(collisionStress[1], walberla::mpi::SUM);
+      walberla::mpi::allReduceInplace(collisionStress[2], walberla::mpi::SUM);
+
+      // collision stress average
+      collisionStress[0] /= real_c(numParticles * tdiff);
+      collisionStress[1] /= real_c(numParticles * tdiff);
+      collisionStress[2] /= real_c(numParticles * tdiff);
+      WALBERLA_LOG_INFO_ON_ROOT("Collision stresses average without normalization is:  "
+                                << collisionStress[0] << " " << collisionStress[1] << " " << collisionStress[2]);
+
+      // hydrodynamic stress reduction
+      walberla::mpi::allReduceInplace(hydrodynamicStress[0], walberla::mpi::SUM);
+      walberla::mpi::allReduceInplace(hydrodynamicStress[1], walberla::mpi::SUM);
+      walberla::mpi::allReduceInplace(hydrodynamicStress[2], walberla::mpi::SUM);
+
+      // hydrodynamic stress average
+      hydrodynamicStress[0] /= real_c(numParticles * tdiff);
+      hydrodynamicStress[1] /= real_c(numParticles * tdiff);
+      hydrodynamicStress[2] /= real_c(numParticles * tdiff);
+      WALBERLA_LOG_INFO_ON_ROOT("Hydrodynamic stresses average without normalization is:  "
+                                << hydrodynamicStress[0] << " " << hydrodynamicStress[1] << " "
+                                << hydrodynamicStress[2]);
    }
-   writeStressesToFile(hydroForceGlobal, collisionForceGlobal);
-
    timeloopTiming.logResultOnRoot();
 
    return EXIT_SUCCESS;
diff --git a/apps/showcases/ChargedParticles/ErrorNorms.h b/apps/showcases/ChargedParticles/ErrorNorms.h
new file mode 100644
index 0000000000000000000000000000000000000000..103c79580c948c2ec5bf7d8e40834deca299ec78
--- /dev/null
+++ b/apps/showcases/ChargedParticles/ErrorNorms.h
@@ -0,0 +1,125 @@
+//
+// Created by avnss on 4/21/2023.
+//
+
+#ifndef WALBERLA_ERRORNORMS_H
+#define WALBERLA_ERRORNORMS_H
+
+namespace walberla
+{
+real_t computeErrorMax(const shared_ptr< StructuredBlockForest >& blocks, BlockDataID& NumericalSol,
+                       const math::AABB& domainAABB, const real_t radius, const real_t epsilon, const real_t charge,
+                       const BlockDataID& potentialErrorID)
+{
+   real_t error   = real_c(0);
+   real_t normmax = real_c(0);
+
+   for (auto block = blocks->begin(); block != blocks->end(); ++block)
+   {
+      ScalarField_T* solutionField = block->getData< ScalarField_T >(NumericalSol);
+
+      real_t blockResult            = real_t(0);
+      real_t blockResult_analytical = real_t(0);
+            WALBERLA_FOR_ALL_CELLS_XYZ_OMP(solutionField, omp parallel for schedule(static) reduction(max: blockResult, max: blockResult_analytical),
+                                           const auto cellAABB = blocks->getBlockLocalCellAABB(*block, Cell(x,y,z));
+                                                   auto cellCenter = cellAABB.center();
+
+                                                       real_t posX = cellCenter[0];
+                                                       real_t posY = cellCenter[1];
+                                                       real_t posZ = cellCenter[2];
+
+                                                       real_t analyticalSol;
+
+                                                       // analytical solution of problem with neumann/dirichlet boundaries
+
+                                                       real_t distance = real_c(
+                                                               sqrt(pow((domainAABB.center()[0] - posX), 2) +
+                                                                    pow((domainAABB.center()[1] - posY), 2) +
+                                                                    pow((domainAABB.center()[2] - posZ), 2)));
+                                                       if (distance >= radius) {
+         analyticalSol = (1 / (4 * math::pi * epsilon)) * (charge / distance);
+                                                       } else {
+         analyticalSol = (1 / (4 * math::pi * epsilon)) * (charge / (2 * radius)) * (3 - pow((distance / radius), 2));
+                                                       }
+
+                                                       real_t currErr = real_c(
+                                                               fabs(solutionField->get(x, y, z) - analyticalSol) /
+                                                               analyticalSol);
+                                                       blockResult = std::max(blockResult, currErr);
+                                                       blockResult_analytical = std::max(blockResult_analytical,
+                                                                                         analyticalSol);
+
+
+            )
+            error   = std::max(error, blockResult);
+            normmax = std::max(normmax, blockResult_analytical);
+   }
+
+   mpi::allReduceInplace(error, mpi::MAX);
+   mpi::allReduceInplace(normmax, mpi::MAX);
+
+   return error;
+}
+
+real_t computeErrorL2(const shared_ptr< StructuredBlockForest >& blocks, BlockDataID& NumericalSol,
+                      const math::AABB& domainAABB, const real_t radius, const real_t epsilon, const real_t charge,
+                      const BlockDataID& potentialErrorID)
+{
+   real_t error  = real_c(0);
+   real_t normL2 = real_c(0);
+   real_t cells  = real_t(0);
+
+   for (auto block = blocks->begin(); block != blocks->end(); ++block)
+   {
+            ScalarField_T* solutionField       = block->getData< ScalarField_T >(NumericalSol);
+            ScalarField_T* potentialErrorField = block->getData< ScalarField_T >(potentialErrorID);
+
+            real_t blockResult            = real_t(0);
+            real_t blockResult_analytical = real_t(0);
+            real_t block_local_cells      = real_t(0);
+            WALBERLA_FOR_ALL_CELLS_XYZ_OMP(solutionField, omp parallel for schedule(static) reduction(+: blockResult, +: blockResult_analytical, +: block_local_cells),
+                                           const auto cellAABB = blocks->getBlockLocalCellAABB(*block, Cell(x,y,z));
+                                                   auto cellCenter = cellAABB.center();
+
+                                                       real_t posX = cellCenter[0];
+                                                       real_t posY = cellCenter[1];
+                                                       real_t posZ = cellCenter[2];
+
+                                                       real_t analyticalSol;
+
+                                                       // analytical solution of problem with neumann/dirichlet boundaries
+
+                                                       real_t distance = real_c(
+                                                               sqrt(pow((domainAABB.center()[0] - posX), 2) +
+                                                                    pow((domainAABB.center()[1] - posY), 2) +
+                                                                    pow((domainAABB.center()[2] - posZ), 2)));
+                                                       if (distance >= radius) {
+         analyticalSol = (1 / (4 * math::pi * epsilon)) * (charge / distance);
+                                                       } else {
+         analyticalSol = (1 / (4 * math::pi * epsilon)) * (charge / (2 * radius)) * (3 - pow((distance / radius), 2));
+                                                       }
+
+                                                       //potentialAnalyticalField->get(x, y, z) = analyticalSol;
+                                                       real_t currErr = (solutionField->get(x, y, z) - analyticalSol);
+                                                       potentialErrorField ->get(x, y, z) = abs(currErr);
+                                                       blockResult += currErr * currErr;
+                                                       blockResult_analytical += analyticalSol * analyticalSol;
+                                                       block_local_cells += 1;
+
+
+            )
+            cells += block_local_cells;
+            error += blockResult;
+            normL2 += blockResult_analytical;
+   }
+   mpi::allReduceInplace(error, mpi::SUM);
+   mpi::allReduceInplace(normL2, mpi::SUM);
+   mpi::allReduceInplace(cells, mpi::SUM);
+
+   WALBERLA_LOG_INFO_ON_ROOT("number of cell is:"
+                             << " " << cells);
+   return (std::sqrt(error / (cells))); //(std::sqrt(error/(cells)))/(std::sqrt(normL2/(cells)));
+}
+} // namespace walberla
+
+#endif // WALBERLA_ERRORNORMS_H
diff --git a/apps/showcases/ChargedParticles/postProcessingUtilities.h b/apps/showcases/ChargedParticles/PostProcessingUtilities.h
similarity index 67%
rename from apps/showcases/ChargedParticles/postProcessingUtilities.h
rename to apps/showcases/ChargedParticles/PostProcessingUtilities.h
index 4f831bac659f763c0f064e419335b613ad9cc4f9..c899e0b6328a0a19ee7f8e845fa3bdd612dbc57e 100644
--- a/apps/showcases/ChargedParticles/postProcessingUtilities.h
+++ b/apps/showcases/ChargedParticles/PostProcessingUtilities.h
@@ -1,27 +1,28 @@
-//
-// Created by avnss on 7/13/2023.
-//
-
 #ifndef WALBERLA_POSTPROCESSINGUTILITIES_H
 #define WALBERLA_POSTPROCESSINGUTILITIES_H
 
 #pragma once
-#include <vector>
-#include <math.h>
 #include "core/DataTypes.h"
-#include "lbm_mesapd_coupling/DataTypes.h"
+
 #include "lbm/field/AddToStorage.h"
-#include <filesystem>
+
+#include "lbm_mesapd_coupling/DataTypes.h"
+
+#include <cmath>
 #include <core/mpi/MPITextFile.h>
+#include <vector>
 
-namespace walberla {
+#include "fstream"
+
+namespace walberla
+{
 
 // class to compute the volume fractions average at each unique height along with xy spatial plane
 template< typename BlockStorage_T, typename FlagField_T >
-class computeSolidVolumeFraction
+class ComputeSolidVolumeFraction
 {
  public:
-   computeSolidVolumeFraction(const shared_ptr< BlockStorage_T >& blocks, const BlockDataID& flagFieldID,
+   ComputeSolidVolumeFraction(const shared_ptr< BlockStorage_T >& blocks, const BlockDataID& flagFieldID,
                               const BlockDataID& particleAndVolumeFractionFieldID,
                               const math::GenericAABB< real_t > simulationDomain)
       : blocks_(blocks), flagFieldID_(flagFieldID), particleAndVolumeFractionFieldID_(particleAndVolumeFractionFieldID),
@@ -70,14 +71,14 @@ class computeSolidVolumeFraction
    math::GenericAABB< real_t > simulationDomain_;
    std::vector< real_t > z_SolidVolumeFraction_;
 
-   void printToSeparateFile(std::string filename)
+   void printToSeparateFile(const std::string& filename)
    {
       WALBERLA_ROOT_SECTION()
       {
          const std::string directoryName = "FractionData";
 
          // Create the directory if it doesn't exist
-         if (!std::filesystem::exists(directoryName)) { std::filesystem::create_directory(directoryName); }
+         if (!filesystem::exists(directoryName)) { filesystem::create_directory(directoryName); }
          const std::string filePath = directoryName + "/" + filename;
 
          // Open the file for writing
@@ -105,15 +106,14 @@ class computeSolidVolumeFraction
 
 // class to compute and write the particle x,y,z velocities and the unique ids into a file
 
-template<typename ParticleAccessor_T >
-void computeParticleProperties(const shared_ptr< ParticleAccessor_T >& ac,
-                               uint_t currentTimeStep)
+template< typename ParticleAccessor_T >
+void ComputeParticleProperties(const shared_ptr< ParticleAccessor_T >& ac, uint_t currentTimeStep)
 {
-   std::string filename            = "ParticleProperties_" + std::to_string(currentTimeStep) + ".txt";
+   std::string const filename      = "ParticleProperties_" + std::to_string(currentTimeStep) + ".txt";
    const std::string directoryName = "ParticleData";
 
    // Create the directory if it doesn't exist
-   if (!std::filesystem::exists(directoryName)) { std::filesystem::create_directory(directoryName); }
+   if (!filesystem::exists(directoryName)) { filesystem::create_directory(directoryName); }
    const std::string filePath = directoryName + "/" + filename;
 
    // Open the file for writing using std::ostringstream
@@ -121,7 +121,6 @@ void computeParticleProperties(const shared_ptr< ParticleAccessor_T >& ac,
 
    WALBERLA_ROOT_SECTION()
    {
-
       outputFile << "Uid"
                  << ","
                  << "Velocity_X"
@@ -155,65 +154,99 @@ void computeParticleProperties(const shared_ptr< ParticleAccessor_T >& ac,
    walberla::mpi::writeMPITextFile(filePath, outputFile.str());
 }
 
-
-// class to compute and write the particle x,y,z velocities and the unique ids into a file
-
-template<typename ParticleAccessor_T >
-void computeParticleStresses(const shared_ptr< ParticleAccessor_T >& ac,
-                             std::vector< real_t >& hydroForceGlobal, std::vector< real_t >& collisionForceGlobal,
-                             std::vector< real_t >& binCount, Vector3<real_t> gravitationForce)
+template< typename ParticleAccessor_T >
+void ComputeParticleStresses(const shared_ptr< ParticleAccessor_T >& ac, std::vector< real_t >& hydroForceGlobal,
+                             std::vector< real_t >& collisionForceGlobal, std::vector< real_t >& binCount,
+                             Vector3< real_t > gravitationForce)
 {
-   real_t diameter = 20;
-
+   real_t const diameter = real_c(20);
 
    for (uint_t i = 0; i < ac->size(); ++i)
    {
       if (isSet(ac->getFlags(i), walberla::mesa_pd::data::particle_flags::GHOST)) continue;
       if (isSet(ac->getFlags(i), walberla::mesa_pd::data::particle_flags::GLOBAL)) continue;
 
-      real_t hydroLubricationForce = (ac->getHydrodynamicForce(i) + gravitationForce).length();   // F'h = Fh - Vp(rhof-rhop)*g
-      real_t collisionForce        = (ac->getForce(i) - (ac->getHydrodynamicForce(i) + gravitationForce)).length(); // Total_force - (F'h)
+      real_t hydroLubricationForce =
+         (ac->getHydrodynamicForce(i) + gravitationForce).length(); // F'h = Fh - Vp(rhof-rhop)*g
+      real_t collisionForce =
+         (ac->getForce(i) - ((ac->getHydrodynamicForce(i) + gravitationForce))).length(); // Total_force - (F'h)
 
-      uint_t bin_index = uint_c((ac->getPosition(i)[2]) / (diameter/2));
-      binCount[bin_index] +=1;
+      uint_t bin_index = uint_c((ac->getPosition(i)[2]) / (diameter / 2));
+      binCount[bin_index] += 1;
       hydroForceGlobal[bin_index] += hydroLubricationForce;
       collisionForceGlobal[bin_index] += collisionForce;
-
    }
-
 }
 
-void writeStressesToFile(std::vector<real_t>& hydroStress, std::vector<real_t>& collisionStress){
+template< typename ParticleAccessor_T >
+void ComputeCollisionFrequency(const shared_ptr< ParticleAccessor_T >& ac, uint_t& collisionCount)
+{
+   real_t const diameter = real_c(20);
 
+   for (uint_t i = 0; i < ac->size(); ++i)
+   {
+      for (uint_t j = 0; j < ac->size(); j++)
+      {
+         if (isSet(ac->getFlags(i), walberla::mesa_pd::data::particle_flags::GLOBAL)) continue;
+         if (isSet(ac->getFlags(j), walberla::mesa_pd::data::particle_flags::GLOBAL)) continue;
+         real_t distance = (ac->getPosition(i) - ac->getPosition(j)).length();
+         if (distance <= diameter) { collisionCount++; }
+      }
+   }
+}
 
+void WriteStressesToFile(std::vector< real_t >& hydroStress, std::vector< real_t >& collisionStress)
+{
    WALBERLA_ROOT_SECTION()
    {
-      std::string filename            = "ParticleStresses.txt";
+      std::string const filename      = "ParticleStresses.txt";
       const std::string directoryName = ".";
 
       // Create the directory if it doesn't exist
-      if (!std::filesystem::exists(directoryName)) { std::filesystem::create_directory(directoryName); }
+      if (!filesystem::exists(directoryName)) { filesystem::create_directory(directoryName); }
       const std::string filePath = directoryName + "/" + filename;
 
       // Open the file for writing using std::ostringstream
       std::ofstream outputFile(filePath);
 
-      outputFile << "Uid"
+      outputFile << "BinId"
                  << ","
                  << "Collision_force"
                  << ","
                  << "Hydro_lub_force" << '\n';
 
-      for (size_t i = 0; i < hydroStress.size(); i++) {
+      for (size_t i = 0; i < hydroStress.size(); i++)
+      {
          outputFile << i << "," << collisionStress[i] << "," << hydroStress[i] << '\n';
       }
       outputFile.close();
    }
 }
 
-}
-
+template< typename ParticleInfo >
+void WriteEnsembledVelocityToFile(const uint_t timestep, ParticleInfo info)
+{
+   WALBERLA_ROOT_SECTION()
+   {
+      std::ofstream file("ensembledVelocity.txt", std::ios_base::app);
+      if (file.is_open())
+      {
+         if (file.tellp() == 0)
+         {
+            file << "timestep"
+                 << ", "
+                 << "ensembledVelocity "
+                 << "\n";
+         }
+         file << timestep << "," << info.ensembledAverageVelocityZ << "\n";
 
+         file.close();
+      }
+      else { std::cerr << "Error opening file for writing.\n"; }
+      file.close();
+   }
+}
 
+} // namespace walberla
 
-#endif //WALBERLA_POSTPROCESSINGUTILITIES_H
+#endif // WALBERLA_POSTPROCESSINGUTILITIES_H
diff --git a/apps/showcases/ChargedParticles/parameterfiles/ForceValidationSingleParticle.prm b/apps/showcases/ChargedParticles/parameterfiles/ForceValidationSingleParticle.prm
index e0f666c2ab5e08f423ca66e0357ebac5285fec87..dec54b1d6a1bdf615140a9bcb113da6b94ba23cc 100644
--- a/apps/showcases/ChargedParticles/parameterfiles/ForceValidationSingleParticle.prm
+++ b/apps/showcases/ChargedParticles/parameterfiles/ForceValidationSingleParticle.prm
@@ -35,6 +35,8 @@ NumericalSetup
     numXBlocks 2;
     numYBlocks 2;
     numZBlocks 2;
+    absResThreshold 1e-19;
+    jacobiIterations 1000;
 
     useLubricationForces true;
     numberOfParticleSubCycles 10;
diff --git a/apps/showcases/ChargedParticles/parameterfiles/singleParticle.prm b/apps/showcases/ChargedParticles/parameterfiles/singleParticle.prm
index ebd2f7320da5ce9a6ae32c9882434db8fb65725c..5a8bbb10d819a9386ab2a2baf6f01a770817f989 100644
--- a/apps/showcases/ChargedParticles/parameterfiles/singleParticle.prm
+++ b/apps/showcases/ChargedParticles/parameterfiles/singleParticle.prm
@@ -35,6 +35,8 @@ NumericalSetup
     numXBlocks 2;
     numYBlocks 2;
     numZBlocks 2;
+    absResThreshold 1e-19;
+    jacobiIterations 1000;
 
     useLubricationForces true;
     numberOfParticleSubCycles 10;
diff --git a/apps/showcases/ChargedParticles/parameterfiles/testIterativeSolvers.prm b/apps/showcases/ChargedParticles/parameterfiles/testIterativeSolvers.prm
index 3c6ead2123b9c7104c7081079952de47767e5c9b..be3631e739f3236995e8aa0f5e12acd683f32c02 100644
--- a/apps/showcases/ChargedParticles/parameterfiles/testIterativeSolvers.prm
+++ b/apps/showcases/ChargedParticles/parameterfiles/testIterativeSolvers.prm
@@ -4,6 +4,7 @@ Setup
     numBlocks    < 2, 2, 2 >;
 
     maxIterations 10000;
+    useAbsResThreshold 1;
     resThreshold 1E-16;
     resCheckFreq 1000;
 }
diff --git a/apps/showcases/ChargedParticles/poisson_solver/DirichletDomainBoundary.h b/apps/showcases/ChargedParticles/poisson_solver/DirichletDomainBoundary.h
index 23de72cd4a39511fd6e888ce0447b5b3e277fe3d..9156d6c1bdc8da47ef46d8b9f986c1b4b291f638 100644
--- a/apps/showcases/ChargedParticles/poisson_solver/DirichletDomainBoundary.h
+++ b/apps/showcases/ChargedParticles/poisson_solver/DirichletDomainBoundary.h
@@ -2,6 +2,7 @@
 #define WALBERLA_DIRICHLETDOMAINBOUNDARY_H
 
 #include "BoundaryCondition.h"
+#include "stencil/D3Q6.h"
 
 namespace walberla
 {
diff --git a/apps/showcases/ChargedParticles/poisson_solver/ParallelCGFixedStencilIteration.h b/apps/showcases/ChargedParticles/poisson_solver/ParallelCGFixedStencilIteration.h
new file mode 100644
index 0000000000000000000000000000000000000000..cad0983b0265682ab5cca90744e32003d0f4566c
--- /dev/null
+++ b/apps/showcases/ChargedParticles/poisson_solver/ParallelCGFixedStencilIteration.h
@@ -0,0 +1,362 @@
+//======================================================================================================================
+//
+//  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 ParallelCGFixedStencilIteration.h
+//! \ingroup pde
+//! \author Richard Angersbach <richard.angersbach@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/Set.h"
+#include "core/logging/Logging.h"
+#include "core/mpi/Reduce.h"
+#include "core/uid/SUID.h"
+
+#include "domain_decomposition/BlockStorage.h"
+
+#include "field/GhostLayerField.h"
+#include "field/iterators/IteratorMacros.h"
+
+#include <functional>
+
+namespace walberla
+{
+namespace pde
+{
+
+template< typename Stencil_T >
+class ParallelCGFixedStencilIteration
+{
+ public:
+   typedef GhostLayerField< real_t, 1 > Field_T;
+
+   ParallelCGFixedStencilIteration(
+      BlockStorage& blocks, const BlockDataID& uId, const BlockDataID& rId, const BlockDataID& dId,
+      const BlockDataID& zId, const BlockDataID& fId, const std::vector< real_t >& weights, const uint_t iterations,
+      const std::function< void() >& synchronizeU, const std::function< void() >& synchronizeD,
+      const std::function< void() >& applyBoundaryHandlingU, const std::function< void() >& applyBoundaryHandlingR,
+      const std::function< void() >& applyBoundaryHandlingD, const real_t residualNormThreshold = real_t(0),
+      const Set< SUID >& requiredSelectors     = Set< SUID >::emptySet(),
+      const Set< SUID >& incompatibleSelectors = Set< SUID >::emptySet());
+
+   void operator()();
+
+ protected:
+   ////////////////////////////
+   // building blocks for CG //
+   ////////////////////////////
+   void calcR();                     // r = f - Au
+   real_t scalarProductRR();         // r*r
+   void copyRToD();                  // d = r
+   void calcAd();                    // z = Ad
+   real_t scalarProductDZ();         // d*z
+   void updateU(const real_t alpha); // u = u + alpha * d
+   void updateR(const real_t alpha); // r = r - alpha * z
+   void updateD(const real_t beta);  // d = r + beta * d
+
+   BlockStorage& blocks_;
+
+   const BlockDataID uId_;
+   const BlockDataID rId_;
+   const BlockDataID dId_;
+   const BlockDataID zId_;
+   const BlockDataID fId_;
+
+   real_t cells_;
+
+   real_t w_[Stencil_T::Size];
+
+   uint_t iterations_;
+   real_t residualNormThreshold_;
+
+   std::function< void() > synchronizeU_;
+   std::function< void() > synchronizeD_;
+
+   std::function< void() > applyBoundaryHandlingU_;
+   std::function< void() > applyBoundaryHandlingR_;
+   std::function< void() > applyBoundaryHandlingD_;
+
+   Set< SUID > requiredSelectors_;
+   Set< SUID > incompatibleSelectors_;
+};
+
+template< typename Stencil_T >
+ParallelCGFixedStencilIteration< Stencil_T >::ParallelCGFixedStencilIteration(
+   BlockStorage& blocks, const BlockDataID& uId, const BlockDataID& rId, const BlockDataID& dId, const BlockDataID& zId,
+   const BlockDataID& fId, const std::vector< real_t >& weights, const uint_t iterations,
+   const std::function< void() >& synchronizeU, const std::function< void() >& synchronizeD,
+   const std::function< void() >& applyBoundaryHandlingU, const std::function< void() >& applyBoundaryHandlingR,
+   const std::function< void() >& applyBoundaryHandlingD, const real_t residualNormThreshold,
+   const Set< SUID >& requiredSelectors, const Set< SUID >& incompatibleSelectors)
+   : blocks_(blocks), uId_(uId), rId_(rId), dId_(dId), zId_(zId), fId_(fId), iterations_(iterations),
+     residualNormThreshold_(residualNormThreshold), synchronizeU_(synchronizeU), synchronizeD_(synchronizeD),
+     applyBoundaryHandlingU_(applyBoundaryHandlingU), applyBoundaryHandlingR_(applyBoundaryHandlingR),
+     applyBoundaryHandlingD_(applyBoundaryHandlingD), requiredSelectors_(requiredSelectors),
+     incompatibleSelectors_(incompatibleSelectors)
+{
+   WALBERLA_ASSERT_EQUAL(weights.size(), Stencil_T::Size);
+   for (uint_t i = uint_t(0); i < Stencil_T::Size; ++i)
+      w_[i] = weights[i];
+
+   uint_t cells(uint_t(0));
+
+   for (auto block = blocks_.begin(requiredSelectors_, incompatibleSelectors_); block != blocks_.end(); ++block)
+   {
+      const Field_T* const u = block->template getData< const Field_T >(uId_);
+      cells += u->xyzSize().numCells();
+   }
+
+   cells_ = real_c(cells);
+   mpi::allReduceInplace(cells_, mpi::SUM);
+}
+
+template< typename Stencil_T >
+void ParallelCGFixedStencilIteration< Stencil_T >::operator()()
+{
+   WALBERLA_LOG_PROGRESS_ON_ROOT("Starting CG iteration with a maximum number of " << iterations_ << " iterations");
+
+   calcR(); // r = f - Au
+
+   real_t rr0          = scalarProductRR(); // r*r
+   real_t residualNorm = std::sqrt(rr0 / cells_);
+
+   if (residualNorm >= residualNormThreshold_)
+   {
+      copyRToD(); // d = r
+
+      uint_t i(uint_t(0));
+      while (i < iterations_)
+      {
+         synchronizeD_();
+
+         calcAd(); // z = Ad
+
+         const real_t alpha = rr0 / scalarProductDZ(); // alpha = r*r / d*z
+         updateU(alpha);                               // u = u + alpha * d
+         updateR(alpha);                               // r = r - alpha * z
+
+         const real_t rr1 = scalarProductRR();
+         residualNorm     = std::sqrt(rr1 / cells_);
+
+         WALBERLA_LOG_PROGRESS_ON_ROOT("Residual norm after CG iteration " << i << " is: " << residualNorm)
+
+         if (residualNorm < residualNormThreshold_)
+         {
+            WALBERLA_LOG_PROGRESS_ON_ROOT("Aborting CG iteration (residual norm threshold reached):"
+                                          "\n  residual norm threshold: "
+                                          << residualNormThreshold_ << "\n  residual norm:           " << residualNorm);
+            break;
+         }
+
+         const real_t beta = rr1 / rr0; // beta = r*r (current) / r*r (previous)
+         updateD(beta);                 // d = r + beta * d
+
+         rr0 = rr1;
+
+         ++i;
+      }
+
+      WALBERLA_LOG_PROGRESS_ON_ROOT("CG iteration finished after " << i << " iterations");
+   }
+   else
+   {
+      WALBERLA_LOG_PROGRESS_ON_ROOT("Aborting CG without a single iteration (residual norm threshold already reached):"
+                                    "\n  residual norm threshold: "
+                                    << residualNormThreshold_ << "\n  residual norm:           " << residualNorm);
+   }
+}
+
+template< typename Stencil_T >
+void ParallelCGFixedStencilIteration< Stencil_T >::calcR() // r = f - Au
+{
+   synchronizeU_();
+
+   for (auto block = blocks_.begin(requiredSelectors_, incompatibleSelectors_); block != blocks_.end(); ++block)
+   {
+      Field_T* rf = block->template getData< Field_T >(rId_);
+      Field_T* ff = block->template getData< Field_T >(fId_);
+      Field_T* uf = block->template getData< Field_T >(uId_);
+
+      WALBERLA_ASSERT_NOT_NULLPTR(rf);
+      WALBERLA_ASSERT_NOT_NULLPTR(ff);
+      WALBERLA_ASSERT_NOT_NULLPTR(uf);
+
+      WALBERLA_ASSERT_EQUAL(rf->xyzSize(), ff->xyzSize());
+      WALBERLA_ASSERT_EQUAL(rf->xyzSize(), uf->xyzSize());
+
+      WALBERLA_ASSERT_GREATER_EQUAL(uf->nrOfGhostLayers(), 1);
+
+      WALBERLA_FOR_ALL_CELLS_XYZ(uf, rf->get(x, y, z) = ff->get(x, y, z);
+
+                                 for (auto dir = Stencil_T::begin(); dir != Stencil_T::end(); ++dir) rf->get(x, y, z) -=
+                                 w_[dir.toIdx()] * uf->getNeighbor(x, y, z, *dir);)
+   }
+
+   applyBoundaryHandlingR_();
+}
+
+template< typename Stencil_T >
+real_t ParallelCGFixedStencilIteration< Stencil_T >::scalarProductRR() // r*r
+{
+   real_t result(real_t(0));
+
+   for (auto block = blocks_.begin(requiredSelectors_, incompatibleSelectors_); block != blocks_.end(); ++block)
+   {
+      Field_T* rf = block->template getData< Field_T >(rId_);
+
+      real_t blockResult(real_t(0));
+
+      WALBERLA_FOR_ALL_CELLS_XYZ_OMP( rf, omp parallel for schedule(static) reduction(+:blockResult),
+         const real_t v = rf->get(x,y,z);
+         blockResult += v * v;
+      )
+
+      result += blockResult;
+   }
+
+   mpi::allReduceInplace(result, mpi::SUM);
+   return result;
+}
+
+template< typename Stencil_T >
+void ParallelCGFixedStencilIteration< Stencil_T >::copyRToD() // d = r
+{
+   for (auto block = blocks_.begin(requiredSelectors_, incompatibleSelectors_); block != blocks_.end(); ++block)
+   {
+      Field_T* rf = block->template getData< Field_T >(rId_);
+      Field_T* df = block->template getData< Field_T >(dId_);
+
+      WALBERLA_ASSERT_NOT_NULLPTR(rf);
+      WALBERLA_ASSERT_NOT_NULLPTR(df);
+
+      WALBERLA_ASSERT_EQUAL(rf->xyzSize(), df->xyzSize());
+
+      WALBERLA_FOR_ALL_CELLS_XYZ(rf, df->get(x, y, z) = rf->get(x, y, z);)
+   }
+
+   applyBoundaryHandlingD_();
+}
+
+template< typename Stencil_T >
+void ParallelCGFixedStencilIteration< Stencil_T >::calcAd() // z = Ad
+{
+   for (auto block = blocks_.begin(requiredSelectors_, incompatibleSelectors_); block != blocks_.end(); ++block)
+   {
+      Field_T* zf = block->template getData< Field_T >(zId_);
+      Field_T* df = block->template getData< Field_T >(dId_);
+
+      WALBERLA_ASSERT_NOT_NULLPTR(zf);
+      WALBERLA_ASSERT_NOT_NULLPTR(df);
+
+      WALBERLA_ASSERT_EQUAL(zf->xyzSize(), df->xyzSize());
+
+      WALBERLA_ASSERT_GREATER_EQUAL(df->nrOfGhostLayers(), 1);
+
+      WALBERLA_FOR_ALL_CELLS_XYZ(df, zf->get(x, y, z) = w_[Stencil_T::idx[stencil::C]] * df->get(x, y, z);
+
+                                 for (auto dir = Stencil_T::beginNoCenter(); dir != Stencil_T::end(); ++dir)
+                                    zf->get(x, y, z) += w_[dir.toIdx()] * df->getNeighbor(x, y, z, *dir);)
+   }
+}
+
+template< typename Stencil_T >
+real_t ParallelCGFixedStencilIteration< Stencil_T >::scalarProductDZ() // d*z
+{
+   real_t result(real_t(0));
+
+   for (auto block = blocks_.begin(requiredSelectors_, incompatibleSelectors_); block != blocks_.end(); ++block)
+   {
+      Field_T* df = block->template getData< Field_T >(dId_);
+      Field_T* zf = block->template getData< Field_T >(zId_);
+
+      WALBERLA_ASSERT_NOT_NULLPTR(df);
+      WALBERLA_ASSERT_NOT_NULLPTR(zf);
+
+      WALBERLA_ASSERT_EQUAL(df->xyzSize(), zf->xyzSize());
+
+      real_t blockResult(real_t(0));
+
+      WALBERLA_FOR_ALL_CELLS_XYZ_OMP( df, omp parallel for schedule(static) reduction(+:blockResult),
+         blockResult += df->get(x,y,z) * zf->get(x,y,z);
+      )
+
+      result += blockResult;
+   }
+
+   mpi::allReduceInplace(result, mpi::SUM);
+   return result;
+}
+
+template< typename Stencil_T >
+void ParallelCGFixedStencilIteration< Stencil_T >::updateU(const real_t alpha) // u = u + alpha * d
+{
+   for (auto block = blocks_.begin(requiredSelectors_, incompatibleSelectors_); block != blocks_.end(); ++block)
+   {
+      Field_T* uf = block->template getData< Field_T >(uId_);
+      Field_T* df = block->template getData< Field_T >(dId_);
+
+      WALBERLA_ASSERT_NOT_NULLPTR(uf);
+      WALBERLA_ASSERT_NOT_NULLPTR(df);
+
+      WALBERLA_ASSERT_EQUAL(uf->xyzSize(), df->xyzSize());
+
+      WALBERLA_FOR_ALL_CELLS_XYZ(uf, uf->get(x, y, z) = uf->get(x, y, z) + alpha * df->get(x, y, z);)
+   }
+
+   applyBoundaryHandlingU_();
+}
+
+template< typename Stencil_T >
+void ParallelCGFixedStencilIteration< Stencil_T >::updateR(const real_t alpha) // r = r - alpha * z
+{
+   for (auto block = blocks_.begin(requiredSelectors_, incompatibleSelectors_); block != blocks_.end(); ++block)
+   {
+      Field_T* rf = block->template getData< Field_T >(rId_);
+      Field_T* zf = block->template getData< Field_T >(zId_);
+
+      WALBERLA_ASSERT_NOT_NULLPTR(rf);
+      WALBERLA_ASSERT_NOT_NULLPTR(zf);
+
+      WALBERLA_ASSERT_EQUAL(rf->xyzSize(), zf->xyzSize());
+
+      WALBERLA_FOR_ALL_CELLS_XYZ(rf, rf->get(x, y, z) = rf->get(x, y, z) - alpha * zf->get(x, y, z);)
+   }
+
+   applyBoundaryHandlingR_();
+}
+
+template< typename Stencil_T >
+void ParallelCGFixedStencilIteration< Stencil_T >::updateD(const real_t beta) // d = r + beta * d
+{
+   for (auto block = blocks_.begin(requiredSelectors_, incompatibleSelectors_); block != blocks_.end(); ++block)
+   {
+      Field_T* df = block->template getData< Field_T >(dId_);
+      Field_T* rf = block->template getData< Field_T >(rId_);
+
+      WALBERLA_ASSERT_NOT_NULLPTR(df);
+      WALBERLA_ASSERT_NOT_NULLPTR(rf);
+
+      WALBERLA_ASSERT_EQUAL(df->xyzSize(), rf->xyzSize());
+
+      WALBERLA_FOR_ALL_CELLS_XYZ(df, df->get(x, y, z) = rf->get(x, y, z) + beta * df->get(x, y, z);)
+   }
+
+   applyBoundaryHandlingD_();
+}
+
+} // namespace pde
+} // namespace walberla
diff --git a/apps/showcases/ChargedParticles/poisson_solver/PoissonSolver.h b/apps/showcases/ChargedParticles/poisson_solver/PoissonSolver.h
index e0c1d5906a0db8fd52ada959cba00cff3437df2a..eaf07d1b33092f89f9d00b2fcf323c65673e01e1 100644
--- a/apps/showcases/ChargedParticles/poisson_solver/PoissonSolver.h
+++ b/apps/showcases/ChargedParticles/poisson_solver/PoissonSolver.h
@@ -3,11 +3,14 @@
 
 #include "pde/all.h"
 
+#include <algorithm>
+
 #include "CustomBoundary.h"
 #include "DirichletDomainBoundary.h"
 #include "Neumann.h"
+#include "ParallelCGFixedStencilIteration.h"
 
-enum Enum { WALBERLA_JACOBI, WALBERLA_SOR, DAMPED_JACOBI };
+enum Enum { WALBERLA_JACOBI, WALBERLA_SOR, WALBERLA_CG, DAMPED_JACOBI };
 
 namespace walberla
 {
@@ -44,9 +47,13 @@ class PoissonSolver
    PoissonSolver(const BlockDataID& src, const BlockDataID& dst, const BlockDataID& rhs,
                  const std::shared_ptr< StructuredBlockForest >& blocks,
                  const std::function< void() >& boundaryHandling,
-                 uint_t iterations = uint_t(1000), real_t residualNormThreshold = real_c(1e-4), uint_t residualCheckFrequency = uint_t(100))
-      : src_(src), dst_(dst), rhs_(rhs), blocks_(blocks), boundaryHandling_(boundaryHandling) {
-
+                 const std::vector< BoundaryCondition >& boundaryConditions, uint_t iterations = uint_t(1000),
+                 bool useAbsResNormThres = false, real_t absResNormThres = real_c(1e-10),
+                 real_t relResNormThres = real_c(1e-4), uint_t resCheckFreq = uint_t(100))
+      : src_(src), dst_(dst), rhs_(rhs), blocks_(blocks), boundaryHandling_(boundaryHandling),
+        boundaryConditions_(boundaryConditions), useRelativeResidualThreshold_(!useAbsResNormThres),
+        relativeResidualReductionFactor_(relResNormThres)
+   {
       // 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()) +
@@ -65,9 +72,33 @@ class PoissonSolver
       commScheme_->addPackInfo(make_shared< field::communication::PackInfo< ScalarField_T > >(src_));
       commScheme_->addPackInfo(make_shared< field::communication::PackInfo< ScalarField_T > >(rhs_));
 
-      // res norm
-
-      residualNorm_ = make_shared< pde::ResidualNorm< Stencil_T > >(blocks_->getBlockStorage(), src_, rhs_, laplaceWeights_);
+      // res norm computation callback
+
+      residualNorm_ =
+         make_shared< pde::ResidualNorm< Stencil_T > >(blocks_->getBlockStorage(), src_, rhs_, laplaceWeights_);
+
+      // handling for absolute/relative residual thresholds
+      // if absolute res thres is used -> simply propagate ctor values for absolute res checks to Jacobi iteration
+      // else
+      //  -> only execute "resCheckFreq" iterations at once and then execute logic for relative exit crit
+      //  -> repeat "numExecutions" times to perform the same amount of total iterations
+      WALBERLA_ASSERT(iterations % resCheckFreq == 0,
+                      "Number of iterations should be divisible by residual check frequency!")
+
+      if (useAbsResNormThres)
+      {
+         residualNormThreshold_  = absResNormThres;
+         residualCheckFrequency_ = resCheckFreq;
+         numIterPerExecution_    = iterations;
+         numExecutions_          = uint_c(1);
+      }
+      else
+      {
+         residualNormThreshold_  = real_c(0);
+         residualCheckFrequency_ = uint_c(0);
+         numIterPerExecution_    = resCheckFreq;
+         numExecutions_          = iterations / resCheckFreq;
+      }
 
       // jacobi
 
@@ -75,14 +106,17 @@ class PoissonSolver
 
       // use custom impl with damping or jacobi from waLBerla
       std::function< void(IBlock*) > jacSweep = {};
-      if (solver == DAMPED_JACOBI) {
+      if (solver == DAMPED_JACOBI)
+      {
          jacSweep = [this](IBlock* block) { dampedJacobiSweep(block); };
-      } else {
+      }
+      else if (solver == WALBERLA_JACOBI) {
          jacSweep = *jacobiFixedSweep_;
       }
 
-      jacobiIteration_ = std::make_unique< pde::JacobiIteration >(blocks_->getBlockStorage(), iterations, *commScheme_,
-                                                                  jacSweep, *residualNorm_, residualNormThreshold, residualCheckFrequency);
+      jacobiIteration_ = std::make_unique< pde::JacobiIteration >(blocks_->getBlockStorage(), numIterPerExecution_,
+                                                                  *commScheme_, jacSweep, *residualNorm_,
+                                                                  residualNormThreshold_, residualCheckFrequency_);
 
       jacobiIteration_->addBoundaryHandling(boundaryHandling_);
 
@@ -93,47 +127,138 @@ class PoissonSolver
       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);
+         blocks_->getBlockStorage(), numIterPerExecution_, *commScheme_, sorFixedSweep_->getRedSweep(),
+         sorFixedSweep_->getBlackSweep(), *residualNorm_, residualNormThreshold_, residualCheckFrequency_);
 
       sorIteration_->addBoundaryHandling(boundaryHandling);
 
+      // CG
+
+      d_ = field::addToStorage< ScalarField_T >(blocks, "d", real_t(0), field::fzyx, uint_t(1));
+      r_ = field::addToStorage< ScalarField_T >(blocks, "r", real_t(0), field::fzyx, uint_t(1));
+      z_ = field::addToStorage< ScalarField_T >(blocks, "z", real_t(0), field::fzyx, uint_t(1));
+
+      syncD_ = make_shared< blockforest::communication::UniformBufferedScheme< Stencil_T > >(blocks_);
+      syncD_->addPackInfo(make_shared< field::communication::PackInfo< ScalarField_T > >(d_));
+
+      // zero-value dirichlet/neumann BCs for CG fields: r and d
+      std::vector< BoundaryCondition > cgBoundaryConditions;
+      for (auto &cond : boundaryConditions_)
+         cgBoundaryConditions.emplace_back(cond.getDirection(), cond.getType(), 0_r);
+
+      applyBoundaryHandlingR_ = CustomBoundary< ScalarField_T >(*blocks, r_, cgBoundaryConditions);
+      applyBoundaryHandlingD_ = CustomBoundary< ScalarField_T >(*blocks, d_, cgBoundaryConditions);
+
+      cgIteration_ = std::make_unique< pde::ParallelCGFixedStencilIteration< Stencil_T > >(
+         blocks->getBlockStorage(), src_, r_, d_, z_, rhs_, laplaceWeights_, uint_t(numIterPerExecution_), *commScheme_,
+         *syncD_, boundaryHandling_, applyBoundaryHandlingR_, applyBoundaryHandlingD_, residualNormThreshold_);
+   }
+
+   // get approximate solution
+   void operator()()
+   {
+      for (uint_t executions = 0; executions < numExecutions_; ++executions)
+      {
+         // execute solver...
+         switch (solver)
+         {
+         case WALBERLA_CG:
+            (*cgIteration_)();
+            break;
+         case WALBERLA_SOR:
+            (*sorIteration_)();
+            break;
+         default:
+            (*jacobiIteration_)();
+            break;
+         }
+
+         // .. and check if (relative) res threshold was reached
+         if (useRelativeResidualThreshold_)
+         {
+            auto curRes = (*residualNorm_)();
+
+            WALBERLA_LOG_INFO_ON_ROOT("Residual norm after " << executions * numIterPerExecution_
+                                                             << " Jacobi iterations: " << curRes);
+
+            if (curRes <= relativeResidualReductionFactor_ * initRes_)
+            {
+               WALBERLA_LOG_INFO_ON_ROOT("Aborting Jacobi iteration (residual norm threshold reached):"
+                                         "\n  residual norm thres: "
+                                         << relativeResidualReductionFactor_ * initRes_
+                                         << "\n  residual norm:           " << curRes);
+
+               break;
+            }
+         }
+      }
+
+      // print residual after solving
+      boundaryHandling_();
+      (*commScheme_)();
+      auto r = (*residualNorm_)();
+
+      WALBERLA_LOG_INFO_ON_ROOT("Residual after solving = " << r);
+   }
+
+   void computeInitialResidual()
+   {
       // compute initial residual and print
       boundaryHandling_();
       (*commScheme_)();
       initRes_ = (*residualNorm_)();
+
       WALBERLA_LOG_INFO_ON_ROOT("Initial residual = " << initRes_);
    }
-
-   // get approximate solution of electric potential
-   void operator()()
-   {
-      if constexpr (solver != WALBERLA_SOR) {
-         (*jacobiIteration_)();
-      } else {
-         (*sorIteration_)();
-      }
+   real_t getResidualNorm(){
+      return (*residualNorm_)();
    }
 
  private:
+   // input fields
    BlockDataID src_;
    BlockDataID dst_;
    BlockDataID rhs_;
+
+   // specialized CG members
+   BlockDataID d_;
+   BlockDataID r_;
+   BlockDataID z_;
+
+   std::shared_ptr< blockforest::communication::UniformBufferedScheme< Stencil_T > > syncD_;
+
+   std::function< void() > applyBoundaryHandlingR_;
+   std::function< void() > applyBoundaryHandlingD_;
+
+   // general solver members
    std::vector< real_t > laplaceWeights_;
    std::shared_ptr< StructuredBlockForest > blocks_;
    std::shared_ptr< blockforest::communication::UniformBufferedScheme< Stencil_T > > commScheme_;
 
    std::function< void() > boundaryHandling_;
+   std::vector< BoundaryCondition > boundaryConditions_;
 
    std::shared_ptr< pde::ResidualNorm< Stencil_T > > residualNorm_;
 
+   // iteration schemes
    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_;
 
+   std::unique_ptr< pde::ParallelCGFixedStencilIteration< Stencil_T > > cgIteration_;
+
+   // general residual variables
+   real_t residualNormThreshold_;
+   uint_t residualCheckFrequency_;
+   uint_t numExecutions_;
+   uint_t numIterPerExecution_;
+
+   // variables for relative residual threshold
    real_t initRes_;
+   bool useRelativeResidualThreshold_;
+   real_t relativeResidualReductionFactor_;
 };
 
 } // namespace walberla
diff --git a/apps/showcases/ChargedParticles/poisson_solver/PotentialValidation.cpp b/apps/showcases/ChargedParticles/poisson_solver/PotentialValidation.cpp
index 45751abc179c7aea5a923c878b7f281b3986bf25..0466bb322ca57c099be80d1613385c265be17cff 100644
--- a/apps/showcases/ChargedParticles/poisson_solver/PotentialValidation.cpp
+++ b/apps/showcases/ChargedParticles/poisson_solver/PotentialValidation.cpp
@@ -88,6 +88,7 @@ void solveElectrostaticPoisson(const shared_ptr< StructuredBlockForest >& blocks
                                BlockDataID& analytical)
 {
    auto numIter = uint_c(100000);
+   auto useAbsResThreshold = true;
    auto resThres = real_c(1e-14);
    auto resCheckFreq = uint_c(1000);
 
@@ -103,6 +104,10 @@ void solveElectrostaticPoisson(const shared_ptr< StructuredBlockForest >& blocks
    geometry::Sphere const sphere(position, R_L);
 
    // set dirichlet function per domain face
+   std::vector< BoundaryCondition > dirichletBoundaryConditions;
+   for (const auto& cond : stencil::D3Q6::dir)
+      dirichletBoundaryConditions.emplace_back(cond, "Dirichlet", 0_r); // dummy value, unused
+
    auto dirichletFunction = DirichletFunctionDomainBoundary< ScalarField_T >(*blocks, solution);
 
 #define GET_POTENTIAL_BCS_LAMBDA(dir) \
@@ -119,9 +124,11 @@ void solveElectrostaticPoisson(const shared_ptr< StructuredBlockForest >& blocks
    dirichletFunction.setFunction(stencil::B, GET_POTENTIAL_BCS_LAMBDA(stencil::B));
    dirichletFunction.setFunction(stencil::T, GET_POTENTIAL_BCS_LAMBDA(stencil::T));
 
-   auto poissonSolverSOR = PoissonSolver< WALBERLA_SOR > (solution, solutionCpy, rhs, blocks, dirichletFunction, numIter, resThres, resCheckFreq);
+   auto poissonSolverSOR = PoissonSolver< WALBERLA_CG > (solution, solutionCpy, rhs, blocks,
+                                                        dirichletFunction, dirichletBoundaryConditions,
+                                                        numIter, useAbsResThreshold, resThres, resThres, resCheckFreq);
 
-   // init rhs with two charged particles
+   // init rhs with one charged particle
 
    for (auto block = blocks->begin(); block != blocks->end(); ++block) {
       ScalarField_T* rhsField = block->getData< ScalarField_T >(rhs);
@@ -152,6 +159,8 @@ void solveElectrostaticPoisson(const shared_ptr< StructuredBlockForest >& blocks
          })
    }
 
+   poissonSolverSOR.computeInitialResidual();
+
    for (auto block = blocks->begin(); block != blocks->end(); ++block) {
       ScalarField_T* analyticalField = block->getData< ScalarField_T >(analytical);
 
@@ -201,8 +210,7 @@ void solveElectrostaticPoisson(const shared_ptr< StructuredBlockForest >& blocks
       }
       mpi::allReduceInplace( error, mpi::SUM );
 
-      return std::sqrt( error / real_c(cells));
-
+      return std::sqrt(error / real_c(cells));
    };
 
    // solve with SOR
@@ -246,7 +254,7 @@ int main(int argc, char** argv)
    BlockDataID solutionCpy = field::addCloneToStorage< ScalarField_T >(blocks, solution, "solution (copy)");
    BlockDataID analytical = field::addCloneToStorage< ScalarField_T >(blocks, solution, "analytical");
 
-   solveElectrostaticPoisson (blocks, dx, radiusScaleFactor, domainAABB, useOverlapFraction, solution, solutionCpy, rhs, analytical);
+   solveElectrostaticPoisson(blocks, dx, radiusScaleFactor, domainAABB, useOverlapFraction, solution, solutionCpy, rhs, analytical);
 
    auto vtkWriter = vtk::createVTKOutput_BlockData(*blocks, "block_data", uint_c(1), uint_c(0), true, "VTK" );
    vtkWriter->addCellDataWriter(make_shared<field::VTKWriter<ScalarField_T> >(solution, "solution"));
diff --git a/apps/showcases/ChargedParticles/poisson_solver/TestIterativeSolvers.cpp b/apps/showcases/ChargedParticles/poisson_solver/TestIterativeSolvers.cpp
index 0d9c26afa49864652d58bef0898c6ecdf4b1105b..73f0cd975ce4a814610602854e5752489abcce58 100644
--- a/apps/showcases/ChargedParticles/poisson_solver/TestIterativeSolvers.cpp
+++ b/apps/showcases/ChargedParticles/poisson_solver/TestIterativeSolvers.cpp
@@ -6,32 +6,33 @@
 #include "core/Environment.h"
 #include "core/grid_generator/SCIterator.h"
 #include "core/logging/all.h"
+#include "core/math/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 {
+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) {
-
+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();
+      auto cellCenter = cellAABB.center();
 
       // snap cell position to actual domain position
       switch (direction) {
@@ -77,22 +78,25 @@ void applyDirichletFunction(const shared_ptr< StructuredBlockStorage > & blocks,
       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);
+            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( sinh ( math::root_two * math::pi * boundaryCoord_z ) );
+            funcVal = real_c(sin(math::pi * boundaryCoord_x)) * real_c(sin(math::pi * boundaryCoord_y)) *
+                      real_c(sinh(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) {
+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);
 
@@ -102,8 +106,10 @@ void resetSolution(const shared_ptr< StructuredBlockStorage > & blocks, BlockDat
    }
 }
 
-void resetRHS(const shared_ptr< StructuredBlockStorage > & blocks, BlockDataID & rhs) {
-   for (auto block = blocks->begin(); block != blocks->end(); ++block) {
+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
@@ -111,25 +117,32 @@ void resetRHS(const shared_ptr< StructuredBlockStorage > & blocks, BlockDataID &
    }
 }
 
-// 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 uint_t numIter, real_t resThres, uint_t resCheckFreq) {
-
+// 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 uint_t numIter, bool useAbsResThres, real_t resThres,
+           uint_t resCheckFreq)
+{
    const bool useDirichlet = testcase == TEST_DIRICHLET_1 || testcase == TEST_DIRICHLET_2;
 
    // set boundary handling depending on scenario
-   std::function< void () > boundaryHandling = {};
+   std::function< void() > boundaryHandling = {};
+   std::vector< BoundaryCondition > boundaryConditions;
 
-   if constexpr (useDirichlet) {
+   if constexpr (useDirichlet)
+   {
       // set dirichlet function per domain face
+      for (const auto& cond : stencil::D3Q6::dir)
+         boundaryConditions.emplace_back(cond, "Dirichlet", 0_r); // dummy value, unused
+
       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); \
-         }
+   [&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));
@@ -139,22 +152,38 @@ void solve(const shared_ptr< StructuredBlockForest > & blocks,
       dirichletFunction.setFunction(stencil::T, GET_BOUNDARY_LAMBDA(stencil::T));
 
       boundaryHandling = dirichletFunction;
-   } else {
-       boundaryHandling = pde::NeumannDomainBoundary< ScalarField_T >(*blocks, solution);
    }
+   else
+   {
+      // set neumann bcs per domain face
+      for (const auto& cond : stencil::D3Q6::dir)
+         boundaryConditions.emplace_back(cond, "Neumann", 0_r); // dummy value, unused
 
-   // solvers: Jacobi and SOR
+      boundaryHandling = NeumannDomainBoundary< ScalarField_T >(*blocks, solution, boundaryConditions);
+   }
 
-   auto poissonSolverJacobi = PoissonSolver< WALBERLA_JACOBI > (solution, solutionCpy, rhs, blocks, boundaryHandling, numIter, resThres, resCheckFreq);
-   auto poissonSolverDampedJac = PoissonSolver< DAMPED_JACOBI > (solution, solutionCpy, rhs, blocks, boundaryHandling, numIter, resThres, resCheckFreq);
-   auto poissonSolverSOR = PoissonSolver< WALBERLA_SOR > (solution, solutionCpy, rhs, blocks, boundaryHandling, numIter, resThres, resCheckFreq);
+   // solvers: (damped) Jacobi, SOR and CG
+
+   auto poissonSolverJacobi =
+      PoissonSolver< WALBERLA_JACOBI >(solution, solutionCpy, rhs, blocks, boundaryHandling, boundaryConditions,
+                                       numIter, useAbsResThres, resThres, resThres, resCheckFreq);
+   auto poissonSolverDampedJac =
+      PoissonSolver< DAMPED_JACOBI >(solution, solutionCpy, rhs, blocks, boundaryHandling, boundaryConditions, numIter,
+                                     useAbsResThres, resThres, resThres, resCheckFreq);
+   auto poissonSolverSOR =
+      PoissonSolver< WALBERLA_SOR >(solution, solutionCpy, rhs, blocks, boundaryHandling, boundaryConditions, numIter,
+                                    useAbsResThres, resThres, resThres, resCheckFreq);
+   auto poissonSolverCG =
+      PoissonSolver< WALBERLA_CG >(solution, solutionCpy, rhs, blocks, boundaryHandling, boundaryConditions, numIter,
+                                   useAbsResThres, resThres, resThres, resCheckFreq);
 
    // 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) {
+      for (auto block = blocks->begin(); block != blocks->end(); ++block)
+      {
          ScalarField_T* solutionField = block->getData< ScalarField_T >(solution);
 
          real_t blockResult = real_t(0);
@@ -163,88 +192,86 @@ void solve(const shared_ptr< StructuredBlockForest > & blocks,
                                     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);
+                                    auto scaleX = real_c(1) / domainAABB.size(0);
+                                    auto scaleY = real_c(1) / domainAABB.size(1);
+                                    auto 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;
+                                    auto posX = cellCenter[0] * scaleX;
+                                    auto posY = cellCenter[1] * scaleY;
+                                    auto 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);
+                                          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 ) );
+                                          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 ) );
+                                          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));
+                                    auto 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 );
+      mpi::allReduceInplace(error, mpi::MAX);
 
       return error;
    };
 
    // init rhs depending on scenario
 
-   for (auto block = blocks->begin(); block != blocks->end(); ++block) {
+   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();
+         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);
+         auto scaleX = real_c(1) / domainAABB.size(0);
+         auto scaleY = real_c(1) / domainAABB.size(1);
+         auto 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;
+         auto posX = cellCenter[0] * scaleX;
+         auto posY = cellCenter[1] * scaleY;
+         auto 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 ) );
+               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) );
+               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();
@@ -269,68 +296,90 @@ void solve(const shared_ptr< StructuredBlockForest > & blocks,
    poissonSolverSOR();
    auto errSOR = computeMaxError();
    WALBERLA_LOG_INFO_ON_ROOT("Error after SOR solver is: " << errSOR);
+
+   // solve with CG
+   WALBERLA_LOG_INFO_ON_ROOT("-- Solve using CG --");
+   resetSolution(blocks, solution, solutionCpy); // reset solutions and solve anew
+   poissonSolverCG();
+   auto errCG = computeMaxError();
+   WALBERLA_LOG_INFO_ON_ROOT("Error after CG solver is: " << errCG);
 }
 
 // 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
+template< bool useDirichlet >
+void solveChargedParticles(const shared_ptr< StructuredBlockForest >& blocks, const math::AABB& domainAABB,
+                           BlockDataID& solution, BlockDataID& solutionCpy, BlockDataID& rhs)
+{
+   // solvers: damped Jacobi, SOR and CG
 
-   auto numIter = uint_c(20000);
-   auto resThres = real_c(1e-5);
-   auto resCheckFreq = uint_c(1000);
+   auto numIter            = uint_c(20000);
+   auto useAbsResThreshold = true;
+   auto resThres           = real_c(1e-5);
+   auto resCheckFreq       = uint_c(1000);
 
    // set boundary handling depending on scenario
-   std::function< void () > boundaryHandling = {};
+   std::function< void() > boundaryHandling = {};
    std::vector< BoundaryCondition > boundaryConditions;
 
-   if constexpr (useDirichlet) {
+   if constexpr (useDirichlet)
+   {
       for (const auto& e : stencil::D3Q6::dir)
          boundaryConditions.emplace_back(e, "Dirichlet", 0_r);
 
       boundaryHandling = DirichletDomainBoundary< ScalarField_T >(*blocks, solution, boundaryConditions);
-   } else {
-      boundaryHandling = pde::NeumannDomainBoundary< ScalarField_T >(*blocks, solution);
    }
+   else {
+      for (const auto& e : stencil::D3Q6::dir)
+         boundaryConditions.emplace_back(e, "Neumann", 0_r);
 
-   auto poissonSolverJacobi = PoissonSolver< DAMPED_JACOBI > (solution, solutionCpy, rhs, blocks, boundaryHandling, numIter, resThres, resCheckFreq);
-   auto poissonSolverSOR = PoissonSolver< WALBERLA_SOR > (solution, solutionCpy, rhs, blocks, boundaryHandling, numIter, resThres, resCheckFreq);
+      boundaryHandling = NeumannDomainBoundary< ScalarField_T >(*blocks, solution, boundaryConditions);
+   }
+
+   auto poissonSolverJacobi =
+      PoissonSolver< DAMPED_JACOBI >(solution, solutionCpy, rhs, blocks, boundaryHandling, boundaryConditions, numIter,
+                                     useAbsResThreshold, resThres, resThres, resCheckFreq);
+   auto poissonSolverSOR =
+      PoissonSolver< WALBERLA_SOR >(solution, solutionCpy, rhs, blocks, boundaryHandling, boundaryConditions, numIter,
+                                    useAbsResThreshold, resThres, resThres, resCheckFreq);
+   auto poissonSolverCG =
+      PoissonSolver< WALBERLA_CG >(solution, solutionCpy, rhs, blocks, boundaryHandling, boundaryConditions, numIter,
+                                   useAbsResThreshold, resThres, resThres, resCheckFreq);
 
    // init rhs with two charged particles
 
-   for (auto block = blocks->begin(); block != blocks->end(); ++block) {
+   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();
+         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 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);
+         const real_t r1 = real_c(0.08) * domainAABB.size(0); const real_t s1 = real_c(1);
+
+         auto posX = cellCenter[0];
+         auto posY = cellCenter[1];
+         auto 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);
          }
@@ -343,6 +392,10 @@ void solveChargedParticles(const shared_ptr< StructuredBlockForest > & blocks,
    // solve with SOR
    resetSolution(blocks, solution, solutionCpy); // reset solutions and solve anew
    poissonSolverSOR();
+
+   // solve with CG
+   resetSolution(blocks, solution, solutionCpy); // reset solutions and solve anew
+   poissonSolverCG();
 }
 
 int main(int argc, char** argv)
@@ -352,7 +405,7 @@ int main(int argc, char** 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" );
+   logging::Logging::instance()->includeLoggingToFile("log");
 
    ///////////////////////////
    // BLOCK STRUCTURE SETUP //
@@ -361,45 +414,46 @@ int main(int argc, char** argv)
    auto cfgFile = env.config();
    if (!cfgFile) WALBERLA_ABORT("Usage: " << argv[0] << " < path-to-configuration-file > \n");
 
-   Config::BlockHandle setup                   = cfgFile->getBlock("Setup");
+   Config::BlockHandle const setup             = cfgFile->getBlock("Setup");
    const Vector3< uint_t > numBlocksPerDim     = setup.getParameter< Vector3< uint_t > >("numBlocks");
    const Vector3< uint_t > numCellsBlockPerDim = setup.getParameter< Vector3< uint_t > >("numCellsPerBlock");
    const uint_t maxIterations                  = setup.getParameter< uint_t >("maxIterations");
+   const bool useAbsResThreshold               = setup.getParameter< bool >("useAbsResThreshold");
    const real_t resThres                       = setup.getParameter< real_t >("resThreshold");
    const uint_t resCheckFreq                   = setup.getParameter< uint_t >("resCheckFreq");
 
    auto domainAABB = math::AABB(0, 0, 0, 1, 1, 1);
-   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,
-      numBlocksPerDim[0], numBlocksPerDim[1], numBlocksPerDim[2],
-      numCellsBlockPerDim[0], numCellsBlockPerDim[1], numCellsBlockPerDim[2],
-      true,
-      false, false, false,
-      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);
+   WALBERLA_LOG_INFO_ON_ROOT("Domain sizes are: x = " << domainAABB.size(0) << ", y = " << domainAABB.size(1)
+                                                      << ", z = " << domainAABB.size(2));
+
+   shared_ptr< StructuredBlockForest > const blocks = blockforest::createUniformBlockGrid(
+      domainAABB, numBlocksPerDim[0], numBlocksPerDim[1], numBlocksPerDim[2], numCellsBlockPerDim[0],
+      numCellsBlockPerDim[1], numCellsBlockPerDim[2], true, false, false, false, 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)");
 
    /* 1. analytical tests */
 
    // 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, maxIterations, resThres, resCheckFreq);
+   WALBERLA_LOG_INFO_ON_ROOT("- Solving analytical Neumann problem with (damped) Jacobi, SOR and CG... -")
+   solve< TEST_NEUMANN >(blocks, domainAABB, solution, solutionCpy, rhs, maxIterations, useAbsResThreshold, resThres,
+                         resCheckFreq);
 
    // ... then solve dirichlet problems
    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, maxIterations, resThres, resCheckFreq);
+   WALBERLA_LOG_INFO_ON_ROOT("- Solving analytical Dirichlet problem (1) with (damped) Jacobi, SOR and CG... -")
+   solve< TEST_DIRICHLET_1 >(blocks, domainAABB, solution, solutionCpy, rhs, maxIterations, useAbsResThreshold,
+                             resThres, resCheckFreq);
 
    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, maxIterations, resThres, resCheckFreq);
+   WALBERLA_LOG_INFO_ON_ROOT("- Solving analytical Dirichlet problem (2) with (damped) Jacobi, SOR and CG... -")
+   solve< TEST_DIRICHLET_2 >(blocks, domainAABB, solution, solutionCpy, rhs, maxIterations, useAbsResThreshold,
+                             resThres, resCheckFreq);
 
    /* 2. experimental charged particle tests */
 
@@ -408,13 +462,13 @@ int main(int argc, char** argv)
    resetSolution(blocks, solution, solutionCpy);
 
    // neumann
-   WALBERLA_LOG_INFO_ON_ROOT("- Run charged particles with Neumann boundaries... -")
-   solveChargedParticles < false > (blocks, domainAABB, solution, solutionCpy, rhs);
+   WALBERLA_LOG_INFO_ON_ROOT("- Run charged particles with Neumann boundaries with damped Jacobi, SOR and CG... -")
+   solveChargedParticles< false >(blocks, domainAABB, solution, solutionCpy, rhs);
 
    // dirichlet
-   WALBERLA_LOG_INFO_ON_ROOT("- Run charged particles with Dirichlet (val = 0) boundaries... -")
+   WALBERLA_LOG_INFO_ON_ROOT("- Run charged particles with Dirichlet (val = 0) boundaries with damped Jacobi, SOR and CG... -")
    resetSolution(blocks, solution, solutionCpy); // reset solutions and solve anew
-   solveChargedParticles < true > (blocks, domainAABB, solution, solutionCpy, rhs);
+   solveChargedParticles< true >(blocks, domainAABB, solution, solutionCpy, rhs);
 
    return EXIT_SUCCESS;
 }
@@ -422,4 +476,3 @@ int main(int argc, char** argv)
 } // namespace walberla
 
 int main(int argc, char** argv) { walberla::main(argc, argv); }
-