diff --git a/.editorconfig b/.editorconfig
new file mode 100644
index 0000000000000000000000000000000000000000..ea4b346c12c70e072c2d17dd430ec55c08385747
--- /dev/null
+++ b/.editorconfig
@@ -0,0 +1,10 @@
+# See https://editorconfig.org/
+root = true # top-most .editorconfig-file
+
+[*]
+tab_width = 3
+indent_style = space
+indent_size = 3
+charset = utf-8
+trim_trailing_whitespace = true
+insert_final_newline = false
\ No newline at end of file
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index cd9f1cd494e09dd2fbdddf55709de77075ab3219..28ddba36c231a862017112b951efa085530fecdd 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -747,7 +747,6 @@ gcc_7_mpionly:
 gcc_7_hybrid:
    <<: *build_definition
    image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:7
-   stage: pretest
    variables:
       <<: *build_hybrid_variables
       WALBERLA_BUILD_WITH_CUDA: "ON"
@@ -755,7 +754,6 @@ gcc_7_hybrid:
    except:
       variables:
          - $DISABLE_PER_COMMIT_BUILDS
-         - $DISABLE_PER_COMMIT_BUILDS
    tags:
       - cuda
       - docker
@@ -770,7 +768,6 @@ gcc_7_serial_dbg:
    except:
       variables:
          - $DISABLE_PER_COMMIT_BUILDS
-         - $DISABLE_PER_COMMIT_BUILDS
    tags:
       - cuda
       - docker
@@ -782,9 +779,9 @@ gcc_7_mpionly_dbg:
       <<: *build_mpionly_dbg_variables
       WALBERLA_BUILD_WITH_CUDA: "ON"
       WALBERLA_ENABLE_GUI: 0
-   except:
+   only:
       variables:
-         - $DISABLE_PER_COMMIT_BUILDS
+         - $ENABLE_NIGHTLY_BUILDS
    tags:
       - cuda
       - docker
@@ -796,9 +793,9 @@ gcc_7_hybrid_dbg:
       <<: *build_hybrid_dbg_variables
       WALBERLA_BUILD_WITH_CUDA: "ON"
       WALBERLA_ENABLE_GUI: 0
-   except:
+   only:
       variables:
-         - $DISABLE_PER_COMMIT_BUILDS
+         - $ENABLE_NIGHTLY_BUILDS
    tags:
       - cuda
       - docker
@@ -810,11 +807,103 @@ gcc_7_hybrid_dbg_sp:
       <<: *build_hybrid_dbg_sp_variables
       WALBERLA_BUILD_WITH_CUDA: "ON"
       WALBERLA_ENABLE_GUI: 0
+   only:
+      variables:
+         - $ENABLE_NIGHTLY_BUILDS
+   tags:
+      - cuda
+      - docker
+
+gcc_8_serial:
+   <<: *build_definition
+   image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:8
+   variables:
+      <<: *build_serial_variables
+      WALBERLA_BUILD_WITH_CUDA: "OFF"
+      WALBERLA_ENABLE_GUI: 0
+   only:
+      variables:
+         - $ENABLE_NIGHTLY_BUILDS
+   tags:
+      - docker
+
+gcc_8_mpionly:
+   <<: *build_definition
+   image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:8
+   variables:
+      <<: *build_mpionly_variables
+      WALBERLA_BUILD_WITH_CUDA: "OFF"
+      WALBERLA_ENABLE_GUI: 0
+   only:
+      variables:
+         - $ENABLE_NIGHTLY_BUILDS
+   tags:
+      - docker
+
+gcc_8_hybrid:
+   <<: *build_definition
+   image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:8
+   stage: pretest
+   variables:
+      <<: *build_hybrid_variables
+      WALBERLA_BUILD_WITH_CUDA: "OFF"
+      WALBERLA_ENABLE_GUI: 0
+   except:
+      variables:
+         - $DISABLE_PER_COMMIT_BUILDS
+   tags:
+      - docker
+
+gcc_8_serial_dbg:
+   <<: *build_definition
+   image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:8
+   variables:
+      <<: *build_serial_dbg_variables
+      WALBERLA_BUILD_WITH_CUDA: "OFF"
+      WALBERLA_ENABLE_GUI: 0
+   except:
+      variables:
+         - $DISABLE_PER_COMMIT_BUILDS
+   tags:
+      - docker
+
+gcc_8_mpionly_dbg:
+   <<: *build_definition
+   image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:8
+   variables:
+      <<: *build_mpionly_dbg_variables
+      WALBERLA_BUILD_WITH_CUDA: "OFF"
+      WALBERLA_ENABLE_GUI: 0
+   except:
+      variables:
+         - $DISABLE_PER_COMMIT_BUILDS
+   tags:
+      - docker
+
+gcc_8_hybrid_dbg:
+   <<: *build_definition
+   image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:8
+   variables:
+      <<: *build_hybrid_dbg_variables
+      WALBERLA_BUILD_WITH_CUDA: "OFF"
+      WALBERLA_ENABLE_GUI: 0
+   except:
+      variables:
+         - $DISABLE_PER_COMMIT_BUILDS
+   tags:
+      - docker
+
+gcc_8_hybrid_dbg_sp:
+   <<: *build_definition
+   image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:8
+   variables:
+      <<: *build_hybrid_dbg_sp_variables
+      WALBERLA_BUILD_WITH_CUDA: "OFF"
+      WALBERLA_ENABLE_GUI: 0
    except:
       variables:
          - $DISABLE_PER_COMMIT_BUILDS
    tags:
-      - cuda
       - docker
 
 clang_3.6_serial:
diff --git a/apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/ForcesOnSphereNearPlaneInShearFlow.cpp b/apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/ForcesOnSphereNearPlaneInShearFlow.cpp
index ccb7cc6318cc9838efae0063809658005c5564f8..de8dd9d31c658d7d4a670acc7c507fd975b56e99 100644
--- a/apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/ForcesOnSphereNearPlaneInShearFlow.cpp
+++ b/apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/ForcesOnSphereNearPlaneInShearFlow.cpp
@@ -361,7 +361,7 @@ private:
 };
 
 void initializeCouetteProfile( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDataID & pdfFieldID, const BlockDataID & boundaryHandlingID,
-                               const real_t & substrateHeight, const real_t & domainHeight, const real_t wallVelocity )
+                               const real_t & domainHeight, const real_t wallVelocity )
 {
    const real_t rho = real_c(1);
 
@@ -377,10 +377,8 @@ void initializeCouetteProfile( const shared_ptr< StructuredBlockStorage > & bloc
 
                                  Vector3< real_t > velocity( real_c(0) );
 
-                                 if( coord[2] >= substrateHeight )
-                                 {
-                                    velocity[0] =  wallVelocity * (coord[2] - substrateHeight) / ( domainHeight - substrateHeight);
-                                 }
+                                 velocity[0] =  wallVelocity * coord[2] / domainHeight;
+
                                  pdfField->setToEquilibrium( x, y, z, velocity, rho );
       )
    }
@@ -632,7 +630,7 @@ int main( int argc, char **argv )
 
 
    // initialize Couette velocity profile in whole domain
-   if( !zeroShearTest ) initializeCouetteProfile(blocks, pdfFieldID, boundaryHandlingID, diameter, domainHeight, wallVelocity);
+   if( !zeroShearTest ) initializeCouetteProfile(blocks, pdfFieldID, boundaryHandlingID, domainHeight, wallVelocity);
 
    ///////////////
    // TIME LOOP //
