Commit b4c848d7 authored by Jakub Klinkovský's avatar Jakub Klinkovský
Browse files

Added tests for DistributedMesh and DistributedMeshSynchronizer

parent 5a51c574
Loading
Loading
Loading
Loading
+25 −12
Original line number Diff line number Diff line
@@ -55,11 +55,32 @@ public:
   template< typename Mesh_ >
   DistributedMesh& operator=( const Mesh_& other )
   {
      getGlobalIndices< 0 >() = other.template getGlobalIndices< 0 >();
      getGlobalIndices< Mesh::getMeshDimension() >() = other.template getGlobalIndices< Mesh::getMeshDimension() >();
      localMesh = other.getLocalMesh();
      group = other.getCommunicationGroup();
      ghostLevels = other.getGhostLevels();
      vtkPointGhostTypesArray = other.vtkPointGhostTypes();
      vtkCellGhostTypesArray = other.vtkCellGhostTypes();
      return *this;
   }

   bool operator==( const DistributedMesh& other ) const
   {
      return ( getGlobalIndices< 0 >() == other.template getGlobalIndices< 0 >() &&
               getGlobalIndices< Mesh::getMeshDimension() >() == other.template getGlobalIndices< Mesh::getMeshDimension() >() &&
               localMesh == other.getLocalMesh() &&
               group == other.getCommunicationGroup() &&
               ghostLevels == other.getGhostLevels() &&
               vtkPointGhostTypesArray == other.vtkPointGhostTypes() &&
               vtkCellGhostTypesArray == other.vtkCellGhostTypes() );
   }

   bool operator!=( const DistributedMesh& other ) const
   {
      return ! operator==( other );
   }

   /**
    * Common methods redirected to the local mesh
    */
