src/lbm_generated/gpu/NonuniformGeneratedGPUPdfPackInfo.impl.h

@@ -55,7 +55,7 @@ std::shared_ptr< NonuniformGeneratedGPUPdfPackInfo< PdfField_T > >
    auto handling = std::make_shared<NonuniformGPUCommDataHandling< LatticeStorageSpecification_T > >(blocks);
    BlockDataID commDataID = sbf->addBlockData(handling, dataIdentifier);
 
-   return std::make_shared<NonuniformGeneratedGPUPdfPackInfo< PdfField_T > >(pdfFieldID, commDataID);
+   return std::make_shared<NonuniformGeneratedGPUPdfPackInfo< PdfField_T > >(sbf->getNumberOfLevels(), pdfFieldID, commDataID);
 }
 
 
@@ -92,6 +92,56 @@ void NonuniformGeneratedGPUPdfPackInfo< PdfField_T >::communicateLocalEqualLevel
    kernels_.localCopyDirection(srcField, srcRegion, dstField, dstRegion, dir, stream);
 }
 
+template< typename PdfField_T>
+void NonuniformGeneratedGPUPdfPackInfo< PdfField_T >::addForLocalEqualLevelComm(
+   const Block* sender, Block* receiver, stencil::Direction dir)
+{
+   const uint_t level = sender->getLevel();
+   auto srcField = const_cast< Block* >(sender)->getData< PdfField_T >(pdfFieldID_);
+   auto dstField = receiver->getData< PdfField_T >(pdfFieldID_);
+
+   CellInterval srcRegion;
+   CellInterval dstRegion;
+   cell_idx_t gls = skipsThroughCoarseBlock(sender, dir) ? 2 : 1;
+   srcField->getSliceBeforeGhostLayer(dir, srcRegion, gls, false);
+   dstField->getGhostRegion(stencil::inverseDir[dir], dstRegion, gls, false);
+
+   strides[0] = int64_t(srcField->xStride());
+   strides[1] = int64_t(srcField->yStride());
+   strides[2] = int64_t(srcField->zStride());
+   strides[3] = int64_t(1 * int64_t(srcField->fStride()));
+
+   value_type* data_pdfs_dst = dstField->dataAt(dstRegion.xMin(), dstRegion.yMin(), dstRegion.zMin(), 0);
+   value_type* data_pdfs_src = srcField->dataAt(srcRegion.xMin(), srcRegion.yMin(), srcRegion.zMin(), 0);
+
+   const uint_t index = level * Stencil::Q + dir;
+   Vector3<int64_t> size(int64_c(srcRegion.xSize()), int64_c(srcRegion.ySize()), int64_c(srcRegion.zSize()));
+
+   equalCommDST[index][size].emplace_back(data_pdfs_dst);
+   equalCommSRC[index][size].emplace_back(data_pdfs_src);
+}
+
+
+template< typename PdfField_T>
+void NonuniformGeneratedGPUPdfPackInfo< PdfField_T >::communicateLocalEqualLevel(uint64_t level, uint8_t timestep, gpuStream_t stream)
+{
+   if(dirty){
+      sync();
+      dirty = false;
+   }
+   for (auto dir = CommunicationStencil::beginNoCenter(); dir != CommunicationStencil::end(); ++dir){
+      const uint_t index = level * Stencil::Q + *dir;
+      for (auto const& x : equalCommSRC[index]){
+         auto key = x.first;
+         value_type** data_pdfs_src_dp = equalCommSRCGPU[index][key];
+         value_type** data_pdfs_dst_dp = equalCommDSTGPU[index][key];
+         std::array< int64_t, 4 > size = { int64_c(equalCommSRC[index][key].size()), key[0], key[1], key[2] };
+
+         kernels_.blockLocalCopyDirection(data_pdfs_src_dp, data_pdfs_dst_dp, *dir, timestep, stream, size, strides);
+      }
+   }
+}
+
 
 template< typename PdfField_T>
 void NonuniformGeneratedGPUPdfPackInfo< PdfField_T >::packDataEqualLevelImpl(
@@ -168,21 +218,13 @@ void NonuniformGeneratedGPUPdfPackInfo< PdfField_T >::communicateLocalCoarseToFi
       Direction const unpackDir      = dstIntervals[index].first;
       CellInterval dstInterval = dstIntervals[index].second;
 
-      uint_t packSize      = kernels_.size(srcInterval);
-
 #ifndef NDEBUG
       Direction const packDir        = srcIntervals[index].first;
       WALBERLA_ASSERT_EQUAL(packDir, stencil::inverseDir[unpackDir])
       uint_t unpackSize = kernels_.redistributeSize(dstInterval);
-      WALBERLA_ASSERT_EQUAL(packSize, unpackSize)
+      WALBERLA_ASSERT_EQUAL(kernels_.size(srcInterval), unpackSize)
 #endif
-
-      // TODO: This is a dirty workaround. Code-generate direct redistribution!
-      unsigned char *buffer;
-      WALBERLA_GPU_CHECK( gpuMalloc( &buffer, packSize))
-      kernels_.packAll(srcField, srcInterval, buffer, stream);
-      kernels_.unpackRedistribute(dstField, dstInterval, buffer, unpackDir, stream);
-      WALBERLA_GPU_CHECK(gpuFree(buffer))
+      kernels_.localCopyRedistribute(srcField, srcInterval, dstField, dstInterval, unpackDir, stream);
    }
 }
 
@@ -190,6 +232,9 @@ template< typename PdfField_T>
 void NonuniformGeneratedGPUPdfPackInfo< PdfField_T >::communicateLocalCoarseToFine(
    const Block* coarseSender, Block* fineReceiver, stencil::Direction dir, GpuBuffer_T & buffer, gpuStream_t stream)
 {
+   // WARNING: This function uses an inplace buffer array.
+   // If possible the direct communicateLocalCoarseToFine without buffer array should be used
+
    auto srcField = const_cast< Block* >(coarseSender)->getData< PdfField_T >(pdfFieldID_);
    auto dstField = fineReceiver->getData< PdfField_T >(pdfFieldID_);
 
@@ -269,22 +314,16 @@ void NonuniformGeneratedGPUPdfPackInfo< PdfField_T >::communicateLocalFineToCoar
 
    CellInterval srcInterval;
    srcField->getGhostRegion(dir, srcInterval, 2);
-   uint_t packSize = kernels_.partialCoalescenceSize(srcInterval, dir);
 
    CellInterval dstInterval = getCoarseBlockCoalescenceInterval(coarseReceiver, fineSender->getId(),
                                                                 invDir, dstField);
 
 #ifndef NDEBUG
    uint_t unpackSize = kernels_.size(dstInterval, invDir);
-   WALBERLA_ASSERT_EQUAL(packSize, unpackSize)
+   WALBERLA_ASSERT_EQUAL(kernels_.partialCoalescenceSize(srcInterval, dir), unpackSize)
 #endif
 
-   // TODO: This is a dirty workaround. Code-generate direct redistribution!
-   unsigned char *buffer;
-   WALBERLA_GPU_CHECK( gpuMalloc( &buffer, packSize))
-   kernels_.packPartialCoalescence(srcField, maskField, srcInterval, buffer, dir, stream);
-   kernels_.unpackCoalescence(dstField, dstInterval, buffer, invDir, stream);
-   WALBERLA_GPU_CHECK(gpuFree(buffer))
+   kernels_.localPartialCoalescence(srcField, maskField, srcInterval, dstField, dstInterval, dir, stream);
 }
 
 

src/lbm_generated/refinement/BasicRecursiveTimeStep.h

@@ -31,6 +31,7 @@
 namespace walberla {
 
 using blockforest::communication::NonUniformBufferedScheme;
+using BlockFunction = std::function<void (const uint_t)>; // parameters: block, level
 
 namespace lbm_generated {
 
@@ -73,11 +74,14 @@ class BasicRecursiveTimeStep
    void operator() () { timestep(0); };
    void addRefinementToTimeLoop(SweepTimeloop & timeloop, uint_t level=0);
 
+   inline void addPostBoundaryHandlingBlockFunction( const BlockFunction & function );
+
  private:
    void timestep(uint_t level);
    void ghostLayerPropagation(Block * block);
-   std::function<void()> executeStreamCollideOnLevel(uint_t level, bool withGhostLayerPropagation=false);
-   std::function<void()> executeBoundaryHandlingOnLevel(uint_t level);
+   std::function< void() > executeStreamCollideOnLevel(uint_t level, bool withGhostLayerPropagation=false);
+   std::function< void() > executeBoundaryHandlingOnLevel(uint_t level);
+   std::function< void() > executePostBoundaryBlockFunctions(uint_t level);
 
    std::shared_ptr< StructuredBlockForest > sbfs_;
    uint_t maxLevel_;
@@ -89,6 +93,8 @@ class BasicRecursiveTimeStep
 
    SweepCollection_T & sweepCollection_;
    BoundaryCollection_T & boundaryCollection_;
+
+   std::vector< BlockFunction >  globalPostBoundaryHandlingBlockFunctions_;
 };
 
 } // namespace lbm_generated

src/lbm_generated/refinement/BasicRecursiveTimeStep.impl.h

@@ -111,6 +111,7 @@ void BasicRecursiveTimeStep< PdfField_T, SweepCollection_T, BoundaryCollection_T
 
    // 1.5 Boundary Handling and Coalescence Preparation
    timeloop.addFuncBeforeTimeStep(executeBoundaryHandlingOnLevel(level), "Refinement Cycle: boundary handling on level " + std::to_string(level));
+   timeloop.addFuncBeforeTimeStep(executePostBoundaryBlockFunctions(level), "Refinement Cycle: post boundary handling block functions on level " + std::to_string(level));
 
    // 1.6 Fine to Coarse Communication, receiving end
    if(level < maxLevel_){
@@ -134,10 +135,12 @@ void BasicRecursiveTimeStep< PdfField_T, SweepCollection_T, BoundaryCollection_T
 
    // 2.5 Boundary Handling and Coalescence Preparation
    timeloop.addFuncBeforeTimeStep(executeBoundaryHandlingOnLevel(level), "Refinement Cycle: boundary handling on level " + std::to_string(level));
+   timeloop.addFuncBeforeTimeStep(executePostBoundaryBlockFunctions(level), "Refinement Cycle: post boundary handling block functions on level " + std::to_string(level));
 
    // 2.6 Fine to Coarse Communication, receiving end
-   if(level < maxLevel_)
+   if(level < maxLevel_){
       timeloop.addFuncBeforeTimeStep(commScheme_->communicateFineToCoarseFunctor(level + 1), "Refinement Cycle: communicate fine to coarse on level " + std::to_string(level + 1));
+   }
 
 }
 
@@ -176,6 +179,16 @@ std::function<void()>  BasicRecursiveTimeStep< PdfField_T, SweepCollection_T, Bo
    };
 }
 
+template< typename PdfField_T, typename SweepCollection_T, typename BoundaryCollection_T >
+std::function<void()>  BasicRecursiveTimeStep< PdfField_T, SweepCollection_T, BoundaryCollection_T >::executePostBoundaryBlockFunctions(uint_t level)
+{
+   return [level, this]() {
+      for( const auto& func : globalPostBoundaryHandlingBlockFunctions_ ){
+         func(level);
+      }
+   };
+}
+
 
 template< typename PdfField_T, typename SweepCollection_T, typename BoundaryCollection_T >
 void BasicRecursiveTimeStep< PdfField_T, SweepCollection_T, BoundaryCollection_T >::ghostLayerPropagation(Block * block)
@@ -193,6 +206,12 @@ void BasicRecursiveTimeStep< PdfField_T, SweepCollection_T, BoundaryCollection_T
    }
 }
 
+template< typename PdfField_T, typename SweepCollection_T, typename BoundaryCollection_T >
+inline void BasicRecursiveTimeStep< PdfField_T, SweepCollection_T, BoundaryCollection_T >::addPostBoundaryHandlingBlockFunction( const BlockFunction & function )
+{
+   globalPostBoundaryHandlingBlockFunctions_.emplace_back( function );
+}
+
 // Refinement Timestep from post collision state:
 //template< typename PdfField_T, typename LbSweep_T >
 //void BasicRecursiveTimeStep< PdfField_T, LbSweep_T >::timestep(uint_t level)

src/lbm_generated/refinement/RefinementScaling.h

@@ -20,44 +20,18 @@
 
 #pragma once
 
-#include "blockforest/BlockDataHandling.h"
+#include "core/DataTypes.h"
 
-#include "domain_decomposition/IBlock.h"
-#include "domain_decomposition/StructuredBlockStorage.h"
-
-namespace walberla
-{
-namespace lbm_generated
+namespace walberla::lbm_generated
 {
 
-class DefaultRefinementScaling : public blockforest::AlwaysInitializeBlockDataHandling< real_t >
+inline real_t relaxationRateScaling( real_t relaxationRate, uint_t refinementLevel )
 {
- public:
-   DefaultRefinementScaling(const weak_ptr< StructuredBlockStorage >& blocks, const real_t parameter)
-      : blocks_(blocks), parameter_(parameter){};
-
-   real_t* initialize(IBlock* const block) override
-   {
-      WALBERLA_ASSERT_NOT_NULLPTR(block)
-      auto blocks = blocks_.lock();
-      WALBERLA_CHECK_NOT_NULLPTR(blocks)
-
-      level_ = block->getBlockStorage().getLevel(*block);
-
-      const real_t level_scale_factor = real_c(uint_t(1) << level_);
-      const real_t one                = real_c(1.0);
-      const real_t half               = real_c(0.5);
-
-      return new real_t(parameter_ / (level_scale_factor * (-parameter_ * half + one) + parameter_ * half));
-   }
-   bool operator==(const DefaultRefinementScaling& other) const { return level_ == other.level_; }
-
- private:
-   const weak_ptr< StructuredBlockStorage > blocks_;
-   const real_t parameter_;
+   const real_t levelScaleFactor = real_c(uint_c(1) << refinementLevel);
+   const real_t one                = real_c(1.0);
+   const real_t half               = real_c(0.5);
 
-   uint_t level_;
-};
+   return real_c(relaxationRate / (levelScaleFactor * (-relaxationRate * half + one) + relaxationRate * half));
+}
 
-} // namespace lbm_generated
-} // namespace walberla
+} // namespace walberla::lbm_generated
\ No newline at end of file
\ No newline at end of file