diff --git a/apps/benchmarks/UniformGrid/UniformGrid.cpp b/apps/benchmarks/UniformGrid/UniformGrid.cpp
index c270b38ca3fcd07869ba7a25842b1dd90f196c0c..0324c7190e3cd1692541a433d33992ab4c16f260 100644
--- a/apps/benchmarks/UniformGrid/UniformGrid.cpp
+++ b/apps/benchmarks/UniformGrid/UniformGrid.cpp
@@ -263,19 +263,28 @@ void createSetupBlockForest( blockforest::SetupBlockForest & sforest, const Conf
 
       WALBERLA_MPI_SECTION()
       {
-         MPIManager::instance()->createCartesianComm( numberOfXProcesses, numberOfYProcesses, numberOfZProcesses, false, false, false );
+         if ( MPIManager::instance()->isCartesianCommValid() )
+         {
+            MPIManager::instance()->createCartesianComm(numberOfXProcesses, numberOfYProcesses, numberOfZProcesses, false, false, false);
 
-         processIdMap = make_shared< std::vector<uint_t> >( numberOfXProcesses * numberOfYProcesses * numberOfZProcesses );
+            processIdMap = make_shared<std::vector<uint_t> >(numberOfXProcesses * numberOfYProcesses * numberOfZProcesses);
 
-         for( uint_t z = 0; z != numberOfZProcesses; ++z ) {
-            for( uint_t y = 0; y != numberOfYProcesses; ++y ) {
-               for( uint_t x = 0; x != numberOfXProcesses; ++x )
-               {
-                  (*processIdMap)[z * numberOfXProcesses * numberOfYProcesses + y * numberOfXProcesses + x] =
-                     uint_c( MPIManager::instance()->cartesianRank( x, y, z ) );
+            for (uint_t z = 0; z != numberOfZProcesses; ++z) {
+               for (uint_t y = 0; y != numberOfYProcesses; ++y) {
+                  for (uint_t x = 0; x != numberOfXProcesses; ++x)
+                  {
+                     (*processIdMap)[z * numberOfXProcesses * numberOfYProcesses + y * numberOfXProcesses + x] =
+                           uint_c(MPIManager::instance()->cartesianRank(x, y, z));
+                  }
                }
             }
          }
+         else {
+            WALBERLA_LOG_WARNING_ON_ROOT( "Your version of OpenMPI contains a bug. See waLBerla issue #73 for more "
+                                          "information. As a workaround, MPI_COMM_WORLD instead of a "
+                                          "Cartesian MPI communicator is used." );
+            MPIManager::instance()->useWorldComm();
+         }
       }
       else
          MPIManager::instance()->useWorldComm();
diff --git a/src/blockforest/Initialization.cpp b/src/blockforest/Initialization.cpp
index acb638accd0ea043a8db2f249eb98a5ea3084e22..c13896c87c7966d3adc9422828d1006f7be7aa6c 100644
--- a/src/blockforest/Initialization.cpp
+++ b/src/blockforest/Initialization.cpp
@@ -229,19 +229,27 @@ createBlockForest(      const AABB& domainAABB,
       //create cartesian communicator only if not yet a cartesian communicator (or other communicator was created)
       if ( ! mpiManager->rankValid() )
       {
-         mpiManager->createCartesianComm( numberOfXProcesses, numberOfYProcesses, numberOfZProcesses, xPeriodic, yPeriodic, zPeriodic );
+         if ( mpiManager->isCartesianCommValid() ) {
+            mpiManager->createCartesianComm( numberOfXProcesses, numberOfYProcesses, numberOfZProcesses, xPeriodic, yPeriodic, zPeriodic );
 
-         processIdMap = new std::vector< uint_t >( numberOfProcesses );
+            processIdMap = new std::vector< uint_t >( numberOfProcesses );
 
-         for( uint_t z = 0; z != numberOfZProcesses; ++z ) {
-            for( uint_t y = 0; y != numberOfYProcesses; ++y ) {
-               for( uint_t x = 0; x != numberOfXProcesses; ++x ) {
+            for( uint_t z = 0; z != numberOfZProcesses; ++z ) {
+               for( uint_t y = 0; y != numberOfYProcesses; ++y ) {
+                  for( uint_t x = 0; x != numberOfXProcesses; ++x ) {
 
-                  (*processIdMap)[ z * numberOfXProcesses * numberOfYProcesses + y * numberOfXProcesses + x ] =
-                        uint_c( MPIManager::instance()->cartesianRank(x,y,z) );
+                     (*processIdMap)[ z * numberOfXProcesses * numberOfYProcesses + y * numberOfXProcesses + x ] =
+                           uint_c( MPIManager::instance()->cartesianRank(x,y,z) );
+                  }
                }
             }
          }
+         else {
+            WALBERLA_LOG_WARNING_ON_ROOT( "Your version of OpenMPI contains a bug. See waLBerla issue #73 for more "
+                                          "information. As a workaround, MPI_COMM_WORLD instead of a "
+                                          "Cartesian MPI communicator is used." );
+            mpiManager->useWorldComm();
+         }
       }
    }
 
diff --git a/src/core/OpenMP.h b/src/core/OpenMP.h
index a7ea139c12e08ca023b186c3c1f948e7b5f6f121..825046dc1570c7106aeb30c0b4f5958617207afa 100644
--- a/src/core/OpenMP.h
+++ b/src/core/OpenMP.h
@@ -14,17 +14,155 @@
 //  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
 //
 //! \file OpenMP.h
-//! \ingroup core
 //! \brief Guarded OpenMP include
 //
 //======================================================================================================================
 
 #pragma once
 
+#include "core/Abort.h"
 
+// MPI SECTION //
+
+#ifdef WALBERLA_BUILD_WITH_OPENMP
+
+#define     WALBERLA_OPENMP_SECTION() if(true)
+#define WALBERLA_NON_OPENMP_SECTION() if(false)
+
+#else
+
+#define     WALBERLA_OPENMP_SECTION() if(false)
+#define WALBERLA_NON_OPENMP_SECTION() if(true)
+
+#endif
 
 #ifdef _OPENMP
+
 #include <omp.h>
+
+#else
+
+namespace walberla {
+
+#define WALBERLA_OPENMP_FUNCTION_ERROR WALBERLA_ABORT( "Invalid OpenMP function call! In case of compiling without OpenMP, OpenMP functions are not available and shouldn't be called!" );
+
+/* schedule kind constants */
+typedef enum omp_sched_t {
+omp_sched_static  = 1,
+omp_sched_dynamic = 2,
+omp_sched_guided  = 3,
+omp_sched_auto    = 4
+} omp_sched_t;
+
+/* set API functions */
+inline void    omp_set_num_threads (int) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline void    omp_set_dynamic     (int) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline void    omp_set_nested      (int) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline void    omp_set_max_active_levels (int) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline void    omp_set_schedule          (omp_sched_t, int) { WALBERLA_OPENMP_FUNCTION_ERROR }
+
+/* query API functions */
+inline int     omp_get_num_threads  (void) { return 1; }
+inline int     omp_get_dynamic      (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int     omp_get_nested       (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int     omp_get_max_threads  (void) { return 1; }
+inline int     omp_get_thread_num   (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int     omp_get_num_procs    (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int     omp_in_parallel      (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int     omp_in_final         (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int     omp_get_active_level        (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int     omp_get_level               (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int     omp_get_ancestor_thread_num (int) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int     omp_get_team_size           (int) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int     omp_get_thread_limit        (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int     omp_get_max_active_levels   (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline void    omp_get_schedule            (omp_sched_t *, int *) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int     omp_get_max_task_priority   (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+
+/* lock API functions */
+typedef struct omp_lock_t {
+    void * _lk;
+} omp_lock_t;
+
+inline void    omp_init_lock    (omp_lock_t *) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline void    omp_set_lock     (omp_lock_t *) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline void    omp_unset_lock   (omp_lock_t *) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline void    omp_destroy_lock (omp_lock_t *) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int     omp_test_lock    (omp_lock_t *) { WALBERLA_OPENMP_FUNCTION_ERROR }
+
+/* nested lock API functions */
+typedef struct omp_nest_lock_t {
+    void * _lk;
+} omp_nest_lock_t;
+
+inline void    omp_init_nest_lock    (omp_nest_lock_t *) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline void    omp_set_nest_lock     (omp_nest_lock_t *) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline void    omp_unset_nest_lock   (omp_nest_lock_t *) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline void    omp_destroy_nest_lock (omp_nest_lock_t *) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int     omp_test_nest_lock    (omp_nest_lock_t *) { WALBERLA_OPENMP_FUNCTION_ERROR }
+
+/* lock hint type for dynamic user lock */
+typedef enum omp_lock_hint_t {
+    omp_lock_hint_none           = 0,
+    omp_lock_hint_uncontended    = 1,
+    omp_lock_hint_contended      = (1<<1 ),
+    omp_lock_hint_nonspeculative = (1<<2 ),
+    omp_lock_hint_speculative    = (1<<3 ),
+    kmp_lock_hint_hle            = (1<<16),
+    kmp_lock_hint_rtm            = (1<<17),
+    kmp_lock_hint_adaptive       = (1<<18)
+} omp_lock_hint_t;
+
+/* hinted lock initializers */
+inline void omp_init_lock_with_hint(omp_lock_t *, omp_lock_hint_t) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline void omp_init_nest_lock_with_hint(omp_nest_lock_t *, omp_lock_hint_t) { WALBERLA_OPENMP_FUNCTION_ERROR }
+
+/* time API functions */
+inline double  omp_get_wtime (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline double  omp_get_wtick (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+
+/* OpenMP 4.0 */
+inline int   omp_get_default_device (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline void  omp_set_default_device (int) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int   omp_is_initial_device (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int   omp_get_num_devices (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int   omp_get_num_teams (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int   omp_get_team_num (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int   omp_get_cancellation (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+
+#   include <stdlib.h>
+/* OpenMP 4.5 */
+inline int    omp_get_initial_device (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline void*  omp_target_alloc(size_t, int) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline void   omp_target_free(void *, int) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int    omp_target_is_present(void *, int) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int    omp_target_memcpy(void *, void *, size_t, size_t, size_t, int, int) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int    omp_target_memcpy_rect(void *, void *, size_t, int, const size_t *,
+                                        const size_t *, const size_t *, const size_t *, const size_t *, int, int) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int    omp_target_associate_ptr(void *, void *, size_t, size_t, int) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int    omp_target_disassociate_ptr(void *, int) { WALBERLA_OPENMP_FUNCTION_ERROR }
+
+/* OpenMP 4.0 affinity API */
+typedef enum omp_proc_bind_t {
+    omp_proc_bind_false = 0,
+    omp_proc_bind_true = 1,
+    omp_proc_bind_master = 2,
+    omp_proc_bind_close = 3,
+    omp_proc_bind_spread = 4
+} omp_proc_bind_t;
+
+inline omp_proc_bind_t omp_get_proc_bind (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+
+/* OpenMP 4.5 affinity API */
+inline int  omp_get_num_places (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int  omp_get_place_num_procs (int) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline void omp_get_place_proc_ids (int, int *) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int  omp_get_place_num (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline int  omp_get_partition_num_places (void) { WALBERLA_OPENMP_FUNCTION_ERROR }
+inline void omp_get_partition_place_nums (int *) { WALBERLA_OPENMP_FUNCTION_ERROR }
+
+} //namespace walberla
+
 #endif
 
 
diff --git a/src/core/debug/CheckFunctions.h b/src/core/debug/CheckFunctions.h
index 4e64bbf358954947fd4fdcd01cf9808d93ed77c4..61808427be29dcb8d358e19b25a4c393c70ccd59 100644
--- a/src/core/debug/CheckFunctions.h
+++ b/src/core/debug/CheckFunctions.h
@@ -96,21 +96,21 @@
 
 
 
-#define WALBERLA_CHECK_2(X,MSG)                             { if( !walberla::debug::check_functions_detail::check              ( (X)      ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check              (           #X,     __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_CHECK_NULLPTR_2(X,MSG)                     { if( !walberla::debug::check_functions_detail::check_nullptr      ( (X)      ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_nullptr      ( (X),      #X,     __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_CHECK_NOT_NULLPTR_2(X,MSG)                 { if( !walberla::debug::check_functions_detail::check_not_nullptr  ( (X)      ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_not_nullptr  (           #X,     __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_CHECK_EQUAL_3(X,Y,MSG)                     { if( !walberla::debug::check_functions_detail::check_equal        ( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_equal        ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_CHECK_UNEQUAL_3(X,Y,MSG)                   { if( !walberla::debug::check_functions_detail::check_unequal      ( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_unequal      ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_CHECK_FLOAT_EQUAL_3(X,Y,MSG)               { if( !walberla::debug::check_functions_detail::check_float_equal  ( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_float_equal  ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_CHECK_FLOAT_UNEQUAL_3(X,Y,MSG)             { if( !walberla::debug::check_functions_detail::check_float_unequal( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_float_unequal( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_CHECK_FLOAT_EQUAL_EPSILON_4(X,Y,EPS,MSG)   { if( !walberla::debug::check_functions_detail::check_float_equal_eps  ( (X), (Y), (EPS) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_float_equal_eps  ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ), (EPS) ); } }
-#define WALBERLA_CHECK_FLOAT_UNEQUAL_EPSILON_4(X,Y,EPS,MSG) { if( !walberla::debug::check_functions_detail::check_float_unequal_eps( (X), (Y), (EPS) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_float_unequal_eps( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ), (EPS) ); } }
-#define WALBERLA_CHECK_IDENTICAL_3(X,Y,MSG)                 { if( !walberla::debug::check_functions_detail::check_identical    ( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_identical    ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_CHECK_NOT_IDENTICAL_3(X,Y,MSG)             { if( !walberla::debug::check_functions_detail::check_not_identical( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_not_identical( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_CHECK_LESS_3(X,Y,MSG)                      { if( !walberla::debug::check_functions_detail::check_less         ( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_less         ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_CHECK_GREATER_3(X,Y,MSG)                   { if( !walberla::debug::check_functions_detail::check_greater      ( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_greater      ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_CHECK_LESS_EQUAL_3(X,Y,MSG)                { if( !walberla::debug::check_functions_detail::check_less_equal   ( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_less_equal   ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_CHECK_GREATER_EQUAL_3(X,Y,MSG)             { if( !walberla::debug::check_functions_detail::check_greater_equal( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_greater_equal( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
+#define WALBERLA_CHECK_2(X,MSG)                             { if( !walberla::debug::check_functions_detail::check              ( (X)      ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check              (           #X,     __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_CHECK_NULLPTR_2(X,MSG)                     { if( !walberla::debug::check_functions_detail::check_nullptr      ( (X)      ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_nullptr      ( (X),      #X,     __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_CHECK_NOT_NULLPTR_2(X,MSG)                 { if( !walberla::debug::check_functions_detail::check_not_nullptr  ( (X)      ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_not_nullptr  (           #X,     __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_CHECK_EQUAL_3(X,Y,MSG)                     { if( !walberla::debug::check_functions_detail::check_equal        ( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_equal        ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_CHECK_UNEQUAL_3(X,Y,MSG)                   { if( !walberla::debug::check_functions_detail::check_unequal      ( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_unequal      ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_CHECK_FLOAT_EQUAL_3(X,Y,MSG)               { if( !walberla::debug::check_functions_detail::check_float_equal  ( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_float_equal  ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_CHECK_FLOAT_UNEQUAL_3(X,Y,MSG)             { if( !walberla::debug::check_functions_detail::check_float_unequal( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_float_unequal( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_CHECK_FLOAT_EQUAL_EPSILON_4(X,Y,EPS,MSG)   { if( !walberla::debug::check_functions_detail::check_float_equal_eps  ( (X), (Y), (EPS) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_float_equal_eps  ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ), (EPS) ); } }
+#define WALBERLA_CHECK_FLOAT_UNEQUAL_EPSILON_4(X,Y,EPS,MSG) { if( !walberla::debug::check_functions_detail::check_float_unequal_eps( (X), (Y), (EPS) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_float_unequal_eps( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ), (EPS) ); } }
+#define WALBERLA_CHECK_IDENTICAL_3(X,Y,MSG)                 { if( !walberla::debug::check_functions_detail::check_identical    ( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_identical    ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_CHECK_NOT_IDENTICAL_3(X,Y,MSG)             { if( !walberla::debug::check_functions_detail::check_not_identical( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_not_identical( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_CHECK_LESS_3(X,Y,MSG)                      { if( !walberla::debug::check_functions_detail::check_less         ( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_less         ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_CHECK_GREATER_3(X,Y,MSG)                   { if( !walberla::debug::check_functions_detail::check_greater      ( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_greater      ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_CHECK_LESS_EQUAL_3(X,Y,MSG)                { if( !walberla::debug::check_functions_detail::check_less_equal   ( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_less_equal   ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_CHECK_GREATER_EQUAL_3(X,Y,MSG)             { if( !walberla::debug::check_functions_detail::check_greater_equal( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_greater_equal( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
 
 
 
diff --git a/src/core/debug/Debug.h b/src/core/debug/Debug.h
index b29d2824e157033614e5cbf09ae4a6a076c70edd..ef717cfc0c95cd4dbcc58b1a63891f65a05a1667 100644
--- a/src/core/debug/Debug.h
+++ b/src/core/debug/Debug.h
@@ -81,24 +81,24 @@
 #define WALBERLA_ASSERT_LESS_EQUAL_2(X,Y)    { if( !walberla::debug::check_functions_detail::check_less_equal   ( (X), (Y) ) ) { walberla::debug::check_functions_detail::check_less_equal   ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler() ); } }
 #define WALBERLA_ASSERT_GREATER_EQUAL_2(X,Y) { if( !walberla::debug::check_functions_detail::check_greater_equal( (X), (Y) ) ) { walberla::debug::check_functions_detail::check_greater_equal( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler() ); } }
 
-#define WALBERLA_ASSERT_2(X,MSG)                         { if( !walberla::debug::check_functions_detail::check                  ( (X)      ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check              (           #X,     __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_ASSERT_NULLPTR_2(X,MSG)                 { if( !walberla::debug::check_functions_detail::check_nullptr          ( (X)      ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_nullptr      ( (X),      #X,     __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_ASSERT_NOT_NULLPTR_2(X,MSG)             { if( !walberla::debug::check_functions_detail::check_not_nullptr      ( (X)      ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_not_nullptr  (           #X,     __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_ASSERT_EQUAL_3(X,Y,MSG)                 { if( !walberla::debug::check_functions_detail::check_equal            ( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_equal        ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_ASSERT_UNEQUAL_3(X,Y,MSG)               { if( !walberla::debug::check_functions_detail::check_unequal          ( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_unequal      ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_ASSERT_FLOAT_EQUAL_3(X,Y,MSG)           { if( !walberla::debug::check_functions_detail::check_float_equal      ( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_float_equal  ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_ASSERT_FLOAT_UNEQUAL_3(X,Y,MSG)         { if( !walberla::debug::check_functions_detail::check_float_unequal    ( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_float_unequal( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
+#define WALBERLA_ASSERT_2(X,MSG)                         { if( !walberla::debug::check_functions_detail::check                  ( (X)      ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check              (           #X,     __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_ASSERT_NULLPTR_2(X,MSG)                 { if( !walberla::debug::check_functions_detail::check_nullptr          ( (X)      ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_nullptr      ( (X),      #X,     __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_ASSERT_NOT_NULLPTR_2(X,MSG)             { if( !walberla::debug::check_functions_detail::check_not_nullptr      ( (X)      ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_not_nullptr  (           #X,     __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_ASSERT_EQUAL_3(X,Y,MSG)                 { if( !walberla::debug::check_functions_detail::check_equal            ( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_equal        ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_ASSERT_UNEQUAL_3(X,Y,MSG)               { if( !walberla::debug::check_functions_detail::check_unequal          ( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_unequal      ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_ASSERT_FLOAT_EQUAL_3(X,Y,MSG)           { if( !walberla::debug::check_functions_detail::check_float_equal      ( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_float_equal  ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_ASSERT_FLOAT_UNEQUAL_3(X,Y,MSG)         { if( !walberla::debug::check_functions_detail::check_float_unequal    ( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_float_unequal( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
 #define WALBERLA_ASSERT_FLOAT_EQUAL_EPSILON_3(X,Y,EPS)   { if( !walberla::debug::check_functions_detail::check_float_equal_eps  ( (X), (Y), (EPS) ) ) { walberla::debug::check_functions_detail::check_float_equal_eps  ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ), (EPS) ); } }
 #define WALBERLA_ASSERT_FLOAT_UNEQUAL_EPSILON_3(X,Y,EPS) { if( !walberla::debug::check_functions_detail::check_float_unequal_eps( (X), (Y), (EPS) ) ) { walberla::debug::check_functions_detail::check_float_unequal_eps( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ), (EPS) ); } }
-#define WALBERLA_ASSERT_IDENTICAL_3(X,Y,MSG)             { if( !walberla::debug::check_functions_detail::check_identical        ( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_identical    ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_ASSERT_NOT_IDENTICAL_3(X,Y,MSG)         { if( !walberla::debug::check_functions_detail::check_not_identical    ( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_not_identical( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_ASSERT_LESS_3(X,Y,MSG)                  { if( !walberla::debug::check_functions_detail::check_less             ( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_less         ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_ASSERT_GREATER_3(X,Y,MSG)               { if( !walberla::debug::check_functions_detail::check_greater          ( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_greater      ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_ASSERT_LESS_EQUAL_3(X,Y,MSG)            { if( !walberla::debug::check_functions_detail::check_less_equal        ( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_less_equal   ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-#define WALBERLA_ASSERT_GREATER_EQUAL_3(X,Y,MSG)         { if( !walberla::debug::check_functions_detail::check_greater_equal    ( (X), (Y) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_greater_equal( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ) ); } }
-
-#define WALBERLA_ASSERT_FLOAT_EQUAL_EPSILON_4(X,Y,EPS,MSG)   { if( !walberla::debug::check_functions_detail::check_float_equal_eps  ( (X), (Y), (EPS) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_float_equal_eps  ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ), (EPS) ); } }
-#define WALBERLA_ASSERT_FLOAT_UNEQUAL_EPSILON_4(X,Y,EPS,MSG) { if( !walberla::debug::check_functions_detail::check_float_unequal_eps( (X), (Y), (EPS) ) ) { std::stringstream ss; ss << MSG; walberla::debug::check_functions_detail::check_float_unequal_eps( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( ss.str() ), (EPS) ); } }
+#define WALBERLA_ASSERT_IDENTICAL_3(X,Y,MSG)             { if( !walberla::debug::check_functions_detail::check_identical        ( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_identical    ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_ASSERT_NOT_IDENTICAL_3(X,Y,MSG)         { if( !walberla::debug::check_functions_detail::check_not_identical    ( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_not_identical( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_ASSERT_LESS_3(X,Y,MSG)                  { if( !walberla::debug::check_functions_detail::check_less             ( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_less         ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_ASSERT_GREATER_3(X,Y,MSG)               { if( !walberla::debug::check_functions_detail::check_greater          ( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_greater      ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_ASSERT_LESS_EQUAL_3(X,Y,MSG)            { if( !walberla::debug::check_functions_detail::check_less_equal        ( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_less_equal   ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+#define WALBERLA_ASSERT_GREATER_EQUAL_3(X,Y,MSG)         { if( !walberla::debug::check_functions_detail::check_greater_equal    ( (X), (Y) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_greater_equal( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ) ); } }
+
+#define WALBERLA_ASSERT_FLOAT_EQUAL_EPSILON_4(X,Y,EPS,MSG)   { if( !walberla::debug::check_functions_detail::check_float_equal_eps  ( (X), (Y), (EPS) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_float_equal_eps  ( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ), (EPS) ); } }
+#define WALBERLA_ASSERT_FLOAT_UNEQUAL_EPSILON_4(X,Y,EPS,MSG) { if( !walberla::debug::check_functions_detail::check_float_unequal_eps( (X), (Y), (EPS) ) ) { std::stringstream _ss_; _ss_ << MSG; walberla::debug::check_functions_detail::check_float_unequal_eps( (X), (Y), #X, #Y, __FILE__, __LINE__, walberla::debug::check_functions_detail::ExitHandler( _ss_.str() ), (EPS) ); } }
 
 #define WALBERLA_ASSERT_3(...) THIS_IS_SUPPOSED_TO_FAIL___YOU_MADE_AN_ERROR_WHEN_USING_AN_ASSERT_MACRO
 #define WALBERLA_ASSERT_4(...) THIS_IS_SUPPOSED_TO_FAIL___YOU_MADE_AN_ERROR_WHEN_USING_AN_ASSERT_MACRO
diff --git a/src/core/mpi/BufferSystem.cpp b/src/core/mpi/BufferSystem.cpp
index 551fce3e60dae2355ff448ab71e380fc9ca77437..a6e64b887be17fe5aaf3a91c160c61a2363b8467 100644
--- a/src/core/mpi/BufferSystem.cpp
+++ b/src/core/mpi/BufferSystem.cpp
@@ -53,6 +53,10 @@ void BufferSystem::iterator::operator++()
    currentRecvBuffer_ = bufferSystem_.waitForNext( currentSenderRank_ );
    if ( ! currentRecvBuffer_ ) {
       WALBERLA_ASSERT_EQUAL( currentSenderRank_, -1 );
+   } else
+   {
+      bufferSystem_.bytesReceived_ += currentRecvBuffer_->size() * sizeof(RecvBuffer::ElementType);
+      bufferSystem_.numberOfReceives_ += 1;
    }
 }
 
@@ -303,7 +307,11 @@ void BufferSystem::sendAll()
       if ( ! iter->second.alreadySent )
       {
          if ( iter->second.buffer.size() > 0 )
+         {
+            bytesSent_     += iter->second.buffer.size() * sizeof(SendBuffer::ElementType);
+            numberOfSends_ += 1;
             currentComm_->send( iter->first, iter->second.buffer );
+         }
 
          iter->second.alreadySent = true;
       }
@@ -331,7 +339,11 @@ void BufferSystem::send( MPIRank rank )
    WALBERLA_ASSERT( ! iter->second.alreadySent ); // this buffer has already been sent
 
    if ( iter->second.buffer.size() > 0 )
+   {
+      bytesSent_     += iter->second.buffer.size() * sizeof(SendBuffer::ElementType);
+      numberOfSends_ += 1;
       currentComm_->send( rank, iter->second.buffer );
+   }
 
    iter->second.alreadySent = true;
 }
@@ -362,6 +374,12 @@ void BufferSystem::startCommunication()
 
    currentComm_->scheduleReceives( recvInfos_ );
    communicationRunning_ = true;
+
+   bytesSent_        = 0;
+   bytesReceived_    = 0;
+
+   numberOfSends_    = 0;
+   numberOfReceives_ = 0;
 }
 
 
diff --git a/src/core/mpi/BufferSystem.h b/src/core/mpi/BufferSystem.h
index a79322d1ab326518156a7114f70cd024d2a9b041..926a7181ffca5cab7ff418cef67057aef20c2066 100644
--- a/src/core/mpi/BufferSystem.h
+++ b/src/core/mpi/BufferSystem.h
@@ -193,6 +193,16 @@ public:
    //@}
    //*******************************************************************************************************************
 
+   ///Bytes sent during the current or last communication
+   int64_t getBytesSent() const { return bytesSent_; }
+   ///Bytes received during the current or last communication
+   int64_t getBytesReceived() const { return bytesReceived_; }
+
+   ///Communication partners during current or last send operation
+   int64_t getNumberOfSends() const { return numberOfSends_; }
+   ///Communication partners during current or last receive operation
+   int64_t getNumberOfReceives() const { return numberOfReceives_; }
+
 
    //* Rank Ranges     *************************************************************************************************
    /*! \name Rank Ranges  */
@@ -240,6 +250,12 @@ protected:
    //stores tags of running communications in debug mode to ensure that
    //each concurrently running communication uses different tags
    static std::set<int> activeTags_;
+
+   int64_t bytesSent_        = 0; ///< number of bytes sent during last communication
+   int64_t bytesReceived_    = 0; ///< number of bytes received during last communication
+
+   int64_t numberOfSends_    = 0; ///< number of communication partners during last send
+   int64_t numberOfReceives_ = 0; ///< number of communication partners during last receive
 };
 
 
diff --git a/src/core/mpi/MPIManager.cpp b/src/core/mpi/MPIManager.cpp
index a48f7304d350523e0be6e8d88fc54f555bfc0928..72861cbd109a91321099f9e95771c88376a8a28c 100644
--- a/src/core/mpi/MPIManager.cpp
+++ b/src/core/mpi/MPIManager.cpp
@@ -22,12 +22,14 @@
 //======================================================================================================================
 
 #include "MPIManager.h"
+#include "core/logging/Logging.h"
 
 #include <boost/exception_ptr.hpp>
 #include <exception>
 #include <iostream>
 #include <stdexcept>
 #include <vector>
+#include <string>
 
 
 namespace walberla{
@@ -156,12 +158,18 @@ void MPIManager::createCartesianComm( int dims[3], int periodicity[3] )
    WALBERLA_ASSERT_GREATER( dims[1], 0 );
    WALBERLA_ASSERT_GREATER( dims[2], 0 );
 
-   MPI_Cart_create( MPI_COMM_WORLD, 3, dims, periodicity, true, &comm_ );
-
-   WALBERLA_ASSERT_UNEQUAL(comm_, MPI_COMM_NULL);
+   if ( ! isCartesianCommValid() ) {
+      WALBERLA_LOG_WARNING_ON_ROOT( "Your version of OpenMPI contains a bug which might lead to a segmentation fault "
+                                    "when generating vtk output. Since the bug only occurs with a 3D Cartesian MPI "
+                                    "communicator, try to use MPI_COMM_WORLD instead. See waLBerla issue #73 for "
+                                    "additional information." );
+   }
 
+   MPI_Cart_create( MPI_COMM_WORLD, 3, dims, periodicity, true, &comm_ );
    MPI_Comm_rank( comm_, &rank_ );
    cartesianSetup_ = true;
+
+   WALBERLA_ASSERT_UNEQUAL(comm_, MPI_COMM_NULL);
 }
 
 
@@ -225,6 +233,25 @@ int MPIManager::cartesianRank( const uint_t x, const uint_t y, const uint_t z )
    return cartesianRank( coords );
 }
 
+bool MPIManager::isCartesianCommValid() const
+{
+   // Avoid using the Cartesian MPI-communicator with certain (buggy) versions of OpenMPI (see waLBerla issue #73)
+   #if defined(OMPI_MAJOR_VERSION) && defined(OMPI_MINOR_VERSION) && defined(OMPI_RELEASE_VERSION)
+      std::string ompi_ver = std::to_string(OMPI_MAJOR_VERSION) + "." + std::to_string(OMPI_MINOR_VERSION) + "." +
+            std::to_string(OMPI_RELEASE_VERSION);
+
+      if (ompi_ver == "2.0.0" || ompi_ver == "2.0.1" || ompi_ver == "2.0.2" || ompi_ver == "2.0.3" ||
+          ompi_ver == "2.1.0" || ompi_ver == "2.1.1") {
+         return false;
+      }
+      else {
+         return true;
+      }
+   #else
+      return true;
+   #endif
+}
+
 std::string MPIManager::getMPIErrorString(int errorCode)
 {
    WALBERLA_NON_MPI_SECTION()
diff --git a/src/core/mpi/MPIManager.h b/src/core/mpi/MPIManager.h
index fe041c2675f72144bb696b225926b689c306bb35..fdbc3e06eab704fa63f10156119a4a5e5eeeeafe 100644
--- a/src/core/mpi/MPIManager.h
+++ b/src/core/mpi/MPIManager.h
@@ -122,6 +122,9 @@ public:
    bool hasCartesianSetup() const { return cartesianSetup_;  }
    /// Rank is valid after calling createCartesianComm() or useWorldComm()
    bool rankValid()         const { return rank_ >= 0;       }
+
+   /// Using a Cartesian MPI communicator is not valid for certain versions of OpenMPI (see waLBerla issue #73)
+   bool isCartesianCommValid() const;
    //@}
    //*******************************************************************************************************************
 
diff --git a/src/field/interpolators/KernelFieldInterpolator.h b/src/field/interpolators/KernelFieldInterpolator.h
index eff9e888a6fa03df0d6b06cb015beb158ee5dd73..50dfde0fc9a370fb1e4b8af95b37ccfe3bae9d4d 100644
--- a/src/field/interpolators/KernelFieldInterpolator.h
+++ b/src/field/interpolators/KernelFieldInterpolator.h
@@ -40,7 +40,7 @@ namespace kernelweights
 
 // corresponds to the smoothed dirac delta function from Roma et al - An Adaptive Version of the Immersed Boundary Method
 // f(r) != 0 for abs(r) < 1.5 -> requires three neighborhood cells
-real_t smoothedDeltaFunction( const real_t & r )
+inline real_t smoothedDeltaFunction( const real_t & r )
 {
    real_t rAbs = std::fabs(r);
    if( rAbs <= real_t(0.5) )
@@ -58,17 +58,17 @@ real_t smoothedDeltaFunction( const real_t & r )
 
 // X: Lagrangian position, x: Eulerian position (usually cell center), global coordinates
 // dx, dy, dz: mesh spacing
-real_t kernelWeightFunction( const real_t & X, const real_t & Y, const real_t & Z,
-                             const real_t & x, const real_t & y, const real_t & z,
-                             const real_t & dx = real_t(1), const real_t & dy = real_t(1), const real_t & dz = real_t(1) )
+inline real_t kernelWeightFunction( const real_t & X, const real_t & Y, const real_t & Z,
+                                    const real_t & x, const real_t & y, const real_t & z,
+                                    const real_t & dx = real_t(1), const real_t & dy = real_t(1), const real_t & dz = real_t(1) )
 {
    return smoothedDeltaFunction( ( X - x ) / dx ) * smoothedDeltaFunction( ( Y - y ) / dy ) * smoothedDeltaFunction( ( Z - z ) / dz );
 }
 
 // X: Lagrangian position, x: Eulerian position (usually cell center), global coordinates
 // dx, dy, dz: mesh spacing
-real_t kernelWeightFunction( const Vector3<real_t> & X, const Vector3<real_t> & x,
-                             const real_t & dx = real_t(1), const real_t & dy = real_t(1), const real_t & dz = real_t(1) )
+inline real_t kernelWeightFunction( const Vector3<real_t> & X, const Vector3<real_t> & x,
+                                    const real_t & dx = real_t(1), const real_t & dy = real_t(1), const real_t & dz = real_t(1) )
 {
    return smoothedDeltaFunction( ( X[0] - x[0] ) / dx ) * smoothedDeltaFunction( ( X[1] - x[1] ) / dy ) * smoothedDeltaFunction( ( X[2] - x[2] ) / dz );
 }
diff --git a/src/pde/CMakeLists.txt b/src/pde/CMakeLists.txt
index 386a98c37d55db966c79de92b58eda5331724fdb..0b013dda8b65a581f5f5b99e300f84b1f2a0924e 100644
--- a/src/pde/CMakeLists.txt
+++ b/src/pde/CMakeLists.txt
@@ -4,7 +4,8 @@
 #
 ###################################################################################################
 
-waLBerla_add_module( DEPENDS blockforest 
+waLBerla_add_module( DEPENDS blockforest
+                             boundary 
                              core
                              domain_decomposition
                              field
diff --git a/src/pde/boundary/Dirichlet.h b/src/pde/boundary/Dirichlet.h
new file mode 100644
index 0000000000000000000000000000000000000000..f5543fdbe95ca04f90bf559604984c1f268a1d3f
--- /dev/null
+++ b/src/pde/boundary/Dirichlet.h
@@ -0,0 +1,339 @@
+//======================================================================================================================
+//
+//  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 Dirichlet.h
+//! \ingroup pde
+//! \author Dominik Bartuschat <dominik.bartuschat@fau.de>
+//! \author Michael Kuron <mkuron@icp.uni-stuttgart.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "field/Field.h"
+
+#include "boundary/Boundary.h"
+
+#include "core/DataTypes.h"
+#include "core/cell/CellInterval.h"
+#include "core/config/Config.h"
+#include "core/debug/Debug.h"
+#include "core/logging/Logging.h"
+
+#include "field/FlagUID.h"
+
+#include "stencil/Directions.h"
+
+#include <vector>
+#include <limits>
+
+
+namespace walberla {
+namespace pde {
+
+
+
+//**********************************************************************************************************************
+/*!
+*   \brief Dirichlet boundary handling for PDEs
+*
+*   This boundary condition imposes a Dirichlet condition with arbitrary values on a PDE.
+*   It does so by modifying the right-hand side and the stencil field.
+*   Anything that has internal copies of the stencil field (e.g. the multigrid V-cycle's coarsened stencils) is
+*   responsible for updating its copies when boundary conditions are changed.
+*
+*   \tparam Stencil_T The stencil used for the discrete operator
+*   \tparam flag_t The integer type used by the flag field
+*/
+//**********************************************************************************************************************
+template< typename Stencil_T, typename flag_t >
+class Dirichlet : public Boundary< flag_t >
+{
+
+   typedef GhostLayerField< real_t, 1 > Field_T;
+   typedef GhostLayerField< real_t, Stencil_T::Size >  StencilField_T;
+
+public:
+
+   static const bool threadsafe = false;
+
+   class DirichletBC : public BoundaryConfiguration {
+   public:
+             DirichletBC( const real_t & _dirichletBC ) : dirichletBC_( _dirichletBC ) {}
+      inline DirichletBC( const Config::BlockHandle & config );
+
+      const real_t & dirichletBC() const { return dirichletBC_; }
+      real_t & dirichletBC() { return dirichletBC_; }
+
+   private:
+
+      real_t dirichletBC_;
+   };
+
+   static shared_ptr<DirichletBC> createConfiguration( const Config::BlockHandle & config ) { return make_shared<DirichletBC>( config ); }
+
+
+
+  //*******************************************************************************************************************
+  /*! Creates a Dirichlet boundary
+   * \param boundaryUID the UID of the boundary condition
+   * \param uid the UID of the flag that marks cells with this boundary condition
+   * \param rhsField pointer to the right-hand side field, which will be adapted by this boundary condition
+   * \param stencilField pointer to the operator stencil field. It should contain the stencil weights that don't take
+   *                     into account the boundary conditions.
+   * \param adaptBCStencilField pointer to the operator stencil field that will be adapted by this boundary condition. 
+   *                            Initially, this field needs to contain the same values as \p stencilField.
+   *                            This is the stencil field that should be passed to the actual PDE solver.
+   * \param flagField pointer to the flag field
+   *******************************************************************************************************************/
+   inline Dirichlet( const BoundaryUID & boundaryUID, const FlagUID & uid, Field_T* const rhsField, const StencilField_T* const stencilField,
+                     StencilField_T* const adaptBCStencilField, FlagField<flag_t> * const flagField );
+
+   void pushFlags( std::vector< FlagUID > & uids ) const { uids.push_back( uid_ ); }
+
+   void beforeBoundaryTreatment() const {}
+   void afterBoundaryTreatment();
+
+   template< typename Buffer_T >
+   inline void packCell( Buffer_T & buffer, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const;
+
+   template< typename Buffer_T >
+   inline void registerCell( Buffer_T & buffer, const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );
+
+   inline void registerCell( const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const BoundaryConfiguration & dirichletBC );
+   inline void registerCells( const flag_t, const CellInterval & cells, const BoundaryConfiguration & dirichletBC );
+   template< typename CellIterator >
+   inline void registerCells( const flag_t, const CellIterator & begin, const CellIterator & end, const BoundaryConfiguration & dirichletBC );
+
+   void unregisterCell( const flag_t, const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz );
+
+   inline void treatDirection( const cell_idx_t  x, const cell_idx_t  y, const cell_idx_t  z, const stencil::Direction dir,
+                               const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz, const flag_t mask );
+
+   inline const real_t & getValue( const cell_idx_t x, cell_idx_t y, cell_idx_t z ) const { return dirichletBC_->get(x,y,z); }
+
+private:
+
+   const FlagUID uid_;
+   const flag_t formerDirichlet_;
+   uint_t numDirtyCells_;
+
+   Field_T* const                rhsField_;
+   shared_ptr< Field_T >         dirichletBC_;
+   const StencilField_T * const  stencilField_;
+   StencilField_T * const        adaptBCStencilField_;
+   FlagField<flag_t> * const     flagField_;
+
+}; // class Dirichlet
+
+
+
+template< typename Stencil_T, typename flag_t >
+inline Dirichlet< Stencil_T, flag_t >::DirichletBC::DirichletBC( const Config::BlockHandle & config  )
+{
+   dirichletBC_ = ( config && config.isDefined( "val" ) ) ? config.getParameter<real_t>( "val" ) : real_c(0.0);
+}
+
+
+
+template< typename Stencil_T, typename flag_t >
+inline Dirichlet< Stencil_T, flag_t >::Dirichlet( const BoundaryUID & boundaryUID, const FlagUID & uid, Field_T* const rhsField, const StencilField_T* const stencilField, StencilField_T* const adaptBCStencilField, FlagField<flag_t> * const flagField ) :
+   Boundary<flag_t>( boundaryUID ), uid_( uid ), formerDirichlet_ (flagField->getOrRegisterFlag("FormerDirichlet")), numDirtyCells_(std::numeric_limits<uint_t>::max()), rhsField_( rhsField ), stencilField_ ( stencilField ), adaptBCStencilField_ ( adaptBCStencilField ), flagField_ (flagField)
+{
+   WALBERLA_ASSERT_NOT_NULLPTR( rhsField_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( stencilField_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( adaptBCStencilField_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( flagField_ );
+
+   WALBERLA_ASSERT_EQUAL( rhsField_->xyzSize(), stencilField_->xyzSize() );
+#ifndef NDEBUG
+   WALBERLA_FOR_ALL_CELLS_XYZ( stencilField_,
+      for( auto dir = Stencil_T::begin(); dir != Stencil_T::end(); ++dir )
+      {
+         WALBERLA_ASSERT_IDENTICAL(stencilField_->get(x,y,z, dir.toIdx()), adaptBCStencilField_->get(x,y,z, dir.toIdx()));
+      }
+   )
+#endif
+
+   dirichletBC_ = make_shared< Field_T >( rhsField_->xSize(), rhsField_->ySize(), rhsField_->zSize(), uint_t(1), field::zyxf );
+
+}
+
+
+template< typename Stencil_T, typename flag_t >
+void Dirichlet< Stencil_T, flag_t >::afterBoundaryTreatment() {
+
+   if (numDirtyCells_>0 && numDirtyCells_ != std::numeric_limits<uint_t>::max()) {
+      WALBERLA_LOG_WARNING("De-registering cells requires re-running Galerkin coarsening");
+   }
+
+   numDirtyCells_=0;
+}
+
+template< typename Stencil_T, typename flag_t >
+template< typename Buffer_T >
+inline void Dirichlet< Stencil_T, flag_t >::packCell( Buffer_T & buffer, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
+{
+   buffer << dirichletBC_->get( x, y, z );
+}
+
+
+
+template< typename Stencil_T, typename flag_t >
+template< typename Buffer_T >
+inline void Dirichlet< Stencil_T, flag_t >::registerCell( Buffer_T & buffer, const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+{
+   buffer >> dirichletBC_->get( x, y, z );
+   if( flagField_->isFlagSet( x, y, z, formerDirichlet_ ) )
+   {
+      flagField_->removeFlag( x, y, z, formerDirichlet_ );
+      --numDirtyCells_;
+   }
+}
+
+
+
+template< typename Stencil_T, typename flag_t >
+inline void Dirichlet< Stencil_T, flag_t >::registerCell( const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+                                                                                       const BoundaryConfiguration & dirichletBC )
+{
+   WALBERLA_ASSERT_EQUAL( dynamic_cast< const DirichletBC * >( &dirichletBC ), &dirichletBC );
+   WALBERLA_ASSERT_NOT_NULLPTR( dirichletBC_ );
+
+   const DirichletBC & val = dynamic_cast< const DirichletBC & >( dirichletBC );
+
+   dirichletBC_->get( x, y, z ) = val.dirichletBC();
+
+   if( flagField_->isFlagSet( x, y, z, formerDirichlet_ ) )
+   {
+      flagField_->removeFlag( x, y, z, formerDirichlet_ );
+      --numDirtyCells_;
+   }
+
+}
+
+
+
+template< typename Stencil_T, typename flag_t >
+inline void Dirichlet< Stencil_T, flag_t >::registerCells( const flag_t, const CellInterval & cells, const BoundaryConfiguration & dirichletBC )
+{
+   WALBERLA_ASSERT_EQUAL( dynamic_cast< const DirichletBC * >( &dirichletBC ), &dirichletBC );
+   WALBERLA_ASSERT_NOT_NULLPTR( dirichletBC_ );
+
+   const DirichletBC & val = dynamic_cast< const DirichletBC & >( dirichletBC );
+
+   for( auto cell = dirichletBC_->beginSliceXYZ( cells ); cell != dirichletBC_->end(); ++cell ) {
+      *cell = val.dirichletBC();
+
+      if( flagField_->isFlagSet( cell.cell(), formerDirichlet_ ) )
+      {
+         flagField_->removeFlag( cell.cell(), formerDirichlet_ );
+         --numDirtyCells_;
+      }
+   }
+}
+
+
+
+template< typename Stencil_T, typename flag_t >
+template< typename CellIterator >
+inline void Dirichlet< Stencil_T, flag_t >::registerCells( const flag_t, const CellIterator & begin, const CellIterator & end,
+                                                                                        const BoundaryConfiguration & dirichletBC )
+{
+   WALBERLA_ASSERT_EQUAL( dynamic_cast< const DirichletBC * >( &dirichletBC ), &dirichletBC );
+   WALBERLA_ASSERT_NOT_NULLPTR( dirichletBC_ );
+
+   const DirichletBC & val = dynamic_cast< const DirichletBC & >( dirichletBC );
+
+   for( auto cell = begin; cell != end; ++cell )
+   {
+      dirichletBC_->get( cell->x(), cell->y(), cell->z() ) = val.dirichletBC();
+
+      if( flagField_->isFlagSet( cell->cell(), formerDirichlet_ ) )
+      {
+         flagField_->removeFlag( cell->cell(), formerDirichlet_ );
+         --numDirtyCells_;
+      }
+   }
+
+}
+
+// Remark: This unregister function works only properly for D3Q7 stencils and convex domains!
+template< typename Stencil_T, typename flag_t >
+void Dirichlet< Stencil_T, flag_t >::unregisterCell( const flag_t, const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz )
+{
+   flagField_->addFlag( nx,ny,nz, formerDirichlet_ );
+   ++numDirtyCells_;
+
+   Cell boundaryCell( nx, ny, nz );
+
+   // Set stencil previously adapted to Dirichlet BC back to un-adapted state
+   for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d )
+   {
+      Cell domainCell = boundaryCell - *d;
+      if( adaptBCStencilField_->isInInnerPart( domainCell ) )
+      {
+//         WALBERLA_LOG_DEVEL("Undo Dirichlet treatment at: " << domainCell );
+
+         // restore original non-center stencil entry of neighboring non-boundary cell
+         adaptBCStencilField_->get( domainCell, d.toIdx() ) = stencilField_->get( domainCell, d.toIdx() );
+
+         // restore original center stencil entry of neighboring non-boundary cell
+         adaptBCStencilField_->get( domainCell, Stencil_T::idx[stencil::C] ) += adaptBCStencilField_->get( domainCell, d.toIdx() );
+      }
+   }
+}
+
+
+//*************************************************************************************************
+/*! \brief Treat one direction by adapting RHS according to Dirichlet boundary value.
+ *
+ */
+
+template< typename Stencil_T, typename flag_t >
+#ifndef NDEBUG
+inline void Dirichlet< Stencil_T, flag_t >::treatDirection( const cell_idx_t  x, const cell_idx_t  y, const cell_idx_t  z, const stencil::Direction dir,
+                                                                                         const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz, const flag_t mask )
+#else
+inline void Dirichlet< Stencil_T, flag_t >::treatDirection( const cell_idx_t  x, const cell_idx_t  y, const cell_idx_t  z, const stencil::Direction dir,
+                                                                                         const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz, const flag_t /*mask*/ )
+#endif
+{
+
+   WALBERLA_ASSERT_EQUAL( nx, x + cell_idx_c( stencil::cx[ dir ] ) );
+   WALBERLA_ASSERT_EQUAL( ny, y + cell_idx_c( stencil::cy[ dir ] ) );
+   WALBERLA_ASSERT_EQUAL( nz, z + cell_idx_c( stencil::cz[ dir ] ) );
+   WALBERLA_ASSERT_UNEQUAL( mask & this->mask_, numeric_cast<flag_t>(0) );
+   WALBERLA_ASSERT_EQUAL( mask & this->mask_, this->mask_ ); // only true if "this->mask_" only contains one single flag, which is the case for the
+                                                             // current implementation of this boundary condition (Dirichlet)
+
+   // Adapt RHS to Dirichlet BC //
+   rhsField_->get( x, y, z ) -= stencilField_->get( x, y, z, Stencil_T::idx[dir] ) * real_t(2) * dirichletBC_->get( nx, ny, nz ); // possibly utilize that off-diagonal entries -1 anyway
+
+   // Adapt Stencils to BCs (former adaptStencilsBC) //
+   // Only required if any new BC cell was added or the BC type of any former BC cell has been changed
+   if (numDirtyCells_>0) {
+
+      // here not thread-safe! (workaround: mutex when adapting central entry)
+      adaptBCStencilField_->get( x, y, z, Stencil_T::idx[stencil::C] ) -= adaptBCStencilField_->get( x, y, z, Stencil_T::idx[dir] );
+      adaptBCStencilField_->get( x, y, z, Stencil_T::idx[dir] ) = 0;
+
+   }
+
+}
+
+
+} // namespace pde
+} // namespace walberla
diff --git a/src/pde/boundary/Dirichlet_withDx.h b/src/pde/boundary/Dirichlet_withDx.h
new file mode 100644
index 0000000000000000000000000000000000000000..c119206926c31891a3d17b73c52e02ec1f92ee3b
--- /dev/null
+++ b/src/pde/boundary/Dirichlet_withDx.h
@@ -0,0 +1,302 @@
+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can 
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT 
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
+//  for more details.
+//  
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file Dirichlet.h
+//! \ingroup pde
+//! \author Dominik Bartuschat <dominik.bartuschat@fau.de>
+//! \author Michael Kuron <mkuron@icp.uni-stuttgart.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "field/Field.h"
+
+#include "boundary/Boundary.h"
+
+#include "core/DataTypes.h"
+#include "core/cell/CellInterval.h"
+#include "core/config/Config.h"
+#include "core/debug/Debug.h"
+#include "core/logging/Logging.h"
+
+#include "field/FlagUID.h"
+
+#include "stencil/Directions.h"
+
+#include <vector>
+#include <array>
+
+
+namespace walberla {
+namespace pde {
+
+
+
+template< typename Stencil_T, typename flag_t >
+class Dirichlet : public Boundary< flag_t >
+{
+
+   typedef GhostLayerField< real_t, 1 > Field_T;
+   typedef GhostLayerField< real_t, Stencil_T::Size >  StencilField_T;
+
+public:
+
+   static const bool threadsafe = false;
+
+   class DirichletBC : public BoundaryConfiguration {
+   public:
+             DirichletBC( const real_t & _dirichletBC ) : dirichletBC_( _dirichletBC ) {}
+      inline DirichletBC( const Config::BlockHandle & config );
+
+      const real_t & dirichletBC() const { return dirichletBC_; }
+      real_t & dirichletBC() { return dirichletBC_; }
+
+   private:
+
+      real_t dirichletBC_;
+   };
+
+   static shared_ptr<DirichletBC> createConfiguration( const Config::BlockHandle & config ) { return make_shared<DirichletBC>( config ); }
+
+
+
+   inline Dirichlet( const BoundaryUID & boundaryUID, const FlagUID & uid, Field_T* const rhsField, const StencilField_T* const stencilField,
+                     StencilField_T* const adaptBCStencilField, FlagField<flag_t> * const flagField, const StructuredBlockStorage& blocks );
+
+   void pushFlags( std::vector< FlagUID > & uids ) const { uids.push_back( uid_ ); }
+
+   void beforeBoundaryTreatment() const {}
+   void afterBoundaryTreatment();
+
+   template< typename Buffer_T >
+   inline void packCell( Buffer_T & buffer, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const;
+
+   template< typename Buffer_T >
+   inline void registerCell( Buffer_T & buffer, const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );
+
+   inline void registerCell( const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const BoundaryConfiguration & dirichletBC );
+   inline void registerCells( const flag_t, const CellInterval & cells, const BoundaryConfiguration & dirichletBC );
+   template< typename CellIterator >
+   inline void registerCells( const flag_t, const CellIterator & begin, const CellIterator & end, const BoundaryConfiguration & dirichletBC );
+
+   void unregisterCell( const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );
+
+   inline void treatDirection( const cell_idx_t  x, const cell_idx_t  y, const cell_idx_t  z, const stencil::Direction dir,
+                               const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz, const flag_t mask );
+
+   inline const real_t & getValue( const cell_idx_t x, cell_idx_t y, cell_idx_t z ) const { return dirichletBC_->get(x,y,z); }
+
+private:
+
+   const FlagUID uid_;
+   const flag_t formerDirichlet_;
+   uint_t numDirtyCells_;
+
+   std::array<real_t, Stencil_T::Q> dx_;
+   Field_T* const                rhsField_;
+   shared_ptr< Field_T >         dirichletBC_;
+   const StencilField_T * const  stencilField_;
+   StencilField_T * const        adaptBCStencilField_;
+   FlagField<flag_t> * const     flagField_;
+
+}; // class Dirichlet
+
+
+
+template< typename Stencil_T, typename flag_t >
+inline Dirichlet< Stencil_T, flag_t >::DirichletBC::DirichletBC( const Config::BlockHandle & config  )
+{
+   dirichletBC_ = ( config && config.isDefined( "val" ) ) ? config.getParameter<real_t>( "val" ) : real_c(0.0);
+}
+
+
+
+template< typename Stencil_T, typename flag_t >
+inline Dirichlet< Stencil_T, flag_t >::Dirichlet( const BoundaryUID & boundaryUID, const FlagUID & uid, Field_T* const rhsField, const StencilField_T* const stencilField, StencilField_T* const adaptBCStencilField, FlagField<flag_t> * const flagField, const StructuredBlockStorage& blocks ) :
+   Boundary<flag_t>( boundaryUID ), uid_( uid ), formerDirichlet_ (flagField->getOrRegisterFlag("FormerDirichlet")), numDirtyCells_(0), rhsField_( rhsField ), stencilField_ ( stencilField ), adaptBCStencilField_ ( adaptBCStencilField ), flagField_ (flagField)
+{
+   WALBERLA_ASSERT_NOT_NULLPTR( rhsField_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( stencilField_ );
+
+   WALBERLA_ASSERT_EQUAL( rhsField_->xyzSize(), stencilField_->xyzSize() );
+
+   dirichletBC_ = make_shared< Field_T >( rhsField_->xSize(), rhsField_->ySize(), rhsField_->zSize(), uint_t(1), field::zyxf );
+
+   for(auto d = Stencil_T::beginNoCenter(); d != Stencil_T::end(); ++d ){
+      dx_[d.toIdx()] = Vector3<real_t>(stencil::cx[d.toIdx()]*blocks.dx(), stencil::cy[d.toIdx()]*blocks.dy(), stencil::cz[d.toIdx()]*blocks.dz() ).sqrLength();
+      WALBERLA_LOG_DEVEL("dx in direction " << d.dirString() << ":" << dx_[d.toIdx()]);
+   }
+
+
+}
+
+
+template< typename Stencil_T, typename flag_t >
+void Dirichlet< Stencil_T, flag_t >::afterBoundaryTreatment() {
+
+   if (numDirtyCells_>0) {
+      WALBERLA_LOG_WARNING("De-registering cells requires re-running Galerkin coarsening");
+   }
+
+   numDirtyCells_=0;
+}
+
+template< typename Stencil_T, typename flag_t >
+template< typename Buffer_T >
+inline void Dirichlet< Stencil_T, flag_t >::packCell( Buffer_T & buffer, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
+{
+   buffer << dirichletBC_->get( x, y, z );
+}
+
+
+
+template< typename Stencil_T, typename flag_t >
+template< typename Buffer_T >
+inline void Dirichlet< Stencil_T, flag_t >::registerCell( Buffer_T & buffer, const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+{
+   buffer >> dirichletBC_->get( x, y, z );
+   if( flagField_->isFlagSet( x, y, z, formerDirichlet_ ) )
+   {
+      flagField_->removeFlag( x, y, z, formerDirichlet_ );
+      --numDirtyCells_;
+   }
+}
+
+
+
+template< typename Stencil_T, typename flag_t >
+inline void Dirichlet< Stencil_T, flag_t >::registerCell( const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+                                                                                       const BoundaryConfiguration & dirichletBC )
+{
+   WALBERLA_ASSERT_EQUAL( dynamic_cast< const DirichletBC * >( &dirichletBC ), &dirichletBC );
+   WALBERLA_ASSERT_NOT_NULLPTR( dirichletBC_ );
+
+   const DirichletBC & val = dynamic_cast< const DirichletBC & >( dirichletBC );
+
+   dirichletBC_->get( x, y, z ) = val.dirichletBC();
+
+   if( flagField_->isFlagSet( x, y, z, formerDirichlet_ ) )
+   {
+      flagField_->removeFlag( x, y, z, formerDirichlet_ );
+      --numDirtyCells_;
+   }
+
+}
+
+
+
+template< typename Stencil_T, typename flag_t >
+inline void Dirichlet< Stencil_T, flag_t >::registerCells( const flag_t, const CellInterval & cells, const BoundaryConfiguration & dirichletBC )
+{
+   WALBERLA_ASSERT_EQUAL( dynamic_cast< const DirichletBC * >( &dirichletBC ), &dirichletBC );
+   WALBERLA_ASSERT_NOT_NULLPTR( dirichletBC_ );
+
+   const DirichletBC & val = dynamic_cast< const DirichletBC & >( dirichletBC );
+
+   for( auto cell = dirichletBC_->beginSliceXYZ( cells ); cell != dirichletBC_->end(); ++cell ) {
+      *cell = val.dirichletBC();
+
+      if( flagField_->isFlagSet( cell.cell(), formerDirichlet_ ) )
+      {
+         flagField_->removeFlag( cell.cell(), formerDirichlet_ );
+         --numDirtyCells_;
+      }
+   }
+}
+
+
+
+template< typename Stencil_T, typename flag_t >
+template< typename CellIterator >
+inline void Dirichlet< Stencil_T, flag_t >::registerCells( const flag_t, const CellIterator & begin, const CellIterator & end,
+                                                                                        const BoundaryConfiguration & dirichletBC )
+{
+   WALBERLA_ASSERT_EQUAL( dynamic_cast< const DirichletBC * >( &dirichletBC ), &dirichletBC );
+   WALBERLA_ASSERT_NOT_NULLPTR( dirichletBC_ );
+
+   const DirichletBC & val = dynamic_cast< const DirichletBC & >( dirichletBC );
+
+   for( auto cell = begin; cell != end; ++cell )
+   {
+      dirichletBC_->get( cell->x(), cell->y(), cell->z() ) = val.dirichletBC();
+
+      if( flagField_->isFlagSet( cell->cell(), formerDirichlet_ ) )
+      {
+         flagField_->removeFlag( cell->cell(), formerDirichlet_ );
+         --numDirtyCells_;
+      }
+   }
+
+}
+
+
+template< typename Stencil_T, typename flag_t >
+void Dirichlet< Stencil_T, flag_t >::unregisterCell( const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+{
+   flagField_->addFlag( x,y,z, formerDirichlet_ );
+   ++numDirtyCells_;
+
+   // Set stencil adapted to BCs back to unadapted state
+   for(auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d ){
+      adaptBCStencilField_->get(x,y,z,d.toIdx()) = stencilField_->get(x,y,z,d.toIdx());
+   }
+
+}
+
+
+//*************************************************************************************************
+/*! \brief Treat one direction by adapting RHS according to Dirichlet boundary value.
+ *
+ */
+
+template< typename Stencil_T, typename flag_t >
+#ifndef NDEBUG
+inline void Dirichlet< Stencil_T, flag_t >::treatDirection( const cell_idx_t  x, const cell_idx_t  y, const cell_idx_t  z, const stencil::Direction dir,
+                                                                                         const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz, const flag_t mask )
+#else
+inline void Dirichlet< Stencil_T, flag_t >::treatDirection( const cell_idx_t  x, const cell_idx_t  y, const cell_idx_t  z, const stencil::Direction dir,
+                                                                                         const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz, const flag_t /*mask*/ )
+#endif
+{
+
+   WALBERLA_ASSERT_EQUAL( nx, x + cell_idx_c( stencil::cx[ dir ] ) );
+   WALBERLA_ASSERT_EQUAL( ny, y + cell_idx_c( stencil::cy[ dir ] ) );
+   WALBERLA_ASSERT_EQUAL( nz, z + cell_idx_c( stencil::cz[ dir ] ) );
+   WALBERLA_ASSERT_UNEQUAL( mask & this->mask_, numeric_cast<flag_t>(0) );
+   WALBERLA_ASSERT_EQUAL( mask & this->mask_, this->mask_ ); // only true if "this->mask_" only contains one single flag, which is the case for the
+                                                             // current implementation of this boundary condition (Dirichlet)
+
+   // Adapt RHS to Dirichlet BC //
+   rhsField_->get( x, y, z ) -= stencilField_->get( x, y, z, Stencil_T::idx[dir] ) * real_t(2) * dx_[ Stencil_T::idx[dir] ] * dirichletBC_->get( nx, ny, nz ); // possibly utilize that off-diagonal entries -1 anyway
+
+   // WALBERLA_LOG_DEVEL("Adapt RHS to Dirichlet value " << dirichletBC_->get( nx, ny, nz ) << " on cell " << Cell(x,y,z) << " for stencil entry " << stencilField_->get( x, y, z, Stencil_T::idx[dir] ) );
+
+   // Adapt Stencils to BCs (former adaptStencilsBC) //
+   // Only required if any new BC cell was added or the BC type of any former BC cell has been changed
+   if (numDirtyCells_>0) {
+
+      // here not thread-safe! (workaround: mutex when adapting central entry)
+      adaptBCStencilField_->get( x, y, z, Stencil_T::idx[stencil::C] ) -= adaptBCStencilField_->get( x, y, z, Stencil_T::idx[dir] );
+      adaptBCStencilField_->get( x, y, z, Stencil_T::idx[dir] ) = 0;
+
+   }
+
+}
+
+
+
+} // namespace pde
+} // namespace walberla
diff --git a/src/pde/boundary/Neumann.h b/src/pde/boundary/Neumann.h
index 42a23295512300ba83b74fe4c38ec19095b534da..ba5fb0e8c7e09ee039a793aee5f1623a10053250 100644
--- a/src/pde/boundary/Neumann.h
+++ b/src/pde/boundary/Neumann.h
@@ -15,16 +15,34 @@
 //
 //! \file Neumann.h
 //! \ingroup pde
+//! \author Dominik Bartuschat <dominik.bartuschat@fau.de>
+//! \author Michael Kuron <mkuron@icp.uni-stuttgart.de>
 //! \author Florian Schornbaum <florian.schornbaum@fau.de>
 //
 //======================================================================================================================
 
 #pragma once
 
+#include "field/Field.h"
+
+#include "boundary/Boundary.h"
+
+#include "core/DataTypes.h"
 #include "core/cell/CellInterval.h"
+#include "core/config/Config.h"
+#include "core/debug/Debug.h"
+#include "core/logging/Logging.h"
+
 #include "domain_decomposition/StructuredBlockStorage.h"
+
+#include "field/FlagUID.h"
+
+#include "stencil/Directions.h"
 #include "stencil/D3Q6.h"
 
+#include <vector>
+#include <limits>
+#include <array>
 
 
 namespace walberla {
@@ -207,5 +225,305 @@ void NeumannDomainBoundary< PdeField >::apply( PdeField * p, const CellInterval
 
 
 
+//**********************************************************************************************************************
+/*!
+*   \brief Neumann boundary handling for PDEs
+*
+*   This boundary condition imposes a Neumann condition with arbitrary values on a PDE.
+*   It does so by modifying the right-hand side and the stencil field.
+*   Anything that has internal copies of the stencil field (e.g. the multigrid V-cycle's coarsened stencils) is
+*   responsible for updating its copies when boundary conditions are changed.
+*
+*   \tparam Stencil_T The stencil used for the discrete operator
+*   \tparam flag_t The integer type used by the flag field
+*/
+//**********************************************************************************************************************
+template< typename Stencil_T, typename flag_t >
+class Neumann : public Boundary< flag_t >
+{
+
+   typedef GhostLayerField< real_t, 1 > Field_T;
+   typedef GhostLayerField< real_t, Stencil_T::Size >  StencilField_T;
+
+public:
+
+   static const bool threadsafe = false;
+
+   class NeumannBC : public BoundaryConfiguration {
+   public:
+             NeumannBC( const real_t & _neumannBC ) : neumannBC_( _neumannBC ) {}
+      inline NeumannBC( const Config::BlockHandle & config );
+
+      const real_t & neumannBC() const { return neumannBC_; }
+      real_t & neumannBC() { return neumannBC_; }
+
+   private:
+
+      real_t neumannBC_;
+   };
+
+   static shared_ptr<NeumannBC> createConfiguration( const Config::BlockHandle & config ) { return make_shared<NeumannBC>( config ); }
+
+
+  //*******************************************************************************************************************
+  /*! Creates a Neumann boundary
+   * \param boundaryUID the UID of the boundary condition
+   * \param uid the UID of the flag that marks cells with this boundary condition
+   * \param rhsField pointer to the right-hand side field, which will be adapted by this boundary condition
+   * \param stencilField pointer to the operator stencil field. It should contain the stencil weights that don't take
+   *                     into account the boundary conditions.
+   * \param adaptBCStencilField pointer to the operator stencil field that will be adapted by this boundary condition. 
+   *                            Initially, this field needs to contain the same values as \p stencilField.
+   *                            This is the stencil field that should be passed to the actual PDE solver.
+   * \param flagField pointer to the flag field
+   * \param blocks
+   *******************************************************************************************************************/
+   inline Neumann( const BoundaryUID & boundaryUID, const FlagUID & uid, Field_T* const rhsField, const StencilField_T* const stencilField,
+                     StencilField_T* const adaptBCStencilField, FlagField<flag_t> * const flagField, const StructuredBlockStorage& blocks );
+
+   void pushFlags( std::vector< FlagUID > & uids ) const { uids.push_back( uid_ ); }
+
+   void beforeBoundaryTreatment() const {}
+   void afterBoundaryTreatment();
+
+   template< typename Buffer_T >
+   inline void packCell( Buffer_T & buffer, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const;
+
+   template< typename Buffer_T >
+   inline void registerCell( Buffer_T & buffer, const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );
+
+   inline void registerCell( const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const BoundaryConfiguration & neumannBC );
+   inline void registerCells( const flag_t, const CellInterval & cells, const BoundaryConfiguration & neumannBC );
+   template< typename CellIterator >
+   inline void registerCells( const flag_t, const CellIterator & begin, const CellIterator & end, const BoundaryConfiguration & neumannBC );
+
+   void unregisterCell( const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );
+
+   inline void treatDirection( const cell_idx_t  x, const cell_idx_t  y, const cell_idx_t  z, const stencil::Direction dir,
+                               const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz, const flag_t mask );
+
+   inline const real_t & getValue( const cell_idx_t x, cell_idx_t y, cell_idx_t z ) const { return neumannBC_->get(x,y,z); }
+
+private:
+
+   const FlagUID uid_;
+   const flag_t formerNeumann_;
+   uint_t numDirtyCells_;
+
+   std::array<real_t, Stencil_T::Q> dx_;
+   Field_T* const                rhsField_;
+   shared_ptr< Field_T >         neumannBC_;
+   const StencilField_T * const  stencilField_;
+   StencilField_T * const        adaptBCStencilField_;
+   FlagField<flag_t> * const     flagField_;
+
+}; // class Neumann
+
+
+
+template< typename Stencil_T, typename flag_t >
+inline Neumann< Stencil_T, flag_t >::NeumannBC::NeumannBC( const Config::BlockHandle & config  )
+{
+   neumannBC_ = ( config && config.isDefined( "val" ) ) ? config.getParameter<real_t>( "val" ) : real_c(0.0);
+}
+
+
+
+template< typename Stencil_T, typename flag_t >
+inline Neumann< Stencil_T, flag_t >::Neumann( const BoundaryUID & boundaryUID, const FlagUID & uid, Field_T* const rhsField, const StencilField_T* const stencilField,
+                                              StencilField_T* const adaptBCStencilField, FlagField<flag_t> * const flagField, const StructuredBlockStorage& blocks  )
+                                              : Boundary<flag_t>( boundaryUID ), uid_( uid ), formerNeumann_ (flagField->getOrRegisterFlag("FormerNeumann")), numDirtyCells_(std::numeric_limits<uint_t>::max()),
+                                                rhsField_( rhsField ), stencilField_ ( stencilField ), adaptBCStencilField_ ( adaptBCStencilField ), flagField_ (flagField)
+{
+   WALBERLA_ASSERT_NOT_NULLPTR( rhsField_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( stencilField_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( adaptBCStencilField_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( flagField_ );
+
+   WALBERLA_ASSERT_EQUAL( rhsField_->xyzSize(), stencilField_->xyzSize() );
+#ifndef NDEBUG
+   WALBERLA_FOR_ALL_CELLS_XYZ( stencilField_,
+      for( auto dir = Stencil_T::begin(); dir != Stencil_T::end(); ++dir )
+      {
+         WALBERLA_ASSERT_IDENTICAL(stencilField_->get(x,y,z, dir.toIdx()), adaptBCStencilField_->get(x,y,z, dir.toIdx()));
+      }
+   )
+#endif
+
+   neumannBC_ = make_shared< Field_T >( rhsField_->xSize(), rhsField_->ySize(), rhsField_->zSize(), uint_t(1), field::zyxf );
+
+   for(auto d = Stencil_T::beginNoCenter(); d != Stencil_T::end(); ++d ){
+      dx_[d.toIdx()] = Vector3<real_t>(real_c(stencil::cx[d.toIdx()])*blocks.dx(), real_c(stencil::cy[d.toIdx()])*blocks.dy(), real_c(stencil::cz[d.toIdx()])*blocks.dz() ).sqrLength();
+      // WALBERLA_LOG_DEVEL("dx in direction " << d.dirString() << ":" << dx_[d.toIdx()]);
+   }
+
+}
+
+
+template< typename Stencil_T, typename flag_t >
+void Neumann< Stencil_T, flag_t >::afterBoundaryTreatment() {
+
+   if (numDirtyCells_>0 && numDirtyCells_ != std::numeric_limits<uint_t>::max()) {
+      WALBERLA_LOG_WARNING("De-registering cells requires re-running Galerkin coarsening");
+   }
+
+   numDirtyCells_=0;
+}
+
+template< typename Stencil_T, typename flag_t >
+template< typename Buffer_T >
+inline void Neumann< Stencil_T, flag_t >::packCell( Buffer_T & buffer, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
+{
+   buffer << neumannBC_->get( x, y, z );
+}
+
+
+
+template< typename Stencil_T, typename flag_t >
+template< typename Buffer_T >
+inline void Neumann< Stencil_T, flag_t >::registerCell( Buffer_T & buffer, const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+{
+   buffer >> neumannBC_->get( x, y, z );
+   if( flagField_->isFlagSet( x, y, z, formerNeumann_ ) )
+   {
+      flagField_->removeFlag( x, y, z, formerNeumann_ );
+      --numDirtyCells_;
+   }
+}
+
+
+
+template< typename Stencil_T, typename flag_t >
+inline void Neumann< Stencil_T, flag_t >::registerCell( const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+                                                                                       const BoundaryConfiguration & neumannBC )
+{
+   WALBERLA_ASSERT_EQUAL( dynamic_cast< const NeumannBC * >( &neumannBC ), &neumannBC );
+   WALBERLA_ASSERT_NOT_NULLPTR( neumannBC_ );
+
+   const NeumannBC & val = dynamic_cast< const NeumannBC & >( neumannBC );
+
+   neumannBC_->get( x, y, z ) = val.neumannBC();
+
+   if( flagField_->isFlagSet( x, y, z, formerNeumann_ ) )
+   {
+      flagField_->removeFlag( x, y, z, formerNeumann_ );
+      --numDirtyCells_;
+   }
+
+}
+
+
+
+template< typename Stencil_T, typename flag_t >
+inline void Neumann< Stencil_T, flag_t >::registerCells( const flag_t, const CellInterval & cells, const BoundaryConfiguration & neumannBC )
+{
+   WALBERLA_ASSERT_EQUAL( dynamic_cast< const NeumannBC * >( &neumannBC ), &neumannBC );
+   WALBERLA_ASSERT_NOT_NULLPTR( neumannBC_ );
+
+   const NeumannBC & val = dynamic_cast< const NeumannBC & >( neumannBC );
+
+   for( auto cell = neumannBC_->beginSliceXYZ( cells ); cell != neumannBC_->end(); ++cell ) {
+      *cell = val.neumannBC();
+
+      if( flagField_->isFlagSet( cell.cell(), formerNeumann_ ) )
+      {
+         flagField_->removeFlag( cell.cell(), formerNeumann_ );
+         --numDirtyCells_;
+      }
+   }
+}
+
+
+
+template< typename Stencil_T, typename flag_t >
+template< typename CellIterator >
+inline void Neumann< Stencil_T, flag_t >::registerCells( const flag_t, const CellIterator & begin, const CellIterator & end,
+                                                                                        const BoundaryConfiguration & neumannBC )
+{
+   WALBERLA_ASSERT_EQUAL( dynamic_cast< const NeumannBC * >( &neumannBC ), &neumannBC );
+   WALBERLA_ASSERT_NOT_NULLPTR( neumannBC_ );
+
+   const NeumannBC & val = dynamic_cast< const NeumannBC & >( neumannBC );
+
+   for( auto cell = begin; cell != end; ++cell )
+   {
+      neumannBC_->get( cell->x(), cell->y(), cell->z() ) = val.neumannBC();
+
+      if( flagField_->isFlagSet( cell->cell(), formerNeumann_ ) )
+      {
+         flagField_->removeFlag( cell->cell(), formerNeumann_ );
+         --numDirtyCells_;
+      }
+   }
+
+}
+
+// Remark: This unregister function works only properly for D3Q7 stencils and convex domains!
+template< typename Stencil_T, typename flag_t >
+void Neumann< Stencil_T, flag_t >::unregisterCell( const flag_t, const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz )
+{
+   flagField_->addFlag( nx, ny, nz, formerNeumann_ );
+   ++numDirtyCells_;
+
+   Cell boundaryCell( nx, ny, nz );
+
+   // Set stencil previously adapted to Neumann BC back to un-adapted state
+   for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d )
+   {
+      Cell domainCell = boundaryCell - *d;
+      if( adaptBCStencilField_->isInInnerPart( domainCell ) )
+      {
+         // restore original non-center stencil entry of neighboring non-boundary cell
+         adaptBCStencilField_->get( domainCell, d.toIdx() ) = stencilField_->get( domainCell, d.toIdx() );
+
+         // restore original center stencil entry of neighboring non-boundary cell
+         adaptBCStencilField_->get( domainCell, Stencil_T::idx[stencil::C] ) -=  adaptBCStencilField_->get( domainCell, d.toIdx() );
+
+      }
+   }
+
+}
+
+
+//*************************************************************************************************
+/*! \brief Treat one direction by adapting RHS according to Neumann boundary value.
+ *
+ */
+
+template< typename Stencil_T, typename flag_t >
+#ifndef NDEBUG
+inline void Neumann< Stencil_T, flag_t >::treatDirection( const cell_idx_t  x, const cell_idx_t  y, const cell_idx_t  z, const stencil::Direction dir,
+                                                                                         const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz, const flag_t mask )
+#else
+inline void Neumann< Stencil_T, flag_t >::treatDirection( const cell_idx_t  x, const cell_idx_t  y, const cell_idx_t  z, const stencil::Direction dir,
+                                                                                         const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz, const flag_t /*mask*/ )
+#endif
+{
+
+   WALBERLA_ASSERT_EQUAL( nx, x + cell_idx_c( stencil::cx[ dir ] ) );
+   WALBERLA_ASSERT_EQUAL( ny, y + cell_idx_c( stencil::cy[ dir ] ) );
+   WALBERLA_ASSERT_EQUAL( nz, z + cell_idx_c( stencil::cz[ dir ] ) );
+   WALBERLA_ASSERT_UNEQUAL( mask & this->mask_, numeric_cast<flag_t>(0) );
+   WALBERLA_ASSERT_EQUAL( mask & this->mask_, this->mask_ ); // only true if "this->mask_" only contains one single flag, which is the case for the
+                                                             // current implementation of this boundary condition (Neumann)
+
+   // WALBERLA_LOG_DEVEL("Neumann treatment at: " << Cell(x,y,z));
+
+   // Adapt RHS to Neumann BC //
+   rhsField_->get( x, y, z ) -= stencilField_->get( x, y, z, Stencil_T::idx[dir] ) * dx_[ Stencil_T::idx[dir] ] * neumannBC_->get( nx, ny, nz ); // possibly utilize that off-diagonal entries -1 anyway
+
+   // Adapt Stencils to BCs (former adaptStencilsBC) //
+   // Only required if any new BC cell was added or the BC type of any former BC cell has been changed
+   if (numDirtyCells_>0) {
+
+      // here not thread-safe! (workaround: mutex when adapting central entry)
+      adaptBCStencilField_->get( x, y, z, Stencil_T::idx[stencil::C] ) += adaptBCStencilField_->get( x, y, z, Stencil_T::idx[dir] );
+      adaptBCStencilField_->get( x, y, z, Stencil_T::idx[dir] ) = 0;
+
+   }
+
+}
+
+
 } // namespace pde
 } // namespace walberla
diff --git a/src/pde/boundary/all.h b/src/pde/boundary/all.h
index 16d35f8642b56d4147c7dd6472234fcf343e433d..9e392663e08bda508eab068879508b8c3579d5b3 100644
--- a/src/pde/boundary/all.h
+++ b/src/pde/boundary/all.h
@@ -22,4 +22,5 @@
 
 #pragma once
 
+#include "Dirichlet.h"
 #include "Neumann.h"
diff --git a/tests/core/CMakeLists.txt b/tests/core/CMakeLists.txt
index ddb9907147e78dbba1d6a9d5fd07be6f2dc033e6..0a000974320fd66a380ac5928e7eafac43d0d64a 100644
--- a/tests/core/CMakeLists.txt
+++ b/tests/core/CMakeLists.txt
@@ -183,6 +183,9 @@ waLBerla_execute_test( NAME DataTypesTest )
 waLBerla_compile_test( FILES GridGeneratorTest.cpp )
 waLBerla_execute_test( NAME GridGeneratorTest )
 
+waLBerla_compile_test( FILES OpenMPWrapperTest.cpp )
+waLBerla_execute_test( NAME OpenMPWrapperTest )
+
 waLBerla_compile_test( FILES SetTest.cpp )
 waLBerla_execute_test( NAME SetTest )
 
diff --git a/tests/core/OpenMPWrapperTest.cpp b/tests/core/OpenMPWrapperTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6da8cbadca525c6121157472d529c6d3c3875e8f
--- /dev/null
+++ b/tests/core/OpenMPWrapperTest.cpp
@@ -0,0 +1,40 @@
+//======================================================================================================================
+//
+//  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 OpenMPWrapperTest.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include <core/debug/TestSubsystem.h>
+#include <core/logging/Logging.h>
+#include <core/OpenMP.h>
+
+using namespace walberla;
+
+int main( int /*argc*/, char** /*argv*/ )
+{
+   debug::enterTestMode();
+   
+   WALBERLA_OPENMP_SECTION()
+   {
+      WALBERLA_LOG_DEVEL_VAR(omp_get_num_threads  ());
+      WALBERLA_LOG_DEVEL_VAR(omp_get_dynamic      ());
+      WALBERLA_LOG_DEVEL_VAR(omp_get_nested       ());
+      WALBERLA_LOG_DEVEL_VAR(omp_get_max_threads  ());
+      WALBERLA_LOG_DEVEL_VAR(omp_get_thread_num   ());
+      WALBERLA_LOG_DEVEL_VAR(omp_get_num_procs    ());
+   }
+}
diff --git a/tests/core/mpi/BufferSystemTest.cpp b/tests/core/mpi/BufferSystemTest.cpp
index a57d24123ed8269d30e71112f4780ce72737b8e6..d84d12da34e0ee84143d5ae1d7f59442f2bec2a6 100644
--- a/tests/core/mpi/BufferSystemTest.cpp
+++ b/tests/core/mpi/BufferSystemTest.cpp
@@ -110,6 +110,9 @@ void symmetricCommunication()
 
       WALBERLA_CHECK_EQUAL( receivedVal, it.rank() );
    }
+
+   WALBERLA_CHECK_EQUAL( bs.getBytesSent(), (MSG_SIZE * sizeof(int) + MSG_SIZE * mpi::BUFFER_DEBUG_OVERHEAD) * 2 );
+   WALBERLA_CHECK_EQUAL( bs.getBytesReceived(), (MSG_SIZE * sizeof(int) + MSG_SIZE * mpi::BUFFER_DEBUG_OVERHEAD) * 2 );
 }
 
 /**
@@ -175,6 +178,9 @@ void asymmetricCommunication()
          WALBERLA_CHECK( it.buffer().isEmpty() );
       }
    }
+
+   WALBERLA_CHECK_EQUAL( bs.getBytesSent(), int64_c(sizeof(int) + mpi::BUFFER_DEBUG_OVERHEAD) * int64_c(rank + rank) );
+   WALBERLA_CHECK_EQUAL( bs.getBytesReceived(), int64_c(sizeof(int) + mpi::BUFFER_DEBUG_OVERHEAD) * int64_c(leftNeighbor + rightNeighbor) );
 }
 
 
diff --git a/tests/pde/BoundaryTest.cpp b/tests/pde/BoundaryTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..17e6b05a3815621c94360ecc2700bd9c5be7a512
--- /dev/null
+++ b/tests/pde/BoundaryTest.cpp
@@ -0,0 +1,445 @@
+//======================================================================================================================
+//
+//  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 BoundaryTest.cpp
+//! \ingroup pde
+//! \author Dominik Bartuschat <dominik.bartuschat@fau.de>
+//! \author Michael Kuron <mkuron@icp.uni-stuttgart.de>
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+#include "boundary/BoundaryHandling.h"
+
+#include "core/Abort.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/mpi/Environment.h"
+#include "core/mpi/MPIManager.h"
+
+#include "field/AddToStorage.h"
+#include "field/GhostLayerField.h"
+#include "field/FlagField.h"
+#include "field/communication/PackInfo.h"
+#include "field/vtk/VTKWriter.h"
+
+#include "pde/iterations/CGFixedStencilIteration.h"
+#include "pde/iterations/CGIteration.h"
+#include "pde/boundary/Dirichlet.h"
+#include "pde/boundary/Neumann.h"
+
+#include "stencil/D2Q5.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+#include "vtk/VTKOutput.h"
+
+#include <cmath>
+
+namespace walberla {
+
+
+
+typedef GhostLayerField< real_t, 1 > Field_T;
+typedef stencil::D2Q5                Stencil_T;
+typedef pde::CGIteration<Stencil_T>::StencilField_T  StencilField_T;
+
+typedef walberla::uint8_t      flag_t;
+typedef FlagField < flag_t >   FlagField_T;
+typedef pde::Dirichlet< Stencil_T, flag_t >  Dirichlet_T;
+typedef pde::Neumann< Stencil_T, flag_t >  Neumann_T;
+
+typedef boost::tuples::tuple< Dirichlet_T, Neumann_T >  BoundaryConditions_T;
+
+typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+
+
+const FlagUID  Domain_Flag( "domain" );
+const FlagUID  Dirichlet_Flag( "dirichlet" );
+const FlagUID  Neumann_Flag( "neumann" );
+
+
+
+///////////////////////
+// BOUNDARY HANDLING //
+///////////////////////
+
+//**********************************************************************************************************************
+/*!
+*   \brief Functor that is used to add a boundary handling object to each block
+*
+*   The member function "BoundaryHandling_T * operator()( IBlock* const block ) const" is called by the framework in
+*   order to add a boundary handling object of type 'BoundaryHandling_T' to each block.
+*/
+//**********************************************************************************************************************
+
+class MyBoundaryHandling
+{
+public:
+
+   MyBoundaryHandling( const BlockDataID & rhsFieldId, const BlockDataID & stencilFieldId, const BlockDataID & adaptBCStencilFieldId, const BlockDataID & flagFieldId, const StructuredBlockStorage& blocks ) :
+      rhsFieldId_( rhsFieldId ), stencilFieldId_ ( stencilFieldId ), adaptBCStencilFieldId_ ( adaptBCStencilFieldId ), flagFieldId_ (flagFieldId), blockStorage_(blocks) {}
+
+   BoundaryHandling_T * operator()( IBlock* const block ) const;
+
+private:
+
+   const BlockDataID rhsFieldId_;
+   const BlockDataID stencilFieldId_;
+   const BlockDataID adaptBCStencilFieldId_;
+   const BlockDataID flagFieldId_;
+   const StructuredBlockStorage & blockStorage_;
+
+}; // class MyBoundaryHandling
+
+BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block ) const
+{
+
+   Field_T * rhsField = block->getData< Field_T >( rhsFieldId_ );
+   StencilField_T * stencilField = block->getData< StencilField_T >( stencilFieldId_ );
+   StencilField_T * adaptStencilField = block->getData< StencilField_T >( adaptBCStencilFieldId_ );
+   FlagField_T * flagField = block->getData< FlagField_T >( flagFieldId_ );
+
+   const flag_t domain = flagField->registerFlag( Domain_Flag ); // register the fluid flag at the flag field
+
+   // A new boundary handling instance that uses the just fetched flag field is created:
+   // Additional, internal flags used by the boundary handling will be stored in this flag field.
+   return new BoundaryHandling_T( "boundary handling", flagField, domain,
+         boost::tuples::make_tuple( Dirichlet_T( "dirichlet", Dirichlet_Flag, rhsField, stencilField, adaptStencilField, flagField ) ,
+                                    Neumann_T( "neumann", Neumann_Flag, rhsField, stencilField, adaptStencilField, flagField, blockStorage_ )
+         )
+   );
+}
+
+
+
+void initRHS( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDataID & rhsId )
+{
+   for( auto block = blocks->begin(); block != blocks->end(); ++block )
+   {
+      Field_T * rhs = block->getData< Field_T >( rhsId );
+      CellInterval xyz = rhs->xyzSize();
+      for( auto cell = xyz.begin(); cell != xyz.end(); ++cell )
+      {
+         const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
+         rhs->get( *cell ) = real_t(4) * math::PI * math::PI * std::sin( real_t(2) * math::PI * p[0] ) * std::sinh( real_t(2) * math::PI * p[1] );
+      }
+   }
+}
+
+
+
+void setRHSConstValue( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDataID & rhsId, const real_t value )
+{
+   for( auto block = blocks->begin(); block != blocks->end(); ++block )
+   {
+      Field_T * rhs = block->getData< Field_T >( rhsId );
+      CellInterval xyz = rhs->xyzSize();
+      for( auto cell = xyz.begin(); cell != xyz.end(); ++cell )
+      {
+         rhs->get( *cell ) = value;
+      }
+   }
+}
+
+
+
+void setBoundaryConditionsDirichl( shared_ptr< StructuredBlockForest > & blocks, const BlockDataID & boundaryHandlingId )
+{
+   CellInterval domainBBInGlobalCellCoordinates = blocks->getDomainCellBB();
+   domainBBInGlobalCellCoordinates.expand(Cell(1,1,0));
+
+   // Iterate through all blocks that are allocated on this process and part of the block structure 'blocks'
+   for( auto block = blocks->begin(); block != blocks->end(); ++block )
+   {
+
+      BoundaryHandling_T * boundaryHandling = block->getData< BoundaryHandling_T >( boundaryHandlingId );
+
+      CellInterval domainBB( domainBBInGlobalCellCoordinates );
+      blocks->transformGlobalToBlockLocalCellInterval( domainBB, *block );
+
+      CellInterval west( domainBB.xMin(), domainBB.yMin(), domainBB.zMin(), domainBB.xMin(), domainBB.yMax(), domainBB.zMax() );
+      CellInterval east( domainBB.xMax(), domainBB.yMin(), domainBB.zMin(), domainBB.xMax(), domainBB.yMax(), domainBB.zMax() );
+      CellInterval south( domainBB.xMin(), domainBB.yMin(), domainBB.zMin(), domainBB.xMax(), domainBB.yMin(), domainBB.zMax() );
+      CellInterval north( domainBB.xMin(), domainBB.yMax(), domainBB.zMin(), domainBB.xMax(), domainBB.yMax(), domainBB.zMax() );
+
+      // Set north boundary to defined function
+      for( auto cell = north.begin(); cell != north.end(); ++cell )
+      {
+
+         const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
+         real_t val = std::sin( real_t( 2 ) * math::PI * p[0] ) * std::sinh( real_t( 2 ) * math::PI * p[1] );
+
+         boundaryHandling->forceBoundary( Dirichlet_Flag, cell->x(), cell->y(), cell->z(), pde::Dirichlet< Stencil_T, flag_t >::DirichletBC( val ) );
+
+      }
+
+      // Set all other boundaries to zero
+      for( auto & ci : { west, east, south } )
+      {
+         boundaryHandling->forceBoundary( Dirichlet_Flag, ci, pde::Dirichlet< Stencil_T, flag_t >::DirichletBC( real_t( 0 ) ) );
+      }
+
+      // All the remaining cells need to be marked as being fluid. The 'fillWithDomain' call does just that.
+      boundaryHandling->fillWithDomain( domainBB );
+   }
+}
+
+
+
+void setBoundaryConditionsMixed( shared_ptr< StructuredBlockForest > & blocks, const BlockDataID & boundaryHandlingId )
+{
+   CellInterval domainBBInGlobalCellCoordinates = blocks->getDomainCellBB();
+   domainBBInGlobalCellCoordinates.expand(Cell(1,1,0));
+
+   // Iterate through all blocks that are allocated on this process and part of the block structure 'blocks'
+   for( auto block = blocks->begin(); block != blocks->end(); ++block )
+   {
+
+      BoundaryHandling_T * boundaryHandling = block->getData< BoundaryHandling_T >( boundaryHandlingId );
+
+      CellInterval domainBB( domainBBInGlobalCellCoordinates );
+      blocks->transformGlobalToBlockLocalCellInterval( domainBB, *block );
+
+      CellInterval west( domainBB.xMin(), domainBB.yMin(), domainBB.zMin(), domainBB.xMin(), domainBB.yMax(), domainBB.zMax() );
+      CellInterval east( domainBB.xMax(), domainBB.yMin(), domainBB.zMin(), domainBB.xMax(), domainBB.yMax(), domainBB.zMax() );
+      CellInterval south( domainBB.xMin(), domainBB.yMin(), domainBB.zMin(), domainBB.xMax(), domainBB.yMin(), domainBB.zMax() );
+      CellInterval north( domainBB.xMin(), domainBB.yMax(), domainBB.zMin(), domainBB.xMax(), domainBB.yMax(), domainBB.zMax() );
+
+      // Set north boundary to large value
+      for( auto cell = north.begin(); cell != north.end(); ++cell )
+      {
+
+         const real_t val = real_t(2);
+         boundaryHandling->forceBoundary( Dirichlet_Flag, cell->x(), cell->y(), cell->z(), pde::Dirichlet< Stencil_T, flag_t >::DirichletBC( val ) );
+
+      }
+
+      // Set south boundary to large value
+      for( auto cell = south.begin(); cell != south.end(); ++cell )
+      {
+
+         const real_t val = real_t(-2);
+         boundaryHandling->forceBoundary( Dirichlet_Flag, cell->x(), cell->y(), cell->z(), pde::Dirichlet< Stencil_T, flag_t >::DirichletBC( val ) );
+
+      }
+
+      // Set all other boundaries to homogeneous Neumann BCs
+      for( auto & ci : { west, east} )
+      {
+         boundaryHandling->forceBoundary( Neumann_Flag, ci, pde::Neumann< Stencil_T, flag_t >::NeumannBC( real_t( 0 ) ) );
+      }
+
+      // All the remaining cells need to be marked as being fluid. The 'fillWithDomain' call does just that.
+      boundaryHandling->fillWithDomain( domainBB );
+   }
+}
+
+
+
+void copyWeightsToStencilField( const shared_ptr< StructuredBlockStorage > & blocks, const std::vector<real_t> & weights, const BlockDataID & stencilId )
+{
+   for( auto block = blocks->begin(); block != blocks->end(); ++block )
+   {
+      StencilField_T * stencil = block->getData< StencilField_T >( stencilId );
+      
+      WALBERLA_FOR_ALL_CELLS_XYZ(stencil,
+         for( auto dir = Stencil_T::begin(); dir != Stencil_T::end(); ++dir )
+            stencil->get(x,y,z,dir.toIdx()) = weights[ dir.toIdx() ];
+      );
+   }
+}
+
+
+
+
+int main( int argc, char** argv )
+{
+   debug::enterTestMode();
+
+   mpi::Environment env( argc, argv );
+
+   const uint_t processes = uint_c( MPIManager::instance()->numProcesses() );
+   if( processes != uint_t(1) && processes != uint_t(4) && processes != uint_t(8) )
+      WALBERLA_ABORT( "The number of processes must be equal to 1, 4, or 8!" );
+
+   logging::Logging::printHeaderOnStream();
+   WALBERLA_ROOT_SECTION() { logging::Logging::instance()->setLogLevel( logging::Logging::PROGRESS ); }
+
+   bool shortrun = false;
+   for( int i = 1; i < argc; ++i )
+      if( std::strcmp( argv[i], "--shortrun" ) == 0 ) shortrun = true;
+
+   const uint_t xBlocks = ( processes == uint_t(1) ) ? uint_t(1) : ( ( processes == uint_t(4) ) ? uint_t(2) : uint_t(4) );
+   const uint_t yBlocks = ( processes == uint_t(1) ) ? uint_t(1) : uint_t(2);
+   const uint_t xCells = ( processes == uint_t(1) ) ? uint_t(200) : ( ( processes == uint_t(4) ) ? uint_t(100) : uint_t(50) );
+   const uint_t yCells = ( processes == uint_t(1) ) ? uint_t(100) : uint_t(50);
+   const real_t xSize = real_t(2);
+   const real_t ySize = real_t(1);
+   const real_t dx = xSize / real_c( xBlocks * xCells + uint_t(1) );
+   const real_t dy = ySize / real_c( yBlocks * yCells + uint_t(1) );
+   auto blocks = blockforest::createUniformBlockGrid( math::AABB( real_t(0.5) * dx, real_t(0.5) * dy, real_t(0),
+                                                                  xSize - real_t(0.5) * dx, ySize - real_t(0.5) * dy, dx ),
+                                                      xBlocks, yBlocks, uint_t(1),
+                                                      xCells, yCells, uint_t(1),
+                                                      true,
+                                                      false, false, false );
+
+   BlockDataID solId = field::addToStorage< Field_T >( blocks, "sol", real_t(0), field::zyxf, uint_t(1) );
+   BlockDataID rId = field::addToStorage< Field_T >( blocks, "r", real_t(0), field::zyxf, uint_t(1) );
+   BlockDataID dId = field::addToStorage< Field_T >( blocks, "d", real_t(0), field::zyxf, uint_t(1) );
+   BlockDataID zId = field::addToStorage< Field_T >( blocks, "z", real_t(0), field::zyxf, uint_t(1) );
+
+   BlockDataID rhsId = field::addToStorage< Field_T >( blocks, "rhs", real_t(0), field::zyxf, uint_t(1) );
+
+   initRHS( blocks, rhsId );
+
+   blockforest::communication::UniformBufferedScheme< Stencil_T > synchronizeD( blocks );
+   synchronizeD.addPackInfo( make_shared< field::communication::PackInfo< Field_T > >( dId ) );
+
+   std::vector< real_t > weights( Stencil_T::Size );
+   weights[ Stencil_T::idx[ stencil::C ] ] = real_t(2) / ( blocks->dx() * blocks->dx() ) + real_t(2) / ( blocks->dy() * blocks->dy() ) + real_t(4) * math::PI * math::PI;
+   weights[ Stencil_T::idx[ stencil::N ] ] = real_t(-1) / ( blocks->dy() * blocks->dy() );
+   weights[ Stencil_T::idx[ stencil::S ] ] = real_t(-1) / ( blocks->dy() * blocks->dy() );
+   weights[ Stencil_T::idx[ stencil::E ] ] = real_t(-1) / ( blocks->dx() * blocks->dx() );
+   weights[ Stencil_T::idx[ stencil::W ] ] = real_t(-1) / ( blocks->dx() * blocks->dx() );
+
+
+   BlockDataID stencilId = field::addToStorage< StencilField_T >( blocks, "w" );
+   BlockDataID stencilBCadaptedId = field::addToStorage< StencilField_T >( blocks, "w" );
+
+   BlockDataID flagId = field::addFlagFieldToStorage< FlagField_T >( blocks, "flagField" );
+
+   copyWeightsToStencilField( blocks, weights, stencilId );
+   copyWeightsToStencilField( blocks, weights, stencilBCadaptedId );
+
+   BlockDataID boundaryHandlingId = blocks->addBlockData< BoundaryHandling_T >( MyBoundaryHandling( rhsId, stencilId, stencilBCadaptedId, flagId, *blocks ),
+                                                                                "boundary handling" );
+
+   // Test Dirichlet BCs //
+   setBoundaryConditionsDirichl( blocks, boundaryHandlingId );
+
+   
+   SweepTimeloop timeloop( blocks, uint_t(2) );
+
+   timeloop.add()
+            << Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingId ), "boundary handling" )
+            << AfterFunction(
+                     pde::CGIteration< Stencil_T >( blocks->getBlockStorage(), solId, rId, dId, zId, rhsId, stencilBCadaptedId,
+                                                    shortrun ? uint_t( 10 ) : uint_t( 10000 ), synchronizeD, real_c( 1e-6 ) ), "CG iteration" );
+
+   timeloop.singleStep();
+
+   // Check values for Dirichlet BCs //
+   if( !shortrun )
+   {
+      Cell cellNearBdry( 75, 2, 0 );
+      real_t solNearBdry( real_c(-0.16347) );
+      Cell cellNearBdryLrg( 24, 95, 0 );
+      real_t solNearBdryLrg( real_c(201.47) );
+      Cell cellDomCentr( 100, 50, 0 );
+      real_t solDomCentr( real_c(0.37587) );
+
+      for( auto block = blocks->begin(); block != blocks->end(); ++block )
+      {
+         Field_T * sol = block->getData < Field_T > ( solId );
+         if( blocks->getBlockCellBB( *block ).contains( cellNearBdry ) )
+         {
+            blocks->transformGlobalToBlockLocalCell(cellNearBdry,*block);
+            // WALBERLA_LOG_DEVEL( "Solution value at cell " << cellNearBdry << ": " << sol->get(cellNearBdry) );
+            WALBERLA_CHECK_LESS( std::fabs( solNearBdry - sol->get( cellNearBdry ) ) / solNearBdry, 0.0001, "Invalid value in cell " << cellNearBdry );
+         }
+         if( blocks->getBlockCellBB( *block ).contains( cellNearBdryLrg ) )
+         {
+            blocks->transformGlobalToBlockLocalCell(cellNearBdryLrg,*block);
+            // WALBERLA_LOG_DEVEL( "Solution value at cell " << cellNearBdryLrg << ": " << sol->get(cellNearBdryLrg) );
+            WALBERLA_CHECK_LESS( std::fabs( solNearBdryLrg - sol->get( cellNearBdryLrg ) ) / solNearBdryLrg, 0.0001, "Invalid value in cell " << cellNearBdryLrg );
+         }
+         if( blocks->getBlockCellBB( *block ).contains( cellDomCentr ) )
+         {
+            blocks->transformGlobalToBlockLocalCell(cellDomCentr,*block);
+            // WALBERLA_LOG_DEVEL( "Solution value at cell " << cellDomCentr << ": " << sol->get(cellDomCentr) );
+            WALBERLA_CHECK_LESS( std::fabs( solDomCentr - sol->get( cellDomCentr ) ) / solDomCentr, 0.0001, "Invalid value in cell " << cellDomCentr );
+         }
+      }
+   }
+   else
+   {
+      Cell cellNearBdry( 75, 2, 0 );
+      real_t solNearBdry( real_c(-0.008355) );
+      Cell cellNearBdryLrg( 24, 95, 0 );
+      real_t solNearBdryLrg( real_c(132.188) );
+      Cell cellDomCentr( 100, 50, 0 );
+      real_t solDomCentr( 0.017603f );
+
+      for( auto block = blocks->begin(); block != blocks->end(); ++block )
+      {
+         Field_T * sol = block->getData < Field_T > ( solId );
+         if( blocks->getBlockCellBB( *block ).contains( cellNearBdry ) )
+         {
+             blocks->transformGlobalToBlockLocalCell(cellNearBdry,*block);
+             // WALBERLA_LOG_DEVEL( "Solution value at cell " << cellNearBdry << ": " << sol->get(cellNearBdry) );
+             WALBERLA_CHECK_LESS( std::fabs( solNearBdry - sol->get( cellNearBdry ) ) / solNearBdry, 0.0001, "Invalid value in cell " << cellNearBdry );
+         }
+         if( blocks->getBlockCellBB( *block ).contains( cellNearBdryLrg ) )
+         {
+             blocks->transformGlobalToBlockLocalCell(cellNearBdryLrg,*block);
+             // WALBERLA_LOG_DEVEL( "Solution value at cell " << cellNearBdryLrg << ": " << sol->get(cellNearBdryLrg) );
+             WALBERLA_CHECK_LESS( std::fabs( solNearBdryLrg - sol->get( cellNearBdryLrg ) ) / solNearBdryLrg, 0.0001, "Invalid value in cell " << cellNearBdryLrg );
+         }
+         if( blocks->getBlockCellBB( *block ).contains( cellDomCentr ) )
+         {
+             blocks->transformGlobalToBlockLocalCell(cellDomCentr,*block);
+             // WALBERLA_LOG_DEVEL( "Solution value at cell " << cellDomCentr << ": " << sol->get(cellDomCentr) );
+             WALBERLA_CHECK_LESS( std::fabs( solDomCentr - sol->get( cellDomCentr ) ) / solDomCentr, 0.0001, "Invalid value in cell " << cellDomCentr );
+         }
+      }
+
+   }
+
+   if( !shortrun )
+   {
+      vtk::writeDomainDecomposition( blocks );
+      field::createVTKOutput< Field_T >( solId, *blocks, "solution_Dirich" )();
+   }
+
+
+   // Test Mixed BCs and re-setting BCs //
+   setBoundaryConditionsMixed( blocks, boundaryHandlingId );
+   // set RHS to zero to get 'ramp' as solution
+   setRHSConstValue( blocks, rhsId, real_t(0));
+
+   timeloop.singleStep();
+
+   // Check values for mixed  BCs //
+   // TODO
+
+   if( !shortrun )
+   {
+      field::createVTKOutput< Field_T >( solId, *blocks, "solution_Mixed" )();
+//      field::createVTKOutput< Field_T >( rhsId, *blocks, "rhs_Mixed" )();
+//      field::createVTKOutput< StencilField_T >( stencilId, *blocks, "originalStencils_Mixed" )();
+//      field::createVTKOutput< StencilField_T >( stencilBCadaptedId, *blocks, "adaptedStencils_Mixed" )();
+   }
+
+   logging::Logging::printFooterOnStream();
+   return EXIT_SUCCESS;
+}
+}
+
+int main( int argc, char* argv[] )
+{
+  return walberla::main( argc, argv );
+}
diff --git a/tests/pde/CMakeLists.txt b/tests/pde/CMakeLists.txt
index 58845042e3ce5173ec1da1baa1a03f03c0279464..5dde027cce3295929ba57314b00a8494a11000a3 100644
--- a/tests/pde/CMakeLists.txt
+++ b/tests/pde/CMakeLists.txt
@@ -22,4 +22,8 @@ waLBerla_execute_test( NAME MGShortTest COMMAND $<TARGET_FILE:MGTest> --shortrun
 waLBerla_execute_test( NAME MGTest COMMAND $<TARGET_FILE:MGTest> PROCESSES 8 CONFIGURATIONS Release RelWithDbgInfo )
 
 waLBerla_compile_test( FILES MGConvergenceTest.cpp DEPENDS blockforest timeloop vtk )
-waLBerla_execute_test( NAME MGConvergenceTest COMMAND $<TARGET_FILE:MGConvergenceTest> PROCESSES 8 CONFIGURATIONS Release RelWithDbgInfo )
\ No newline at end of file
+waLBerla_execute_test( NAME MGConvergenceTest COMMAND $<TARGET_FILE:MGConvergenceTest> PROCESSES 8 CONFIGURATIONS Release RelWithDbgInfo )
+
+waLBerla_compile_test( FILES BoundaryTest.cpp DEPENDS blockforest timeloop vtk boundary )
+waLBerla_execute_test( NAME BoundaryShortTest COMMAND $<TARGET_FILE:BoundaryTest> --shortrun PROCESSES 8 )
+waLBerla_execute_test( NAME BoundaryTest COMMAND $<TARGET_FILE:BoundaryTest> PROCESSES 8 CONFIGURATIONS Release RelWithDbgInfo )
diff --git a/tests/pe_coupling/discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp b/tests/pe_coupling/discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp
index 80d152379546eb88a70eca1ed23cead805b32466..1a4293fec9670eabbe9b0669f24e37b61d674ee3 100644
--- a/tests/pe_coupling/discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp
+++ b/tests/pe_coupling/discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp
@@ -62,7 +62,6 @@ namespace hindered_settling_dynamics_dpm
 using namespace walberla;
 using walberla::uint_t;
 
-
 ///////////////
 // CONSTANTS //
 ///////////////
diff --git a/tests/pe_coupling/discrete_particle_methods/SphereWallCollisionBehaviorDPM.cpp b/tests/pe_coupling/discrete_particle_methods/SphereWallCollisionBehaviorDPM.cpp
index 2cbdd76dc50ff5df1d5ae15dfe98a02fcbeb4c6a..f40e6156de97586ebb7029e269caa4862a62e5b3 100644
--- a/tests/pe_coupling/discrete_particle_methods/SphereWallCollisionBehaviorDPM.cpp
+++ b/tests/pe_coupling/discrete_particle_methods/SphereWallCollisionBehaviorDPM.cpp
@@ -61,7 +61,6 @@ namespace sphere_wall_collision_behavior_dpm
 using namespace walberla;
 using walberla::uint_t;
 
-
 ///////////////
 // CONSTANTS //
 ///////////////
diff --git a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
index caa0987ed52d9bb2e03b6ffd6618dbe7063e81d1..2e449157fc7bf338a321fcb794e917d63e227e84 100644
--- a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
@@ -73,11 +73,10 @@ namespace drag_force_sphere_mem
 
 using namespace walberla;
 using walberla::uint_t;
-using lbm::force_model::SimpleConstant;
 
 // PDF field, flag field & body field
 using ForceModel_T = lbm::force_model::LuoConstant;
-typedef lbm::D3Q19< lbm::collision_model::TRT, false, ForceModel_T, 1>  LatticeModel_T;
+typedef lbm::D3Q19< lbm::collision_model::TRT, false, ForceModel_T >  LatticeModel_T;
 
 using Stencil_T = LatticeModel_T::Stencil;
 using PdfField_T = lbm::PdfField<LatticeModel_T>;
diff --git a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp
index f940ae03303270c35832b1cd2c4ef16ddc95515f..358d165bae76cd446dc8709fd7059b2eb9f6a0ae 100644
--- a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp
@@ -83,10 +83,9 @@ namespace drag_force_sphere_mem_refinement
 
 using namespace walberla;
 using walberla::uint_t;
-using lbm::force_model::SimpleConstant;
 
 // PDF field, flag field & body field
-typedef lbm::D3Q19< lbm::collision_model::TRT, false, lbm::force_model::SimpleConstant, 1>  LatticeModel_T;
+typedef lbm::D3Q19< lbm::collision_model::TRT, false, lbm::force_model::SimpleConstant >  LatticeModel_T;
 
 using Stencil_T = LatticeModel_T::Stencil;
 using PdfField_T = lbm::PdfField<LatticeModel_T>;
@@ -531,7 +530,8 @@ int main( int argc, char **argv )
    ////////////////////////
 
    // create the lattice model
-   LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( omega ), SimpleConstant( Vector3<real_t> ( setup.extForce, 0, 0 ) ) );
+   LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( omega ),
+                                                 lbm::force_model::SimpleConstant( Vector3<real_t> ( setup.extForce, 0, 0 ) ) );
 
    // add PDF field ( uInit = <0,0,0>, rhoInit = 1 )
    BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (zyxf)", latticeModel,
diff --git a/tests/pe_coupling/momentum_exchange_method/GlobalBodyAsBoundaryMEMStaticRefinement.cpp b/tests/pe_coupling/momentum_exchange_method/GlobalBodyAsBoundaryMEMStaticRefinement.cpp
index 2afeebd2d7d60a9653250e6f8bd38c2d4ae57382..b494339bbeaa0c1c8c7af8e0e3185a741d5b70ea 100644
--- a/tests/pe_coupling/momentum_exchange_method/GlobalBodyAsBoundaryMEMStaticRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/GlobalBodyAsBoundaryMEMStaticRefinement.cpp
@@ -79,7 +79,7 @@ using walberla::uint_t;
 //////////////
 
 // PDF field, flag field & body field
-typedef lbm::D3Q19< lbm::collision_model::TRT, false >  LatticeModel_T;
+typedef lbm::D3Q19< lbm::collision_model::TRT >  LatticeModel_T;
 using Stencil_T = LatticeModel_T::Stencil;
 using PdfField_T = lbm::PdfField<LatticeModel_T>;
 
diff --git a/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp b/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
index 0f4dbd3fecab43b90f7a5c12601bdb9beb3f16e0..2ebe21f515c607ccc96698d66d137dc0ce3f2b4e 100644
--- a/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
@@ -77,8 +77,7 @@ using walberla::uint_t;
 //////////////
 
 // pdf field & flag field
-
-typedef lbm::D3Q19< lbm::collision_model::TRT, false >  LatticeModel_T;
+typedef lbm::D3Q19< lbm::collision_model::TRT >  LatticeModel_T;
 using Stencil_T = LatticeModel_T::Stencil;
 using PdfField_T = lbm::PdfField<LatticeModel_T>;
 
diff --git a/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp b/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
index 2ecfd788de6d83a3dc1e9677f7bdba9a437095b9..678aa5822d5525967d045d2065f0e54aff873d2f 100644
--- a/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
@@ -78,14 +78,12 @@ namespace periodic_particle_channel_mem
 using namespace walberla;
 using walberla::uint_t;
 
-
-
 //////////////
 // TYPEDEFS //
 //////////////
 
 // pdf field & flag field
-typedef lbm::D3Q19< lbm::collision_model::TRT, false >  LatticeModel_T;
+typedef lbm::D3Q19< lbm::collision_model::TRT >  LatticeModel_T;
 using Stencil_T = LatticeModel_T::Stencil;
 using PdfField_T = lbm::PdfField<LatticeModel_T>;
 
diff --git a/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp b/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
index a238ce167bc2e2b17f7ad38d5faf1d4667b8bbdc..ca31d93920e212d15884fc0c77a152f55805fbfb 100644
--- a/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
@@ -76,14 +76,13 @@ namespace segre_silberberg_mem
 
 using namespace walberla;
 using walberla::uint_t;
-using lbm::force_model::SimpleConstant;
 
 //////////////
 // TYPEDEFS //
 //////////////
 
 // PDF field, flag field & body field
-typedef lbm::D3Q19< lbm::collision_model::TRT, false, SimpleConstant >  LatticeModel_T;
+typedef lbm::D3Q19< lbm::collision_model::TRT, false, lbm::force_model::SimpleConstant >  LatticeModel_T;
 using Stencil_T = LatticeModel_T::Stencil;
 using PdfField_T = lbm::PdfField<LatticeModel_T>;
 
@@ -582,7 +581,8 @@ int main( int argc, char **argv )
    ////////////////////////
 
    // create the lattice model
-   LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( omega ), SimpleConstant( Vector3<real_t> ( setup.forcing, real_t(0), real_t(0) ) ) );
+   LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( omega ),
+                                                 lbm::force_model::SimpleConstant( Vector3<real_t> ( setup.forcing, real_t(0), real_t(0) ) ) );
 
    // add PDF field
    BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (zyxf)", latticeModel,
diff --git a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp
index b90dcf0056803b867190b27ccae6dd2e640620ec..482bc1883b9335798770053d4d4e6889d1849779 100644
--- a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp
@@ -80,7 +80,7 @@ using walberla::uint_t;
 //////////////
 
 // PDF field, flag field & body field
-typedef lbm::D3Q19< lbm::collision_model::TRT, false >  LatticeModel_T;
+typedef lbm::D3Q19< lbm::collision_model::TRT >  LatticeModel_T;
 using Stencil_T = LatticeModel_T::Stencil;
 using PdfField_T = lbm::PdfField<LatticeModel_T>;
 
diff --git a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp
index 885e9ed5f5e9e075800de8ccaabcd9e7a4d96203..558747d90711f0a7e5f01c600e3621c3dadb849e 100644
--- a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp
@@ -86,7 +86,7 @@ using walberla::uint_t;
 //////////////
 
 // PDF field, flag field & body field
-typedef lbm::D3Q19< lbm::collision_model::TRT, false >  LatticeModel_T;
+typedef lbm::D3Q19< lbm::collision_model::TRT >  LatticeModel_T;
 using Stencil_T = LatticeModel_T::Stencil;
 using PdfField_T = lbm::PdfField<LatticeModel_T>;
 
diff --git a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp
index 3cf8bf5fd8b8ee09e99f7455045ed437a287a376..f0863496bfc2b1a6c311447e1661de8e37cd1cd2 100644
--- a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp
@@ -82,7 +82,7 @@ using walberla::uint_t;
 //////////////
 
 // PDF field, flag field & body field
-typedef lbm::D3Q19< lbm::collision_model::TRT, false >  LatticeModel_T;
+typedef lbm::D3Q19< lbm::collision_model::TRT >  LatticeModel_T;
 using Stencil_T = LatticeModel_T::Stencil;
 using PdfField_T = lbm::PdfField<LatticeModel_T>;
 
diff --git a/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp b/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
index adca4508693a8a5802c3ea3a5fa883ad00bef304..1d3bab503558cede6f0d84773ec7b717869157c2 100644
--- a/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
@@ -80,8 +80,7 @@ using namespace walberla;
 using walberla::uint_t;
 
 // PDF field, flag field & body field
-using ForceModel_T = lbm::force_model::None;
-typedef lbm::D3Q19<lbm::collision_model::TRT, false, ForceModel_T> LatticeModel_T;
+typedef lbm::D3Q19<lbm::collision_model::TRT> LatticeModel_T;
 
 using Stencil_T = LatticeModel_T::Stencil;
 using PdfField_T = lbm::PdfField<LatticeModel_T>;
diff --git a/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp b/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp
index 0205b4695a0eeff177c2d6a00f5acac3ac8c6315..b654d5dfaacbad39af9e10ebd3e5d803b9c44d33 100644
--- a/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp
@@ -69,14 +69,12 @@ namespace taylor_coette_flow_mem
 using namespace walberla;
 using walberla::uint_t;
 
-
-
 //////////////
 // TYPEDEFS //
 //////////////
 
 // pdf field & flag field
-typedef lbm::D3Q19< lbm::collision_model::TRT, false >  LatticeModel_T;
+typedef lbm::D3Q19< lbm::collision_model::TRT >  LatticeModel_T;
 using Stencil_T = LatticeModel_T::Stencil;
 using PdfField_T = lbm::PdfField<LatticeModel_T>;
 
diff --git a/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp b/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
index 7f43c9fd26fd49fc1453fcb677542cab9064c5b4..68a74cde6dda8d148e7928429a5ee9b4f8dfe168 100644
--- a/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
@@ -76,7 +76,7 @@ using namespace walberla;
 using walberla::uint_t;
 
 // PDF field, flag field & body field
-typedef lbm::D3Q19< lbm::collision_model::TRT, false, lbm::force_model::None, 1>  LatticeModel_T;
+typedef lbm::D3Q19< lbm::collision_model::TRT >  LatticeModel_T;
 
 using Stencil_T = LatticeModel_T::Stencil;
 using PdfField_T = lbm::PdfField<LatticeModel_T>;
diff --git a/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSM.cpp b/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSM.cpp
index fb92b0c68975877bf00a5066acd6ce08735113ba..b5e53a23fc1856a18548a4929634e94c5b320179 100644
--- a/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSM.cpp
+++ b/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSM.cpp
@@ -73,11 +73,10 @@ namespace drag_force_sphere_psm
 
 using namespace walberla;
 using walberla::uint_t;
-using lbm::force_model::SimpleConstant;
 
 // PDF field, flag field & body field
 using ForceModel_T = lbm::force_model::LuoConstant;
-typedef lbm::D3Q19< lbm::collision_model::SRT, false, ForceModel_T, 1>  LatticeModel_T;
+typedef lbm::D3Q19< lbm::collision_model::SRT, false, ForceModel_T>  LatticeModel_T;
 
 using Stencil_T = LatticeModel_T::Stencil;
 using PdfField_T = lbm::PdfField<LatticeModel_T>;
diff --git a/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSMRefinement.cpp b/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSMRefinement.cpp
index 58add3e2404c162f19f653552033e1be81a322ea..b0ca02da90f9b912399af654658655dd8f8219b1 100644
--- a/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSMRefinement.cpp
+++ b/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSMRefinement.cpp
@@ -83,10 +83,9 @@ namespace drag_force_sphere_psm_refinement
 
 using namespace walberla;
 using walberla::uint_t;
-using lbm::force_model::SimpleConstant;
 
 // PDF field, flag field & body field
-typedef lbm::D3Q19< lbm::collision_model::SRT, false, lbm::force_model::SimpleConstant, 1>  LatticeModel_T;
+typedef lbm::D3Q19< lbm::collision_model::SRT, false, lbm::force_model::SimpleConstant >  LatticeModel_T;
 
 using Stencil_T = LatticeModel_T::Stencil;
 using PdfField_T = lbm::PdfField<LatticeModel_T>;
@@ -536,7 +535,7 @@ int main( int argc, char **argv )
    ////////////////////////
 
    // create the lattice model
-   LatticeModel_T latticeModel = LatticeModel_T( omega, SimpleConstant( Vector3<real_t> ( setup.extForce, 0, 0 ) ) );
+   LatticeModel_T latticeModel = LatticeModel_T( omega, lbm::force_model::SimpleConstant( Vector3<real_t> ( setup.extForce, 0, 0 ) ) );
 
    // add PDF field ( uInit = <0,0,0>, rhoInit = 1 )
    BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (zyxf)", latticeModel,
diff --git a/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp b/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp
index 46f2a10e4fbc83242bd90a2b1510a4425543a090..935a0c1d51ef9dfbc175772cb1a7540e617dcaf1 100644
--- a/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp
+++ b/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp
@@ -76,14 +76,13 @@ namespace segre_silberberg_psm
 
 using namespace walberla;
 using walberla::uint_t;
-using lbm::force_model::SimpleConstant;
 
 //////////////
 // TYPEDEFS //
 //////////////
 
 // PDF field, flag field & body field
-typedef lbm::D3Q19< lbm::collision_model::SRT, false, SimpleConstant >  LatticeModel_T;
+typedef lbm::D3Q19< lbm::collision_model::SRT, false, lbm::force_model::SimpleConstant >  LatticeModel_T;
 using Stencil_T = LatticeModel_T::Stencil;
 using PdfField_T = lbm::PdfField<LatticeModel_T>;
 
@@ -532,7 +531,7 @@ int main( int argc, char **argv )
    ////////////////////////
 
    // create the lattice model
-   LatticeModel_T latticeModel = LatticeModel_T( omega, SimpleConstant( Vector3<real_t> ( setup.forcing, real_c(0), real_c(0) ) ) );
+   LatticeModel_T latticeModel = LatticeModel_T( omega, lbm::force_model::SimpleConstant( Vector3<real_t> ( setup.forcing, real_c(0), real_c(0) ) ) );
 
    // add PDF field
    BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (zyxf)", latticeModel,
diff --git a/tests/pe_coupling/partially_saturated_cells_method/TorqueSpherePSM.cpp b/tests/pe_coupling/partially_saturated_cells_method/TorqueSpherePSM.cpp
index 64114c43e0f76a36da85962123e6ca71074fb144..7cd33d47ed694d99b08525e161837bb85b4b8b7a 100644
--- a/tests/pe_coupling/partially_saturated_cells_method/TorqueSpherePSM.cpp
+++ b/tests/pe_coupling/partially_saturated_cells_method/TorqueSpherePSM.cpp
@@ -75,7 +75,7 @@ using namespace walberla;
 using walberla::uint_t;
 
 // PDF field, flag field & body field
-typedef lbm::D3Q19< lbm::collision_model::SRT, false, lbm::force_model::None, 1>  LatticeModel_T;
+typedef lbm::D3Q19< lbm::collision_model::SRT, false >  LatticeModel_T;
 
 using Stencil_T = LatticeModel_T::Stencil;
 using PdfField_T = lbm::PdfField<LatticeModel_T>;