Commit 6a13ce6f authored by Jakub Klinkovský's avatar Jakub Klinkovský
Browse files

Refactored functions in tnl-decompose-mesh 

Despite the ugly diff, the code did not change much - only the two
static class member functions became plain functions, changed order and
indentation level, and the struct DecomposeMesh was removed.
parent e183940d

src/Tools/tnl-decompose-mesh.cpp

+446 −448
Original line number Diff line number Diff line
@@ -277,156 +277,7 @@ void setMETISoptions( idx_t options[METIS_NOPTIONS], const Config::ParameterCont 
}



template< typename Mesh >

struct DecomposeMesh

{

   static void run( const Mesh& mesh, const Config::ParameterContainer& parameters )

   {

      using Index = typename Mesh::GlobalIndexType;



      // warn if input mesh has 64-bit indices, but METIS uses only 32-bit indices

      if( IDXTYPEWIDTH == 32 && sizeof(Index) > 4 )

         std::cerr << "Warning: the input mesh uses 64-bit indices, but METIS was compiled only with 32-bit indices. "

                      "Decomposition may not work correctly if the index values overflow the 32-bit type." << std::endl;



      // get the mesh connectivity information in a format suitable for METIS. Actually, the same

      // format is used by the XML-based VTK formats - the only difference is that METIS requires

      // `offsets` to start with 0.

      std::vector< idx_t > connectivity, offsets;

      offsets.push_back(0);

      const Index cellsCount = mesh.template getEntitiesCount< typename Mesh::Cell >();

      for( Index i = 0; i < cellsCount; i++ ) {

         const auto& entity = mesh.template getEntity< typename Mesh::Cell >( i );

         const Index subvertices = entity.template getSubentitiesCount< 0 >();

         for( Index j = 0; j < subvertices; j++ )

            connectivity.push_back( entity.template getSubentityIndex< 0 >( j ) );

         offsets.push_back( connectivity.size() );

      }



      // number of elements (cells)

      idx_t ne = mesh.template getEntitiesCount< typename Mesh::Cell >();

      // number of nodes (vertices)

      idx_t nn = mesh.template getEntitiesCount< 0 >();

      // pointers to arrays storing the mesh in a CSR-like format

      idx_t* eptr = offsets.data();

      idx_t* eind = connectivity.data();

      // Specifies the number of common nodes that two elements must have in order to put an edge

      // between them in the dual graph.

      idx_t ncommon = Mesh::getMeshDimension();

      if( parameters.checkParameter( "min-common-vertices" ) )

         ncommon = parameters.template getParameter< unsigned >( "min-common-vertices" );

      // numbering scheme for the adjacency structure of a graph or element-node structure of a mesh

      // (0 is C-style, 1 is Fortran-style)

      idx_t numflag = 0;

      // These arrays store the adjacency structure of the generated dual graph. The format is of

      // the adjacency structure is described in Section 5.5 of the METIS manual. Memory for these

      // arrays is allocated by METIS’ API using the standard malloc function. It is the

      // responsibility of the application to free this memory by calling free. METIS provides the

      // METIS_Free that is a wrapper to C’s free function.

      idx_t* xadj = nullptr;

      idx_t* adjncy = nullptr;



      // We could use METIS_PartMeshDual directly instead of METIS_MeshToDual + METIS_PartGraph*,

      // but we need to reuse the dual graph for the generation of ghost cells.

      std::cout << "Running METIS_MeshToDual..." << std::endl;

      int status = METIS_MeshToDual(&ne, &nn, eptr, eind, &ncommon, &numflag, &xadj, &adjncy);



      // wrap xadj and adjncy with shared_ptr

      std::shared_ptr<idx_t> shared_xadj {xadj, METIS_Free};

      std::shared_ptr<idx_t> shared_adjncy {adjncy, METIS_Free};



      switch( status )

      {

         case METIS_OK: break;

         case METIS_ERROR_INPUT:

            throw std::runtime_error( "METIS_MeshToDual failed due to an input error." );

         case METIS_ERROR_MEMORY:

            throw std::runtime_error( "METIS_MeshToDual failed due to a memory allocation error." );

         case METIS_ERROR:

         default:

            throw std::runtime_error( "METIS_MeshToDual failed with an unspecified error." );

      }



      // The number of vertices in the graph.

      idx_t nvtxs = ne;

      // The number of balancing constraints. It should be at least 1.

      idx_t ncon = 1;

      // An array of size `ne` specifying the weights of the elements. A NULL value can be passed

      // to indicate that all elements have an equal weight.

      idx_t* vwgt = nullptr;

      // An array of size `ne` specifying the size of the elements that is used for computing the

      // total communication volume as described in Section 5.7 of the METIS manual. A NULL value

      // can be passed when the objective is cut or when all elements have an equal size.

      idx_t* vsize = nullptr;

      // The weights of the edges as described in Section 5.5 of the METIS manual.

      idx_t* adjwgt = nullptr;  // METIS_PartMeshDual uses NULL too

      // The number of parts to partition the mesh.

      idx_t nparts = parameters.template getParameter< unsigned >( "subdomains" );

      // This is an array of size nparts that specifies the desired weight for each partition. The

      // target partition weight for the i-th partition is specified at tpwgts[i] (the numbering for

      // the partitions starts from 0). The sum of the tpwgts[] entries must be 1.0. A NULL value

      // can be passed to indicate that the graph should be equally divided among the partitions.

      real_t* tpwgts = nullptr;

      // This is an array of size ncon that specifies the allowed load imbalance tolerance for each

      // constraint. For the i-th partition and j-th constraint the allowed weight is the

      // ubvec[j]*tpwgts[i*ncon+j] fraction of the j-th’s constraint total weight. The load

      // imbalances must be greater than 1.0. A NULL value can be passed indicating that the load

      // imbalance tolerance for each constraint should be 1.001 (for ncon=1) or 1.01 (for ncon>1).

      real_t* ubvec = nullptr;  // METIS_PartMeshDual uses NULL too

      // Upon successful completion, this variable stores either the edgecut or the total

      // communication volume of the dual graph’s partitioning.

      idx_t objval = 0;

      // Array of size `ne` that upon successful completion stores the partition array for the

      // elements of the mesh.

      MetisIndexArray part_array( nvtxs );

      idx_t* part = part_array.getData();



      // Array of METIS options as described in Section 5.4 of the METIS manual.

      idx_t options[METIS_NOPTIONS];

      // future-proof (or just in case we forgot to set some options explicitly)

      METIS_SetDefaultOptions(options);

      // set METIS options from parameters

      setMETISoptions(options, parameters);



      if( nparts == 1 ) {

         // k-way partitioning from Metis fails for nparts == 1 (segfault),

         // RB succeeds but produces nonsense

         part_array.setValue( 0 );

      }

      else {

         if( options[METIS_OPTION_PTYPE] == METIS_PTYPE_KWAY ) {

            std::cout << "Running METIS_PartGraphKway..." << std::endl;

            status = METIS_PartGraphKway(&nvtxs, &ncon, xadj, adjncy, vwgt, vsize, adjwgt, &nparts, tpwgts, ubvec, options, &objval, part);

         }

         else {

            std::cout << "Running METIS_PartGraphRecursive..." << std::endl;

            status = METIS_PartGraphRecursive(&nvtxs, &ncon, xadj, adjncy, vwgt, vsize, adjwgt, &nparts, tpwgts, ubvec, options, &objval, part);

         }



         switch( status )

         {

            case METIS_OK: break;

            case METIS_ERROR_INPUT:

               throw std::runtime_error( "METIS_PartGraph failed due to an input error." );

            case METIS_ERROR_MEMORY:

               throw std::runtime_error( "METIS_PartGraph failed due to a memory allocation error." );

            case METIS_ERROR:

            default:

               throw std::runtime_error( "METIS_PartGraph failed with an unspecified error." );

         }

      }



      // deallocate auxiliary vectors

      connectivity.clear();

      connectivity.shrink_to_fit();

      offsets.clear();

      offsets.shrink_to_fit();



      const unsigned ghost_levels = parameters.getParameter< unsigned >( "ghost-levels" );

      const std::string pvtuFileName = parameters.template getParameter< String >( "output-file" );

      decompose_and_save( mesh, nparts, part_array, shared_xadj, shared_adjncy, ncommon, ghost_levels, pvtuFileName );

   }



