SF_numbering.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) NumeriCor GmbH
// SPDX-License-Identifier: Apache-2.0
 
#ifndef _SF_NUMBERING_H
#define _SF_NUMBERING_H
 
#include <cassert>
#include <cstddef>
#include <iostream>
#include <mpi.h>
 
#include "hashmap.hpp"
 
#include "SF_container.h"
#include "SF_vector.h"
 
namespace SF {
 
template<class T, class S>
class numbering {
public:
 
  virtual void operator() (meshdata<T, S> & mesh) = 0;
};
 
 
template<class T, class S>
class submesh_numbering : public numbering<T, S>
{
public:
  inline void renumber_sorted_ascending(meshdata<T, S> & submesh, SF_nbr in_nbr, SF_nbr out_nbr)
  {
    MPI_Comm comm = submesh.comm;
 
    int size, rank;
    MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
 
    // determine global min and max indices
    const vector<T> & rnod = submesh.get_numbering(in_nbr);
    T gmax = global_max(rnod, comm);
    T gmin = global_min(rnod, comm);
 
    // block size
    T bsize = (gmax - gmin) / size + 1;
 
    // distribute indices uniquely and linearly ascending (w.r.t. rank index) --------------------
    vector<T> dest(rnod.size()), sbuff(rnod);
    for(size_t i=0; i<dest.size(); i++)
      dest[i] = (rnod[i] - gmin) / bsize;
    binary_sort_copy(dest, sbuff);
 
    commgraph<size_t> grph;
    grph.configure(dest, comm);
 
    size_t rsize = sum(grph.rcnt);
    vector<T> rbuff(rsize);
    MPI_Exchange(grph, sbuff, rbuff, comm);
 
    // we need a data structure to store where we got the indices from
    vector<T> proc(rsize), acc_cnt(rsize, T(1));
    grph.source_ranks(proc);
 
    // compute local node set
    binary_sort_copy(rbuff, proc);
    unique_accumulate(rbuff, acc_cnt);
    rsize = rbuff.size();
 
    // communicate local sizes
    MPI_Allgather(&rsize, sizeof(size_t), MPI_BYTE, grph.rcnt.data(), sizeof(size_t), MPI_BYTE, comm);
    dsp_from_cnt(grph.rcnt, grph.rdsp);
 
    // now the new numbering can be derived from the data layout
    vector<T> new_nbr(rsize);
    interval(new_nbr, grph.rdsp[rank], grph.rdsp[rank+1]);
 
    // we send the old and new numberings back -------------------------------------------------
    sbuff.resize(proc.size());
 
    // expand the old numbering into the send buffer
    for(size_t i=0, idx=0; i<acc_cnt.size(); i++)
      for(T j=0; j<acc_cnt[i]; j++, idx++) sbuff[idx] = rbuff[i];
    dest.assign(proc.begin(), proc.end());
    binary_sort_copy(dest, sbuff); // permute for sending
    grph.configure(dest, comm);    // reconfigure comm graph
    rbuff.resize(sum(grph.rcnt));  // resize receive buffer
    MPI_Exchange(grph, sbuff, rbuff, comm);  // communicate
 
    vector<T> old_idx(rbuff);
 
    // expand the new numbering into the send buffer
    for(size_t i=0, idx=0; i<acc_cnt.size(); i++)
      for(T j=0; j<acc_cnt[i]; j++, idx++) sbuff[idx] = new_nbr[i];
    dest.assign(proc.begin(), proc.end());
    binary_sort_copy(dest, sbuff); // permute for sending
    MPI_Exchange(grph, sbuff, rbuff, comm);  // communicate
 
    vector<T> new_idx(rbuff);
 
    // add a new numbering to submesh -----------------------------------------
    // create mapping between original and new numbering
    index_mapping<T> map(old_idx, new_idx);
    vector<T> & snod = submesh.register_numbering(out_nbr);
    snod.resize(rnod.size());
    // use mapping for new numbering
    for(size_t i=0; i < rnod.size(); i++) snod[i] = map.forward_map(rnod[i]);
  }
 
