diff --git a/src/field/all.h b/src/field/all.h index 12a69dfe7ba2ba592cd30ccb6636ffcd93ce422d..05ea39cc7534d32cc00dd84dca8cd72ded61119a 100644 --- a/src/field/all.h +++ b/src/field/all.h @@ -45,6 +45,7 @@ #include "allocation/all.h" #include "blockforest/all.h" #include "communication/all.h" +#include "distributors/all.h" #include "interpolators/all.h" #include "iterators/all.h" #include "refinement/all.h" diff --git a/src/field/distributors/DistributorCreators.h b/src/field/distributors/DistributorCreators.h new file mode 100644 index 0000000000000000000000000000000000000000..78d790c626ffff7acbd354cec9fab8ae4a7b7e02 --- /dev/null +++ b/src/field/distributors/DistributorCreators.h @@ -0,0 +1,108 @@ +//====================================================================================================================== +// +// 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 DistributorCreators.h +//! \ingroup field +//! \author Christoph Rettinger <christoph.rettinger@fau.de> +// +//====================================================================================================================== + +#pragma once + +#include "core/DataTypes.h" +#include "core/Set.h" + +#include "blockforest/BlockDataHandling.h" + +#include "domain_decomposition/StructuredBlockStorage.h" + +namespace walberla { +namespace field { + +//********************************************************************************************************************** +/*! DistributorCreators + * + * \ingroup field + * + * Distributor_T: A distributor that has a constructor + * ( const shared_ptr<StructuredBlockStorage> & blockStorage, const IBlock & block, const BaseField_T & baseField, + * const FlagField_T & flagField, const flag_t & evaluationMask ) + * and distribution functions: + * template< typename ForwardIterator_T > inline void distribute( const Vector3<real_t> & position, ForwardIterator_T distributeValueBegin ) + * template< typename ForwardIterator_T > inline void distribute( const real_t & x, const real_t & y, const real_t & z, ForwardIterator_T distributeValueBegin ) + * + * See NearestNeighborDistributor for an example implementation. + * + * A distributor is aware of the flag field (FlagField_T) and distributes values only to cells flagged by a given mask. + * + * Distributors are used to spread a given value to the corresponding destination field. + * E.g. if a certain force has to be applied at some specific position onto the fluid, a distributor can be used + * to do so by distributing this force value (and conservation fo this force value is ensured) onto the force field. + * + */ +//********************************************************************************************************************** +template< typename Distributor_T, typename FlagField_T > +class DistributorHandling : public blockforest::AlwaysInitializeBlockDataHandling< Distributor_T > +{ +public: + + DistributorHandling( const shared_ptr<StructuredBlockStorage> & blockStorage, + const BlockDataID & distributionDestinationFieldID, + const ConstBlockDataID & flagFieldID, + const Set< FlagUID > & cellsToEvaluate ) : + blockStorage_( blockStorage ), distributionDestinationFieldID_( distributionDestinationFieldID ), flagFieldID_( flagFieldID ), cellsToEvaluate_( cellsToEvaluate ) + {} + + Distributor_T * initialize( IBlock * const block ) + { + typedef typename Distributor_T::BaseField_T DistributionDestinationField_T; + typedef typename FlagField_T::flag_t flag_t; + DistributionDestinationField_T * distributionDestinationField = block->getData< DistributionDestinationField_T >( distributionDestinationFieldID_ ); + const FlagField_T * flagField = block->getData< FlagField_T >( flagFieldID_ ); + + WALBERLA_ASSERT_NOT_NULLPTR( distributionDestinationField ); + WALBERLA_ASSERT_NOT_NULLPTR( flagField ); + + const flag_t evaluationMask = flagField->getMask( cellsToEvaluate_ ); + + return new Distributor_T( blockStorage_, *block, *distributionDestinationField, *flagField, evaluationMask ); + } + +private: + + shared_ptr<StructuredBlockStorage> blockStorage_; + BlockDataID distributionDestinationFieldID_; + ConstBlockDataID flagFieldID_; + Set< FlagUID > cellsToEvaluate_; + +}; // class DistributorHandling + + + +template< typename Distributor_T, typename FlagField_T > +inline BlockDataID addDistributor( const shared_ptr< StructuredBlockStorage > & blocks, + const BlockDataID & distributionDestinationFieldID, + const ConstBlockDataID & flagFieldID, + const Set< FlagUID > & cellsToEvaluate, + const std::string & identifier = std::string(), + const Set<SUID> & requiredSelectors = Set<SUID>::emptySet(), + const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet() ) +{ + return blocks->addBlockData( make_shared< DistributorHandling< Distributor_T, FlagField_T > >( blocks, distributionDestinationFieldID, flagFieldID, cellsToEvaluate ), identifier, requiredSelectors, incompatibleSelectors ); +} + + +} // namespace field +} // namespace walberla diff --git a/src/field/distributors/KernelDistributor.h b/src/field/distributors/KernelDistributor.h new file mode 100644 index 0000000000000000000000000000000000000000..e08369330afbf0ecc2abe476f89394e2494dbf4f --- /dev/null +++ b/src/field/distributors/KernelDistributor.h @@ -0,0 +1,155 @@ +//====================================================================================================================== +// +// 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 KernelDistributor.h +//! \ingroup field +//! \author Christoph Rettinger <christoph.rettinger@fau.de> +// +//====================================================================================================================== + +#pragma once + +#include "core/debug/Debug.h" +#include "core/math/Vector3.h" + +#include "domain_decomposition/StructuredBlockStorage.h" + +#include "field/interpolators/KernelFieldInterpolator.h" +#include "field/GhostLayerField.h" + +#include <vector> + +namespace walberla { +namespace field { + +/*! Distributor for walberla::field::GhostLayerField with kernel strategy + * + * \ingroup field + * + * This distributor uses a smoothed dirac kernel function to distribute values to the field at the given position. + * The applied weights are given in the namespace field::kernelweights. + * Needs the full neighborhood of the containing cell, i.e. 27 cells. + * Never construct this distributor directly, but use the functionality from DistributorCreators.h instead. + */ +template< typename Field_T, typename FlagField_T > +class KernelDistributor +{ +public: + + static const uint_t F_SIZE = Field_T::F_SIZE; + + typedef Field_T BaseField_T; + typedef typename FlagField_T::flag_t flag_t; + typedef KernelDistributor<Field_T,FlagField_T> OwnType; + + KernelDistributor( const shared_ptr<StructuredBlockStorage> & blockStorage, const IBlock & block, + BaseField_T & baseField, const FlagField_T & flagField, + const flag_t & evaluationMask ) + : blockStorage_( blockStorage ), block_( block ), baseField_( baseField ), flagField_( flagField ), evaluationMask_( evaluationMask ) + { + WALBERLA_ASSERT(baseField.nrOfGhostLayers() > uint_t(0), "field for kernel distribution needs at least one ghost layer"); + } + + inline bool operator==( const OwnType & other ){ return baseField_ == other.baseField_; } + + template< typename ForwardIterator_T > + inline void distribute( const Vector3<real_t> & position, ForwardIterator_T distributeValueBegin ) + { + distribute( position[0], position[1], position[2], distributeValueBegin ); + } + + template< typename ForwardIterator_T > + inline void distribute( const real_t & x, const real_t & y, const real_t & z, ForwardIterator_T distributeValueBegin ) + { + WALBERLA_ASSERT(block_.getAABB().contains(x,y,z), + "Distribution position <" << x << ", " << y << ", " << z << "> is not contained inside the block of this distributor with AABB " << block_.getAABB() << " !"); + + Cell centerCell = blockStorage_->getBlockLocalCell( block_, x, y, z ); + + const real_t dx = blockStorage_->dx( blockStorage_->getLevel( block_ ) ); + const real_t dy = blockStorage_->dy( blockStorage_->getLevel( block_ ) ); + const real_t dz = blockStorage_->dz( blockStorage_->getLevel( block_ ) ); + + const uint_t neighborhoodSize = cell_idx_t(1); + + CellInterval cellNeighborhood( centerCell[0] - cell_idx_c(neighborhoodSize), centerCell[1] - cell_idx_c(neighborhoodSize), centerCell[2] - cell_idx_c(neighborhoodSize), + centerCell[0] + cell_idx_c(neighborhoodSize), centerCell[1] + cell_idx_c(neighborhoodSize), centerCell[2] + cell_idx_c(neighborhoodSize) ); + + const uint_t kernelSizeOneDirection = uint_t(2) * neighborhoodSize + uint_t(1); + std::vector<real_t> weights( kernelSizeOneDirection * kernelSizeOneDirection * kernelSizeOneDirection, real_t(0) ); + + uint_t counter = uint_t(0); + real_t sumOfWeights = real_t(0); + real_t sumOfWeightsUnavailable = real_t(0); + + // get distribution weights and count available cells in surrounding cells + for( auto cellIt = cellNeighborhood.begin(); cellIt != cellNeighborhood.end(); ++cellIt ) + { + Vector3<real_t> curCellCenter = blockStorage_->getBlockLocalCellCenter( block_, *cellIt ); + if( flagField_.isPartOfMaskSet( *cellIt, evaluationMask_ ) ) + { + weights[counter] = kernelweights::kernelWeightFunction( x, y, z, curCellCenter[0], curCellCenter[1], curCellCenter[2], dx, dy, dz ); + sumOfWeights += weights[counter]; + } + else + { + weights[counter] = real_t(0); + sumOfWeightsUnavailable += kernelweights::kernelWeightFunction( x, y, z, curCellCenter[0], curCellCenter[1], curCellCenter[2], dx, dy, dz ); + } + ++counter; + } + + // check if at least one cell was available, to prevent division by 0 + if( sumOfWeights <= real_t(0) ) + return; + + // scale domain weights if some non-domain cells are in neighborhood + const real_t scalingFactor = real_t(1) + sumOfWeightsUnavailable / sumOfWeights; + + // distribute the values to the neighboring domain cells with the corresponding (scaled) weighting + counter = uint_t(0); + for( auto cellIt = cellNeighborhood.begin(); cellIt != cellNeighborhood.end(); ++cellIt ) + { + if ( weights[counter] > real_t(0) ) + { + addWeightedCellValue( distributeValueBegin, *cellIt, scalingFactor * weights[counter] ); + } + ++counter; + } + } + +private: + + template< typename ForwardIterator_T > + void addWeightedCellValue( ForwardIterator_T distributeValueBegin, const Cell & curCell, const real_t & weighting ) + { + for( uint_t f = uint_t(0); f < F_SIZE; ++f ) + { + baseField_( curCell, f) += weighting * (*distributeValueBegin); + ++distributeValueBegin; + } + } + + shared_ptr<StructuredBlockStorage> blockStorage_; + const IBlock & block_; + BaseField_T & baseField_; + const FlagField_T & flagField_; + flag_t evaluationMask_; + +}; + + +} // namespace field +} // namespace walberla diff --git a/src/field/distributors/NearestNeighborDistributor.h b/src/field/distributors/NearestNeighborDistributor.h new file mode 100644 index 0000000000000000000000000000000000000000..7c8f57092cabe0daf61fe57bc589d6728624dfc3 --- /dev/null +++ b/src/field/distributors/NearestNeighborDistributor.h @@ -0,0 +1,138 @@ +//====================================================================================================================== +// +// 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 NearestNeighborDistributor.h +//! \ingroup field +//! \author Christoph Rettinger <christoph.rettinger@fau.de> +// +//====================================================================================================================== + +#pragma once + +#include "core/math/Vector3.h" + +#include "domain_decomposition/StructuredBlockStorage.h" + +#include "field/GhostLayerField.h" + +#include <numeric> +#include <vector> + +namespace walberla { +namespace field { + +/*! Distributor for walberla::field::Field with nearest neighbor strategy + * + * \ingroup field + * + * This distributor distributes the given value to the nearest cell, flagged in the evaluation mask, of the given position. + * This is usually the cell that contains the distribution position. + * If this cell is a cell that is not within the evaluation mask, the direct neighborhood (8 cells) will be searched for an available cell. + * Never construct this distributor directly, but use the functionality from DistributorCreators.h instead. + */ +template< typename Field_T, typename FlagField_T > +class NearestNeighborDistributor +{ +public: + + static const uint_t F_SIZE = Field_T::F_SIZE; + + typedef Field_T BaseField_T; + typedef typename FlagField_T::flag_t flag_t; + typedef NearestNeighborDistributor<Field_T,FlagField_T> OwnType; + + NearestNeighborDistributor( const shared_ptr<StructuredBlockStorage> & blockStorage, const IBlock & block, + BaseField_T & baseField, const FlagField_T & flagField, + const flag_t & evaluationMask ) + : blockStorage_( blockStorage ), block_( block ), baseField_( baseField ), flagField_( flagField ), evaluationMask_( evaluationMask ) + {} + + inline bool operator==( const OwnType & other ){ return baseField_ == other.baseField_; } + + template< typename ForwardIterator_T > + inline void distribute( const Vector3<real_t> & position, ForwardIterator_T distributeValueBegin ) + { + distribute( position[0], position[1], position[2], distributeValueBegin ); + } + + template< typename ForwardIterator_T > + inline void distribute( const real_t & x, const real_t & y, const real_t & z, ForwardIterator_T distributeValueBegin ) + { + WALBERLA_ASSERT(block_.getAABB().contains(x,y,z), + "Distribution position <" << x << ", " << y << ", " << z << "> is not contained inside the block of this distributor with AABB " << block_.getAABB() << " !"); + + Cell nearestCell = blockStorage_->getBlockLocalCell( block_, x, y, z ); + + if( flagField_.isPartOfMaskSet( nearestCell, evaluationMask_ ) ) + { + for( uint_t f = uint_t(0); f < F_SIZE; ++f ) + { + baseField_( nearestCell, f) += *distributeValueBegin; + ++distributeValueBegin; + } + return; + } + else + { + // look for available cell in direct neighborhood + + CellInterval fieldXYZSize = baseField_.xyzSize(); + + Vector3<real_t> nearestCellCenter = blockStorage_->getBlockLocalCellCenter( block_, nearestCell ); + const cell_idx_t xNeighbor = cell_idx_c( floor( x - nearestCellCenter[0] ) ); + const cell_idx_t yNeighbor = cell_idx_c( floor( y - nearestCellCenter[1] ) ); + const cell_idx_t zNeighbor = cell_idx_c( floor( z - nearestCellCenter[2] ) ); + + const cell_idx_t xMin = nearestCell.x() + xNeighbor; + const cell_idx_t yMin = nearestCell.y() + yNeighbor; + const cell_idx_t zMin = nearestCell.z() + zNeighbor; + + for( cell_idx_t zC = zMin; zC <= zMin + cell_idx_t(1); ++zC) + { + for( cell_idx_t yC = yMin; yC <= yMin + cell_idx_t(1); ++yC) + { + for( cell_idx_t xC = xMin; xC <= xMin + cell_idx_t(1); ++xC) + { + Cell curCell(xC,yC,zC); + if( fieldXYZSize.contains(curCell) ) + { + if (flagField_.isPartOfMaskSet(curCell, evaluationMask_)) + { + for (uint_t f = uint_t(0); f < F_SIZE; ++f) + { + baseField_(curCell, f) += *distributeValueBegin; + ++distributeValueBegin; + } + return; + } + } + } + } + } + } + } + +private: + + shared_ptr<StructuredBlockStorage> blockStorage_; + const IBlock & block_; + BaseField_T & baseField_; + const FlagField_T & flagField_; + flag_t evaluationMask_; + +}; + +} // namespace field +} // namespace walberla diff --git a/src/field/distributors/all.h b/src/field/distributors/all.h new file mode 100644 index 0000000000000000000000000000000000000000..8a31a22b6e70b8fd6d7647a6a0840220f8d02386 --- /dev/null +++ b/src/field/distributors/all.h @@ -0,0 +1,26 @@ +//====================================================================================================================== +// +// 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 field +//! \author Christoph Rettinger <christoph.rettinger@fau.de> +// +//====================================================================================================================== + +#pragma once + +#include "DistributorCreators.h" +#include "KernelDistributor.h" +#include "NearestNeighborDistributor.h" diff --git a/tests/field/CMakeLists.txt b/tests/field/CMakeLists.txt index 2e490ea32b22d0266e8a7aac742be338758057be..d7534a81375d2e85498a6c8fae5b51c0c5b0a68d 100644 --- a/tests/field/CMakeLists.txt +++ b/tests/field/CMakeLists.txt @@ -10,6 +10,9 @@ waLBerla_execute_test( NAME AccuracyEvaluationTest4 COMMAND $<TARGET_FILE:Accura waLBerla_compile_test( FILES communication/FieldPackInfoTest.cpp DEPENDS blockforest ) waLBerla_execute_test( NAME FieldPackInfoTest ) +waLBerla_compile_test( FILES distributors/DistributionTest.cpp) +waLBerla_execute_test( NAME DistributionTest ) + waLBerla_compile_test( FILES FieldTest.cpp ) waLBerla_execute_test( NAME FieldTest ) diff --git a/tests/field/distributors/DistributionTest.cpp b/tests/field/distributors/DistributionTest.cpp new file mode 100644 index 0000000000000000000000000000000000000000..4d0e4f644d153c055c9d61c0648568a13df94eb0 --- /dev/null +++ b/tests/field/distributors/DistributionTest.cpp @@ -0,0 +1,536 @@ +//====================================================================================================================== +// +// 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 DistributionTest.cpp +//! \ingroup field +//! \author Christoph Rettinger <christoph.rettinger@fau.de> +// +//====================================================================================================================== + +#include "blockforest/all.h" + +#include "core/debug/TestSubsystem.h" +#include "core/mpi/Environment.h" +#include "core/math/all.h" + +#include "field/AddToStorage.h" +#include "field/FlagField.h" +#include "field/GhostLayerField.h" +#include "field/distributors/all.h" + +#include <vector> + +using namespace walberla; + +namespace distribution_tests { + +const uint_t FieldGhostLayers( 1 ); + +typedef walberla::uint8_t flag_t; +typedef FlagField< flag_t > FlagField_T; + +typedef GhostLayerField< real_t, 1> ScalarField_T; +typedef GhostLayerField< Vector3<real_t>, 1> Vec3Field_T; +typedef GhostLayerField< real_t, 3> MultiComponentField_T; + + +const FlagUID Domain_Flag ( "domain" ); +const FlagUID Boundary_Flag ( "boundary" ); + +void initFlagField( FlagField_T * field, IBlock * const /*block*/ ) +{ + auto domainFlag = field->getOrRegisterFlag( Domain_Flag ); + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ( field, field->addFlag( x, y, z, domainFlag ); ); +} + +void resetScalarField( const shared_ptr<StructuredBlockStorage> & blocks, + const BlockDataID & scalarFieldID ) +{ + for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) + { + auto sField = blockIt->getData<ScalarField_T>( scalarFieldID ); + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(sField, + sField->get(x,y,z) = real_t(0); + ); + } +} + +void resetVectorField( const shared_ptr<StructuredBlockStorage> & blocks, + const BlockDataID & vectorFieldID ) +{ + for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) + { + auto vField = blockIt->getData<Vec3Field_T>( vectorFieldID ); + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(vField, + vField->get(x,y,z) = Vector3<real_t>(real_t(0)); + ); + } +} + +void resetMultiCompField( const shared_ptr<StructuredBlockStorage> & blocks, + const BlockDataID & multiComponentFieldID ) +{ + for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) + { + auto mField = blockIt->getData<MultiComponentField_T>( multiComponentFieldID ); + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(mField, + mField->get(x,y,z,0) = real_t(0); + mField->get(x,y,z,1) = real_t(0); + mField->get(x,y,z,2) = real_t(0); + ); + } +} + + +void setBoundaryFlags( const shared_ptr<StructuredBlockStorage> & blocks, + const BlockDataID & flagFieldID ) +{ + for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) + { + auto flagField = blockIt->getData<FlagField_T>( flagFieldID ); + auto domainFlag = flagField->getOrRegisterFlag( Domain_Flag ); + auto boundaryFlag = flagField->getOrRegisterFlag( Boundary_Flag ); + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, + if( x == 2) + { + flagField->removeFlag(x,y,z,domainFlag); + flagField->addFlag(x,y,z,boundaryFlag); + } + ); + } +} + +void getScalarFieldQuantities( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & flagFieldID, + const BlockDataID & fieldID, + real_t & summarizedValue, real_t & minValue, real_t & maxValue) +{ + real_t sum( real_t(0) ); + real_t min( real_t(0) ); + real_t max( real_t(0) ); + for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) + { + auto field = blockIt->getData<ScalarField_T>( fieldID ); + auto flagField = blockIt->getData<FlagField_T>( flagFieldID ); + auto domainFlag = flagField->getFlag( Domain_Flag ); + + CellInterval xyzSizeWithGhostLayers = field->xyzSizeWithGhostLayer(); + for( auto cellIt = xyzSizeWithGhostLayers.begin(); cellIt != xyzSizeWithGhostLayers.end(); ++cellIt ) + { + if( flagField->isFlagSet(*cellIt,domainFlag)) + { + real_t value = field->get(*cellIt); + sum += value; + min = std::min(min, value); + max = std::max(max, value); + } + } + } + WALBERLA_MPI_SECTION() + { + mpi::allReduceInplace( sum, mpi::SUM ); + mpi::allReduceInplace( min, mpi::MIN ); + mpi::allReduceInplace( max, mpi::MAX ); + } + + summarizedValue = sum; + minValue = min; + maxValue = max; +} + +void getVectorFieldQuantities( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & flagFieldID, + const BlockDataID & fieldID, + Vector3<real_t> & summarizedValue, Vector3<real_t> & minValue, Vector3<real_t> & maxValue) +{ + Vector3<real_t> sum( real_t(0) ); + Vector3<real_t> min( real_t(0) ); + Vector3<real_t> max( real_t(0) ); + for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) + { + auto field = blockIt->getData<Vec3Field_T>( fieldID ); + auto flagField = blockIt->getData<FlagField_T>( flagFieldID ); + auto domainFlag = flagField->getFlag( Domain_Flag ); + CellInterval xyzSizeWithGhostLayers = field->xyzSizeWithGhostLayer(); + for( auto cellIt = xyzSizeWithGhostLayers.begin(); cellIt != xyzSizeWithGhostLayers.end(); ++cellIt ) + { + if( flagField->isFlagSet(*cellIt,domainFlag)) + { + Vector3<real_t> value = field->get(*cellIt); + sum += value; + for (size_t i = 0; i < 3; ++i) { + min[i] = std::min(min[i], value[i]); + max[i] = std::max(max[i], value[i]); + } + } + } + } + WALBERLA_MPI_SECTION() + { + for( size_t i = 0; i < 3; ++i) { + sum[i] = mpi::allReduce( sum[i], mpi::SUM ); + min[i] = mpi::allReduce( min[i], mpi::MIN ); + max[i] = mpi::allReduce( max[i], mpi::MAX ); + } + } + + summarizedValue = sum; + minValue = min; + maxValue = max; +} + +void getMultiCompFieldQuantities( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & flagFieldID, + const BlockDataID & fieldID, + std::vector<real_t> & summarizedValue, std::vector<real_t> & minValue, std::vector<real_t> & maxValue ) +{ + std::vector<real_t> sum( 3 ); + std::vector<real_t> min( 3 ); + std::vector<real_t> max( 3 ); + for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) + { + auto field = blockIt->getData<MultiComponentField_T>( fieldID ); + auto flagField = blockIt->getData<FlagField_T>( flagFieldID ); + auto domainFlag = flagField->getFlag( Domain_Flag ); + CellInterval xyzSizeWithGhostLayers = field->xyzSizeWithGhostLayer(); + for( auto cellIt = xyzSizeWithGhostLayers.begin(); cellIt != xyzSizeWithGhostLayers.end(); ++cellIt ) + { + if( flagField->isFlagSet(*cellIt,domainFlag)) + { + real_t value0 = field->get(*cellIt, 0); + sum[0] += value0; + min[0] = std::min(min[0], value0); + max[0] = std::max(max[0], value0); + real_t value1 = field->get(*cellIt, 1); + sum[1] += value1; + min[1] = std::min(min[1], value1); + max[1] = std::max(max[1], value1); + real_t value2 = field->get(*cellIt, 2); + sum[2] += value2; + min[2] = std::min(min[2], value2); + max[2] = std::max(max[2], value2); + } + } + } + WALBERLA_MPI_SECTION() + { + mpi::allReduceInplace( sum, mpi::SUM ); + mpi::allReduceInplace( min, mpi::MIN ); + mpi::allReduceInplace( max, mpi::MAX ); + } + + summarizedValue = sum; + minValue = min; + maxValue = max; +} + +void testNearestNeighborDistributor( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & flagFieldID, + const BlockDataID & scalarFieldID, const BlockDataID & vectorFieldID, const BlockDataID & multiComponentFieldID ) +{ + // distributors + typedef field::NearestNeighborDistributor<ScalarField_T, FlagField_T> ScalarDistributor_T; + typedef field::NearestNeighborDistributor<Vec3Field_T, FlagField_T> Vec3Distributor_T; + typedef field::NearestNeighborDistributor<MultiComponentField_T, FlagField_T> MultiComponentDistributor_T; + BlockDataID scalarDistributorID = field::addDistributor< ScalarDistributor_T, FlagField_T >( blocks, scalarFieldID, flagFieldID, Domain_Flag ); + BlockDataID vectorDistributorID = field::addDistributor< Vec3Distributor_T, FlagField_T >( blocks, vectorFieldID, flagFieldID, Domain_Flag ); + BlockDataID multiComponentDistributorID = field::addDistributor< MultiComponentDistributor_T, FlagField_T >( blocks, multiComponentFieldID, flagFieldID, Domain_Flag ); + + // check scalar distribution + { + Vector3<real_t> distributionPoint(real_t(1.9), real_t(2.1), real_t(2.1)); + real_t distributionValue(real_t(100)); + auto containingBlockID = blocks->getBlock(distributionPoint); + if( containingBlockID != NULL ) + { + auto distPtr = containingBlockID->getData<ScalarDistributor_T>(scalarDistributorID); + distPtr->distribute(distributionPoint, &distributionValue); + } + real_t sum(real_t(0)), min(real_t(0)), max(real_t(0)); + getScalarFieldQuantities(blocks, flagFieldID, scalarFieldID, sum, min, max); + WALBERLA_CHECK_FLOAT_EQUAL(sum, distributionValue, "NearestNeighborDistributor: Sum of scalar distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(min, real_t(0), "NearestNeighborDistributor: Min of scalar distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(max, distributionValue, "NearestNeighborDistributor: Max of scalar distribution failed!" ); + + resetScalarField(blocks, scalarFieldID); + } + + // check vector distribution + { + Vector3<real_t> distributionPoint(real_t(5.4),real_t(2.1),real_t(3.2)); + Vector3<real_t> distributionValue(real_t(100), real_t(-10), real_t(1)); + auto containingBlockID = blocks->getBlock(distributionPoint); + if( containingBlockID != NULL ) + { + auto distPtr = containingBlockID->getData<Vec3Distributor_T>(vectorDistributorID); + distPtr->distribute(distributionPoint, &distributionValue); + } + Vector3<real_t> sum(real_t(0)), min(real_t(0)), max(real_t(0)); + getVectorFieldQuantities(blocks, flagFieldID, vectorFieldID, sum, min, max); + WALBERLA_CHECK_FLOAT_EQUAL(sum[0], distributionValue[0], "NearestNeighborDistributor: Sum of vec[0] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(sum[1], distributionValue[1], "NearestNeighborDistributor: Sum of vec[1] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(sum[2], distributionValue[2], "NearestNeighborDistributor: Sum of vec[2] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(min[0], real_t(0), "NearestNeighborDistributor: Min of vec[0] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(min[1], distributionValue[1], "NearestNeighborDistributor: Min of vec[1] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(min[2], real_t(0), "NearestNeighborDistributor: Min of vec[2] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(max[0], distributionValue[0], "NearestNeighborDistributor: Max of vec[0] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(max[1], real_t(0), "NearestNeighborDistributor: Max of vec[1] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(max[2], distributionValue[2], "NearestNeighborDistributor: Max of vec[2] distribution failed!" ); + + resetVectorField(blocks, vectorFieldID); + } + + // check multi component distribution + { + Vector3<real_t> distributionPoint(real_t(4.4),real_t(2.1),real_t(3.2)); + std::vector<real_t> distributionValue(3); + distributionValue[0] = real_t(100); + distributionValue[1] = real_t(-10); + distributionValue[2] = real_t(1); + auto containingBlockID = blocks->getBlock(distributionPoint); + if( containingBlockID != NULL ) + { + auto distPtr = containingBlockID->getData<MultiComponentDistributor_T>(multiComponentDistributorID); + distPtr->distribute(distributionPoint, distributionValue.begin()); + } + std::vector<real_t> sum(3), min(3), max(3); + getMultiCompFieldQuantities(blocks, flagFieldID, multiComponentFieldID, sum, min, max); + WALBERLA_CHECK_FLOAT_EQUAL(sum[0], distributionValue[0], "NearestNeighborDistributor: Sum of Multi Component[0] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(sum[1], distributionValue[1], "NearestNeighborDistributor: Sum of Multi Component[1] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(sum[2], distributionValue[2], "NearestNeighborDistributor: Sum of Multi Component[2] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(min[0], real_t(0), "NearestNeighborDistributor: Min of Multi Component[0] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(min[1], distributionValue[1], "NearestNeighborDistributor: Min of Multi Component[1] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(min[2], real_t(0), "NearestNeighborDistributor: Min of Multi Component[2] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(max[0], distributionValue[0], "NearestNeighborDistributor: Max of Multi Component[0] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(max[1], real_t(0), "NearestNeighborDistributor: Max of Multi Component[1] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(max[2], distributionValue[2], "NearestNeighborDistributor: Max of Multi Component[2] distribution failed!" ); + + resetMultiCompField(blocks, multiComponentFieldID); + } +} + + +void testKernelDistributor( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & flagFieldID, + const BlockDataID & scalarFieldID, const BlockDataID & vectorFieldID, const BlockDataID & multiComponentFieldID ) +{ + // distributors + typedef field::KernelDistributor<ScalarField_T, FlagField_T> ScalarDistributor_T; + typedef field::KernelDistributor<Vec3Field_T, FlagField_T> Vec3Distributor_T; + typedef field::KernelDistributor<MultiComponentField_T, FlagField_T> MultiComponentDistributor_T; + BlockDataID scalarDistributorID = field::addDistributor< ScalarDistributor_T, FlagField_T >( blocks, scalarFieldID, flagFieldID, Domain_Flag ); + BlockDataID vectorDistributorID = field::addDistributor< Vec3Distributor_T, FlagField_T >( blocks, vectorFieldID, flagFieldID, Domain_Flag ); + BlockDataID multiComponentDistributorID = field::addDistributor< MultiComponentDistributor_T, FlagField_T >( blocks, multiComponentFieldID, flagFieldID, Domain_Flag ); + + // check scalar distribution + { + Vector3<real_t> distributionPoint(real_t(1.9), real_t(2.1), real_t(2.1)); + real_t distributionValue(real_t(100)); + auto containingBlockID = blocks->getBlock(distributionPoint); + if( containingBlockID != NULL ) + { + auto distPtr = containingBlockID->getData<ScalarDistributor_T>(scalarDistributorID); + distPtr->distribute(distributionPoint, &distributionValue); + } + real_t sum(real_t(0)), min(real_t(0)), max(real_t(0)); + getScalarFieldQuantities(blocks, flagFieldID, scalarFieldID, sum, min, max); + WALBERLA_CHECK_FLOAT_EQUAL(sum, distributionValue, "KernelDistributor: Sum of scalar distribution failed!" ); + WALBERLA_CHECK(min >= real_t(0), "KernelDistributor: Min of scalar distribution failed!" ); + WALBERLA_CHECK(max <= distributionValue, "KernelDistributor: Max of scalar distribution failed!" ); + + resetScalarField(blocks, scalarFieldID); + } + + // check vector distribution + { + Vector3<real_t> distributionPoint(real_t(5.4),real_t(2.1),real_t(3.2)); + Vector3<real_t> distributionValue(real_t(100), real_t(-10), real_t(1)); + auto containingBlockID = blocks->getBlock(distributionPoint); + if( containingBlockID != NULL ) + { + auto distPtr = containingBlockID->getData<Vec3Distributor_T>(vectorDistributorID); + distPtr->distribute(distributionPoint, &distributionValue); + } + Vector3<real_t> sum(real_t(0)), min(real_t(0)), max(real_t(0)); + getVectorFieldQuantities(blocks, flagFieldID, vectorFieldID, sum, min, max); + WALBERLA_CHECK_FLOAT_EQUAL(sum[0], distributionValue[0], "KernelDistributor: Sum of vec[0] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(sum[1], distributionValue[1], "KernelDistributor: Sum of vec[1] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(sum[2], distributionValue[2], "KernelDistributor: Sum of vec[2] distribution failed!" ); + WALBERLA_CHECK(min[0] >= real_t(0), "KernelDistributor: Min of vec[0] distribution failed!" ); + WALBERLA_CHECK(min[1] >= distributionValue[1], "KernelDistributor: Min of vec[1] distribution failed!" ); + WALBERLA_CHECK(min[2] >= real_t(0), "KernelDistributor: Min of vec[2] distribution failed!" ); + WALBERLA_CHECK(max[0] <= distributionValue[0], "KernelDistributor: Max of vec[0] distribution failed!" ); + WALBERLA_CHECK(max[1] <= real_t(0), "KernelDistributor: Max of vec[1] distribution failed!" ); + WALBERLA_CHECK(max[2] <= distributionValue[2], "KernelDistributor: Max of vec[2] distribution failed!" ); + + resetVectorField(blocks, vectorFieldID); + } + + // check multi component distribution + { + Vector3<real_t> distributionPoint(real_t(4.4),real_t(2.1),real_t(3.2)); + std::vector<real_t> distributionValue(3); + distributionValue[0] = real_t(100); + distributionValue[1] = real_t(-10); + distributionValue[2] = real_t(1); + auto containingBlockID = blocks->getBlock(distributionPoint); + if( containingBlockID != NULL ) + { + auto distPtr = containingBlockID->getData<MultiComponentDistributor_T>(multiComponentDistributorID); + distPtr->distribute(distributionPoint, distributionValue.begin()); + } + std::vector<real_t> sum(3), min(3), max(3); + getMultiCompFieldQuantities(blocks, flagFieldID, multiComponentFieldID, sum, min, max); + WALBERLA_CHECK_FLOAT_EQUAL(sum[0], distributionValue[0], "KernelDistributor: Sum of Multi Component[0] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(sum[1], distributionValue[1], "KernelDistributor: Sum of Multi Component[1] distribution failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(sum[2], distributionValue[2], "KernelDistributor: Sum of Multi Component[2] distribution failed!" ); + WALBERLA_CHECK(min[0] >= real_t(0), "KernelDistributor: Min of Multi Component[0] distribution failed!" ); + WALBERLA_CHECK(min[1] >= distributionValue[1], "KernelDistributor: Min of Multi Component[1] distribution failed!" ); + WALBERLA_CHECK(min[2] >= real_t(0), "KernelDistributor: Min of Multi Component[2] distribution failed!" ); + WALBERLA_CHECK(max[0] <= distributionValue[0], "KernelDistributor: Max of Multi Component[0] distribution failed!" ); + WALBERLA_CHECK(max[1] <= real_t(0), "KernelDistributor: Max of Multi Component[1] distribution failed!" ); + WALBERLA_CHECK(max[2] <= distributionValue[2], "KernelDistributor: Max of Multi Component[2] distribution failed!" ); + + resetMultiCompField(blocks, multiComponentFieldID); + } +} + + +void testNearestNeighborDistributorAtBoundary( const shared_ptr<StructuredBlockStorage> & blocks, + const BlockDataID & flagFieldID, const BlockDataID & scalarFieldID ) +{ + // distributor + typedef field::NearestNeighborDistributor<ScalarField_T, FlagField_T> ScalarDistributor_T; + BlockDataID scalarDistributorID = field::addDistributor<ScalarDistributor_T, FlagField_T>(blocks, scalarFieldID, flagFieldID, Domain_Flag); + + // check scalar interpolation close to boundary + { + Vector3<real_t> distributionPoint(real_t(1.9), real_t(2.1), real_t(2.1)); + real_t distributionValue(real_t(100)); + auto containingBlockID = blocks->getBlock(distributionPoint); + if (containingBlockID != NULL) { + auto distPtr = containingBlockID->getData<ScalarDistributor_T>(scalarDistributorID); + distPtr->distribute(distributionPoint, &distributionValue); + } + real_t sum(real_t(0)), min(real_t(0)), max(real_t(0)); + getScalarFieldQuantities(blocks, flagFieldID, scalarFieldID, sum, min, max); + WALBERLA_CHECK_FLOAT_EQUAL(sum, distributionValue, "NearestNeighborDistributor: Sum of scalar distribution near boundary failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(min, real_t(0), "NearestNeighborDistributor: Min of scalar distribution near boundary failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(max, distributionValue, "NearestNeighborDistributor: Max of scalar distribution near boundary failed!" ); + + resetScalarField(blocks, scalarFieldID); + } + + // check scalar interpolation inside boundary + { + Vector3<real_t> distributionPoint(real_t(2.7), real_t(2.1), real_t(1.1)); + real_t distributionValue(real_t(100)); + auto containingBlockID = blocks->getBlock(distributionPoint); + if (containingBlockID != NULL) { + auto distPtr = containingBlockID->getData<ScalarDistributor_T>(scalarDistributorID); + distPtr->distribute(distributionPoint, &distributionValue); + } + real_t sum(real_t(0)), min(real_t(0)), max(real_t(0)); + getScalarFieldQuantities(blocks, flagFieldID, scalarFieldID, sum, min, max); + WALBERLA_CHECK_FLOAT_EQUAL(sum, distributionValue, "NearestNeighborDistributor: Sum of scalar distribution inside boundary failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(min, real_t(0), "NearestNeighborDistributor: Min of scalar distribution inside boundary failed!" ); + WALBERLA_CHECK_FLOAT_EQUAL(max, distributionValue, "NearestNeighborDistributor: Max of scalar distribution inside boundary failed!" ); + + resetScalarField(blocks, scalarFieldID); + } +} + +void testKernelDistributorAtBoundary( const shared_ptr<StructuredBlockStorage> & blocks, + const BlockDataID & flagFieldID, const BlockDataID & scalarFieldID ) +{ + // distributor + typedef field::KernelDistributor<ScalarField_T, FlagField_T> ScalarDistributor_T; + BlockDataID scalarDistributorID = field::addDistributor<ScalarDistributor_T, FlagField_T>(blocks, scalarFieldID, flagFieldID, Domain_Flag); + + // check scalar interpolation close to boundary + { + Vector3<real_t> distributionPoint(real_t(1.9), real_t(2.1), real_t(2.1)); + real_t distributionValue(real_t(100)); + auto containingBlockID = blocks->getBlock(distributionPoint); + if (containingBlockID != NULL) { + auto distPtr = containingBlockID->getData<ScalarDistributor_T>(scalarDistributorID); + distPtr->distribute(distributionPoint, &distributionValue); + } + real_t sum(real_t(0)), min(real_t(0)), max(real_t(0)); + getScalarFieldQuantities(blocks, flagFieldID, scalarFieldID, sum, min, max); + WALBERLA_CHECK_FLOAT_EQUAL(sum, distributionValue, "KernelDistributor: Sum of scalar distribution near boundary failed!" ); + WALBERLA_CHECK(min >= real_t(0), "KernelDistributor: Min of scalar distribution near boundary failed!" ); + WALBERLA_CHECK(max <= distributionValue, "KernelDistributor: Max of scalar distribution near boundary failed!" ); + + resetScalarField(blocks, scalarFieldID); + } + + // check scalar interpolation inside boundary + { + Vector3<real_t> distributionPoint(real_t(2.7), real_t(2.1), real_t(1.1)); + real_t distributionValue(real_t(100)); + auto containingBlockID = blocks->getBlock(distributionPoint); + if (containingBlockID != NULL) { + auto distPtr = containingBlockID->getData<ScalarDistributor_T>(scalarDistributorID); + distPtr->distribute(distributionPoint, &distributionValue); + } + real_t sum(real_t(0)), min(real_t(0)), max(real_t(0)); + getScalarFieldQuantities(blocks, flagFieldID, scalarFieldID, sum, min, max); + WALBERLA_CHECK_FLOAT_EQUAL(sum, distributionValue, "KernelDistributor: Sum of scalar distribution inside boundary failed!" ); + WALBERLA_CHECK(min >= real_t(0), "KernelDistributor: Min of scalar distribution inside boundary failed!" ); + WALBERLA_CHECK(max <= distributionValue, "KernelDistributor: Max of scalar distribution inside boundary failed!" ); + + resetScalarField(blocks, scalarFieldID); + } +} + + +int main(int argc, char **argv) { + + mpi::Environment mpiEnv(argc, argv); + debug::enterTestMode(); + + const uint_t numberOfBlocksInDirection = 2; + const uint_t numberOfCellsPerBlockInDirection = 4; + const real_t dx = real_t(1); + + // block storage + auto blocks = blockforest::createUniformBlockGrid( numberOfBlocksInDirection, numberOfBlocksInDirection, numberOfBlocksInDirection, + numberOfCellsPerBlockInDirection, numberOfCellsPerBlockInDirection, numberOfCellsPerBlockInDirection, + dx, 0, false, false, + false, false, false, + false ); + // flag field + BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>( blocks, "flag field", FieldGhostLayers, false, initFlagField ); + + // data fields + BlockDataID scalarFieldID = field::addToStorage< ScalarField_T >( blocks, "scalar field", real_t(0), field::zyxf, FieldGhostLayers ); + BlockDataID vectorFieldID = field::addToStorage< Vec3Field_T >( blocks, "vec3 field", Vector3<real_t>(real_t(0)), field::zyxf, FieldGhostLayers ); + BlockDataID multiComponentFieldID = field::addToStorage< MultiComponentField_T >( blocks, "multi component field", real_t(0), field::zyxf, FieldGhostLayers ); + + // test all distributors with domain flags everywhere, i.e. without special boundary treatment necessary + testNearestNeighborDistributor(blocks, flagFieldID, scalarFieldID, vectorFieldID, multiComponentFieldID); + testKernelDistributor(blocks, flagFieldID, scalarFieldID, vectorFieldID, multiComponentFieldID); + + // set some boundary flags in flag field and invalidate the corresponding scalar field values + setBoundaryFlags(blocks, flagFieldID ); + + // test all distributors' behavior close to boundary cells + testNearestNeighborDistributorAtBoundary(blocks, flagFieldID, scalarFieldID); + testKernelDistributorAtBoundary(blocks, flagFieldID, scalarFieldID); + + return 0; +} + +} // namespace field_distribution_tests + +int main( int argc, char **argv ){ + distribution_tests::main(argc, argv); +} \ No newline at end of file