IndexDistribution.inc.hpp

#pragma once

#include <utility>

#include <dune/common/timer.hh>

#if HAVE_MPI
#include <amdis/common/parallel/Collective.hpp>
#include <amdis/common/parallel/RequestOperations.hpp>
#endif

namespace AMDiS {

#if HAVE_MPI

template <class GlobalId, class LI>
  template <class Basis>
void PetscParallelIndexDistribution<GlobalId,LI>::
update(Basis const& basis)
{
  Dune::Timer t;

  // clear all vectors and reset sizes
  reset();
  buildParallelIndexSet(basis, *indexSet_);

  remoteIndices_->setIndexSets(*indexSet_, *indexSet_, world_);
  remoteIndices_->template rebuild<true>();

  // 1. insert and number owner DOFs
  globalIndices_.resize(indexSet_->size(), 0);
  owner_.resize(indexSet_->size(), false);
  localSize_ = 0;
  ghostSize_ = 0;
  for (auto const& ip : *indexSet_) {
    if (ip.local().attribute() == Attribute::owner) {
      size_type idx = ip.local();
      globalIndices_[idx] = localSize_++;
      owner_[idx] = true;
    } else {
      ghostSize_++;
    }
  }

  // communicate the local sizes from all processors
  Mpi::all_gather(world_, localSize_, sizes_);

  // at which global index do the local partitions start
  starts_.resize(world_.size() + 1);
  starts_[0] = 0;
  std::partial_sum(sizes_.begin(), sizes_.end(), starts_.begin()+1);
  globalSize_ = starts_.back();

  // update the global index for all local indices by shifting by global start position
  for (auto& i : globalIndices_)
    i += starts_[world_.rank()];

  // build up the communication of overlap DOFs that do not yet have a global index
  // assigned. Therefore send the global index for all already computed owner DOFs
  // to the neighboring remote processors. And receive from those their owner DOFs
  // global indices.
  using GlobalAssoc = std::pair<typename IndexSet::GlobalIndex, size_type>; // {globalId, globalIndex}
  std::vector<std::vector<GlobalAssoc>> sendList(world_.size());
  std::vector<std::size_t> receiveList(world_.size(), 0);

  // Communicate attributes at the interface
  for (auto const& rim : *remoteIndices_) {
    int p = rim.first;

    auto* sourceRemoteIndexList = rim.second.first;
    auto* targetRemoteIndexList = rim.second.second;

    // send to overlap
    for (auto const& ri : *sourceRemoteIndexList) {
      auto const& lip = ri.localIndexPair();
      Attribute remoteAttr = ri.attribute();
      Attribute myAttr = lip.local().attribute();
      if (myAttr == Attribute::owner && remoteAttr != Attribute::owner) {
        size_type globalIndex = globalIndices_[size_type(lip.local())];
        sendList[p].push_back({lip.global(), globalIndex});
      }
    }

    // receive from owner
    for (auto const& ri : *targetRemoteIndexList) {
      auto const& lip = ri.localIndexPair();
      Attribute remoteAttr = ri.attribute();
      Attribute myAttr = lip.local().attribute();
      if (myAttr != Attribute::owner && remoteAttr == Attribute::owner) {
        receiveList[p]++;
      }
    }
  }
  // all ghostDOFs must be communicated!
  assert(ghostSize_ == std::accumulate(receiveList.begin(), receiveList.end(), 0u));

  // send {globalId, globalIndex} to remote processors
  std::vector<Mpi::Request> sendRequests;
  for (int p = 0; p < world_.size(); ++p) {
    if (!sendList[p].empty()) {
      sendRequests.emplace_back( world_.isend(sendList[p], p, tag_) );
    }
  }

  // receive {globalID, globalIndex} from remote processors
  std::vector<Mpi::Request> recvRequests;
  std::vector<std::vector<GlobalAssoc>> recvData(world_.size());
  for (int p = 0; p < world_.size(); ++p) {
    if (receiveList[p] > 0)
      recvRequests.emplace_back( world_.irecv(recvData[p], p, tag_) );
  }

  Mpi::wait_all(recvRequests.begin(), recvRequests.end());

  ghostIndices_.reserve(ghostSize_);
  ghostLocalIndices_.resize(indexSet_->size(), LocalIndex(-1));

  // insert all remote global indices into the map
  std::size_t counter = 0;
  for (int p = 0; p < world_.size(); ++p) {
    auto const& data = recvData[p];
    assert(data.size() == receiveList[p]);
    for (auto const& d : data) {
      typename IndexSet::IndexPair const& l = indexSet_->at(d.first);
      assert(!owner_[size_type(l.local())]);

      globalIndices_[size_type(l.local())] = d.second;
      ghostIndices_.push_back(d.second);
      ghostLocalIndices_[size_type(l.local())] = counter++;
    }
  }
  assert(counter == ghostSize_);
  assert(ghostSize_ + localSize_ == indexSet_->size());

  Mpi::wait_all(sendRequests.begin(), sendRequests.end());
  msg("update IndexDistribution need {} sec", t.elapsed());
}


template <class GlobalId, class LI>
void PetscParallelIndexDistribution<GlobalId,LI>::
debug() const
{
  int p = world_.rank();
  std::cout << "[" << p << "]  sizes_.size()=" << sizes_.size() << ", my_size=" << sizes_[p] << std::endl;
  std::cout << "[" << p << "]  starts_.size()=" << starts_.size() << ", my_start=" << starts_[p] << std::endl;
  std::cout << "[" << p << "]  localSize_=" << localSize_ << ", globalSize_=" << globalSize_ << ", ghostSize_=" << ghostSize_ << std::endl;
  std::cout << "[" << p << "]  globalIndices_.size()=" << globalIndices_.size() << std::endl;
  std::cout << "[" << p << "]  ghostIndices_.size()=" << ghostIndices_.size() << std::endl;
  std::cout << "[" << p << "]  ghostLocalIndices_.size()=" << ghostLocalIndices_.size() << std::endl;
  std::cout << "[" << p << "]  owner_.size()=" << owner_.size() << std::endl;
}
#endif

} // end namespace AMDiS