  inline void operator() (meshdata<T, S> & submesh)
  {
    renumber_sorted_ascending(submesh, NBR_REF, NBR_SUBMESH);
    renumber_sorted_ascending(submesh, NBR_ELEM_REF, NBR_ELEM_SUBMESH);
  }
 
};
 
template<class T> inline
void reindex_cuthill_mckee(const vector<T> & n2n_cnt,
                           const vector<T> & n2n_dsp,
                           const vector<T> & n2n_con,
                           const bool reverse,
                           hashmap::unordered_map<T,T> & old2new)
{
  /*
   * We use a Cuthill-McKee style breadth-first traversal to generate a permutation
   * of the original indices. This permutation is then renumbered
   * either forward or reverse.
   *
   * The permutation is generated as following:
   * Initially, an arbitrary node is selected. It is inserted in the permutation
   * and in the imaginary set R (represented by a boolean vector).
   * Then we loop over each node in perm, adding all connected nodes, that are not yet
   * in R, in ascending order to perm and R.
   *
   * Note that this is a simplified variant of canonical Cuthill-McKee: the traversal
   * seeds at node 0 and orders neighbours by index, instead of using the degree and
   * pseudo-peripheral node heuristics. It therefore reduces the bandwidth less than a
   * canonical implementation would.
   */
  size_t nnod = n2n_cnt.size();
  if(nnod == 0) return;   // empty partition: nothing to renumber
 
  vector<bool>   inR(nnod, false);
  vector<T> perm(nnod);
  perm.resize(1); perm[0] = 0; inR[0] = true;
  T pidx=0;
  // nodes below this were all visited, so seeding never has to look at them again
  size_t seed = 0;
  size_t ncomp = 1;
 
  while(perm.size() < nnod)
  {
    T nidx;
    if (pidx < T(perm.size())) nidx = perm[pidx++];
    else {
      // the traversal has exhausted the current connected component. Seed the next
      // one with the first unvisited node and consume it as if it had been reached.
      while(inR[seed] == true) seed++;
      nidx = T(seed);
      inR[seed] = true;
      perm.push_back(nidx);
      pidx++;
      ncomp++;
    }
    T start = n2n_dsp[nidx], stop = start + n2n_cnt[nidx];
    vector<T> adj(stop - start);
 
    T widx=0;
    for(T j = start; j<stop; j++) {
      T cidx = n2n_con[j];
      if( ! inR[cidx] ) {
        adj[widx++] = cidx;
        inR[cidx] = true;
      }
    }
 
    adj.resize(widx);
    binary_sort(adj);
    perm.append(adj.begin(), adj.end());
  }
 
  // Reported once per rank rather than per component. Several components are legitimate
  // for meshes that are decoupled by design, such as the cells of an EMI mesh, but they
  // also weaken the ordering: the traversal restarts and cannot relate one component to
  // the next, so the bandwidth reduction only holds within each of them.
  if(ncomp > 1)
    std::cerr << "Note: the mesh graph splits into " << ncomp << " connected components; "
              << "renumbering treats each of them separately." << std::endl;
 
  if(!reverse) {
    // Cuthill-McKee
    for(size_t i=0; i<perm.size(); i++) old2new[perm[i]] = i;
  } else {
    // Reverse Cuthill-McKee
    for(size_t i=0, j=perm.size()-1; i<perm.size(); i++, j--) old2new[perm[j]] = i;
  }
}
 
 
template<class T, class S>
class petsc_numbering : public numbering<T, S>
{
private:
  overlapping_layout<T> & _pl;
  bool _use_rcm;
public:
  petsc_numbering(overlapping_layout<T> & pl, bool use_rcm = false) : _pl(pl), _use_rcm(use_rcm)
  {}
 
