Domain.h

Go to the documentation of this file.

/***************************************************************************
 *  ThunderEgg, a library for solvers on adaptively refined block-structured
 *  Cartesian grids.
 *
 *  Copyright (c) 2019-2021 Scott Aiton
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 ***************************************************************************/
 
#ifndef THUNDEREGG_DOMAIN_H
#define THUNDEREGG_DOMAIN_H
 
#include <ThunderEgg/PatchInfo.h>
#include <ThunderEgg/Timer.h>
#include <map>
#include <set>
#include <vector>
 
namespace ThunderEgg {
template<int D>
class Domain
{
private:
  Communicator comm;
  int id = -1;
  std::array<int, D> ns;
  int num_ghost_cells;
  int num_cells_in_patch;
  int num_cells_in_patch_with_ghost;
  std::vector<PatchInfo<D>> pinfos;
  int global_num_patches = 1;
  mutable std::shared_ptr<Timer> timer;
 
  void indexPatchesLocal()
  {
    // index patches
    int curr_index = 0;
    std::map<int, int> id_to_local_index;
    for (auto& pinfo : pinfos) {
      pinfo.local_index = curr_index;
      id_to_local_index[pinfo.id] = pinfo.local_index;
      curr_index++;
    }
 
    // set local index in nbrinfo objects
    for (auto& pinfo : pinfos) {
      pinfo.setNeighborLocalIndexes(id_to_local_index);
    }
  }
  void indexPatchesGlobal()
  {
    // get starting global index
    int num_local_patches = (int)pinfos.size();
    int curr_global_index;
    MPI_Scan(&num_local_patches, &curr_global_index, 1, MPI_INT, MPI_SUM, comm.getMPIComm());
    curr_global_index -= num_local_patches;
 
    // index the patches
    std::map<int, int> id_to_global_index;
    for (auto& pinfo : pinfos) {
      pinfo.global_index = curr_global_index;
      id_to_global_index[pinfo.id] = pinfo.global_index;
      curr_global_index++;
    }
 
    std::map<int, std::set<std::pair<int, int>>> ranks_to_ids_and_global_indexes_outgoing;
    std::map<int, std::set<int>> ranks_to_ids_incoming;
    for (auto& pinfo : pinfos) {
      auto ranks = pinfo.getNbrRanks();
      auto ids = pinfo.getNbrIds();
      for (size_t idx = 0; idx < ranks.size(); idx++) {
        int nbr_id = ids[idx];
        int nbr_rank = ranks[idx];
        if (nbr_rank != comm.getRank()) {
          ranks_to_ids_and_global_indexes_outgoing[nbr_rank].insert(
            std::make_pair(pinfo.id, pinfo.global_index));
          ranks_to_ids_incoming[nbr_rank].insert(nbr_id);
        }
      }
    }
 
    // prepare to recieve data
    std::vector<MPI_Request> recv_requests;
 
    // allocate incoming vectors and post recvs
    std::map<int, std::vector<int>> rank_to_incoming_data;
    for (const auto& pair : ranks_to_ids_incoming) {
      int source_rank = pair.first;
      std::vector<int>& incoming_data = rank_to_incoming_data[source_rank];
      incoming_data.resize(pair.second.size());
 
      MPI_Request request;
      MPI_Irecv(incoming_data.data(),
                (int)incoming_data.size(),
                MPI_INT,
                source_rank,
                0,
                comm.getMPIComm(),
                &request);
      recv_requests.push_back(request);
    }
 
    // prepare outgoing vector of data and send it
    std::vector<MPI_Request> send_requests;
 
    // post sends
    std::map<int, std::vector<int>> rank_to_outgoing_data;
    for (const auto& pair : ranks_to_ids_and_global_indexes_outgoing) {
      int dest_rank = pair.first;
      // allocate and fill vector
      std::vector<int>& data = rank_to_outgoing_data[dest_rank];
      data.reserve(pair.second.size());
      for (const auto& id_and_global_index : pair.second) {
        data.push_back(id_and_global_index.second);
      }
      MPI_Request request;
      MPI_Isend(data.data(), (int)data.size(), MPI_INT, dest_rank, 0, comm.getMPIComm(), &request);
      send_requests.push_back(request);
    }
 
    // add global indexes to map as recvs come in
    for (size_t i = 0; i < recv_requests.size(); i++) {
      MPI_Status status;
      int request_index;
      MPI_Waitany((int)recv_requests.size(), recv_requests.data(), &request_index, &status);
 
      int source_rank = status.MPI_SOURCE;
      const std::set<int>& incoming_ids = ranks_to_ids_incoming[source_rank];
      const std::vector<int>& data = rank_to_incoming_data[source_rank];
 
      auto curr_id = incoming_ids.cbegin();
      auto curr_global_index_iter = data.cbegin();
      while (curr_id != incoming_ids.cend()) {
        id_to_global_index[*curr_id] = *curr_global_index_iter;
        curr_id++;
        curr_global_index_iter++;
      }
    }
 
    // wait for all the sends to finsh
    MPI_Waitall((int)send_requests.size(), send_requests.data(), MPI_STATUSES_IGNORE);
 
    // update global indexes in nbrinfo objects
    for (auto& pinfo : pinfos) {
      pinfo.setNeighborGlobalIndexes(id_to_global_index);
    }
  }
 
public:
  template<class InputIterator>
  Domain(Communicator comm,
         int id,
         std::array<int, D> ns,
         int num_ghost_cells,
         InputIterator first_pinfo,
         InputIterator last_pinfo)
    : comm(comm)
    , id(id)
    , ns(ns)
    , num_ghost_cells(num_ghost_cells)
    , pinfos(first_pinfo, last_pinfo)
  {
    num_cells_in_patch = 1;
    num_cells_in_patch_with_ghost = 1;
    for (size_t i = 0; i < D; i++) {
      num_cells_in_patch *= ns[i];
      num_cells_in_patch_with_ghost *= (ns[i] + 2 * num_ghost_cells);
    }
 
    int num_local_domains = pinfos.size();
    MPI_Allreduce(&num_local_domains, &global_num_patches, 1, MPI_INT, MPI_SUM, comm.getMPIComm());
 
    indexPatchesLocal();
    indexPatchesGlobal();
  }
  const Communicator& getCommunicator() const { return comm; }
  const std::vector<PatchInfo<D>>& getPatchInfoVector() const { return pinfos; }
  const std::array<int, D>& getNs() const { return ns; }
  int getNumGlobalPatches() const { return global_num_patches; }
  int getNumLocalPatches() const { return (int)pinfos.size(); }
  int getNumGlobalCells() const { return global_num_patches * num_cells_in_patch; }
  int getNumLocalCells() const { return ((int)pinfos.size()) * num_cells_in_patch; }
  int getNumLocalCellsWithGhost() const
  {
    return ((int)pinfos.size()) * num_cells_in_patch_with_ghost;
  }
  int getNumCellsInPatch() const { return num_cells_in_patch; }
  int getNumGhostCells() const { return num_ghost_cells; }
  double volume() const
  {
    double sum = 0;
    for (auto& pinfo : pinfos) {
      double patch_vol = 1;
      for (size_t i = 0; i < D; i++) {
        patch_vol *= pinfo.spacings[i] * pinfo.ns[i];
      }
      sum += patch_vol;
    }
    double retval;
    MPI_Allreduce(&sum, &retval, 1, MPI_DOUBLE, MPI_SUM, comm.getMPIComm());
    return retval;
  }
  void setTimer(std::shared_ptr<Timer> timer) const
  {
    this->timer = timer;
    timer->addDomain(id, *this);
  }
  std::shared_ptr<Timer> getTimer() const { return timer; }
  bool hasTimer() const { return timer != nullptr; }
  int getId() const { return id; }
};
 
template<int D>
void
to_json(tpl::nlohmann::json& j, const Domain<D>& domain)
{
  for (auto pinfo : domain.getPatchInfoVector()) {
    j.push_back(pinfo);
  }
}
 
extern template class Domain<2>;
extern template class Domain<3>;
extern template void
to_json<2>(tpl::nlohmann::json& j, const Domain<2>& domain);
extern template void
to_json<3>(tpl::nlohmann::json& j, const Domain<3>& domain);
} // namespace ThunderEgg
#endif