torus.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (C) DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_GRID_YASPGRID_TORUS_HH
#define DUNE_GRID_YASPGRID_TORUS_HH
 
#include <array>
#include <bitset>
#include <cmath>
#include <deque>
#include <iostream>
#include <vector>
 
#if HAVE_MPI
#include <mpi.h>
#endif
 
#include <dune/common/binaryfunctions.hh>
#include <dune/common/streamoperators.hh>
#include <dune/grid/common/exceptions.hh>
 
#include "partitioning.hh"
 
namespace Dune
{
 
  template<class Communication, int d>
  class Torus {
  public:
    typedef std::array<int, d> iTupel;
 
 
  private:
    struct CommPartner {
      int rank;
      iTupel delta;
      int index;
    };
 
    struct CommTask {
      int rank;      // process to send to / receive from
      int size;      // size of buffer
      void *buffer;  // buffer to send / receive
    };
 
  public:
    Torus ()
    {}
 
    Torus (Communication comm, int tag, iTupel size, int overlap, const Yasp::Partitioning<d>* partitioner)
      : _comm(comm), _tag(tag)
    {
      // determine dimensions
      partitioner->partition(size, _comm.size(), _dims, overlap);
 
      // compute increments for lexicographic ordering
      int inc = 1;
      for (int i=0; i<d; i++)
      {
        _increment[i] = inc;
        inc *= _dims[i];
      }
 
      // check whether the load balancing matches the size of the communicator
      if (inc != _comm.size())
        DUNE_THROW(Dune::Exception, "Communicator size and result of the given load balancer do not match!");
 
      // make full schedule
      proclists();
    }
 
    int rank () const
    {
      return _comm.rank();
    }
 
    iTupel coord () const
    {
      return rank_to_coord(_comm.rank());
    }
 
    int procs () const
    {
      return _comm.size();
    }
 
    const iTupel & dims () const
    {
      return _dims;
    }
 
    int dims (int i) const
    {
      return _dims[i];
    }
 
    Communication comm () const
    {
      return _comm;
    }
 
    int tag () const
    {
      return _tag;
    }
 
    bool inside (iTupel c) const
    {
      for (int i=d-1; i>=0; i--)
        if (c[i]<0 || c[i]>=_dims[i]) return false;
      return true;
    }
 
    iTupel rank_to_coord (int rank) const
    {
      iTupel coord;
      rank = rank%_comm.size();
      for (int i=d-1; i>=0; i--)
      {
        coord[i] = rank/_increment[i];
        rank = rank%_increment[i];
      }
      return coord;
    }
 
    int coord_to_rank (iTupel coord) const
    {
      for (int i=0; i<d; i++) coord[i] = coord[i]%_dims[i];
      int rank = 0;
      for (int i=0; i<d; i++) rank += coord[i]*_increment[i];
      return rank;
    }
 
    int rank_relative (int rank, int dir, int cnt) const
    {
      iTupel coord = rank_to_coord(rank);
      coord[dir] = (coord[dir]+_dims[dir]+cnt)%_dims[dir];
      return coord_to_rank(coord);
    }
 
    int color (const iTupel & coord) const
    {
      int c = 0;
      int power = 1;
 
      // interior coloring
      for (int i=0; i<d; i++)
      {
        if (coord[i]%2==1) c += power;
        power *= 2;
      }
 
      // extra colors for boundary processes
      for (int i=0; i<d; i++)
      {
        if (_dims[i]>1 && coord[i]==_dims[i]-1) c += power;
        power *= 2;
      }
 
      return c;
    }
 
    int color (int rank) const
    {
      return color(rank_to_coord(rank));
    }
 
    int neighbors () const
    {
      int n=1;
      for (int i=0; i<d; ++i)
        n *= 3;
      return n-1;
    }
 
