Commit 5921ada4 authored by Dominik Schuster's avatar Dominik Schuster
Browse files

Merge remote-tracking branch 'origin/master' into fluidizedbed_showcase

parents ad067f02 b6e20d58
......@@ -64,7 +64,7 @@ SCIterator::SCIterator(const AABB& domain, const Vector3<real_t>& pointOfReferen
auto min = domain.min() - pointOfReference_;
iReturn_ = int_c( ceil( min[0] / spacing[0] ) );
i_ = iReturn_;
jReturn_ = int_c(ceil( min[1] / spacing[1] ) + real_c(0.1));
jReturn_ = int_c(ceil( min[1] / spacing[1] ));
j_ = jReturn_;
k_ = int_c (ceil( min[2] / spacing[2] ));
......
//======================================================================================================================
//
// 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 Martin Bauer <martin.bauer@fau.de>
//! \author Florian Schornbaum <florian.schornbaum@fau.de>
//! \brief Packs only certain components of a field
//
//======================================================================================================================
#pragma once
#include "communication/UniformPackInfo.h"
#include "core/debug/Debug.h"
#include "stencil/Directions.h"
namespace walberla {
namespace field {
namespace communication {
/**
* \brief PackInfo that packs only the components that point to the neighbor.
*
* see also documentation for FieldPackInfo
*
* \warning For fields with nrOfGhostLayers > 1: only the components pointing towards
* the boundary are communicated, which may not be the desired behavior
* for the 'inner' ghost layers
*/
template< typename GhostLayerField_T, typename Stencil_T >
class StencilRestrictedPackInfo : public walberla::communication::UniformPackInfo
{
public:
StencilRestrictedPackInfo( const BlockDataID & fieldId ) : fieldId_( fieldId ) {}
virtual ~StencilRestrictedPackInfo() {}
bool constantDataExchange() const { return true; }
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;
const BlockDataID fieldId_;
static_assert(GhostLayerField_T::F_SIZE == Stencil_T::Size, "Size of stencil and f size of field have to be equal");
};
template< typename GhostLayerField_T, typename Stencil >
void StencilRestrictedPackInfo<GhostLayerField_T, Stencil>::unpackData( IBlock * receiver, stencil::Direction dir, mpi::RecvBuffer & buffer )
{
if( Stencil::idx[ stencil::inverseDir[dir] ] >= Stencil::Size )
return;
GhostLayerField_T * pdfField = receiver->getData< GhostLayerField_T >( fieldId_ );
WALBERLA_ASSERT_NOT_NULLPTR( pdfField );
WALBERLA_ASSERT_EQUAL( pdfField->nrOfGhostLayers(), 1 );
stencil::Direction packerDirection = stencil::inverseDir[dir];
for(auto i = pdfField->beginGhostLayerOnlyXYZ(dir); i != pdfField->end(); ++i )
for(uint_t f = 0; f < Stencil::d_per_d_length[packerDirection]; ++f)
buffer >> i.getF( Stencil::idx[ Stencil::d_per_d[packerDirection][f] ] );
}
template< typename GhostLayerField_T, typename Stencil >
void StencilRestrictedPackInfo<GhostLayerField_T, Stencil>::communicateLocal( const IBlock * sender, IBlock * receiver, stencil::Direction dir )
{
if( Stencil::idx[dir] >= Stencil::Size )
return;
const GhostLayerField_T * sf = sender ->getData< GhostLayerField_T >( fieldId_ );
GhostLayerField_T * rf = receiver->getData< GhostLayerField_T >( fieldId_ );
WALBERLA_ASSERT_EQUAL( sf->xyzSize(), rf->xyzSize() );
typename GhostLayerField_T::const_iterator srcIter = sf->beginSliceBeforeGhostLayerXYZ(dir);
typename GhostLayerField_T::iterator dstIter = rf->beginGhostLayerOnlyXYZ(stencil::inverseDir[dir]);
while( srcIter != sf->end() )
{
for( uint_t f = 0; f < Stencil::d_per_d_length[dir]; ++f )
dstIter.getF( Stencil::idx[ Stencil::d_per_d[dir][f] ] ) = srcIter.getF( Stencil::idx[ Stencil::d_per_d[dir][f] ] );
++srcIter;
++dstIter;
}
WALBERLA_ASSERT( srcIter == sf->end() );
WALBERLA_ASSERT( dstIter == rf->end() );
}
template< typename GhostLayerField_T, typename Stencil >
void StencilRestrictedPackInfo<GhostLayerField_T, Stencil>::packDataImpl( const IBlock * sender, stencil::Direction dir, mpi::SendBuffer & outBuffer ) const
{
if( Stencil::idx[dir] >= Stencil::Size )
return;
const GhostLayerField_T * pdfField = sender->getData< GhostLayerField_T >( fieldId_ );
WALBERLA_ASSERT_NOT_NULLPTR( pdfField );
WALBERLA_ASSERT_EQUAL( pdfField->nrOfGhostLayers(), 1 );
for( auto i = pdfField->beginSliceBeforeGhostLayerXYZ(dir); i != pdfField->end(); ++i )
for(uint_t f = 0; f < Stencil::d_per_d_length[dir]; ++f)
outBuffer << i.getF( Stencil::idx[ Stencil::d_per_d[dir][f] ] );
}
} // namespace communication
} // namespace field
} // namespace walberla
......@@ -26,12 +26,19 @@
#include "field/communication/PackInfo.h"
#include "field/communication/StencilRestrictedPackInfo.h"
#include "field/communication/UniformMPIDatatypeInfo.h"
#include "python_coupling/helper/MplHelpers.h"
#include "python_coupling/helper/BoostPythonHelpers.h"
#include "python_coupling/helper/MplHelpers.h"
#include "stencil/D2Q9.h"
#include "stencil/D3Q7.h"
#include "stencil/D3Q15.h"
#include "stencil/D3Q19.h"
#include "stencil/D3Q27.h"
namespace walberla {
namespace field {
......@@ -39,6 +46,86 @@ namespace field {
namespace internal {
//===================================================================================================================
//
// createStencilRestrictedPackInfo Export
//
//===================================================================================================================
template< typename FieldType >
typename std::enable_if<FieldType::F_SIZE == 27, boost::python::object>::type
createStencilRestrictedPackInfoObject( BlockDataID bdId )
{
typedef GhostLayerField<typename FieldType::value_type, 27> GlField_T;
using field::communication::StencilRestrictedPackInfo;
return boost::python::object( make_shared< StencilRestrictedPackInfo<GlField_T, stencil::D3Q27> >( bdId) );
}
template< typename FieldType >
typename std::enable_if<FieldType::F_SIZE == 19, boost::python::object>::type
createStencilRestrictedPackInfoObject( BlockDataID bdId )
{
typedef GhostLayerField<typename FieldType::value_type, 19> GlField_T;
using field::communication::StencilRestrictedPackInfo;
return boost::python::object( make_shared< StencilRestrictedPackInfo<GlField_T, stencil::D3Q19> >( bdId) );
}
template< typename FieldType >
typename std::enable_if<FieldType::F_SIZE == 15, boost::python::object>::type
createStencilRestrictedPackInfoObject( BlockDataID bdId )
{
typedef GhostLayerField<typename FieldType::value_type, 15> GlField_T;
using field::communication::StencilRestrictedPackInfo;
return boost::python::object( make_shared< StencilRestrictedPackInfo<GlField_T, stencil::D3Q15> >( bdId) );
}
template< typename FieldType >
typename std::enable_if<FieldType::F_SIZE == 7, boost::python::object>::type
createStencilRestrictedPackInfoObject( BlockDataID bdId )
{
typedef GhostLayerField<typename FieldType::value_type, 7> GlField_T;
using field::communication::StencilRestrictedPackInfo;
return boost::python::object( make_shared< StencilRestrictedPackInfo<GlField_T, stencil::D3Q7> >( bdId) );
}
template< typename FieldType >
typename std::enable_if<FieldType::F_SIZE == 9, boost::python::object>::type
createStencilRestrictedPackInfoObject( BlockDataID bdId )
{
typedef GhostLayerField<typename FieldType::value_type, 9> GlField_T;
using field::communication::StencilRestrictedPackInfo;
return boost::python::object( make_shared< StencilRestrictedPackInfo<GlField_T, stencil::D2Q9> >( bdId) );
}
template< typename FieldType >
typename std::enable_if<!(FieldType::F_SIZE == 9 ||
FieldType::F_SIZE == 7 ||
FieldType::F_SIZE == 15 ||
FieldType::F_SIZE == 19 ||
FieldType::F_SIZE == 27), boost::python::object>::type
createStencilRestrictedPackInfoObject( BlockDataID )
{
PyErr_SetString( PyExc_ValueError, "This works only for fields with fSize in 7, 9, 15, 19 or 27" );
throw boost::python::error_already_set();
}
FunctionExporterClass( createStencilRestrictedPackInfoObject, boost::python::object( BlockDataID ) );
template< typename FieldVector>
boost::python::object createStencilRestrictedPackInfo( const shared_ptr<StructuredBlockStorage> & bs,
const std::string & blockDataName )
{
auto bdId = python_coupling::blockDataIDFromString( *bs, blockDataName );
if ( bs->begin() == bs->end() ) {
// if no blocks are on this field an arbitrary PackInfo can be returned
return createStencilRestrictedPackInfoObject< GhostLayerField<real_t,1> > ( bdId );
}
IBlock * firstBlock = & ( * bs->begin() );
python_coupling::Dispatcher<FieldVector, Exporter_createStencilRestrictedPackInfoObject > dispatcher( firstBlock );
return dispatcher( bdId )( bdId ) ;
}
//===================================================================================================================
//
// createPackInfo Export
......@@ -114,6 +201,34 @@ namespace internal {
return dispatcher( bdId )( bdId, numberOfGhostLayers );
}
template< typename T>
void exportStencilRestrictedPackInfo()
{
using field::communication::StencilRestrictedPackInfo;
using namespace boost::python;
{
typedef StencilRestrictedPackInfo<GhostLayerField<T, 9>, stencil::D2Q9> Pi;
class_< Pi, shared_ptr<Pi>, bases<walberla::communication::UniformPackInfo>, boost::noncopyable >( "StencilRestrictedPackInfo", no_init );
}
{
typedef StencilRestrictedPackInfo<GhostLayerField<T, 7>, stencil::D3Q7> Pi;
class_< Pi, shared_ptr<Pi>, bases<walberla::communication::UniformPackInfo>, boost::noncopyable >( "StencilRestrictedPackInfo", no_init );
}
{
typedef StencilRestrictedPackInfo<GhostLayerField<T, 15>, stencil::D3Q15> Pi;
class_< Pi, shared_ptr<Pi>, bases<walberla::communication::UniformPackInfo>, boost::noncopyable >( "StencilRestrictedPackInfo", no_init );
}
{
typedef StencilRestrictedPackInfo<GhostLayerField<T, 19>, stencil::D3Q19> Pi;
class_< Pi, shared_ptr<Pi>, bases<walberla::communication::UniformPackInfo>, boost::noncopyable >( "StencilRestrictedPackInfo", no_init );
}
{
typedef StencilRestrictedPackInfo<GhostLayerField<T, 27>, stencil::D3Q27> Pi;
class_< Pi, shared_ptr<Pi>, bases<walberla::communication::UniformPackInfo>, boost::noncopyable >( "StencilRestrictedPackInfo", no_init );
}
}
} // namespace internal
......@@ -126,8 +241,13 @@ void exportCommunicationClasses()
{
using namespace boost::python;
internal::exportStencilRestrictedPackInfo<float>();
internal::exportStencilRestrictedPackInfo<double>();
def( "createMPIDatatypeInfo",&internal::createMPIDatatypeInfo<FieldTypes>, ( arg("blocks"), arg("blockDataName"), arg("numberOfGhostLayers" ) =0 ) );
def( "createPackInfo", &internal::createPackInfo<FieldTypes>, ( arg("blocks"), arg("blockDataName"), arg("numberOfGhostLayers" ) =0 ) );
def( "createStencilRestrictedPackInfo", &internal::createStencilRestrictedPackInfo<FieldTypes>,
(arg("blocks"), arg("blockDataName") ));
}
......
//======================================================================================================================
//
// 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 BinaryRawFileReader.cpp
//! \ingroup geometry
//! \author Christian Godenschwager <christian.godenschwager@fau.de>
//
//======================================================================================================================
#include "BinaryRawFile.h"
namespace walberla {
namespace geometry {
BinaryRawFile::BinaryRawFile( const std::string & filename, const Vector3< uint_t > & size, const std::string & datatype )
: size_( size )
{
init( filename, datatype );
}
BinaryRawFile::BinaryRawFile( const Config::BlockHandle & configBlock )
{
size_ = configBlock.getParameter< Vector3< uint_t > >( "size" );
const std::string filename = configBlock.getParameter< std::string >( "filename" );
const std::string datatype = configBlock.getParameter< std::string >( "datatype" );
init( filename, datatype );
}
void BinaryRawFile::init( const std::string & filename, const std::string & datatype )
{
if (datatype == "uint8")
init<uint8_t>( filename );
else if (datatype == "int8")
init<int8_t>( filename );
else if (datatype == "uint16")
init<uint16_t>( filename );
else if (datatype == "int16")
init<int16_t>( filename );
else if (datatype == "uint32")
init<uint32_t>( filename );
else if (datatype == "int32")
init<int32_t>( filename );
else if (datatype == "uint64")
init<uint64_t>( filename );
else if (datatype == "int64")
init<int64_t>( filename );
else if (datatype == "float")
init<float>( filename );
else if (datatype == "double")
init<double>( filename );
else
WALBERLA_ABORT( "Unknown datatype!" );
}
} // namespace geometry
} // namespace walberla
\ No newline at end of file
//======================================================================================================================
//
// 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 BinaryRawFileReader.h
//! \ingroup geometry
//! \author Christian Godenschwager <christian.godenschwager@fau.de>
//
//======================================================================================================================
#pragma once
#include "core/DataTypes.h"
#include "core/mpi/Broadcast.h"
#include "core/mpi/MPIManager.h"
#include "core/mpi/BufferDataTypeExtensions.h"
#include "core/math/AABB.h"
#include "core/math/Vector3.h"
#include "core/config/Config.h"
#include <fstream>
#include <string>
#include <vector>
#include <iterator>
namespace walberla {
namespace geometry {
class BinaryRawFile
{
public:
BinaryRawFile( const std::string & filename, const Vector3< uint_t > & size, const std::string & datatype );
template< typename T >
static BinaryRawFile loadFile( const std::string & filename, const Vector3< uint_t > & size );
BinaryRawFile( const Config::BlockHandle & configBlock );
inline bool get( const Vector3< uint_t > & pos ) const;
inline bool get( const uint_t x, const uint_t y, const uint_t z ) const;
const Vector3< uint_t> & size() const { return size_; }
private:
template< typename T >
BinaryRawFile( const std::string & filename, const Vector3< uint_t > & size, const T );
void init( const std::string & filename, const std::string & datatype );
template< typename T >
void init( const std::string & filename );
Vector3< uint_t > size_;
std::vector< bool > data_;
};
class BinaryRawFileInterpolator
{
public:
enum Interpolator { NEAREST_NEIGHBOR };
BinaryRawFileInterpolator( const AABB & aabb, const BinaryRawFile & binaryRawFile, Interpolator interpolator )
: aabb_(aabb), interpolator_(interpolator), binaryRawFile_( binaryRawFile ) {}
inline bool get( const Vector3< real_t > & pos ) const;
inline bool get( const real_t x, const real_t y, const real_t z ) const;
const AABB & aabb() const { return aabb_; }
Interpolator interpolator() const { return interpolator_; }
private:
inline bool getNearestNeighbor( const real_t x, const real_t y, const real_t z ) const;
AABB aabb_;
Interpolator interpolator_;
const BinaryRawFile & binaryRawFile_;
};
template<typename T>
BinaryRawFile::BinaryRawFile( const std::string & filename, const Vector3<uint_t> & size, const T )
: size_(size)
{
init<T>( filename );
}
template<typename T>
void BinaryRawFile::init( const std::string & filename )
{
const uint_t numElements = size_[0] * size_[1] * size_[2];
data_.reserve( numElements );
WALBERLA_ROOT_SECTION()
{
std::ifstream ifs( filename, std::ifstream::in | std::ifstream::binary );
std::transform( std::istream_iterator<T>( ifs ), std::istream_iterator<T>(),
std::back_inserter( data_ ), []( const T v ) { return v > T( 0 ); } );
WALBERLA_CHECK_EQUAL( numElements, data_.size(), "Error reading file \"" << filename << "\"!" );
}
mpi::broadcastObject( data_ );
}
template< typename T >
static BinaryRawFile loadFile( const std::string & filename, const Vector3< uint_t > & size )
{
return BinaryRawFileReader( filename, size, T() );
}
bool BinaryRawFile::get( const Vector3< uint_t > & pos ) const
{
return get( pos[0], pos[1], pos[2] );
}
inline bool BinaryRawFile::get( const uint_t x, const uint_t y, const uint_t z ) const
{
WALBERLA_ASSERT_LESS( x, size_[0] );
WALBERLA_ASSERT_LESS( y, size_[1] );
WALBERLA_ASSERT_LESS( z, size_[2] );
const uint_t i = z * size_[0] * size_[1] + y * size_[0] + x;
WALBERLA_ASSERT_LESS( i, data_.size() );
return data_[i];
}
bool BinaryRawFileInterpolator::get( const Vector3< real_t > & pos ) const
{
return get( pos[0], pos[1], pos[2] );
}
bool BinaryRawFileInterpolator::get( const real_t x, const real_t y, const real_t z ) const
{
WALBERLA_ASSERT( aabb_.contains( x, y, z ) );
switch (interpolator_)
{
case NEAREST_NEIGHBOR:
return getNearestNeighbor( x, y, z );
default:
WALBERLA_ABORT( "Unknown Interpolator!" );
}
}
bool BinaryRawFileInterpolator::getNearestNeighbor( const real_t x, const real_t y, const real_t z ) const
{
uint_t xInt = uint_c( (x - aabb_.xMin()) / aabb_.xSize() * real_t( binaryRawFile_.size()[0] ) + real_t( 0.5 ) );
uint_t yInt = uint_c( (y - aabb_.yMin()) / aabb_.ySize() * real_t( binaryRawFile_.size()[1] ) + real_t( 0.5 ) );
uint_t zInt = uint_c( (z - aabb_.zMin()) / aabb_.zSize() * real_t( binaryRawFile_.size()[2] ) + real_t( 0.5 ) );
return binaryRawFile_.get( xInt, yInt, zInt );
}
} // namespace geometry
} // namespace walberla
\ No newline at end of file
......@@ -24,6 +24,7 @@
#pragma once
#include "lbm/field/PdfField.h"
#include "field/communication/StencilRestrictedPackInfo.h"
#include "communication/UniformPackInfo.h"
#include "core/debug/Debug.h"
#include "stencil/Directions.h"
......@@ -55,96 +56,8 @@ namespace lbm {
* \ingroup lbm
*/
template< typename LatticeModel_T >
class PdfFieldPackInfo : public walberla::communication::UniformPackInfo
{
public:
typedef PdfField<LatticeModel_T> PdfField_T;
typedef typename LatticeModel_T::Stencil Stencil;
PdfFieldPackInfo( const BlockDataID & pdfFieldId ) : pdfFieldId_( pdfFieldId ) {}
virtual ~PdfFieldPackInfo() {}
bool constantDataExchange() const { return true; }
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;
const BlockDataID pdfFieldId_;
};
template< typename LatticeModel_T >
void PdfFieldPackInfo< LatticeModel_T >::unpackData( IBlock * receiver, stencil::Direction dir, mpi::RecvBuffer & buffer )
{
if( Stencil::idx[ stencil::inverseDir[dir] ] >= Stencil::Size )
return;
PdfField_T * pdfField = receiver->getData< PdfField_T >( pdfFieldId_ );