[BUGFIX] fixed uninitialized vertex coordinates (parmetis)

src/blockforest/loadbalancing/DynamicParMetis.cpp

@@ -124,7 +124,7 @@ bool DynamicParMetis::operator()( std::vector< std::pair< const PhantomBlock *,
    globalTimer.start();
 
    //create new communicator which excludes processes which do not have blocks
-   MPI_Comm subComm;
+   MPI_Comm subComm = MPIManager::instance()->comm();
    MPI_Group allGroup, subGroup;
    MPI_Comm_group( MPIManager::instance()->comm(), &allGroup );
    std::vector<int> ranks;
@@ -142,7 +142,7 @@ bool DynamicParMetis::operator()( std::vector< std::pair< const PhantomBlock *,
    if (subComm != MPI_COMM_NULL)
    {
       const std::pair<uint_t, uint_t> blockSequenceRange = getBlockSequenceRange( phantomForest, subComm );
-      const std::map< blockforest::BlockID, uint_t > mapping = getBlockIdToSequenceMapping( phantomForest, blockSequenceRange, subComm );
+      const std::map< blockforest::BlockID, uint_t > mapping = getBlockIdToSequenceMapping( phantomForest, blockSequenceRange, subComm ); //blockid to vertex id
 
       std::vector<int64_t> vtxdist = mpi::allGather( int64_c( blockSequenceRange.second ), subComm );
       vtxdist.insert( vtxdist.begin(), uint_t( 0 ) );

src/blockforest/loadbalancing/DynamicParMetis.h

@@ -23,10 +23,12 @@
 
 #include "blockforest/PhantomBlockForest.h"
 
+#include "core/debug/Debug.h"
 #include "core/DataTypes.h"
 #include "core/math/Vector3.h"
 #include "core/mpi/MPIWrapper.h"
 
+#include <cmath>
 #include <map>
 
 namespace walberla {
@@ -100,7 +102,7 @@ public:
    vsize_t getVertexSize() const { return vertexSize_; }
    void setVertexSize( const vsize_t size ) { vertexSize_ = size; }
 
-   const Vector3<double> & getVertexCoords() const { return vertexCoords_; }
+   const Vector3<double> & getVertexCoords() const { WALBERLA_ASSERT( !std::isnan(vertexCoords_[0]) && !std::isnan(vertexCoords_[1]) && !std::isnan(vertexCoords_[2]) ); return vertexCoords_; }
    void setVertexCoords( const Vector3<double> & p ) { vertexCoords_ = p; }
 
    void setEdgeWeight( const blockforest::BlockID & blockId, const weight_t edgeWeight ) { edgeWeights_[blockId] = edgeWeight; }
@@ -117,7 +119,7 @@ private:
 
    weight_t vertexWeight_; /// Defines the weight of a block
    vsize_t vertexSize_;    /// Defines the cost of rebalancing a block
-   Vector3<double> vertexCoords_; /// Defines where the block is located in space. Needed by some ParMetis algorithms
+   Vector3<double> vertexCoords_ = Vector3<double>(std::numeric_limits<double>::signaling_NaN()); /// Defines where the block is located in space. Needed by some ParMetis algorithms
    std::map< blockforest::BlockID, weight_t > edgeWeights_; /// Defines the cost of communication with other blocks
 };
 

tests/pe/CMakeLists.txt

@@ -69,6 +69,11 @@ waLBerla_execute_test( NAME   PE_OVERLAP )
 waLBerla_compile_test( NAME   PE_PARALLELEQUIVALENCE FILES ParallelEquivalence.cpp DEPENDS core blockforest  )
 waLBerla_execute_test( NAME   PE_PARALLELEQUIVALENCE PROCESSES 4 )
 
+if( WALBERLA_BUILD_WITH_PARMETIS )
+   waLBerla_compile_test( NAME   PE_PARMETIS FILES ParMetis.cpp DEPENDS core blockforest  )
+   waLBerla_execute_test( NAME   PE_PARMETIS PROCESSES 64 )
+endif()
+
 waLBerla_compile_test( NAME   PE_PARSEMESSAGE FILES ParseMessage.cpp DEPENDS core  )
 waLBerla_execute_test( NAME   PE_PARSEMESSAGE )
 

tests/pe/ParMetis.cpp

+//======================================================================================================================
+//
+//  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 ParMetis.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "pe/utility/CreateWorld.h"
+
+#include <blockforest/loadbalancing/DynamicParMetis.h>
+#include <core/Environment.h>
+#include <core/logging/Logging.h>
+#include <core/math/Sample.h>
+
+using namespace walberla;
+using namespace walberla::pe;
+
+class ReGrid
+{
+public:
+   void operator()( std::vector< std::pair< const Block *, uint_t > > & minTargetLevels,
+                    std::vector< const Block * > &, const BlockForest & /*forest*/ )
+   {
+      std::for_each( minTargetLevels.begin(),
+                     minTargetLevels.end(),
+                     [](auto& pair){pair.second = 1;} );
+   }
+
+};
+
+class MetisAssignmentFunctor
+{
+public:
+
+   typedef blockforest::DynamicParMetisBlockInfo           PhantomBlockWeight;
+   typedef blockforest::DynamicParMetisBlockInfoPackUnpack PhantomBlockWeightPackUnpackFunctor;
+
+   void operator()( std::vector< std::pair< const PhantomBlock *, walberla::any > > & blockData, const PhantomBlockForest & )
+   {
+      for( auto it = blockData.begin(); it != blockData.end(); ++it )
+      {
+         const auto&   corner = it->first->getAABB().maxCorner();
+         const int weight     = int_c( corner[0] + corner[1] + corner[2] );
+         blockforest::DynamicParMetisBlockInfo info(0);
+         info.setVertexWeight( weight );
+         info.setVertexSize( weight );
+         info.setVertexCoords( it->first->getAABB().center() );
+         for( uint_t nb = uint_t(0); nb < it->first->getNeighborhoodSize(); ++nb )
+         {
+            info.setEdgeWeight(it->first->getNeighborId(nb), int64_c(10) );
+         }
+         it->second = info;
+      }
+   }
+};
+
+int parmetisTest(const std::string& algorithm,
+                 const std::string& weightsToUse,
+                 const std::string& edgeSource)
+{
+   walberla::MPIManager::instance()->resetMPI();
+
+   WALBERLA_LOG_INFO_ON_ROOT("****** " << algorithm << " | " << weightsToUse << " | " << edgeSource);
+
+   // create forest
+   shared_ptr< BlockForest > forest   = createBlockForest( math::AABB(0,0,0,40,40,40),
+                                                           Vector3<uint_t>(4, 4, 4),
+                                                           Vector3<bool>(false, false, false),
+                                                           64,
+                                                           0);
+   if (!forest)
+   {
+      WALBERLA_LOG_INFO_ON_ROOT( "No BlockForest created ... exiting!");
+      return EXIT_SUCCESS;
+   }
+
+   forest->setRefreshMinTargetLevelDeterminationFunction( ReGrid() );
+
+   auto assFunc = MetisAssignmentFunctor();
+   forest->setRefreshPhantomBlockDataAssignmentFunction( assFunc );
+   forest->setRefreshPhantomBlockDataPackFunction( MetisAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor() );
+   forest->setRefreshPhantomBlockDataUnpackFunction( MetisAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor() );
+
+   auto alg     = blockforest::DynamicParMetis::stringToAlgorithm(    algorithm );
+   auto vWeight = blockforest::DynamicParMetis::stringToWeightsToUse( weightsToUse );
+   auto eWeight = blockforest::DynamicParMetis::stringToEdgeSource(   edgeSource );
+
+   auto prepFunc = blockforest::DynamicParMetis( alg, vWeight, eWeight );
+   forest->setRefreshPhantomBlockMigrationPreparationFunction( prepFunc );
+
+   forest->refresh();
+
+   math::Sample numBlocks;
+   numBlocks.castToRealAndInsert( forest->size() );
+   numBlocks.mpiGatherRoot();
+   WALBERLA_LOG_INFO_ON_ROOT("#blocks: " << numBlocks.format() );
+
+   int weight = 0;
+   for (const auto& block : *forest)
+   {
+      const auto&   corner = block.getAABB().maxCorner();
+      weight    += int_c( corner[0] + corner[1] + corner[2] );
+   }
+   math::Sample weightSample;
+   weightSample.castToRealAndInsert( weight );
+   weightSample.mpiGatherRoot();
+   WALBERLA_LOG_INFO_ON_ROOT("#weights: " << weightSample.format() );
+
+   return EXIT_SUCCESS;
+}
+
+int main( int argc, char ** argv )
+{
+   walberla::debug::enterTestMode();
+   walberla::MPIManager::instance()->initializeMPI( &argc, &argv );
+
+   std::vector<std::string> algs = {"PART_GEOM_KWAY", "PART_KWAY", "PART_ADAPTIVE_REPART", "REFINE_KWAY"};
+   std::vector<std::string> wtu  = {"NO_WEIGHTS", "EDGE_WEIGHTS", "VERTEX_WEIGHTS", "BOTH_WEIGHTS"};
+   std::vector<std::string> es   = {"EDGES_FROM_FOREST", "EDGES_FROM_EDGE_WEIGHTS"};
+
+   for (const auto& a : algs)
+      for (const auto& w : wtu)
+         for (const auto& e : es)
+         {
+            parmetisTest(a,w,e);
+         }
+
+   return EXIT_SUCCESS;
+}