    bool is_neighbor (iTupel delta, std::bitset<d> periodic) const
    {
      iTupel coord = rank_to_coord(_comm.rank()); // my own coordinate with 0 <= c_i < dims_i
 
      for (int i=0; i<d; i++)
      {
        if (delta[i]<0)
        {
          // if I am on the boundary and domain is not periodic => no neighbor
          if (coord[i]==0 && periodic[i]==false) return false;
        }
        if (delta[i]>0)
        {
          // if I am on the boundary and domain is not periodic => no neighbor
          if (coord[i]==_dims[i]-1 && periodic[i]==false) return false;
        }
      }
      return true;
    }
 
    double partition (int rank, iTupel origin_in, iTupel size_in, iTupel& origin_out, iTupel& size_out) const
    {
      iTupel coord = rank_to_coord(rank);
      double maxsize = 1;
      double sz = 1;
 
      // make a tensor product partition
      for (int i=0; i<d; i++)
      {
        // determine
        int m = size_in[i]/_dims[i];
        int r = size_in[i]%_dims[i];
 
        sz *= size_in[i];
 
        if (coord[i]<_dims[i]-r)
        {
          origin_out[i] = origin_in[i] + coord[i]*m;
          size_out[i] = m;
          maxsize *= m;
        }
        else
        {
          origin_out[i] = origin_in[i] + (_dims[i]-r)*m + (coord[i]-(_dims[i]-r))*(m+1);
          size_out[i] = m+1;
          maxsize *= m+1;
        }
      }
      return maxsize/(sz/_comm.size());
    }
 
    class ProcListIterator {
    public:
      ProcListIterator (typename std::deque<CommPartner>::const_iterator iter)
      {
        i = iter;
      }
 
      int rank () const
      {
        return i->rank;
      }
 
      iTupel delta () const
      {
        return i->delta;
      }
 
      int index () const
      {
        return i->index;
      }
 
      int distance () const
      {
        int dist = 0;
        iTupel delta=i->delta;
        for (int j=0; j<d; ++j)
          dist += std::abs(delta[j]);
        return dist;
      }
 
      bool operator== (const ProcListIterator& iter) const
      {
        return i == iter.i;
      }
 
 
      bool operator!= (const ProcListIterator& iter) const
      {
        return i != iter.i;
      }
 
      ProcListIterator& operator++ ()
      {
        ++i;
        return *this;
      }
 
    private:
      typename std::deque<CommPartner>::const_iterator i;
    };
 
    ProcListIterator sendbegin () const
    {
      return ProcListIterator(_sendlist.begin());
    }
 
    ProcListIterator sendend () const
    {
      return ProcListIterator(_sendlist.end());
    }
 
    ProcListIterator recvbegin () const
    {
      return ProcListIterator(_recvlist.begin());
    }
 
    ProcListIterator recvend () const
    {
      return ProcListIterator(_recvlist.end());
    }
 
    void send (int rank, void* buffer, int size) const
    {
      CommTask task;
      task.rank = rank;
      task.buffer = buffer;
      task.size = size;
      if (rank!=_comm.rank())
        _sendrequests.push_back(task);
      else
        _localsendrequests.push_back(task);
    }
 
    void recv (int rank, void* buffer, int size) const
    {
      CommTask task;
      task.rank = rank;
      task.buffer = buffer;
      task.size = size;
      if (rank!=_comm.rank())
        _recvrequests.push_back(task);
      else
        _localrecvrequests.push_back(task);
    }
 
