File mpicartesiangroup.h

File List > cartesian > mpicartesiangroup.h

Go to the documentation of this file

#ifndef TEMPLAT_PARALLEL_MPI_COMM_MPICARTESIANGROUP_H
#define TEMPLAT_PARALLEL_MPI_COMM_MPICARTESIANGROUP_H

/* This file is part of TempLat, available at https://cosmolattice.github.io/templat .
   Copyright 2021-2026 The TempLat authors, see AUTHORS.md.
   Released under the MIT license, see LICENSE.md. */

// File info: Main contributor(s): Wessel Valkenburg, Year: 2019

#include "TempLat/util/exception.h"

#include "TempLat/parallel/mpi/comm/mpicommreference.h"
#include "TempLat/parallel/mpi/comm/mpidomainsplit.h"

namespace TempLat
{

  MakeException(MPICartesianGroupException);

  class MPICartesianGroup
  {
  public:
    // Put public methods here. These should change very little over time.
    MPICartesianGroup(MPICommReference topologyGroup, ptrdiff_t nDimensions, std::vector<int> decomposition,
                      std::vector<int> expectedCoords = std::vector<int>(),
                      std::vector<int> periodic = std::vector<int>())
        : mBaseGroup(topologyGroup), mCartesianGroup(mBaseGroup), mNDimensions(nDimensions),
          mDecomposition(decomposition), mPeriodic(periodic), mExpectedCoords(expectedCoords)
    {
      verifyInput();

      createGroups();

      verifyCoordsMatchTopology();
    }

    MPICartesianGroup(ptrdiff_t nDimensions, std::vector<int> decomposition)
        : MPICartesianGroup(MPICommReference(), nDimensions, decomposition)
    {
    }

    MPI_Comm getComm() const { return mCartesianGroup; }
    MPICommReference getBaseComm() const { return mBaseGroup; }

    int getRank() { return mCartesianGroup.getRank(); }

    ptrdiff_t getNDimensions() const { return mNDimensions; }

    const std::vector<int> &getPosition() const { return mSelfPosition; }

    const std::vector<int> &getDecomposition() const { return mDecomposition; }

    auto size() const { return mBaseGroup.size(); }

    friend std::ostream &operator<<(std::ostream &ostream, MPICartesianGroup mGr)
    {
      ostream << "MPICartesianGroup:\n  decomposition: " << mGr.getDecomposition() << "\n  this rank: " << mGr.getRank()
              << "\n  this position: " << mGr.getPosition() << "\n";
      return ostream;
    }

  private:
    /* Put all member variables and private methods here. These may change arbitrarily. */

    MPICommReference mBaseGroup;
    MPICommReference mCartesianGroup;

    const ptrdiff_t mNDimensions;
    std::vector<int> mDecomposition;
    std::vector<int> mPeriodic;
    std::vector<int> mExpectedCoords;
    std::vector<int> mSelfPosition;

    std::vector<int> fetchPosition(int ofRank)
    {
      std::vector<int> result(mDecomposition.size(), 0);
#ifdef HAVE_MPI
      // sayMPI << this << "Getting coords: " << ofRank << ", " << result << " nd: " <<mNDimensions<< "\n";
      if (0 != MPI_Cart_coords(mCartesianGroup, ofRank, (int)result.size(), result.data())) {
        throw MPICartesianGroupException("Could not get result from MPI_Cart_coords.");
      }
#endif
      return result;
    }

    void createGroups()
    {
#ifdef HAVE_MPI
      /* No reordering. The rank order of mBaseGroup is the topology, so the Cartesian
       * communicator must preserve it: the layout's local starts are derived from that order,
       * and ghost neighbours are derived from MPI_Cart_shift on the communicator built here.
       * Letting MPI permute ranks lets those two disagree. */
      mCartesianGroup = createOneGroup(mBaseGroup, mNDimensions, false);
#endif

      mSelfPosition = fetchPosition(mCartesianGroup.rank());
    }

    MPICommReference createOneGroup(MPICommReference fromGroup, ptrdiff_t nDim, bool reorder)
    {
      MPI_Comm newComm;

#ifdef HAVE_MPI
      /* Build from fromGroup, not mBaseGroup: this used to ignore its own argument. */
      MPI_Cart_create(fromGroup, nDim, mDecomposition.data(), mPeriodic.data(), reorder, &newComm);
#else
      (void)fromGroup;
      (void)nDim;
      (void)reorder;
      newComm = MPI_COMM_WORLD;
#endif
      return MPICommReference(newComm);
    }

    void verifyCoordsMatchTopology() const
    {
      if (mExpectedCoords.empty()) return;

      if ((ptrdiff_t)mExpectedCoords.size() != mNDimensions)
        throw MPICartesianGroupException("Expected coordinates have ", mExpectedCoords.size(), " entries, but this is a ",
                                         mNDimensions, "-dimensional group.");

      for (ptrdiff_t i = 0; i < mNDimensions; ++i) {
        if (mSelfPosition[i] != mExpectedCoords[i])
          throw MPICartesianGroupException(
              "Cartesian coordinates disagree with the topology that produced this group: MPI places this rank at ",
              mSelfPosition, " but the topology's owner reports ", mExpectedCoords, " (decomposition ", mDecomposition,
              "). The layout's local starts follow the latter while ghost exchange follows the former, so continuing "
              "would silently exchange the wrong data at subdomain boundaries.");
      }
    }

    void verifyInput()
    {
      while ((ptrdiff_t)mDecomposition.size() < mNDimensions)
        mDecomposition.push_back(1);
      while (mPeriodic.size() < mDecomposition.size())
        mPeriodic.push_back(1);

      if ((ptrdiff_t)mDecomposition.size() > mNDimensions) mDecomposition.resize(mNDimensions);
      if ((ptrdiff_t)mPeriodic.size() > mNDimensions) mPeriodic.resize(mNDimensions);

      /* The decomposition must tile the communicator exactly. This used to "repair" a mismatch
       * by overwriting mDecomposition[0] and then zeroing trailing entries from the right, which
       * silently mangled an explicitly-pinned backend grid into something the backend would not
       * use — precisely the disagreement this class now exists to prevent. */
      ptrdiff_t groupSize = mBaseGroup.size();
      ptrdiff_t product = 1;
      for (auto &&it : mDecomposition)
        product *= it;

      if (product != groupSize)
        throw MPICartesianGroupException("Decomposition ", mDecomposition, " describes ", product,
                                         " ranks but the communicator has ", groupSize,
                                         ". The decomposition must tile the communicator exactly; build the group via "
                                         "FFTMPIDomainSplit::makeMPIGroup so it matches the FFT backend.");
    }
  };
} // namespace TempLat

#endif