@@ -160,16 +181,6 @@ public:
      const GlobalIndexType pointsCount = localMesh.template getEntitiesCount< 0 >();
      const GlobalIndexType cellsCount = localMesh.template getEntitiesCount< Mesh::getMeshDimension() >();

      // TODO: the mesh should explicitly store ghost counts (or offsets) - useful for efficient iteration
      GlobalIndexType ghostPoints = 0;
      for( GlobalIndexType i = 0; i < pointsCount; i++ )
         if( localMesh.template isGhostEntity< 0 >( i ) )
            ghostPoints++;
      GlobalIndexType ghostCells = 0;
      for( GlobalIndexType i = 0; i < cellsCount; i++ )
         if( localMesh.template isGhostEntity< Mesh::getMeshDimension() >( i ) )
            ghostCells++;

      CommunicatorType::Barrier();
      for( int i = 0; i < CommunicatorType::GetSize(); i++ ) {
         if( i == CommunicatorType::GetRank() ) {
@@ -179,8 +190,10 @@ public:
                << "\tCells count:\t" << cellsCount << "\n"
                << "\tPoints count:\t" << pointsCount << "\n"
                << "\tGhost levels:\t" << getGhostLevels() << "\n"
                << "\tGhost cells:\t" << ghostCells << "\n"
                << "\tGhost points:\t" << ghostPoints << "\n";
                << "\tGhost cells count:\t" << localMesh.template getGhostEntitiesCount< Mesh::getMeshDimension() >() << "\n"
                << "\tGhost points count:\t" << localMesh.template getGhostEntitiesCount< 0 >() << "\n"
                << "\tBoundary cells count:\t" << localMesh.template getBoundaryEntitiesCount< Mesh::getMeshDimension() >() << "\n"
                << "\tBoundary points count:\t" << localMesh.template getBoundaryEntitiesCount< 0 >() << "\n";
            const GlobalIndexType globalPointIndices = getGlobalIndices< 0 >().getSize();
            const GlobalIndexType globalCellIndices = getGlobalIndices< Mesh::getMeshDimension() >().getSize();
            if( getGhostLevels() > 0 ) {
+41 −32
Original line number Diff line number Diff line
@@ -38,10 +38,15 @@ public:

   DistributedMeshSynchronizer() = default;

   void initialize( const Mesh& mesh )
   template< typename DistributedMesh >
   void initialize( const DistributedMesh& mesh )
   {
      static_assert( std::is_same< typename Mesh::CommunicatorType, CommunicatorType >::value,
                     "Mesh::CommunnicatorType does not match" );
      static_assert( std::is_same< typename DistributedMesh::MeshType, typename MeshFunction::MeshType >::value,
                     "The type of the local mesh does not match the type of the mesh function's mesh." );
      static_assert( std::is_same< typename DistributedMesh::CommunicatorType, CommunicatorType >::value,
                     "DistributedMesh::CommunnicatorType does not match" );
      if( mesh.getGhostLevels() <= 0 )
         throw std::logic_error( "There are no ghost levels on the distributed mesh." );
      localEntitiesCount = mesh.getLocalMesh().template getEntitiesCount< EntityDimension >();

      // GOTCHA: https://devblogs.nvidia.com/cuda-pro-tip-always-set-current-device-avoid-multithreading-bugs/
@@ -56,10 +61,10 @@ public:

      // exchange the global index offsets so that each rank can determine the
      // owner of every entity by its global index
      const IndexType ownStart = mesh.template getGlobalIndices< EntityDimension >()[ 0 ];
      Containers::Array< IndexType, DeviceType > offsets( nproc );
      const IndexType ownStart = mesh.template getGlobalIndices< EntityDimension >().getElement( 0 );
      Containers::Array< IndexType, Devices::Host, int > offsets( nproc );
      {
         Containers::Array< IndexType, DeviceType > sendbuf( nproc );
         Containers::Array< IndexType, Devices::Host, int > sendbuf( nproc );
         sendbuf.setValue( ownStart );
         CommunicatorType::Alltoall( sendbuf.getData(), 1,
                                     offsets.getData(), 1,
@@ -75,31 +80,27 @@ public:
         return nproc - 1;
      };

      // TODO: initialization of the distributed mesh should set the ranges for
      //       local and ghost entities so we can iterate just over the ghost
      //       entities
      // FIXME: the local vertices are not contiguous !!!
      // count local ghost entities for each rank
      Containers::Array< IndexType, DeviceType > localGhostCounts( nproc );
      Containers::Array< IndexType, Devices::Host, int > localGhostCounts( nproc );
      localGhostCounts.setValue( 0 );
      IndexType firstGhost = 0;
      bool ghost_found = false;
      for( IndexType local_idx = 0; local_idx < localEntitiesCount; local_idx++ ) {
         const IndexType global_idx = mesh.template getGlobalIndices< EntityDimension >()[ local_idx ];
         const bool isNeighbor = mesh.getLocalMesh().template isGhostEntity< EntityDimension >( local_idx );
         if( ! isNeighbor ) {
            // TODO: this is just for testing/debugging
            if( ghost_found )
               std::cerr << "ERROR: global indices of local entities are not contiguous (local index " << local_idx << " is after a ghost entity has been found)" << std::endl;
            ++firstGhost;
         }
         else {
      Containers::Array< IndexType, Devices::Host, IndexType > hostGlobalIndices;
      hostGlobalIndices = mesh.template getGlobalIndices< EntityDimension >();
      IndexType prev_global_idx = 0;
      for( IndexType local_idx = mesh.getLocalMesh().template getGhostEntitiesOffset< EntityDimension >();
           local_idx < localEntitiesCount;
           local_idx++ )
      {
         if( ! std::is_same< DeviceType, Devices::Cuda >::value )
         if( ! mesh.getLocalMesh().template isGhostEntity< EntityDimension >( local_idx ) )
            throw std::runtime_error( "encountered local entity while iterating over ghost entities - the mesh is probably inconsistent or there is a bug in the DistributedMeshSynchronizer" );
         const IndexType global_idx = hostGlobalIndices[ local_idx ];
         if( global_idx < prev_global_idx )
            throw std::runtime_error( "ghost indices are not sorted - the mesh is probably inconsistent or there is a bug in the DistributedMeshSynchronizer" );
         prev_global_idx = global_idx;
         const int owner = getOwner( global_idx );
            if( owner != rank ) {
         if( owner == rank )
            throw std::runtime_error( "the owner of a ghost entity cannot be the local rank - the mesh is probably inconsistent or there is a bug in the DistributedMeshSynchronizer" );
         ++localGhostCounts[ owner ];
               ghost_found = true;
            }
         }
      }

      // exchange the ghost counts
@@ -135,6 +136,7 @@ public:
         std::vector< typename CommunicatorType::Request > requests;

         // send our ghost indices to the neighboring ranks
         IndexType firstGhost = mesh.getLocalMesh().template getGhostEntitiesOffset< EntityDimension >();
         for( int i = 0; i < nproc; i++ ) {
            if( ghostEntitiesCounts( rank, i ) > 0 ) {
               requests.push_back( CommunicatorType::ISend(
@@ -211,12 +213,19 @@ public:
         }
      }

      auto sendBuffersView = sendBuffers.getView();
      const auto ghostNeighborsView = ghostNeighbors.getConstView();
      const auto functionDataView = function.getData().getConstView();
      auto copy_kernel = [sendBuffersView, functionDataView, ghostNeighborsView] __cuda_callable__ ( IndexType k, IndexType offset ) mutable
      {
         sendBuffersView[ offset + k ] = functionDataView[ ghostNeighborsView[ offset + k ] ];
      };

      for( int i = 0; i < nproc; i++ ) {
         if( ghostEntitiesCounts( i, rank ) > 0 ) {
            const IndexType offset = ghostNeighborOffsets[ i ];
            // copy data to send buffers
            for( IndexType k = 0; k < ghostEntitiesCounts( i, rank ); k++ )
               sendBuffers[ offset + k ] = function.getData()[ ghostNeighbors[ offset + k ] ];
            Algorithms::ParallelFor< DeviceType >::exec( (IndexType) 0, ghostEntitiesCounts( i, rank ), copy_kernel, offset );

            // issue async send operation
            requests.push_back( CommunicatorType::ISend(
@@ -258,7 +267,7 @@ protected:
    * entity owned by the i-th rank. All ghost entities owned by the i-th
    * rank are assumed to be indexed contiguously in the local mesh.
    */
   Containers::Array< IndexType, DeviceType, IndexType > ghostOffsets;
   Containers::Array< IndexType, Devices::Host, int > ghostOffsets;

   /**
    * Ghost neighbor offsets: array of size nproc + 1 where the i-th value is
+3 −3
Original line number Diff line number Diff line
@@ -103,11 +103,11 @@ protected:
   };

   template< int Dimension >
   struct UpdateBoundaryIndices
   struct UpdateEntityTagsLayer
   {
      static void exec( Mesh& mesh )
      {
         mesh.template updateBoundaryIndices< Dimension >();
         mesh.template updateEntityTagsLayer< Dimension >();
      }
   };

@@ -148,7 +148,7 @@ public:
                                                      kernel,
                                                      &meshPointer.template modifyData< DeviceType >() );

         Algorithms::TemplateStaticFor< int, 0, Mesh::getMeshDimension() + 1, UpdateBoundaryIndices >::execHost( mesh );
         Algorithms::TemplateStaticFor< int, 0, Mesh::getMeshDimension() + 1, UpdateEntityTagsLayer >::execHost( mesh );
      }
   };

+44 −33
Original line number Diff line number Diff line
@@ -57,6 +57,7 @@ public:
      tags = other.tags;
      boundaryIndices = other.boundaryIndices;
      interiorIndices = other.interiorIndices;
      ghostsOffset = other.ghostsOffset;
      return *this;
   }

@@ -69,6 +70,7 @@ public:
      tags = other.tags;
      boundaryIndices = other.boundaryIndices;
      interiorIndices = other.interiorIndices;
      ghostsOffset = other.ghostsOffset;
      return *this;
   }

@@ -76,6 +78,7 @@ public:
   void setEntitiesCount( DimensionTag, const GlobalIndexType& entitiesCount )
   {
      tags.setSize( entitiesCount );
      ghostsOffset = entitiesCount;
   }

   void resetEntityTags( DimensionTag )
@@ -113,41 +116,36 @@ public:
      return tags[ entityIndex ] & EntityTags::GhostEntity;
   }

   void updateBoundaryIndices( DimensionTag )
   void updateEntityTagsLayer( DimensionTag )
   {
      // count boundary entities - custom reduction because expression templates don't support filtering bits this way:
      //    const GlobalIndexType boundaryEntities = sum(cast< GlobalIndexType >( _tagsVector & EntityTags::BoundaryEntity ));
      // NOTE: boundary and ghost entities may overlap for vertices, so we count all categories separately
      // NOTE: boundary/interior entities may overlap with ghost entities, so we count all categories separately
      const auto tags_view = tags.getConstView();
      auto is_interior = [=] __cuda_callable__ ( GlobalIndexType entityIndex ) -> GlobalIndexType
      {
         return ! (tags_view[ entityIndex ] & ( EntityTags::BoundaryEntity | EntityTags::GhostEntity ));
      };
      auto is_boundary = [=] __cuda_callable__ ( GlobalIndexType entityIndex ) -> GlobalIndexType
      {
         return tags_view[ entityIndex ] & EntityTags::BoundaryEntity;
         return bool(tags_view[ entityIndex ] & EntityTags::BoundaryEntity);
      };
      auto is_ghost = [=] __cuda_callable__ ( GlobalIndexType entityIndex ) -> GlobalIndexType
      {
         return tags_view[ entityIndex ] & EntityTags::GhostEntity;
         return bool(tags_view[ entityIndex ] & EntityTags::GhostEntity);
      };
      const GlobalIndexType interiorEntities = Algorithms::Reduction< Device >::reduce( (GlobalIndexType) 0, tags.getSize(), std::plus<>{}, is_interior, (GlobalIndexType) 0 );
      const GlobalIndexType boundaryEntities = Algorithms::Reduction< Device >::reduce( (GlobalIndexType) 0, tags.getSize(), std::plus<>{}, is_boundary, (GlobalIndexType) 0 );
      const GlobalIndexType ghostEntities = Algorithms::Reduction< Device >::reduce( (GlobalIndexType) 0, tags.getSize(), std::plus<>{}, is_ghost, (GlobalIndexType) 0 );

      interiorIndices.setSize( interiorEntities );
      interiorIndices.setSize( tags.getSize() - boundaryEntities );
      boundaryIndices.setSize( boundaryEntities );
      ghostIndices.setSize( ghostEntities );
      ghostsOffset = tags.getSize() - ghostEntities;

      if( ! std::is_same< Device, Devices::Cuda >::value ) {
         GlobalIndexType i = 0, b = 0, g = 0;
         GlobalIndexType i = 0, b = 0;
         for( GlobalIndexType e = 0; e < tags.getSize(); e++ ) {
            if( ! (tags[ e ] & ( EntityTags::BoundaryEntity | EntityTags::GhostEntity )) )
               interiorIndices[ i++ ] = e;
            if( tags[ e ] & EntityTags::BoundaryEntity )
               boundaryIndices[ b++ ] = e;
            if( tags[ e ] & EntityTags::GhostEntity )
               ghostIndices[ g++ ] = e;
            else
               interiorIndices[ i++ ] = e;
            if( tags[ e ] & EntityTags::GhostEntity && ghostEntities > 0 && e < ghostsOffset )
               throw std::runtime_error( "The mesh is inconsistent - ghost entities of dimension " + std::to_string(DimensionTag::value) + " are not ordered after local entities." );
         }
      }
      // TODO: parallelize directly on the device
@@ -156,28 +154,26 @@ public:
         using OrderingHostArray   = typename OrderingArray::template Self< typename OrderingArray::ValueType, Devices::Host >;

         EntityTagsHostArray hostTags;
         OrderingHostArray hostBoundaryIndices, hostInteriorIndices, hostGhostIndices;
         OrderingHostArray hostBoundaryIndices, hostInteriorIndices;

         hostTags.setLike( tags );
         hostInteriorIndices.setLike( interiorIndices );
         hostBoundaryIndices.setLike( boundaryIndices );
         hostGhostIndices.setLike( ghostIndices );

         hostTags = tags;

         GlobalIndexType i = 0, b = 0, g = 0;
         GlobalIndexType i = 0, b = 0;
         for( GlobalIndexType e = 0; e < tags.getSize(); e++ ) {
            if( ! (hostTags[ e ] & ( EntityTags::BoundaryEntity | EntityTags::GhostEntity )) )
               hostInteriorIndices[ i++ ] = e;
            if( hostTags[ e ] & EntityTags::BoundaryEntity )
               hostBoundaryIndices[ b++ ] = e;
            if( hostTags[ e ] & EntityTags::GhostEntity )
               hostGhostIndices[ g++ ] = e;
            else
               hostInteriorIndices[ i++ ] = e;
            if( hostTags[ e ] & EntityTags::GhostEntity && ghostEntities > 0 && e < ghostsOffset )
               throw std::runtime_error( "The mesh is inconsistent - ghost entities of dimension " + std::to_string(DimensionTag::value) + " are not ordered after local entities." );
         }

         interiorIndices = hostInteriorIndices;
         boundaryIndices = hostBoundaryIndices;
         ghostIndices = hostGhostIndices;
      }
   }

@@ -205,6 +201,18 @@ public:
      return interiorIndices[ i ];
   }

   __cuda_callable__
   GlobalIndexType getGhostEntitiesCount( DimensionTag ) const
   {
      return tags.getSize() - ghostsOffset;
   }

   __cuda_callable__
   GlobalIndexType getGhostEntitiesOffset( DimensionTag ) const
   {
      return ghostsOffset;
   }

   void save( File& file ) const
   {
      file << tags;
@@ -213,7 +221,7 @@ public:
   void load( File& file )
   {
      file >> tags;
      updateBoundaryIndices( DimensionTag() );
      updateEntityTagsLayer( DimensionTag() );
   }

   void print( std::ostream& str ) const
@@ -224,14 +232,14 @@ public:
      str << interiorIndices << std::endl;
      str << "Indices of the boundary entities of dimension " << DimensionTag::value << " are: ";
      str << boundaryIndices << std::endl;
      str << "Indices of the ghost entities of dimension " << DimensionTag::value << " are: ";
      str << ghostIndices << std::endl;
      str << "Index of the first ghost entity of dimension " << DimensionTag::value << " is: ";
      str << ghostsOffset << std::endl;
   }

   bool operator==( const Layer& layer ) const
   {
      TNL_ASSERT( ( tags == layer.tags && interiorIndices == layer.interiorIndices && boundaryIndices == layer.boundaryIndices && ghostIndices == layer.ghostIndices ) ||
                  ( tags != layer.tags && interiorIndices != layer.interiorIndices && boundaryIndices != layer.boundaryIndices && ghostIndices != layer.ghostIndices ),
      TNL_ASSERT( ( tags == layer.tags && interiorIndices == layer.interiorIndices && boundaryIndices == layer.boundaryIndices && ghostsOffset == layer.ghostsOffset ) ||
                  ( tags != layer.tags && interiorIndices != layer.interiorIndices && boundaryIndices != layer.boundaryIndices && ghostsOffset != layer.ghostsOffset ),
                  std::cerr << "The EntityTags layer is in inconsistent state - this is probably a bug in the boundary tags initializer." << std::endl
                            << "tags                  = " << tags << std::endl
                            << "layer.tags            = " << layer.tags << std::endl
@@ -239,14 +247,15 @@ public:
                            << "layer.interiorIndices = " << layer.interiorIndices << std::endl
                            << "boundaryIndices       = " << boundaryIndices << std::endl
                            << "layer.boundaryIndices = " << layer.boundaryIndices << std::endl
                            << "ghostIndices          = " << ghostIndices << std::endl
                            << "layer.ghostIndices    = " << layer.ghostIndices << std::endl; );
                            << "ghostsOffset          = " << ghostsOffset << std::endl
                            << "layer.ghostsOffset    = " << layer.ghostsOffset << std::endl; );
      return tags == layer.tags;
   }

private:
   EntityTagsArrayType tags;
   OrderingArray interiorIndices, boundaryIndices, ghostIndices;
   OrderingArray interiorIndices, boundaryIndices;
   GlobalIndexType ghostsOffset = 0;

   // friend class is needed for templated assignment operators
   template< typename MeshConfig_, typename Device_, typename DimensionTag_, bool TagStorage_ >
@@ -277,11 +286,13 @@ protected:
   void removeEntityTag( DimensionTag, const GlobalIndexType&, TagType ) const {}
   void isBoundaryEntity( DimensionTag, const GlobalIndexType& ) const {}
   void isGhostEntity( DimensionTag, const GlobalIndexType& ) const {}
   void updateBoundaryIndices( DimensionTag ) {}
   void updateEntityTagsLayer( DimensionTag ) {}
   void getBoundaryEntitiesCount( DimensionTag ) const {}
   void getBoundaryEntityIndex( DimensionTag, const GlobalIndexType& i ) const {}
   void getInteriorEntitiesCount( DimensionTag ) const {}
   void getInteriorEntityIndex( DimensionTag, const GlobalIndexType& i ) const {}
   void getGhostEntitiesCount() const;
   void getGhostEntitiesOffset() const;

   void save( File& file ) const {}
   void load( File& file ) {}
+38 −16
Original line number Diff line number Diff line
@@ -34,11 +34,13 @@ protected:
   using LayerType::removeEntityTag;
   using LayerType::isBoundaryEntity;
   using LayerType::isGhostEntity;
   using LayerType::updateBoundaryIndices;
   using LayerType::updateEntityTagsLayer;
   using LayerType::getBoundaryEntitiesCount;
   using LayerType::getBoundaryEntityIndex;
   using LayerType::getInteriorEntitiesCount;
   using LayerType::getInteriorEntityIndex;
   using LayerType::getGhostEntitiesCount;
   using LayerType::getGhostEntitiesOffset;

   using BaseType::setEntitiesCount;
   using BaseType::resetEntityTags;
@@ -47,11 +49,13 @@ protected:
   using BaseType::removeEntityTag;
   using BaseType::isBoundaryEntity;
   using BaseType::isGhostEntity;
   using BaseType::updateBoundaryIndices;
   using BaseType::updateEntityTagsLayer;
   using BaseType::getBoundaryEntitiesCount;
   using BaseType::getBoundaryEntityIndex;
   using BaseType::getInteriorEntitiesCount;
   using BaseType::getInteriorEntityIndex;
   using BaseType::getGhostEntitiesCount;
   using BaseType::getGhostEntitiesOffset;


   LayerInheritor() = default;
@@ -114,16 +118,18 @@ class LayerInheritor< MeshConfig, Device, DimensionTag< MeshConfig::meshDimensio
protected:
   void setEntitiesCount();
   void resetEntityTags();
   void getEntityTag();
   void getEntityTag() const;
   void addEntityTag();
   void removeEntityTag();
   void isBoundaryEntity();
   void isGhostEntity();
   void updateBoundaryIndices();
   void getBoundaryEntitiesCount();
   void getBoundaryEntityIndex();
   void getInteriorEntitiesCount();
   void getInteriorEntityIndex();
   void isBoundaryEntity() const;
   void isGhostEntity() const;
   void updateEntityTagsLayer();
   void getBoundaryEntitiesCount() const;
   void getBoundaryEntityIndex() const;
   void getInteriorEntitiesCount() const;
   void getInteriorEntityIndex() const;
   void getGhostEntitiesCount() const;
   void getGhostEntitiesOffset() const;

   LayerInheritor() = default;
   explicit LayerInheritor( const LayerInheritor& other ) {}
@@ -240,6 +246,28 @@ public:
      return BaseType::getInteriorEntityIndex( DimensionTag< Dimension >(), i );
   }

   template< int Dimension >
   __cuda_callable__
   GlobalIndexType getGhostEntitiesCount() const
   {
      static_assert( WeakTrait< Dimension >::entityTagsEnabled, "You try to access entity tags which are not configured for storage." );
      return BaseType::getGhostEntitiesCount( DimensionTag< Dimension >() );
   }

   template< int Dimension >
   __cuda_callable__
   GlobalIndexType getGhostEntitiesOffset() const
   {
      static_assert( WeakTrait< Dimension >::entityTagsEnabled, "You try to access entity tags which are not configured for storage." );
      return BaseType::getGhostEntitiesOffset( DimensionTag< Dimension >() );
   }

   template< int Dimension >
   void updateEntityTagsLayer()
   {
      BaseType::updateEntityTagsLayer( DimensionTag< Dimension >() );
   }

protected:
   template< int Dimension >
   void entityTagsSetEntitiesCount( const GlobalIndexType& entitiesCount )
@@ -252,12 +280,6 @@ protected:
   {
      BaseType::resetEntityTags( DimensionTag< Dimension >() );
   }

   template< int Dimension >
   void updateBoundaryIndices()
   {
      BaseType::updateBoundaryIndices( DimensionTag< Dimension >() );
   }
};

} // namespace EntityTags
Loading