Commit 5298876a authored by Christoph Rettinger's avatar Christoph Rettinger
Browse files

Merge branch 'new_field_utilities' into 'master'

New field utilities

See merge request walberla/walberla!61
parents 4e57552a 1e46802c
......@@ -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"
......
//======================================================================================================================
//
// 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 UniformPullReductionPackInfo.h
//! \ingroup field
//! \author Tobias Schruff <schruff@iww.rwth-aachen.de>
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#pragma once
#include "communication/UniformPackInfo.h"
#include "core/debug/Debug.h"
#include "field/GhostLayerField.h"
#include "stencil/Directions.h"
namespace walberla {
namespace field {
namespace communication {
/**
* Data packing/unpacking for ghost layer based communication of a single walberla::field::Field
*
* \ingroup comm
*
* template ReduceOperation is e.g. std::plus
*
* This pack info is used to apply a given reduce operation (e.g. +) to the values in the interior of a ghost layer field
* together with the values coming from the sender's ghost layer.
*/
template< template<typename> class ReduceOperation, typename GhostLayerField_T >
class UniformPullReductionPackInfo : public walberla::communication::UniformPackInfo
{
public:
typedef typename GhostLayerField_T::value_type T;
UniformPullReductionPackInfo( const BlockDataID & bdID ) : bdID_( bdID ), communicateAllGhostLayers_( true ),
numberOfGhostLayers_( 0 ) {}
UniformPullReductionPackInfo( const BlockDataID & bdID, const uint_t numberOfGhostLayers ) : bdID_( bdID ),
communicateAllGhostLayers_( false ), numberOfGhostLayers_( numberOfGhostLayers ) {}
virtual ~UniformPullReductionPackInfo() {}
bool constantDataExchange() const { return mpi::BufferSizeTrait<T>::constantSize; }
bool threadsafeReceiving() const { return true; }
void unpackData(IBlock * receiver, stencil::Direction dir, mpi::RecvBuffer & buffer);
void communicateLocal(const IBlock * sender, IBlock * receiver, stencil::Direction dir);
protected:
void packDataImpl(const IBlock * sender, stencil::Direction dir, mpi::SendBuffer & outBuffer) const;
uint_t numberOfGhostLayersToCommunicate( const GhostLayerField_T * const field ) const;
const BlockDataID bdID_;
bool communicateAllGhostLayers_;
uint_t numberOfGhostLayers_;
ReduceOperation<T> op_;
};
template< template<typename> class ReduceOperation, typename GhostLayerField_T >
void UniformPullReductionPackInfo<ReduceOperation, GhostLayerField_T>::unpackData( IBlock * receiver, stencil::Direction dir, mpi::RecvBuffer & buffer )
{
GhostLayerField_T * f = receiver->getData< GhostLayerField_T >( bdID_ );
WALBERLA_ASSERT_NOT_NULLPTR(f);
cell_idx_t nrOfGhostLayers = cell_idx_c( numberOfGhostLayersToCommunicate( f ) );
T buf( 0 );
for( auto i = f->beginSliceBeforeGhostLayer( dir, nrOfGhostLayers ); i != f->end(); ++i ) {
buffer >> buf;
*i = op_( *i, buf );
}
}
template< template<typename> class ReduceOperation, typename GhostLayerField_T>
void UniformPullReductionPackInfo<ReduceOperation, GhostLayerField_T>::communicateLocal( const IBlock * sender, IBlock * receiver, stencil::Direction dir )
{
const GhostLayerField_T * sf = sender ->getData< GhostLayerField_T >( bdID_ );
GhostLayerField_T * rf = receiver->getData< GhostLayerField_T >( bdID_ );
WALBERLA_ASSERT_EQUAL(sf->xSize(), rf->xSize());
WALBERLA_ASSERT_EQUAL(sf->ySize(), rf->ySize());
WALBERLA_ASSERT_EQUAL(sf->zSize(), rf->zSize());
uint_t nrOfGhostLayers = numberOfGhostLayersToCommunicate( sf );
auto srcIter = sf->beginGhostLayerOnly( nrOfGhostLayers, dir );
auto dstIter = rf->beginSliceBeforeGhostLayer( stencil::inverseDir[dir], cell_idx_c( nrOfGhostLayers ) );
while( srcIter != sf->end() ) {
*dstIter = op_( *srcIter, *dstIter );
++srcIter;
++dstIter;
}
WALBERLA_ASSERT(srcIter == sf->end() && dstIter == rf->end());
}
template< template<typename> class ReduceOperation, typename GhostLayerField_T>
void UniformPullReductionPackInfo<ReduceOperation, GhostLayerField_T>::packDataImpl( const IBlock * sender, stencil::Direction dir, mpi::SendBuffer & outBuffer ) const
{
const GhostLayerField_T * f = sender->getData< GhostLayerField_T >( bdID_ );
WALBERLA_ASSERT_NOT_NULLPTR(f);
uint_t nrOfGhostLayers = numberOfGhostLayersToCommunicate( f );
for( auto i = f->beginGhostLayerOnly( nrOfGhostLayers, dir ); i != f->end(); ++i )
outBuffer << *i;
}
template< template<typename> class ReduceOperation, typename GhostLayerField_T>
uint_t UniformPullReductionPackInfo<ReduceOperation, GhostLayerField_T>::numberOfGhostLayersToCommunicate( const GhostLayerField_T * const field ) const
{
if( communicateAllGhostLayers_ )
{
return field->nrOfGhostLayers();
}
else
{
WALBERLA_ASSERT_LESS_EQUAL( numberOfGhostLayers_, field->nrOfGhostLayers() );
return numberOfGhostLayers_;
}
}
} // namespace communication
} // namespace field
} // namespace walberla
......@@ -27,3 +27,4 @@
#include "MPIDatatypes.h"
#include "ReducePackInfo.h"
#include "UniformMPIDatatypeInfo.h"
#include "UniformPullReductionPackInfo.h"
//======================================================================================================================
//
// This file is part of waLBerla. waLBerla is free software: you can
// redistribute it and/or modify it under the terms of the GNU General Public
// License as published by the Free Software Foundation, either version 3 of
// the License, or (at your option) any later version.
//
// waLBerla is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
// for more details.
//
// You should have received a copy of the GNU General Public License along
// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
//! \file 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 weak_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 weak_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:
weak_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
//======================================================================================================================
//
// 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 weak_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() << " !");
WALBERLA_CHECK( !blockStorage_.expired() );
auto blockStorage = blockStorage_.lock();
WALBERLA_CHECK_NOT_NULLPTR(blockStorage);
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;
}
}
weak_ptr<StructuredBlockStorage> blockStorage_;
const IBlock & block_;
BaseField_T & baseField_;
const FlagField_T & flagField_;
flag_t evaluationMask_;
};
} // namespace field
} // namespace walberla
//======================================================================================================================
//
// 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 weak_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() << " !");
WALBERLA_CHECK( !blockStorage_.expired() );
auto blockStorage = blockStorage_.lock();
WALBERLA_CHECK_NOT_NULLPTR(blockStorage);
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:
weak_ptr<StructuredBlockStorage> blockStorage_;
const IBlock & block_;
BaseField_T & baseField_;
const FlagField_T & flagField_;
flag_t evaluationMask_;
};
} // namespace field
} // namespace walberla
//======================================================================================================================
//
// 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.