oomph::METIS Namespace Reference

Namespace for METIS graph partitioning routines. More...

Typedefs
typedef void(*	ErrorToWeightFctPt) (const double &spatial_error, const double &max_error, const double &min_error, int &weight)

Functions
void	default_error_to_weight_fct (const double &spatial_error, const double &max_error, const double &min_error, int &weight)
void	uniform_partition_mesh (Problem *problem_pt, const unsigned &ndomain, Vector< unsigned > &element_domain)
void	partition_mesh (Problem *problem_pt, const unsigned &ndomain, const unsigned &objective, Vector< unsigned > &element_domain)
void	partition_mesh (OomphCommunicator *comm_pt, Mesh *mesh_pt, const unsigned &ndomain, const unsigned &objective, Vector< unsigned > &element_domain)

Variables
ErrorToWeightFctPt	Error_to_weight_fct_pt = &default_error_to_weight_fct

Detailed Description

Namespace for METIS graph partitioning routines.

Typedef Documentation

ErrorToWeightFctPt

typedef void(* oomph::METIS::ErrorToWeightFctPt) (const double &spatial_error, const double &max_error, const double &min_error, int &weight)

Typedef for function pointer to to function that translates spatial error into weight for METIS partitioning.

Function Documentation

default_error_to_weight_fct()

void oomph::METIS::default_error_to_weight_fct	(	const double &	spatial_error,
		const double &	max_error,
		const double &	min_error,
		int &	weight
	)

Default function that translates spatial error into weight for METIS partitioning (unit weight regardless of input).

Default function that translates spatial error into weight for METIS partitioning (unit weight regardless of input)

    {
      weight = 1;
    }

Referenced by partition_mesh().

partition_mesh() [1/2]

void oomph::METIS::partition_mesh	(	OomphCommunicator *	comm_pt,
		Mesh *	mesh_pt,
		const unsigned &	ndomain,
		const unsigned &	objective,
		Vector< unsigned > &	element_domain
	)

Use METIS to assign each element to a domain. On return, element_domain[ielem] contains the number of the domain [0,1,...,ndomain-1] to which element ielem has been assigned.

• objective=0: minimise edgecut.
• objective=1: minimise total communications volume.

Partioning is based on nodal graph of mesh.

partition_mesh() [2/2]

void oomph::METIS::partition_mesh	(	Problem *	problem_pt,
		const unsigned &	ndomain,
		const unsigned &	objective,
		Vector< unsigned > &	element_domain
	)

Use METIS to assign each element to a domain. On return, element_domain[ielem] contains the number of the domain [0,1,...,ndomain-1] to which element ielem has been assigned.

• objective=0: minimise edgecut.
• objective=1: minimise total communications volume.

Partioning is based on nodal graph of mesh.

Use METIS to assign each element to a domain. On return, element_domain[ielem] contains the number of the domain [0,1,...,ndomain-1] to which element ielem has been assigned.

• objective=0: minimise edgecut.
• objective=1: minimise total communications volume.