  inline void operator() (meshdata<T, S> & mesh)
  {
    vector<T> & petsc_idx = mesh.register_numbering(NBR_PETSC);
    petsc_idx.assign(_pl.num_local_idx(), -1);
 
    int size, rank;
    MPI_Comm comm = mesh.comm;
    MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
 
 
    const vector<T> & alg_dsp = _pl.algebraic_layout();
    const vector<T> & alg_nod = _pl.algebraic_nodes();
    const T alg_start = alg_dsp[rank], alg_end = alg_dsp[rank+1];
 
    if(!_use_rcm) {
      // trivial ascending numbering
      for(size_t i=0; i<alg_nod.size(); i++)
        petsc_idx[alg_nod[i]] = alg_start + i;
    }
    // a rank may own no algebraic nodes at all. The restricted graph is then empty and
    // the connectivity kernels would index into zero-sized buffers, so skip straight to
    // the communication step below.
    else if(alg_nod.size() > 0) {
      // reversed breadth-first renumbering on the restricted submesh of local algebraic nodes.
      meshdata<T,S> restr_mesh;
      vector<T> cnt(mesh.dsp.size() - 1, 0);
      hashmap::unordered_map<T,T> nodal2restr;
      restr_mesh.con.reserve(mesh.con.size());
 
      for(size_t i=0; i<alg_nod.size(); i++) {
        nodal2restr[alg_nod[i]] = i;
      }
 
      for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++)
      {
        T start = mesh.dsp[eidx], stop = mesh.dsp[eidx+1];
        for(T i = start; i<stop; i++) {
          if(nodal2restr.count(mesh.con[i])) {
            cnt[eidx]++;
            restr_mesh.con.push_back(nodal2restr[mesh.con[i]]);
          }
        }
      }
 
      dsp_from_cnt(cnt, restr_mesh.dsp);
 
      vector<T> n2n_cnt, n2n_dsp, n2n_con;
      hashmap::unordered_map<T,T> old2new;
 
      nodal_connectivity_graph(restr_mesh, n2n_cnt, n2n_con);
      dsp_from_cnt(n2n_cnt, n2n_dsp);
 
      reindex_cuthill_mckee(n2n_cnt, n2n_dsp, n2n_con, true, old2new);
 
      assert(old2new.size() == alg_nod.size());
 
      // Apply the reordering to the algebraic node list, not to the indices handed out.
      // Both give a node the same PETSc index, but only this way does the position of a
      // node in algebraic_nodes() keep matching its slot in the PETSc vector, which is
      // how the ionic model masks and the asynchronous output address the local arrays.
      // at() rather than operator[]: a node missing from the graph would otherwise be
      // silently mapped to 0, leaving perm_nod partly unwritten and corrupting the node
      // list the layout is about to adopt.
      vector<T> perm_nod(alg_nod.size());
      for(size_t i=0; i<alg_nod.size(); i++)
        perm_nod[old2new.at(i)] = alg_nod[i];
 
      _pl.permute_algebraic_nodes(perm_nod);
 
      // alg_nod refers to the layout's list, so it reflects the new order here
      for(size_t i=0; i<alg_nod.size(); i++)
        petsc_idx[alg_nod[i]] = alg_start + i;
    }
    // communicate the new numbering so that all nodes in the local DD domain
    // -- also the ones not in the algebraic range -- have a new index
 
    _pl.reduce(petsc_idx, "max");
 
    bool print_indexing = false;
    if(print_indexing) {
      for(int pid = 0; pid < size; pid++) {
        if(pid == rank) {
          for(size_t eidx = 0; eidx < 100; eidx++) {
            printf("rank %d elem %d: ", rank, int(eidx));
            for(T j = mesh.dsp[eidx]; j < mesh.dsp[eidx+1] - 1; j++) {
              const T c = mesh.con[j];
              printf("%lld ", static_cast<long long>(petsc_idx[c]));
            }
            printf("%lld\n",
                   static_cast<long long>(petsc_idx[mesh.con[mesh.dsp[eidx+1] - 1]]));
          }
        }
        MPI_Barrier(PETSC_COMM_WORLD);
      }
    }
  }
};
 
}
 
#endif