    void exchange () const
    {
      // handle local requests first
      if (_localsendrequests.size()!=_localrecvrequests.size())
      {
        std::cout << "[" << rank() << "]: ERROR: local sends/receives do not match in exchange!" << std::endl;
        return;
      }
      for (unsigned int i=0; i<_localsendrequests.size(); i++)
      {
        if (_localsendrequests[i].size!=_localrecvrequests[i].size)
        {
          std::cout << "[" << rank() << "]: ERROR: size in local sends/receive does not match in exchange!" << std::endl;
          return;
        }
        memcpy(_localrecvrequests[i].buffer,_localsendrequests[i].buffer,_localsendrequests[i].size);
      }
      _localsendrequests.clear();
      _localrecvrequests.clear();
 
#if HAVE_MPI
      // handle foreign requests
 
      std::vector<MPI_Request> requests(_sendrequests.size() + _recvrequests.size());
      MPI_Request* req = requests.data();
 
      // issue sends to foreign processes
      for (unsigned int i=0; i<_sendrequests.size(); i++)
        if (_sendrequests[i].rank!=rank())
        {
          //          std::cout << "[" << rank() << "]" << " send " << _sendrequests[i].size << " bytes "
          //                    << "to " << _sendrequests[i].rank << " p=" << _sendrequests[i].buffer << std::endl;
          MPI_Isend(_sendrequests[i].buffer, _sendrequests[i].size, MPI_BYTE,
                    _sendrequests[i].rank, _tag, _comm, req++);
        }
 
      // issue receives from foreign processes
      for (unsigned int i=0; i<_recvrequests.size(); i++)
        if (_recvrequests[i].rank!=rank())
        {
          //          std::cout << "[" << rank() << "]"  << " recv " << _recvrequests[i].size << " bytes "
          //                    << "fm " << _recvrequests[i].rank << " p=" << _recvrequests[i].buffer << std::endl;
          MPI_Irecv(_recvrequests[i].buffer, _recvrequests[i].size, MPI_BYTE,
                    _recvrequests[i].rank, _tag, _comm, req++);
        }
 
      // Wait for communication to complete
      MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
 
      // clear request buffers
      _sendrequests.clear();
      _recvrequests.clear();
#endif
    }
 
    double global_max (double x) const
    {
      double res = 0.0;
      _comm.template allreduce<Dune::Max<double>,double>(&x, &res, 1);
      return res;
    }
 
    void print (std::ostream& s) const
    {
      s << "[" << rank() <<  "]: Torus " << procs() << " processor(s) arranged as " << dims() << std::endl;
      for (ProcListIterator i=sendbegin(); i!=sendend(); ++i)
      {
        s << "[" << rank() <<  "]: send to   "
          << "rank=" << i.rank()
          << " index=" << i.index()
          << " delta=" << i.delta() << " dist=" << i.distance() << std::endl;
      }
      for (ProcListIterator i=recvbegin(); i!=recvend(); ++i)
      {
        s << "[" << rank() <<  "]: recv from "
          << "rank=" << i.rank()
          << " index=" << i.index()
          << " delta=" << i.delta() << " dist=" << i.distance() << std::endl;
      }
    }
 
  private:
 
    void proclists ()
    {
      // compile the full neighbor list
      CommPartner cp;
      iTupel delta;
 
      std::fill(delta.begin(), delta.end(), -1);
      bool ready = false;
      iTupel me, nb;
      me = rank_to_coord(_comm.rank());
      int index = 0;
      int last = neighbors()-1;
      while (!ready)
      {
        // find neighbors coordinates
        for (int i=0; i<d; i++)
          nb[i] = ( me[i]+_dims[i]+delta[i] ) % _dims[i];
 
        // find neighbors rank
        int nbrank = coord_to_rank(nb);
 
        // check if delta is not zero
        for (int i=0; i<d; i++)
          if (delta[i]!=0)
          {
            cp.rank = nbrank;
            cp.delta = delta;
            cp.index = index;
            _recvlist.push_back(cp);
            cp.index = last-index;
            _sendlist.push_front(cp);
            index++;
            break;
          }
 
        // next neighbor
        ready = true;
        for (int i=0; i<d; i++)
          if (delta[i]<1)
          {
            (delta[i])++;
            ready=false;
            break;
          }
          else
          {
            delta[i] = -1;
          }
      }
 
    }
 
    Communication _comm;
 
    iTupel _dims;
    iTupel _increment;
    int _tag;
    std::deque<CommPartner> _sendlist;
    std::deque<CommPartner> _recvlist;
 
    mutable std::vector<CommTask> _sendrequests;
    mutable std::vector<CommTask> _recvrequests;
    mutable std::vector<CommTask> _localsendrequests;
    mutable std::vector<CommTask> _localrecvrequests;
 
  };
 
  template <class Communication, int d>
  inline std::ostream& operator<< (std::ostream& s, const Torus<Communication, d> & t)
  {
    t.print(s);
    return s;
  }
}
 
#endif