Skip to content

GitLab

Commit 601c690b authored by Dominik Bartuschat's avatar Dominik Bartuschat Committed by Michael Kuron
Browse files

Galerkin coarsening for multigrid 

unified prolongate and restrict operators to work for GCA and DCA
parent fa80e3b2
Hide whitespace changes
Inline Side-by-side
src/pde/iterations/VCycles.h
View file @ 601c690b
......@@ -32,6 +32,7 @@
#include "pde/sweeps/RBGSFixedStencil.h"
#include "pde/sweeps/RBGS.h"
#include "pde/sweeps/Multigrid.h"
#include <functional>
......@@ -41,7 +42,21 @@
namespace walberla {
namespace pde {
template< typename Stencil_T >
//**********************************************************************************************************************
/*!
* \brief Class for multigrid V-cycle
*
* \tparam Stencil_T The stencil used for the discrete operator
* \tparam OperatorCoarsening_T The coarsening operator to use, defaults to direct coarsening
* \tparam Restrict_T The restriction operator to use
* \tparam ProlongateAndCorrect_T The prolongation and correction operator to use
*/
//**********************************************************************************************************************
template< typename Stencil_T,
typename OperatorCoarsening_T = CoarsenStencilFieldsDCA<Stencil_T>,
typename Restrict_T = Restrict<Stencil_T>,
typename ProlongateAndCorrect_T = ProlongateAndCorrect<Stencil_T>
>
class VCycles
{
public:
......@@ -50,6 +65,21 @@ public:
typedef std::vector< real_t > Weight_T;
typedef GhostLayerField< real_t, Stencil_T::Size > StencilField_T;
//*******************************************************************************************************************
/*! Creates a multigrid V-cycle with a fixed stencil
* \param blocks the block storage where the fields are stored
* \param uFieldId the block data id of the solution field on the finest level
* \param fFieldId the block data id of the right-hand side field on the finest level
* \param weights vector of stencil weights for the discrete operator
* \param iterations maximum number of V-cycles to perform
* \param numLvl number of grid levels to use (including the finest level)
* \param preSmoothingIters number of Gauss-Seidel iterations before restriction
* \param postSmoothingIters number of Gauss-Seidel iterations after prolongation
* \param coarseIters number of Conjugate Gradient iterations on coarsest grid
* \param residualNorm function that returns the norm of the current residuum
* \param residualNormThreshold norm threshold below which the iteration is terminated
* \param residualCheckFrequency how often to check whether the threshold has been reached
*******************************************************************************************************************/
VCycles( shared_ptr< StructuredBlockForest > blocks, const BlockDataID & uFieldId, const BlockDataID & fFieldId,
const Weight_T weights,
const uint_t iterations, const uint_t numLvl,
......@@ -59,8 +89,25 @@ public:
const Set<SUID> & requiredSelectors = Set<SUID>::emptySet(),
const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet() );
//*******************************************************************************************************************
/*! Creates a multigrid V-cycle with a stencil field
* \param blocks the block storage where the fields are stored
* \param uFieldId the block data id of the solution field on the finest level
* \param fFieldId the block data id of the right-hand side field on the finest level
* \param stencilFieldId the block data id of the stencil field for the finest level.
* The values stored in the field must not change after this class has been constructed.
* \param operatorCoarsening function that performs the stencil coarsening
* \param iterations maximum number of V-cycles to perform
* \param numLvl number of grid levels to use (including the finest level)
* \param preSmoothingIters number of Gauss-Seidel iterations before restriction
* \param postSmoothingIters number of Gauss-Seidel iterations after prolongation
* \param coarseIters number of Conjugate Gradient iterations on coarsest grid
* \param residualNorm function that returns the norm of the current residuum
* \param residualNormThreshold norm threshold below which the iteration is terminated
* \param residualCheckFrequency how often to check whether the threshold has been reached
*******************************************************************************************************************/
VCycles( shared_ptr< StructuredBlockForest > blocks, const BlockDataID & uFieldId, const BlockDataID & fFieldId,
const BlockDataID & stencilFieldId,
const BlockDataID & stencilFieldId, const OperatorCoarsening_T & operatorCoarsening,
const uint_t iterations, const uint_t numLvl,
const uint_t preSmoothingIters, const uint_t postSmoothingIters,
const uint_t coarseIters, const std::function< real_t () > & residualNorm,
......@@ -76,9 +123,7 @@ public:
const std::vector<real_t> & convergenceRate() { return convergenceRate_; }
static Vector3<uint_t> getSizeForLevel( const uint_t level, const shared_ptr< StructuredBlockStorage > & blocks, IBlock * const block );
protected:
void coarsenStencilFields();
private:
StructuredBlockForest & blocks_;
std::vector< Weight_T > weights_;
......
src/pde/iterations/VCycles.impl.h
View file @ 601c690b
......@@ -39,10 +39,10 @@
namespace walberla {
namespace pde {
template< typename Stencil_T >
VCycles< Stencil_T >::VCycles( shared_ptr< StructuredBlockForest > blocks, const BlockDataID & uFieldId, const BlockDataID & fFieldId,
const Weight_T weights,
const uint_t iterations, const uint_t numLvl,
template< typename Stencil_T, typename OperatorCoarsening_T, typename Restrict_T, typename ProlongateAndCorrect_T >
VCycles< Stencil_T, OperatorCoarsening_T, Restrict_T, ProlongateAndCorrect_T >::VCycles(
shared_ptr< StructuredBlockForest > blocks, const BlockDataID & uFieldId, const BlockDataID & fFieldId,
const Weight_T weights, const uint_t iterations, const uint_t numLvl,
const uint_t preSmoothingIters, const uint_t postSmoothingIters,
const uint_t coarseIters, const std::function< real_t () > & residualNorm,
const real_t residualNormThreshold, const uint_t residualCheckFrequency,
......@@ -54,6 +54,9 @@ VCycles< Stencil_T >::VCycles( shared_ptr< StructuredBlockForest > blocks, const
iterationsPerformed_( uint_t(0) ), thresholdReached_( false ), residualNorm_( residualNorm ), convergenceRate_(), stencilId_(),
requiredSelectors_( requiredSelectors ), incompatibleSelectors_( incompatibleSelectors )
{
static_assert(std::is_same<OperatorCoarsening_T, CoarsenStencilFieldsDCA<Stencil_T>>::value, "Use of weight requires DCA, use constructor with stencil field if you want to employ GCA");
// Set up fields for finest level
uId_.push_back( uFieldId );
fId_.push_back( fFieldId );
......@@ -87,7 +90,7 @@ VCycles< Stencil_T >::VCycles( shared_ptr< StructuredBlockForest > blocks, const
// Set up fields for coarser levels
for ( uint_t lvl = 1; lvl < numLvl; ++lvl )
{
auto getSize = std::bind(VCycles<Stencil_T>::getSizeForLevel, lvl, std::placeholders::_1, std::placeholders::_2);
auto getSize = std::bind(VCycles<Stencil_T, OperatorCoarsening_T, Restrict_T, ProlongateAndCorrect_T>::getSizeForLevel, lvl, std::placeholders::_1, std::placeholders::_2);
uId_.push_back( field::addToStorage< PdeField_T >( blocks, "u_"+boost::lexical_cast<std::string>(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
fId_.push_back( field::addToStorage< PdeField_T >( blocks, "f_"+boost::lexical_cast<std::string>(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
rId_.push_back( field::addToStorage< PdeField_T >( blocks, "r_"+boost::lexical_cast<std::string>(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
......@@ -100,7 +103,7 @@ VCycles< Stencil_T >::VCycles( shared_ptr< StructuredBlockForest > blocks, const
}
// Set up fields for CG on coarsest level
auto getFineSize = std::bind(VCycles<Stencil_T>::getSizeForLevel, numLvl-1, std::placeholders::_1, std::placeholders::_2);
auto getFineSize = std::bind(VCycles<Stencil_T, OperatorCoarsening_T, Restrict_T, ProlongateAndCorrect_T>::getSizeForLevel, numLvl-1, std::placeholders::_1, std::placeholders::_2);
dId_ = field::addToStorage< PdeField_T >( blocks, "d", getFineSize, real_t(0), field::zyxf, uint_t(1) );
zId_ = field::addToStorage< PdeField_T >( blocks, "z", getFineSize, real_t(0), field::zyxf, uint_t(1) );
......@@ -132,9 +135,9 @@ VCycles< Stencil_T >::VCycles( shared_ptr< StructuredBlockForest > blocks, const
for (uint_t lvl = 0; lvl < numLvl-1; ++lvl)
{
computeResidual_.push_back( ComputeResidualFixedStencil< Stencil_T >( blocks, uId_[lvl], fId_[lvl], weights_[lvl], rId_[lvl] ) );
restrict_.push_back( Restrict< Stencil_T >(blocks, rId_[lvl], fId_[lvl+1] ) );
restrict_.push_back( Restrict_T(blocks, rId_[lvl], fId_[lvl+1] ) );
zeroize_.push_back( Zeroize(blocks, uId_[lvl+1]) );
prolongateAndCorrect_.push_back( ProlongateAndCorrect< Stencil_T >(blocks, uId_[lvl+1], uId_[lvl]) );
prolongateAndCorrect_.push_back( ProlongateAndCorrect_T(blocks, uId_[lvl+1], uId_[lvl]) );
}
// Set up CG coarse-grid iteration
......@@ -145,9 +148,10 @@ VCycles< Stencil_T >::VCycles( shared_ptr< StructuredBlockForest > blocks, const
template< typename Stencil_T >
VCycles< Stencil_T >::VCycles( shared_ptr< StructuredBlockForest > blocks, const BlockDataID & uFieldId, const BlockDataID & fFieldId,
const BlockDataID & stencilFieldId,
template< typename Stencil_T, typename OperatorCoarsening_T, typename Restrict_T, typename ProlongateAndCorrect_T >
VCycles< Stencil_T, OperatorCoarsening_T, Restrict_T, ProlongateAndCorrect_T >::VCycles(
shared_ptr< StructuredBlockForest > blocks, const BlockDataID & uFieldId, const BlockDataID & fFieldId,
const BlockDataID & stencilFieldId, const OperatorCoarsening_T & operatorCoarsening,
const uint_t iterations, const uint_t numLvl,
const uint_t preSmoothingIters, const uint_t postSmoothingIters,
const uint_t coarseIters, const std::function< real_t () > & residualNorm,
......@@ -192,17 +196,19 @@ VCycles< Stencil_T >::VCycles( shared_ptr< StructuredBlockForest > blocks, const
// Set up fields for coarser levels
for ( uint_t lvl = 1; lvl < numLvl; ++lvl )
{
auto getSize = std::bind(VCycles<Stencil_T>::getSizeForLevel, lvl, std::placeholders::_1, std::placeholders::_2);
auto getSize = std::bind(VCycles<Stencil_T, OperatorCoarsening_T, Restrict_T, ProlongateAndCorrect_T>::getSizeForLevel, lvl, std::placeholders::_1, std::placeholders::_2);
uId_.push_back( field::addToStorage< PdeField_T >( blocks, "u_"+boost::lexical_cast<std::string>(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
fId_.push_back( field::addToStorage< PdeField_T >( blocks, "f_"+boost::lexical_cast<std::string>(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
rId_.push_back( field::addToStorage< PdeField_T >( blocks, "r_"+boost::lexical_cast<std::string>(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
stencilId_.push_back( field::addToStorage< StencilField_T >( blocks, "w_"+boost::lexical_cast<std::string>(lvl), getSize, real_t(0), field::zyxf, uint_t(1) ) );
}
coarsenStencilFields(); // scaling by ( 1/h^2 )^lvl
// CoarsenStencilFieldsDCA<Stencil_T>( blocks, stencilId_, numLvl, uint_t(2)) (); // scaling by ( 1/h^2 )^lvl
// CoarsenStencilFieldsGCA<Stencil_T>( blocks, stencilId_, numLvl, real_t(2)) ();
operatorCoarsening(stencilId_);
// Set up fields for CG on coarsest level
auto getFineSize = std::bind(VCycles<Stencil_T>::getSizeForLevel, numLvl-1, std::placeholders::_1, std::placeholders::_2);
auto getFineSize = std::bind(VCycles<Stencil_T, OperatorCoarsening_T, Restrict_T, ProlongateAndCorrect_T>::getSizeForLevel, numLvl-1, std::placeholders::_1, std::placeholders::_2);
dId_ = field::addToStorage< PdeField_T >( blocks, "d", getFineSize, real_t(0), field::zyxf, uint_t(1) );
zId_ = field::addToStorage< PdeField_T >( blocks, "z", getFineSize, real_t(0), field::zyxf, uint_t(1) );
......@@ -234,9 +240,9 @@ VCycles< Stencil_T >::VCycles( shared_ptr< StructuredBlockForest > blocks, const
for (uint_t lvl = 0; lvl < numLvl-1; ++lvl)
{
computeResidual_.push_back( ComputeResidual< Stencil_T >( blocks, uId_[lvl], fId_[lvl], stencilId_[lvl], rId_[lvl] ) );
restrict_.push_back( Restrict< Stencil_T >(blocks, rId_[lvl], fId_[lvl+1] ) );
restrict_.push_back( Restrict_T(blocks, rId_[lvl], fId_[lvl+1] ) );
zeroize_.push_back( Zeroize(blocks, uId_[lvl+1]) );
prolongateAndCorrect_.push_back( ProlongateAndCorrect< Stencil_T >(blocks, uId_[lvl+1], uId_[lvl]) );
prolongateAndCorrect_.push_back( ProlongateAndCorrect_T(blocks, uId_[lvl+1], uId_[lvl]) );
}
// Set up CG coarse-grid iteration
......@@ -247,8 +253,8 @@ VCycles< Stencil_T >::VCycles( shared_ptr< StructuredBlockForest > blocks, const
template< typename Stencil_T >
void VCycles< Stencil_T >::operator()()
template< typename Stencil_T, typename OperatorCoarsening_T, typename Restrict_T, typename ProlongateAndCorrect_T >
void VCycles< Stencil_T, OperatorCoarsening_T, Restrict_T, ProlongateAndCorrect_T >::operator()()
{
WALBERLA_LOG_PROGRESS_ON_ROOT( "Starting VCycle iteration with a maximum number of " << iterations_ << " cycles and " << numLvl_ << " levels" );
thresholdReached_ = false;
......@@ -300,32 +306,8 @@ void VCycles< Stencil_T >::operator()()
template< typename Stencil_T >
void VCycles< Stencil_T >::coarsenStencilFields()
{
const real_t scalingFactor = real_t(0.25); // scaling by ( 1/h^2 )^lvl
WALBERLA_ASSERT_EQUAL(numLvl_, stencilId_.size(), "This function can only be called when operating with stencil fields!");
for ( uint_t lvl = 1; lvl < numLvl_; ++lvl )
{
for( auto block = blocks_.begin( requiredSelectors_, incompatibleSelectors_ ); block != blocks_.end(); ++block )
{
StencilField_T * fine = block->template getData< StencilField_T >( stencilId_[lvl-1] );
StencilField_T * coarse = block->template getData< StencilField_T >( stencilId_[lvl] );
WALBERLA_FOR_ALL_CELLS_XYZ(coarse,
for( auto dir = Stencil_T::begin(); dir != Stencil_T::end(); ++dir )
coarse->get(x,y,z, dir.toIdx()) = scalingFactor * fine->get(x,y,z, dir.toIdx());
)
}
}
}
template< typename Stencil_T >
void VCycles< Stencil_T >::VCycle()
template< typename Stencil_T, typename OperatorCoarsening_T, typename Restrict_T, typename ProlongateAndCorrect_T >
void VCycles< Stencil_T, OperatorCoarsening_T, Restrict_T, ProlongateAndCorrect_T >::VCycle()
{
// pre-smoothen -- go from fine to coarse
for (uint_t l = 0; l < numLvl_-1; ++l)
......@@ -376,8 +358,8 @@ void VCycles< Stencil_T >::VCycle()
template< typename Stencil_T >
Vector3<uint_t> VCycles< Stencil_T >::getSizeForLevel( const uint_t level, const shared_ptr< StructuredBlockStorage > & blocks, IBlock * const block )
template< typename Stencil_T, typename OperatorCoarsening_T, typename Restrict_T, typename ProlongateAndCorrect_T >
Vector3<uint_t> VCycles< Stencil_T, OperatorCoarsening_T, Restrict_T, ProlongateAndCorrect_T >::getSizeForLevel( const uint_t level, const shared_ptr< StructuredBlockStorage > & blocks, IBlock * const block )
{
Vector3<uint_t> cells( blocks->getNumberOfXCells( *block ), blocks->getNumberOfYCells( *block ), blocks->getNumberOfZCells( *block ) );
......
src/pde/sweeps/Multigrid.h
View file @ 601c690b
......@@ -31,6 +31,13 @@ namespace pde {
//**********************************************************************************************************************
/*!
* \brief Restriction sweep for multigrid
*
* \tparam Stencil_T The stencil used for the discrete operator
*/
//**********************************************************************************************************************
template< typename Stencil_T >
class Restrict
{
......@@ -38,6 +45,11 @@ public:
typedef GhostLayerField< real_t, 1 > Field_T;
//*******************************************************************************************************************
/* \param blocks the block storage where the fields are stored
* \param fineFieldId the block data id of the fine field
* \param coarseFieldId the block data id of the coarse field
*******************************************************************************************************************/
Restrict( const shared_ptr< domain_decomposition::StructuredBlockStorage > & blocks,
const BlockDataID & fineFieldId, const BlockDataID & coarseFieldId ) :
blocks_( blocks ), fineFieldId_( fineFieldId ), coarseFieldId_( coarseFieldId ) {}
......@@ -53,6 +65,13 @@ private:
//**********************************************************************************************************************
/*!
* \brief Prolongation and correction sweep for multigrid
*
* \tparam Stencil_T The stencil used for the discrete operator
*/
//**********************************************************************************************************************
template< typename Stencil_T >
class ProlongateAndCorrect
{
......@@ -60,6 +79,11 @@ public:
typedef GhostLayerField< real_t, 1 > Field_T;
//*******************************************************************************************************************
/* \param blocks the block storage where the fields are stored
* \param coarseFieldId the block data id of the coarse field
* \param fineFieldId the block data id of the fine field
*******************************************************************************************************************/
ProlongateAndCorrect( const shared_ptr< domain_decomposition::StructuredBlockStorage > & blocks,
const BlockDataID & coarseFieldId, const BlockDataID & fineFieldId ) :
blocks_( blocks ), fineFieldId_( fineFieldId ), coarseFieldId_( coarseFieldId ) {}
......@@ -75,6 +99,13 @@ private:
//**********************************************************************************************************************
/*!
* \brief Residual calculation sweep for multigrid with stencil field
*
* \tparam Stencil_T The stencil used for the discrete operator
*/
//**********************************************************************************************************************
template< typename Stencil_T >
class ComputeResidual
{
......@@ -83,6 +114,13 @@ public:
typedef GhostLayerField< real_t, 1 > Field_T;
typedef GhostLayerField< real_t, Stencil_T::Size > StencilField_T;
//*******************************************************************************************************************
/* \param blocks the block storage where the fields are stored
* \param uId the block data id of the solution field
* \param fId the block data id of the right-hand side field
* \param stencilId the block data id of the stencil field
* \param rId the block data id of the residual field
*******************************************************************************************************************/
ComputeResidual( const shared_ptr< domain_decomposition::StructuredBlockStorage > & blocks, const BlockDataID & uId,
const BlockDataID & fId, const BlockDataID & stencilId, const BlockDataID & rId )
: blocks_( blocks ), uId_( uId ), fId_( fId ), stencilId_( stencilId ), rId_( rId )
......@@ -101,6 +139,13 @@ private:
//**********************************************************************************************************************
/*!
* \brief Residual calculation sweep for multigrid with fixed stencil weights
*
* \tparam Stencil_T The stencil used for the discrete operator
*/
//**********************************************************************************************************************
template< typename Stencil_T >
class ComputeResidualFixedStencil
{
......@@ -109,6 +154,13 @@ public:
typedef GhostLayerField< real_t, 1 > Field_T;
typedef std::vector< real_t > Weight_T;
//*******************************************************************************************************************
/* \param blocks the block storage where the fields are stored
* \param uId the block data id of the solution field
* \param fId the block data id of the right-hand side field
* \param weights vector of stencil weights for the discrete operator
* \param rId the block data id of the residual field
*******************************************************************************************************************/
ComputeResidualFixedStencil( const shared_ptr< domain_decomposition::StructuredBlockStorage > & blocks, const BlockDataID & uId,
const BlockDataID & fId, const std::vector <real_t> & weights, const BlockDataID & rId )
: blocks_( blocks ), uId_( uId ), fId_( fId ), weights_( weights ), rId_( rId )
......@@ -127,12 +179,21 @@ private:
//**********************************************************************************************************************
/*!
* \brief Sweep that sets all values in a field to zero
*/
//**********************************************************************************************************************
class Zeroize
{
public:
typedef GhostLayerField< real_t, 1 > Field_T;
//*******************************************************************************************************************
/* \param blocks the block storage where the field is stored
* \param fieldId the block data id of the field
*******************************************************************************************************************/
Zeroize( const shared_ptr< domain_decomposition::StructuredBlockStorage > & blocks, const BlockDataID & fieldId )
: blocks_( blocks ), fieldId_( fieldId )
{ }
......@@ -150,6 +211,97 @@ private:
//**********************************************************************************************************************
/*!
* \brief Direct Coarsening Approach for the stencil field
*
* \tparam Stencil_T The stencil used for the discrete operator
*/
//**********************************************************************************************************************
template< typename Stencil_T >
class CoarsenStencilFieldsDCA
{
public:
typedef GhostLayerField< real_t, Stencil_T::Size > StencilField_T;
//*******************************************************************************************************************
/* \param blocks the block storage where the fields are stored
* \param numLvl number of grid levels to use (including the finest level)
* \param operatorOrder the order of the (continuum) differential operator, e.g. 2 for Laplace
*******************************************************************************************************************/
CoarsenStencilFieldsDCA( shared_ptr< StructuredBlockForest > blocks,
const uint_t numLvl, const uint_t operatorOrder,
const Set<SUID> & requiredSelectors = Set<SUID>::emptySet(),
const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet() )
: blocks_( blocks ), numLvl_(numLvl), operatorOrder_(operatorOrder),
requiredSelectors_( requiredSelectors ), incompatibleSelectors_( incompatibleSelectors )
{ }
//*******************************************************************************************************************
/* \param stencilFieldId a vector of the block data ids of the stencil field for all levels (finest first)
*******************************************************************************************************************/
void operator()(const std::vector<BlockDataID> & stencilFieldId) const;
private:
shared_ptr< domain_decomposition::StructuredBlockStorage > blocks_;
uint_t numLvl_;
uint_t operatorOrder_;
Set< SUID > requiredSelectors_;
Set< SUID > incompatibleSelectors_;
};
//**********************************************************************************************************************
/*!
* \brief Galerkin Coarsening Approach for the stencil field
*
* \tparam Stencil_T The stencil used for the discrete operator
*/
//**********************************************************************************************************************
template< typename Stencil_T >
class CoarsenStencilFieldsGCA
{
public:
typedef GhostLayerField< real_t, Stencil_T::Size > StencilField_T;
//*******************************************************************************************************************
/* \param blocks the block storage where the fields are stored
* \param numLvl number of grid levels to use (including the finest level)
* \param overrelaxFact overrelaxation factor, e.g. 2 for the Poisson equation
*******************************************************************************************************************/
CoarsenStencilFieldsGCA( shared_ptr< StructuredBlockForest > blocks,
const uint_t numLvl, const real_t overrelaxFact = real_t(1),
const Set<SUID> & requiredSelectors = Set<SUID>::emptySet(),
const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet() )
: blocks_( blocks ), numLvl_(numLvl), overrelaxFact_(overrelaxFact),
requiredSelectors_( requiredSelectors ), incompatibleSelectors_( incompatibleSelectors )
{ }
//*******************************************************************************************************************
/* \param stencilFieldId a vector of the block data ids of the stencil field for all levels (finest first)
*******************************************************************************************************************/
void operator()(const std::vector<BlockDataID> & stencilFieldId) const;
private:
shared_ptr< domain_decomposition::StructuredBlockStorage > blocks_;
uint_t numLvl_;
real_t overrelaxFact_;
Set< SUID > requiredSelectors_;
Set< SUID > incompatibleSelectors_;
};
} // namespace pde
} // namespace walberla
......
src/pde/sweeps/Multigrid.impl.h
View file @ 601c690b
......@@ -24,6 +24,19 @@
#include "Multigrid.h"
#include "stencil/D3Q7.h"
#ifdef _MSC_VER
# pragma warning(push)
// disable warning for multi_array: "declaration of 'extents' hides global declaration"
# pragma warning( disable : 4459 )
#endif //_MSC_VER
#include <boost/multi_array.hpp>
#ifdef _MSC_VER
# pragma warning(pop)
#endif //_MSC_VER
namespace walberla {
......@@ -53,6 +66,7 @@ void Restrict< Stencil_T >::operator()( IBlock * const block ) const
else
{
WALBERLA_ASSERT_EQUAL(z, 0);
coarse->get(x,y,z) = real_t(0);
}
coarse->get(x,y,z) += fine->get(fx , fy , fz ) + fine->get(fx+1, fy , fz )
......@@ -140,5 +154,120 @@ void ComputeResidualFixedStencil< Stencil_T >::operator()( IBlock * const block
template< typename Stencil_T >
void CoarsenStencilFieldsDCA<Stencil_T >::operator()( const std::vector<BlockDataID> & stencilFieldId ) const
{
const real_t scalingFactor = real_t(1)/real_c(2<< (operatorOrder_-1) ); // scaling by ( 1/h^operatorOrder )^lvl
WALBERLA_ASSERT_EQUAL(numLvl_, stencilFieldId.size(), "This function can only be called when operating with stencil fields!");
for ( uint_t lvl = 1; lvl < numLvl_; ++lvl )
{
for( auto block = blocks_->begin( requiredSelectors_, incompatibleSelectors_ ); block != blocks_->end(); ++block )
{
StencilField_T * fine = block->template getData< StencilField_T >( stencilFieldId[lvl-1] );
StencilField_T * coarse = block->template getData< StencilField_T >( stencilFieldId[lvl] );
WALBERLA_FOR_ALL_CELLS_XYZ(coarse,
for( auto dir = Stencil_T::begin(); dir != Stencil_T::end(); ++dir )
coarse->get(x,y,z, dir.toIdx()) = scalingFactor * fine->get(x,y,z, dir.toIdx());
)
}
}
}
template< >
void CoarsenStencilFieldsGCA< stencil::D3Q7 >::operator()( const std::vector<BlockDataID> & stencilFieldId ) const
{
WALBERLA_ASSERT_EQUAL(numLvl_, stencilFieldId.size(), "This function can only be called when operating with stencil fields!");
// Apply Galerkin coarsening to each level //
// currently only implemented for CCMG with constant restriction and prolongation
for ( uint_t lvl = 1; lvl < numLvl_; ++lvl )
{
for( auto block = blocks_->begin( requiredSelectors_, incompatibleSelectors_ ); block != blocks_->end(); ++block )
{
StencilField_T * fine = block->getData< StencilField_T >( stencilFieldId[lvl-1] );
StencilField_T * coarse = block->getData< StencilField_T >( stencilFieldId[lvl] );
typedef boost::multi_array<real_t, 3> Array3D;
Array3D p(boost::extents[7][7][7]);
Array3D r(boost::extents[2][2][2]);
// Set to restriction weights from constant restrict operator
for(Array3D::index z=0; z<2; ++z) {
for(Array3D::index y=0; y<2; ++y) {
for(Array3D::index x=0; x<2; ++x) {
r[x][y][z] = real_t(1);
}
}
}
// Init to 0 //
for(Array3D::index k=0; k<7; ++k) {
for(Array3D::index j=0; j<7; ++j) {
for(Array3D::index i=0; i<7; ++i) {
p[i][j][k] = real_t(0);
}
}
}
// Set center to prolongation weights, including overrelaxation factor (latter therefore no longer needed in ProlongateAndCorrect)
for(Array3D::index k=0; k<2; ++k) {
for(Array3D::index j=0; j<2; ++j) {
for(Array3D::index i=0; i<2; ++i) {
p[i+2][j+2][k+2] = real_t(0.125)/overrelaxFact_; // Factor 0.125 such that same prolongation operator for DCA and GCA
}
}
}
WALBERLA_FOR_ALL_CELLS_XYZ(coarse,
Array3D ap(boost::extents[7][7][7]);
const cell_idx_t fx = 2*x;
const cell_idx_t fy = 2*y;
const cell_idx_t fz = 2*z;
for(Array3D::index k=0; k<7; ++k)
for(Array3D::index j=0; j<7; ++j)
for(Array3D::index i=0; i<7; ++i)
ap[i][j][k] = real_t(0);
// Tested for spatially varying stencils! //
for(Array3D::index k=1; k<5; ++k)
for(Array3D::index j=1; j<5; ++j)
for(Array3D::index i=1; i<5; ++i) {
ap[i][j][k] = real_t(0);
for(auto d = stencil::D3Q7::begin(); d != stencil::D3Q7::end(); ++d ){