src/mesh/CMakeLists.txt

@@ -12,6 +12,7 @@ endif()
 target_sources( mesh
       PRIVATE
       MeshConversion.h
+      Utility.h
       )
 
 add_subdirectory( blockforest )

src/mesh/Utility.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file Utility.h
+//! \ingroup mesh
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#pragma once
+
+#include "core/debug/CheckFunctions.h"
+#include "core/debug/Debug.h"
+
+
+namespace walberla::mesh
+{
+
+/**
+* \brief Color the faces of a mesh according to its vertices
+*
+* Iterate over all faces and colors them in their vertex color.
+* If no uniform coloring of the vertices is given, a default color is taken.
+*
+* \tparam MeshType The type of the Mesh
+*
+* \param mesh The Mesh source mesh
+* \param defaultColor Default color if no uniform coloring is given
+ */
+template< typename MeshType >
+void vertexToFaceColor(MeshType& mesh, const typename MeshType::Color& defaultColor)
+{
+   WALBERLA_CHECK(mesh.has_vertex_colors())
+   mesh.request_face_colors();
+
+   for (auto faceIt = mesh.faces_begin(); faceIt != mesh.faces_end(); ++faceIt)
+   {
+      typename MeshType::Color vertexColor;
+
+      bool useVertexColor = true;
+
+      auto vertexIt = mesh.fv_iter(*faceIt);
+      WALBERLA_ASSERT(vertexIt.is_valid())
+
+      vertexColor = mesh.color(*vertexIt);
+
+      ++vertexIt;
+      while (vertexIt.is_valid() && useVertexColor)
+      {
+         if (vertexColor != mesh.color(*vertexIt)) useVertexColor = false;
+         ++vertexIt;
+      }
+
+      mesh.set_color(*faceIt, useVertexColor ? vertexColor : defaultColor);
+   }
+}
+}
\ No newline at end of file

src/vtk/VTKOutput.cpp

@@ -153,7 +153,7 @@ void VTKOutput::init( const std::string & identifier )
    std::string uniqueIdentifier = uniqueIdentifierStrStream.str();
 
    if( VTKUID::exists( uniqueIdentifier ) )
-      WALBERLA_ABORT( "VTKOutput with identifier \"" << uniqueIdentifier << "\" already exists. VTKOutput objects must have unique identifiers!" );
+      WALBERLA_ABORT( "VTKOutput with identifier \"" << uniqueIdentifier << "\" already exists. VTKOutput objects must have unique identifiers!" )
 
    VTKUID( uniqueIdentifier, true );
 
@@ -179,7 +179,7 @@ void VTKOutput::init( const std::string & identifier )
          filesystem::create_directories( path );
    }
 
-   WALBERLA_MPI_WORLD_BARRIER();
+   WALBERLA_MPI_WORLD_BARRIER()
 }
 
 
@@ -211,7 +211,7 @@ void VTKOutput::forceWrite( uint_t number, const bool immediatelyWriteCollectors
             filesystem::create_directory( tpath );
       }
    }
-   WALBERLA_MPI_WORLD_BARRIER();
+   WALBERLA_MPI_WORLD_BARRIER()
 
 
    if( useMPIIO_ )
