MPIGhostFiller.h

Go to the documentation of this file.

/***************************************************************************
 *  ThunderEgg, a library for solvers on adaptively refined block-structured
 *  Cartesian grids.
 *
 *  Copyright (c) 2020-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_MPIGHOSTFILLER_H
#define THUNDEREGG_MPIGHOSTFILLER_H
 
#include <ThunderEgg/DimensionalArray.h>
#include <ThunderEgg/Domain.h>
#include <ThunderEgg/GhostFiller.h>
#include <ThunderEgg/GhostFillingType.h>
 
namespace ThunderEgg {
template<int D>
class MPIGhostFiller : public GhostFiller<D>
{
private:
  template<int M>
  class RemoteCallPrototype
  {
  public:
    int id;
    Face<D, M> face;
    NbrType nbr_type;
    Orthant<M> orthant;
    int local_index;
    RemoteCallPrototype(int id,
                        Face<D, M> face,
                        NbrType nbr_type,
                        Orthant<M> orthant,
                        int local_index)
      : id(id)
      , face(face)
      , nbr_type(nbr_type)
      , orthant(orthant)
      , local_index(local_index)
    {}
    bool operator<(const RemoteCallPrototype& other) const
    {
      return std::forward_as_tuple(id, face.opposite(), nbr_type, orthant, local_index) <
             std::forward_as_tuple(
               other.id, other.face.opposite(), other.nbr_type, other.orthant, other.local_index);
    }
  };
 
  template<int M>
  class IncomingGhostPrototype
  {
  public:
    int id;
    Face<D, M> face;
    int local_index;
    IncomingGhostPrototype(int id, Face<D, M> face, int local_index)
      : id(id)
      , face(face)
      , local_index(local_index)
    {}
    bool operator<(const IncomingGhostPrototype& other) const
    {
      return std::forward_as_tuple(id, face, local_index) <
             std::forward_as_tuple(other.id, other.face, other.local_index);
    }
  };
 
  template<int M>
  struct RemoteCall
  {
    Face<D, M> face;
    NbrType nbr_type;
    Orthant<M> orthant;
    int local_index;
    size_t offset;
    RemoteCall(const RemoteCallPrototype<M>& prototype, size_t offset)
      : face(prototype.face)
      , nbr_type(prototype.nbr_type)
      , orthant(prototype.orthant)
      , local_index(prototype.local_index)
      , offset(offset)
    {}
  };
 
  template<int M>
  struct IncomingGhost
  {
    Face<D, M> face;
    int local_index;
    size_t offset;
    IncomingGhost(const IncomingGhostPrototype<M>& prototype, size_t offset)
      : face(prototype.face)
      , local_index(prototype.local_index)
      , offset(offset)
    {}
  };
 
  template<int M>
  struct LocalCall
  {
    Face<D, M> face;
    NbrType nbr_type;
    Orthant<M> orthant;
    int local_index;
    int nbr_local_index;
    LocalCall(Face<D, M> face,
              NbrType nbr_type,
              Orthant<M> orthant,
              int local_index,
              size_t nbr_local_index)
      : face(face)
      , nbr_type(nbr_type)
      , orthant(orthant)
      , local_index(local_index)
      , nbr_local_index(nbr_local_index)
    {}
  };
 
  struct RemoteCallSet
  {
    int rank;
    size_t send_buffer_length = 0;
    template<int M>
    using RemoteCallDeque = std::deque<RemoteCall<M>>;
    DimensionalArray<D, RemoteCallDeque> remote_calls;
    size_t recv_buffer_length = 0;
    template<int M>
    using IncomingGhostDeque = std::deque<IncomingGhost<M>>;
    DimensionalArray<D, IncomingGhostDeque> incoming_ghosts;
 
    explicit RemoteCallSet(int rank)
      : rank(rank)
    {}
  };
 
  std::vector<RemoteCallSet> remote_call_sets;
 
  template<int M>
  using LocalCallDeque = std::deque<LocalCall<M>>;
  DimensionalArray<D, LocalCallDeque> local_calls;
 
  template<int M>
  class GhostViewInfo
  {
  private:
    std::array<std::array<int, D + 1>, Face<D, M>::number_of> strides;
    std::array<std::array<int, D + 1>, Face<D, M>::number_of> ghost_start;
    std::array<std::array<int, D + 1>, Face<D, M>::number_of> start;
    std::array<std::array<int, D + 1>, Face<D, M>::number_of> end;
    std::array<std::array<int, D + 1>, Face<D, M>::number_of> ghost_end;
    std::array<size_t, Face<D, M>::number_of> sizes;
 
  public:
    GhostViewInfo() = default;
    explicit GhostViewInfo(const Domain<D>& domain)
    {
      std::array<int, D + 1> ns;
      for (int i = 0; i < D; i++) {
        ns[i] = domain.getNs()[i];
      }
      ns[D] = 1;
      int num_ghost_cells = domain.getNumGhostCells();
      for (Face<D, M> face : Face<D, M>::getValues()) {
        std::array<int, D + 1> ghost_ns = ns;
 
        std::array<int, D + 1> face_ghost_start;
        face_ghost_start.fill(0);
 
        std::array<int, D + 1> face_start;
        face_start.fill(0);
 
        std::array<int, D + 1> face_end;
        face_end = ns;
        for (int& v : face_end) {
          v--;
        }
 
        std::array<int, D + 1> face_ghost_end;
        face_ghost_end = ns;
        for (int& v : face_ghost_end) {
          v--;
        }
 
        for (Side<D> side : face.getSides()) {
          ghost_ns[side.getAxisIndex()] = num_ghost_cells;
          if (side.isLowerOnAxis()) {
            face_ghost_start[side.getAxisIndex()] = -num_ghost_cells;
            face_ghost_end[side.getAxisIndex()] = -1;
            face_end[side.getAxisIndex()] = -1;
          } else {
            face_ghost_start[side.getAxisIndex()] = ns[side.getAxisIndex()];
            face_start[side.getAxisIndex()] = ns[side.getAxisIndex()];
            face_ghost_end[side.getAxisIndex()] = ns[side.getAxisIndex()] + num_ghost_cells - 1;
          }
        }
        strides[face.getIndex()][0] = 1;
        for (size_t i = 1; i < D + 1; i++) {
          strides[face.getIndex()][i] = ghost_ns[i - 1] * strides[face.getIndex()][i - 1];
        }
        sizes[face.getIndex()] = ghost_ns[D - 1] * strides[face.getIndex()][D - 1];
        int offset = 0;
        for (Side<D> side : face.getSides()) {
          if (side.isHigherOnAxis()) {
            offset += -(ns[side.getAxisIndex()] + num_ghost_cells) *
                      strides[face.getIndex()][side.getAxisIndex()];
          }
        }
        ghost_start[face.getIndex()] = face_ghost_start;
        start[face.getIndex()] = face_start;
        end[face.getIndex()] = face_end;
        ghost_end[face.getIndex()] = face_ghost_end;
      }
    }
 
    size_t getSize(Face<D, M> face) const { return sizes[face.getIndex()]; }
 
    PatchView<const double, D> getPatchView(double* buffer_ptr,
                                            Face<D, M> face,
                                            int num_components) const
    {
      std::array<int, D + 1> my_end = end[face.getIndex()];
      std::array<int, D + 1> my_ghost_end = ghost_end[face.getIndex()];
      my_end[D] = num_components - 1;
      my_ghost_end[D] = num_components - 1;
      return PatchView<const double, D>(buffer_ptr,
                                        strides[face.getIndex()],
                                        ghost_start[face.getIndex()],
                                        start[face.getIndex()],
                                        my_end,
                                        my_ghost_end);
    }
  };
 
  DimensionalArray<D, GhostViewInfo> ghost_local_data_infos;
 
  Domain<D> domain;
 
  GhostFillingType fill_type;
 
  template<int M>
  PatchView<const double, D> getPatchViewForBuffer(double* buffer_ptr,
                                                   Face<D, M> face,
                                                   int num_components) const
  {
    const GhostViewInfo<M>& gld_info = ghost_local_data_infos.template get<M>();
    return gld_info.getPatchView(buffer_ptr, face, num_components);
  }
  std::vector<MPI_Request> postRecvs(std::vector<std::vector<double>>& buffers) const
  {
    std::vector<MPI_Request> recv_requests(buffers.size());
    for (size_t i = 0; i < recv_requests.size(); i++) {
      MPI_Irecv(buffers[i].data(),
                buffers[i].size(),
                MPI_DOUBLE,
                remote_call_sets[i].rank,
                0,
                MPI_COMM_WORLD,
                &recv_requests[i]);
    }
    return recv_requests;
  }
 
  template<int M>
  void addRecvBufferToGhost(const RemoteCallSet& remote_call_set,
                            std::vector<double>& buffer,
                            const Vector<D>& u) const
  {
    for (const IncomingGhost<M>& incoming_ghost :
         remote_call_set.incoming_ghosts.template get<M>()) {
      int local_index = incoming_ghost.local_index;
      Face<D, M> face = incoming_ghost.face;
      size_t buffer_offset = incoming_ghost.offset;
 
      PatchView<const double, D> local_view = u.getPatchView(local_index);
      double* buffer_ptr = buffer.data() + buffer_offset * u.getNumComponents();
      PatchView<const double, D> buffer_view =
        getPatchViewForBuffer(buffer_ptr, face, u.getNumComponents());
      std::array<size_t, D - M> start;
      start.fill(0);
      std::array<size_t, D - M> end;
      end.fill(domain.getNumGhostCells() - 1);
      Loop::Nested<D - M>(start, end, [&](const std::array<size_t, D - M>& offset) {
        View<double, M + 1> local_slice = local_view.getGhostSliceOn(face, offset);
        View<double, M + 1> buffer_slice = buffer_view.getGhostSliceOn(face, offset);
        Loop::OverInteriorIndexes<M + 1>(local_slice, [&](const std::array<int, M + 1>& coord) {
          local_slice[coord] += buffer_slice[coord];
        });
      });
    }
  }
  void processRecvs(std::vector<MPI_Request>& requests,
                    std::vector<std::vector<double>>& buffers,
                    const Vector<D>& u) const
  {
    size_t num_requests = requests.size();
    for (size_t i = 0; i < num_requests; i++) {
      int finished_index;
      MPI_Waitany(requests.size(), requests.data(), &finished_index, MPI_STATUS_IGNORE);
      switch (fill_type) {
        case GhostFillingType::Corners:
          if constexpr (D >= 2) {
            addRecvBufferToGhost<0>(remote_call_sets[finished_index], buffers[finished_index], u);
          }
        case GhostFillingType::Edges:
          if constexpr (D == 3) {
            addRecvBufferToGhost<1>(remote_call_sets[finished_index], buffers[finished_index], u);
          }
        case GhostFillingType::Faces:
          addRecvBufferToGhost<D - 1>(remote_call_sets[finished_index], buffers[finished_index], u);
          break;
        default:
          throw RuntimeError("Unsupported GhostFilling Type");
      }
    }
  }
 
  template<int M>
  void fillSendBuffer(const RemoteCallSet& remote_call_set,
                      std::vector<double>& buffer,
                      const Vector<D>& u) const
  {
    for (const RemoteCall<M>& call : remote_call_set.remote_calls.template get<M>()) {
      const PatchInfo<D>& pinfo = domain.getPatchInfoVector()[call.local_index];
      Face<D, M> face = call.face.opposite();
      PatchView<const double, D> local_view = u.getPatchView(call.local_index);
      double* buffer_ptr = buffer.data() + call.offset * u.getNumComponents();
 
      // create View objects for the buffer
      PatchView<const double, D> buffer_view =
        getPatchViewForBuffer(buffer_ptr, face, u.getNumComponents());
 
      // make the call
      fillGhostCellsForNbrPatchPriv(
        pinfo, local_view, buffer_view, call.face, call.nbr_type, call.orthant);
    }
  }
 
  std::vector<MPI_Request> postSends(std::vector<std::vector<double>>& buffers,
                                     const Vector<D>& u) const
  {
    std::vector<MPI_Request> send_requests(remote_call_sets.size());
    for (size_t i = 0; i < remote_call_sets.size(); i++) {
      switch (fill_type) {
        case GhostFillingType::Corners:
          if constexpr (D >= 2) {
            fillSendBuffer<0>(remote_call_sets[i], buffers[i], u);
          }
        case GhostFillingType::Edges:
          if constexpr (D == 3) {
            fillSendBuffer<1>(remote_call_sets[i], buffers[i], u);
          }
        case GhostFillingType::Faces:
          fillSendBuffer<D - 1>(remote_call_sets[i], buffers[i], u);
          break;
        default:
          throw RuntimeError("Unsupported GhostFilling Type");
      }
      MPI_Isend(buffers[i].data(),
                buffers[i].size(),
                MPI_DOUBLE,
                remote_call_sets[i].rank,
                0,
                MPI_COMM_WORLD,
                &send_requests[i]);
    }
    return send_requests;
  }
 
  template<int M>
  void fillGhostCellsForNbrPatchPriv(const PatchInfo<D>& pinfo,
                                     const PatchView<const double, D>& local_view,
                                     const PatchView<const double, D>& nbr_view,
                                     Face<D, M> face,
                                     NbrType nbr_type,
                                     Orthant<M> orthant_on_coarse) const
  {
    if constexpr (M == D - 1) {
      fillGhostCellsForNbrPatch(pinfo, local_view, nbr_view, face, nbr_type, orthant_on_coarse);
    } else if constexpr (M > 0) {
      fillGhostCellsForEdgeNbrPatch(pinfo, local_view, nbr_view, face, nbr_type, orthant_on_coarse);
    } else {
      fillGhostCellsForCornerNbrPatch(pinfo, local_view, nbr_view, face, nbr_type);
    }
  }
 
  template<int M>
  void processLocalFills(const Vector<D>& u) const
  {
    for (const LocalCall<M>& call : local_calls.template get<M>()) {
      const PatchInfo<D>& pinfo = domain.getPatchInfoVector()[call.local_index];
      PatchView<const double, D> local_view = u.getPatchView(call.local_index);
      PatchView<const double, D> nbr_view = u.getPatchView(call.nbr_local_index);
      fillGhostCellsForNbrPatchPriv(
        pinfo, local_view, nbr_view, call.face, call.nbr_type, call.orthant);
    }
    if constexpr (M > 0) {
      processLocalFills<M - 1>(u);
    }
  }
 
  template<int M>
  static void addNormalNbrCalls(
    const PatchInfo<D>& pinfo,
    Face<D, M> f,
    std::deque<LocalCall<M>>& my_local_calls,
    std::map<int, std::set<RemoteCallPrototype<M>>>& rank_to_remote_call_prototypes,
    std::map<int, std::set<IncomingGhostPrototype<M>>>& rank_to_incoming_ghost_prototypes)
  {
    int rank;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    auto nbrinfo = pinfo.getNormalNbrInfo(f);
    if (nbrinfo.rank == rank) {
      my_local_calls.emplace_back(
        f, NbrType::Normal, Orthant<M>::null(), pinfo.local_index, nbrinfo.local_index);
    } else {
      rank_to_remote_call_prototypes[nbrinfo.rank].emplace(
        nbrinfo.id, f, NbrType::Normal, Orthant<M>::null(), pinfo.local_index);
      rank_to_incoming_ghost_prototypes[nbrinfo.rank].emplace(pinfo.id, f, pinfo.local_index);
    }
  }
 
  template<int M>
  static void addFineNbrCalls(
    const PatchInfo<D>& pinfo,
    Face<D, M> f,
    std::deque<LocalCall<M>>& my_local_calls,
    std::map<int, std::set<RemoteCallPrototype<M>>>& rank_to_remote_call_prototypes,
    std::map<int, std::set<IncomingGhostPrototype<M>>>& rank_to_incoming_ghost_prototypes)
  {
    int rank;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    auto nbrinfo = pinfo.getFineNbrInfo(f);
    for (size_t i = 0; i < Orthant<M>::num_orthants; i++) {
      if (nbrinfo.ranks[i] == rank) {
        my_local_calls.emplace_back(
          f, NbrType::Fine, Orthant<M>(i), pinfo.local_index, nbrinfo.local_indexes[i]);
      } else {
        rank_to_remote_call_prototypes[nbrinfo.ranks[i]].emplace(
          nbrinfo.ids[i], f, NbrType::Fine, Orthant<M>(i), pinfo.local_index);
        rank_to_incoming_ghost_prototypes[nbrinfo.ranks[i]].emplace(pinfo.id, f, pinfo.local_index);
      }
    }
  }
 
  template<int M>
  static void addCoarseNbrCalls(
    const PatchInfo<D>& pinfo,
    Face<D, M> f,
    std::deque<LocalCall<M>>& my_local_calls,
    std::map<int, std::set<RemoteCallPrototype<M>>>& rank_to_remote_call_prototypes,
    std::map<int, std::set<IncomingGhostPrototype<M>>>& rank_to_incoming_ghost_prototypes)
  {
    int rank;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    auto nbrinfo = pinfo.getCoarseNbrInfo(f);
    if (nbrinfo.rank == rank) {
      my_local_calls.emplace_back(
        f, NbrType::Coarse, nbrinfo.orth_on_coarse, pinfo.local_index, nbrinfo.local_index);
    } else {
      rank_to_remote_call_prototypes[nbrinfo.rank].emplace(
        nbrinfo.id, f, NbrType::Coarse, nbrinfo.orth_on_coarse, pinfo.local_index);
      rank_to_incoming_ghost_prototypes[nbrinfo.rank].emplace(pinfo.id, f, pinfo.local_index);
    }
  }
 
  template<int M>
  void enumerateCalls(std::deque<LocalCall<M>>& my_local_calls,
                      std::map<int, RemoteCallSet>& rank_to_remote_call_sets) const
  {
    int rank;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    std::map<int, std::set<RemoteCallPrototype<M>>> rank_to_remote_call_prototypes;
    std::map<int, std::set<IncomingGhostPrototype<M>>> rank_to_incoming_ghost_prototypes;
 
    for (const PatchInfo<D>& pinfo : domain.getPatchInfoVector()) {
      for (Face<D, M> f : Face<D, M>::getValues()) {
        if (pinfo.hasNbr(f)) {
          switch (pinfo.getNbrType(f)) {
            case NbrType::Normal:
              addNormalNbrCalls(pinfo,
                                f,
                                my_local_calls,
                                rank_to_remote_call_prototypes,
                                rank_to_incoming_ghost_prototypes);
              break;
            case NbrType::Fine:
              addFineNbrCalls(pinfo,
                              f,
                              my_local_calls,
                              rank_to_remote_call_prototypes,
                              rank_to_incoming_ghost_prototypes);
              break;
            case NbrType::Coarse:
              addCoarseNbrCalls(pinfo,
                                f,
                                my_local_calls,
                                rank_to_remote_call_prototypes,
                                rank_to_incoming_ghost_prototypes);
              break;
            default:
              throw RuntimeError("Unsupported NbrType");
          }
        }
      }
    }
    const GhostViewInfo<M>& ghost_local_data_info = ghost_local_data_infos.template get<M>();
    for (const auto& pair : rank_to_remote_call_prototypes) {
      int rank = pair.first;
      rank_to_remote_call_sets.emplace(rank, rank);
      RemoteCallSet& remote_call_set = rank_to_remote_call_sets.at(rank);
      std::tuple<int, Face<D, M>> prev_id_side;
      for (const RemoteCallPrototype<M>& call : pair.second) {
        size_t offset = remote_call_set.send_buffer_length;
 
        // calculate length in buffer need for ghost cells
        size_t length = ghost_local_data_info.getSize(call.face);
        // add length to buffer length
        // if its the same side of the patch, resuse the previous buffer space
        if (std::make_tuple(call.id, call.face) == prev_id_side) {
          offset -= length;
        } else {
          remote_call_set.send_buffer_length += length;
        }
        remote_call_set.remote_calls.template get<M>().emplace_back(call, offset);
        prev_id_side = std::make_tuple(call.id, call.face);
      }
    }
    for (const auto& pair : rank_to_incoming_ghost_prototypes) {
      RemoteCallSet& remote_call_set = rank_to_remote_call_sets.at(pair.first);
      for (const IncomingGhostPrototype<M>& prototype : pair.second) {
        // add length for ghosts to buffer length
        size_t length = ghost_local_data_info.getSize(prototype.face);
        size_t offset = remote_call_set.recv_buffer_length;
        remote_call_set.recv_buffer_length += length;
 
        // add ghost to incoming ghosts
        remote_call_set.incoming_ghosts.template get<M>().emplace_back(prototype, offset);
      }
    }
  }
 
  template<int M>
  void zeroGhostCellsOnAllFaces(const PatchInfo<D>& pinfo, PatchView<const double, D>& view) const
  {
    for (Face<D, M> f : Face<D, M>::getValues()) {
      if (pinfo.hasNbr(f)) {
        std::array<size_t, D - M> start;
        start.fill(0);
        std::array<size_t, D - M> end;
        end.fill(domain.getNumGhostCells() - 1);
        Loop::Nested<D - M>(start, end, [&](const std::array<size_t, D - M>& offset) {
          View<double, M + 1> this_ghost = view.getGhostSliceOn(f, offset);
          Loop::OverInteriorIndexes<M + 1>(
            this_ghost, [&](const std::array<int, M + 1>& coord) { this_ghost[coord] = 0; });
        });
      }
    }
  }
 
  void zeroGhostCells(const Vector<D>& u) const
  {
    for (const PatchInfo<D>& pinfo : domain.getPatchInfoVector()) {
      PatchView<const double, D> this_patch = u.getPatchView(pinfo.local_index);
      switch (fill_type) {
        case GhostFillingType::Corners:
          if constexpr (D >= 2) {
            zeroGhostCellsOnAllFaces<0>(pinfo, this_patch);
          }
        case GhostFillingType::Edges:
          if constexpr (D == 3) {
            zeroGhostCellsOnAllFaces<1>(pinfo, this_patch);
          }
        case GhostFillingType::Faces:
          zeroGhostCellsOnAllFaces<D - 1>(pinfo, this_patch);
          break;
        default:
          throw RuntimeError("Unsupported GhostFilling Type");
      }
    }
  }
 
public:
  MPIGhostFiller(const Domain<D>& domain, GhostFillingType fill_type)
    : domain(domain)
    , fill_type(fill_type)
  {
    std::map<int, RemoteCallSet> rank_to_remote_call_sets;
 
    switch (fill_type) {
      case GhostFillingType::Corners:
        if constexpr (D >= 2) {
          ghost_local_data_infos.template get<0>() = GhostViewInfo<0>(domain);
          enumerateCalls<0>(local_calls.template get<0>(), rank_to_remote_call_sets);
        }
      case GhostFillingType::Edges:
        if constexpr (D == 3) {
          ghost_local_data_infos.template get<1>() = GhostViewInfo<1>(domain);
          enumerateCalls<1>(local_calls.template get<1>(), rank_to_remote_call_sets);
        }
      case GhostFillingType::Faces:
        ghost_local_data_infos.template get<D - 1>() = GhostViewInfo<D - 1>(domain);
        enumerateCalls<D - 1>(local_calls.template get<D - 1>(), rank_to_remote_call_sets);
        break;
      default:
        throw RuntimeError("Unsupported GhostFilling Type");
    }
 
    remote_call_sets.reserve(rank_to_remote_call_sets.size());
    for (const auto& pair : rank_to_remote_call_sets) {
      remote_call_sets.push_back(pair.second);
    }
  }
  virtual void fillGhostCellsForNbrPatch(const PatchInfo<D>& pinfo,
                                         const PatchView<const double, D>& local_view,
                                         const PatchView<const double, D>& nbr_view,
                                         Side<D> side,
                                         NbrType nbr_type,
                                         Orthant<D - 1> orthant_on_coarse) const = 0;
  virtual void fillGhostCellsForEdgeNbrPatch(const PatchInfo<D>& pinfo,
                                             const PatchView<const double, D>& local_view,
                                             const PatchView<const double, D>& nbr_view,
                                             Edge edge,
                                             NbrType nbr_type,
                                             Orthant<1> orthant_on_coarse) const = 0;
  virtual void fillGhostCellsForCornerNbrPatch(const PatchInfo<D>& pinfo,
                                               const PatchView<const double, D>& local_view,
                                               const PatchView<const double, D>& nbr_view,
                                               Corner<D> corner,
                                               NbrType nbr_type) const = 0;
 
  virtual void fillGhostCellsForLocalPatch(const PatchInfo<D>& pinfo,
                                           const PatchView<const double, D>& view) const = 0;
 
  void fillGhost(const Vector<D>& u) const override
  {
    if constexpr (ENABLE_DEBUG) {
      if (u.getNumLocalPatches() != domain.getNumLocalPatches()) {
        throw RuntimeError("u vector is incorrect length. Expected Lenght of " +
                           std::to_string(domain.getNumLocalPatches()) + " but vector was length " +
                           std::to_string(u.getNumLocalPatches()));
      }
    }
    // zero out ghost cells
    zeroGhostCells(u);
 
    // allocate recv buffers and post recvs
    std::vector<std::vector<double>> recv_buffers(remote_call_sets.size());
    for (size_t i = 0; i < remote_call_sets.size(); i++) {
      recv_buffers[i].resize(remote_call_sets[i].recv_buffer_length * u.getNumComponents());
    }
    std::vector<MPI_Request> recv_requests = postRecvs(recv_buffers);
 
    // allocate send buffers
    std::vector<std::vector<double>> out_buffers(remote_call_sets.size());
    for (size_t i = 0; i < remote_call_sets.size(); i++) {
      out_buffers[i].resize(remote_call_sets[i].send_buffer_length * u.getNumComponents());
    }
    std::vector<MPI_Request> send_requests = postSends(out_buffers, u);
 
    // perform local operations
    for (const PatchInfo<D>& pinfo : domain.getPatchInfoVector()) {
      PatchView<const double, D> view = u.getPatchView(pinfo.local_index);
      fillGhostCellsForLocalPatch(pinfo, view);
    }
    processLocalFills<D - 1>(u);
 
    processRecvs(recv_requests, recv_buffers, u);
 
    // wait for sends for finish
    MPI_Waitall(send_requests.size(), send_requests.data(), MPI_STATUS_IGNORE);
  }
 
  GhostFillingType getFillType() const { return fill_type; }
 
  const Domain<D>& getDomain() const { return domain; }
};
extern template class MPIGhostFiller<1>;
extern template class MPIGhostFiller<2>;
extern template class MPIGhostFiller<3>;
} // namespace ThunderEgg
#endif