   static void

void

decompose_and_save( const Mesh& mesh,

                    const unsigned nparts,

                    const MetisIndexArray& part,
@@ -757,7 +608,154 @@ struct DecomposeMesh 
      }

   }

}

};



template< typename Mesh >

void run( const Mesh& mesh, const Config::ParameterContainer& parameters )

{

   using Index = typename Mesh::GlobalIndexType;



   // warn if input mesh has 64-bit indices, but METIS uses only 32-bit indices

   if( IDXTYPEWIDTH == 32 && sizeof(Index) > 4 )

      std::cerr << "Warning: the input mesh uses 64-bit indices, but METIS was compiled only with 32-bit indices. "

                   "Decomposition may not work correctly if the index values overflow the 32-bit type." << std::endl;



   // get the mesh connectivity information in a format suitable for METIS. Actually, the same

   // format is used by the XML-based VTK formats - the only difference is that METIS requires

   // `offsets` to start with 0.

   std::vector< idx_t > connectivity, offsets;

   offsets.push_back(0);

   const Index cellsCount = mesh.template getEntitiesCount< typename Mesh::Cell >();

   for( Index i = 0; i < cellsCount; i++ ) {

      const auto& entity = mesh.template getEntity< typename Mesh::Cell >( i );

      const Index subvertices = entity.template getSubentitiesCount< 0 >();

      for( Index j = 0; j < subvertices; j++ )

         connectivity.push_back( entity.template getSubentityIndex< 0 >( j ) );

      offsets.push_back( connectivity.size() );

   }



   // number of elements (cells)

   idx_t ne = mesh.template getEntitiesCount< typename Mesh::Cell >();

   // number of nodes (vertices)

   idx_t nn = mesh.template getEntitiesCount< 0 >();

   // pointers to arrays storing the mesh in a CSR-like format

   idx_t* eptr = offsets.data();

   idx_t* eind = connectivity.data();

   // Specifies the number of common nodes that two elements must have in order to put an edge

   // between them in the dual graph.

   idx_t ncommon = Mesh::getMeshDimension();

   if( parameters.checkParameter( "min-common-vertices" ) )

      ncommon = parameters.template getParameter< unsigned >( "min-common-vertices" );

   // numbering scheme for the adjacency structure of a graph or element-node structure of a mesh

   // (0 is C-style, 1 is Fortran-style)

   idx_t numflag = 0;

   // These arrays store the adjacency structure of the generated dual graph. The format is of

   // the adjacency structure is described in Section 5.5 of the METIS manual. Memory for these

   // arrays is allocated by METIS’ API using the standard malloc function. It is the

   // responsibility of the application to free this memory by calling free. METIS provides the

   // METIS_Free that is a wrapper to C’s free function.

   idx_t* xadj = nullptr;

   idx_t* adjncy = nullptr;



   // We could use METIS_PartMeshDual directly instead of METIS_MeshToDual + METIS_PartGraph*,

   // but we need to reuse the dual graph for the generation of ghost cells.

   std::cout << "Running METIS_MeshToDual..." << std::endl;

   int status = METIS_MeshToDual(&ne, &nn, eptr, eind, &ncommon, &numflag, &xadj, &adjncy);



   // wrap xadj and adjncy with shared_ptr

   std::shared_ptr<idx_t> shared_xadj {xadj, METIS_Free};

   std::shared_ptr<idx_t> shared_adjncy {adjncy, METIS_Free};



   switch( status )

   {

      case METIS_OK: break;

      case METIS_ERROR_INPUT:

         throw std::runtime_error( "METIS_MeshToDual failed due to an input error." );

      case METIS_ERROR_MEMORY:

         throw std::runtime_error( "METIS_MeshToDual failed due to a memory allocation error." );

      case METIS_ERROR:

      default:

         throw std::runtime_error( "METIS_MeshToDual failed with an unspecified error." );

   }



   // The number of vertices in the graph.

   idx_t nvtxs = ne;

   // The number of balancing constraints. It should be at least 1.

   idx_t ncon = 1;

   // An array of size `ne` specifying the weights of the elements. A NULL value can be passed

   // to indicate that all elements have an equal weight.

   idx_t* vwgt = nullptr;

   // An array of size `ne` specifying the size of the elements that is used for computing the

   // total communication volume as described in Section 5.7 of the METIS manual. A NULL value

   // can be passed when the objective is cut or when all elements have an equal size.

   idx_t* vsize = nullptr;

   // The weights of the edges as described in Section 5.5 of the METIS manual.

   idx_t* adjwgt = nullptr;  // METIS_PartMeshDual uses NULL too

   // The number of parts to partition the mesh.

   idx_t nparts = parameters.template getParameter< unsigned >( "subdomains" );

   // This is an array of size nparts that specifies the desired weight for each partition. The

   // target partition weight for the i-th partition is specified at tpwgts[i] (the numbering for

   // the partitions starts from 0). The sum of the tpwgts[] entries must be 1.0. A NULL value

   // can be passed to indicate that the graph should be equally divided among the partitions.

   real_t* tpwgts = nullptr;

   // This is an array of size ncon that specifies the allowed load imbalance tolerance for each

   // constraint. For the i-th partition and j-th constraint the allowed weight is the

   // ubvec[j]*tpwgts[i*ncon+j] fraction of the j-th’s constraint total weight. The load

   // imbalances must be greater than 1.0. A NULL value can be passed indicating that the load

   // imbalance tolerance for each constraint should be 1.001 (for ncon=1) or 1.01 (for ncon>1).

   real_t* ubvec = nullptr;  // METIS_PartMeshDual uses NULL too

   // Upon successful completion, this variable stores either the edgecut or the total

   // communication volume of the dual graph’s partitioning.

   idx_t objval = 0;

   // Array of size `ne` that upon successful completion stores the partition array for the

   // elements of the mesh.

   MetisIndexArray part_array( nvtxs );

   idx_t* part = part_array.getData();



   // Array of METIS options as described in Section 5.4 of the METIS manual.

   idx_t options[METIS_NOPTIONS];

   // future-proof (or just in case we forgot to set some options explicitly)

   METIS_SetDefaultOptions(options);

   // set METIS options from parameters

   setMETISoptions(options, parameters);



   if( nparts == 1 ) {

      // k-way partitioning from Metis fails for nparts == 1 (segfault),

      // RB succeeds but produces nonsense

      part_array.setValue( 0 );

   }

   else {

      if( options[METIS_OPTION_PTYPE] == METIS_PTYPE_KWAY ) {

         std::cout << "Running METIS_PartGraphKway..." << std::endl;

         status = METIS_PartGraphKway(&nvtxs, &ncon, xadj, adjncy, vwgt, vsize, adjwgt, &nparts, tpwgts, ubvec, options, &objval, part);

      }

      else {

         std::cout << "Running METIS_PartGraphRecursive..." << std::endl;

         status = METIS_PartGraphRecursive(&nvtxs, &ncon, xadj, adjncy, vwgt, vsize, adjwgt, &nparts, tpwgts, ubvec, options, &objval, part);

      }



      switch( status )

      {

         case METIS_OK: break;

         case METIS_ERROR_INPUT:

            throw std::runtime_error( "METIS_PartGraph failed due to an input error." );

         case METIS_ERROR_MEMORY:

            throw std::runtime_error( "METIS_PartGraph failed due to a memory allocation error." );

         case METIS_ERROR:

         default:

            throw std::runtime_error( "METIS_PartGraph failed with an unspecified error." );

      }

   }



   // deallocate auxiliary vectors

   connectivity.clear();

   connectivity.shrink_to_fit();

   offsets.clear();

   offsets.shrink_to_fit();



   const unsigned ghost_levels = parameters.getParameter< unsigned >( "ghost-levels" );

   const std::string pvtuFileName = parameters.template getParameter< String >( "output-file" );

   decompose_and_save( mesh, nparts, part_array, shared_xadj, shared_adjncy, ncommon, ghost_levels, pvtuFileName );

}



int main( int argc, char* argv[] )

{
@@ -780,7 +778,7 @@ int main( int argc, char* argv[] ) 
   auto wrapper = [&] ( const auto& reader, auto&& mesh )

   {

      using MeshType = std::decay_t< decltype(mesh) >;

      DecomposeMesh< MeshType >::run( std::forward<MeshType>(mesh), parameters );

      run( std::forward<MeshType>(mesh), parameters );

      return true;

   };

   return ! Meshes::resolveAndLoadMesh< DecomposeMeshConfigTag, Devices::Host >( wrapper, inputFileName, inputFileFormat );