Partioning is based on dual graph of mesh.

  {
    // Global mesh
    Mesh* mesh_pt = problem_pt->mesh_pt();
 
    // Number of elements
    unsigned nelem = mesh_pt->nelement();
 
#ifdef PARANOID
    if (nelem != element_domain.size())
    {
      std::ostringstream error_stream;
      error_stream << "element_domain Vector has wrong length " << nelem << " "
                   << element_domain.size() << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Setup dual graph
    //------------------
 
    // Start timer
    clock_t cpu_start = clock();
 
    // Container to collect all elements associated with given global eqn number
    std::map<unsigned, std::set<unsigned>> elements_connected_with_global_eqn;
 
    // Container for all unique global eqn numbers
    std::set<unsigned> all_global_eqns;
 
    // Loop over all elements
    for (unsigned e = 0; e < nelem; e++)
    {
      GeneralisedElement* el_pt = mesh_pt->element_pt(e);
 
      // Add all global eqn numbers
      unsigned ndof = el_pt->ndof();
      for (unsigned j = 0; j < ndof; j++)
      {
        // Get global eqn number
        unsigned eqn_number = el_pt->eqn_number(j);
        elements_connected_with_global_eqn[eqn_number].insert(e);
        all_global_eqns.insert(eqn_number);
      }
    }
 
    // Now reverse the lookup scheme to find out all elements
    // that are connected because they share the same global eqn
    Vector<std::set<unsigned>> connected_elements(nelem);
 
    // Counter for total number of entries in connected_elements structure
    unsigned count = 0;
 
    // Loop over all global eqns
    for (std::set<unsigned>::iterator it = all_global_eqns.begin();
         it != all_global_eqns.end();
         it++)
    {
      // Get set of elements connected with this data item
      std::set<unsigned> elements = elements_connected_with_global_eqn[*it];
 
      // Double loop over connnected elements: Everybody's connected to
      // everybody
      for (std::set<unsigned>::iterator it1 = elements.begin();
           it1 != elements.end();
           it1++)
      {
        for (std::set<unsigned>::iterator it2 = elements.begin();
             it2 != elements.end();
             it2++)
        {
          if ((*it1) != (*it2))
          {
            connected_elements[(*it1)].insert(*it2);
          }
        }
      }
    }
 
 
    // Now convert into C-style packed array for interface with METIS
    int* xadj = new int[nelem + 1];
    Vector<int> adjacency_vector;
 
    // Reserve (too much) space
    adjacency_vector.reserve(count);
 
    // Initialise counters
    unsigned ientry = 0;
 
    // Loop over all elements
    for (unsigned e = 0; e < nelem; e++)
    {
      // First entry for current element
      xadj[e] = ientry;
 
      // Loop over elements that are connected to current element
      typedef std::set<unsigned>::iterator IT;
      for (IT it = connected_elements[e].begin();
           it != connected_elements[e].end();
           it++)
      {
        // Copy into adjacency array
        adjacency_vector.push_back(*it);
 
        // We've just made another entry
        ientry++;
      }
 
      // Entry after last entry for current element:
      xadj[e + 1] = ientry;
    }
 
    // End timer
    clock_t cpu_end = clock();
 
    // Doc
    double cpu0 = double(cpu_end - cpu_start) / CLOCKS_PER_SEC;
    oomph_info
      << "CPU time for setup of METIS data structures            [nelem="
      << nelem << "]: " << cpu0 << " sec" << std::endl;
 
 
    // Call METIS graph partitioner
    //-----------------------------
 
    // Start timer
    cpu_start = clock();
 
    // Number of vertices in graph
    int nvertex = nelem;
 
    // No vertex weights
    int* vwgt = 0;
 
    // No edge weights
    int* adjwgt = 0;
 
    // Flag indicating that graph isn't weighted: 0; vertex weights only: 2
    // Note that wgtflag==2 requires nodal weights to be stored in vwgt.
    int wgtflag = 0;
 
    // Use C-style numbering (first array entry is zero)
    int numflag = 0;
 
    // Number of desired partitions
    int nparts = ndomain;
 
    // Use default options
    int* options = new int[10];
    options[0] = 0;
 
#ifdef OOMPH_TRANSITION_TO_VERSION_3
    switch (objective)
    {
      case 0:
        // Edge-cut minimization
        options[0] = METIS_OBJTYPE_CUT;
        break;
 
      case 1:
        // communication volume minimisation
        options[0] = METIS_OBJTYPE_VOL;
        break;
 
      default:
        std::ostringstream error_stream;
        error_stream << "Wrong objective for METIS. objective = " << objective
                     << std::endl;
 
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Number of cut edges in graph
    int* edgecut = new int[nelem];
 
    // Array containing the partition information
    int* part = new int[nelem];
 
    // Can we get an error estimate?
 
    unsigned n_mesh = problem_pt->nsub_mesh();
 
    if (n_mesh == 0)
    {
      RefineableMeshBase* mmesh_pt = dynamic_cast<RefineableMeshBase*>(mesh_pt);
      if (mmesh_pt != 0)
      {
        // Bias distribution?
        if (Error_to_weight_fct_pt != &default_error_to_weight_fct)
        {
          oomph_info
            << "Biasing element distribution via spatial error estimate\n";
 
          // Adjust flag and provide storage for weights
          wgtflag = 2;
          vwgt = new int[nelem];
 
          // Get error for all elements
          Vector<double> elemental_error(nelem);
          mmesh_pt->spatial_error_estimator_pt()->get_element_errors(
            mesh_pt, elemental_error);
 
          double max_error =
            *(std::max_element(elemental_error.begin(), elemental_error.end()));
          double min_error =
            *(std::min_element(elemental_error.begin(), elemental_error.end()));
 
          // Bias weights
          int weight = 1;
          for (unsigned e = 0; e < nelem; e++)
          {
            // Translate error into weight
            Error_to_weight_fct_pt(
              elemental_error[e], max_error, min_error, weight);
            vwgt[e] = weight;
          }
        }
      }
    }
    else // There are submeshes
    {
      // Are any of the submeshes refineable?
      bool refineable_submesh_exists = false;
      // Vector to store "start and end point" for loops in submeshes
      Vector<unsigned> loop_helper(n_mesh + 1);
      loop_helper[0] = 0;
 
      // Loop over submeshes
      for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
      {
        // Store the end of the loop
        loop_helper[i_mesh + 1] =
          problem_pt->mesh_pt(i_mesh)->nelement() + loop_helper[i_mesh];
 
        RefineableMeshBase* mmesh_pt =
          dynamic_cast<RefineableMeshBase*>(problem_pt->mesh_pt(i_mesh));
        if (mmesh_pt != 0)
        {
          refineable_submesh_exists = true;
        }
      }
 
      // If a refineable submesh exists
      if (refineable_submesh_exists)
      {
        // Bias distribution?
        if (Error_to_weight_fct_pt != &default_error_to_weight_fct)
        {
          oomph_info
            << "Biasing element distribution via spatial error estimate\n";
 
          // Adjust flag and provide storage for weights
          wgtflag = 2;
          vwgt = new int[nelem];
 
          // Loop over submeshes
          for (unsigned i_mesh = 0; i_mesh < n_mesh; i_mesh++)
          {
            RefineableMeshBase* mmesh_pt =
              dynamic_cast<RefineableMeshBase*>(problem_pt->mesh_pt(i_mesh));
            if (mmesh_pt != 0)
            {
              // Get error for all elements
              unsigned nsub_elem =
                loop_helper[i_mesh + 1] - loop_helper[i_mesh];
              Vector<double> elemental_error(nsub_elem);
              mmesh_pt->spatial_error_estimator_pt()->get_element_errors(
                problem_pt->mesh_pt(i_mesh), elemental_error);
 
              double max_error = *(std::max_element(elemental_error.begin(),
                                                    elemental_error.end()));
              double min_error = *(std::min_element(elemental_error.begin(),
                                                    elemental_error.end()));
 
              // Bias weights
              int weight = 1;
              unsigned start = loop_helper[i_mesh];
              unsigned end = loop_helper[i_mesh + 1];
              for (unsigned e = start; e < end; e++)
              {
                unsigned error_index = e - start;
                // Translate error into weight
                Error_to_weight_fct_pt(
                  elemental_error[error_index], max_error, min_error, weight);
                vwgt[e] = weight;
              }
            }
            else // This mesh is not refineable
            {
              // There's no error estimator, so use the default weight
              int weight = 1;
              unsigned start = loop_helper[i_mesh];
              unsigned end = loop_helper[i_mesh + 1];
              for (unsigned e = start; e < end; e++)
              {
                vwgt[e] = weight;
              }
            }
          }
        }
      }
    }
 
#ifdef OOMPH_TRANSITION_TO_VERSION_3
 
    // Call partitioner
    METIS_PartGraphKway(&nvertex,
                        xadj,
                        &adjacency_vector[0],
                        vwgt,
                        adjwgt,
                        &wgtflag,
                        &numflag,
                        &nparts,
                        options,
                        edgecut,
                        part);
#else
    // original code to delete in version 3
 
    // Call partitioner
    if (objective == 0)
    {
      // Partition with the objective of minimising the edge cut
      METIS_PartGraphKway(&nvertex,
                          xadj,
                          &adjacency_vector[0],
                          vwgt,
                          adjwgt,
                          &wgtflag,
                          &numflag,
                          &nparts,
                          options,
                          edgecut,
                          part);
    }
    else if (objective == 1)
    {
      // Partition with the objective of minimising the total communication
      // volume
      METIS_PartGraphVKway(&nvertex,
                           xadj,
                           &adjacency_vector[0],
                           vwgt,
                           adjwgt,
                           &wgtflag,
                           &numflag,
                           &nparts,
                           options,
                           edgecut,
                           part);
    }
    else
    {
      std::ostringstream error_stream;
      error_stream << "Wrong objective for METIS. objective = " << objective
                   << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
#ifdef PARANOID
    std::vector<bool> done(nparts, false);
#endif
 
    // Copy across
    for (unsigned e = 0; e < nelem; e++)
    {
      element_domain[e] = part[e];
#ifdef PARANOID
      done[part[e]] = true;
#endif
    }
 
 
#ifdef PARANOID
    // Check
    std::ostringstream error_stream;
    bool shout = false;
    for (int p = 0; p < nparts; p++)
    {
      if (!done[p])
      {
        shout = true;
        error_stream << "No elements on processor " << p
                     << "when trying to partition " << nelem << "elements over "
                     << nparts << " processors!\n";
      }
    }
    if (shout)
    {
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
 
    // End timer
    cpu_end = clock();
 
    // Doc
    double cpu1 = double(cpu_end - cpu_start) / CLOCKS_PER_SEC;
    oomph_info
      << "CPU time for METIS mesh partitioning                   [nelem="
      << nelem << "]: " << cpu1 << " sec" << std::endl;
 
 
    // Cleanup
    delete[] xadj;
    delete[] part;
    delete[] edgecut;
    delete[] options;
  }

References default_error_to_weight_fct(), e(), oomph::Mesh::element_pt(), Eigen::placeholders::end, oomph::GeneralisedElement::eqn_number(), Error_to_weight_fct_pt, oomph::ErrorEstimator::get_element_errors(), j, MeshRefinement::max_error, oomph::Problem::mesh_pt(), oomph::METIS_PartGraphKway(), oomph::METIS_PartGraphVKway(), MeshRefinement::min_error, oomph::GeneralisedElement::ndof(), oomph::Mesh::nelement(), oomph::Problem::nsub_mesh(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, p, oomph::RefineableMeshBase::spatial_error_estimator_pt(), and oomph::CumulativeTimings::start().

uniform_partition_mesh()

void oomph::METIS::uniform_partition_mesh	(	Problem *	problem_pt,
		const unsigned &	ndomain,
		Vector< unsigned > &	element_domain
	)

Partition mesh uniformly by dividing elements equally over the partitions, in the order in which they are returned by problem. On return, element_domain[ielem] contains the number of the domain [0,1,...,ndomain-1] to which element ielem has been assigned.

  {
    // Number of elements
    unsigned nelem = problem_pt->mesh_pt()->nelement();
 
#ifdef PARANOID
    if (nelem != element_domain.size())
    {
      std::ostringstream error_stream;
      error_stream << "element_domain Vector has wrong length " << nelem << " "
                   << element_domain.size() << std::endl;
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
 
    // Uniform partitioning
    unsigned nel_per_domain = int(float(nelem) / float(ndomain));
    for (unsigned ielem = 0; ielem < nelem; ielem++)
    {
      unsigned idomain = unsigned(float(ielem) / float(nel_per_domain));
      element_domain[ielem] = idomain;
    }
  }

References int(), oomph::Problem::mesh_pt(), oomph::Mesh::nelement(), OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

Variable Documentation

Error_to_weight_fct_pt

ErrorToWeightFctPt oomph::METIS::Error_to_weight_fct_pt = &default_error_to_weight_fct

Function pointer to to function that translates spatial error into weight for METIS partitioning.

Referenced by partition_mesh().