From 45719cceb070e53717f729cb5fa44e0e2c186055 Mon Sep 17 00:00:00 2001
From: markus holzer <markus.holzer@fau.de>
Date: Mon, 12 Oct 2020 20:44:37 +0200
Subject: [PATCH] Started rebasing
---
src/lbm/all.h | 1 +
src/lbm/refinement_rebase/BoundarySetup.h | 206 ++
.../EqualLevelBorderStreamCorrection.h | 171 ++
src/lbm/refinement_rebase/LinearExplosion.h | 498 +++++
src/lbm/refinement_rebase/PdfFieldPackInfo.h | 1452 ++++++++++++++
.../refinement_rebase/PdfFieldSyncPackInfo.h | 44 +
.../RefinementFunctorWrapper.h | 67 +
src/lbm/refinement_rebase/TimeStep.h | 1751 +++++++++++++++++
.../refinement_rebase/TimeStepPdfPackInfo.h | 44 +
src/lbm/refinement_rebase/TimeTracker.h | 109 +
.../VorticityBasedLevelDetermination.h | 161 ++
src/lbm/refinement_rebase/all.h | 32 +
tests/lbm/CMakeLists.txt | 3 +-
.../lbm/codegen/CodeGenerationRefinement.cpp | 63 +-
14 files changed, 4570 insertions(+), 32 deletions(-)
create mode 100644 src/lbm/refinement_rebase/BoundarySetup.h
create mode 100644 src/lbm/refinement_rebase/EqualLevelBorderStreamCorrection.h
create mode 100644 src/lbm/refinement_rebase/LinearExplosion.h
create mode 100644 src/lbm/refinement_rebase/PdfFieldPackInfo.h
create mode 100644 src/lbm/refinement_rebase/PdfFieldSyncPackInfo.h
create mode 100644 src/lbm/refinement_rebase/RefinementFunctorWrapper.h
create mode 100644 src/lbm/refinement_rebase/TimeStep.h
create mode 100644 src/lbm/refinement_rebase/TimeStepPdfPackInfo.h
create mode 100644 src/lbm/refinement_rebase/TimeTracker.h
create mode 100644 src/lbm/refinement_rebase/VorticityBasedLevelDetermination.h
create mode 100644 src/lbm/refinement_rebase/all.h
diff --git a/src/lbm/all.h b/src/lbm/all.h
index d4b0a4f15..91b01b6b0 100644
--- a/src/lbm/all.h
+++ b/src/lbm/all.h
@@ -35,5 +35,6 @@
#include "gui/all.h"
#include "lattice_model/all.h"
#include "refinement/all.h"
+#include "refinement_rebase/all.h"
#include "sweeps/all.h"
#include "vtk/all.h"
diff --git a/src/lbm/refinement_rebase/BoundarySetup.h b/src/lbm/refinement_rebase/BoundarySetup.h
new file mode 100644
index 000000000..4271b9984
--- /dev/null
+++ b/src/lbm/refinement_rebase/BoundarySetup.h
@@ -0,0 +1,206 @@
+//======================================================================================================================
+//
+// 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 LinearExplosion.h
+//! \ingroup lbm
+//! \author Florian Schornbaum <florian.schornbaum@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "blockforest/Block.h"
+#include "blockforest/BlockNeighborhoodSection.h"
+
+#include "boundary/Boundary.h"
+
+#include "core/cell/CellInterval.h"
+#include "core/debug/Debug.h"
+
+#include "domain_decomposition/StructuredBlockStorage.h"
+
+#include "field/FlagUID.h"
+
+#include "stencil/D3Q27.h"
+
+
+
+namespace walberla {
+namespace lbm {
+namespace refinement {
+
+
+namespace internal {
+inline shared_ptr< BoundaryConfiguration > defaultBoundaryConfiguration( const Cell &, const Vector3< real_t > & )
+{
+ return shared_ptr< BoundaryConfiguration >( new BoundaryConfiguration );
+}
+}
+
+
+//**********************************************************************************************************************
+/*!
+* \brief Function for consistently setting boundaries across block borders (even if there is a level transition)
+*
+* Boundaries are set not just in the inside of the block that is passed to this function, but also in the block's
+* ghost layers.
+* The user must provide a callback function "isBoundary": For any point in the 3D simulation space, this function must
+* either return true or false. Returning true means at this point in space, a boundary of type "flag" must be set.
+* If the boundary requires a corresponding boundary configuration parameter, a second callback function "parameter"
+* must be provided: If a boundary is set at a certain cell, this cell is passed to the callback function together
+* with the point in the 3D simulation space that, prior to that, was passed to "isBoundary". The function then must
+* return a shared pointer to a BoundaryConfiguration object.
+*
+* \param blocks block storage data structure
+* \param block block for which boundaries of type "flag" will be set
+* \param boundaryHandlingId block data ID for the boundary handling object
+* \param flag flag UID that corresponds to the boundary which will be set
+* \param isBoundary callback function for determining if for a certain point in the 3D simulation space a
+* boundary corresponding to "flag" must be set
+* \param parameter callback function for retrieving boundary configuration parameters
+*/
+//**********************************************************************************************************************
+template< typename BoundaryHandling_T, bool ForceBoundary >
+void consistentlySetBoundary( StructuredBlockStorage & blocks, Block & block, const BlockDataID & boundaryHandlingId, const FlagUID & flag,
+ const std::function< bool ( const Vector3< real_t > & ) > & isBoundary,
+ const std::function< shared_ptr< BoundaryConfiguration > ( const Cell &, const Vector3< real_t > & ) > & parameter = internal::defaultBoundaryConfiguration );
+
+
+
+/// For documentation see function "template< typename BoundaryHandling_T, bool ForceBoundary > void consistentlySetBoundary(...)"
+template< typename BoundaryHandling_T >
+void consistentlySetBoundary( StructuredBlockStorage & blocks, Block & block, const BlockDataID & boundaryHandlingId, const FlagUID & flag,
+ const std::function< bool ( const Vector3< real_t > & ) > & isBoundary,
+ const std::function< shared_ptr< BoundaryConfiguration > ( const Cell &, const Vector3< real_t > & ) > & parameter = internal::defaultBoundaryConfiguration )
+{
+ consistentlySetBoundary< BoundaryHandling_T, false >( blocks, block, boundaryHandlingId, flag, isBoundary, parameter );
+}
+
+
+
+/// For documentation see function "template< typename BoundaryHandling_T, bool ForceBoundary > void consistentlySetBoundary(...)"
+template< typename BoundaryHandling_T >
+void consistentlyForceBoundary( StructuredBlockStorage & blocks, Block & block, const BlockDataID & boundaryHandlingId, const FlagUID & flag,
+ const std::function< bool ( const Vector3< real_t > & ) > & isBoundary,
+ const std::function< shared_ptr< BoundaryConfiguration > ( const Cell &, const Vector3< real_t > & ) > & parameter = internal::defaultBoundaryConfiguration )
+{
+ consistentlySetBoundary< BoundaryHandling_T, true >( blocks, block, boundaryHandlingId, flag, isBoundary, parameter );
+}
+
+
+
+template< typename BoundaryHandling_T, bool ForceBoundary >
+void consistentlySetBoundary( StructuredBlockStorage & blocks, Block & block, const BlockDataID & boundaryHandlingId, const FlagUID & flag,
+ const std::function< bool ( const Vector3< real_t > & ) > & isBoundary,
+ const std::function< shared_ptr< BoundaryConfiguration > ( const Cell &, const Vector3< real_t > & ) > & parameter )
+{
+ const uint_t level = blocks.getLevel(block);
+
+ BoundaryHandling_T * boundaryHandling = block.template getData< BoundaryHandling_T >( boundaryHandlingId );
+
+ const auto * const flagField = boundaryHandling->getFlagField();
+ int ghostLayers = int_c( flagField->nrOfGhostLayers() );
+
+ CellInterval cells = flagField->xyzSize();
+ cells.expand( cell_idx_c(ghostLayers) );
+
+ std::vector< CellInterval > coarseRegions;
+ for( auto dir = stencil::D3Q27::beginNoCenter(); dir != stencil::D3Q27::end(); ++dir )
+ {
+ const auto index = blockforest::getBlockNeighborhoodSectionIndex( dir.cx(), dir.cy(), dir.cz() );
+ if( block.neighborhoodSectionHasLargerBlock( index ) )
+ {
+ CellInterval coarseRegion( cells );
+ for( uint_t i = 0; i != 3; ++i )
+ {
+ const auto c = stencil::c[i][*dir];
+
+ if( c == -1 ) coarseRegion.max()[i] = coarseRegion.min()[i] + cell_idx_c( 2 * ghostLayers - 1 );
+ else if( c == 1 ) coarseRegion.min()[i] = coarseRegion.max()[i] - cell_idx_c( 2 * ghostLayers - 1 );
+ }
+ coarseRegions.push_back( coarseRegion );
+ }
+ }
+
+ for( auto cell = cells.begin(); cell != cells.end(); ++cell )
+ {
+ bool inCoarseRegion( false );
+ for( auto region = coarseRegions.begin(); region != coarseRegions.end() && !inCoarseRegion; ++region )
+ inCoarseRegion = region->contains( *cell );
+
+ if( !inCoarseRegion )
+ {
+ Vector3< real_t > center;
+ blocks.getBlockLocalCellCenter( block, *cell, center );
+ blocks.mapToPeriodicDomain( center );
+
+ if( isBoundary(center) )
+ {
+ if( ForceBoundary )
+ {
+ auto p = parameter( *cell, center );
+ boundaryHandling->forceBoundary( flag, cell->x(), cell->y(), cell->z(), *p );
+ }
+ else
+ {
+ auto p = parameter( *cell, center );
+ boundaryHandling->setBoundary( flag, cell->x(), cell->y(), cell->z(), *p );
+ }
+ }
+ }
+ else
+ {
+ Cell globalCell( *cell );
+ blocks.transformBlockLocalToGlobalCell( globalCell, block );
+
+ Cell coarseCell( globalCell );
+ for( uint_t i = 0; i < 3; ++i )
+ {
+ if( coarseCell[i] < cell_idx_t(0) )
+ {
+ coarseCell[i] = -( ( cell_idx_t(1) - coarseCell[i] ) >> 1 );
+ }
+ else
+ {
+ coarseCell[i] >>= 1;
+ }
+ }
+
+ Vector3< real_t > coarseCenter;
+ blocks.getCellCenter( coarseCenter, coarseCell, level - uint_t(1) );
+ blocks.mapToPeriodicDomain( coarseCenter );
+
+ if( isBoundary(coarseCenter) )
+ {
+ if( ForceBoundary )
+ {
+ auto p = parameter( *cell, coarseCenter );
+ boundaryHandling->forceBoundary( flag, cell->x(), cell->y(), cell->z(), *p );
+ }
+ else
+ {
+ auto p = parameter( *cell, coarseCenter );
+ boundaryHandling->setBoundary( flag, cell->x(), cell->y(), cell->z(), *p );
+ }
+ }
+ }
+ }
+}
+
+
+
+} // namespace refinement
+} // namespace lbm
+} // namespace walberla
diff --git a/src/lbm/refinement_rebase/EqualLevelBorderStreamCorrection.h b/src/lbm/refinement_rebase/EqualLevelBorderStreamCorrection.h
new file mode 100644
index 000000000..2cc428b98
--- /dev/null
+++ b/src/lbm/refinement_rebase/EqualLevelBorderStreamCorrection.h
@@ -0,0 +1,171 @@
+//======================================================================================================================
+//
+// 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 EqualLevelBorderStreamCorrection.h
+//! \ingroup lbm
+//! \author Florian Schornbaum <florian.schornbaum@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "lbm/field/PdfField.h"
+
+#include "blockforest/Block.h"
+#include "blockforest/BlockNeighborhoodSection.h"
+
+#include "core/cell/CellInterval.h"
+
+#include "stencil/D2CornerStencil.h"
+#include "stencil/D3EdgeCornerStencil.h"
+
+#include <type_traits>
+
+
+
+namespace walberla {
+namespace lbm {
+namespace refinement {
+
+
+
+namespace internal {
+template< typename LatticeModel_T, class Enable = void >
+struct EdgeCornerStencil;
+template< typename LatticeModel_T >
+struct EdgeCornerStencil< LatticeModel_T, typename std::enable_if< LatticeModel_T::Stencil::D == 2 >::type >
+{
+ typedef stencil::D2CornerStencil type;
+};
+template< typename LatticeModel_T >
+struct EdgeCornerStencil< LatticeModel_T, typename std::enable_if< LatticeModel_T::Stencil::D == 3 >::type >
+{
+ typedef stencil::D3EdgeCornerStencil type;
+};
+}
+
+
+
+template< typename LatticeModel_T >
+class EqualLevelBorderStreamCorrection
+{
+public:
+
+ typedef PdfField< LatticeModel_T > PdfField_T;
+ typedef typename internal::EdgeCornerStencil< LatticeModel_T >::type EdgeCornerStencil_T;
+
+ EqualLevelBorderStreamCorrection( const BlockDataID & pdfFieldId ) : pdfFieldId_( pdfFieldId ) {}
+
+ void prepare( Block * block ) const { helper( block, true ); }
+ void apply( Block * block ) const { helper( block, false ); }
+
+protected:
+
+ bool finerNeighborExistsInVicinity( Block * block, stencil::Direction dir ) const;
+ void helper( Block * block, const bool prep ) const;
+
+ BlockDataID pdfFieldId_;
+};
+
+
+
+template< typename LatticeModel_T >
+bool EqualLevelBorderStreamCorrection< LatticeModel_T >::finerNeighborExistsInVicinity( Block * block, stencil::Direction dir ) const
+{
+ Vector3<int> min( -1 );
+ Vector3<int> max( 1 );
+
+ for( uint_t i = 0; i != 3; ++i )
+ {
+ if( stencil::c[i][dir] == -1 ) max[i] = 0;
+ if( stencil::c[i][dir] == 1 ) min[i] = 0;
+ }
+ if( LatticeModel_T::Stencil::D == uint_t(2) )
+ min[2] = max[2] = 0;
+
+ for( int z = min[2]; z <= max[2]; ++z ) {
+ for( int y = min[1]; y <= max[1]; ++y ) {
+ for( int x = min[0]; x <= max[0]; ++x )
+ {
+ if( x == 0 && y == 0 && z == 0 )
+ continue;
+ if( block->neighborhoodSectionHasSmallerBlocks( blockforest::getBlockNeighborhoodSectionIndex(x,y,z) ) )
+ return true;
+ }
+ }
+ }
+
+ return false;
+}
+
+
+
+template< typename LatticeModel_T >
+void EqualLevelBorderStreamCorrection< LatticeModel_T >::helper( Block * block, const bool prep ) const
+{
+ PdfField_T * pdfField = block->template getData< PdfField_T >( pdfFieldId_ );
+ const CellInterval xyzSize = pdfField->xyzSize();
+
+ for( auto dir = EdgeCornerStencil_T::beginNoCenter(); dir != EdgeCornerStencil_T::end(); ++dir )
+ {
+ if( block->neighborhoodSectionHasEquallySizedBlock( blockforest::getBlockNeighborhoodSectionIndex(*dir) ) &&
+ finerNeighborExistsInVicinity( block, *dir ) )
+ {
+ CellInterval region( xyzSize );
+ region.expand( cell_idx_t(1) );
+ for( uint_t i = 0; i < 3; ++i )
+ {
+ if( stencil::c[i][*dir] == -1 ) region.max()[i] = region.min()[i];
+ else if( stencil::c[i][*dir] == 1 ) region.min()[i] = region.max()[i];
+ else
+ {
+ region.min()[i] += cell_idx_t(1);
+ region.max()[i] -= cell_idx_t(1);
+ }
+ }
+
+ auto invDir = stencil::inverseDir[ *dir ];
+
+ if( prep )
+ {
+ for( auto cell = pdfField->beginSliceXYZ( region ); cell != pdfField->end(); ++cell ) {
+ for( uint_t n = 0; n < PdfField_T::Stencil::d_per_d_length[invDir]; ++n )
+ {
+ auto d = PdfField_T::Stencil::d_per_d[invDir][n];
+ auto idx = PdfField_T::Stencil::idx[ d ];
+ cell[idx] = cell.neighbor( d, idx );
+ }
+ }
+ }
+ else
+ {
+ for( auto cell = pdfField->beginSliceXYZ( region ); cell != pdfField->end(); ++cell ) {
+ for( uint_t n = 0; n < PdfField_T::Stencil::d_per_d_length[invDir]; ++n )
+ {
+ auto d = PdfField_T::Stencil::d_per_d[invDir][n];
+ auto idx = PdfField_T::Stencil::idx[ d ];
+ cell.neighbor( d, idx ) = cell[idx];
+ }
+ }
+ }
+ }
+ }
+}
+
+
+
+} // namespace refinement
+} // namespace lbm
+} // namespace walberla
diff --git a/src/lbm/refinement_rebase/LinearExplosion.h b/src/lbm/refinement_rebase/LinearExplosion.h
new file mode 100644
index 000000000..c70b7befe
--- /dev/null
+++ b/src/lbm/refinement_rebase/LinearExplosion.h
@@ -0,0 +1,498 @@
+//======================================================================================================================
+//
+// 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 LinearExplosion.h
+//! \ingroup lbm
+//! \author Florian Schornbaum <florian.schornbaum@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "lbm/field/PdfField.h"
+#include "lbm/lattice_model/NeighborsStencil.h"
+
+#include "blockforest/Block.h"
+#include "blockforest/BlockNeighborhoodSection.h"
+
+#include "core/OpenMP.h"
+#include "core/cell/CellInterval.h"
+#include "core/debug/Debug.h"
+
+#include "stencil/D3Q27.h"
+
+#include <map>
+
+
+
+namespace walberla {
+namespace lbm {
+namespace refinement {
+
+
+
+namespace internal {
+
+// These two functions have been (static) member functions of class LinearExplosion. Sadly, there have been strange compile errors with
+// the IBM compiler -> possible bug in the compiler. As a workaround, these two functions are now two free functions ...
+
+template< typename PdfField_T, typename BoundaryHandling_T, typename CoarseField >
+void fillTemporaryCoarseField( const cell_idx_t y, const cell_idx_t z, const CellInterval & coarse,
+ const PdfField_T * const pdfField, const BoundaryHandling_T * const boundaryHandling,
+ const shared_ptr< CoarseField > & tmpField );
+
+template< typename PdfField_T, typename BoundaryHandling_T, typename CoarseField, typename BoolField >
+void linearInterpolation( const cell_idx_t y, const cell_idx_t z, const CellInterval & fine,
+ PdfField_T * const pdfField, const BoundaryHandling_T * const boundaryHandling,
+ const shared_ptr< CoarseField > & tmpField, const shared_ptr< BoolField > & boolField );
+
+} // namespace internal
+
+
+
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+class LinearExplosion
+{
+private:
+
+ typedef typename LatticeModel_T::Stencil Stencil;
+
+ typedef PdfField< LatticeModel_T > PdfField_T;
+
+ typedef field::GhostLayerField< typename PdfField_T::value_type, PdfField_T::F_SIZE > CoarseField;
+ typedef field::GhostLayerField< bool, 1 > BoolField;
+
+public:
+
+ LinearExplosion( const BlockDataID & pdfFieldId, const ConstBlockDataID & boundaryHandlingId ) :
+ pdfFieldId_( pdfFieldId ), boundaryHandlingId_( boundaryHandlingId ) {}
+
+ void operator()( Block * block );
+
+protected:
+
+
+
+ BlockDataID pdfFieldId_;
+ ConstBlockDataID boundaryHandlingId_;
+
+ std::map< Cell, std::pair< shared_ptr< CoarseField >, shared_ptr< BoolField > > > tmpFields_;
+};
+
+
+
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+void LinearExplosion< LatticeModel_T, BoundaryHandling_T >::operator()( Block * block )
+{
+ PdfField_T * pdfField = block->template getData< PdfField_T >( pdfFieldId_ );
+ const BoundaryHandling_T * boundaryHandling = block->template getData< BoundaryHandling_T >( boundaryHandlingId_ );
+
+ CellInterval fineInterval = pdfField->xyzSize();
+ CellInterval coarseInterval = fineInterval;
+
+ fineInterval.expand( cell_idx_t(4) );
+ coarseInterval.xMax() = ( coarseInterval.xMax() >> 1 );
+ coarseInterval.yMax() = ( coarseInterval.yMax() >> 1 );
+ coarseInterval.zMax() = ( Stencil::D == uint_t(3) ) ? ( coarseInterval.zMax() >> 1 ) : cell_idx_t(0);
+ coarseInterval.expand( cell_idx_t(2) );
+
+ std::vector< CellInterval > coarseIntervals;
+ std::vector< CellInterval > fineIntervals;
+
+ // determine regions that need to be interpolated
+
+ for( auto d = NeighborsStencil< LatticeModel_T >::type::beginNoCenter(); d != NeighborsStencil< LatticeModel_T >::type::end(); ++d )
+ {
+ if( block->neighborhoodSectionHasLargerBlock( blockforest::getBlockNeighborhoodSectionIndex(*d) ) )
+ {
+ CellInterval coarse = coarseInterval;
+ CellInterval fine = fineInterval;
+
+ for( uint_t i = 0; i < 3; ++i )
+ {
+ if( stencil::c[i][*d] == -1 )
+ {
+ coarse.max()[i] = coarse.min()[i] + cell_idx_t(1);
+ fine.max()[i] = fine.min()[i] + cell_idx_t(3);
+ }
+ else if( stencil::c[i][*d] == 1 )
+ {
+ coarse.min()[i] = coarse.max()[i] - cell_idx_t(1);
+ fine.min()[i] = fine.max()[i] - cell_idx_t(3);
+ }
+ else
+ {
+ WALBERLA_ASSERT_EQUAL( stencil::c[i][*d], 0 );
+ coarse.min()[i] += cell_idx_t(2);
+ coarse.max()[i] -= cell_idx_t(2);
+ fine.min()[i] += cell_idx_t(4);
+ fine.max()[i] -= cell_idx_t(4);
+ }
+ }
+
+ coarseIntervals.push_back( coarse );
+ fineIntervals.push_back( fine );
+ }
+ }
+
+ // it there are no coarse neighbors -> return
+
+ if( coarseIntervals.empty() )
+ {
+ WALBERLA_ASSERT( fineIntervals.empty() );
+ return;
+ }
+
+ // fetch temporary fields from temporary field storage "tmpFields_"
+
+ if( tmpFields_.find( coarseInterval.max() ) == tmpFields_.end() )
+ {
+ tmpFields_[ coarseInterval.max() ] = std::make_pair( make_shared< CoarseField >( coarseInterval.xSize() - uint_t(4), coarseInterval.ySize() - uint_t(4),
+ coarseInterval.zSize() - uint_t(4), uint_t(2), pdfField->layout() ),
+ make_shared< BoolField >( coarseInterval.xSize() - uint_t(4), coarseInterval.ySize() - uint_t(4),
+ coarseInterval.zSize() - uint_t(4), uint_t(3), pdfField->layout() ) );
+ }
+ auto fields = tmpFields_[ coarseInterval.max() ];
+ auto tmpField = fields.first;
+ auto boolField = fields.second;
+
+#ifndef NDEBUG
+ std::fill( tmpField->beginWithGhostLayer(), tmpField->end(), std::numeric_limits< typename PdfField_T::value_type >::quiet_NaN() );
+ std::fill( boolField->beginWithGhostLayer(), boolField->end(), false );
+#endif
+
+ // fill temporary coarse field with values (only regions that are required for the interpolation!)
+
+ for( auto coarse = coarseIntervals.begin(); coarse != coarseIntervals.end(); ++coarse )
+ {
+#ifdef _OPENMP
+
+ if( coarse->zSize() >= coarse->ySize() )
+ {
+ int zMin = int_c( coarse->zMin() );
+ int zMax = int_c( coarse->zMax() ) + 1;
+
+ #pragma omp parallel for schedule(static)
+ for( int zi = zMin; zi < zMax; ++zi )
+ {
+ const cell_idx_t z = cell_idx_c(zi);
+ for( cell_idx_t y = coarse->yMin(); y <= coarse->yMax(); ++y )
+ internal::fillTemporaryCoarseField< PdfField_T, BoundaryHandling_T, CoarseField >( y, z, *coarse, pdfField, boundaryHandling, tmpField );
+ }
+ }
+ else
+ {
+ int yMin = int_c( coarse->yMin() );
+ int yMax = int_c( coarse->yMax() ) + 1;
+
+ #pragma omp parallel for schedule(static)
+ for( int yi = yMin; yi < yMax; ++yi )
+ {
+ const cell_idx_t y = cell_idx_c(yi);
+ for( cell_idx_t z = coarse->zMin(); z <= coarse->zMax(); ++z )
+ internal::fillTemporaryCoarseField< PdfField_T, BoundaryHandling_T, CoarseField >( y, z, *coarse, pdfField, boundaryHandling, tmpField );
+ }
+ }
+#else
+ for( cell_idx_t z = coarse->zMin(); z <= coarse->zMax(); ++z ) {
+ for( cell_idx_t y = coarse->yMin(); y <= coarse->yMax(); ++y )
+ {
+ internal::fillTemporaryCoarseField< PdfField_T, BoundaryHandling_T, CoarseField >( y, z, *coarse, pdfField, boundaryHandling, tmpField );
+ }
+ }
+#endif
+
+ CellInterval expanded( *coarse );
+ expanded.expand( cell_idx_t(1) );
+
+ for( auto cell = boolField->beginSliceXYZ( expanded ); cell != boolField->end(); ++cell )
+ *cell = false;
+ }
+
+ for( auto coarse = coarseIntervals.begin(); coarse != coarseIntervals.end(); ++coarse ) {
+ for( auto z = coarse->zMin(); z <= coarse->zMax(); ++z ) {
+ for( auto y = coarse->yMin(); y <= coarse->yMax(); ++y ) {
+ for( auto x = coarse->xMin(); x <= coarse->xMax(); ++x )
+ {
+ if( boundaryHandling->isDomain( cell_idx_t(2) * x, cell_idx_t(2) * y, cell_idx_t(2) * z ) )
+ {
+ WALBERLA_ASSERT( boundaryHandling->isDomain( cell_idx_t(2) * x + cell_idx_t(1), cell_idx_t(2) * y, cell_idx_t(2) * z ) );
+ WALBERLA_ASSERT( boundaryHandling->isDomain( cell_idx_t(2) * x, cell_idx_t(2) * y + cell_idx_t(1), cell_idx_t(2) * z ) );
+ WALBERLA_ASSERT( boundaryHandling->isDomain( cell_idx_t(2) * x + cell_idx_t(1), cell_idx_t(2) * y + cell_idx_t(1), cell_idx_t(2) * z ) );
+ if( Stencil::D == uint_t(3) )
+ {
+ WALBERLA_ASSERT( boundaryHandling->isDomain( cell_idx_t(2) * x, cell_idx_t(2) * y, cell_idx_t(2) * z + cell_idx_t(1) ) );
+ WALBERLA_ASSERT( boundaryHandling->isDomain( cell_idx_t(2) * x + cell_idx_t(1), cell_idx_t(2) * y, cell_idx_t(2) * z + cell_idx_t(1) ) );
+ WALBERLA_ASSERT( boundaryHandling->isDomain( cell_idx_t(2) * x, cell_idx_t(2) * y + cell_idx_t(1), cell_idx_t(2) * z + cell_idx_t(1) ) );
+ WALBERLA_ASSERT( boundaryHandling->isDomain( cell_idx_t(2) * x + cell_idx_t(1), cell_idx_t(2) * y + cell_idx_t(1), cell_idx_t(2) * z + cell_idx_t(1) ) );
+ }
+
+ boolField->get(x,y,z) = true;
+ }
+ else
+ {
+ WALBERLA_ASSERT( !boundaryHandling->isDomain( cell_idx_t(2) * x + cell_idx_t(1), cell_idx_t(2) * y, cell_idx_t(2) * z ) );
+ WALBERLA_ASSERT( !boundaryHandling->isDomain( cell_idx_t(2) * x, cell_idx_t(2) * y + cell_idx_t(1), cell_idx_t(2) * z ) );
+ WALBERLA_ASSERT( !boundaryHandling->isDomain( cell_idx_t(2) * x + cell_idx_t(1), cell_idx_t(2) * y + cell_idx_t(1), cell_idx_t(2) * z ) );
+ if( Stencil::D == uint_t(3) )
+ {
+ WALBERLA_ASSERT( !boundaryHandling->isDomain( cell_idx_t(2) * x, cell_idx_t(2) * y, cell_idx_t(2) * z + cell_idx_t(1) ) );
+ WALBERLA_ASSERT( !boundaryHandling->isDomain( cell_idx_t(2) * x + cell_idx_t(1), cell_idx_t(2) * y, cell_idx_t(2) * z + cell_idx_t(1) ) );
+ WALBERLA_ASSERT( !boundaryHandling->isDomain( cell_idx_t(2) * x, cell_idx_t(2) * y + cell_idx_t(1), cell_idx_t(2) * z + cell_idx_t(1) ) );
+ WALBERLA_ASSERT( !boundaryHandling->isDomain( cell_idx_t(2) * x + cell_idx_t(1), cell_idx_t(2) * y + cell_idx_t(1), cell_idx_t(2) * z + cell_idx_t(1) ) );
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // linear interpolation
+
+ for( auto fine = fineIntervals.begin(); fine != fineIntervals.end(); ++fine )
+ {
+ WALBERLA_ASSERT( (fine->xSize() & uint_t(1)) == uint_t(0) );
+ WALBERLA_ASSERT( (fine->ySize() & uint_t(1)) == uint_t(0) );
+ WALBERLA_ASSERT( ( Stencil::D == uint_t(2) && fine->zSize() == uint_t(1) ) || ( Stencil::D == uint_t(3) && (fine->zSize() & uint_t(1)) == uint_t(0) ) );
+
+#ifdef _OPENMP
+
+ if( fine->zSize() >= fine->ySize() )
+ {
+ int zSize = (Stencil::D == uint_t(3)) ? ( int_c( fine->zSize() ) / 2 ) : 1;
+
+ #pragma omp parallel for schedule(static)
+ for( int zi = 0; zi < zSize; ++zi )
+ {
+ const cell_idx_t z = fine->zMin() + cell_idx_c(zi) * cell_idx_t(2);
+ for( cell_idx_t y = fine->yMin(); y <= fine->yMax(); y += cell_idx_t(2) )
+ internal::linearInterpolation< PdfField_T, BoundaryHandling_T, CoarseField, BoolField >( y, z, *fine, pdfField, boundaryHandling, tmpField, boolField );
+ }
+ }
+ else
+ {
+ int ySize = int_c( fine->ySize() ) / 2;
+
+ #pragma omp parallel for schedule(static)
+ for( int yi = 0; yi < ySize; ++yi )
+ {
+ const cell_idx_t y = fine->yMin() + cell_idx_c(yi) * cell_idx_t(2);
+ for( cell_idx_t z = fine->zMin(); z <= fine->zMax(); z += cell_idx_t(2) )
+ internal::linearInterpolation< PdfField_T, BoundaryHandling_T, CoarseField, BoolField >( y, z, *fine, pdfField, boundaryHandling, tmpField, boolField );
+ }
+ }
+#else
+ for( cell_idx_t z = fine->zMin(); z <= fine->zMax(); z += cell_idx_t(2) ) {
+ for( cell_idx_t y = fine->yMin(); y <= fine->yMax(); y += cell_idx_t(2) )
+ {
+ internal::linearInterpolation< PdfField_T, BoundaryHandling_T, CoarseField, BoolField >( y, z, *fine, pdfField, boundaryHandling, tmpField, boolField );
+ }
+ }
+#endif
+ }
+}
+
+
+
+namespace internal {
+
+template< typename PdfField_T, typename BoundaryHandling_T, typename CoarseField >
+void fillTemporaryCoarseField( const cell_idx_t y, const cell_idx_t z, const CellInterval & coarse,
+ const PdfField_T * const pdfField, const BoundaryHandling_T * const boundaryHandling,
+ const shared_ptr< CoarseField > & tmpField )
+{
+ for( cell_idx_t x = coarse.xMin(); x <= coarse.xMax(); ++x ) {
+ for( uint_t f = 0; f < PdfField_T::F_SIZE; ++f )
+ {
+ const auto fx = cell_idx_t(2) * x;
+ const auto fy = cell_idx_t(2) * y;
+ const auto fz = cell_idx_t(2) * z;
+
+ const auto value = pdfField->get( fx, fy, fz, f );
+
+ if( boundaryHandling->isDomain(fx,fy,fz) )
+ {
+ WALBERLA_ASSERT( !math::isnan( value ) );
+ }
+
+ tmpField->get(x,y,z,f) = value;
+
+ /*
+ WALBERLA_ASSERT( !math::isnan( pdfField->get( fx , fy , fz , f ) ) );
+ WALBERLA_ASSERT( !math::isnan( pdfField->get( fx , fy , fz + cell_idx_t(1), f ) ) );
+ WALBERLA_ASSERT( !math::isnan( pdfField->get( fx , fy + cell_idx_t(1), fz , f ) ) );
+ WALBERLA_ASSERT( !math::isnan( pdfField->get( fx , fy + cell_idx_t(1), fz + cell_idx_t(1), f ) ) );
+ WALBERLA_ASSERT( !math::isnan( pdfField->get( fx + cell_idx_t(1), fy , fz , f ) ) );
+ WALBERLA_ASSERT( !math::isnan( pdfField->get( fx + cell_idx_t(1), fy , fz + cell_idx_t(1), f ) ) );
+ WALBERLA_ASSERT( !math::isnan( pdfField->get( fx + cell_idx_t(1), fy + cell_idx_t(1), fz , f ) ) );
+ WALBERLA_ASSERT( !math::isnan( pdfField->get( fx + cell_idx_t(1), fy + cell_idx_t(1), fz + cell_idx_t(1), f ) ) );
+
+ auto value = pdfField->get( fx , fy , fz , f );
+ value += pdfField->get( fx , fy , fz + cell_idx_t(1), f );
+ value += pdfField->get( fx , fy + cell_idx_t(1), fz , f );
+ value += pdfField->get( fx , fy + cell_idx_t(1), fz + cell_idx_t(1), f );
+ value += pdfField->get( fx + cell_idx_t(1), fy , fz , f );
+ value += pdfField->get( fx + cell_idx_t(1), fy , fz + cell_idx_t(1), f );
+ value += pdfField->get( fx + cell_idx_t(1), fy + cell_idx_t(1), fz , f );
+ value += pdfField->get( fx + cell_idx_t(1), fy + cell_idx_t(1), fz + cell_idx_t(1), f );
+
+ tmpField->get(x,y,z,f) = real_t(0.125) * value;
+ */
+ }
+ }
+}
+
+
+
+template< typename PdfField_T, typename BoundaryHandling_T, typename CoarseField, typename BoolField >
+void linearInterpolation( const cell_idx_t y, const cell_idx_t z, const CellInterval & fine,
+ PdfField_T * const pdfField, const BoundaryHandling_T * const boundaryHandling,
+ const shared_ptr< CoarseField > & tmpField, const shared_ptr< BoolField > & boolField )
+{
+ // see: Chen et al., "Grid refinement in lattice Boltzmann methods based on volumetric formulation", 2006
+
+ static const real_t weights[8][3] = { { -0.25, -0.25, -0.25 },
+ { -0.25, -0.25, 0.25 },
+ { -0.25, 0.25, -0.25 },
+ { -0.25, 0.25, 0.25 },
+ { 0.25, -0.25, -0.25 },
+ { 0.25, -0.25, 0.25 },
+ { 0.25, 0.25, -0.25 },
+ { 0.25, 0.25, 0.25 } };
+
+ for( cell_idx_t x = fine.xMin(); x < fine.xMax(); x += cell_idx_t(2) )
+ {
+ if( boundaryHandling->isDomain(x,y,z) )
+ {
+ WALBERLA_ASSERT( boundaryHandling->isDomain( x + cell_idx_t(1), y, z ) );
+ WALBERLA_ASSERT( boundaryHandling->isDomain( x, y + cell_idx_t(1), z ) );
+ WALBERLA_ASSERT( boundaryHandling->isDomain( x + cell_idx_t(1), y + cell_idx_t(1), z ) );
+ if( PdfField_T::Stencil::D == uint_t(3) )
+ {
+ WALBERLA_ASSERT( boundaryHandling->isDomain( x, y, z + cell_idx_t(1) ) );
+ WALBERLA_ASSERT( boundaryHandling->isDomain( x + cell_idx_t(1), y, z + cell_idx_t(1) ) );
+ WALBERLA_ASSERT( boundaryHandling->isDomain( x, y + cell_idx_t(1), z + cell_idx_t(1) ) );
+ WALBERLA_ASSERT( boundaryHandling->isDomain( x + cell_idx_t(1), y + cell_idx_t(1), z + cell_idx_t(1) ) );
+ }
+
+ Cell cell( ( ( x + cell_idx_t(4) ) >> 1 ) - cell_idx_t(2),
+ ( ( y + cell_idx_t(4) ) >> 1 ) - cell_idx_t(2),
+ (PdfField_T::Stencil::D == uint_t(3)) ? ( ( ( z + cell_idx_t(4) ) >> 1 ) - cell_idx_t(2) ) : cell_idx_t(0) );
+
+ Cell min[3], max[3];
+
+ min[0][0] = cell[0] - cell_idx_t(1); min[0][1] = cell[1]; min[0][2] = cell[2];
+ max[0][0] = cell[0] + cell_idx_t(1); max[0][1] = cell[1]; max[0][2] = cell[2];
+ min[1][0] = cell[0]; min[1][1] = cell[1] - cell_idx_t(1); min[1][2] = cell[2];
+ max[1][0] = cell[0]; max[1][1] = cell[1] + cell_idx_t(1); max[1][2] = cell[2];
+ min[2][0] = cell[0]; min[2][1] = cell[1]; min[2][2] = cell[2] - cell_idx_t(1);
+ max[2][0] = cell[0]; max[2][1] = cell[1]; max[2][2] = cell[2] + cell_idx_t(1);
+
+ for( uint_t f = 0; f < PdfField_T::F_SIZE; ++f )
+ {
+ WALBERLA_ASSERT( !math::isnan( tmpField->get( cell, f ) ) );
+
+ const auto v = tmpField->get( cell, f );
+
+ Vector3< real_t > grad( real_t(0) );
+
+ for( uint_t i = 0; i < PdfField_T::Stencil::D; ++i )
+ {
+
+#define WALBERLA_LBM_REFINEMENT_EXPLOSION_EXCLUDE_EXTRAPOLATION
+#ifdef WALBERLA_LBM_REFINEMENT_EXPLOSION_EXCLUDE_EXTRAPOLATION
+
+ if( boolField->get( max[i] ) && boolField->get( min[i] ) )
+ {
+ WALBERLA_ASSERT( !math::isnan( tmpField->get( max[i], f ) ) );
+ WALBERLA_ASSERT( !math::isnan( tmpField->get( min[i], f ) ) );
+
+ grad[i] = real_t(0.5) * ( tmpField->get( max[i], f ) - tmpField->get( min[i], f ) );
+ }
+
+#else
+ if( boolField->get( max[i] ) )
+ {
+ WALBERLA_ASSERT( !math::isnan( tmpField->get( max[i], f ) ) );
+ if( boolField->get( min[i] ) )
+ {
+ WALBERLA_ASSERT( !math::isnan( tmpField->get( min[i], f ) ) );
+ grad[i] = real_t(0.5) * ( tmpField->get( max[i], f ) - tmpField->get( min[i], f ) );
+ }
+ else
+ {
+ grad[i] = tmpField->get( max[i], f ) - v;
+ }
+
+ }
+ else if( boolField->get( min[i] ) )
+ {
+ WALBERLA_ASSERT( !math::isnan( tmpField->get( min[i], f ) ) );
+ grad[i] = v - tmpField->get( min[i], f );
+ }
+#endif
+#undef WALBERLA_LBM_REFINEMENT_EXPLOSION_EXCLUDE_EXTRAPOLATION
+ }
+
+#define WALBERLA_LBM_REFINEMENT_EXPLOSION_CHEN_CORRECTION
+#ifdef WALBERLA_LBM_REFINEMENT_EXPLOSION_CHEN_CORRECTION
+
+ const auto direction = PdfField_T::Stencil::dir[f];
+ Vector3< real_t > cNorm( stencil::cNorm[0][direction], stencil::cNorm[1][direction], stencil::cNorm[2][direction] );
+
+ grad = grad - cNorm * ( cNorm * grad );
+#endif
+#undef WALBERLA_LBM_REFINEMENT_EXPLOSION_CHEN_CORRECTION
+
+ const auto xx = x + cell_idx_t(1);
+ const auto yy = y + cell_idx_t(1);
+ pdfField->get( x , y , z , f ) = v + ( weights[0][0] * grad[0] + weights[0][1] * grad[1] + weights[0][2] * grad[2] );
+ pdfField->get( x , yy, z , f ) = v + ( weights[2][0] * grad[0] + weights[2][1] * grad[1] + weights[2][2] * grad[2] );
+ pdfField->get( xx, y , z , f ) = v + ( weights[4][0] * grad[0] + weights[4][1] * grad[1] + weights[4][2] * grad[2] );
+ pdfField->get( xx, yy, z , f ) = v + ( weights[6][0] * grad[0] + weights[6][1] * grad[1] + weights[6][2] * grad[2] );
+
+ if( PdfField_T::Stencil::D == uint_t(3) )
+ {
+ const auto zz = z + cell_idx_t(1);
+ pdfField->get( x , y , zz, f ) = v + ( weights[1][0] * grad[0] + weights[1][1] * grad[1] + weights[1][2] * grad[2] );
+ pdfField->get( x , yy, zz, f ) = v + ( weights[3][0] * grad[0] + weights[3][1] * grad[1] + weights[3][2] * grad[2] );
+ pdfField->get( xx, y , zz, f ) = v + ( weights[5][0] * grad[0] + weights[5][1] * grad[1] + weights[5][2] * grad[2] );
+ pdfField->get( xx, yy, zz, f ) = v + ( weights[7][0] * grad[0] + weights[7][1] * grad[1] + weights[7][2] * grad[2] );
+ }
+ }
+ }
+ else
+ {
+ WALBERLA_ASSERT( !boundaryHandling->isDomain( x + cell_idx_t(1), y, z ) );
+ WALBERLA_ASSERT( !boundaryHandling->isDomain( x, y + cell_idx_t(1), z ) );
+ WALBERLA_ASSERT( !boundaryHandling->isDomain( x + cell_idx_t(1), y + cell_idx_t(1), z ) );
+ if( PdfField_T::Stencil::D == uint_t(3) )
+ {
+ WALBERLA_ASSERT( !boundaryHandling->isDomain( x, y, z + cell_idx_t(1) ) );
+ WALBERLA_ASSERT( !boundaryHandling->isDomain( x + cell_idx_t(1), y, z + cell_idx_t(1) ) );
+ WALBERLA_ASSERT( !boundaryHandling->isDomain( x, y + cell_idx_t(1), z + cell_idx_t(1) ) );
+ WALBERLA_ASSERT( !boundaryHandling->isDomain( x + cell_idx_t(1), y + cell_idx_t(1), z + cell_idx_t(1) ) );
+ }
+ }
+ }
+}
+
+} // namespace internal
+
+
+
+} // namespace refinement
+} // namespace lbm
+} // namespace walberla
diff --git a/src/lbm/refinement_rebase/PdfFieldPackInfo.h b/src/lbm/refinement_rebase/PdfFieldPackInfo.h
new file mode 100644
index 000000000..6fbfb9221
--- /dev/null
+++ b/src/lbm/refinement_rebase/PdfFieldPackInfo.h
@@ -0,0 +1,1452 @@
+//======================================================================================================================
+//
+// 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 PdfFieldPackInfo.h
+//! \ingroup lbm
+//! \author Florian Schornbaum <florian.schornbaum@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "TimeStepPdfPackInfo.h"
+
+#include "lbm/field/PdfField.h"
+#include "blockforest/BlockNeighborhoodSection.h"
+#include "core/cell/CellInterval.h"
+#include "stencil/Directions.h"
+
+
+namespace walberla {
+namespace lbm {
+namespace refinement {
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+#endif
+class PdfFieldPackInfo : public TimeStepPdfPackInfo
+{
+public:
+
+ typedef PdfField<LatticeModel_T> PdfField_T;
+ typedef typename LatticeModel_T::Stencil Stencil;
+
+#ifdef NDEBUG
+ PdfFieldPackInfo( const BlockDataID & pdfFieldId, const bool _optimizedEqualLevelCommunication = true,
+ const bool _optimizedForLinearExplosion = true ) :
+ pdfFieldId_( pdfFieldId ), optimizedEqualLevelCommunication_( _optimizedEqualLevelCommunication ),
+ optimizedForLinearExplosion_( _optimizedForLinearExplosion ), equalLevelCells_( equalLevelCells() ) {}
+#else
+ PdfFieldPackInfo( const BlockDataID & pdfFieldId, const ConstBlockDataID & boundaryHandlingId,
+ const bool _optimizedEqualLevelCommunication = true, const bool _optimizedForLinearExplosion = true ) :
+ pdfFieldId_( pdfFieldId ), boundaryHandlingId_( boundaryHandlingId ),
+ optimizedEqualLevelCommunication_( _optimizedEqualLevelCommunication ), optimizedForLinearExplosion_( _optimizedForLinearExplosion ),
+ equalLevelCells_( equalLevelCells() ) {}
+#endif
+
+ virtual ~PdfFieldPackInfo() {}
+
+ bool optimizedEqualLevelCommunication() const { return optimizedEqualLevelCommunication_; }
+ void optimizeEqualLevelCommunication( const bool value = true ) { optimizedEqualLevelCommunication_ = value; }
+
+ bool optimizedForLinearExplosion() const { return optimizedForLinearExplosion_; }
+ void optimizeForLinearExplosion( const bool value = true ) { optimizedForLinearExplosion_ = value; }
+
+ bool constantDataExchange() const { return true; }
+ bool threadsafeReceiving() const { return true; }
+
+ void unpackDataEqualLevel( Block * receiver, stencil::Direction dir, mpi::RecvBuffer & buffer );
+ void communicateLocalEqualLevel( const Block * sender, Block * receiver, stencil::Direction dir );
+
+ void unpackDataCoarseToFine( Block * fineReceiver, const BlockID & coarseSender, stencil::Direction dir, mpi::RecvBuffer & buffer );
+ void communicateLocalCoarseToFine( const Block * coarseSender, Block * fineReceiver, stencil::Direction dir );
+
+ void unpackDataFineToCoarse( Block * coarseReceiver, const BlockID & fineSender, stencil::Direction dir, mpi::RecvBuffer & buffer );
+ void communicateLocalFineToCoarse( const Block * fineSender, Block * coarseReceiver, stencil::Direction dir );
+
+protected:
+
+ void packDataEqualLevelImpl( const Block * sender, stencil::Direction dir, mpi::SendBuffer & buffer ) const;
+ void packDataCoarseToFineImpl( const Block * coarseSender, const BlockID & fineReceiver, stencil::Direction dir, mpi::SendBuffer & buffer ) const;
+ void packDataFineToCoarseImpl( const Block * fineSender, const BlockID & coarseReceiver, stencil::Direction dir, mpi::SendBuffer & buffer ) const;
+
+ ///////////////////////////////////////////////////////////////////////
+ // Helper functions for determining packing/unpacking cell intervals //
+ ///////////////////////////////////////////////////////////////////////
+
+ static inline CellInterval equalLevelPackInterval ( stencil::Direction dir, const CellInterval & cellBB, const uint_t numberOfLayers );
+ static CellInterval equalLevelUnpackInterval( stencil::Direction dir, const CellInterval & cellBB, const uint_t numberOfLayers );
+
+ inline bool equalLevelFaceIntervalSplitable( const CellInterval & interval, stencil::Direction dir ) const;
+ inline std::vector< CellInterval > splitEqualLevelFaceInterval( const CellInterval & interval, stencil::Direction dir ) const; // dir: from sender to receiver
+
+ static inline Vector3< cell_idx_t > getNeighborShift( const BlockID & smallBlock, stencil::Direction dir ); // dir: direction from big to small block
+
+ static CellInterval coarseToFinePackInterval ( stencil::Direction dir, const CellInterval & cellBB, const BlockID & smallBlock );
+ static CellInterval coarseToFineUnpackInterval( stencil::Direction dir, const CellInterval & cellBB, const BlockID & smallBlock );
+
+ static inline CellInterval fineToCoarsePackInterval ( stencil::Direction dir, const CellInterval & cellBB );
+ static CellInterval fineToCoarseUnpackInterval( stencil::Direction dir, const CellInterval & cellBB, const BlockID & smallBlock );
+
+ //////////////////////////////
+ // General helper functions //
+ //////////////////////////////
+
+ static inline uint_t equalLevelCells();
+
+ static inline bool isFaceDirection( stencil::Direction dir );
+ static inline bool isEdgeDirection( stencil::Direction dir );
+ static inline bool isCornerDirection( stencil::Direction dir );
+
+ static inline bool blocksConnectedByFaces( const Block * block, const BlockID & neighbor );
+ static inline bool blocksConnectedByEdges( const Block * block, const BlockID & neighbor );
+
+ static inline bool divisibleByTwo( const CellInterval & cellBB );
+
+ static bool coarserNeighborExistsInVicinity( const Block * block, stencil::Direction dir );
+
+
+
+ BlockDataID pdfFieldId_;
+#ifndef NDEBUG
+ ConstBlockDataID boundaryHandlingId_;
+#endif
+
+ bool optimizedEqualLevelCommunication_;
+ bool optimizedForLinearExplosion_;
+
+ uint_t equalLevelCells_;
+};
+
+
+
+/////////////////
+// Equal level //
+/////////////////
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+void PdfFieldPackInfo< LatticeModel_T >::packDataEqualLevelImpl( const Block * sender, stencil::Direction dir, mpi::SendBuffer & buffer ) const
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+void PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::packDataEqualLevelImpl( const Block * sender, stencil::Direction dir,
+ mpi::SendBuffer & buffer ) const
+#endif
+{
+#ifndef NDEBUG
+ if( Stencil::D == uint_t(2) )
+ WALBERLA_ASSERT_EQUAL( stencil::cz[dir], 0 );
+#endif
+
+ const bool optimizedCommunication = ( optimizedEqualLevelCommunication_ && !coarserNeighborExistsInVicinity( sender, dir ) );
+
+ if( optimizedCommunication && Stencil::d_per_d_length[dir] == uint_t(0) )
+ return;
+
+ const PdfField_T * field = sender->getData< PdfField_T >( pdfFieldId_ );
+
+ if( optimizedCommunication )
+ {
+ CellInterval packingInterval = equalLevelPackInterval( dir, field->xyzSize(), uint_t(1) );
+
+ for( auto cell = field->beginSliceXYZ( packingInterval ); cell != field->end(); ++cell )
+ for( uint_t d = 0; d < Stencil::d_per_d_length[dir]; ++d )
+ buffer << cell.getF( Stencil::idx[ Stencil::d_per_d[dir][d] ] );
+
+ /*
+ * --> fzyx (sadly, seems to be slower ... ?)
+ *
+ for( auto z = packingInterval.zMin(); z <= packingInterval.zMax(); ++z ) {
+ for( auto y = packingInterval.yMin(); y <= packingInterval.yMax(); ++y ) {
+ for( uint_t d = 0; d < Stencil::d_per_d_length[dir]; ++d ) {
+
+ auto * s = &( field->get( packingInterval.xMin(), y, z, Stencil::d_per_d[dir][d] ) );
+
+ for( uint_t x = 0; x < packingInterval.xSize(); ++x )
+ buffer << s[x];
+ }
+ }
+ }
+ */
+ }
+ else
+ {
+ CellInterval packingInterval = equalLevelPackInterval( dir, field->xyzSize(), equalLevelCells_ );
+
+ if( optimizedEqualLevelCommunication_ && isFaceDirection( dir ) && Stencil::D == uint_t(3) )
+ {
+ WALBERLA_ASSERT( equalLevelFaceIntervalSplitable( packingInterval, dir ) );
+
+ std::vector< CellInterval > intervals = splitEqualLevelFaceInterval( packingInterval, dir );
+ WALBERLA_ASSERT_EQUAL( intervals.size(), 5 );
+
+ for( uint_t i = 0; i < 4; ++i )
+ {
+ for( auto cell = field->beginSliceXYZ( intervals[i] ); cell != field->end(); ++cell )
+ for( uint_t d = 0; d < Stencil::Size; ++d )
+ buffer << cell.getF( d );
+ }
+
+ if( Stencil::d_per_d_length[dir] > uint_t(0) )
+ for( auto cell = field->beginSliceXYZ( intervals[4] ); cell != field->end(); ++cell )
+ for( uint_t d = 0; d < Stencil::d_per_d_length[dir]; ++d )
+ buffer << cell.getF( Stencil::idx[ Stencil::d_per_d[dir][d] ] );
+ }
+ else
+ {
+ for( auto cell = field->beginSliceXYZ( packingInterval ); cell != field->end(); ++cell )
+ for( uint_t d = 0; d < Stencil::Size; ++d )
+ buffer << cell.getF( d );
+
+ /*
+ * --> fzyx (sadly, seems to be slower ... ?)
+ *
+ for( auto z = packingInterval.zMin(); z <= packingInterval.zMax(); ++z ) {
+ for( auto y = packingInterval.yMin(); y <= packingInterval.yMax(); ++y ) {
+ for( uint_t d = 0; d < Stencil::Size; ++d ) {
+
+ auto * s = &( field->get( packingInterval.xMin(), y, z, d ) );
+
+ for( uint_t x = 0; x < packingInterval.xSize(); ++x )
+ buffer << s[x];
+ }
+ }
+ }
+ */
+ }
+ }
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+void PdfFieldPackInfo< LatticeModel_T >::unpackDataEqualLevel( Block * receiver, stencil::Direction dir, mpi::RecvBuffer & buffer )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+void PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::unpackDataEqualLevel( Block * receiver, stencil::Direction dir,
+ mpi::RecvBuffer & buffer )
+#endif
+{
+#ifndef NDEBUG
+ if( Stencil::D == uint_t(2) )
+ WALBERLA_ASSERT_EQUAL( stencil::cz[dir], 0 );
+#endif
+
+ const bool optimizedCommunication = ( optimizedEqualLevelCommunication_ && !coarserNeighborExistsInVicinity( receiver, dir ) );
+
+ const auto invDir = stencil::inverseDir[dir];
+
+ if( optimizedCommunication && Stencil::d_per_d_length[invDir] == uint_t(0) )
+ return;
+
+ PdfField_T * field = receiver->getData< PdfField_T >( pdfFieldId_ );
+
+ if( optimizedCommunication )
+ {
+ CellInterval unpackingInterval = equalLevelUnpackInterval( dir, field->xyzSize(), uint_t(1) );
+
+ for( auto cell = field->beginSliceXYZ( unpackingInterval ); cell != field->end(); ++cell )
+ for( uint_t d = 0; d < Stencil::d_per_d_length[invDir]; ++d )
+ buffer >> cell.getF( Stencil::idx[ Stencil::d_per_d[invDir][d] ] );
+
+ /*
+ * --> fzyx (sadly, seems to be slower ... ?)
+ *
+ for( auto z = unpackingInterval.zMin(); z <= unpackingInterval.zMax(); ++z ) {
+ for( auto y = unpackingInterval.yMin(); y <= unpackingInterval.yMax(); ++y ) {
+ for( uint_t d = 0; d < Stencil::d_per_d_length[invDir]; ++d ) {
+
+ auto * r = &( field->get( unpackingInterval.xMin(), y, z, Stencil::d_per_d[invDir][d] ) );
+
+ for( uint_t x = 0; x < unpackingInterval.xSize(); ++x )
+ buffer >> r[x];
+ }
+ }
+ }
+ */
+ }
+ else
+ {
+ CellInterval unpackingInterval = equalLevelUnpackInterval( dir, field->xyzSize(), equalLevelCells_ );
+
+ if( optimizedEqualLevelCommunication_ && isFaceDirection( dir ) && Stencil::D == uint_t(3) )
+ {
+ WALBERLA_ASSERT( equalLevelFaceIntervalSplitable( unpackingInterval, invDir ) );
+
+ std::vector< CellInterval > intervals = splitEqualLevelFaceInterval( unpackingInterval, invDir );
+ WALBERLA_ASSERT_EQUAL( intervals.size(), 5 );
+
+ for( uint_t i = 0; i < 4; ++i )
+ {
+ for( auto cell = field->beginSliceXYZ( intervals[i] ); cell != field->end(); ++cell )
+ for( uint_t d = 0; d < Stencil::Size; ++d )
+ buffer >> cell.getF( d );
+ }
+
+ if( Stencil::d_per_d_length[invDir] > uint_t(0) )
+ for( auto cell = field->beginSliceXYZ( intervals[4] ); cell != field->end(); ++cell )
+ for( uint_t d = 0; d < Stencil::d_per_d_length[invDir]; ++d )
+ buffer >> cell.getF( Stencil::idx[ Stencil::d_per_d[invDir][d] ] );
+ }
+ else
+ {
+ for( auto cell = field->beginSliceXYZ( unpackingInterval ); cell != field->end(); ++cell )
+ for( uint_t d = 0; d < Stencil::Size; ++d )
+ buffer >> cell.getF( d );
+
+ /*
+ * --> fzyx (sadly, seems to be slower ... ?)
+ *
+ for( auto z = unpackingInterval.zMin(); z <= unpackingInterval.zMax(); ++z ) {
+ for( auto y = unpackingInterval.yMin(); y <= unpackingInterval.yMax(); ++y ) {
+ for( uint_t d = 0; d < Stencil::Size; ++d ) {
+
+ auto * r = &( field->get( unpackingInterval.xMin(), y, z, d ) );
+
+ for( uint_t x = 0; x < unpackingInterval.xSize(); ++x )
+ buffer >> r[x];
+ }
+ }
+ }
+ */
+ }
+ }
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+void PdfFieldPackInfo< LatticeModel_T >::communicateLocalEqualLevel( const Block * sender, Block * receiver, stencil::Direction dir )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+void PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::communicateLocalEqualLevel( const Block * sender, Block * receiver, stencil::Direction dir )
+#endif
+{
+#ifndef NDEBUG
+ if( Stencil::D == uint_t(2) )
+ WALBERLA_ASSERT_EQUAL( stencil::cz[dir], 0 );
+#endif
+
+ const bool optimizedCommunication = ( optimizedEqualLevelCommunication_ && !coarserNeighborExistsInVicinity( sender, dir ) );
+
+ if( optimizedCommunication && Stencil::d_per_d_length[dir] == uint_t(0) )
+ return;
+
+ const PdfField_T * sf = sender->getData< PdfField_T >( pdfFieldId_ );
+ PdfField_T * rf = receiver->getData< PdfField_T >( pdfFieldId_ );
+
+ WALBERLA_ASSERT_EQUAL( sf->xyzSize(), rf->xyzSize() );
+
+ if( optimizedCommunication )
+ {
+ CellInterval packingInterval = equalLevelPackInterval( dir, sf->xyzSize(), uint_t(1) );
+ CellInterval unpackingInterval = equalLevelUnpackInterval( stencil::inverseDir[dir], rf->xyzSize(), uint_t(1) );
+
+ WALBERLA_ASSERT_EQUAL( packingInterval.xSize(), unpackingInterval.xSize() );
+ WALBERLA_ASSERT_EQUAL( packingInterval.ySize(), unpackingInterval.ySize() );
+ WALBERLA_ASSERT_EQUAL( packingInterval.zSize(), unpackingInterval.zSize() );
+
+ auto sCell = sf->beginSliceXYZ( packingInterval );
+ auto rCell = rf->beginSliceXYZ( unpackingInterval );
+ while( sCell != sf->end() )
+ {
+ WALBERLA_ASSERT( rCell != rf->end() );
+
+ for( uint_t d = 0; d < Stencil::d_per_d_length[dir]; ++d )
+ {
+ const auto idx = Stencil::idx[ Stencil::d_per_d[dir][d] ];
+ rCell.getF( idx ) = sCell.getF( idx );
+ }
+
+ ++sCell;
+ ++rCell;
+ }
+ WALBERLA_ASSERT( rCell == rf->end() );
+
+ /*
+ * --> fzyx (sadly, seems to be slower ... ?)
+ *
+ for( auto zp = packingInterval.zMin(), zu = unpackingInterval.zMin(); zp <= packingInterval.zMax(); ++zp, ++zu ) {
+ for( auto yp = packingInterval.yMin(), yu = unpackingInterval.yMin(); yp <= packingInterval.yMax(); ++yp, ++yu ) {
+ for( uint_t d = 0; d < Stencil::d_per_d_length[dir]; ++d ) {
+
+ auto * s = &( sf->get( packingInterval.xMin(), yp, zp, Stencil::d_per_d[dir][d] ) );
+ auto * r = &( rf->get( unpackingInterval.xMin(), yu, zu, Stencil::d_per_d[dir][d] ) );
+
+ std::copy( s, s+packingInterval.xSize(), r );
+ }
+ }
+ }
+ */
+ }
+ else
+ {
+ CellInterval packingInterval = equalLevelPackInterval( dir, sf->xyzSize(), equalLevelCells_ );
+ CellInterval unpackingInterval = equalLevelUnpackInterval( stencil::inverseDir[dir], rf->xyzSize(), equalLevelCells_ );
+
+ WALBERLA_ASSERT_EQUAL( packingInterval.xSize(), unpackingInterval.xSize() );
+ WALBERLA_ASSERT_EQUAL( packingInterval.ySize(), unpackingInterval.ySize() );
+ WALBERLA_ASSERT_EQUAL( packingInterval.zSize(), unpackingInterval.zSize() );
+
+ if( optimizedEqualLevelCommunication_ && isFaceDirection( dir ) && Stencil::D == uint_t(3) )
+ {
+ WALBERLA_ASSERT( equalLevelFaceIntervalSplitable( packingInterval, dir ) );
+ WALBERLA_ASSERT( equalLevelFaceIntervalSplitable( unpackingInterval, dir ) );
+
+ std::vector< CellInterval > pIntervals = splitEqualLevelFaceInterval( packingInterval, dir );
+ std::vector< CellInterval > uIntervals = splitEqualLevelFaceInterval( unpackingInterval, dir );
+ WALBERLA_ASSERT_EQUAL( pIntervals.size(), 5 );
+ WALBERLA_ASSERT_EQUAL( uIntervals.size(), 5 );
+
+ for( uint_t i = 0; i < 4; ++i )
+ {
+ auto sCell = sf->beginSliceXYZ( pIntervals[i] );
+ auto rCell = rf->beginSliceXYZ( uIntervals[i] );
+ while( sCell != sf->end() )
+ {
+ WALBERLA_ASSERT( rCell != rf->end() );
+
+ for( uint_t d = 0; d < Stencil::Size; ++d )
+ rCell.getF( d ) = sCell.getF( d );
+
+ ++sCell;
+ ++rCell;
+ }
+ WALBERLA_ASSERT( rCell == rf->end() );
+ }
+
+ if( Stencil::d_per_d_length[dir] > uint_t(0) )
+ {
+ auto sCell = sf->beginSliceXYZ( pIntervals[4] );
+ auto rCell = rf->beginSliceXYZ( uIntervals[4] );
+ while( sCell != sf->end() )
+ {
+ WALBERLA_ASSERT( rCell != rf->end() );
+
+ for( uint_t d = 0; d < Stencil::d_per_d_length[dir]; ++d )
+ {
+ const auto idx = Stencil::idx[ Stencil::d_per_d[dir][d] ];
+ rCell.getF( idx ) = sCell.getF( idx );
+ }
+
+ ++sCell;
+ ++rCell;
+ }
+ WALBERLA_ASSERT( rCell == rf->end() );
+ }
+ }
+ else
+ {
+ auto sCell = sf->beginSliceXYZ( packingInterval );
+ auto rCell = rf->beginSliceXYZ( unpackingInterval );
+ while( sCell != sf->end() )
+ {
+ WALBERLA_ASSERT( rCell != rf->end() );
+
+ for( uint_t d = 0; d < Stencil::Size; ++d )
+ rCell.getF( d ) = sCell.getF( d );
+
+ ++sCell;
+ ++rCell;
+ }
+ WALBERLA_ASSERT( rCell == rf->end() );
+
+ /*
+ * --> fzyx (sadly, seems to be slower ... ?)
+ *
+ for( auto zp = packingInterval.zMin(), zu = unpackingInterval.zMin(); zp <= packingInterval.zMax(); ++zp, ++zu ) {
+ for( auto yp = packingInterval.yMin(), yu = unpackingInterval.yMin(); yp <= packingInterval.yMax(); ++yp, ++yu ) {
+ for( uint_t d = 0; d < Stencil::Size; ++d ) {
+
+ auto * s = &( sf->get( packingInterval.xMin(), yp, zp, d ) );
+ auto * r = &( rf->get( unpackingInterval.xMin(), yu, zu, d ) );
+
+ std::copy( s, s+packingInterval.xSize(), r );
+ }
+ }
+ }
+ */
+ }
+ }
+}
+
+
+
+////////////////////
+// Coarse to fine //
+////////////////////
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+void PdfFieldPackInfo< LatticeModel_T >::packDataCoarseToFineImpl( const Block * coarseSender, const BlockID & fineReceiver,
+ stencil::Direction dir, mpi::SendBuffer & buffer ) const
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+void PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::packDataCoarseToFineImpl( const Block * coarseSender, const BlockID & fineReceiver,
+ stencil::Direction dir, mpi::SendBuffer & buffer ) const
+#endif
+{
+#ifndef NDEBUG
+ if( Stencil::D == uint_t(2) )
+ WALBERLA_ASSERT_EQUAL( stencil::cz[dir], 0 );
+#endif
+
+ const PdfField_T * field = coarseSender->getData< PdfField_T >( pdfFieldId_ );
+
+#ifndef NDEBUG
+ const BoundaryHandling_T * boundaryHandling = coarseSender->template getData< BoundaryHandling_T >( boundaryHandlingId_ );
+#endif
+
+ CellInterval packingInterval = coarseToFinePackInterval( dir, field->xyzSize(), fineReceiver );
+
+ for( cell_idx_t z = packingInterval.zMin(); z <= packingInterval.zMax(); ++z ) {
+ for( cell_idx_t y = packingInterval.yMin(); y <= packingInterval.yMax(); ++y ) {
+ for( cell_idx_t x = packingInterval.xMin(); x <= packingInterval.xMax(); ++x ) {
+ for( uint_t d = 0; d < Stencil::Size; ++d )
+ //for( uint_t d = 0; d < Stencil::d_per_d_length[dir]; ++d )
+ {
+ buffer << field->get( x, y, z, d );
+ //buffer << field->get( x, y, z, Stencil::idx[ Stencil::d_per_d[dir][d] ] );
+
+#ifndef NDEBUG
+ if( boundaryHandling->isDomain(x,y,z) )
+ WALBERLA_ASSERT( !math::isnan( field->get( x, y, z, d ) ) );
+#endif
+ }
+ }
+ }
+ }
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+void PdfFieldPackInfo< LatticeModel_T >::unpackDataCoarseToFine( Block * fineReceiver, const BlockID & /*coarseSender*/,
+ stencil::Direction dir, mpi::RecvBuffer & buffer )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+void PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::unpackDataCoarseToFine( Block * fineReceiver, const BlockID & /*coarseSender*/,
+ stencil::Direction dir, mpi::RecvBuffer & buffer )
+#endif
+{
+#ifndef NDEBUG
+ if( Stencil::D == uint_t(2) )
+ WALBERLA_ASSERT_EQUAL( stencil::cz[dir], 0 );
+#endif
+
+ PdfField_T * field = fineReceiver->getData< PdfField_T >( pdfFieldId_ );
+
+ CellInterval unpackingInterval = coarseToFineUnpackInterval( dir, field->xyzSize(), fineReceiver->getId() );
+
+ if( optimizedForLinearExplosion_ )
+ {
+ for( cell_idx_t z = unpackingInterval.zMin(); z <= unpackingInterval.zMax(); z += cell_idx_t(2) ) {
+ for( cell_idx_t y = unpackingInterval.yMin(); y <= unpackingInterval.yMax(); y += cell_idx_t(2) ) {
+ for( cell_idx_t x = unpackingInterval.xMin(); x <= unpackingInterval.xMax(); x += cell_idx_t(2) ) {
+ for( uint_t idx = 0; idx < Stencil::Size; ++idx )
+ {
+ typename PdfField_T::value_type value;
+ buffer >> value;
+ field->get( x, y, z, idx ) = value;
+ }
+ }
+ }
+ }
+ }
+ else
+ {
+ for( cell_idx_t z = unpackingInterval.zMin(); z <= unpackingInterval.zMax(); z += cell_idx_t(2) ) {
+ for( cell_idx_t y = unpackingInterval.yMin(); y <= unpackingInterval.yMax(); y += cell_idx_t(2) ) {
+ for( cell_idx_t x = unpackingInterval.xMin(); x <= unpackingInterval.xMax(); x += cell_idx_t(2) ) {
+ for( uint_t idx = 0; idx < Stencil::Size; ++idx )
+ //for( uint_t d = 0; d < Stencil::d_per_d_length[ stencil::inverseDir[dir] ]; ++d )
+ {
+ typename PdfField_T::value_type value;
+ buffer >> value;
+
+ //const auto idx = Stencil::idx[ Stencil::d_per_d[ stencil::inverseDir[dir] ][d] ];
+
+ field->get( x, y, z, idx ) = value;
+ field->get( x + cell_idx_t(1), y, z, idx ) = value;
+ field->get( x, y + cell_idx_t(1), z, idx ) = value;
+ field->get( x + cell_idx_t(1), y + cell_idx_t(1), z, idx ) = value;
+ if( Stencil::D == uint_t(3) )
+ {
+ field->get( x, y, z + cell_idx_t(1), idx ) = value;
+ field->get( x + cell_idx_t(1), y, z + cell_idx_t(1), idx ) = value;
+ field->get( x, y + cell_idx_t(1), z + cell_idx_t(1), idx ) = value;
+ field->get( x + cell_idx_t(1), y + cell_idx_t(1), z + cell_idx_t(1), idx ) = value;
+ }
+ }
+ }
+ }
+ }
+ }
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+void PdfFieldPackInfo< LatticeModel_T >::communicateLocalCoarseToFine( const Block * coarseSender, Block * fineReceiver, stencil::Direction dir )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+void PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::communicateLocalCoarseToFine( const Block * coarseSender, Block * fineReceiver,
+ stencil::Direction dir )
+#endif
+{
+#ifndef NDEBUG
+ if( Stencil::D == uint_t(2) )
+ WALBERLA_ASSERT_EQUAL( stencil::cz[dir], 0 );
+#endif
+
+ const PdfField_T * sf = coarseSender->getData< PdfField_T >( pdfFieldId_ );
+ PdfField_T * rf = fineReceiver->getData< PdfField_T >( pdfFieldId_ );
+
+#ifndef NDEBUG
+ const BoundaryHandling_T * boundaryHandling = coarseSender->template getData< BoundaryHandling_T >( boundaryHandlingId_ );
+#endif
+
+ WALBERLA_ASSERT_EQUAL( sf->xyzSize(), rf->xyzSize() );
+
+ CellInterval packingInterval = coarseToFinePackInterval( dir, sf->xyzSize(), fineReceiver->getId() );
+ CellInterval unpackingInterval = coarseToFineUnpackInterval( stencil::inverseDir[dir], rf->xyzSize(), fineReceiver->getId() );
+
+ WALBERLA_ASSERT_EQUAL( packingInterval.xSize() * uint_t(2), unpackingInterval.xSize() );
+ WALBERLA_ASSERT_EQUAL( packingInterval.ySize() * uint_t(2), unpackingInterval.ySize() );
+ if( Stencil::D == uint_t(3) )
+ {
+ WALBERLA_ASSERT_EQUAL( packingInterval.zSize() * uint_t(2), unpackingInterval.zSize() );
+ }
+ else
+ {
+ WALBERLA_ASSERT_EQUAL( packingInterval.zSize(), unpackingInterval.zSize() );
+ }
+
+ if( optimizedForLinearExplosion_ )
+ {
+ cell_idx_t rz = unpackingInterval.zMin();
+ for( cell_idx_t sz = packingInterval.zMin(); sz <= packingInterval.zMax(); ++sz )
+ {
+ cell_idx_t ry = unpackingInterval.yMin();
+ for( cell_idx_t sy = packingInterval.yMin(); sy <= packingInterval.yMax(); ++sy )
+ {
+ cell_idx_t rx = unpackingInterval.xMin();
+ for( cell_idx_t sx = packingInterval.xMin(); sx <= packingInterval.xMax(); ++sx ) {
+ for( uint_t idx = 0; idx < Stencil::Size; ++idx )
+ {
+ rf->get( rx, ry, rz, idx ) = sf->get( sx, sy, sz, idx );
+
+#ifndef NDEBUG
+ if( boundaryHandling->isDomain(sx,sy,sz) )
+ WALBERLA_ASSERT( !math::isnan( sf->get( sx, sy, sz, idx ) ), sx << ", " << sy << ", " << sz << ", " << idx << " coarse sender block = " << coarseSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> " );
+#endif
+ }
+ rx += cell_idx_t(2);
+ }
+ ry += cell_idx_t(2);
+ WALBERLA_ASSERT_EQUAL( rx, unpackingInterval.xMax() + cell_idx_t(1) );
+ }
+ rz += cell_idx_t(2);
+ WALBERLA_ASSERT_EQUAL( ry, unpackingInterval.yMax() + cell_idx_t(1) );
+ }
+ if( Stencil::D == uint_t(3) )
+ {
+ WALBERLA_ASSERT_EQUAL( rz, unpackingInterval.zMax() + cell_idx_t(1) );
+ }
+ else
+ {
+ WALBERLA_ASSERT_EQUAL( rz, unpackingInterval.zMax() + cell_idx_t(2) );
+ }
+ }
+ else
+ {
+ cell_idx_t rz = unpackingInterval.zMin();
+ for( cell_idx_t sz = packingInterval.zMin(); sz <= packingInterval.zMax(); ++sz )
+ {
+ cell_idx_t ry = unpackingInterval.yMin();
+ for( cell_idx_t sy = packingInterval.yMin(); sy <= packingInterval.yMax(); ++sy )
+ {
+ cell_idx_t rx = unpackingInterval.xMin();
+ for( cell_idx_t sx = packingInterval.xMin(); sx <= packingInterval.xMax(); ++sx ) {
+ for( uint_t idx = 0; idx < Stencil::Size; ++idx )
+ //for( uint_t d = 0; d < Stencil::d_per_d_length[dir]; ++d )
+ {
+ //const auto idx = Stencil::idx[ Stencil::d_per_d[dir][d] ];
+
+ const auto & value = sf->get(sx,sy,sz,idx);
+
+#ifndef NDEBUG
+ if( boundaryHandling->isDomain(sx,sy,sz) )
+ WALBERLA_ASSERT( !math::isnan( value ), "value at " << sx << ", " << sy << ", " << sz << ", " << idx << " coarse sender block = " << coarseSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> " );
+#endif
+
+ rf->get( rx, ry, rz, idx ) = value;
+ rf->get( rx + cell_idx_t(1), ry, rz, idx ) = value;
+ rf->get( rx, ry + cell_idx_t(1), rz, idx ) = value;
+ rf->get( rx + cell_idx_t(1), ry + cell_idx_t(1), rz, idx ) = value;
+ if( Stencil::D == uint_t(3) )
+ {
+ rf->get( rx, ry, rz + cell_idx_t(1), idx ) = value;
+ rf->get( rx + cell_idx_t(1), ry, rz + cell_idx_t(1), idx ) = value;
+ rf->get( rx, ry + cell_idx_t(1), rz + cell_idx_t(1), idx ) = value;
+ rf->get( rx + cell_idx_t(1), ry + cell_idx_t(1), rz + cell_idx_t(1), idx ) = value;
+ }
+ }
+ rx += cell_idx_t(2);
+ }
+ ry += cell_idx_t(2);
+ WALBERLA_ASSERT_EQUAL( rx, unpackingInterval.xMax() + cell_idx_t(1) );
+ }
+ rz += cell_idx_t(2);
+ WALBERLA_ASSERT_EQUAL( ry, unpackingInterval.yMax() + cell_idx_t(1) );
+ }
+ if( Stencil::D == uint_t(3) )
+ {
+ WALBERLA_ASSERT_EQUAL( rz, unpackingInterval.zMax() + cell_idx_t(1) );
+ }
+ else
+ {
+ WALBERLA_ASSERT_EQUAL( rz, unpackingInterval.zMax() + cell_idx_t(2) );
+ }
+ }
+}
+
+
+
+////////////////////
+// Fine to coarse //
+////////////////////
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+void PdfFieldPackInfo< LatticeModel_T >::packDataFineToCoarseImpl( const Block * fineSender, const BlockID & coarseReceiver,
+ stencil::Direction dir, mpi::SendBuffer & buffer ) const
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+void PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::packDataFineToCoarseImpl( const Block * fineSender, const BlockID & coarseReceiver,
+ stencil::Direction dir, mpi::SendBuffer & buffer ) const
+#endif
+{
+#ifndef NDEBUG
+ if( Stencil::D == uint_t(2) )
+ WALBERLA_ASSERT_EQUAL( stencil::cz[dir], 0 );
+#endif
+
+ if( Stencil::d_per_d_length[dir] == uint_t(0) )
+ return;
+
+ if( ( ( isEdgeDirection(dir) || isCornerDirection(dir) ) && blocksConnectedByFaces( fineSender, coarseReceiver ) ) ||
+ ( isCornerDirection(dir) && blocksConnectedByEdges( fineSender, coarseReceiver ) ) )
+ return;
+
+ const PdfField_T * field = fineSender->getData< PdfField_T >( pdfFieldId_ );
+
+#ifndef NDEBUG
+ const BoundaryHandling_T * boundaryHandling = fineSender->template getData< BoundaryHandling_T >( boundaryHandlingId_ );
+#endif
+
+ CellInterval packingInterval = fineToCoarsePackInterval( dir, field->xyzSize() );
+
+ const real_t factor = ( Stencil::D == uint_t(3) ) ? real_t( 0.125 ) : real_t( 0.25 );
+
+ for( cell_idx_t z = packingInterval.zMin(); z <= packingInterval.zMax(); z += cell_idx_t(2) ) {
+ for( cell_idx_t y = packingInterval.yMin(); y <= packingInterval.yMax(); y += cell_idx_t(2) ) {
+ for( cell_idx_t x = packingInterval.xMin(); x <= packingInterval.xMax(); x += cell_idx_t(2) ) {
+ for( uint_t d = 0; d < Stencil::d_per_d_length[dir]; ++d )
+ {
+ const auto idx = Stencil::idx[ Stencil::d_per_d[dir][d] ];
+
+#ifndef NDEBUG
+ if( boundaryHandling->isDomain(x,y,z) )
+ {
+ WALBERLA_ASSERT( !math::isnan( field->get( x, y, z, idx ) ), x << ", " << y << ", " << z << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> ");
+ WALBERLA_ASSERT( !math::isnan( field->get( x + cell_idx_t(1), y, z, idx ) ), x + cell_idx_t(1) << ", " << y << ", " << z << ", " << idx <<" fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> " );
+ WALBERLA_ASSERT( !math::isnan( field->get( x, y + cell_idx_t(1), z, idx ) ), x << ", " << y + cell_idx_t(1) << ", " << z << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> " );
+ WALBERLA_ASSERT( !math::isnan( field->get( x + cell_idx_t(1), y + cell_idx_t(1), z, idx ) ), x + cell_idx_t(1) << ", " << y + cell_idx_t(1) << ", " << z << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> " );
+ if( Stencil::D == uint_t(3) )
+ {
+ WALBERLA_ASSERT( !math::isnan( field->get( x, y, z + cell_idx_t(1), idx ) ), x << ", " << y << ", " << z + cell_idx_t(1) << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> " );
+ WALBERLA_ASSERT( !math::isnan( field->get( x + cell_idx_t(1), y, z + cell_idx_t(1), idx ) ), x + cell_idx_t(1) << ", " << y << ", " << z + cell_idx_t(1) << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> " );
+ WALBERLA_ASSERT( !math::isnan( field->get( x, y + cell_idx_t(1), z + cell_idx_t(1), idx ) ), x << ", " << y + cell_idx_t(1) << ", " << z + cell_idx_t(1) << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> " );
+ WALBERLA_ASSERT( !math::isnan( field->get( x + cell_idx_t(1), y + cell_idx_t(1), z + cell_idx_t(1), idx ) ), x + cell_idx_t(1) << ", " << y + cell_idx_t(1) << ", " << z + cell_idx_t(1) << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> " );
+ }
+ }
+#endif
+
+ typename PdfField_T::value_type value = field->get( x, y, z, idx );
+ value += field->get( x + cell_idx_t(1), y, z, idx );
+ value += field->get( x, y + cell_idx_t(1), z, idx );
+ value += field->get( x + cell_idx_t(1), y + cell_idx_t(1), z, idx );
+ if( Stencil::D == uint_t(3) )
+ {
+ value += field->get( x, y, z + cell_idx_t(1), idx );
+ value += field->get( x + cell_idx_t(1), y, z + cell_idx_t(1), idx );
+ value += field->get( x, y + cell_idx_t(1), z + cell_idx_t(1), idx );
+ value += field->get( x + cell_idx_t(1), y + cell_idx_t(1), z + cell_idx_t(1), idx );
+ }
+
+ buffer << ( factor * value );
+ }
+ }
+ }
+ }
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+void PdfFieldPackInfo< LatticeModel_T >::unpackDataFineToCoarse( Block * coarseReceiver, const BlockID & fineSender,
+ stencil::Direction dir, mpi::RecvBuffer & buffer )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+void PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::unpackDataFineToCoarse( Block * coarseReceiver, const BlockID & fineSender,
+ stencil::Direction dir, mpi::RecvBuffer & buffer )
+#endif
+{
+#ifndef NDEBUG
+ if( Stencil::D == uint_t(2) )
+ WALBERLA_ASSERT_EQUAL( stencil::cz[dir], 0 );
+#endif
+
+ auto invDir = stencil::inverseDir[dir];
+ if( Stencil::d_per_d_length[invDir] == uint_t(0) )
+ return;
+
+ if( ( ( isEdgeDirection(dir) || isCornerDirection(dir) ) && blocksConnectedByFaces( coarseReceiver, fineSender ) ) ||
+ ( isCornerDirection(dir) && blocksConnectedByEdges( coarseReceiver, fineSender ) ) )
+ return;
+
+ PdfField_T * field = coarseReceiver->getData< PdfField_T >( pdfFieldId_ );
+
+#ifndef NDEBUG
+ const BoundaryHandling_T * boundaryHandling = coarseReceiver->template getData< BoundaryHandling_T >( boundaryHandlingId_ );
+#endif
+
+ CellInterval unpackingInterval = fineToCoarseUnpackInterval( dir, field->xyzSize(), fineSender );
+
+ for( cell_idx_t z = unpackingInterval.zMin(); z <= unpackingInterval.zMax(); ++z ) {
+ for( cell_idx_t y = unpackingInterval.yMin(); y <= unpackingInterval.yMax(); ++y ) {
+ for( cell_idx_t x = unpackingInterval.xMin(); x <= unpackingInterval.xMax(); ++x ) {
+ for( uint_t d = 0; d < Stencil::d_per_d_length[invDir]; ++d )
+ {
+ buffer >> field->get( x, y, z, Stencil::idx[ Stencil::d_per_d[invDir][d] ] );
+
+#ifndef NDEBUG
+ if( boundaryHandling->isDomain(x,y,z) )
+ WALBERLA_ASSERT( !math::isnan( field->get( x, y, z, Stencil::idx[ Stencil::d_per_d[invDir][d] ] ) ) );
+#endif
+ }
+ }
+ }
+ }
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+void PdfFieldPackInfo< LatticeModel_T >::communicateLocalFineToCoarse( const Block * fineSender, Block * coarseReceiver, stencil::Direction dir )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+void PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::communicateLocalFineToCoarse( const Block * fineSender, Block * coarseReceiver,
+ stencil::Direction dir )
+#endif
+{
+#ifndef NDEBUG
+ if( Stencil::D == uint_t(2) )
+ WALBERLA_ASSERT_EQUAL( stencil::cz[dir], 0 );
+#endif
+
+ if( Stencil::d_per_d_length[dir] == uint_t(0) )
+ return;
+
+ if( ( ( isEdgeDirection(dir) || isCornerDirection(dir) ) && blocksConnectedByFaces( fineSender, coarseReceiver->getId() ) ) ||
+ ( isCornerDirection(dir) && blocksConnectedByEdges( fineSender, coarseReceiver->getId() ) ) )
+ return;
+
+ const PdfField_T * sf = fineSender->getData< PdfField_T >( pdfFieldId_ );
+ PdfField_T * rf = coarseReceiver->getData< PdfField_T >( pdfFieldId_ );
+
+#ifndef NDEBUG
+ const BoundaryHandling_T * boundaryHandling = fineSender->template getData< BoundaryHandling_T >( boundaryHandlingId_ );
+#endif
+
+ WALBERLA_ASSERT_EQUAL( sf->xyzSize(), rf->xyzSize() );
+
+ CellInterval packingInterval = fineToCoarsePackInterval( dir, sf->xyzSize() );
+ CellInterval unpackingInterval = fineToCoarseUnpackInterval( stencil::inverseDir[dir], rf->xyzSize(), fineSender->getId() );
+
+ WALBERLA_ASSERT_EQUAL( packingInterval.xSize(), unpackingInterval.xSize() * uint_t(2) );
+ WALBERLA_ASSERT_EQUAL( packingInterval.ySize(), unpackingInterval.ySize() * uint_t(2) );
+ if( Stencil::D == uint_t(3) )
+ {
+ WALBERLA_ASSERT_EQUAL( packingInterval.zSize(), unpackingInterval.zSize() * uint_t(2) );
+ }
+ else
+ {
+ WALBERLA_ASSERT_EQUAL( packingInterval.zSize(), unpackingInterval.zSize() );
+ }
+
+ const real_t factor = ( Stencil::D == uint_t(3) ) ? real_t( 0.125 ) : real_t( 0.25 );
+
+ cell_idx_t sz = packingInterval.zMin();
+ for( cell_idx_t rz = unpackingInterval.zMin(); rz <= unpackingInterval.zMax(); ++rz )
+ {
+ cell_idx_t sy = packingInterval.yMin();
+ for( cell_idx_t ry = unpackingInterval.yMin(); ry <= unpackingInterval.yMax(); ++ry )
+ {
+ cell_idx_t sx = packingInterval.xMin();
+ for( cell_idx_t rx = unpackingInterval.xMin(); rx <= unpackingInterval.xMax(); ++rx ) {
+ for( uint_t d = 0; d < Stencil::d_per_d_length[dir]; ++d )
+ {
+ const auto idx = Stencil::idx[ Stencil::d_per_d[dir][d] ];
+
+#ifndef NDEBUG
+ if( boundaryHandling->isDomain(sx,sy,sz) )
+ {
+ WALBERLA_ASSERT( !math::isnan( sf->get( sx, sy, sz, idx ) ), sx << ", " << sy << ", " << sz << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> ");
+ WALBERLA_ASSERT( !math::isnan( sf->get( sx + cell_idx_t(1), sy, sz, idx ) ), sx + cell_idx_t(1) << ", " << sy << ", " << sz << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> ");
+ WALBERLA_ASSERT( !math::isnan( sf->get( sx, sy + cell_idx_t(1), sz, idx ) ), sx << ", " << sy + cell_idx_t(1) << ", " << sz << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> ");
+ WALBERLA_ASSERT( !math::isnan( sf->get( sx + cell_idx_t(1), sy + cell_idx_t(1), sz, idx ) ), sx + cell_idx_t(1) << ", " << sy + cell_idx_t(1) << ", " << sz << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> ");
+ if( Stencil::D == uint_t(3) )
+ {
+ WALBERLA_ASSERT( !math::isnan( sf->get( sx, sy, sz + cell_idx_t(1), idx ) ), sx << ", " << sy << ", " << sz + cell_idx_t(1) << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> ");
+ WALBERLA_ASSERT( !math::isnan( sf->get( sx + cell_idx_t(1), sy, sz + cell_idx_t(1), idx ) ), sx + cell_idx_t(1) << ", " << sy << ", " << sz + cell_idx_t(1) << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> ");
+ WALBERLA_ASSERT( !math::isnan( sf->get( sx, sy + cell_idx_t(1), sz + cell_idx_t(1), idx ) ), sx << ", " << sy + cell_idx_t(1) << ", " << sz + cell_idx_t(1) << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> ");
+ WALBERLA_ASSERT( !math::isnan( sf->get( sx + cell_idx_t(1), sy + cell_idx_t(1), sz + cell_idx_t(1), idx ) ), sx + cell_idx_t(1) << ", " << sy + cell_idx_t(1) << ", " << sz + cell_idx_t(1) << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> ");
+ }
+ }
+#endif
+
+ typename PdfField_T::value_type value = sf->get( sx, sy, sz, idx );
+ value += sf->get( sx + cell_idx_t(1), sy, sz, idx );
+ value += sf->get( sx, sy + cell_idx_t(1), sz, idx );
+ value += sf->get( sx + cell_idx_t(1), sy + cell_idx_t(1), sz, idx );
+ if( Stencil::D == uint_t(3) )
+ {
+ value += sf->get( sx, sy, sz + cell_idx_t(1), idx );
+ value += sf->get( sx + cell_idx_t(1), sy, sz + cell_idx_t(1), idx );
+ value += sf->get( sx, sy + cell_idx_t(1), sz + cell_idx_t(1), idx );
+ value += sf->get( sx + cell_idx_t(1), sy + cell_idx_t(1), sz + cell_idx_t(1), idx );
+ }
+
+ rf->get( rx, ry, rz, idx ) = factor * value;
+ }
+ sx += cell_idx_t(2);
+ }
+ sy += cell_idx_t(2);
+ WALBERLA_ASSERT_EQUAL( sx, packingInterval.xMax() + cell_idx_t(1) );
+ }
+ sz += cell_idx_t(2);
+ WALBERLA_ASSERT_EQUAL( sy, packingInterval.yMax() + cell_idx_t(1) );
+ }
+ if( Stencil::D == uint_t(3) )
+ {
+ WALBERLA_ASSERT_EQUAL( sz, packingInterval.zMax() + cell_idx_t(1) );
+ }
+ else
+ {
+ WALBERLA_ASSERT_EQUAL( sz, packingInterval.zMax() + cell_idx_t(2) );
+ }
+}
+
+
+
+///////////////////////////////////////////////////////////////////////
+// Helper functions for determining packing/unpacking cell intervals //
+///////////////////////////////////////////////////////////////////////
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+inline CellInterval PdfFieldPackInfo< LatticeModel_T >::equalLevelPackInterval( stencil::Direction dir, const CellInterval & cellBB,
+ const uint_t numberOfLayers )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+inline CellInterval PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::equalLevelPackInterval( stencil::Direction dir, const CellInterval & cellBB,
+ const uint_t numberOfLayers )
+#endif
+{
+ CellInterval interval = equalLevelUnpackInterval( dir, cellBB, numberOfLayers );
+
+ for( uint_t i = 0; i != Stencil::D; ++i )
+ {
+ const auto offset = cell_idx_c( stencil::c[i][dir] ) * cell_idx_c( numberOfLayers );
+ interval.min()[i] -= offset;
+ interval.max()[i] -= offset;
+ }
+
+ return interval;
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+CellInterval PdfFieldPackInfo< LatticeModel_T >::equalLevelUnpackInterval( stencil::Direction dir, const CellInterval & cellBB,
+ const uint_t numberOfLayers )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+CellInterval PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::equalLevelUnpackInterval( stencil::Direction dir, const CellInterval & cellBB,
+ const uint_t numberOfLayers )
+#endif
+{
+ CellInterval interval( cellBB );
+ interval.expand( cell_idx_c(numberOfLayers) );
+
+ for( uint_t i = 0; i != 3; ++i )
+ {
+ const auto c = cell_idx_c( stencil::c[i][dir] );
+
+ if( c == -1 ) interval.max()[i] = interval.min()[i] + cell_idx_c( numberOfLayers - 1 );
+ else if( c == 1 ) interval.min()[i] = interval.max()[i] - cell_idx_c( numberOfLayers - 1 );
+ else
+ {
+ WALBERLA_ASSERT_EQUAL( c, cell_idx_t(0) );
+ interval.min()[i] += cell_idx_c( numberOfLayers );
+ interval.max()[i] -= cell_idx_c( numberOfLayers );
+ }
+ }
+
+ return interval;
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+inline bool PdfFieldPackInfo< LatticeModel_T >::equalLevelFaceIntervalSplitable( const CellInterval & interval, stencil::Direction dir ) const
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+inline bool PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::equalLevelFaceIntervalSplitable( const CellInterval & interval, stencil::Direction dir ) const
+#endif
+{
+ if( stencil::cx[dir] != 0 )
+ {
+ WALBERLA_ASSERT_EQUAL( stencil::cy[dir], 0 );
+ WALBERLA_ASSERT_EQUAL( stencil::cz[dir], 0 );
+
+ WALBERLA_ASSERT_EQUAL( interval.xSize(), equalLevelCells_ );
+
+ return interval.ySize() > uint_t(2) && interval.zSize() > uint_t(2);
+ }
+ else if( stencil::cy[dir] != 0 )
+ {
+ WALBERLA_ASSERT_EQUAL( stencil::cx[dir], 0 );
+ WALBERLA_ASSERT_EQUAL( stencil::cz[dir], 0 );
+
+ WALBERLA_ASSERT_EQUAL( interval.ySize(), equalLevelCells_ );
+
+ return interval.xSize() > uint_t(2) && interval.zSize() > uint_t(2);
+ }
+
+ WALBERLA_ASSERT_UNEQUAL( stencil::cz[dir], 0 );
+ WALBERLA_ASSERT_EQUAL( stencil::cx[dir], 0 );
+ WALBERLA_ASSERT_EQUAL( stencil::cy[dir], 0 );
+
+ WALBERLA_ASSERT_EQUAL( interval.zSize(), equalLevelCells_ );
+
+ return interval.xSize() > uint_t(2) && interval.ySize() > uint_t(2);
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+inline std::vector< CellInterval > PdfFieldPackInfo< LatticeModel_T >::splitEqualLevelFaceInterval( const CellInterval & interval,
+ stencil::Direction dir ) const // dir: from sender to receiver
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+inline std::vector< CellInterval > PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::splitEqualLevelFaceInterval( const CellInterval & interval,
+ stencil::Direction dir ) const // dir: from sender to receiver
+#endif
+{
+ std::vector< CellInterval > intervals;
+
+ if( stencil::cx[dir] != 0 )
+ {
+ WALBERLA_ASSERT_EQUAL( stencil::cy[dir], 0 );
+ WALBERLA_ASSERT_EQUAL( stencil::cz[dir], 0 );
+
+ WALBERLA_ASSERT_EQUAL( interval.xSize(), equalLevelCells_ );
+
+ intervals.push_back( CellInterval( interval.xMin(), interval.yMin(), interval.zMax(),
+ interval.xMax(), interval.yMax(), interval.zMax() ) );
+ intervals.push_back( CellInterval( interval.xMin(), interval.yMin(), interval.zMin() + cell_idx_t(1),
+ interval.xMax(), interval.yMin(), interval.zMax() - cell_idx_t(1) ) );
+ intervals.push_back( CellInterval( interval.xMin(), interval.yMax(), interval.zMin() + cell_idx_t(1),
+ interval.xMax(), interval.yMax(), interval.zMax() - cell_idx_t(1) ) );
+ intervals.push_back( CellInterval( interval.xMin(), interval.yMin(), interval.zMin(),
+ interval.xMax(), interval.yMax(), interval.zMin() ) );
+
+ const cell_idx_t x = ( stencil::cx[dir] > 0 ) ? interval.xMax() : interval.xMin();
+ intervals.push_back( CellInterval( x, interval.yMin() + cell_idx_t(1), interval.zMin() + cell_idx_t(1),
+ x, interval.yMax() - cell_idx_t(1), interval.zMax() - cell_idx_t(1) ) );
+ }
+ else if( stencil::cy[dir] != 0 )
+ {
+ WALBERLA_ASSERT_EQUAL( stencil::cx[dir], 0 );
+ WALBERLA_ASSERT_EQUAL( stencil::cz[dir], 0 );
+
+ WALBERLA_ASSERT_EQUAL( interval.ySize(), equalLevelCells_ );
+
+ intervals.push_back( CellInterval( interval.xMin(), interval.yMin(), interval.zMax(),
+ interval.xMax(), interval.yMax(), interval.zMax() ) );
+ intervals.push_back( CellInterval( interval.xMin(), interval.yMin(), interval.zMin() + cell_idx_t(1),
+ interval.xMin(), interval.yMax(), interval.zMax() - cell_idx_t(1) ) );
+ intervals.push_back( CellInterval( interval.xMax(), interval.yMin(), interval.zMin() + cell_idx_t(1),
+ interval.xMax(), interval.yMax(), interval.zMax() - cell_idx_t(1) ) );
+ intervals.push_back( CellInterval( interval.xMin(), interval.yMin(), interval.zMin(),
+ interval.xMax(), interval.yMax(), interval.zMin() ) );
+
+ const cell_idx_t y = ( stencil::cy[dir] > 0 ) ? interval.yMax() : interval.yMin();
+ intervals.push_back( CellInterval( interval.xMin() + cell_idx_t(1), y, interval.zMin() + cell_idx_t(1),
+ interval.xMax() - cell_idx_t(1), y, interval.zMax() - cell_idx_t(1) ) );
+ }
+ else
+ {
+ WALBERLA_ASSERT_UNEQUAL( stencil::cz[dir], 0 );
+ WALBERLA_ASSERT_EQUAL( stencil::cx[dir], 0 );
+ WALBERLA_ASSERT_EQUAL( stencil::cy[dir], 0 );
+
+ WALBERLA_ASSERT_EQUAL( interval.zSize(), equalLevelCells_ );
+
+ intervals.push_back( CellInterval( interval.xMin(), interval.yMax(), interval.zMin(),
+ interval.xMax(), interval.yMax(), interval.zMax() ) );
+ intervals.push_back( CellInterval( interval.xMin(), interval.yMin() + cell_idx_t(1), interval.zMin(),
+ interval.xMin(), interval.yMax() - cell_idx_t(1), interval.zMax() ) );
+ intervals.push_back( CellInterval( interval.xMax(), interval.yMin() + cell_idx_t(1), interval.zMin(),
+ interval.xMax(), interval.yMax() - cell_idx_t(1), interval.zMax() ) );
+ intervals.push_back( CellInterval( interval.xMin(), interval.yMin(), interval.zMin(),
+ interval.xMax(), interval.yMin(), interval.zMax() ) );
+
+ const cell_idx_t z = ( stencil::cz[dir] > 0 ) ? interval.zMax() : interval.zMin();
+ intervals.push_back( CellInterval( interval.xMin() + cell_idx_t(1), interval.yMin() + cell_idx_t(1), z,
+ interval.xMax() - cell_idx_t(1), interval.yMax() - cell_idx_t(1), z ) );
+ }
+
+ return intervals;
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+inline Vector3< cell_idx_t > PdfFieldPackInfo< LatticeModel_T >::getNeighborShift( const BlockID & smallBlock, stencil::Direction dir ) // dir: direction from big to small block
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+inline Vector3< cell_idx_t > PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::getNeighborShift( const BlockID & smallBlock,
+ stencil::Direction dir ) // dir: direction from big to small block
+#endif
+{
+ Vector3< cell_idx_t > shift;
+
+ uint_t branchId = smallBlock.getBranchId();
+
+ shift[0] = ( stencil::cx[dir] == 0 ) ? ( ( ( branchId & uint_t(1) ) == uint_t(0) ) ? cell_idx_t(-1) : cell_idx_t(1) ) : cell_idx_t(0);
+ shift[1] = ( stencil::cy[dir] == 0 ) ? ( ( ( branchId & uint_t(2) ) == uint_t(0) ) ? cell_idx_t(-1) : cell_idx_t(1) ) : cell_idx_t(0);
+ shift[2] = ( Stencil::D == uint_t(3) ) ?
+ ( ( stencil::cz[dir] == 0 ) ? ( ( ( branchId & uint_t(4) ) == uint_t(0) ) ? cell_idx_t(-1) : cell_idx_t(1) ) : cell_idx_t(0) ) : cell_idx_t(0);
+
+ return shift;
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+CellInterval PdfFieldPackInfo< LatticeModel_T >::coarseToFinePackInterval( stencil::Direction dir, const CellInterval & cellBB,
+ const BlockID & smallBlock )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+CellInterval PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::coarseToFinePackInterval( stencil::Direction dir, const CellInterval & cellBB,
+ const BlockID & smallBlock )
+#endif
+{
+ CellInterval interval = equalLevelPackInterval( dir, cellBB, uint_t(2) );
+ //CellInterval interval = equalLevelPackInterval( dir, cellBB, uint_t(1) );
+
+ Vector3< cell_idx_t > shift = getNeighborShift( smallBlock, dir );
+
+ WALBERLA_ASSERT( divisibleByTwo( cellBB ) );
+
+ for( uint_t i = 0; i != Stencil::D; ++i )
+ {
+ if( shift[i] == cell_idx_t(-1) )
+ interval.max()[i] = interval.min()[i] + cell_idx_c( cellBB.size(i) / uint_t(2) ) + cell_idx_t(1);
+ //interval.max()[i] = interval.min()[i] + cell_idx_c( cellBB.size(i) / uint_t(2) );
+ if( shift[i] == cell_idx_t( 1) )
+ interval.min()[i] = interval.max()[i] - cell_idx_c( cellBB.size(i) / uint_t(2) ) - cell_idx_t(1);
+ //interval.min()[i] = interval.max()[i] - cell_idx_c( cellBB.size(i) / uint_t(2) );
+ }
+
+ WALBERLA_ASSERT( cellBB.contains( interval ) );
+
+ return interval;
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+CellInterval PdfFieldPackInfo< LatticeModel_T >::coarseToFineUnpackInterval( stencil::Direction dir, const CellInterval & cellBB,
+ const BlockID & smallBlock )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+CellInterval PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::coarseToFineUnpackInterval( stencil::Direction dir, const CellInterval & cellBB,
+ const BlockID & smallBlock )
+#endif
+{
+ CellInterval interval = equalLevelUnpackInterval( dir, cellBB, uint_t(4) );
+ //CellInterval interval = equalLevelUnpackInterval( dir, cellBB, uint_t(2) );
+
+ Vector3< cell_idx_t > shift = getNeighborShift( smallBlock, dir );
+
+ for( uint_t i = 0; i != Stencil::D; ++i )
+ {
+ if( shift[i] == cell_idx_t(-1) )
+ interval.max()[i] += cell_idx_t(4);
+ //interval.max()[i] += cell_idx_t(2);
+ if( shift[i] == cell_idx_t( 1) )
+ interval.min()[i] -= cell_idx_t(4);
+ //interval.min()[i] -= cell_idx_t(2);
+ }
+
+#ifndef NDEBUG
+ CellInterval expandedCellBB( cellBB );
+ expandedCellBB.expand( cell_idx_t(4) );
+ WALBERLA_ASSERT( expandedCellBB.contains( interval ) );
+#endif
+
+ return interval;
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+inline CellInterval PdfFieldPackInfo< LatticeModel_T >::fineToCoarsePackInterval( stencil::Direction dir, const CellInterval & cellBB )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+inline CellInterval PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::fineToCoarsePackInterval( stencil::Direction dir, const CellInterval & cellBB )
+#endif
+{
+ return equalLevelUnpackInterval( dir, cellBB, uint_t(2) );
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+CellInterval PdfFieldPackInfo< LatticeModel_T >::fineToCoarseUnpackInterval( stencil::Direction dir, const CellInterval & cellBB,
+ const BlockID & smallBlock )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+CellInterval PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::fineToCoarseUnpackInterval( stencil::Direction dir, const CellInterval & cellBB,
+ const BlockID & smallBlock )
+#endif
+{
+ CellInterval interval = equalLevelPackInterval( dir, cellBB, uint_t(1) );
+ Vector3< cell_idx_t > shift = getNeighborShift( smallBlock, dir );
+
+ WALBERLA_ASSERT( divisibleByTwo( cellBB ) );
+
+ for( uint_t i = 0; i != Stencil::D; ++i )
+ {
+ if( shift[i] == cell_idx_t(-1) )
+ interval.max()[i] = interval.min()[i] + cell_idx_c( cellBB.size(i) / uint_t(2) ) - cell_idx_t(1);
+ if( shift[i] == cell_idx_t( 1) )
+ interval.min()[i] = interval.max()[i] - cell_idx_c( cellBB.size(i) / uint_t(2) ) + cell_idx_t(1);
+ }
+
+ WALBERLA_ASSERT( cellBB.contains( interval ) );
+
+ return interval;
+}
+
+
+
+//////////////////////////////
+// General helper functions //
+//////////////////////////////
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+inline uint_t PdfFieldPackInfo< LatticeModel_T >::equalLevelCells()
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+inline uint_t PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::equalLevelCells()
+#endif
+{
+ if( Stencil::containsDir( stencil::TNE ) || Stencil::containsDir( stencil::TNW ) || Stencil::containsDir( stencil::TSE ) ||
+ Stencil::containsDir( stencil::TSW ) || Stencil::containsDir( stencil::BNE ) || Stencil::containsDir( stencil::BNW ) ||
+ Stencil::containsDir( stencil::BSE ) || Stencil::containsDir( stencil::BSW ) )
+ return uint_t(3);
+
+ return uint_t(2);
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+inline bool PdfFieldPackInfo< LatticeModel_T >::isFaceDirection( stencil::Direction dir )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+inline bool PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::isFaceDirection( stencil::Direction dir )
+#endif
+{
+ return ( dir == stencil::N ) || ( dir == stencil::S ) || ( dir == stencil::W ) ||
+ ( dir == stencil::E ) || ( dir == stencil::T ) || ( dir == stencil::B );
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+inline bool PdfFieldPackInfo< LatticeModel_T >::isEdgeDirection( stencil::Direction dir )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+inline bool PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::isEdgeDirection( stencil::Direction dir )
+#endif
+{
+ return ( dir == stencil::NW ) || ( dir == stencil::NE ) || ( dir == stencil::SW ) || ( dir == stencil::SE ) ||
+ ( dir == stencil::TN ) || ( dir == stencil::TS ) || ( dir == stencil::TW ) || ( dir == stencil::TE ) ||
+ ( dir == stencil::BN ) || ( dir == stencil::BS ) || ( dir == stencil::BW ) || ( dir == stencil::BE );
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+inline bool PdfFieldPackInfo< LatticeModel_T >::isCornerDirection( stencil::Direction dir )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+inline bool PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::isCornerDirection( stencil::Direction dir )
+#endif
+{
+ return ( dir == stencil::TNE ) || ( dir == stencil::TNW ) || ( dir == stencil::TSE ) || ( dir == stencil::TSW ) ||
+ ( dir == stencil::BNE ) || ( dir == stencil::BNW ) || ( dir == stencil::BSE ) || ( dir == stencil::BSW );
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+inline bool PdfFieldPackInfo< LatticeModel_T >::blocksConnectedByFaces( const Block * block, const BlockID & neighbor )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+inline bool PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::blocksConnectedByFaces( const Block * block, const BlockID & neighbor )
+#endif
+{
+ const uint_t face[] = { uint_t(4), uint_t(10), uint_t(12), uint_t(13), uint_t(15), uint_t(21) };
+ for( int i = 0; i != 6; ++i )
+ {
+ for( uint_t n = 0; n != block->getNeighborhoodSectionSize( face[i] ); ++n )
+ if( block->getNeighborId( face[i], n ) == neighbor )
+ return true;
+ }
+ return false;
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+inline bool PdfFieldPackInfo< LatticeModel_T >::blocksConnectedByEdges( const Block * block, const BlockID & neighbor )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+inline bool PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::blocksConnectedByEdges( const Block * block, const BlockID & neighbor )
+#endif
+{
+ const uint_t face[] = { uint_t( 1), uint_t( 3), uint_t( 5), uint_t( 7), uint_t( 9), uint_t(11),
+ uint_t(14), uint_t(16), uint_t(18), uint_t(20), uint_t(22), uint_t(24) };
+
+ for( int i = 0; i != 12; ++i )
+ {
+ for( uint_t n = 0; n != block->getNeighborhoodSectionSize( face[i] ); ++n )
+ if( block->getNeighborId( face[i], n ) == neighbor )
+ return true;
+ }
+ return false;
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+inline bool PdfFieldPackInfo< LatticeModel_T >::divisibleByTwo( const CellInterval & cellBB )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+inline bool PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::divisibleByTwo( const CellInterval & cellBB )
+#endif
+{
+ return ( ( cellBB.xSize() & uint_t(1) ) == uint_t(0) ) &&
+ ( ( cellBB.ySize() & uint_t(1) ) == uint_t(0) ) &&
+ ( ( Stencil::D == uint_t(2) && cellBB.zSize() == uint_t(1) ) || ( Stencil::D == uint_t(3) && ( cellBB.zSize() & uint_t(1) ) == uint_t(0) ) );
+}
+
+
+
+#ifdef NDEBUG
+template< typename LatticeModel_T >
+bool PdfFieldPackInfo< LatticeModel_T >::coarserNeighborExistsInVicinity( const Block * block, stencil::Direction dir )
+#else
+template< typename LatticeModel_T, typename BoundaryHandling_T >
+bool PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::coarserNeighborExistsInVicinity( const Block * block, stencil::Direction dir )
+#endif
+{
+ if( block->getLevel() == uint_t(0) )
+ return false;
+
+ Vector3<int> min( -1 );
+ Vector3<int> max( 1 );
+
+ for( uint_t i = 0; i != 3; ++i )
+ {
+ if( stencil::c[i][dir] == -1 ) max[i] = 0;
+ if( stencil::c[i][dir] == 1 ) min[i] = 0;
+ }
+ if( LatticeModel_T::Stencil::D == uint_t(2) )
+ min[2] = max[2] = 0;
+
+ for( int z = min[2]; z <= max[2]; ++z ) {
+ for( int y = min[1]; y <= max[1]; ++y ) {
+ for( int x = min[0]; x <= max[0]; ++x )
+ {
+ if( x == 0 && y == 0 && z == 0 )
+ continue;
+ if( block->neighborhoodSectionHasLargerBlock( blockforest::getBlockNeighborhoodSectionIndex(x,y,z) ) )
+ return true;
+ }
+ }
+ }
+
+ return false;
+}
+
+
+
+} // namespace refinement
+} // namespace lbm
+} // namespace walberla
diff --git a/src/lbm/refinement_rebase/PdfFieldSyncPackInfo.h b/src/lbm/refinement_rebase/PdfFieldSyncPackInfo.h
new file mode 100644
index 000000000..0b68676b7
--- /dev/null
+++ b/src/lbm/refinement_rebase/PdfFieldSyncPackInfo.h
@@ -0,0 +1,44 @@
+//======================================================================================================================
+//
+// 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 PdfFieldSyncPackInfo.h
+//! \ingroup lbm
+//! \author Florian Schornbaum <florian.schornbaum@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "lbm/field/PdfField.h"
+#include "field/refinement/PackInfo.h"
+
+
+namespace walberla {
+namespace lbm {
+namespace refinement {
+
+template <typename LatticeModel_T>
+class PdfFieldSyncPackInfo : public field::refinement::PackInfo< PdfField<LatticeModel_T>, typename LatticeModel_T::Stencil >
+{
+public:
+ PdfFieldSyncPackInfo( const BlockDataID & fieldId )
+ : field::refinement::PackInfo< PdfField<LatticeModel_T>, typename LatticeModel_T::Stencil >( fieldId )
+ {
+ }
+};
+
+} // namespace refinement
+} // namespace lbm
+} // namespace walberla
diff --git a/src/lbm/refinement_rebase/RefinementFunctorWrapper.h b/src/lbm/refinement_rebase/RefinementFunctorWrapper.h
new file mode 100644
index 000000000..9ea0b99e1
--- /dev/null
+++ b/src/lbm/refinement_rebase/RefinementFunctorWrapper.h
@@ -0,0 +1,67 @@
+//======================================================================================================================
+//
+// 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 RefinementFunctorWrapper.h
+//! \ingroup lbm
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "domain_decomposition/BlockStorage.h"
+
+#include <functional>
+
+namespace walberla {
+namespace lbm {
+namespace refinement {
+
+class FunctorWrapper {
+public:
+ FunctorWrapper( const std::function<void()> & fct)
+ : fct_(fct) {
+ }
+
+ void operator()(const uint_t /*level*/, const uint_t /*executionCounter*/) {
+ fct_();
+ }
+
+private:
+ std::function<void(void)> fct_;
+};
+
+class SweepAsFunctorWrapper {
+public:
+ SweepAsFunctorWrapper( const std::function<void(IBlock * )> & fct,
+ const shared_ptr <StructuredBlockStorage> & blockStorage )
+ : fct_(fct), blockStorage_(blockStorage) {
+ }
+
+ void operator()(const uint_t level, const uint_t /*executionCounter*/) {
+ for (auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt) {
+ if (blockStorage_->getLevel(*blockIt) != level) continue;
+ fct_(blockIt.get());
+ }
+ }
+
+private:
+ std::function<void(IBlock * )> fct_;
+ shared_ptr <StructuredBlockStorage> blockStorage_;
+};
+
+} // namespace refinement
+} // namespace lbm
+} // namespace walberla
diff --git a/src/lbm/refinement_rebase/TimeStep.h b/src/lbm/refinement_rebase/TimeStep.h
new file mode 100644
index 000000000..1be6a470d
--- /dev/null
+++ b/src/lbm/refinement_rebase/TimeStep.h
@@ -0,0 +1,1751 @@
+//======================================================================================================================
+//
+// 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 TimeStep.h
+//! \ingroup lbm
+//! \author Florian Schornbaum <florian.schornbaum@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "EqualLevelBorderStreamCorrection.h"
+#include "LinearExplosion.h"
+#include "PdfFieldPackInfo.h"
+#include "TimeStepPdfPackInfo.h"
+
+#include "blockforest/StructuredBlockForest.h"
+#include "blockforest/communication/NonUniformBufferedScheme.h"
+
+#include "core/debug/CheckFunctions.h"
+#include "core/debug/Debug.h"
+#include "core/math/Uint.h"
+#include "core/selectable/IsSetSelected.h"
+#include "core/timing/TimingPool.h"
+
+#include "lbm/lattice_model/NeighborsStencil.h"
+
+#include "stencil/D2Q9.h"
+#include "stencil/D3Q27.h"
+
+
+
+namespace walberla {
+namespace lbm {
+namespace refinement {
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+class TimeStep
+{
+public:
+
+ typedef typename NeighborsStencil<LatticeModel_T>::type CommunicationStencil_T;
+
+ typedef std::function< void ( const uint_t, const uint_t ) > VoidFunction; // parameters: level, execution count
+ typedef std::function< void ( IBlock *, const uint_t, const uint_t ) > BlockFunction; // parameters: block, level, execution count
+
+private:
+
+ class VoidFunctionWrappper
+ {
+ public:
+ VoidFunctionWrappper( std::vector< std::pair< VoidFunction, std::string > > & functions, const uint_t index ) :
+ functions_( functions ), index_( index ) {}
+ void operator()( const uint_t level, const uint_t executionCount )
+ {
+ functions_[index_].first( level, executionCount );
+ }
+ private:
+ std::vector< std::pair< VoidFunction, std::string > > & functions_;
+ uint_t index_;
+ };
+
+ class BlockFunctionWrappper
+ {
+ public:
+ BlockFunctionWrappper( std::vector< std::pair< BlockFunction, std::string > > & functions, const uint_t index ) :
+ functions_( functions ), index_( index ) {}
+ void operator()( IBlock * block, const uint_t level, const uint_t executionCount )
+ {
+ functions_[index_].first( block, level, executionCount );
+ }
+ private:
+ std::vector< std::pair< BlockFunction, std::string > > & functions_;
+ uint_t index_;
+ };
+
+ template< typename F >
+ class SharedVoidFunctor
+ {
+ public:
+ SharedVoidFunctor( const shared_ptr<F> & functorPtr ) : functorPtr_( functorPtr ) { }
+ void operator()( const uint_t level, const uint_t executionCount ){ (*functorPtr_)( level, executionCount ); }
+ private:
+ shared_ptr<F> functorPtr_;
+ };
+
+ template< typename F >
+ class SharedBlockFunctor
+ {
+ public:
+ SharedBlockFunctor( const shared_ptr<F> & functorPtr ) : functorPtr_( functorPtr ) { }
+ void operator()( IBlock * block, const uint_t level, const uint_t executionCount ){ (*functorPtr_)( block, level, executionCount ); }
+ private:
+ shared_ptr<F> functorPtr_;
+ };
+
+public:
+
+ typedef typename LatticeModel_T::Stencil Stencil;
+
+ static const uint_t StreamIncludedGhostLayers = 2;
+
+ // constructors
+
+ TimeStep( weak_ptr< StructuredBlockForest > blocks, shared_ptr< Sweep_T > & sweep,
+ const BlockDataID & pdfFieldId, const BlockDataID & boundaryHandlingId,
+ const Set<SUID> & requiredBlockSelectors = Set<SUID>::emptySet(),
+ const Set<SUID> & incompatibleBlockSelectors = Set<SUID>::emptySet() );
+
+ TimeStep( weak_ptr< StructuredBlockForest > blocks, shared_ptr< Sweep_T > & sweep,
+ const BlockDataID & pdfFieldId, const BlockDataID & boundaryHandlingId,
+ const shared_ptr< TimeStepPdfPackInfo > & pdfPackInfo,
+ const Set<SUID> & requiredBlockSelectors = Set<SUID>::emptySet(),
+ const Set<SUID> & incompatibleBlockSelectors = Set<SUID>::emptySet() );
+
+ // member functions
+
+ bool asynchronousCommunicationIsUsed() const { return asynchronousCommunication_; }
+ void asynchronousCommunication( const bool value = true ) { asynchronousCommunication_ = value; }
+
+ bool optimizedCommunicationIsUsed() const { return optimizedCommunication_; }
+ void optimizeCommunication( const bool value = true )
+ {
+ optimizedCommunication_ = value;
+ pdfPackInfo_->optimizeEqualLevelCommunication( optimizedCommunication_ );
+ pdfPackInfo_->optimizeForLinearExplosion( optimizedCommunication_ && performLinearExplosion_ );
+ }
+
+ bool equalLevelBorderStreamCorrectionIsPerformed() const { return performEqualLevelBorderStreamCorrection_; }
+ void performEqualLevelBorderStreamCorrection( const bool value = true )
+ {
+ performEqualLevelBorderStreamCorrection_ = value;
+ }
+
+ bool linearExplosionIsPerformed() const { return performLinearExplosion_; }
+ void performLinearExplosion( const bool value = true )
+ {
+ performLinearExplosion_ = value;
+ pdfPackInfo_->optimizeForLinearExplosion( optimizedCommunication_ && performLinearExplosion_ );
+ }
+
+ void deactivateTiming() { timing_ = false; }
+ void enableTiming( const shared_ptr<WcTimingPool> & timingPool, const shared_ptr<WcTimingPool> & levelwiseTimingPool )
+ {
+ timing_ = true;
+ timingPool_ = timingPool;
+ levelwiseTimingPool_ = levelwiseTimingPool;
+ refresh( postCollideVoidFunctions_.size() );
+ }
+
+ void enableTiming()
+ {
+ timing_ = true;
+ if( !timingPool_ )
+ timingPool_ = make_shared<WcTimingPool>();
+ if( !levelwiseTimingPool_ )
+ levelwiseTimingPool_ = make_shared<WcTimingPool>();
+ refresh( postCollideVoidFunctions_.size() );
+ }
+
+ const shared_ptr<WcTimingPool> & getTimingPool() const { return timingPool_; }
+ const shared_ptr<WcTimingPool> & getLevelWiseTimingPool() const { return levelwiseTimingPool_; }
+
+ void operator()()
+ {
+ auto blocks = blocks_.lock();
+ WALBERLA_CHECK_NOT_NULLPTR( blocks, "Trying to access 'TimeStep' (refinement) for a block storage object that doesn't exist anymore" );
+ if( blocks->getNumberOfLevels() > postCollideVoidFunctions_.size() )
+ refresh( blocks->getNumberOfLevels() );
+ recursiveStep( uint_t(0), uint_t(0) );
+ }
+
+ void addPackInfo( const typename blockforest::communication::NonUniformBufferedScheme< CommunicationStencil_T >::PackInfo & packInfo )
+ {
+ communication_.addPackInfo( packInfo );
+ }
+
+ inline void addPostCollideVoidFunction ( const VoidFunction & function, const std::string & identifier );
+ inline void addPostCollideBlockFunction( const BlockFunction & function, const std::string & identifier );
+ inline void addPostCollideVoidFunction ( const VoidFunction & function, const std::string & identifier, const uint_t level );
+ inline void addPostCollideBlockFunction( const BlockFunction & function, const std::string & identifier, const uint_t level );
+
+ inline void addPostBoundaryHandlingVoidFunction ( const VoidFunction & function, const std::string & identifier );
+ inline void addPostBoundaryHandlingBlockFunction( const BlockFunction & function, const std::string & identifier );
+ inline void addPostBoundaryHandlingVoidFunction ( const VoidFunction & function, const std::string & identifier, const uint_t level );
+ inline void addPostBoundaryHandlingBlockFunction( const BlockFunction & function, const std::string & identifier, const uint_t level );
+
+ inline void addPostStreamVoidFunction ( const VoidFunction & function, const std::string & identifier );
+ inline void addPostStreamBlockFunction( const BlockFunction & function, const std::string & identifier );
+ inline void addPostStreamVoidFunction ( const VoidFunction & function, const std::string & identifier, const uint_t level );
+ inline void addPostStreamBlockFunction( const BlockFunction & function, const std::string & identifier, const uint_t level );
+
+ template< typename F >
+ inline void addPostCollideVoidFunction ( const shared_ptr<F> & functorPtr, const std::string & identifier );
+ template< typename F >
+ inline void addPostCollideBlockFunction( const shared_ptr<F> & functorPtr, const std::string & identifier );
+ template< typename F >
+ inline void addPostCollideVoidFunction ( const shared_ptr<F> & functorPtr, const std::string & identifier, const uint_t level );
+ template< typename F >
+ inline void addPostCollideBlockFunction( const shared_ptr<F> & functorPtr, const std::string & identifier, const uint_t level );
+
+ template< typename F >
+ inline void addPostBoundaryHandlingVoidFunction ( const shared_ptr<F> & functorPtr, const std::string & identifier );
+ template< typename F >
+ inline void addPostBoundaryHandlingBlockFunction( const shared_ptr<F> & functorPtr, const std::string & identifier );
+ template< typename F >
+ inline void addPostBoundaryHandlingVoidFunction ( const shared_ptr<F> & functorPtr, const std::string & identifier, const uint_t level );
+ template< typename F >
+ inline void addPostBoundaryHandlingBlockFunction( const shared_ptr<F> & functorPtr, const std::string & identifier, const uint_t level );
+
+ template< typename F >
+ inline void addPostStreamVoidFunction ( const shared_ptr<F> & functorPtr, const std::string & identifier );
+ template< typename F >
+ inline void addPostStreamBlockFunction( const shared_ptr<F> & functorPtr, const std::string & identifier );
+ template< typename F >
+ inline void addPostStreamVoidFunction ( const shared_ptr<F> & functorPtr, const std::string & identifier, const uint_t level );
+ template< typename F >
+ inline void addPostStreamBlockFunction( const shared_ptr<F> & functorPtr, const std::string & identifier, const uint_t level );
+
+private:
+
+ void init( const BlockDataID & pdfFieldId, const BlockDataID & boundaryHandlingId );
+ void consistencyChecks( const BlockDataID & pdfFieldId, const BlockDataID & boundaryHandlingId ) const;
+ void refresh( const uint_t levels );
+
+ std::string getTimingPoolString( const std::string & name, const std::string & suffix = std::string("") ) const;
+ std::string getLevelwiseTimingPoolString( const std::string & name, const uint_t level, const std::string & suffix = std::string("") ) const;
+
+ void createTimers( const uint_t levels );
+
+ template< typename Function >
+ void addFunction( std::vector< std::vector< std::pair< Function, std::string > > > & functions,
+ const Function & function, const std::string & identifier );
+
+ template< typename Function >
+ void addFunction( std::vector< std::vector< std::pair< Function, std::string > > > & functions,
+ const Function & function, const std::string & identifier, const uint_t level );
+
+ std::vector< Block * > selectedBlocks( const uint_t level ) const;
+
+ inline void startTiming( const std::string & name, const uint_t level, const std::string & suffix = std::string("") );
+ inline void stopTiming( const std::string & name, const uint_t level, const std::string & suffix = std::string("") );
+
+ void collide( std::vector< Block * > & blocks, const uint_t level, const uint_t executionCount );
+ void stream( std::vector< Block * > & blocks, const uint_t level, const uint_t executionCount );
+ void finishStream( std::vector< Block * > & blocks, const uint_t level, const uint_t executionCount );
+ void streamCollide( std::vector< Block * > & blocks, const uint_t level, const uint_t executionCount );
+
+ void startCommunicationEqualLevel( const uint_t level );
+ void endCommunicationEqualLevel( const uint_t level );
+ void startCommunicationCoarseToFine( const uint_t level );
+ void endCommunicationCoarseToFine( const uint_t level );
+ void startCommunicationFineToCoarse( const uint_t level );
+ void endCommunicationFineToCoarse( const uint_t level );
+
+ void performLinearExplosion( std::vector< Block * > & blocks, const uint_t level );
+
+ void recursiveStep( const uint_t level, const uint_t executionCount );
+
+
+
+ weak_ptr< StructuredBlockForest > blocks_;
+
+ shared_ptr< Sweep_T > sweep_; // stream & collide
+ typename BoundaryHandling_T::BlockSweep boundarySweep_; // pre-stream boundary treatment
+ typename BoundaryHandling_T::BlockSweep boundarySweepWithLayers_; // pre-stream boundary treatment (including ghost layers)
+
+ bool asynchronousCommunication_;
+ bool optimizedCommunication_;
+
+ shared_ptr< TimeStepPdfPackInfo > pdfPackInfo_;
+
+ blockforest::communication::NonUniformBufferedScheme< CommunicationStencil_T > communication_;
+
+ bool performEqualLevelBorderStreamCorrection_;
+ EqualLevelBorderStreamCorrection< LatticeModel_T > equalLevelBorderStreamCorrection_;
+
+ bool performLinearExplosion_;
+ LinearExplosion< LatticeModel_T, BoundaryHandling_T > linearExplosion_;
+
+ bool timing_;
+ shared_ptr<WcTimingPool> timingPool_;
+ shared_ptr<WcTimingPool> levelwiseTimingPool_;
+
+ Set<SUID> requiredBlockSelectors_;
+ Set<SUID> incompatibleBlockSelectors_;
+
+ std::vector< std::pair< VoidFunction, std::string > > globalPostCollideVoidFunctions_;
+ std::vector< std::pair< BlockFunction, std::string > > globalPostCollideBlockFunctions_;
+
+ std::vector< std::pair< VoidFunction, std::string > > globalPostBoundaryHandlingVoidFunctions_;
+ std::vector< std::pair< BlockFunction, std::string > > globalPostBoundaryHandlingBlockFunctions_;
+
+ std::vector< std::pair< VoidFunction, std::string > > globalPostStreamVoidFunctions_;
+ std::vector< std::pair< BlockFunction, std::string > > globalPostStreamBlockFunctions_;
+
+ std::vector< std::vector< std::pair< VoidFunction, std::string > > > postCollideVoidFunctions_;
+ std::vector< std::vector< std::pair< BlockFunction, std::string > > > postCollideBlockFunctions_;
+
+ std::vector< std::vector< std::pair< VoidFunction, std::string > > > postBoundaryHandlingVoidFunctions_;
+ std::vector< std::vector< std::pair< BlockFunction, std::string > > > postBoundaryHandlingBlockFunctions_;
+
+ std::vector< std::vector< std::pair< VoidFunction, std::string > > > postStreamVoidFunctions_;
+ std::vector< std::vector< std::pair< BlockFunction, std::string > > > postStreamBlockFunctions_;
+
+}; // class TimeStep
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+const uint_t TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::StreamIncludedGhostLayers;
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::TimeStep( weak_ptr< StructuredBlockForest > blocks, shared_ptr< Sweep_T > & sweep,
+ const BlockDataID & pdfFieldId, const BlockDataID & boundaryHandlingId,
+ const Set<SUID> & requiredBlockSelectors /*= Set<SUID>::emptySet()*/,
+ const Set<SUID> & incompatibleBlockSelectors /*= Set<SUID>::emptySet()*/ ) :
+ blocks_( blocks ), sweep_( sweep ),
+ boundarySweep_( BoundaryHandling_T::getBlockSweep( boundaryHandlingId, uint_t(0) ) ),
+ boundarySweepWithLayers_( BoundaryHandling_T::getBlockSweep( boundaryHandlingId, StreamIncludedGhostLayers ) ),
+ asynchronousCommunication_( true ), optimizedCommunication_( true ),
+#ifdef NDEBUG
+ pdfPackInfo_( make_shared< lbm::refinement::PdfFieldPackInfo< LatticeModel_T > >( pdfFieldId, true, true ) ),
+#else
+ pdfPackInfo_( make_shared< lbm::refinement::PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T > >( pdfFieldId, boundaryHandlingId, true, true ) ),
+#endif
+ communication_( blocks, requiredBlockSelectors, incompatibleBlockSelectors ),
+ performEqualLevelBorderStreamCorrection_( true ), equalLevelBorderStreamCorrection_( pdfFieldId ),
+ performLinearExplosion_( true ), linearExplosion_( pdfFieldId, boundaryHandlingId ),
+ timing_( false ),
+ requiredBlockSelectors_( requiredBlockSelectors ), incompatibleBlockSelectors_( incompatibleBlockSelectors )
+{
+ init( pdfFieldId, boundaryHandlingId );
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::TimeStep( weak_ptr< StructuredBlockForest > blocks, shared_ptr< Sweep_T > & sweep,
+ const BlockDataID & pdfFieldId, const BlockDataID & boundaryHandlingId,
+ const shared_ptr< TimeStepPdfPackInfo > & pdfPackInfo,
+ const Set<SUID> & requiredBlockSelectors /*= Set<SUID>::emptySet()*/,
+ const Set<SUID> & incompatibleBlockSelectors /*= Set<SUID>::emptySet()*/ ) :
+ blocks_( blocks ), sweep_( sweep ),
+ boundarySweep_( BoundaryHandling_T::getBlockSweep( boundaryHandlingId, uint_t(0) ) ),
+ boundarySweepWithLayers_( BoundaryHandling_T::getBlockSweep( boundaryHandlingId, StreamIncludedGhostLayers ) ),
+ asynchronousCommunication_( true ), optimizedCommunication_( true ),
+ pdfPackInfo_( pdfPackInfo ),
+ communication_( blocks, requiredBlockSelectors, incompatibleBlockSelectors ),
+ performEqualLevelBorderStreamCorrection_( true ), equalLevelBorderStreamCorrection_( pdfFieldId ),
+ performLinearExplosion_( true ), linearExplosion_( pdfFieldId, boundaryHandlingId ),
+ timing_( false ),
+ requiredBlockSelectors_( requiredBlockSelectors ), incompatibleBlockSelectors_( incompatibleBlockSelectors )
+{
+ init( pdfFieldId, boundaryHandlingId );
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::init( const BlockDataID & pdfFieldId, const BlockDataID & boundaryHandlingId )
+{
+ communication_.addPackInfo( pdfPackInfo_ );
+ communication_.setLocalMode( blockforest::WAIT );
+ consistencyChecks( pdfFieldId, boundaryHandlingId );
+ auto blocks = blocks_.lock();
+ WALBERLA_CHECK_NOT_NULLPTR( blocks, "Trying to access 'TimeStep' (refinement) for a block storage object that doesn't exist anymore" );
+ refresh( blocks->getNumberOfLevels() );
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::consistencyChecks( const BlockDataID & pdfFieldId, const BlockDataID & boundaryHandlingId ) const
+{
+ // consistency checks / field validity checks
+
+ auto blocks = blocks_.lock();
+ WALBERLA_CHECK_NOT_NULLPTR( blocks, "Trying to access 'TimeStep' (refinement) for a block storage object that doesn't exist anymore" );
+ for( auto block = blocks->begin(); block != blocks->end(); ++block )
+ {
+ if( !selectable::isSetSelected( block->getState(), requiredBlockSelectors_, incompatibleBlockSelectors_ ) )
+ continue;
+
+ auto * pdfField = block->template getData< PdfField< LatticeModel_T > >( pdfFieldId );
+ if( pdfField == NULL )
+ {
+ WALBERLA_ABORT( "Could not get the PDF field from block " << block->getId() << ". Check if it is allocated on "
+ "the block and if the LatticeModel matches!" );
+ }
+
+ auto * boundaryHandling = block->template getData< BoundaryHandling_T >( boundaryHandlingId );
+ if( boundaryHandling == NULL )
+ {
+ WALBERLA_ABORT( "Could not get the boundary handling from block " << block->getId() << ". Check if it is "
+ "allocated on the block and if its type matches!" );
+ }
+
+ auto * flagField = boundaryHandling->getFlagField();
+ WALBERLA_ASSERT_NOT_NULLPTR( flagField );
+
+ if( LatticeModel_T::Stencil::D == uint_t(3) )
+ {
+ if( ( pdfField->xSize() & uint_t(1) ) == uint_t(1) || ( pdfField->ySize() & uint_t(1) ) == uint_t(1) ||
+ ( pdfField->zSize() & uint_t(1) ) == uint_t(1) )
+ WALBERLA_ABORT( "The x-, y-, and z-size of the PDF field on each block must be divisible by two!\n"
+ "(PDF field size on block " << block->getId() <<
+ " is " << pdfField->xSize() << " x " << pdfField->ySize() << " x " << pdfField->zSize() << ")" );
+
+ if( pdfField->xSize() < uint_t(4) || pdfField->ySize() < uint_t(4) || pdfField->zSize() < uint_t(4) )
+ WALBERLA_ABORT( "The size of the PDF field on each block must be at least 4x4x4 cells!\n"
+ "(PDF field size on block " << block->getId() <<
+ " is " << pdfField->xSize() << " x " << pdfField->ySize() << " x " << pdfField->zSize() << ")" );
+ }
+ else
+ {
+ WALBERLA_CHECK_EQUAL( LatticeModel_T::Stencil::D, uint_t(2) );
+
+ if( ( pdfField->xSize() & uint_t(1) ) == uint_t(1) || ( pdfField->ySize() & uint_t(1) ) == uint_t(1) )
+ WALBERLA_ABORT( "The x- and y-size of the PDF field on each block must be divisible by two!\n"
+ "(PDF field size on block " << block->getId() <<
+ " is " << pdfField->xSize() << " x " << pdfField->ySize() << " x " << pdfField->zSize() << ")" );
+
+ if( pdfField->xSize() < uint_t(4) || pdfField->ySize() < uint_t(4) )
+ WALBERLA_ABORT( "The size of the PDF field on each block must be at least 4x4x1 cells!\n"
+ "(PDF field size on block " << block->getId() <<
+ " is " << pdfField->xSize() << " x " << pdfField->ySize() << " x " << pdfField->zSize() << ")" );
+
+ if( pdfField->zSize() != uint_t(1) )
+ WALBERLA_ABORT( "The z-size of the PDF field on each block must be 1!\n"
+ "(PDF field size on block " << block->getId() <<
+ " is " << pdfField->xSize() << " x " << pdfField->ySize() << " x " << pdfField->zSize() << ")" );
+ }
+
+ if( pdfField->xSize() != flagField->xSize() || pdfField->ySize() != flagField->ySize() || pdfField->zSize() != flagField->zSize() )
+ WALBERLA_ABORT( "The size of the PDF field must be identical to the size of the flag field!\n"
+ "(PDF field [" << pdfField->xSize() << " x " << pdfField->ySize() << " x " << pdfField->zSize() <<
+ "] vs. flag field [" << flagField->xSize() << " x " << flagField->ySize() << " x " << flagField->zSize() << "])" );
+
+ if( pdfField->nrOfGhostLayers() < uint_t(4) )
+ WALBERLA_ABORT( "The PDF field of each block needs at least 4 ghost layers!\n"
+ "(currently only possesses " << pdfField->nrOfGhostLayers() << " on block " << block->getId() << ")" );
+
+ if( flagField->nrOfGhostLayers() < uint_t(4) )
+ WALBERLA_ABORT( "The flag field of each block needs at least 4 ghost layers!\n"
+ "(currently only possesses " << flagField->nrOfGhostLayers() << " on block " << block->getId() << ")" );
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::refresh( const uint_t levels )
+{
+ createTimers( levels );
+
+ if( levels > postCollideVoidFunctions_.size() )
+ {
+ const uint_t previousLevels = postCollideVoidFunctions_.size();
+
+ postCollideVoidFunctions_.resize( levels );
+ postCollideBlockFunctions_.resize( levels );
+
+ postBoundaryHandlingVoidFunctions_.resize( levels );
+ postBoundaryHandlingBlockFunctions_.resize( levels );
+
+ postStreamVoidFunctions_.resize( levels );
+ postStreamBlockFunctions_.resize( levels );
+
+ for( uint_t f = uint_t(0); f != globalPostCollideVoidFunctions_.size(); ++f )
+ {
+ VoidFunctionWrappper wrappedFunction( globalPostCollideVoidFunctions_, f );
+ const std::string identifier = globalPostCollideVoidFunctions_[f].second;
+ for( uint_t i = previousLevels; i != levels; ++i )
+ {
+ postCollideVoidFunctions_[i].push_back( std::make_pair( wrappedFunction, identifier ) );
+ if( timing_ && ! levelwiseTimingPool_->timerExists( getLevelwiseTimingPoolString( identifier, i ) ) )
+ levelwiseTimingPool_->registerTimer( getLevelwiseTimingPoolString( identifier, i ) );
+ }
+ }
+
+ for( uint_t f = uint_t(0); f != globalPostCollideBlockFunctions_.size(); ++f )
+ {
+ BlockFunctionWrappper wrappedFunction( globalPostCollideBlockFunctions_, f );
+ const std::string identifier = globalPostCollideBlockFunctions_[f].second;
+ for( uint_t i = previousLevels; i != levels; ++i )
+ {
+ postCollideBlockFunctions_[i].push_back( std::make_pair( wrappedFunction, identifier ) );
+ if( timing_ && ! levelwiseTimingPool_->timerExists( getLevelwiseTimingPoolString( identifier, i ) ) )
+ levelwiseTimingPool_->registerTimer( getLevelwiseTimingPoolString( identifier, i ) );
+ }
+ }
+
+ for( uint_t f = uint_t(0); f != globalPostBoundaryHandlingVoidFunctions_.size(); ++f )
+ {
+ VoidFunctionWrappper wrappedFunction( globalPostBoundaryHandlingVoidFunctions_, f );
+ const std::string identifier = globalPostBoundaryHandlingVoidFunctions_[f].second;
+ for( uint_t i = previousLevels; i != levels; ++i )
+ {
+ postBoundaryHandlingVoidFunctions_[i].push_back( std::make_pair( wrappedFunction, identifier ) );
+ if( timing_ && ! levelwiseTimingPool_->timerExists( getLevelwiseTimingPoolString( identifier, i ) ) )
+ levelwiseTimingPool_->registerTimer( getLevelwiseTimingPoolString( identifier, i ) );
+ }
+ }
+
+ for( uint_t f = uint_t(0); f != globalPostBoundaryHandlingBlockFunctions_.size(); ++f )
+ {
+ BlockFunctionWrappper wrappedFunction( globalPostBoundaryHandlingBlockFunctions_, f );
+ const std::string identifier = globalPostBoundaryHandlingBlockFunctions_[f].second;
+ for( uint_t i = previousLevels; i != levels; ++i )
+ {
+ postBoundaryHandlingBlockFunctions_[i].push_back( std::make_pair( wrappedFunction, identifier ) );
+ if( timing_ && ! levelwiseTimingPool_->timerExists( getLevelwiseTimingPoolString( identifier, i ) ) )
+ levelwiseTimingPool_->registerTimer( getLevelwiseTimingPoolString( identifier, i ) );
+ }
+ }
+
+ for( uint_t f = uint_t(0); f != globalPostStreamVoidFunctions_.size(); ++f )
+ {
+ VoidFunctionWrappper wrappedFunction( globalPostStreamVoidFunctions_, f );
+ const std::string identifier = globalPostStreamVoidFunctions_[f].second;
+ for( uint_t i = previousLevels; i != levels; ++i )
+ {
+ postStreamVoidFunctions_[i].push_back( std::make_pair( wrappedFunction, identifier ) );
+ if( timing_ && ! levelwiseTimingPool_->timerExists( getLevelwiseTimingPoolString( identifier, i ) ) )
+ levelwiseTimingPool_->registerTimer( getLevelwiseTimingPoolString( identifier, i ) );
+ }
+ }
+
+ for( uint_t f = uint_t(0); f != globalPostStreamBlockFunctions_.size(); ++f )
+ {
+ BlockFunctionWrappper wrappedFunction( globalPostStreamBlockFunctions_, f );
+ const std::string identifier = globalPostStreamBlockFunctions_[f].second;
+ for( uint_t i = previousLevels; i != levels; ++i )
+ {
+ postStreamBlockFunctions_[i].push_back( std::make_pair( wrappedFunction, identifier ) );
+ if( timing_ && ! levelwiseTimingPool_->timerExists( getLevelwiseTimingPoolString( identifier, i ) ) )
+ levelwiseTimingPool_->registerTimer( getLevelwiseTimingPoolString( identifier, i ) );
+ }
+ }
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+template< typename Function >
+void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addFunction( std::vector< std::vector< std::pair< Function, std::string > > > & functions,
+ const Function & function, const std::string & identifier )
+{
+ for( uint_t i = uint_t(0); i < functions.size(); ++i )
+ {
+ functions[i].push_back( std::make_pair( function, identifier ) );
+ if( timing_ && ! levelwiseTimingPool_->timerExists( getLevelwiseTimingPoolString( identifier, i ) ) )
+ levelwiseTimingPool_->registerTimer( getLevelwiseTimingPoolString( identifier, i ) );
+ }
+ if( timing_ && ! timingPool_->timerExists( getTimingPoolString( identifier ) ) )
+ timingPool_->registerTimer( getTimingPoolString( identifier ) );
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+template< typename Function >
+void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addFunction( std::vector< std::vector< std::pair< Function, std::string > > > & functions,
+ const Function & function, const std::string & identifier,
+ const uint_t level )
+{
+ if( level > ( postCollideVoidFunctions_.size() - uint_t(1) ) )
+ refresh( level + uint_t(1) );
+ functions[level].push_back( std::make_pair( function, identifier ) );
+ if( timing_ && ! levelwiseTimingPool_->timerExists( getLevelwiseTimingPoolString( identifier, level ) ) )
+ levelwiseTimingPool_->registerTimer( getLevelwiseTimingPoolString( identifier, level ) );
+ if( timing_ && ! timingPool_->timerExists( getTimingPoolString( identifier ) ) )
+ timingPool_->registerTimer( getTimingPoolString( identifier ) );
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostCollideVoidFunction( const VoidFunction & function, const std::string & identifier )
+{
+ VoidFunctionWrappper wrappedFunction( globalPostCollideVoidFunctions_, globalPostCollideVoidFunctions_.size() );
+ globalPostCollideVoidFunctions_.push_back( std::make_pair( function, identifier ) );
+ addFunction< VoidFunction >( postCollideVoidFunctions_, wrappedFunction, identifier );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostCollideBlockFunction( const BlockFunction & function, const std::string & identifier )
+{
+ BlockFunctionWrappper wrappedFunction( globalPostCollideBlockFunctions_, globalPostCollideBlockFunctions_.size() );
+ globalPostCollideBlockFunctions_.push_back( std::make_pair( function, identifier ) );
+ addFunction< BlockFunction >( postCollideBlockFunctions_, wrappedFunction, identifier );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostCollideVoidFunction( const VoidFunction & function, const std::string & identifier, const uint_t level )
+{
+ addFunction< VoidFunction >( postCollideVoidFunctions_, function, identifier, level );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostCollideBlockFunction( const BlockFunction & function, const std::string & identifier, const uint_t level )
+{
+ addFunction< BlockFunction >( postCollideBlockFunctions_, function, identifier, level );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostBoundaryHandlingVoidFunction( const VoidFunction & function, const std::string & identifier )
+{
+ VoidFunctionWrappper wrappedFunction( globalPostBoundaryHandlingVoidFunctions_, globalPostBoundaryHandlingVoidFunctions_.size() );
+ globalPostBoundaryHandlingVoidFunctions_.push_back( std::make_pair( function, identifier ) );
+ addFunction< VoidFunction >( postBoundaryHandlingVoidFunctions_, wrappedFunction, identifier );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostBoundaryHandlingBlockFunction( const BlockFunction & function, const std::string & identifier )
+{
+ BlockFunctionWrappper wrappedFunction( globalPostBoundaryHandlingBlockFunctions_, globalPostBoundaryHandlingBlockFunctions_.size() );
+ globalPostBoundaryHandlingBlockFunctions_.push_back( std::make_pair( function, identifier ) );
+ addFunction< BlockFunction >( postBoundaryHandlingBlockFunctions_, wrappedFunction, identifier );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostBoundaryHandlingVoidFunction( const VoidFunction & function, const std::string & identifier, const uint_t level )
+{
+ addFunction< VoidFunction >( postBoundaryHandlingVoidFunctions_, function, identifier, level );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostBoundaryHandlingBlockFunction( const BlockFunction & function, const std::string & identifier, const uint_t level )
+{
+ addFunction< BlockFunction >( postBoundaryHandlingBlockFunctions_, function, identifier, level );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostStreamVoidFunction( const VoidFunction & function, const std::string & identifier )
+{
+ VoidFunctionWrappper wrappedFunction( globalPostStreamVoidFunctions_, globalPostStreamVoidFunctions_.size() );
+ globalPostStreamVoidFunctions_.push_back( std::make_pair( function, identifier ) );
+ addFunction< VoidFunction >( postStreamVoidFunctions_, wrappedFunction, identifier );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostStreamBlockFunction( const BlockFunction & function, const std::string & identifier )
+{
+ BlockFunctionWrappper wrappedFunction( globalPostStreamBlockFunctions_, globalPostStreamBlockFunctions_.size() );
+ globalPostStreamBlockFunctions_.push_back( std::make_pair( function, identifier ) );
+ addFunction< BlockFunction >( postStreamBlockFunctions_, wrappedFunction, identifier );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostStreamVoidFunction( const VoidFunction & function, const std::string & identifier, const uint_t level )
+{
+ addFunction< VoidFunction >( postStreamVoidFunctions_, function, identifier, level );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostStreamBlockFunction( const BlockFunction & function, const std::string & identifier, const uint_t level )
+{
+ addFunction< BlockFunction >( postStreamBlockFunctions_, function, identifier, level );
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+template< typename F >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostCollideVoidFunction( const shared_ptr<F> & functorPtr, const std::string & identifier )
+{
+ addPostCollideVoidFunction( SharedVoidFunctor<F>(functorPtr), identifier );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+template< typename F >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostCollideBlockFunction( const shared_ptr<F> & functorPtr, const std::string & identifier )
+{
+ addPostCollideBlockFunction( SharedBlockFunctor<F>(functorPtr), identifier );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+template< typename F >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostCollideVoidFunction( const shared_ptr<F> & functorPtr, const std::string & identifier, const uint_t level )
+{
+ addPostCollideVoidFunction( SharedVoidFunctor<F>(functorPtr), identifier, level );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+template< typename F >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostCollideBlockFunction( const shared_ptr<F> & functorPtr, const std::string & identifier, const uint_t level )
+{
+ addPostCollideBlockFunction( SharedBlockFunctor<F>(functorPtr), identifier, level );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+template< typename F >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostBoundaryHandlingVoidFunction( const shared_ptr<F> & functorPtr, const std::string & identifier )
+{
+ addPostBoundaryHandlingVoidFunction( SharedVoidFunctor<F>(functorPtr), identifier );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+template< typename F >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostBoundaryHandlingBlockFunction( const shared_ptr<F> & functorPtr, const std::string & identifier )
+{
+ addPostBoundaryHandlingBlockFunction( SharedBlockFunctor<F>(functorPtr), identifier );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+template< typename F >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostBoundaryHandlingVoidFunction( const shared_ptr<F> & functorPtr, const std::string & identifier, const uint_t level )
+{
+ addPostBoundaryHandlingVoidFunction( SharedVoidFunctor<F>(functorPtr), identifier, level );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+template< typename F >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostBoundaryHandlingBlockFunction( const shared_ptr<F> & functorPtr, const std::string & identifier, const uint_t level )
+{
+ addPostBoundaryHandlingBlockFunction( SharedBlockFunctor<F>(functorPtr), identifier, level );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+template< typename F >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostStreamVoidFunction( const shared_ptr<F> & functorPtr, const std::string & identifier )
+{
+ addPostStreamVoidFunction( SharedVoidFunctor<F>(functorPtr), identifier );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+template< typename F >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostStreamBlockFunction( const shared_ptr<F> & functorPtr, const std::string & identifier )
+{
+ addPostStreamBlockFunction( SharedBlockFunctor<F>(functorPtr), identifier );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+template< typename F >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostStreamVoidFunction( const shared_ptr<F> & functorPtr, const std::string & identifier, const uint_t level )
+{
+ addPostStreamVoidFunction( SharedVoidFunctor<F>(functorPtr), identifier, level );
+}
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+template< typename F >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::addPostStreamBlockFunction( const shared_ptr<F> & functorPtr, const std::string & identifier, const uint_t level )
+{
+ addPostStreamBlockFunction( SharedBlockFunctor<F>(functorPtr), identifier, level );
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+std::string TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::getTimingPoolString( const std::string & name,
+ const std::string & suffix ) const
+{
+ std::ostringstream oss;
+ oss << name;
+ if( !suffix.empty() )
+ oss << " " << suffix;
+ return oss.str();
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+std::string TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::getLevelwiseTimingPoolString( const std::string & name, const uint_t level,
+ const std::string & suffix ) const
+{
+ std::ostringstream oss;
+ oss << name << " (" << level << ")";
+ if( !suffix.empty() )
+ oss << " " << suffix;
+ return oss.str();
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::createTimers( const uint_t levels )
+{
+ // unify registered timers across all processes
+
+ if( timing_ )
+ {
+ std::vector< std::string > timers;
+ timers.push_back( getTimingPoolString( "boundary handling" ) );
+ timers.push_back( getTimingPoolString( "collide" ) );
+ timers.push_back( getTimingPoolString( "communication coarse to fine", "[pack & send]" ) );
+ timers.push_back( getTimingPoolString( "communication coarse to fine", "[wait & unpack]" ) );
+ timers.push_back( getTimingPoolString( "communication equal level", "[pack & send]" ) );
+ timers.push_back( getTimingPoolString( "communication equal level", "[wait & unpack]" ) );
+ timers.push_back( getTimingPoolString( "communication fine to coarse", "[pack & send]" ) );
+ timers.push_back( getTimingPoolString( "communication fine to coarse", "[wait & unpack]" ) );
+ timers.push_back( getTimingPoolString( "equal level border stream correction", "[prepare]" ) );
+ timers.push_back( getTimingPoolString( "equal level border stream correction", "[apply]" ) );
+ timers.push_back( getTimingPoolString( "linear explosion" ) );
+ timers.push_back( getTimingPoolString( "stream" ) );
+ timers.push_back( getTimingPoolString( "stream & collide" ) );
+
+ for( auto timer = timers.begin(); timer != timers.end(); ++timer )
+ if( ! timingPool_->timerExists( *timer ) )
+ timingPool_->registerTimer( *timer );
+
+ timers.clear();
+ timers.push_back( getLevelwiseTimingPoolString( "stream & collide", levels - uint_t(1) ) );
+ for( uint_t i = uint_t(0); i < levels; ++i )
+ {
+ timers.push_back( getLevelwiseTimingPoolString( "boundary handling", i ) );
+ timers.push_back( getLevelwiseTimingPoolString( "collide", i ) );
+ timers.push_back( getLevelwiseTimingPoolString( "communication equal level", i, "[pack & send]" ) );
+ timers.push_back( getLevelwiseTimingPoolString( "communication equal level", i, "[wait & unpack]" ) );
+ timers.push_back( getLevelwiseTimingPoolString( "linear explosion", i ) );
+ timers.push_back( getLevelwiseTimingPoolString( "stream", i ) );
+ if( i != uint_t(0) )
+ {
+ timers.push_back( getLevelwiseTimingPoolString( "communication coarse to fine", i, "[pack & send]" ) );
+ timers.push_back( getLevelwiseTimingPoolString( "communication coarse to fine", i, "[wait & unpack]" ) );
+ timers.push_back( getLevelwiseTimingPoolString( "communication fine to coarse", i, "[pack & send]" ) );
+ timers.push_back( getLevelwiseTimingPoolString( "communication fine to coarse", i, "[wait & unpack]" ) );
+ }
+ if( i != ( levels - uint_t(1) ) )
+ {
+ timers.push_back( getLevelwiseTimingPoolString( "equal level border stream correction", i, "[prepare]" ) );
+ timers.push_back( getLevelwiseTimingPoolString( "equal level border stream correction", i, "[apply]" ) );
+ }
+ }
+
+ for( auto timer = timers.begin(); timer != timers.end(); ++timer )
+ if( ! levelwiseTimingPool_->timerExists( *timer ) )
+ levelwiseTimingPool_->registerTimer( *timer );
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+std::vector< Block * > TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::selectedBlocks( const uint_t level ) const
+{
+ std::vector< Block * > levelBlocks;
+ auto blocks = blocks_.lock();
+ WALBERLA_CHECK_NOT_NULLPTR( blocks, "Trying to access 'TimeStep' (refinement) for a block storage object that doesn't exist anymore" );
+ blocks->getBlocks( levelBlocks, level );
+
+ std::vector< Block * > sBlocks;
+
+ for( auto block = levelBlocks.begin(); block != levelBlocks.end(); ++block )
+ if( selectable::isSetSelected( (*block)->getState(), requiredBlockSelectors_, incompatibleBlockSelectors_ ) )
+ sBlocks.push_back( *block );
+
+ return sBlocks;
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::startTiming( const std::string & name, const uint_t level,
+ const std::string & suffix )
+{
+ (*timingPool_)[ getTimingPoolString( name, suffix ) ].start();
+ (*levelwiseTimingPool_)[ getLevelwiseTimingPoolString( name, level, suffix ) ].start();
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+inline void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::stopTiming( const std::string & name, const uint_t level,
+ const std::string & suffix )
+{
+ (*timingPool_)[ getTimingPoolString( name, suffix ) ].end();
+ (*levelwiseTimingPool_)[ getLevelwiseTimingPoolString( name, level, suffix ) ].end();
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::collide( std::vector< Block * > & blocks, const uint_t level, const uint_t executionCount )
+{
+ if( timing_ )
+ {
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ startTiming( "collide", level );
+ sweep_->collide( *block );
+ stopTiming( "collide", level );
+
+ for( auto func = postCollideBlockFunctions_[level].begin(); func != postCollideBlockFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( *block, level, executionCount );
+ stopTiming( func->second, level );
+ }
+ }
+ for( auto func = postCollideVoidFunctions_[level].begin(); func != postCollideVoidFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( level, executionCount );
+ stopTiming( func->second, level );
+ }
+ }
+ else
+ {
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ sweep_->collide( *block );
+
+ for( auto func = postCollideBlockFunctions_[level].begin(); func != postCollideBlockFunctions_[level].end(); ++func )
+ (func->first)( *block, level, executionCount );
+ }
+ for( auto func = postCollideVoidFunctions_[level].begin(); func != postCollideVoidFunctions_[level].end(); ++func )
+ (func->first)( level, executionCount );
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::stream( std::vector< Block * > & blocks, const uint_t level, const uint_t executionCount )
+{
+ auto _blocks = blocks_.lock();
+ WALBERLA_CHECK_NOT_NULLPTR( _blocks, "Trying to access 'TimeStep' (refinement) for a block storage object that doesn't exist anymore" );
+
+ const uint_t finestLevel = _blocks->getNumberOfLevels() - uint_t(1);
+ const bool doEqualLevelBorderStreamCorrection = ( performEqualLevelBorderStreamCorrection_ && level != finestLevel );
+
+ if( timing_ )
+ {
+ if( postBoundaryHandlingVoidFunctions_[level].empty() )
+ {
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ bool coarseNeighbors = false;
+ for( uint_t i = uint_t(0); i < uint_t(26) && !coarseNeighbors; ++i )
+ coarseNeighbors = (*block)->neighborhoodSectionHasLargerBlock(i);
+
+ if( coarseNeighbors )
+ {
+ startTiming( "boundary handling", level );
+ boundarySweepWithLayers_( *block );
+ stopTiming( "boundary handling", level );
+
+ for( auto func = postBoundaryHandlingBlockFunctions_[level].begin(); func != postBoundaryHandlingBlockFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( *block, level, executionCount );
+ stopTiming( func->second, level );
+ }
+
+ startTiming( "stream", level );
+ sweep_->stream( *block, StreamIncludedGhostLayers );
+ stopTiming( "stream", level );
+ }
+ else
+ {
+ startTiming( "boundary handling", level );
+ boundarySweep_( *block );
+ stopTiming( "boundary handling", level );
+
+ for( auto func = postBoundaryHandlingBlockFunctions_[level].begin(); func != postBoundaryHandlingBlockFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( *block, level, executionCount );
+ stopTiming( func->second, level );
+ }
+
+ startTiming( "stream", level );
+ sweep_->stream( *block, uint_t(0) );
+ stopTiming( "stream", level );
+ }
+
+ if( doEqualLevelBorderStreamCorrection )
+ {
+ startTiming( "equal level border stream correction", level, "[prepare]" );
+ equalLevelBorderStreamCorrection_.prepare( *block );
+ stopTiming( "equal level border stream correction", level, "[prepare]" );
+ }
+ }
+ }
+ else
+ {
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ bool coarseNeighbors = false;
+ for( uint_t i = uint_t(0); i < uint_t(26) && !coarseNeighbors; ++i )
+ coarseNeighbors = (*block)->neighborhoodSectionHasLargerBlock(i);
+
+ if( coarseNeighbors )
+ {
+ startTiming( "boundary handling", level );
+ boundarySweepWithLayers_( *block );
+ stopTiming( "boundary handling", level );
+ }
+ else
+ {
+ startTiming( "boundary handling", level );
+ boundarySweep_( *block );
+ stopTiming( "boundary handling", level );
+ }
+ for( auto func = postBoundaryHandlingBlockFunctions_[level].begin(); func != postBoundaryHandlingBlockFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( *block, level, executionCount );
+ stopTiming( func->second, level );
+ }
+ }
+ for( auto func = postBoundaryHandlingVoidFunctions_[level].begin(); func != postBoundaryHandlingVoidFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( level, executionCount );
+ stopTiming( func->second, level );
+ }
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ bool coarseNeighbors = false;
+ for( uint_t i = uint_t(0); i < uint_t(26) && !coarseNeighbors; ++i )
+ coarseNeighbors = (*block)->neighborhoodSectionHasLargerBlock(i);
+
+ if( coarseNeighbors )
+ {
+ startTiming( "stream", level );
+ sweep_->stream( *block, StreamIncludedGhostLayers );
+ stopTiming( "stream", level );
+ }
+ else
+ {
+ startTiming( "stream", level );
+ sweep_->stream( *block, uint_t(0) );
+ stopTiming( "stream", level );
+ }
+
+ if( doEqualLevelBorderStreamCorrection )
+ {
+ startTiming( "equal level border stream correction", level, "[prepare]" );
+ equalLevelBorderStreamCorrection_.prepare( *block );
+ stopTiming( "equal level border stream correction", level, "[prepare]" );
+ }
+ }
+ }
+ }
+ else
+ {
+ if( postBoundaryHandlingVoidFunctions_[level].empty() )
+ {
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ bool coarseNeighbors = false;
+ for( uint_t i = uint_t(0); i < uint_t(26) && !coarseNeighbors; ++i )
+ coarseNeighbors = (*block)->neighborhoodSectionHasLargerBlock(i);
+
+ if( coarseNeighbors )
+ {
+ boundarySweepWithLayers_( *block );
+ for( auto func = postBoundaryHandlingBlockFunctions_[level].begin(); func != postBoundaryHandlingBlockFunctions_[level].end(); ++func )
+ (func->first)( *block, level, executionCount );
+ sweep_->stream( *block, StreamIncludedGhostLayers );
+ }
+ else
+ {
+ boundarySweep_( *block );
+ for( auto func = postBoundaryHandlingBlockFunctions_[level].begin(); func != postBoundaryHandlingBlockFunctions_[level].end(); ++func )
+ (func->first)( *block, level, executionCount );
+ sweep_->stream( *block, uint_t(0) );
+ }
+
+ if( doEqualLevelBorderStreamCorrection )
+ equalLevelBorderStreamCorrection_.prepare( *block );
+ }
+ }
+ else
+ {
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ bool coarseNeighbors = false;
+ for( uint_t i = uint_t(0); i < uint_t(26) && !coarseNeighbors; ++i )
+ coarseNeighbors = (*block)->neighborhoodSectionHasLargerBlock(i);
+
+ if( coarseNeighbors )
+ {
+ boundarySweepWithLayers_( *block );
+ }
+ else
+ {
+ boundarySweep_( *block );
+ }
+ for( auto func = postBoundaryHandlingBlockFunctions_[level].begin(); func != postBoundaryHandlingBlockFunctions_[level].end(); ++func )
+ (func->first)( *block, level, executionCount );
+ }
+ for( auto func = postBoundaryHandlingVoidFunctions_[level].begin(); func != postBoundaryHandlingVoidFunctions_[level].end(); ++func )
+ (func->first)( level, executionCount );
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ bool coarseNeighbors = false;
+ for( uint_t i = uint_t(0); i < uint_t(26) && !coarseNeighbors; ++i )
+ coarseNeighbors = (*block)->neighborhoodSectionHasLargerBlock(i);
+
+ if( coarseNeighbors )
+ {
+ sweep_->stream( *block, StreamIncludedGhostLayers );
+ }
+ else
+ {
+ sweep_->stream( *block, uint_t(0) );
+ }
+
+ if( doEqualLevelBorderStreamCorrection )
+ equalLevelBorderStreamCorrection_.prepare( *block );
+ }
+ }
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::finishStream( std::vector< Block * > & blocks, const uint_t level, const uint_t executionCount )
+{
+ auto _blocks = blocks_.lock();
+ WALBERLA_CHECK_NOT_NULLPTR( _blocks, "Trying to access 'TimeStep' (refinement) for a block storage object that doesn't exist anymore" );
+
+ const uint_t finestLevel = _blocks->getNumberOfLevels() - uint_t(1);
+ const bool doEqualLevelBorderStreamCorrection = ( performEqualLevelBorderStreamCorrection_ && level != finestLevel );
+
+ if( timing_ )
+ {
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ if( doEqualLevelBorderStreamCorrection )
+ {
+ startTiming( "equal level border stream correction", level, "[apply]" );
+ equalLevelBorderStreamCorrection_.apply( *block );
+ stopTiming( "equal level border stream correction", level, "[apply]" );
+ }
+ for( auto func = postStreamBlockFunctions_[level].begin(); func != postStreamBlockFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( *block, level, executionCount );
+ stopTiming( func->second, level );
+ }
+ }
+ for( auto func = postStreamVoidFunctions_[level].begin(); func != postStreamVoidFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( level, executionCount );
+ stopTiming( func->second, level );
+ }
+ }
+ else
+ {
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ if( doEqualLevelBorderStreamCorrection )
+ equalLevelBorderStreamCorrection_.apply( *block );
+
+ for( auto func = postStreamBlockFunctions_[level].begin(); func != postStreamBlockFunctions_[level].end(); ++func )
+ (func->first)( *block, level, executionCount );
+ }
+ for( auto func = postStreamVoidFunctions_[level].begin(); func != postStreamVoidFunctions_[level].end(); ++func )
+ (func->first)( level, executionCount );
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::streamCollide( std::vector< Block * > & blocks, const uint_t level, const uint_t executionCount )
+{
+#ifndef NDEBUG
+ auto _blocks = blocks_.lock();
+ WALBERLA_CHECK_NOT_NULLPTR( _blocks, "Trying to access 'TimeStep' (refinement) for a block storage object that doesn't exist anymore" );
+ WALBERLA_ASSERT_EQUAL( level, _blocks->getNumberOfLevels() - uint_t(1) );
+#endif
+
+ if( postStreamVoidFunctions_[level].empty() )
+ {
+ if( timing_ )
+ {
+ if( postBoundaryHandlingVoidFunctions_[level].empty() )
+ {
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ bool coarseNeighborsOrPostStreamFunctions = !(postStreamBlockFunctions_[level].empty());
+ for( uint_t i = uint_t(0); i < uint_t(26) && !coarseNeighborsOrPostStreamFunctions; ++i )
+ coarseNeighborsOrPostStreamFunctions = (*block)->neighborhoodSectionHasLargerBlock(i);
+
+ if( coarseNeighborsOrPostStreamFunctions )
+ {
+ startTiming( "boundary handling", level );
+ boundarySweepWithLayers_( *block );
+ stopTiming( "boundary handling", level );
+
+ for( auto func = postBoundaryHandlingBlockFunctions_[level].begin(); func != postBoundaryHandlingBlockFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( *block, level, executionCount );
+ stopTiming( func->second, level );
+ }
+
+ startTiming( "stream", level );
+ sweep_->stream( *block, StreamIncludedGhostLayers );
+ stopTiming( "stream", level );
+
+ for( auto func = postStreamBlockFunctions_[level].begin(); func != postStreamBlockFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( *block, level, executionCount );
+ stopTiming( func->second, level );
+ }
+
+ startTiming( "collide", level );
+ sweep_->collide( *block );
+ stopTiming( "collide", level );
+ }
+ else
+ {
+ startTiming( "boundary handling", level );
+ boundarySweep_( *block );
+ stopTiming( "boundary handling", level );
+
+ for( auto func = postBoundaryHandlingBlockFunctions_[level].begin(); func != postBoundaryHandlingBlockFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( *block, level, executionCount );
+ stopTiming( func->second, level );
+ }
+
+ startTiming( "stream & collide", level );
+ (*sweep_)( *block );
+ stopTiming( "stream & collide", level );
+ }
+ for( auto func = postCollideBlockFunctions_[level].begin(); func != postCollideBlockFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( *block, level, executionCount + uint_t(1) );
+ stopTiming( func->second, level );
+ }
+ }
+ for( auto func = postCollideVoidFunctions_[level].begin(); func != postCollideVoidFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( level, executionCount + uint_t(1) );
+ stopTiming( func->second, level );
+ }
+ }
+ else
+ {
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ bool coarseNeighbors = false;
+ for( uint_t i = uint_t(0); i < uint_t(26) && !coarseNeighbors; ++i )
+ coarseNeighbors = (*block)->neighborhoodSectionHasLargerBlock(i);
+
+ if( coarseNeighbors )
+ {
+ startTiming( "boundary handling", level );
+ boundarySweepWithLayers_( *block );
+ stopTiming( "boundary handling", level );
+ }
+ else
+ {
+ startTiming( "boundary handling", level );
+ boundarySweep_( *block );
+ stopTiming( "boundary handling", level );
+ }
+ for( auto func = postBoundaryHandlingBlockFunctions_[level].begin(); func != postBoundaryHandlingBlockFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( *block, level, executionCount );
+ stopTiming( func->second, level );
+ }
+ }
+ for( auto func = postBoundaryHandlingVoidFunctions_[level].begin(); func != postBoundaryHandlingVoidFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( level, executionCount );
+ stopTiming( func->second, level );
+ }
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ bool coarseNeighborsOrPostStreamFunctions = !(postStreamBlockFunctions_[level].empty());
+ for( uint_t i = uint_t(0); i < uint_t(26) && !coarseNeighborsOrPostStreamFunctions; ++i )
+ coarseNeighborsOrPostStreamFunctions = (*block)->neighborhoodSectionHasLargerBlock(i);
+
+ if( coarseNeighborsOrPostStreamFunctions )
+ {
+ startTiming( "stream", level );
+ sweep_->stream( *block, StreamIncludedGhostLayers );
+ stopTiming( "stream", level );
+
+ for( auto func = postStreamBlockFunctions_[level].begin(); func != postStreamBlockFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( *block, level, executionCount );
+ stopTiming( func->second, level );
+ }
+
+ startTiming( "collide", level );
+ sweep_->collide( *block );
+ stopTiming( "collide", level );
+ }
+ else
+ {
+ startTiming( "stream & collide", level );
+ (*sweep_)( *block );
+ stopTiming( "stream & collide", level );
+ }
+ for( auto func = postCollideBlockFunctions_[level].begin(); func != postCollideBlockFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( *block, level, executionCount + uint_t(1) );
+ stopTiming( func->second, level );
+ }
+ }
+ for( auto func = postCollideVoidFunctions_[level].begin(); func != postCollideVoidFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( level, executionCount + uint_t(1) );
+ stopTiming( func->second, level );
+ }
+ }
+ }
+ else
+ {
+ if( postBoundaryHandlingVoidFunctions_[level].empty() )
+ {
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ bool coarseNeighborsOrPostStreamFunctions = !(postStreamBlockFunctions_[level].empty());
+ for( uint_t i = uint_t(0); i < uint_t(26) && !coarseNeighborsOrPostStreamFunctions; ++i )
+ coarseNeighborsOrPostStreamFunctions = (*block)->neighborhoodSectionHasLargerBlock(i);
+
+ if( coarseNeighborsOrPostStreamFunctions )
+ {
+ boundarySweepWithLayers_( *block );
+ for( auto func = postBoundaryHandlingBlockFunctions_[level].begin(); func != postBoundaryHandlingBlockFunctions_[level].end(); ++func )
+ (func->first)( *block, level, executionCount );
+ sweep_->stream( *block, StreamIncludedGhostLayers );
+ for( auto func = postStreamBlockFunctions_[level].begin(); func != postStreamBlockFunctions_[level].end(); ++func )
+ (func->first)( *block, level, executionCount );
+ sweep_->collide( *block );
+ }
+ else
+ {
+ boundarySweep_( *block );
+ for( auto func = postBoundaryHandlingBlockFunctions_[level].begin(); func != postBoundaryHandlingBlockFunctions_[level].end(); ++func )
+ (func->first)( *block, level, executionCount );
+ (*sweep_)( *block );
+ }
+ for( auto func = postCollideBlockFunctions_[level].begin(); func != postCollideBlockFunctions_[level].end(); ++func )
+ (func->first)( *block, level, executionCount + uint_t(1) );
+ }
+ for( auto func = postCollideVoidFunctions_[level].begin(); func != postCollideVoidFunctions_[level].end(); ++func )
+ (func->first)( level, executionCount + uint_t(1) );
+ }
+ else
+ {
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ bool coarseNeighbors = false;
+ for( uint_t i = uint_t(0); i < uint_t(26) && !coarseNeighbors; ++i )
+ coarseNeighbors = (*block)->neighborhoodSectionHasLargerBlock(i);
+
+ if( coarseNeighbors )
+ {
+ boundarySweepWithLayers_( *block );
+ }
+ else
+ {
+ boundarySweep_( *block );
+ }
+ for( auto func = postBoundaryHandlingBlockFunctions_[level].begin(); func != postBoundaryHandlingBlockFunctions_[level].end(); ++func )
+ (func->first)( *block, level, executionCount );
+ }
+ for( auto func = postBoundaryHandlingVoidFunctions_[level].begin(); func != postBoundaryHandlingVoidFunctions_[level].end(); ++func )
+ (func->first)( level, executionCount );
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ bool coarseNeighborsOrPostStreamFunctions = !(postStreamBlockFunctions_[level].empty());
+ for( uint_t i = uint_t(0); i < uint_t(26) && !coarseNeighborsOrPostStreamFunctions; ++i )
+ coarseNeighborsOrPostStreamFunctions = (*block)->neighborhoodSectionHasLargerBlock(i);
+
+ if( coarseNeighborsOrPostStreamFunctions )
+ {
+ sweep_->stream( *block, StreamIncludedGhostLayers );
+ for( auto func = postStreamBlockFunctions_[level].begin(); func != postStreamBlockFunctions_[level].end(); ++func )
+ (func->first)( *block, level, executionCount );
+ sweep_->collide( *block );
+ }
+ else
+ {
+ (*sweep_)( *block );
+ }
+ for( auto func = postCollideBlockFunctions_[level].begin(); func != postCollideBlockFunctions_[level].end(); ++func )
+ (func->first)( *block, level, executionCount + uint_t(1) );
+ }
+ for( auto func = postCollideVoidFunctions_[level].begin(); func != postCollideVoidFunctions_[level].end(); ++func )
+ (func->first)( level, executionCount + uint_t(1) );
+ }
+ }
+ }
+ else
+ {
+ stream( blocks, level, executionCount );
+
+ if( timing_ )
+ {
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ for( auto func = postStreamBlockFunctions_[level].begin(); func != postStreamBlockFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( *block, level, executionCount );
+ stopTiming( func->second, level );
+ }
+ }
+ for( auto func = postStreamVoidFunctions_[level].begin(); func != postStreamVoidFunctions_[level].end(); ++func )
+ {
+ startTiming( func->second, level );
+ (func->first)( level, executionCount );
+ stopTiming( func->second, level );
+ }
+ }
+ else
+ {
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ for( auto func = postStreamBlockFunctions_[level].begin(); func != postStreamBlockFunctions_[level].end(); ++func )
+ (func->first)( *block, level, executionCount );
+ }
+ for( auto func = postStreamVoidFunctions_[level].begin(); func != postStreamVoidFunctions_[level].end(); ++func )
+ (func->first)( level, executionCount );
+ }
+
+ collide( blocks, level, executionCount + uint_t(1) );
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::startCommunicationEqualLevel( const uint_t level )
+{
+ if( timing_ )
+ {
+ startTiming( "communication equal level", level, "[pack & send]" );
+ communication_.startCommunicateEqualLevel( level );
+ stopTiming( "communication equal level", level, "[pack & send]" );
+ }
+ else
+ {
+ communication_.startCommunicateEqualLevel( level );
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::endCommunicationEqualLevel( const uint_t level )
+{
+ if( timing_ )
+ {
+ startTiming( "communication equal level", level, "[wait & unpack]" );
+ communication_.waitCommunicateEqualLevel( level );
+ stopTiming( "communication equal level", level, "[wait & unpack]" );
+ }
+ else
+ {
+ communication_.waitCommunicateEqualLevel( level );
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::startCommunicationCoarseToFine( const uint_t level )
+{
+ if( timing_ )
+ {
+ startTiming( "communication coarse to fine", level, "[pack & send]" );
+ communication_.startCommunicateCoarseToFine( level );
+ stopTiming( "communication coarse to fine", level, "[pack & send]" );
+ }
+ else
+ {
+ communication_.startCommunicateCoarseToFine( level );
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::endCommunicationCoarseToFine( const uint_t level )
+{
+ if( timing_ )
+ {
+ startTiming( "communication coarse to fine", level, "[wait & unpack]" );
+ communication_.waitCommunicateCoarseToFine( level );
+ stopTiming( "communication coarse to fine", level, "[wait & unpack]" );
+ }
+ else
+ {
+ communication_.waitCommunicateCoarseToFine( level );
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::startCommunicationFineToCoarse( const uint_t level )
+{
+ if( timing_ )
+ {
+ startTiming( "communication fine to coarse", level, "[pack & send]" );
+ communication_.startCommunicateFineToCoarse( level );
+ stopTiming( "communication fine to coarse", level, "[pack & send]" );
+ }
+ else
+ {
+ communication_.startCommunicateFineToCoarse( level );
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::endCommunicationFineToCoarse( const uint_t level )
+{
+ if( timing_ )
+ {
+ startTiming( "communication fine to coarse", level, "[wait & unpack]" );
+ communication_.waitCommunicateFineToCoarse( level );
+ stopTiming( "communication fine to coarse", level, "[wait & unpack]" );
+ }
+ else
+ {
+ communication_.waitCommunicateFineToCoarse( level );
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::performLinearExplosion( std::vector< Block * > & blocks, const uint_t level )
+{
+ if( performLinearExplosion_ )
+ {
+ if( timing_ )
+ {
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ {
+ startTiming( "linear explosion", level );
+ linearExplosion_( *block );
+ stopTiming( "linear explosion", level );
+ }
+ }
+ else
+ {
+ for( auto block = blocks.begin(); block != blocks.end(); ++block )
+ linearExplosion_( *block );
+ }
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename Sweep_T, typename BoundaryHandling_T >
+void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::recursiveStep( const uint_t level, const uint_t executionCount )
+{
+ // -> "A generic, mass conservative local grid refinement technique for lattice-Boltzmann schemes", Rhode et al., 2005
+
+ auto _blocks = blocks_.lock();
+ WALBERLA_CHECK_NOT_NULLPTR( _blocks, "Trying to access 'TimeStep' (refinement) for a block storage object that doesn't exist anymore" );
+
+ const uint_t coarsestLevel = uint_t(0);
+ const uint_t finestLevel = _blocks->getNumberOfLevels() - uint_t(1);
+
+ const uint_t executionCount1st = (executionCount + uint_t(1)) * uint_t(2) - uint_t(2);
+ const uint_t executionCount2nd = (executionCount + uint_t(1)) * uint_t(2) - uint_t(1);
+
+ std::vector< Block * > blocks = selectedBlocks( level );
+
+ WALBERLA_LOG_DETAIL("Starting recursive step with level " << level << " and execution count " << executionCount);
+
+ if( asynchronousCommunication_ && level != coarsestLevel )
+ {
+ WALBERLA_LOG_DETAIL("Start communication coarse to fine, initiated by fine level " << level );
+ startCommunicationCoarseToFine( level ); // [start] explosion (initiated by fine level, involves "level" and "level-1")
+ }
+
+ WALBERLA_LOG_DETAIL("Colliding on level " << level);
+ collide( blocks, level, executionCount1st );
+
+ if( asynchronousCommunication_ )
+ {
+ WALBERLA_LOG_DETAIL("Start communication equal level, initiated by level " << level );
+ startCommunicationEqualLevel( level ); // [start] equal level communication
+ }
+
+ if( level != finestLevel )
+ {
+ WALBERLA_LOG_DETAIL("Calling recursive step with level " << level + uint_t(1) << " and execution count " << executionCount1st );
+ recursiveStep( level + uint_t(1), executionCount1st );
+
+ if( asynchronousCommunication_ ) {
+ WALBERLA_LOG_DETAIL("Start communication fine to coarse, initiated by coarse level " << level );
+ startCommunicationFineToCoarse(level + uint_t(1)); // [start] coalescence (initiated by coarse level)
+ }
+ }
+
+ if( level != coarsestLevel )
+ {
+ if( !asynchronousCommunication_ ) {
+ WALBERLA_LOG_DETAIL("Start communication coarse to fine, initiated by fine level " << level );
+ startCommunicationCoarseToFine(level); // [start] explosion (initiated by fine level, involves "level" and "level-1")
+ }
+ WALBERLA_LOG_DETAIL("End communication coarse to fine, initiated by fine level " << level );
+ endCommunicationCoarseToFine( level ); // [end] explosion (initiated by fine level, involves "level" and "level-1")
+ WALBERLA_LOG_DETAIL("Perform linear explosion on level " << level );
+ performLinearExplosion( blocks, level );
+ }
+
+ if( !asynchronousCommunication_ ) {
+ WALBERLA_LOG_DETAIL("Start communication equal level, initiated by level " << level );
+ startCommunicationEqualLevel(level); // [start] equal level communication
+ }
+ WALBERLA_LOG_DETAIL("End communication equal level, initiated by level " << level );
+ endCommunicationEqualLevel( level ); // [end] equal level communication
+
+ // performLinearExplosion( blocks, level ); // if equal level neighbor values are needed, linear explosion should be performed here
+
+ if( level == finestLevel && level != coarsestLevel )
+ {
+ WALBERLA_LOG_DETAIL("Stream + collide on level " << level );
+ streamCollide( blocks, level, executionCount1st );
+ }
+ else
+ {
+ WALBERLA_LOG_DETAIL("Stream on level " << level );
+ stream( blocks, level, executionCount1st );
+
+ if( level != finestLevel )
+ {
+ if( !asynchronousCommunication_ ) {
+ WALBERLA_LOG_DETAIL("Start communication fine to coarse, initiated by coarse level " << level );
+ startCommunicationFineToCoarse(level + uint_t(1)); // [start] coalescence (initiated by coarse level)
+ }
+ WALBERLA_LOG_DETAIL("End communication fine to coarse, initiated by coarse level " << level );
+ endCommunicationFineToCoarse( level + uint_t(1) ); // [end] coalescence (initiated by coarse level)
+ }
+
+ WALBERLA_LOG_DETAIL("Finish stream on level " << level );
+ finishStream( blocks, level, executionCount1st );
+
+ if( level == coarsestLevel )
+ {
+ WALBERLA_LOG_DETAIL("End recursive step on level " << level );
+ return;
+ }
+
+ WALBERLA_LOG_DETAIL("Colliding on level " << level);
+ collide( blocks, level, executionCount2nd );
+ }
+
+ if( asynchronousCommunication_ ) {
+ WALBERLA_LOG_DETAIL("Start communication equal level, initiated by level " << level );
+ startCommunicationEqualLevel(level); // [start] equal level communication
+ }
+
+ if( level != finestLevel )
+ {
+ WALBERLA_LOG_DETAIL("Calling recursive step with level " << level + uint_t(1) << " and execution count " << executionCount2nd );
+ recursiveStep( level + uint_t(1), executionCount2nd );
+
+ if( asynchronousCommunication_ ) {
+ WALBERLA_LOG_DETAIL("Start communication fine to coarse, initiated by coarse level " << level );
+ startCommunicationFineToCoarse(level + uint_t(1)); // [start] coalescence (initiated by coarse level)
+ }
+ }
+
+ if( !asynchronousCommunication_ ) {
+ WALBERLA_LOG_DETAIL("Start communication equal level, initiated by level " << level );
+ startCommunicationEqualLevel(level); // [start] equal level communication
+ }
+ WALBERLA_LOG_DETAIL("End communication equal level, initiated by level " << level );
+ endCommunicationEqualLevel( level ); // [end] equal level communication
+
+ WALBERLA_LOG_DETAIL("Stream on level " << level );
+ stream( blocks, level, executionCount2nd );
+
+ if( level != finestLevel )
+ {
+ if( !asynchronousCommunication_ )
+ {
+ WALBERLA_LOG_DETAIL("Start communication fine to coarse, initiated by coarse level " << level );
+ startCommunicationFineToCoarse( level + uint_t(1) ); // [start] coalescence (initiated by coarse level)
+ }
+ WALBERLA_LOG_DETAIL("End communication fine to coarse, initiated by coarse level " << level );
+ endCommunicationFineToCoarse( level + uint_t(1) ); // [end] coalescence (initiated by coarse level)
+ }
+
+ WALBERLA_LOG_DETAIL("Finish stream on level " << level );
+ finishStream( blocks, level, executionCount2nd );
+
+ WALBERLA_LOG_DETAIL("End recursive step on level " << level );
+}
+
+
+
+
+template< typename LatticeModel_T, typename BoundaryHandling_T, typename Sweep_T >
+shared_ptr< TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T > >
+makeTimeStep( weak_ptr< StructuredBlockForest > blocks, shared_ptr< Sweep_T > & sweep,
+ const BlockDataID & pdfFieldId, const BlockDataID & boundaryHandlingId,
+ const Set<SUID> & requiredBlockSelectors = Set<SUID>::emptySet(),
+ const Set<SUID> & incompatibleBlockSelectors = Set<SUID>::emptySet() )
+{
+ typedef TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T > TS_T;
+ return shared_ptr< TS_T >( new TS_T( blocks, sweep, pdfFieldId, boundaryHandlingId, requiredBlockSelectors, incompatibleBlockSelectors ) );
+}
+
+
+
+template< typename LatticeModel_T, typename BoundaryHandling_T, typename Sweep_T >
+shared_ptr< TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T > >
+makeTimeStep( weak_ptr< StructuredBlockForest > blocks, shared_ptr< Sweep_T > & sweep,
+ const BlockDataID & pdfFieldId, const BlockDataID & boundaryHandlingId,
+ const shared_ptr< TimeStepPdfPackInfo > & pdfPackInfo,
+ const Set<SUID> & requiredBlockSelectors = Set<SUID>::emptySet(),
+ const Set<SUID> & incompatibleBlockSelectors = Set<SUID>::emptySet() )
+{
+ typedef TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T > TS_T;
+ return shared_ptr< TS_T >( new TS_T( blocks, sweep, pdfFieldId, boundaryHandlingId, pdfPackInfo, requiredBlockSelectors, incompatibleBlockSelectors ) );
+}
+
+
+
+} // namespace refinement
+} // namespace lbm
+} // namespace walberla
diff --git a/src/lbm/refinement_rebase/TimeStepPdfPackInfo.h b/src/lbm/refinement_rebase/TimeStepPdfPackInfo.h
new file mode 100644
index 000000000..d4749b882
--- /dev/null
+++ b/src/lbm/refinement_rebase/TimeStepPdfPackInfo.h
@@ -0,0 +1,44 @@
+//======================================================================================================================
+//
+// 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 TimeStep.h
+//! \ingroup lbm
+//! \author Christian Godenschwager <christian.godenschwager@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "blockforest/communication/NonUniformPackInfo.h"
+
+namespace walberla {
+namespace lbm {
+namespace refinement {
+
+
+class TimeStepPdfPackInfo : public blockforest::communication::NonUniformPackInfo
+{
+public:
+ virtual bool optimizedEqualLevelCommunication() const = 0;
+ virtual void optimizeEqualLevelCommunication( const bool value = true ) = 0;
+
+ virtual bool optimizedForLinearExplosion() const = 0;
+ virtual void optimizeForLinearExplosion( const bool value = true ) = 0;
+};
+
+
+} // namespace refinement
+} // namespace lbm
+} // namespace walberla
diff --git a/src/lbm/refinement_rebase/TimeTracker.h b/src/lbm/refinement_rebase/TimeTracker.h
new file mode 100644
index 000000000..3cba8dd52
--- /dev/null
+++ b/src/lbm/refinement_rebase/TimeTracker.h
@@ -0,0 +1,109 @@
+//======================================================================================================================
+//
+// 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 BlockSweepWrapper.h
+//! \ingroup lbm
+//! \author Florian Schornbaum <florian.schornbaum@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/DataTypes.h"
+#include "core/debug/Debug.h"
+
+
+namespace walberla {
+namespace lbm {
+
+
+
+/// Keeps track of the simulation time (One step on the coarsest level counts as '1', one step on the
+/// next finer level counts as '0.5', one step on the next level as '0.25' etc.)
+class TimeTracker
+{
+public:
+
+ TimeTracker( const uint_t initialTime = uint_t(0) )
+ {
+ fractionCounter_.push_back( initialTime );
+ fractionBase_.push_back( uint_t(1) );
+ time_.push_back( real_c( initialTime ) );
+ }
+
+ // for use without refinement time step
+ void operator()()
+ {
+ WALBERLA_ASSERT_GREATER( fractionCounter_.size(), 0 );
+ WALBERLA_ASSERT_GREATER( fractionBase_.size(), 0 );
+ WALBERLA_ASSERT_GREATER( time_.size(), 0 );
+ fractionCounter_[0] += uint_t(1);
+ time_[0] = real_c( fractionCounter_[0] );
+ }
+
+ // for use with refinement time step (addPost*VoidFunction)
+#ifdef NDEBUG
+ void operator()( const uint_t level, const uint_t )
+#else
+ void operator()( const uint_t level, const uint_t executionCount )
+#endif
+ {
+ checkLevel( level );
+
+ WALBERLA_ASSERT_EQUAL( fractionCounter_[ level ] % fractionBase_[ level ], executionCount );
+ fractionCounter_[ level ] += uint_t(1);
+
+ time_[ level ] = real_c( fractionCounter_[ level ] ) / real_c( fractionBase_[ level ] );
+ }
+
+ real_t getTime( const uint_t level = uint_t(0) )
+ {
+ checkLevel( level );
+
+ return time_[ level ];
+ }
+
+private:
+
+ void checkLevel( const uint_t level )
+ {
+ if( level >= time_.size() )
+ {
+ const uint_t previousSize = time_.size();
+
+ fractionCounter_.resize( level + uint_t(1) );
+ fractionBase_.resize( level + uint_t(1) );
+ time_.resize( level + uint_t(1) );
+
+ for( uint_t i = previousSize; i <= level; ++i )
+ {
+ fractionCounter_[ i ] = fractionCounter_[0] * math::uintPow2( i );
+ fractionBase_[ i ] = math::uintPow2( i );
+ time_[ i ] = time_[0];
+ }
+ }
+ }
+
+ std::vector< uint_t > fractionCounter_;
+ std::vector< uint_t > fractionBase_;
+
+ std::vector< real_t > time_;
+
+}; // class TimeTracker
+
+
+
+} // namespace lbm
+} // namespace walberla
diff --git a/src/lbm/refinement_rebase/VorticityBasedLevelDetermination.h b/src/lbm/refinement_rebase/VorticityBasedLevelDetermination.h
new file mode 100644
index 000000000..c3f4ee20d
--- /dev/null
+++ b/src/lbm/refinement_rebase/VorticityBasedLevelDetermination.h
@@ -0,0 +1,161 @@
+//======================================================================================================================
+//
+// 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 VorticityBasedLevelDetermination.h
+//! \ingroup lbm
+//! \author Florian Schornbaum <florian.schornbaum@fau.de>
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "blockforest/BlockForest.h"
+#include "core/math/Vector3.h"
+#include "domain_decomposition/BlockDataID.h"
+#include "field/GhostLayerField.h"
+
+#include <vector>
+
+namespace walberla {
+namespace lbm {
+namespace refinement {
+
+
+/*!\brief Level determination for refinement check based on (scaled) vorticity values
+ *
+ * If (scaled) vorticity magnitude is below lowerLimit in all cells of a block, that block could be coarsened.
+ * If the (scaled) vorticity value is above the upperLimit for at least one cell, that block gets marked for refinement.
+ * Else, the block remains on the current level.
+ *
+ * The scaling originates from neglecting the actual mesh size on the block to obtain different vorticity values for
+ * different mesh sizes.
+ */
+template< typename Filter_T, bool Pseudo2D = false >
+class VorticityBasedLevelDetermination // used as a 'BlockForest::RefreshMinTargetLevelDeterminationFunction'
+{
+
+public:
+
+ typedef GhostLayerField< Vector3<real_t>, 1 > VectorField_T;
+
+ VorticityBasedLevelDetermination( const ConstBlockDataID & fieldID, const Filter_T & filter,
+ const real_t upperLimit, const real_t lowerLimit, const uint_t maxLevel ) :
+ fieldID_( fieldID ), filter_( filter ),
+ upperLimit_( upperLimit * upperLimit ), lowerLimit_( lowerLimit * lowerLimit ), maxLevel_( maxLevel )
+ {}
+
+ void operator()( std::vector< std::pair< const Block *, uint_t > > & minTargetLevels,
+ std::vector< const Block * > & blocksAlreadyMarkedForRefinement,
+ const BlockForest & forest );
+
+private:
+
+ ConstBlockDataID fieldID_;
+
+ Filter_T filter_;
+
+ real_t upperLimit_;
+ real_t lowerLimit_;
+
+ uint_t maxLevel_;
+
+}; // class VorticityRefinement
+
+template< typename Filter_T, bool Pseudo2D >
+void VorticityBasedLevelDetermination< Filter_T, Pseudo2D >::operator()( std::vector< std::pair< const Block *, uint_t > > & minTargetLevels,
+ std::vector< const Block * > &, const BlockForest & )
+{
+ for( auto it = minTargetLevels.begin(); it != minTargetLevels.end(); ++it )
+ {
+ const Block * const block = it->first;
+ const VectorField_T * u = block->template getData< VectorField_T >( fieldID_ );
+
+ if( u == nullptr )
+ {
+ it->second = uint_t(0);
+ continue;
+ }
+
+ WALBERLA_ASSERT_GREATER_EQUAL(u->nrOfGhostLayers(), uint_t(1));
+
+ CellInterval interval = u->xyzSize();
+ Cell expand( cell_idx_c(-1), cell_idx_c(-1), Pseudo2D ? cell_idx_t(0) : cell_idx_c(-1) );
+ interval.expand( expand );
+
+ const cell_idx_t one( cell_idx_t(1) );
+ const real_t half( real_c(0.5) );
+
+ bool refine( false );
+ bool coarsen( true );
+
+ filter_( *block );
+
+ const cell_idx_t xSize = cell_idx_c( interval.xSize() );
+ const cell_idx_t ySize = cell_idx_c( interval.ySize() );
+ const cell_idx_t zSize = cell_idx_c( interval.zSize() );
+
+ for( cell_idx_t z = cell_idx_t(0); z < zSize; ++z ) {
+ for (cell_idx_t y = cell_idx_t(0); y < ySize; ++y) {
+ for (cell_idx_t x = cell_idx_t(0); x < xSize; ++x) {
+ if( filter_(x,y,z) && filter_(x+one,y,z) && filter_(x-one,y,z) && filter_(x,y+one,z) && filter_(x,y-one,z) &&
+ ( Pseudo2D || (filter_(x,y,z+one) && filter_(x,y,z-one)) ) )
+ {
+ const Vector3< real_t > xa = u->get(x+one,y,z);
+ const Vector3< real_t > xb = u->get(x-one,y,z);
+ const Vector3< real_t > ya = u->get(x,y+one,z);
+ const Vector3< real_t > yb = u->get(x,y-one,z);
+ const Vector3< real_t > za = Pseudo2D ? Vector3< real_t >(0) : u->get(x,y,z+one);
+ const Vector3< real_t > zb = Pseudo2D ? Vector3< real_t >(0) : u->get(x,y,z-one);
+
+ const real_t duzdy = half * ( ya[2] - yb[2] );
+ const real_t duydz = half * ( za[1] - zb[1] );
+ const real_t duxdz = half * ( za[0] - zb[0] );
+ const real_t duzdx = half * ( xa[2] - xb[2] );
+ const real_t duydx = half * ( xa[1] - xb[1] );
+ const real_t duxdy = half * ( ya[0] - yb[0] );
+
+ // since no dx was used in the gradient calculation, this value is actually curl * dx
+ // this is needed to have changing curl values on different grid levels
+ const Vector3< real_t > curl( duzdy - duydz, duxdz - duzdx, duydx - duxdy );
+ const auto curlSqr = curl.sqrLength();
+
+ if( curlSqr > lowerLimit_ )
+ {
+ coarsen = false;
+ if( curlSqr > upperLimit_ )
+ refine = true;
+ }
+ }
+ }
+ }
+ }
+
+ if( refine && block->getLevel() < maxLevel_ )
+ {
+ WALBERLA_ASSERT( !coarsen );
+ it->second = block->getLevel() + uint_t(1);
+ }
+ if( coarsen && block->getLevel() > uint_t(0) )
+ {
+ WALBERLA_ASSERT( !refine );
+ it->second = block->getLevel() - uint_t(1);
+ }
+ }
+}
+
+} // namespace refinement
+} // namespace lbm
+} // namespace walberla
diff --git a/src/lbm/refinement_rebase/all.h b/src/lbm/refinement_rebase/all.h
new file mode 100644
index 000000000..a63a1ba38
--- /dev/null
+++ b/src/lbm/refinement_rebase/all.h
@@ -0,0 +1,32 @@
+//======================================================================================================================
+//
+// 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 all.h
+//! \ingroup lbm
+//! \author Christian Godenschwager <christian.godenschwager@fau.de>
+//! \brief Collective header file for module lbm
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "BoundarySetup.h"
+#include "PdfFieldPackInfo.h"
+#include "PdfFieldSyncPackInfo.h"
+#include "RefinementFunctorWrapper.h"
+#include "TimeStep.h"
+#include "TimeStepPdfPackInfo.h"
+#include "TimeTracker.h"
+#include "VorticityBasedLevelDetermination.h"
diff --git a/tests/lbm/CMakeLists.txt b/tests/lbm/CMakeLists.txt
index f66bd62d0..4cc28ab9a 100644
--- a/tests/lbm/CMakeLists.txt
+++ b/tests/lbm/CMakeLists.txt
@@ -92,7 +92,8 @@ waLBerla_generate_target_from_python(NAME FieldLayoutAndVectorizationTestGenerat
waLBerla_compile_test( FILES codegen/FieldLayoutAndVectorizationTest.cpp DEPENDS FieldLayoutAndVectorizationTestGenerated)
-set( configs D3Q19_SRT_INCOMP D3Q19_SRT_COMP D3Q27_SRT_INCOMP D3Q27_SRT_COMP D3Q27_TRT_INCOMP D3Q19_MRT_COMP )
+# set( configs D3Q19_SRT_INCOMP D3Q19_SRT_COMP D3Q27_SRT_INCOMP D3Q27_SRT_COMP D3Q27_TRT_INCOMP D3Q19_MRT_COMP )
+set( configs D3Q19_SRT_INCOMP )
list( APPEND RefinementTargets "" )
foreach( config ${configs} )
diff --git a/tests/lbm/codegen/CodeGenerationRefinement.cpp b/tests/lbm/codegen/CodeGenerationRefinement.cpp
index 47508ab70..fa728c748 100644
--- a/tests/lbm/codegen/CodeGenerationRefinement.cpp
+++ b/tests/lbm/codegen/CodeGenerationRefinement.cpp
@@ -42,19 +42,20 @@
#include <lbm/field/PdfField.h>
#include <lbm/lattice_model/D3Q19.h>
#include <lbm/lattice_model/D3Q27.h>
-#include <lbm/refinement/TimeStep.h>
#include <lbm/sweeps/CellwiseSweep.h>
#include <timeloop/SweepTimeloop.h>
#include <type_traits>
-#include "CodeGenerationRefinement_D3Q19_SRT_COMP_LatticeModel.h"
+#include <lbm/refinement_rebase/TimeStep.h>
+
+// #include "CodeGenerationRefinement_D3Q19_SRT_COMP_LatticeModel.h"
#include "CodeGenerationRefinement_D3Q19_SRT_INCOMP_LatticeModel.h"
-#include "CodeGenerationRefinement_D3Q27_SRT_COMP_LatticeModel.h"
-#include "CodeGenerationRefinement_D3Q27_SRT_INCOMP_LatticeModel.h"
-#include "CodeGenerationRefinement_D3Q27_TRT_INCOMP_LatticeModel.h"
-#include "CodeGenerationRefinement_D3Q19_MRT_COMP_LatticeModel.h"
+// #include "CodeGenerationRefinement_D3Q27_SRT_COMP_LatticeModel.h"
+// #include "CodeGenerationRefinement_D3Q27_SRT_INCOMP_LatticeModel.h"
+// #include "CodeGenerationRefinement_D3Q27_TRT_INCOMP_LatticeModel.h"
+// #include "CodeGenerationRefinement_D3Q19_MRT_COMP_LatticeModel.h"
-//#define TEST_USES_VTK_OUTPUT
+#define TEST_USES_VTK_OUTPUT
#ifdef TEST_USES_VTK_OUTPUT
#include "lbm/vtk/Density.h"
#include "lbm/vtk/Velocity.h"
@@ -503,13 +504,13 @@ int main( int argc, char ** argv )
fieldIds.emplace_back();
- typedef lbm::D3Q19<SRT, true> walberalLM_D3Q19_SRT_COMP;
- typedef lbm::CodeGenerationRefinement_D3Q19_SRT_COMP_LatticeModel lbmpyLM_D3Q19_SRT_COMP;
-
- AddRefinementTimeStep< walberalLM_D3Q19_SRT_COMP >::add(blocks, timeloop, walberalLM_D3Q19_SRT_COMP(SRT(GlobalOmega)),
- fieldIds.back(), field::fzyx, flagFieldId, velocity, "( waLBerla D3Q19 SRT COMP )");
- AddRefinementTimeStep< lbmpyLM_D3Q19_SRT_COMP >::add(blocks, timeloop, lbmpyLM_D3Q19_SRT_COMP(GlobalOmega),
- fieldIds.back(), field::fzyx, flagFieldId, velocity, "( lbmpy D3Q19 SRT COMP )");
+// typedef lbm::D3Q19<SRT, true> walberalLM_D3Q19_SRT_COMP;
+// typedef lbm::CodeGenerationRefinement_D3Q19_SRT_COMP_LatticeModel lbmpyLM_D3Q19_SRT_COMP;
+//
+// AddRefinementTimeStep< walberalLM_D3Q19_SRT_COMP >::add(blocks, timeloop, walberalLM_D3Q19_SRT_COMP(SRT(GlobalOmega)),
+// fieldIds.back(), field::fzyx, flagFieldId, velocity, "( waLBerla D3Q19 SRT COMP )");
+// AddRefinementTimeStep< lbmpyLM_D3Q19_SRT_COMP >::add(blocks, timeloop, lbmpyLM_D3Q19_SRT_COMP(GlobalOmega),
+// fieldIds.back(), field::fzyx, flagFieldId, velocity, "( lbmpy D3Q19 SRT COMP )");
////////////////////////////////
/// D3Q27 SRT incompressible ///
@@ -517,13 +518,13 @@ int main( int argc, char ** argv )
fieldIds.emplace_back();
- typedef lbm::D3Q27<SRT, false> walberalLM_D3Q27_SRT_INCOMP;
- typedef lbm::CodeGenerationRefinement_D3Q27_SRT_INCOMP_LatticeModel lbmpyLM_D3Q27_SRT_INCOMP;
-
- AddRefinementTimeStep< walberalLM_D3Q27_SRT_INCOMP >::add(blocks, timeloop, walberalLM_D3Q27_SRT_INCOMP(SRT(GlobalOmega)),
- fieldIds.back(), field::fzyx, flagFieldId, velocity, "( waLBerla D3Q27 SRT INCOMP )");
- AddRefinementTimeStep< lbmpyLM_D3Q27_SRT_INCOMP >::add(blocks, timeloop, lbmpyLM_D3Q27_SRT_INCOMP(GlobalOmega),
- fieldIds.back(), field::fzyx, flagFieldId, velocity, "( lbmpy D3Q27 SRT INCOMP )");
+// typedef lbm::D3Q27<SRT, false> walberalLM_D3Q27_SRT_INCOMP;
+// typedef lbm::CodeGenerationRefinement_D3Q27_SRT_INCOMP_LatticeModel lbmpyLM_D3Q27_SRT_INCOMP;
+//
+// AddRefinementTimeStep< walberalLM_D3Q27_SRT_INCOMP >::add(blocks, timeloop, walberalLM_D3Q27_SRT_INCOMP(SRT(GlobalOmega)),
+// fieldIds.back(), field::fzyx, flagFieldId, velocity, "( waLBerla D3Q27 SRT INCOMP )");
+// AddRefinementTimeStep< lbmpyLM_D3Q27_SRT_INCOMP >::add(blocks, timeloop, lbmpyLM_D3Q27_SRT_INCOMP(GlobalOmega),
+// fieldIds.back(), field::fzyx, flagFieldId, velocity, "( lbmpy D3Q27 SRT INCOMP )");
////////////////////////////////
/// D3Q27 SRT compressible ///
@@ -531,13 +532,13 @@ int main( int argc, char ** argv )
fieldIds.emplace_back();
- typedef lbm::D3Q27<SRT, true> walberalLM_D3Q27_SRT_COMP;
- typedef lbm::CodeGenerationRefinement_D3Q27_SRT_COMP_LatticeModel lbmpyLM_D3Q27_SRT_COMP;
-
- AddRefinementTimeStep< walberalLM_D3Q27_SRT_COMP >::add(blocks, timeloop, walberalLM_D3Q27_SRT_COMP(SRT(GlobalOmega)),
- fieldIds.back(), field::fzyx, flagFieldId, velocity, "( waLBerla D3Q27 SRT COMP )");
- AddRefinementTimeStep< lbmpyLM_D3Q27_SRT_COMP >::add(blocks, timeloop, lbmpyLM_D3Q27_SRT_COMP(GlobalOmega),
- fieldIds.back(), field::fzyx, flagFieldId, velocity, "( lbmpy D3Q27 SRT COMP )");
+// typedef lbm::D3Q27<SRT, true> walberalLM_D3Q27_SRT_COMP;
+// typedef lbm::CodeGenerationRefinement_D3Q27_SRT_COMP_LatticeModel lbmpyLM_D3Q27_SRT_COMP;
+//
+// AddRefinementTimeStep< walberalLM_D3Q27_SRT_COMP >::add(blocks, timeloop, walberalLM_D3Q27_SRT_COMP(SRT(GlobalOmega)),
+// fieldIds.back(), field::fzyx, flagFieldId, velocity, "( waLBerla D3Q27 SRT COMP )");
+// AddRefinementTimeStep< lbmpyLM_D3Q27_SRT_COMP >::add(blocks, timeloop, lbmpyLM_D3Q27_SRT_COMP(GlobalOmega),
+// fieldIds.back(), field::fzyx, flagFieldId, velocity, "( lbmpy D3Q27 SRT COMP )");
/////////////////////
/// RUN TEST RUNS ///
@@ -554,14 +555,14 @@ int main( int argc, char ** argv )
//////////////////
check<walberalLM_D3Q19_SRT_INCOMP, lbmpyLM_D3Q19_SRT_INCOMP>( blocks, fieldIds[0][0], fieldIds[0][1] );
- check<walberalLM_D3Q19_SRT_COMP, lbmpyLM_D3Q19_SRT_COMP> ( blocks, fieldIds[1][0], fieldIds[1][1] );
+ // check<walberalLM_D3Q19_SRT_COMP, lbmpyLM_D3Q19_SRT_COMP> ( blocks, fieldIds[1][0], fieldIds[1][1] );
//////////////////
/// D3Q27, SRT ///
//////////////////
- check<walberalLM_D3Q27_SRT_INCOMP, lbmpyLM_D3Q27_SRT_INCOMP>( blocks, fieldIds[2][0], fieldIds[2][1] );
- check<walberalLM_D3Q27_SRT_COMP, lbmpyLM_D3Q27_SRT_COMP> ( blocks, fieldIds[3][0], fieldIds[3][1] );
+ // check<walberalLM_D3Q27_SRT_INCOMP, lbmpyLM_D3Q27_SRT_INCOMP>( blocks, fieldIds[2][0], fieldIds[2][1] );
+ // check<walberalLM_D3Q27_SRT_COMP, lbmpyLM_D3Q27_SRT_COMP> ( blocks, fieldIds[3][0], fieldIds[3][1] );
return EXIT_SUCCESS;
}
--
GitLab