@@ -308,7 +308,7 @@ void VTKOutput::write( const bool immediatelyWriteCollectors, const int simultan
 
 void VTKOutput::writeDomainDecomposition( const std::string& path, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates ) const
 {
-   WALBERLA_ASSERT_NOT_NULLPTR( unstructuredBlockStorage_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( unstructuredBlockStorage_ )
 
    const uint_t numberOfBlocks = unstructuredBlockStorage_->getNumberOfBlocks();
 
@@ -328,7 +328,7 @@ void VTKOutput::writeDomainDecomposition( const std::string& path, const Set<SUI
    writeDomainDecompositionPieces( ofs, requiredStates, incompatibleStates );
 
    ofs << " </UnstructuredGrid>\n"
-       << "</VTKFile>" << std::endl;
+       << "</VTKFile>" << '\n';
 
    ofs.close();
 }
@@ -474,7 +474,7 @@ void VTKOutput::writeDomainDecompositionPieces( std::ostream& ofs, const Set<SUI
 void VTKOutput::computeOutputPoints( std::vector<Vector3<real_t> > & points, std::vector<bool> & outputPoint,
                                      uint_t & numberOfPoints ) const
 {
-   WALBERLA_ASSERT_NOT_NULLPTR( pointDataSource_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( pointDataSource_ )
 
    pointDataSource_->configure();
    points = pointDataSource_->getPoints();
@@ -574,7 +574,7 @@ void VTKOutput::writePointDataPieceHelper( const std::vector<Vector3<real_t> > &
    for( auto dataArray = attributes.begin(); dataArray != attributes.end(); ++dataArray, ++dataIndex )
    {
       const uint_t components = dataArray->components;
-      WALBERLA_ASSERT_GREATER( components, 0 );
+      WALBERLA_ASSERT_GREATER( components, 0 )
 
       ofs << "    <DataArray type=\"" << dataArray->type << "\" Name=\"" << dataArray->name << "\" NumberOfComponents=\""
          << components << "\" format=\"" << format_ << "\">\n";
@@ -637,7 +637,7 @@ void VTKOutput::writePointData( const std::string& path )
    writePointDataPieceHelper( points, outputPoint, numberOfPoints, ofs );
 
    ofs << " </UnstructuredGrid>\n"
-       << "</VTKFile>" << std::endl;
+       << "</VTKFile>" << '\n';
 
    ofs.close();
 }
@@ -664,7 +664,7 @@ void VTKOutput::computeOutputPolylines( std::vector< std::vector< Vector3< real_
    std::vector< std::vector< bool > > & outputPolylinePoint, std::vector< size_t > & polylineSize,
    uint_t & numberOfPolylines, uint_t & numberOfPolylinePoints ) const
 {
-   WALBERLA_ASSERT_NOT_NULLPTR( polylineDataSource_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( polylineDataSource_ )
 
    polylineDataSource_->configure();
    lines = polylineDataSource_->getPolyLines();
@@ -861,7 +861,7 @@ void VTKOutput::writePolylineDataPieceHelper( const std::vector< std::vector< Ve
    for( auto dataArray = attributes.begin(); dataArray != attributes.end(); ++dataArray, ++dataIndex )
    {
       const uint_t components = dataArray->components;
-      WALBERLA_ASSERT_GREATER( components, 0 );
+      WALBERLA_ASSERT_GREATER( components, 0 )
 
       ofs << "    <DataArray type=\"" << dataArray->type << "\" Name=\"" << dataArray->name << "\" NumberOfComponents=\""
          << components << "\" format=\"" << format_ << "\">\n";
@@ -932,7 +932,7 @@ void VTKOutput::writePolylineData( const std::string& path )
    writePolylineDataPieceHelper( lines, outputPolylinePoint, polylineSize, numberOfPolylines, numberOfPolylinePoints, ofs );
 
    ofs << " </UnstructuredGrid>\n"
-       << "</VTKFile>" << std::endl;
+       << "</VTKFile>" << '\n';
 
    ofs.close();
 }
@@ -1007,19 +1007,18 @@ void VTKOutput::computeVTUCells( const IBlock& block, CellVector & cellsOut ) co
 
 void VTKOutput::writeBlocks( const std::string& path, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates )
 {
-   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ )
 
    if( !configured_ ) {
       if( !forcePVTU_ && cellInclusionFunctions_.empty() && cellExclusionFunctions_.empty() &&
           blockStorage_->getNumberOfLevels() == 1 && ghostLayers_ == 0 ) // uniform data -> vti
       {
-            uniformGrid_ = true;
+         uniformGrid_ = true;
       }
       configured_ = true;
    }
 
-   if(!uniformGrid_ && oneFilePerProcess_)
-   {
+   if(!uniformGrid_ && oneFilePerProcess_){
       const int rank         = MPIManager::instance()->rank();
       std::ostringstream file;
       file << path << "/dataRank[" << rank << "].vtu";
@@ -1027,49 +1026,18 @@ void VTKOutput::writeBlocks( const std::string& path, const Set<SUID>& requiredS
       writeParallelVTU( ofs, requiredStates, incompatibleStates );
       ofs.close();
    }
-   else
-   {
-      std::vector< const IBlock* > blocks;
-      for( auto block = blockStorage_->begin(); block != blockStorage_->end(); ++block )
-      {
-         if( selectable::isSetSelected( uid::globalState() + block->getState(), requiredStates, incompatibleStates ) )
-            blocks.push_back( block.get() );
+   else{
+      if (uniformGrid_ || amrFileFormat_){ // uniform data -> vti  amr data -> vti
+         if (samplingDx_ <= real_c(0) || samplingDy_ <= real_c(0) || samplingDz_ <= real_c(0))
+            writeVTI(path, requiredStates, incompatibleStates);
+         else
+            writeVTI_sampling(path, requiredStates, incompatibleStates);
       }
-      for( auto it = blocks.begin(); it != blocks.end(); ++it )
-      {
-         WALBERLA_ASSERT_NOT_NULLPTR( *it );
-         const IBlock& block = **it;
-         const uint_t level = blockStorage_->getLevel(block);
-
-         std::ostringstream file;
-         file << path << "/block [" << block.getId() << "] level[" << level << "].";
-
-         if( uniformGrid_ || amrFileFormat_ ) // uniform data -> vti  amr data -> vti
-         {
-            file << "vti";
-            std::ofstream ofs( file.str().c_str()  );
-            if( samplingDx_ <= real_c(0) || samplingDy_ <= real_c(0) || samplingDz_ <= real_c(0) )
-               writeVTI( ofs, block );
-            else
-               writeVTI_sampling( ofs, block );
-            ofs.close();
-         }
-         else // unstructured data -> vtu
-         {
-            CellVector cells; // cells to be written to file
-            computeVTUCells( block, cells );
-
-            if( !cells.empty() )
-            {
-               file << "vtu";
-               std::ofstream ofs( file.str().c_str()  );
-               if( samplingDx_ <= real_c(0) || samplingDy_ <= real_c(0) || samplingDz_ <= real_c(0) )
-                  writeVTU( ofs, block, cells );
-               else
-                  writeVTU_sampling( ofs, block, cells );
-               ofs.close();
-            }
-         }
+      else{ // unstructured data -> vtu
+         if (samplingDx_ <= real_c(0) || samplingDy_ <= real_c(0) || samplingDz_ <= real_c(0))
+            writeVTU(path, requiredStates, incompatibleStates);
+         else
+            writeVTU_sampling(path, requiredStates, incompatibleStates);
       }
    }
 }
@@ -1078,7 +1046,7 @@ void VTKOutput::writeBlocks( const std::string& path, const Set<SUID>& requiredS
 
 void VTKOutput::writeBlockPieces( std::ostream & oss, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates )
 {
-   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ )
 
    std::vector< const IBlock* > blocks;
    for( auto block = blockStorage_->begin(); block != blockStorage_->end(); ++block )
@@ -1098,7 +1066,7 @@ void VTKOutput::writeBlockPieces( std::ostream & oss, const Set<SUID>& requiredS
 
    for( auto it = blocks.begin(); it != blocks.end(); ++it )
    {
-      WALBERLA_ASSERT_NOT_NULLPTR( *it );
+      WALBERLA_ASSERT_NOT_NULLPTR( *it )
       const IBlock& block = **it;
 
       if( uniformGrid_ ) // uniform data -> vti
@@ -1124,32 +1092,71 @@ void VTKOutput::writeBlockPieces( std::ostream & oss, const Set<SUID>& requiredS
    }
 }
 
+void VTKOutput::writeVTI(const std::string& path, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates) const
+{
+   const AABB& domain = blockStorage_->getDomain();
+   std::vector< const IBlock* > blocks;
+
+   for( auto block = blockStorage_->begin(); block != blockStorage_->end(); ++block ){
+      if( selectable::isSetSelected( uid::globalState() + block->getState(), requiredStates, incompatibleStates ) )
+         blocks.push_back( block.get() );
+   }
+   for( auto it = blocks.begin(); it != blocks.end(); ++it ){
+      WALBERLA_ASSERT_NOT_NULLPTR(*it)
+      const IBlock& block = **it;
+      const uint_t level  = blockStorage_->getLevel(block);
 
+      CellVector cells; // cells to be written to file
+      computeVTUCells( block, cells );
+      if(cells.empty())
+         continue;
 
-void VTKOutput::writeVTI( std::ostream& ofs, const IBlock& block ) const
+      CellInterval localCellInterval = cells.boundingBox();
+      CellInterval globalCellInterval;
+      blockStorage_->transformBlockLocalToGlobalCellInterval(globalCellInterval, block, localCellInterval);
+
+      std::ostringstream file;
+      file << path << "/block [" << block.getId() << "] level[" << level << "].";
+      file << "vti";
+      std::ofstream ofs(file.str().c_str());
+
+      ofs << "<?xml version=\"1.0\"?>\n"
+          << "<VTKFile type=\"ImageData\" version=\"0.1\" byte_order=\"" << endianness_ << "\">\n"
+          << " <ImageData WholeExtent=\"" << globalCellInterval.xMin() << " " << (globalCellInterval.xMax() + 1) << " " << globalCellInterval.yMin() << " "
+          << (globalCellInterval.yMax() + 1) << " " << globalCellInterval.zMin() << " " << (globalCellInterval.zMax() + 1) << "\"" << " Origin=\""
+          << domain.xMin() << " " << domain.yMin() << " " << domain.zMin() << "\"" << " Spacing=\""
+          << blockStorage_->dx(level) << " " << blockStorage_->dy(level) << " " << blockStorage_->dz(level) << "\">\n";
+
+      writeVTIPiece(ofs, block, globalCellInterval, localCellInterval);
+
+      ofs << " </ImageData>\n"
+          << "</VTKFile>" << '\n';
+      ofs.close();
+   }
+}
+
+void VTKOutput::writeVTIPiece( std::ostream& ofs, const IBlock& block, const CellInterval& globalCellInterval, CellInterval& localCellInterval) const
 {
-   const CellInterval& cellBB = blockStorage_->getBlockCellBB( block );
-   const AABB&         domain = blockStorage_->getDomain();
-   const uint_t        level  = blockStorage_->getLevel( block );
+   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ )
+   WALBERLA_ASSERT_EQUAL( ghostLayers_, 0 )
 
-   ofs << "<?xml version=\"1.0\"?>\n"
-       << "<VTKFile type=\"ImageData\" version=\"0.1\" byte_order=\"" << endianness_ << "\">\n"
-       << " <ImageData WholeExtent=\"" << cellBB.xMin() << " " << ( cellBB.xMax() + 1 ) << " "
-       << cellBB.yMin() << " " << ( cellBB.yMax() + 1 ) << " "
-       << cellBB.zMin() << " " << ( cellBB.zMax() + 1 ) << "\""
-       << " Origin=\"" << domain.xMin() << " " << domain.yMin() << " " << domain.zMin() << "\""
-       << " Spacing=\"" << blockStorage_->dx(level) << " " << blockStorage_->dy(level) << " " << blockStorage_->dz(level) << "\">\n";
+   ofs << "  <Piece Extent=\"" << globalCellInterval.xMin() << " " << ( globalCellInterval.xMax() + 1 ) << " "
+       << globalCellInterval.yMin() << " " << ( globalCellInterval.yMax() + 1 ) << " "
+       << globalCellInterval.zMin() << " " << ( globalCellInterval.zMax() + 1 ) << "\">\n"
+       << "   <CellData>\n";
 
-   writeVTIPiece( ofs, block );
+   writeCellData( ofs, block, localCellInterval );
 
-   ofs << " </ImageData>\n"
-       << "</VTKFile>" << std::endl;
+   ofs << "   </CellData>\n"
+       << "  </Piece>\n";
 }
 
+
+
 void VTKOutput::writeVTIPiece( std::ostream& ofs, const IBlock& block ) const
 {
-   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ );
-   WALBERLA_ASSERT_EQUAL( ghostLayers_, 0 );
+   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ )
+   WALBERLA_ASSERT_EQUAL( ghostLayers_, 0 )
 
    const CellInterval& cellBB = blockStorage_->getBlockCellBB( block );
 
@@ -1165,41 +1172,60 @@ void VTKOutput::writeVTIPiece( std::ostream& ofs, const IBlock& block ) const
          for( uint_t x = 0; x < blockStorage_->getNumberOfXCells( block ); ++x )
             cells.push_back( x, y, z );
 
+   WALBERLA_LOG_DEVEL_VAR(cells)
+
    writeCellData( ofs, block, cells );
 
    ofs << "   </CellData>\n"
       << "  </Piece>\n";
 }
 
-void VTKOutput::writeVTI_sampling( std::ostream& ofs, const IBlock& block ) const
+void VTKOutput::writeVTI_sampling(const std::string& path, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates) const
 {
-   const AABB&  blockBB = block.getAABB();
-   const AABB&  domain  = blockStorage_->getDomain();
+   std::vector< const IBlock* > blocks;
+   for( auto block = blockStorage_->begin(); block != blockStorage_->end(); ++block ){
+      if( selectable::isSetSelected( uid::globalState() + block->getState(), requiredStates, incompatibleStates ) )
+         blocks.push_back( block.get() );
+   }
+   for( auto it = blocks.begin(); it != blocks.end(); ++it )
+   {
+      WALBERLA_ASSERT_NOT_NULLPTR(*it)
+      const IBlock& block = **it;
+      const uint_t level  = blockStorage_->getLevel(block);
 
-   CellInterval cellBB = getSampledCellInterval( blockBB );
+      std::ostringstream file;
+      file << path << "/block [" << block.getId() << "] level[" << level << "].";
+      file << "vti";
+      std::ofstream ofs(file.str().c_str());
 
-   ofs << "<?xml version=\"1.0\"?>\n"
-       << "<VTKFile type=\"ImageData\" version=\"0.1\" byte_order=\"" << endianness_ << "\">\n"
-       << " <ImageData WholeExtent=\"" << cellBB.xMin() << " " << ( cellBB.xMax() + 1 ) << " "
-       << cellBB.yMin() << " " << ( cellBB.yMax() + 1 ) << " "
-       << cellBB.zMin() << " " << ( cellBB.zMax() + 1 ) << "\""
-       << " Origin=\"" << domain.xMin() << " " << domain.yMin() << " " << domain.zMin() << "\""
-       << " Spacing=\"" << samplingDx_ << " " << samplingDy_ << " " << samplingDz_ << "\">\n";
+      const AABB& blockBB = block.getAABB();
+      const AABB& domain  = blockStorage_->getDomain();
+
+      CellInterval cellBB = getSampledCellInterval(blockBB);
 
-   writeVTIPiece_sampling( ofs, block );
+      ofs << "<?xml version=\"1.0\"?>\n"
+          << "<VTKFile type=\"ImageData\" version=\"0.1\" byte_order=\"" << endianness_ << "\">\n"
+          << " <ImageData WholeExtent=\"" << cellBB.xMin() << " " << (cellBB.xMax() + 1) << " " << cellBB.yMin() << " "
+          << (cellBB.yMax() + 1) << " " << cellBB.zMin() << " " << (cellBB.zMax() + 1) << "\"" << " Origin=\""
+          << domain.xMin() << " " << domain.yMin() << " " << domain.zMin() << "\"" << " Spacing=\"" << samplingDx_
+          << " " << samplingDy_ << " " << samplingDz_ << "\">\n";
+
+      writeVTIPiece_sampling(ofs, block);
 
-   ofs << " </ImageData>\n"
-       << "</VTKFile>" << std::endl;
+      ofs << " </ImageData>\n"
+          << "</VTKFile>" << '\n';
+      ofs.close();
+   }
 }
 
 void VTKOutput::writeVTIPiece_sampling( std::ostream& ofs, const IBlock& block ) const
 {
-   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ );
-   WALBERLA_ASSERT_EQUAL( ghostLayers_, 0 );
+   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ )
+   WALBERLA_ASSERT_EQUAL( ghostLayers_, 0 )
 
-   WALBERLA_ASSERT_GREATER( samplingDx_, real_t(0) );
-   WALBERLA_ASSERT_GREATER( samplingDy_, real_t(0) );
-   WALBERLA_ASSERT_GREATER( samplingDz_, real_t(0) );
+   WALBERLA_ASSERT_GREATER( samplingDx_, real_t(0) )
+   WALBERLA_ASSERT_GREATER( samplingDy_, real_t(0) )
+   WALBERLA_ASSERT_GREATER( samplingDz_, real_t(0) )
 
    const AABB&  blockBB = block.getAABB();
    const AABB&  domain = blockStorage_->getDomain();
@@ -1234,16 +1260,40 @@ void VTKOutput::writeVTIPiece_sampling( std::ostream& ofs, const IBlock& block )
 
 
 
-void VTKOutput::writeVTU( std::ostream& ofs, const IBlock& block, const CellVector& cells ) const
+void VTKOutput::writeVTU(const std::string& path, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates) const
 {
-   ofs << "<?xml version=\"1.0\"?>\n"
-       << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"" << endianness_ << "\">\n"
-       << " <UnstructuredGrid>\n";
+   std::vector< const IBlock* > blocks;
+   for( auto block = blockStorage_->begin(); block != blockStorage_->end(); ++block ){
+      if( selectable::isSetSelected( uid::globalState() + block->getState(), requiredStates, incompatibleStates ) )
+         blocks.push_back( block.get() );
+   }
+   for( auto it = blocks.begin(); it != blocks.end(); ++it ){
+      WALBERLA_ASSERT_NOT_NULLPTR(*it)
+      const IBlock& block = **it;
+      const uint_t level  = blockStorage_->getLevel(block);
 
-   writeVTUPiece( ofs, block, cells );
+      std::ostringstream file;
+      file << path << "/block [" << block.getId() << "] level[" << level << "].";
 
-   ofs << " </UnstructuredGrid>\n"
-       << "</VTKFile>" << std::endl;
+      CellVector cells; // cells to be written to file
+      computeVTUCells( block, cells );
+
+      if( !cells.empty() ){
+         file << "vtu";
+         std::ofstream ofs(file.str().c_str());
+
+         ofs << "<?xml version=\"1.0\"?>\n"
+             << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"" << endianness_ << "\">\n"
+             << " <UnstructuredGrid>\n";
+
+         writeVTUPiece(ofs, block, cells);
+
+         ofs << " </UnstructuredGrid>\n"
+             << "</VTKFile>" << '\n';
+
+         ofs.close();
+      }
+   }
 }
 
 void VTKOutput::writeParallelVTU( std::ostream& ofs, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates  ) const
@@ -1252,7 +1302,7 @@ void VTKOutput::writeParallelVTU( std::ostream& ofs, const Set<SUID>& requiredSt
        << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"" << endianness_ << "\">\n"
        << " <UnstructuredGrid>\n";
 
-   WALBERLA_ASSERT_NOT_NULLPTR(blockStorage_);
+   WALBERLA_ASSERT_NOT_NULLPTR(blockStorage_)
    const uint_t finestLevel = blockStorage_->getNumberOfLevels() - 1;
 
    std::map< Vertex, Index, VertexCompare > vimap; // vertex<->index map
@@ -1260,34 +1310,30 @@ void VTKOutput::writeParallelVTU( std::ostream& ofs, const Set<SUID>& requiredSt
    std::vector< Index >                     ci;    // ci[0] to ci[7]: indices for cell number one, ci[8] to ci[15]: ...
    uint_t numberOfCells = 0;
 
-   for( auto block = blockStorage_->begin(); block != blockStorage_->end(); ++block )
-   {
+   for( auto block = blockStorage_->begin(); block != blockStorage_->end(); ++block ){
       if( !selectable::isSetSelected( uid::globalState() + block->getState(), requiredStates, incompatibleStates ) )
          continue;
 
-      // CellVector cells; // cells to be written to file
-      // computeVTUCells( *block, cells );
+      CellVector cells; // cells to be written to file
+      computeVTUCells( *block, cells );
 
       const uint_t level = blockStorage_->getLevel(*block);
       const cell_idx_t factorToFinest = 1 << (finestLevel - level);
-      const CellInterval cells = blockStorage_->getBlockCellBB(*block); //  These are global cells
 
-      for (auto cell = cells.begin(); cell != cells.end(); ++cell)
-      {
+      for (auto cell = cells.begin(); cell != cells.end(); ++cell){
          numberOfCells++;
-         const AABB aabb = blockStorage_->getCellAABB(*cell, level);
+         Cell globalCell;
+         blockStorage_->transformBlockLocalToGlobalCell(globalCell, *block, *cell);
+         const AABB aabb = blockStorage_->getCellAABB(globalCell, level);
          for (cell_idx_t z = 0; z != 2; ++z) {
             for (cell_idx_t y = 0; y != 2; ++y) {
-               for (cell_idx_t x = 0; x != 2; ++x)
-               {
-                  const Vertex v((cell->x() + x) * factorToFinest, (cell->y() + y) * factorToFinest, (cell->z() + z) * factorToFinest);
+               for (cell_idx_t x = 0; x != 2; ++x){
+                  const Vertex v((globalCell.x() + x) * factorToFinest, (globalCell.y() + y) * factorToFinest, (globalCell.z() + z) * factorToFinest);
                   auto mapping = vimap.find(v);
-                  if (mapping != vimap.end()) // vertex already exists
-                  {
+                  if (mapping != vimap.end()){ // vertex already exists
                      ci.push_back(mapping->second);
                   }
-                  else // new vertex
-                  {
+                  else{ // new vertex
                      vimap[v] = numeric_cast< Index >(vc.size());
                      ci.push_back(numeric_cast< Index >(vc.size()));
                      vc.emplace_back((x == 0) ? aabb.xMin() : aabb.xMax(),
@@ -1310,15 +1356,13 @@ void VTKOutput::writeParallelVTU( std::ostream& ofs, const Set<SUID>& requiredSt
        << "  </Piece>\n";
 
    ofs << " </UnstructuredGrid>\n"
-       << "</VTKFile>" << std::endl;
-
+       << "</VTKFile>" << '\n';
 }
 
 
-
 void VTKOutput::writeVTUPiece( std::ostream& ofs, const IBlock& block, const CellVector& cells ) const
 {
-   WALBERLA_ASSERT_NOT_NULLPTR(blockStorage_);
+   WALBERLA_ASSERT_NOT_NULLPTR(blockStorage_)
 
    // setting up vertex-index mapping --->
 
@@ -1355,7 +1399,7 @@ void VTKOutput::writeVTUPiece( std::ostream& ofs, const IBlock& block, const Cel
       }
    }
 
-   WALBERLA_ASSERT_EQUAL(ci.size(), 8 * cells.size());
+   WALBERLA_ASSERT_EQUAL(ci.size(), 8 * cells.size())
 
    // <--- setting up vertex-index mapping
 
@@ -1394,17 +1438,42 @@ void VTKOutput::writeVTUPiece( std::ostream& ofs, const IBlock& block, const Cel
 
 
 
-void VTKOutput::writeVTU_sampling( std::ostream& ofs, const IBlock& block, const CellVector& localCells ) const
+void VTKOutput::writeVTU_sampling(const std::string& path, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates) const
 {
-   ofs << "<?xml version=\"1.0\"?>\n"
-       << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"" << endianness_ << "\">\n"
-       << " <UnstructuredGrid>\n";
+   std::vector< const IBlock* > blocks;
+   for( auto block = blockStorage_->begin(); block != blockStorage_->end(); ++block ){
+      if( selectable::isSetSelected( uid::globalState() + block->getState(), requiredStates, incompatibleStates ) )
+         blocks.push_back( block.get() );
+   }
+   for( auto it = blocks.begin(); it != blocks.end(); ++it )
+   {
+      WALBERLA_ASSERT_NOT_NULLPTR(*it)
+      const IBlock& block = **it;
+      const uint_t level  = blockStorage_->getLevel(block);
 
-   writeVTUPiece_sampling( ofs, block, localCells );
+      std::ostringstream file;
+      file << path << "/block [" << block.getId() << "] level[" << level << "].";
 
-   ofs << " </UnstructuredGrid>\n"
-       << "</VTKFile>" << std::endl;
+      CellVector cells; // cells to be written to file
+      computeVTUCells(block, cells);
+
+      if (!cells.empty())
+      {
+         file << "vtu";
+         std::ofstream ofs(file.str().c_str());
+
+         ofs << "<?xml version=\"1.0\"?>\n"
+             << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"" << endianness_ << "\">\n"
+             << " <UnstructuredGrid>\n";
+
+         writeVTUPiece_sampling(ofs, block, cells);
 
+         ofs << " </UnstructuredGrid>\n"
+             << "</VTKFile>" << '\n';
+
+         ofs.close();
+      }
+   }
 }
 
 
@@ -1444,7 +1513,7 @@ void VTKOutput::writeVTUPiece_sampling(std::ostream& ofs, const IBlock& block, c
       }
    }
 
-   WALBERLA_ASSERT_EQUAL(ci.size(), 8 * cells.size());
+   WALBERLA_ASSERT_EQUAL(ci.size(), 8 * cells.size())
 
    // <--- setting up vertex-index mapping
 
@@ -1541,7 +1610,7 @@ void VTKOutput::writeVTUHeaderPiece( std::ostream& ofs, const uint_t numberOfCel
 
 uint8_t VTKOutput::ghostLayerNr( const IBlock& block, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
 {
-   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ )
 
    const cell_idx_t xMax = cell_idx_c( blockStorage_->getNumberOfXCells(block) ) - cell_idx_c(1);
    const cell_idx_t yMax = cell_idx_c( blockStorage_->getNumberOfYCells(block) ) - cell_idx_c(1);
@@ -1550,17 +1619,17 @@ uint8_t VTKOutput::ghostLayerNr( const IBlock& block, const cell_idx_t x, const
    cell_idx_t xDistance = 0;
    if( x < 0 ) xDistance = -x;
    else if( x > xMax ) xDistance = x - xMax;
-   WALBERLA_ASSERT_GREATER_EQUAL( xDistance, 0 );
+   WALBERLA_ASSERT_GREATER_EQUAL( xDistance, 0 )
 
    cell_idx_t yDistance = 0;
    if( y < 0 ) yDistance = -y;
    else if( y > yMax ) yDistance = y - yMax;
-   WALBERLA_ASSERT_GREATER_EQUAL( xDistance, 0 );
+   WALBERLA_ASSERT_GREATER_EQUAL( xDistance, 0 )
 
    cell_idx_t zDistance = 0;
    if( z < 0 ) zDistance = -z;
    else if( z > zMax ) zDistance = z - zMax;
-   WALBERLA_ASSERT_GREATER_EQUAL( xDistance, 0 );
+   WALBERLA_ASSERT_GREATER_EQUAL( xDistance, 0 )
 
    return uint8_c( std::max( xDistance, std::max( yDistance, zDistance )) );
 }
@@ -1569,7 +1638,7 @@ uint8_t VTKOutput::ghostLayerNr( const IBlock& block, const cell_idx_t x, const
 
 std::vector< VTKOutput::SamplingCell > VTKOutput::getSamplingCells( const IBlock& block, const CellVector& cells ) const
 {
-   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ )
 
    std::vector< SamplingCell > samplingCells;
    CellSet cellSet( cells );
@@ -1630,10 +1699,45 @@ std::vector< VTKOutput::SamplingCell > VTKOutput::getSamplingCells( const IBlock
 }
 
 
+void VTKOutput::writeCellData( std::ostream& ofs, const IBlock& block, const CellInterval& cells ) const
+{
+   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ )
+
+   for( auto writer = cellDataWriter_.begin(); writer != cellDataWriter_.end(); ++writer )
+   {
+      (*writer)->configure( block, *blockStorage_ );
+
+      ofs << "    <DataArray type=\"" << (*writer)->typeString() << "\" Name=\"" << (*writer)->identifier()
+          << "\" NumberOfComponents=\"" << (*writer)->fSize() << "\" format=\"" << format_ << "\">\n";
+
+      if( binary_ )
+      {
+         Base64Writer base64;
+         for( auto cell = cells.begin(); cell != cells.end(); ++cell )
+            for( uint_t f = 0; f != (*writer)->fSize(); ++f )
+               (*writer)->push( base64, cell->x(), cell->y(), cell->z(), cell_idx_c(f) );
+         ofs << "     "; base64.toStream( ofs );
+      }
+      else
+      {
+         for( auto cell = cells.begin(); cell != cells.end(); ++cell ) {
+            ofs << "     ";
+            for( uint_t f = 0; f != (*writer)->fSize(); ++f )
+            {
+               (*writer)->push( ofs, cell->x(), cell->y(), cell->z(), cell_idx_c(f) );
+               ofs << ( ( f == (*writer)->fSize() - 1 ) ? "\n" : " " );
+            }
+         }
+      }
+
+      ofs << "    </DataArray>\n";
+   }
+}
+
 
 void VTKOutput::writeCellData( std::ostream& ofs, const IBlock& block, const CellVector& cells ) const
 {
-   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ )
 
    for( auto writer = cellDataWriter_.begin(); writer != cellDataWriter_.end(); ++writer )
    {
@@ -1669,7 +1773,7 @@ void VTKOutput::writeCellData( std::ostream& ofs, const IBlock& block, const Cel
 
 void VTKOutput::writeCellData( std::ostream& ofs, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates ) const
 {
-   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ )
 
    for( auto writer = cellDataWriter_.begin(); writer != cellDataWriter_.end(); ++writer )
    {
@@ -1712,7 +1816,7 @@ void VTKOutput::writeCellData( std::ostream& ofs, const Set<SUID>& requiredState
 
 void VTKOutput::writeCellData( std::ostream& ofs, const IBlock& block, const std::vector< SamplingCell >& cells ) const
 {
-   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ )
 
    for( auto writer = cellDataWriter_.begin(); writer != cellDataWriter_.end(); ++writer )
    {
@@ -1756,18 +1860,14 @@ void VTKOutput::writeCellData( std::ostream& ofs, const IBlock& block, const std
 void VTKOutput::writeCollectors( const bool barrier )
 {
    if( barrier )
-      WALBERLA_MPI_WORLD_BARRIER();
+      WALBERLA_MPI_WORLD_BARRIER()
 
    WALBERLA_NON_ROOT_SECTION() { return; }
-
-   WALBERLA_ASSERT_EQUAL( MPIManager::instance()->worldRank(), 0 );
+   WALBERLA_ASSERT_EQUAL( MPIManager::instance()->worldRank(), 0 )
 
    if(!amrFileFormat_)
       writePVD();
 
-
-
-
    for( auto collector = collectorsToWrite_.begin(); collector != collectorsToWrite_.end(); ++collector )
    {
       if( uniformGrid_ )
@@ -1816,20 +1916,20 @@ void VTKOutput::writePVD()
       {
          if( line.find("</Collection>") != std::string::npos )
          {
-            WALBERLA_ASSERT_GREATER( ofs.tellg(), 0 );
+            WALBERLA_ASSERT_GREATER( ofs.tellg(), 0 )
             pvdEnd_ = ofs.tellg();
             pvdEnd_ -= int_c(line.size()) + 1;
             break;
          }
       }
-      WALBERLA_ASSERT_GREATER( pvdEnd_, 0 );
+      WALBERLA_ASSERT_GREATER( pvdEnd_, 0 )
       ofs.seekp(pvdEnd_);
    }
    else
    {
       ofs.seekp(pvdEnd_);
    }
-   WALBERLA_ASSERT_GREATER(ofs.tellp(), 0);
+   WALBERLA_ASSERT_GREATER(ofs.tellp(), 0)
 
    std::string ending;
    if( useMPIIO_ )
@@ -1854,7 +1954,7 @@ void VTKOutput::writePVD()
    allCollectors_.clear();
 
    pvdEnd_ = ofs.tellp();
-   WALBERLA_ASSERT_GREATER( pvdEnd_, 0 );
+   WALBERLA_ASSERT_GREATER( pvdEnd_, 0 )
    ofs << " </Collection>\n"
        << "</VTKFile>\n";
 
@@ -1890,14 +1990,14 @@ void VTKOutput::writeVTHBSeries()
             break;
          }
       }
-      WALBERLA_ASSERT_GREATER( pvdEnd_, 0 );
+      WALBERLA_ASSERT_GREATER( pvdEnd_, 0 )
       ofs.seekp(pvdEnd_);
    }
    else
    {
       ofs.seekp(pvdEnd_);
    }
-   WALBERLA_ASSERT_GREATER(ofs.tellp(), 0);
+   WALBERLA_ASSERT_GREATER(ofs.tellp(), 0)
 
    std::string ending = ".vthb";
 
@@ -1910,7 +2010,7 @@ void VTKOutput::writeVTHBSeries()
    allCollectors_.clear();
 
    pvdEnd_ = ofs.tellp();
-   WALBERLA_ASSERT_GREATER( pvdEnd_, 0 );
+   WALBERLA_ASSERT_GREATER( pvdEnd_, 0 )
    ofs << "   ]\n"
        << "}\n";
 
@@ -1921,7 +2021,7 @@ void VTKOutput::writeVTHBSeries()
 
 void VTKOutput::writePVTI( const uint_t collector ) const
 {
-   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ )
 
    std::ostringstream collection;
    collection << baseFolder_ << "/" << identifier_ << "/" << executionFolder_ << "_" << collector << ".pvti";
@@ -1969,7 +2069,7 @@ void VTKOutput::writePVTI( const uint_t collector ) const
 
 void VTKOutput::writeVTHB( const uint_t collector ) const
 {
-   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ )
 
    std::ostringstream collection;
    collection << baseFolder_ << "/" << identifier_ << "/" << executionFolder_ << "_" << collector << ".vthb";
@@ -2002,10 +2102,10 @@ void VTKOutput::writeVTHB( const uint_t collector ) const
 
 void VTKOutput::writePVTI_sampled( const uint_t collector ) const
 {
-   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ );
-   WALBERLA_ASSERT_GREATER( samplingDx_, real_t(0) );
-   WALBERLA_ASSERT_GREATER( samplingDy_, real_t(0) );
-   WALBERLA_ASSERT_GREATER( samplingDz_, real_t(0) );
+   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ )
+   WALBERLA_ASSERT_GREATER( samplingDx_, real_t(0) )
+   WALBERLA_ASSERT_GREATER( samplingDy_, real_t(0) )
+   WALBERLA_ASSERT_GREATER( samplingDz_, real_t(0) )
 
    std::ostringstream collection;
    collection << baseFolder_ << "/" << identifier_ << "/" << executionFolder_ << "_" << collector << ".pvti";
@@ -2052,7 +2152,7 @@ void VTKOutput::writePVTI_sampled( const uint_t collector ) const
 
 bool VTKOutput::writeCombinedVTI( std::string localPart, const uint_t collector ) const
 {
-   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ );
+   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ )
 
    const MPI_Comm comm    = MPIManager::instance()->comm();
    const int rank         = MPIManager::instance()->rank();
@@ -2098,10 +2198,10 @@ bool VTKOutput::writeCombinedVTI( std::string localPart, const uint_t collector
 
 bool VTKOutput::writeCombinedVTI_sampled( std::string localPart, const uint_t collector ) const
 {
-   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ );
-   WALBERLA_ASSERT_GREATER( samplingDx_, real_t( 0 ) );
-   WALBERLA_ASSERT_GREATER( samplingDy_, real_t( 0 ) );
-   WALBERLA_ASSERT_GREATER( samplingDz_, real_t( 0 ) );
+   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ )
+   WALBERLA_ASSERT_GREATER( samplingDx_, real_t( 0 ) )
+   WALBERLA_ASSERT_GREATER( samplingDy_, real_t( 0 ) )
+   WALBERLA_ASSERT_GREATER( samplingDz_, real_t( 0 ) )
 
    const MPI_Comm comm    = MPIManager::instance()->comm();
    const int rank         = MPIManager::instance()->rank();
@@ -2233,11 +2333,11 @@ void VTKOutput::getFilenames( std::vector< filesystem::path >& files, const uint
    path << baseFolder_ << "/" << identifier_ << "/" << executionFolder_ << "_" << collector;
    filesystem::path directory( path.str() );
 
-   WALBERLA_ASSERT( filesystem::exists( directory ) );
+   WALBERLA_ASSERT( filesystem::exists( directory ) )
 
    for( filesystem::directory_iterator file( directory ); file != filesystem::directory_iterator(); ++file )
    {
-      WALBERLA_ASSERT( filesystem::is_regular_file( file->path() ) && !filesystem::is_directory( file->path() ) );
+      WALBERLA_ASSERT( filesystem::is_regular_file( file->path() ) && !filesystem::is_directory( file->path() ) )
       files.push_back( file->path() );
    }
 }
@@ -2249,12 +2349,12 @@ void VTKOutput::getFilenamesSortedByLevel( std::vector< std::vector< filesystem:
    path << baseFolder_ << "/" << identifier_ << "/" << executionFolder_ << "_" << collector;
    filesystem::path directory( path.str() );
 
-   WALBERLA_ASSERT( filesystem::exists( directory ) );
+   WALBERLA_ASSERT( filesystem::exists( directory ) )
 
    for( filesystem::directory_iterator file( directory ); file != filesystem::directory_iterator(); ++file )
    {
       std::string fileName = file->path().string();
-      WALBERLA_ASSERT( filesystem::is_regular_file( file->path() ) && !filesystem::is_directory( file->path() ) );
+      WALBERLA_ASSERT( filesystem::is_regular_file( file->path() ) && !filesystem::is_directory( file->path() ) )
 
       std::size_t pos1 = fileName.find("level[");
       WALBERLA_ASSERT_UNEQUAL(pos1, std::string::npos, "file names of the block data must contain the block level for AMR data files")
@@ -2271,7 +2371,7 @@ void VTKOutput::getFilenamesSortedByLevel( std::vector< std::vector< filesystem:
 
 void VTKOutput::writePPointData( std::ofstream& ofs ) const
 {
-   WALBERLA_ASSERT( !( pointDataSource_ && polylineDataSource_ ) );
+   WALBERLA_ASSERT( !( pointDataSource_ && polylineDataSource_ ) )
 
    if( pointDataSource_ )
    {
@@ -2310,9 +2410,9 @@ void VTKOutput::writePCellData( std::ofstream& ofs ) const
 
 CellInterval VTKOutput::getSampledCellInterval( const AABB & aabb ) const
 {
-   WALBERLA_ASSERT_GREATER( samplingDx_, real_t(0) );
-   WALBERLA_ASSERT_GREATER( samplingDy_, real_t(0) );
-   WALBERLA_ASSERT_GREATER( samplingDz_, real_t(0) );
+   WALBERLA_ASSERT_GREATER( samplingDx_, real_t(0) )
+   WALBERLA_ASSERT_GREATER( samplingDy_, real_t(0) )
+   WALBERLA_ASSERT_GREATER( samplingDz_, real_t(0) )
 
    const AABB& domain = blockStorage_->getDomain();
 

src/vtk/VTKOutput.h

@@ -236,14 +236,15 @@ private:
    void writeBlocks( const std::string& path, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates );
    void writeBlockPieces( std::ostream & oss, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates );
 
-   void writeVTI( std::ostream& ofs, const IBlock& block ) const;
-   void writeVTI_sampling( std::ostream& ofs, const IBlock& block ) const;
+   void writeVTI(const std::string& path, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates) const;
+   void writeVTI_sampling(const std::string& path, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates) const;
 
+   void writeVTIPiece( std::ostream& ofs, const IBlock& block, const CellInterval& globalCellInterval, CellInterval& localCellInterval) const;
    void writeVTIPiece( std::ostream& ofs, const IBlock& block ) const;
    void writeVTIPiece_sampling( std::ostream& ofs, const IBlock& block ) const;
 
-   void writeVTU( std::ostream& ofs, const IBlock& block, const CellVector& cells ) const;
-   void writeVTU_sampling( std::ostream& ofs, const IBlock& block, const CellVector& cells ) const;
+   void writeVTU(const std::string& path, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates) const;
+   void writeVTU_sampling(const std::string& path, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates) const;
    void writeParallelVTU( std::ostream& ofs, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates ) const;
 
    void writeVTUPiece(std::ostream& ofs, const IBlock& block, const CellVector& cells) const;
@@ -256,6 +257,7 @@ private:
 
    std::vector< SamplingCell > getSamplingCells( const IBlock& block, const CellVector& cells ) const;
 
+   void writeCellData( std::ostream& ofs, const IBlock& block, const CellInterval& cells ) const;
    void writeCellData( std::ostream& ofs, const IBlock& block, const CellVector& cells ) const;
    void writeCellData( std::ostream& ofs, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates ) const;
    void writeCellData( std::ostream& ofs, const IBlock& block, const std::vector< SamplingCell >& cells ) const;
@@ -339,7 +341,7 @@ private:
 inline void VTKOutput::setInitialWriteCallsToSkip( const uint_t initialWriteCallsToSkip )
 {
    if( executionCounter_ > 0 )
-      WALBERLA_ABORT( "Setting the number of initial write calls to skip is only possible as long as \"write\" has not yet been called!" );
+      WALBERLA_ABORT( "Setting the number of initial write calls to skip is only possible as long as \"write\" has not yet been called!" )
 
    initialWriteCallsToSkip_ = initialWriteCallsToSkip;
 }
@@ -350,10 +352,10 @@ inline void VTKOutput::addAABBInclusionFilter( const AABB & aabb )
 {
    if( outputDomainDecomposition_ )
       WALBERLA_ABORT( "You are trying to add an AABB inclusion filter to VTKOutput \"" << identifier_ << "\", "
-                      "but this VTKOutput is intended for outputting the domain decomposition, AABB filters won't have any effect." );
+                      "but this VTKOutput is intended for outputting the domain decomposition, AABB filters won't have any effect." )
    if( configured_ )
       WALBERLA_ABORT( "You are trying to add an AABB inclusion filter to VTKOutput \"" << identifier_ << "\", "
-                      "but this VTKOutput is already configured.\nYou can only add filters as long as no output has been written." );
+                      "but this VTKOutput is already configured.\nYou can only add filters as long as no output has been written." )
 
    if( pointDataSource_ || polylineDataSource_ )
       aabbInclusionFilters_.push_back( aabb );
@@ -367,10 +369,10 @@ inline void VTKOutput::addAABBExclusionFilter( const AABB & aabb )
 {
    if( outputDomainDecomposition_ )
       WALBERLA_ABORT( "You are trying to add an AABB exclusion filter to VTKOutput \"" << identifier_ << "\", "
-                      "but this VTKOutput is intended for outputting the domain decomposition, AABB filters won't have any effect." );
+                      "but this VTKOutput is intended for outputting the domain decomposition, AABB filters won't have any effect." )
    if( configured_ )
       WALBERLA_ABORT( "You are trying to add an AABB exclusion filter to VTKOutput \"" << identifier_ << "\", "
-                      "but this VTKOutput is already configured.\nYou can only add filters as long as no output has been written." );
+                      "but this VTKOutput is already configured.\nYou can only add filters as long as no output has been written." )
 
    if( pointDataSource_ || polylineDataSource_ )
       aabbExclusionFilters_.push_back( aabb );
@@ -384,18 +386,18 @@ inline void VTKOutput::addCellInclusionFilter( CellFilter f )
 {
    if( outputDomainDecomposition_ )
       WALBERLA_ABORT( "You are trying to add a cell inclusion filter to VTKOutput \"" << identifier_ << "\", "
-                      "but this VTKOutput is intended for outputting the domain decomposition, cell filters won't have any effect." );
+                      "but this VTKOutput is intended for outputting the domain decomposition, cell filters won't have any effect." )
    if( pointDataSource_ )
       WALBERLA_ABORT( "You are trying to add a cell inclusion filter to VTKOutput \"" << identifier_ << "\", "
-                      "but this VTKOutput is intended for outputting arbitrary point data, cell filters won't have any effect." );
+                      "but this VTKOutput is intended for outputting arbitrary point data, cell filters won't have any effect." )
 
    if( polylineDataSource_ )
       WALBERLA_ABORT( "You are trying to add a cell inclusion filter to VTKOutput \"" << identifier_ << "\", "
-                      "but this VTKOutput is intended for outputting arbitrary polyline data, cell filters won't have any effect." );
+                      "but this VTKOutput is intended for outputting arbitrary polyline data, cell filters won't have any effect." )
 
    if( configured_ )
       WALBERLA_ABORT( "You are trying to add a cell inclusion filter to VTKOutput \"" << identifier_ << "\", "
-                      "but this VTKOutput is already configured.\nYou can only add filters as long as no output has been written." );
+                      "but this VTKOutput is already configured.\nYou can only add filters as long as no output has been written." )
 
    cellInclusionFunctions_.push_back(f);
 }
@@ -406,18 +408,18 @@ inline void VTKOutput::addCellExclusionFilter( CellFilter f )
 {
    if( outputDomainDecomposition_ )
       WALBERLA_ABORT( "You are trying to add a cell exclusion filter to VTKOutput \"" << identifier_ << "\", "
-                      "but this VTKOutput is intended for outputting the domain decomposition, cell filters won't have any effect." );
+                      "but this VTKOutput is intended for outputting the domain decomposition, cell filters won't have any effect." )
    if( pointDataSource_ )
       WALBERLA_ABORT( "You are trying to add a cell exclusion filter to VTKOutput \"" << identifier_ << "\", "
-                      "but this VTKOutput is intended for outputting arbitrary point data, cell filters won't have any effect." );
+                      "but this VTKOutput is intended for outputting arbitrary point data, cell filters won't have any effect." )
 
    if( polylineDataSource_ )
       WALBERLA_ABORT( "You are trying to add a cell exclusion filter to VTKOutput \"" << identifier_ << "\", "
-                      "but this VTKOutput is intended for outputting arbitrary polyline data, cell filters won't have any effect." );
+                      "but this VTKOutput is intended for outputting arbitrary polyline data, cell filters won't have any effect." )
 
    if( configured_ )
       WALBERLA_ABORT( "You are trying to add a cell exclusion filter to VTKOutput \"" << identifier_ << "\", "
-                      "but this VTKOutput is already configured.\nYou can only add filters as long as no output has been written." );
+                      "but this VTKOutput is already configured.\nYou can only add filters as long as no output has been written." )
 
    cellExclusionFunctions_.push_back(f);
 }
@@ -428,18 +430,18 @@ inline void VTKOutput::addCellDataWriter( const shared_ptr< BlockCellDataWriterI
 {
    if( outputDomainDecomposition_ )
       WALBERLA_ABORT( "You are trying to add a block cell data writer to VTKOutput \"" << identifier_ << "\", "
-                      "but this VTKOutput is intended for outputting the domain decomposition, cell data writers won't have any effect." );
+                      "but this VTKOutput is intended for outputting the domain decomposition, cell data writers won't have any effect." )
    if( pointDataSource_ )
       WALBERLA_ABORT( "You are trying to add a block cell data writer to VTKOutput \"" << identifier_ << "\", "
-                      "but this VTKOutput is intended for outputting arbitrary point data, cell data writers won't have any effect." );
+                      "but this VTKOutput is intended for outputting arbitrary point data, cell data writers won't have any effect." )
 
    if( polylineDataSource_ )
       WALBERLA_ABORT( "You are trying to add a block cell data writer to VTKOutput \"" << identifier_ << "\", "
-                      "but this VTKOutput is intended for outputting arbitrary polyline data, cell data writers won't have any effect." );
+                      "but this VTKOutput is intended for outputting arbitrary polyline data, cell data writers won't have any effect." )
 
    if( configured_ )
       WALBERLA_ABORT( "You are trying to add a block cell data writer to VTKOutput \"" << identifier_ << "\", "
-                      "but this VTKOutput is already configured.\nYou can only add cell data writers as long as no output has been written." );
+                      "but this VTKOutput is already configured.\nYou can only add cell data writers as long as no output has been written." )
 
    cellDataWriter_.push_back( writer );
 }
@@ -452,22 +454,22 @@ inline void VTKOutput::setSamplingResolution( const real_t spacing )
    {
       if( outputDomainDecomposition_ )
          WALBERLA_ABORT( "You are trying to set a sampling resolution of " << spacing << " to VTKOutput \"" << identifier_ << "\", "
-                         "but this VTKOutput is intended for outputting the domain decomposition, setting a sampling resolution won't have any effect." );
+                         "but this VTKOutput is intended for outputting the domain decomposition, setting a sampling resolution won't have any effect." )
       if( pointDataSource_ )
          WALBERLA_ABORT( "You are trying to set a sampling resolution of " << spacing << " to VTKOutput \"" << identifier_ << "\", "
-                         "but this VTKOutput is intended for outputting arbitrary point data, setting a sampling resolution won't have any effect." );
+                         "but this VTKOutput is intended for outputting arbitrary point data, setting a sampling resolution won't have any effect." )
 
       if( polylineDataSource_ )
          WALBERLA_ABORT( "You are trying to set a sampling resolution of " << spacing << " to VTKOutput \"" << identifier_ << "\", "
-                         "but this VTKOutput is intended for outputting arbitrary polyline data, setting a sampling resolution won't have any effect." );
+                         "but this VTKOutput is intended for outputting arbitrary polyline data, setting a sampling resolution won't have any effect." )
 
       if( configured_ )
          WALBERLA_ABORT( "You are trying to set a sampling resolution of " << spacing << " to VTKOutput \"" << identifier_ << "\", "
-                         "but this VTKOutput is already configured.\nYou can only change the sampling resolution as long as no output has been written." );
+                         "but this VTKOutput is already configured.\nYou can only change the sampling resolution as long as no output has been written." )
       if( ghostLayers_ > 0 )
          WALBERLA_ABORT( "You are trying to set a sampling resolution of " << spacing << " to VTKOutput \"" << identifier_ << "\", "
                          "but this VTKOutput is configured to output " << ghostLayers_ << " ghost layer(s).\n"
-                         "Writing ghost layers and using sampling at the same time is not supported!" );
+                         "Writing ghost layers and using sampling at the same time is not supported!" )
    }
 
    samplingDx_ = spacing;
@@ -484,25 +486,25 @@ inline void VTKOutput::setSamplingResolution( const real_t dx, const real_t dy,
       if( outputDomainDecomposition_ )
          WALBERLA_ABORT( "You are trying to set a sampling resolution of ( " << dx << ", " << dy << ", " << dz << " ) "
                          "to VTKOutput \"" << identifier_ << "\", "
-                         "but this VTKOutput is intended for outputting the domain decomposition, setting a sampling resolution won't have any effect." );
+                         "but this VTKOutput is intended for outputting the domain decomposition, setting a sampling resolution won't have any effect." )
       if( pointDataSource_ )
          WALBERLA_ABORT( "You are trying to set a sampling resolution of ( " << dx << ", " << dy << ", " << dz << " ) "
                          "to VTKOutput \"" << identifier_ << "\", "
-                         "but this VTKOutput is intended for outputting arbitrary point data, setting a sampling resolution won't have any effect." );
+                         "but this VTKOutput is intended for outputting arbitrary point data, setting a sampling resolution won't have any effect." )
 
       if( polylineDataSource_ )
          WALBERLA_ABORT( "You are trying to set a sampling resolution of ( " << dx << ", " << dy << ", " << dz << " ) "
                          "to VTKOutput \"" << identifier_ << "\", "
-                         "but this VTKOutput is intended for outputting arbitrary polyline data, setting a sampling resolution won't have any effect." );
+                         "but this VTKOutput is intended for outputting arbitrary polyline data, setting a sampling resolution won't have any effect." )
 
       if( configured_ )
          WALBERLA_ABORT( "You are trying to set a sampling resolution of ( " << dx << ", " << dy << ", " << dz << " ) "
                          "to VTKOutput \"" << identifier_ << "\", "
-                         "but this VTKOutput is already configured.\nYou can only change the sampling resolution as long as no output has been written." );
+                         "but this VTKOutput is already configured.\nYou can only change the sampling resolution as long as no output has been written." )
       if( ghostLayers_ > 0 )
          WALBERLA_ABORT( "You are trying to set a sampling resolution of ( " << dx << ", " << dy << ", " << dz << " ) "
                          "to VTKOutput \"" << identifier_ << "\", but this VTKOutput is configured to output " << ghostLayers_ << " ghost layer(s).\n"
-                         "Writing ghost layers and using sampling at the same time is not supported!" );
+                         "Writing ghost layers and using sampling at the same time is not supported!" )
    }
 
    samplingDx_ = dx;
@@ -563,7 +565,7 @@ inline shared_ptr<VTKOutput> createVTKOutput_DomainDecomposition( const shared_p
                                                                   const uint_t initialExecutionCount = 0 )
 {
    if( !sbs )
-      WALBERLA_ABORT( "creating VTK output for domain decomposition failed (StructuredBlockStorage shared pointer is NULL)" );
+      WALBERLA_ABORT( "creating VTK output for domain decomposition failed (StructuredBlockStorage shared pointer is NULL)" )
 
    return createVTKOutput_DomainDecomposition( *sbs, identifier, writeFrequency, baseFolder, executionFolder,
                                                continuousNumbering, binary, littleEndian, useMPIIO, initialExecutionCount );
@@ -578,7 +580,7 @@ inline shared_ptr<VTKOutput> createVTKOutput_DomainDecomposition( const shared_p
                                                                   const uint_t initialExecutionCount = 0 )
 {
    if( !sbs )
-      WALBERLA_ABORT( "creating VTK output for domain decomposition failed (StructuredBlockStorage shared pointer is NULL)" );
+      WALBERLA_ABORT( "creating VTK output for domain decomposition failed (StructuredBlockStorage shared pointer is NULL)" )
 
    return createVTKOutput_DomainDecomposition( sbs->getBlockStorage(), identifier, writeFrequency, baseFolder, executionFolder,
                                                continuousNumbering, binary, littleEndian, useMPIIO, initialExecutionCount );
@@ -609,7 +611,7 @@ inline shared_ptr<VTKOutput> createVTKOutput_BlockData( const shared_ptr< const
                                                         const uint_t initialExecutionCount = 0, const bool amrFileFormat = false )
 {
    if( !sbs )
-      WALBERLA_ABORT( "creating VTK output for block data failed (StructuredBlockStorage shared pointer is NULL)" );
+      WALBERLA_ABORT( "creating VTK output for block data failed (StructuredBlockStorage shared pointer is NULL)" )
 
    return createVTKOutput_BlockData( *sbs, identifier, writeFrequency, ghostLayers, forcePVTU, baseFolder, executionFolder,
                                                continuousNumbering, binary, littleEndian, useMPIIO, initialExecutionCount, amrFileFormat );
@@ -684,7 +686,7 @@ inline void writeDomainDecomposition( const shared_ptr< const BlockStorage > & b
                                       bool useMPIIO = true )
 {
    if( !bs )
-      WALBERLA_ABORT( "Writing domain decomposition failed (StructuredBlockStorage shared pointer is NULL)" );
+      WALBERLA_ABORT( "Writing domain decomposition failed (StructuredBlockStorage shared pointer is NULL)" )
 
    writeDomainDecomposition( *bs, identifier, baseFolder, executionFolder, binary, littleEndian,
                              simultaneousIOOperations, requiredStates, incompatibleStates, useMPIIO );
@@ -700,7 +702,7 @@ inline void writeDomainDecomposition( const shared_ptr< const StructuredBlockSto
                                       bool useMPIIO = true )
 {
    if( !sbs )
-      WALBERLA_ABORT( "Writing domain decomposition failed (StructuredBlockStorage shared pointer is NULL)" );
+      WALBERLA_ABORT( "Writing domain decomposition failed (StructuredBlockStorage shared pointer is NULL)" )
 
    writeDomainDecomposition( *sbs, identifier, baseFolder, executionFolder, binary, littleEndian,
                              simultaneousIOOperations, requiredStates, incompatibleStates, useMPIIO );

tests/field/CMakeLists.txt

@@ -71,6 +71,11 @@ waLBerla_generate_target_from_python(NAME CodegenJacobiCPUGeneratedJacobiKernel
 waLBerla_compile_test( FILES codegen/CodegenJacobiCPU.cpp DEPENDS gui timeloop CodegenJacobiCPUGeneratedJacobiKernel)
 waLBerla_execute_test( NAME CodegenJacobiCPU )
 
+waLBerla_generate_target_from_python(NAME CodegenADEGenerated FILE codegen/ADEKernel.py
+        OUT_FILES ADE.cpp ADE.h Dirichlet.cpp Dirichlet.h Neumann.cpp Neumann.h ADEPackInfo.cpp ADEPackInfo.h)
+waLBerla_compile_test( FILES codegen/ADE.cpp DEPENDS blockforest core field stencil timeloop vtk CodegenADEGenerated)
+waLBerla_execute_test( NAME CodegenADE )
+
 waLBerla_generate_target_from_python(NAME SweepCollectionKernel FILE codegen/SweepCollection.py
         OUT_FILES SweepCollection.h SweepCollection.cpp)
 waLBerla_compile_test( FILES codegen/SweepCollection.cpp DEPENDS timeloop SweepCollectionKernel)

tests/field/codegen/ADE.cpp

+//======================================================================================================================
+//
+//  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 ADE.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "core/Environment.h"
+#include "core/debug/TestSubsystem.h"
+
+#include "field/AddToStorage.h"
+#include "field/communication/PackInfo.h"
+#include "field/vtk/VTKWriter.h"
+
+#include "stencil/D2Q5.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+#include "vtk/VTKOutput.h"
+
+#include "ADE.h"
+#include "Dirichlet.h"
+#include "Neumann.h"
+#include "ADEPackInfo.h"
+
+
+using namespace walberla;
+
+using ScalarField = GhostLayerField<real_t, 1>;
+using VectorField = GhostLayerField<real_t, 2>;
+using flag_t      = walberla::uint8_t;
+using FlagField_T = FlagField< flag_t >;
+
+using CommScheme = blockforest::communication::UniformBufferedScheme<stencil::D2Q5>;
+
+void setupBoundary(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID,
+                   const FlagUID& fluidFlagUID, const FlagUID& dirichletZeroFlagUID, const FlagUID& dirichletOneFlagUID, const FlagUID& neumannFlagUID)
+{
+
+   const AABB & domain = sbfs->getDomain();
+   const AABB inflowAABB( 10, 25, 0, 30, 39, domain.zMax());
+
+   for (auto bIt = sbfs->begin(); bIt != sbfs->end(); ++bIt)
+   {
+
+      Block& b             = dynamic_cast< Block& >(*bIt);
+      auto flagField       = b.getData< FlagField_T >(flagFieldID);
+
+      uint8_t fluidFlag = flagField->registerFlag(fluidFlagUID);
+      uint8_t dirichletZeroFlag = flagField->registerFlag(dirichletZeroFlagUID);
+      uint8_t dirichletOneFlag = flagField->registerFlag(dirichletOneFlagUID);
+      uint8_t neumannFlag = flagField->registerFlag(neumannFlagUID);
+
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, {
+         Cell globalCell = Cell(x, y, z);
+         sbfs->transformBlockLocalToGlobalCell(globalCell, *bIt);
+         const Vector3< real_t > globalPoint = sbfs->getCellCenter(globalCell);
+
+         if(globalPoint[0] < 0 || globalPoint[1] < 0 || globalPoint[1] > domain.yMax())
+         {
+            addFlag(flagField->get(x, y, z), dirichletZeroFlag);
+         }
+         else if (globalPoint[0] > domain.xMax())
+         {
+            addFlag(flagField->get(x, y, z), neumannFlag);
+         }
+         else if (inflowAABB.contains(globalPoint))
+         {
+            addFlag(flagField->get(x, y, z), dirichletOneFlag);
+         }
+         else
+         {
+            addFlag(flagField->get(x, y, z), fluidFlag);
+         }
+      }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
+
+   }
+}
+
+int main( int argc, char ** argv )
+{
+   mpi::Environment env( argc, argv );
+   debug::enterTestMode();
+
+   logging::Logging::instance()->setLogLevel( logging::Logging::INFO );
+
+   uint_t xSize = 128;
+   uint_t ySize = 64;
+   // Create blocks
+   shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid (
+      uint_t(1) , uint_t(1),  uint_t(1),  // number of blocks in x,y,z direction
+      xSize, ySize, uint_t(1),            // how many cells per block (x,y,z)
+      real_c(1.0),                          // dx: length of one cell in physical coordinates
+      false,                                // one block per process - "false" means all blocks to one process
+      false, false, false );                 // full periodicity
+
+   real_t diffusivity = real_c(0.01);
+   real_t dx = real_c(1.0);
+   real_t dt = real_c(0.01);
+
+   BlockDataID src = field::addToStorage<ScalarField>(blocks, "src", real_c(0.0));
+   BlockDataID vel = field::addToStorage<VectorField>(blocks, "vel", real_c(0.0));
+   BlockDataID flagFieldID = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
+
+   for (auto& block : *blocks)
+   {
+      const AABB & domain = blocks->getDomain();
+      const AABB inflowAABB( 10, 25, 0, 30, 39, domain.zMax());
+
+      auto * velField = block.getData< VectorField > ( vel );
+      auto * srcField = block.getData< ScalarField > ( src );
+
+      for( auto it = velField->beginWithGhostLayer(1); it != velField->end(); ++it )
+      {
+         velField->get(it.cell(), 0) = real_c(1.0);
+         if (inflowAABB.contains( it.x(), it.y(), it.z() ))
+         {
+            srcField->get(it.cell()) = real_c(1.0);
+         }
+      }
+   }
+
+   const FlagUID fluidFlagUID("Fluid");
+   const FlagUID dirichletZeroFlagUID("dirichletZero");
+   const FlagUID dirichletOneFlagUID("dirichletOne");
+   const FlagUID neumannFlagUID("neumann");
+
+   setupBoundary(blocks, flagFieldID, fluidFlagUID, dirichletZeroFlagUID, dirichletOneFlagUID, neumannFlagUID);
+
+   pystencils::Dirichlet dirichletZero(blocks, src, real_c(0.0));
+   pystencils::Dirichlet dirichletOne(blocks, src, real_c(1.0));
+   pystencils::Neumann neumann(blocks, src);
+
+   dirichletZero.fillFromFlagField< FlagField_T >(blocks, flagFieldID, dirichletZeroFlagUID, fluidFlagUID);
+   dirichletOne.fillFromFlagField< FlagField_T >(blocks, flagFieldID, dirichletOneFlagUID, fluidFlagUID);
+   neumann.fillFromFlagField< FlagField_T >(blocks, flagFieldID, neumannFlagUID, fluidFlagUID);
+
+   CommScheme commScheme(blocks);
+   commScheme.addDataToCommunicate( make_shared<pystencils::ADEPackInfo>(src, vel) );
+
+   // Create Timeloop
+   const uint_t numberOfTimesteps = uint_t(100); // number of timesteps for non-gui runs
+   SweepTimeloop timeloop ( blocks, numberOfTimesteps );
+
+   // Registering the sweep
+   timeloop.add() << BeforeFunction(  commScheme, "Communication" )
+                  << Sweep(dirichletZero, "dirichlet zero");
+   timeloop.add() << Sweep(dirichletOne, "dirichlet one" );
+   timeloop.add() << Sweep(neumann, "neumann" );
+   timeloop.add() << Sweep( pystencils::ADE(src, vel, diffusivity, dt, dx), "advection diffusion kernel" );
+
+   const uint_t vtkWriteFrequency = 0;
+   if (vtkWriteFrequency > 0)
+   {
+      auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "vtk", vtkWriteFrequency, 0, false,
+                                                      "vtk_ADE", "simulation_step", false, true, true,
+                                                      false, 0);
+
+      auto concentrationWriter = make_shared< field::VTKWriter< ScalarField, float32 > >(src, "concentration");
+      vtkOutput->addCellDataWriter(concentrationWriter);
+
+      auto flagWriter = make_shared< field::VTKWriter< FlagField_T > >(flagFieldID, "flag");
+      vtkOutput->addCellDataWriter(flagWriter);
+
+
+      timeloop.addFuncAfterTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
+   }
+   timeloop.run();
+
+
+   return EXIT_SUCCESS;
+}

tests/field/codegen/ADEKernel.py

+import sympy as sp
+import pystencils as ps
+import pystencils_walberla as psw
+
+with psw.CodeGeneration() as ctx:
+    dim = 2
+    c, c_tmp, u = ps.fields(f"c, c_tmp, u({dim}): [{dim}d]")
+    D, dx, dt = sp.symbols("D dx dt")
+
+    pde = ps.fd.transient(c) - ps.fd.diffusion(c, D) + ps.fd.advection(c, u)
+
+    discretize = ps.fd.Discretization2ndOrder(dx, dt)
+    assign = [ps.Assignment(c_tmp.center(), discretize(pde)),]
+
+    target = ps.Target.GPU if ctx.gpu else ps.Target.CPU
+
+    psw.generate_sweep(ctx, 'ADE', assign, field_swaps=[(c, c_tmp)], target=target)
+
+    dirichlet = ps.Dirichlet(sp.Symbol("concentration"))
+    neumann = ps.Neumann()
+
+    psw.generate_boundary(ctx, 'Dirichlet', dirichlet, field=c, target=target)
+    psw.generate_boundary(ctx, 'Neumann', neumann, field=c, target=target)
+
+    psw.generate_pack_info_from_kernel(ctx, 'ADEPackInfo', assign, target=target)
\ No newline at end of file

tests/lbm_generated/CMakeLists.txt

@@ -35,6 +35,21 @@ waLBerla_compile_test( FILES InterpolationNoSlip.cpp DEPENDS InterpolationNoSlip
 # waLBerla_execute_test( NAME InterpolationNoSlip1 COMMAND $<TARGET_FILE:InterpolationNoSlip> ${CMAKE_CURRENT_SOURCE_DIR}/InterpolationNoSlip.prm -Parameters.distanceWall=0.1 )
 # waLBerla_execute_test( NAME InterpolationNoSlip2 COMMAND $<TARGET_FILE:InterpolationNoSlip> ${CMAKE_CURRENT_SOURCE_DIR}/InterpolationNoSlip.prm -Parameters.distanceWall=0.5 )
 waLBerla_execute_test( NAME InterpolationNoSlip3 COMMAND $<TARGET_FILE:InterpolationNoSlip> ${CMAKE_CURRENT_SOURCE_DIR}/InterpolationNoSlip.prm )
+
+
+if( WALBERLA_BUILD_WITH_CUDA )
+    waLBerla_generate_target_from_python(NAME RefinementCPUGPUComparisonGenerated
+            FILE RefinementCPUGPUComparison.py
+            CODEGEN_CFG refinement_comparison_codegen
+            OUT_FILES CPUVersionStorageSpecification.h CPUVersionStorageSpecification.cpp
+            CPUVersionSweepCollection.h CPUVersionSweepCollection.cpp
+            GPUVersionStorageSpecification.h GPUVersionStorageSpecification.cu
+            GPUVersionSweepCollection.h GPUVersionSweepCollection.cu
+            RefinementCPUGPUComparisonHeader.h)
+    waLBerla_compile_test( FILES RefinementCPUGPUComparison.cpp DEPENDS RefinementCPUGPUComparisonGenerated blockforest gpu field lbm_generated timeloop )
+endif()
+
+
 endif()
 
 waLBerla_generate_target_from_python(NAME FreeSlipRefinementGenerated
@@ -52,3 +67,5 @@ waLBerla_compile_test( FILES FreeSlipRefinement.cpp DEPENDS FreeSlipRefinementGe
 if( WALBERLA_DOUBLE_ACCURACY )
 waLBerla_compile_test( FILES LDC.cpp DEPENDS blockforest field lbm_generated timeloop )
 endif()
+
+

tests/lbm_generated/RefinementCPUGPUComparison.cpp

+//======================================================================================================================
+//
+//  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 RefinementCPUGPUComparison.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//! \brief Shear flow that compares GPU and CPU refinement
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/NonUniformBufferedScheme.h"
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/logging/Initialization.h"
+#include "core/math/Vector3.h"
+#include "core/timing/RemainingTimeLogger.h"
+
+#include "field/AddToStorage.h"
+#include "field/FlagField.h"
+#include "field/GhostLayerField.h"
+#include "field/vtk/VTKWriter.h"
+
+#include "geometry/InitBoundaryHandling.h"
+
+#include "gpu/GPUWrapper.h"
+#include "gpu/communication/NonUniformGPUScheme.h"
+#include "gpu/AddGPUFieldToStorage.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+#include "lbm_generated/communication/NonuniformGeneratedPdfPackInfo.h"
+#include "lbm_generated/field/AddToStorage.h"
+#include "lbm_generated/field/PdfField.h"
+#include "lbm_generated/refinement/BasicRecursiveTimeStep.h"
+#include "lbm_generated/gpu/AddToStorage.h"
+#include "lbm_generated/gpu/BasicRecursiveTimeStepGPU.h"
+#include "lbm_generated/gpu/GPUPdfField.h"
+#include "lbm_generated/gpu/NonuniformGeneratedGPUPdfPackInfo.h"
+
+// include the generated header file. It includes all generated classes
+#include "RefinementCPUGPUComparisonHeader.h"
+
+using namespace walberla;
+using namespace std::placeholders;
+
+using StorageSpecificationCPU = lbm::CPUVersionStorageSpecification;
+using Stencil_T              = StorageSpecificationCPU::Stencil;
+using CommunicationStencil_T = StorageSpecificationCPU::CommunicationStencil;
+using PdfFieldCPU             = lbm_generated::PdfField< StorageSpecificationCPU >;
+using SweepCollectionCPU = lbm::CPUVersionSweepCollection;
+
+
+using StorageSpecificationGPU = lbm::GPUVersionStorageSpecification;
+using PdfFieldCPU2             = lbm_generated::PdfField< StorageSpecificationGPU >;
+using SweepCollectionGPU = lbm::GPUVersionSweepCollection;
+using GPUPdfField_T = lbm_generated::GPUPdfField< StorageSpecificationGPU >;
+
+using VectorField_T = GhostLayerField< real_t, StorageSpecificationCPU::Stencil::D >;
+using ScalarField_T = GhostLayerField< real_t, 1 >;
+
+
+using RefinementSelectionFunctor = SetupBlockForest::RefinementSelectionFunction;
+
+struct IDs
+{
+   BlockDataID pdfField;
+   BlockDataID velocityField;
+   BlockDataID densityField;
+
+   // CPU Fields for the GPU simulation (just for result checking on the CPU, not for actual calculations)
+   BlockDataID pdfField2;
+   BlockDataID velocityField2;
+   BlockDataID densityField2;
+
+   BlockDataID pdfFieldGPU;
+   BlockDataID velocityFieldGPU;
+   BlockDataID densityFieldGPU;
+};
+
+class Refinement
+{
+ public:
+   Refinement(const uint_t depth) : refinementDepth_(depth){};
+
+   void operator()(SetupBlockForest& forest)
+   {
+      std::vector< SetupBlock* > blocks;
+      forest.getBlocks(blocks);
+      auto center = forest.getDomain().center();
+      AABB refinementAABB = AABB(center[0] - 1, center[1] - 1, center[2] - 1,
+                                 center[0] + 1, center[1] + 1, center[2] + 1);
+
+      for (auto b : blocks)
+      {
+         if (refinementAABB.intersects(b->getAABB()) && b->getLevel() < refinementDepth_)
+         {
+            b->setMarker(true);
+         }
+      }
+   }
+
+ private:
+   const uint_t refinementDepth_;
+};
+
+namespace
+{
+void createSetupBlockForest(SetupBlockForest& setupBfs, const Config::BlockHandle& domainSetup)
+{
+   Vector3< real_t > domainSize = domainSetup.getParameter< Vector3< real_t > >("domainSize");
+   Vector3< uint_t > rootBlocks = domainSetup.getParameter< Vector3< uint_t > >("rootBlocks");
+   Vector3< bool > periodic     = domainSetup.getParameter< Vector3< bool > >("periodic");
+   const uint_t refinementDepth = domainSetup.getParameter< uint_t >("refinementDepth", uint_c(1));
+
+   auto refSelection = Refinement(refinementDepth);
+   setupBfs.addRefinementSelectionFunction(std::function< void(SetupBlockForest&) >(refSelection));
+   AABB domain(real_t(0.0), real_t(0.0), real_t(0.0), domainSize[0], domainSize[1], domainSize[2]);
+   setupBfs.init(domain, rootBlocks[0], rootBlocks[1], rootBlocks[2], periodic[0], periodic[1], periodic[2]);
+   setupBfs.balanceLoad(blockforest::StaticLevelwiseCurveBalance(true), uint_c(MPIManager::instance()->numProcesses()));
+}
+
+void initShearVelocity(const shared_ptr< StructuredBlockForest >& blocks, BlockDataID velFieldID,
+                       const real_t xMagnitude = real_c(0.005), const real_t fluctuationMagnitude = real_c(0.05))
+{
+   math::seedRandomGenerator(42);
+   auto halfY = blocks->getDomainCellBB().yMax() / 2;
+   for (auto& block : *blocks)
+   {
+      Block& b           = dynamic_cast< Block& >(block);
+      const uint_t level = blocks->getLevel(b);
+
+      auto velField = block.getData< VelocityField_T >(velFieldID);
+      for (auto it = velField->beginWithGhostLayer(2); it != velField->end(); ++it)
+      {
+         Cell globalCell;
+         blocks->transformBlockLocalToGlobalCell(globalCell, block, it.cell());
+         Vector3< real_t > cellCenter = blocks->getCellCenter(globalCell, level);
+         const real_t randomReal = xMagnitude * math::realRandom< real_t >(-fluctuationMagnitude, fluctuationMagnitude);
+         velField->get(it.cell(), 1) = randomReal;
+         velField->get(it.cell(), 2) = real_c(0.0);
+
+         if (cellCenter[1] >= halfY) { velField->get(it.cell(), 0) = xMagnitude; }
+         else { velField->get(it.cell(), 0) = -xMagnitude; }
+      }
+   }
+}
+
+void initialiseAll(const shared_ptr< StructuredBlockForest >& blocks, IDs& ids, SweepCollectionCPU& sweepCollection,
+                   SweepCollectionGPU& sweepCollectionGPU, const real_t xMagnitude = real_c(0.005))
+{
+   initShearVelocity(blocks, ids.velocityField, xMagnitude);
+   for (auto& block : *blocks)
+   {
+      sweepCollection.initialise(&block, cell_idx_c(1));
+   }
+
+   initShearVelocity(blocks, ids.velocityField2, xMagnitude);
+   gpu::fieldCpy< gpu::GPUField< real_t >, VectorField_T >(blocks, ids.velocityFieldGPU, ids.velocityField2);
+   for (auto& block : *blocks)
+   {
+      sweepCollectionGPU.initialise(&block, cell_idx_c(1));
+   }
+   WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
+   WALBERLA_GPU_CHECK(gpuPeekAtLastError())
+   WALBERLA_MPI_BARRIER()
+}
+
+void checkResult(const shared_ptr< StructuredBlockForest >& blocks, IDs& ids)
+{
+   gpu::fieldCpy< PdfFieldCPU2, GPUPdfField_T >(blocks, ids.pdfField2, ids.pdfFieldGPU);
+   WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
+   WALBERLA_GPU_CHECK(gpuPeekAtLastError())
+   WALBERLA_MPI_BARRIER()
+
+   WALBERLA_LOG_INFO_ON_ROOT("Checking results ....")
+
+   for (auto& block : *blocks)
+   {
+      Block& b       = dynamic_cast< Block& >(block);
+      uint_t level   = b.getLevel();
+      auto blockID   = b.getId();
+      auto blockAABB = b.getAABB();
+
+      auto pdfs1 = b.getData< PdfFieldCPU >(ids.pdfField);
+      auto pdfs2 = b.getData< PdfFieldCPU2 >(ids.pdfField2);
+
+      auto ci = pdfs1->xyzSizeWithGhostLayer();
+      for (auto it = ci.begin(); it != ci.end(); ++it)
+      {
+         for (cell_idx_t i = 0; i < cell_idx_c(Stencil_T::Q); ++i)
+         {
+            WALBERLA_CHECK_FLOAT_EQUAL(pdfs1->get(*it, i), pdfs2->get(*it, i),
+                                       "Error occurred on block " << blockID << " with AABB: " << blockAABB
+                                                                  << " on level: " << level << " at index " << *it
+                                                                  << " on pdf: " << i)
+         }
+      }
+   }
+   WALBERLA_LOG_INFO_ON_ROOT("CPU and GPU version are equal")
+}
+
+} // close anonymous namespace
+
+
+class BoundaryCollectionCPU
+{
+ public:
+   void operator() (IBlock * /*block*/){}
+};
+
+class BoundaryCollectionGPU
+{
+ public:
+   void operator() (IBlock * /*block*/, gpuStream_t stream = nullptr){}
+};
+
+int main(int argc, char** argv)
+{
+   walberla::Environment walberlaEnv(argc, argv);
+   mpi::MPIManager::instance()->useWorldComm();
+
+   logging::configureLogging(walberlaEnv.config());
+
+   // read parameters
+   auto domainSetup = walberlaEnv.config()->getOneBlock("DomainSetup");
+   auto parameters  = walberlaEnv.config()->getOneBlock("Parameters");
+
+   const real_t omega           = parameters.getParameter< real_t >("omega");
+   const real_t xMagnitude      = parameters.getParameter< real_t >("xMagnitude");
+   const uint_t timesteps       = parameters.getParameter< uint_t >("timesteps", uint_c(10)) + uint_c(1);
+   const uint_t refinementDepth = parameters.getParameter< uint_t >("refinementDepth", uint_c(1));
+   const std::string check      = parameters.getParameter< std::string >("check");
+
+   auto loggingParameters = walberlaEnv.config()->getOneBlock("Logging");
+   bool writeSetupForestAndReturn = loggingParameters.getParameter<bool>("writeSetupForestAndReturn", false);
+
+   auto remainingTimeLoggerFrequency =
+      parameters.getParameter< real_t >("remainingTimeLoggerFrequency", real_c(3.0)); // in seconds
+
+   SetupBlockForest setupBfs;
+   WALBERLA_LOG_INFO_ON_ROOT("Generating SetupBlockForest...")
+   createSetupBlockForest(setupBfs, domainSetup);
+
+   // Create structured block forest
+   Vector3< uint_t > cellsPerBlock = domainSetup.getParameter< Vector3< uint_t > >("cellsPerBlock");
+
+   WALBERLA_LOG_INFO_ON_ROOT("Creating structured block forest...")
+   auto bfs    = std::make_shared< BlockForest >(uint_c(MPIManager::instance()->worldRank()), setupBfs);
+   auto blocks = std::make_shared< StructuredBlockForest >(bfs, cellsPerBlock[0], cellsPerBlock[1], cellsPerBlock[2]);
+   blocks->createCellBoundingBoxes();
+
+   if (writeSetupForestAndReturn)
+   {
+      WALBERLA_LOG_INFO_ON_ROOT("Writing SetupBlockForest to VTK file")
+      WALBERLA_ROOT_SECTION() { vtk::writeDomainDecomposition(blocks, "FreeSlipRefinementDomainDecomposition", "vtk_out"); }
+
+      WALBERLA_LOG_INFO_ON_ROOT("Blocks created: " << setupBfs.getNumberOfBlocks())
+      for (uint_t level = 0; level <= refinementDepth; level++)
+      {
+         WALBERLA_LOG_INFO_ON_ROOT("Level " << level << " Blocks: " << setupBfs.getNumberOfBlocks(level))
+      }
+      WALBERLA_LOG_INFO_ON_ROOT("Ending program")
+      return EXIT_SUCCESS;
+   }
+
+   WALBERLA_LOG_INFO_ON_ROOT("Blocks created: " << setupBfs.getNumberOfBlocks())
+   for (uint_t level = 0; level <= refinementDepth; level++)
+   {
+      WALBERLA_LOG_INFO_ON_ROOT("Level " << level << " Blocks: " << setupBfs.getNumberOfBlocks(level))
+   }
+
+   StorageSpecificationCPU StorageSpec = StorageSpecificationCPU();
+   IDs ids;
+
+   ids.pdfField      = lbm_generated::addPdfFieldToStorage(blocks, "pdf field", StorageSpec, uint_c(2));
+   ids.velocityField = field::addToStorage< VectorField_T >(blocks, "Velocity", real_c(0.0), field::fzyx, uint_c(2));
+   ids.densityField  = field::addToStorage< ScalarField_T >(blocks, "density", real_c(1.0), field::fzyx, uint_c(2));
+
+   StorageSpecificationGPU StorageSpecGPU = StorageSpecificationGPU();
+   ids.pdfField2      = lbm_generated::addPdfFieldToStorage(blocks, "pdf field 2", StorageSpecGPU, uint_c(2));
+   ids.velocityField2 = field::addToStorage< VectorField_T >(blocks, "Velocity 2", real_c(0.0), field::fzyx, uint_c(2));
+   ids.densityField2  = field::addToStorage< ScalarField_T >(blocks, "density 2", real_c(1.0), field::fzyx, uint_c(2));
+
+   ids.pdfFieldGPU = lbm_generated::addGPUPdfFieldToStorage< PdfFieldCPU2 >(blocks, ids.pdfField2, StorageSpecGPU, "pdfs on GPU", true);
+   ids.velocityFieldGPU = gpu::addGPUFieldToStorage< VelocityField_T >(blocks, ids.velocityField2, "velocity on GPU", true);
+   ids.densityFieldGPU = gpu::addGPUFieldToStorage< ScalarField_T >(blocks, ids.densityField2, "density on GPU", true);
+
+   SweepCollectionCPU sweepCollectionCPU(blocks, ids.pdfField, ids.densityField, ids.velocityField, omega);
+   SweepCollectionGPU sweepCollectionGPU(blocks, ids.pdfFieldGPU, ids.densityFieldGPU, ids.velocityFieldGPU, omega);
+
+   WALBERLA_LOG_INFO_ON_ROOT("Setting up communication...")
+   auto commCPU = std::make_shared< blockforest::communication::NonUniformBufferedScheme< CommunicationStencil_T > >(blocks);
+   auto packInfoCPU = lbm_generated::setupNonuniformPdfCommunication< PdfFieldCPU >(blocks, ids.pdfField);
+   commCPU->addPackInfo(packInfoCPU);
+
+   auto commGPU = std::make_shared< NonUniformGPUScheme< CommunicationStencil_T > >(blocks, false);
+   auto packInfoGPU      = lbm_generated::setupNonuniformGPUPdfCommunication< GPUPdfField_T >(blocks, ids.pdfFieldGPU);
+   commGPU->addPackInfo(packInfoGPU);
+
+   BoundaryCollectionCPU bcCPU = BoundaryCollectionCPU();
+   BoundaryCollectionGPU bcGPU = BoundaryCollectionGPU();
+
+   lbm_generated::BasicRecursiveTimeStep< PdfFieldCPU , SweepCollectionCPU , BoundaryCollectionCPU > LBMRefinementCPU( blocks, ids.pdfField, sweepCollectionCPU, bcCPU, commCPU, packInfoCPU);
+   lbm_generated::BasicRecursiveTimeStepGPU< GPUPdfField_T, SweepCollectionGPU , BoundaryCollectionGPU > LBMRefinementGPU(blocks, ids.pdfFieldGPU, sweepCollectionGPU, bcGPU, commGPU, packInfoGPU);
+
+   if(check == "coarseToFine"){
+      WALBERLA_LOG_INFO_ON_ROOT("Checking only coarse to fine communication")
+      initialiseAll(blocks, ids, sweepCollectionCPU, sweepCollectionGPU, xMagnitude);
+
+      for (auto& block : *blocks){
+         sweepCollectionCPU.streamCollide(&block);
+         sweepCollectionGPU.streamCollide(&block);
+      }
+      WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
+      WALBERLA_GPU_CHECK(gpuPeekAtLastError())
+      WALBERLA_MPI_BARRIER()
+
+      for (uint_t level = 1; level <= refinementDepth; level++){
+         commCPU->communicateCoarseToFine(level);
+      }
+
+      for (uint_t level = 1; level <= refinementDepth; level++){
+         commGPU->communicateCoarseToFine(level);
+      }
+
+      WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
+      WALBERLA_GPU_CHECK(gpuPeekAtLastError())
+      WALBERLA_MPI_BARRIER()
+
+      WALBERLA_LOG_INFO_ON_ROOT("Finished execution of the communication")
+   }
+   else if (check == "full"){
+      WALBERLA_LOG_INFO_ON_ROOT("Checking full refinement step")
+      initialiseAll(blocks, ids, sweepCollectionCPU, sweepCollectionGPU, xMagnitude);
+
+      SweepTimeloop timeloopCPU(blocks->getBlockStorage(), timesteps);
+      SweepTimeloop timeloopGPU(blocks->getBlockStorage(), timesteps);
+
+      uint_t vtkWriteFrequency = parameters.getParameter< uint_t >("vtkWriteFrequency", 0);
+      if (vtkWriteFrequency > 0)
+      {
+         auto vtkOutputCPU = vtk::createVTKOutput_BlockData(*blocks, "RefinementComparisonCPU", vtkWriteFrequency, 0, false, "vtk_out",
+                                                         "simulation_step", false, true, true, false, 0);
+
+         auto velWriterCPU = make_shared< field::VTKWriter< VectorField_T > >(ids.velocityField, "velocity");
+         vtkOutputCPU->addBeforeFunction([&]() {
+            for (auto& block : *blocks)
+            {
+               sweepCollectionCPU.calculateMacroscopicParameters(&block);
+            }
+         });
+
+         vtkOutputCPU->addCellDataWriter(velWriterCPU);
+         timeloopCPU.addFuncBeforeTimeStep(vtk::writeFiles(vtkOutputCPU), "VTK Output");
+
+
+         auto vtkOutputGPU = vtk::createVTKOutput_BlockData(*blocks, "RefinementComparisonGPU", vtkWriteFrequency, 0, false, "vtk_out",
+                                                            "simulation_step", false, true, true, false, 0);
+
+         auto velWriterGPU = make_shared< field::VTKWriter< VectorField_T > >(ids.velocityField2, "velocity");
+         vtkOutputGPU->addBeforeFunction([&]() {
+            for (auto& block : *blocks)
+            {
+               sweepCollectionGPU.calculateMacroscopicParameters(&block);
+            }
+            gpu::fieldCpy< VelocityField_T, gpu::GPUField< real_t > >(blocks, ids.velocityField2, ids.velocityFieldGPU);
+         });
+
+         vtkOutputGPU->addCellDataWriter(velWriterGPU);
+         timeloopGPU.addFuncBeforeTimeStep(vtk::writeFiles(vtkOutputGPU), "VTK Output");
+      }
+
+      LBMRefinementCPU.addRefinementToTimeLoop(timeloopCPU);
+      LBMRefinementGPU.addRefinementToTimeLoop(timeloopGPU, uint_c(0));
+
+      if (remainingTimeLoggerFrequency > 1.0)
+      {
+         timeloopCPU.addFuncAfterTimeStep(timing::RemainingTimeLogger(timeloopCPU.getNrOfTimeSteps(),
+                                                                      remainingTimeLoggerFrequency), "remaining time logger");
+
+         timeloopGPU.addFuncAfterTimeStep(timing::RemainingTimeLogger(timeloopGPU.getNrOfTimeSteps(),
+                                                                      remainingTimeLoggerFrequency), "remaining time logger");
+      }
+
+      WALBERLA_LOG_INFO_ON_ROOT("Starting CPU Simulation with " << timesteps << " timesteps")
+      timeloopCPU.run();
+      WALBERLA_MPI_BARRIER()
+      WALBERLA_LOG_INFO_ON_ROOT("Finished CPU Simulation")
+
+      WALBERLA_LOG_INFO_ON_ROOT("Starting GPU Simulation with " << timesteps << " timesteps")
+      timeloopGPU.run();
+      WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
+      WALBERLA_GPU_CHECK(gpuPeekAtLastError())
+      WALBERLA_MPI_BARRIER()
+      WALBERLA_LOG_INFO_ON_ROOT("Finished GPU Simulation")
+   }
+
+   checkResult(blocks, ids);
+   return EXIT_SUCCESS;
+}

tests/lbm_generated/RefinementCPUGPUComparison.prm

+Parameters
+{
+	omega           1.8;
+	xMagnitude   0.05;
+	timesteps       10;
+
+	remainingTimeLoggerFrequency 5; // in seconds
+	vtkWriteFrequency 0;
+
+	check full;
+}
+
+DomainSetup
+{
+   domainSize    <32, 32, 32>;
+   rootBlocks    <4, 4, 4>;
+
+   cellsPerBlock <  8, 8, 8 >;
+   periodic      <  1, 1, 1 >;
+   refinementDepth 1;
+}
+
+Logging
+{
+    logLevel info;  // info progress detail tracing
+    writeSetupForestAndReturn false;
+}
+

tests/lbm_generated/RefinementCPUGPUComparison.py

+import sympy as sp
+from pystencils import Target
+
+from pystencils.field import fields
+from pystencils.simp import insert_aliases, insert_constants
+
+from lbmpy import LBStencil, LBMConfig, LBMOptimisation
+from lbmpy.creationfunctions import create_lb_collision_rule
+from lbmpy.enums import Method, Stencil
+
+from pystencils_walberla import CodeGeneration, generate_info_header
+from lbmpy_walberla import generate_lbm_package, lbm_boundary_generator
+
+omega = sp.symbols("omega")
+inlet_velocity = sp.symbols("u_x")
+max_threads = 256
+
+with CodeGeneration() as ctx:
+    sweep_params = {'block_size': (128, 1, 1)} if ctx.gpu else {}
+
+    dtype = 'float64'
+    pdf_dtype = 'float64'
+
+    stencil = LBStencil(Stencil.D3Q19)
+    q = stencil.Q
+    dim = stencil.D
+
+    streaming_pattern = 'pull'
+
+    pdfs, pdfs_tmp = fields(f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): {pdf_dtype}[3D]", layout='fzyx')
+    velocity_field, density_field = fields(f"velocity({dim}), density(1) : {dtype}[{dim}D]", layout='fzyx')
+
+    macroscopic_fields = {'density': density_field, 'velocity': velocity_field}
+
+    method_enum = Method.SRT
+    lbm_config = LBMConfig(
+        method=method_enum,
+        stencil=stencil,
+        relaxation_rate=omega,
+        compressible=True,
+        # subgrid_scale_model=SubgridScaleModel.QR,
+        # fourth_order_correction=1e-4 if method_enum == Method.CUMULANT and stencil.Q == 27 else False,
+        field_name='pdfs',
+        streaming_pattern=streaming_pattern,
+    )
+
+    lbm_opt = LBMOptimisation(cse_global=False, cse_pdfs=False)
+
+    collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+    if lbm_config.method == Method.CUMULANT:
+        collision_rule = insert_constants(collision_rule)
+        collision_rule = insert_aliases(collision_rule)
+    lb_method = collision_rule.method
+
+    generate_lbm_package(ctx, name="CPUVersion", collision_rule=collision_rule,
+                         lbm_config=lbm_config, lbm_optimisation=lbm_opt,
+                         nonuniform=True,
+                         macroscopic_fields=macroscopic_fields,
+                         target=Target.CPU, data_type=dtype, pdfs_data_type=pdf_dtype)
+
+    generate_lbm_package(ctx, name="GPUVersion", collision_rule=collision_rule,
+                         lbm_config=lbm_config, lbm_optimisation=lbm_opt,
+                         nonuniform=True,
+                         macroscopic_fields=macroscopic_fields, gpu_indexing_params=sweep_params,
+                         target=Target.GPU, data_type=dtype, pdfs_data_type=pdf_dtype,
+                         max_threads=max_threads)
+
+    field_typedefs = {'VelocityField_T': velocity_field,
+                      'ScalarField_T': density_field}
+
+    # Info header containing correct template definitions for stencil and field
+    generate_info_header(ctx, 'RefinementCPUGPUComparisonHeader',
+                         field_typedefs=field_typedefs)