source: src/comm/comm_mpi.cpp@ 97c25dd

/**
 * @file   comm_mpi.cpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Wed Jun 16 13:21:06 2010
 *
 * @brief  Class for MPI-based communication.
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#ifdef HAVE_MPI

#include <mpi.h>

#ifdef HAVE_BOOST_FILESYSTEM
#include <boost/filesystem.hpp>
namespace fs = boost::filesystem;
#endif

#include <cstring>
#include <sstream>

#include "base/helper.hpp"
#include "base/linked_cell_list.hpp"
#include "base/timer.hpp"
#include "comm/comm_mpi.hpp"
#include "grid/grid.hpp"
#include "grid/multigrid.hpp"
#include "grid/tempgrid.hpp"
#include "mg.hpp"

static char print_buffer[512];

using namespace VMG;

void CommMPI::CommSubgrid(Grid& grid_old, Grid& grid_new)
{
  Timer::Start("CommSubgrid");

  MPI_Comm comm = comm_info.GetUnionCommunicator(grid_old, grid_new);

  if (comm != MPI_COMM_NULL) {

    int comm_rank, comm_size;
    MPI_Comm_rank(comm, &comm_rank);
    MPI_Comm_size(comm, &comm_size);

    Index begin_local, end_local;
    Index begin_remote, end_remote;
    Index sizes, subsizes, starts;
    Index offset_new, offset_old;
    int* buffer = new int[6*comm_size];

    for (int i=0; i<3; ++i) {

      offset_old[i] = (grid_old.Local().HaloEnd1()[i] > 0 ? grid_old.Local().Begin()[i] : 0);
      offset_new[i] = (grid_new.Local().HaloEnd1()[i] > 0 ? grid_new.Local().Begin()[i] : 0);

      buffer[6*comm_rank + i  ] = grid_new.Global().BeginLocal()[i];
      buffer[6*comm_rank + i+3] = grid_new.Global().EndLocal()[i];

    }

    //TODO: Create LUT to avoid global communication here
    MPI_Allgather(MPI_IN_PLACE, 6, MPI_INT, buffer, 6, MPI_INT, comm);

    for (int i=0; i<comm_size; ++i) {

      for (int j=0; j<3; ++j) {

        begin_remote[j] = buffer[6*i + j];
        end_remote[j] = buffer[6*i + j+3];

      }

      begin_remote = begin_remote.Clamp(grid_old.Global().BeginLocal(), grid_old.Global().EndLocal());
      end_remote = end_remote.Clamp(grid_old.Global().BeginLocal(), grid_old.Global().EndLocal());

      if ((end_remote - begin_remote).Product() > 0) {

        MPI_Datatype dt_temp;

        sizes = grid_old.Local().SizeTotal();
        subsizes = end_remote - begin_remote;
        starts = begin_remote - grid_old.Global().BeginLocal() + offset_old;

        MPI_Type_create_subarray(3, sizes.vec(), subsizes.vec(), starts.vec(), MPI_ORDER_C, MPI_DOUBLE, &dt_temp);

        MPI_Type_commit(&dt_temp);

        MPI_Isend(&grid_old(0), 1, dt_temp, i, 0, comm, &Request());

        MPI_Type_free(&dt_temp);

      }

    }

    for (int i=0; i<3; ++i) {

      buffer[6*comm_rank + i  ] = grid_old.Global().BeginLocal()[i];
      buffer[6*comm_rank + i+3] = grid_old.Global().EndLocal()[i];

    }

    //TODO: Create LUT to avoid global communication here
    MPI_Allgather(MPI_IN_PLACE, 6, MPI_INT, buffer, 6, MPI_INT, comm);

    for (int i=0; i<comm_size; ++i) {

      for (int j=0; j<3; ++j) {

        begin_remote[j] = buffer[6*i + j];
        end_remote[j] = buffer[6*i + j+3];

      }

      begin_remote = begin_remote.Clamp(grid_new.Global().BeginLocal(), grid_new.Global().EndLocal());
      end_remote = end_remote.Clamp(grid_new.Global().BeginLocal(), grid_new.Global().EndLocal());

      if ((end_remote - begin_remote).Product() > 0) {

        MPI_Datatype dt_temp;

        sizes = grid_new.Local().SizeTotal();
        subsizes = end_remote - begin_remote;
        starts = begin_remote - grid_new.Global().BeginLocal() + offset_new;

        MPI_Type_create_subarray(3, sizes.vec(), subsizes.vec(), starts.vec(), MPI_ORDER_C, MPI_DOUBLE, &dt_temp);

        MPI_Type_commit(&dt_temp);

        MPI_Irecv(&grid_new(0), 1, dt_temp, i, 0, comm, &Request());

        MPI_Type_free(&dt_temp);

      }

    }

    ClearRequests();

    delete [] buffer;

  }

  CommToGhosts(grid_new);

  Timer::Stop("CommSubgrid");
}

void CommMPI::CommAddSubgrid(Grid& grid_old, Grid& grid_new)
{
  Timer::Start("CommAddSubgrid");

  MPI_Comm comm = comm_info.GetUnionCommunicator(grid_old, grid_new);

  if (comm != MPI_COMM_NULL) {

    int comm_rank, comm_size;
    MPI_Comm_rank(comm, &comm_rank);
    MPI_Comm_size(comm, &comm_size);

    Index begin_local, end_local;
    Index begin_remote, end_remote;
    Index sizes, subsizes, starts;
    int* buffer = new int[6*comm_size];

    for (int i=0; i<3; ++i) {

      buffer[6*comm_rank + i  ] = grid_new.Global().BeginLocal()[i];
      buffer[6*comm_rank + i+3] = grid_new.Global().EndLocal()[i];

    }

    //TODO: Create LUT to avoid global communication here
    MPI_Allgather(MPI_IN_PLACE, 6, MPI_INT, buffer, 6, MPI_INT, comm);

    for (int i=0; i<comm_size; ++i) {
      for (int j=0; j<3; ++j) {

        begin_remote[j] = buffer[6*i + j];
        end_remote[j] = buffer[6*i + j+3];

      }

      begin_remote = begin_remote.Clamp(grid_old.Global().BeginLocal(), grid_old.Global().EndLocal());
      end_remote = end_remote.Clamp(grid_old.Global().BeginLocal(), grid_old.Global().EndLocal());

      if ((end_remote - begin_remote).Product() > 0) {

        MPI_Datatype dt_temp;

        sizes = grid_old.Local().SizeTotal();
        subsizes = end_remote - begin_remote;
        starts = begin_remote - grid_old.Global().BeginLocal();

        MPI_Type_create_subarray(3, sizes.vec(), subsizes.vec(), starts.vec(), MPI_ORDER_C, MPI_DOUBLE, &dt_temp);

        MPI_Type_commit(&dt_temp);

        MPI_Isend(&grid_old(0), 1, dt_temp, i, 0, comm, &Request());

        MPI_Type_free(&dt_temp);

      }

    }

    for (int i=0; i<3; ++i) {

      buffer[6*comm_rank + i  ] = grid_old.Global().BeginLocal()[i];
      buffer[6*comm_rank + i+3] = grid_old.Global().EndLocal()[i];

    }

    //TODO: Create LUT to avoid global communication here
    MPI_Allgather(MPI_IN_PLACE, 6, MPI_INT, buffer, 6, MPI_INT, comm);

    for (int i=0; i<comm_size; ++i) {

      for (int j=0; j<3; ++j) {

        begin_remote[j] = buffer[6*i + j];
        end_remote[j] = buffer[6*i + j+3];

      }

      begin_remote = begin_remote.Clamp(grid_new.Global().BeginLocal(), grid_new.Global().EndLocal());
      end_remote = end_remote.Clamp(grid_new.Global().BeginLocal(), grid_new.Global().EndLocal());

      if ((end_remote - begin_remote).Product() > 0) {

        MPI_Datatype dt_temp;

        sizes = grid_new.Local().SizeTotal();
        subsizes = end_remote - begin_remote;
        starts = begin_remote - grid_old.Global().BeginLocal();

        MPI_Type_create_subarray(3, sizes.vec(), subsizes.vec(), starts.vec(), MPI_ORDER_C, MPI_DOUBLE, &dt_temp);

        MPI_Type_commit(&dt_temp);

        ReceiveAndAddSubgrid(grid_new, comm, dt_temp, i, 0);

        MPI_Type_free(&dt_temp);

      }

    }

    ClearRequests();

    delete [] buffer;

  }

  CommToGhosts(grid_new);

  Timer::Stop("CommAddSubgrid");
}

void CommMPI::CommFromGhosts(Grid& grid)
{
  Timer::Start("CommFromGhosts");

  MPI_Comm comm = comm_info.GetCommunicator(grid);

  if (comm != MPI_COMM_NULL) {

    const Index has_halo_1 = grid.Local().HasHalo1();
    const Index has_halo_2 = grid.Local().HasHalo2();
    Tuple<int> neighbors;

    for (int i=0; i<3; ++i) {

      MPI_Cart_shift(comm, i, 1, &neighbors.Left(), &neighbors.Right());

      if (has_halo_1[i]) {
        MPI_Datatype dts1 = comm_info.GetDatatypeSubarray(grid, grid.Iterators().Halo1()[i]);
        MPI_Isend(&grid(0), 1, dts1, neighbors.Left(),  1, comm, &Request());
      }

      if (has_halo_2[i]) {

        MPI_Datatype dts2 = comm_info.GetDatatypeSubarray(grid, grid.Iterators().Halo2()[i]);
        MPI_Isend(&grid(0), 1, dts2, neighbors.Right(), 2, comm, &Request());

        MPI_Datatype dtr1 = comm_info.GetDatatypeSubarray(grid, grid.Iterators().NearBoundary2()[i]);
        ReceiveAndAddSubgrid(grid, comm, dtr1, neighbors.Right(), 1);

      }

      if (has_halo_1[i]) {
        MPI_Datatype dtr2 = comm_info.GetDatatypeSubarray(grid, grid.Iterators().NearBoundary1()[i]);
        ReceiveAndAddSubgrid(grid, comm, dtr2, neighbors.Left(), 2);
      }

    }

    ClearRequests();

  }

  Timer::Stop("CommFromGhosts");
}

void CommMPI::CommToGhosts(Grid& grid)
{
  Timer::Start("CommToGhosts");

  Timer::Start("CommToGhostsGetCommunicator");
  MPI_Comm comm = comm_info.GetCommunicator(grid);
  Timer::Stop("CommToGhostsGetCommunicator");

  if (comm != MPI_COMM_NULL) {

    const Index has_halo_1 = grid.Local().HasHalo1();
    const Index has_halo_2 = grid.Local().HasHalo2();
    Tuple<int> neighbors;

    for (int i=2; i>=0; --i) {

      Timer::Start("CommToGhostsCartShift");
      MPI_Cart_shift(comm, i, 1, &neighbors.Left(), &neighbors.Right());
      Timer::Stop("CommToGhostsCartShift");

      MPI_Datatype dts_left, dts_right;
      MPI_Datatype dtr_left, dtr_right;
      int num_left, num_right;

      if (has_halo_1[i]) {
        Timer::Start("CommToGhostsGetDatatypeSubarray");
        dts_left = comm_info.GetDatatypeSubarray(grid, grid.Iterators().NearBoundary1()[i]);
        dtr_left = comm_info.GetDatatypeSubarray(grid, grid.Iterators().Halo1()[i]);
        Timer::Stop("CommToGhostsGetDatatypeSubarray");
        num_left = 1;
      }else {
        dts_left = MPI_INT;
        dtr_left = MPI_INT;
        num_left = 0;
      }

      if (has_halo_2[i]) {
        Timer::Start("CommToGhostsGetDatatypeSubarray");
        dts_right = comm_info.GetDatatypeSubarray(grid, grid.Iterators().NearBoundary2()[i]);
        dtr_right = comm_info.GetDatatypeSubarray(grid, grid.Iterators().Halo2()[i]);
        Timer::Stop("CommToGhostsGetDatatypeSubarray");
        num_right = 1;
      }else {
        dts_right = MPI_INT;
        dtr_right = MPI_INT;
        num_right = 0;
      }

      Timer::Start("CommToGhostsCommunication");
      MPI_Sendrecv(&(grid(0)), num_right, dts_right, neighbors.Right(), 3,
                   &(grid(0)), num_left,  dtr_left,  neighbors.Left(),  3,
                   comm, MPI_STATUS_IGNORE);
      MPI_Sendrecv(&(grid(0)), num_left,  dts_left,  neighbors.Left(),  4,
                   &(grid(0)), num_right, dtr_right, neighbors.Right(), 4,
                   comm, MPI_STATUS_IGNORE);
      Timer::Stop("CommToGhostsCommunication");

    }
  }

  Timer::Stop("CommToGhosts");
}

void CommMPI::CommToGhostsRB(Grid& grid, const int& offset)
{
  MPI_Comm comm = comm_info.GetCommunicator(grid);

  if (comm != MPI_COMM_NULL) {

    Index i;
    int num_left, num_right;
    MPI_Datatype dts_left, dts_right;
    MPI_Datatype dtr_left, dtr_right;
    Tuple<int> neighbors;
    std::vector<vmg_float> buffer_send_left, buffer_send_right;
    std::vector<vmg_float> buffer_recv_left, buffer_recv_right;
    std::vector<vmg_float>::const_iterator iter;

    const Index halo_size_1 = grid.Local().HaloEnd1() - grid.Local().HaloBegin1();
    const Index halo_size_2 = grid.Local().HaloEnd2() - grid.Local().HaloBegin2();

    buffer_send_left.clear();
    buffer_send_right.clear();
    buffer_recv_left.clear();
    buffer_recv_right.clear();

    MPI_Cart_shift(comm, 2, 1, &neighbors.Left(), &neighbors.Right());

    if (halo_size_1.Z()) {

      buffer_send_left.reserve(grid.Local().Size().X() *
                               grid.Local().Size().Y() *
                               halo_size_1.Z() / 2);

      for (i.X()=grid.Local().Begin().X(); i.X()<grid.Local().End().X(); ++i.X())
        for (i.Y()=grid.Local().Begin().Y(); i.Y()<grid.Local().End().Y(); ++i.Y())
          for (i.Z()=grid.Local().Begin().Z() + (i.X()+i.Y()+offset)%2;
               i.Z()<grid.Local().Begin().Z() + halo_size_1.Z();
               i.Z() += 2)
            buffer_send_left.push_back(grid.GetVal(i));
    }

    if (halo_size_2.Z()) {

      buffer_send_right.reserve(grid.Local().Size().X() *
                                grid.Local().Size().Y() *
                                halo_size_2.Z() / 2);

      for (i.X()=grid.Local().Begin().X(); i.X()<grid.Local().End().X(); ++i.X())
        for (i.Y()=grid.Local().Begin().Y(); i.Y()<grid.Local().End().Y(); ++i.Y())
          for (i.Z()=grid.Local().End().Z() - halo_size_2.Z() + (i.X()+i.Y()+offset)%2;
               i.Z()<grid.Local().End().Z();
               i.Z() += 2)
            buffer_send_right.push_back(grid.GetVal(i));

    }

    buffer_recv_left.resize(buffer_send_left.size());
    buffer_recv_right.resize(buffer_send_right.size());

    MPI_Sendrecv(&buffer_send_right.front(), buffer_send_right.size(), MPI_DOUBLE, neighbors.Right(), 7,
                 &buffer_recv_left.front(), buffer_recv_left.size(), MPI_DOUBLE, neighbors.Left(), 7,
                 comm, MPI_STATUS_IGNORE);

    MPI_Sendrecv(&buffer_send_left.front(), buffer_send_left.size(), MPI_DOUBLE, neighbors.Left(), 8,
                 & buffer_recv_right.front(), buffer_recv_right.size(), MPI_DOUBLE, neighbors.Right(), 8,
                 comm, MPI_STATUS_IGNORE);

    if (halo_size_1.Z()) {

      iter = buffer_recv_left.begin();
      for (i.X()=grid.Local().Begin().X(); i.X()<grid.Local().End().X(); ++i.X())
        for (i.Y()=grid.Local().Begin().Y(); i.Y()<grid.Local().End().Y(); ++i.Y())
          for (i.Z() = (i.X()+i.Y()+offset)%2; i.Z() < grid.Local().Begin().Z(); i.Z() += 2)
            grid(i) = *iter++;

    }

    if (halo_size_2.Z()) {

      iter = buffer_recv_right.begin();
      for (i.X()=grid.Local().Begin().X(); i.X()<grid.Local().End().X(); ++i.X())
        for (i.Y()=grid.Local().Begin().Y(); i.Y()<grid.Local().End().Y(); ++i.Y())
          for (i.Z()=grid.Local().End().Z() + (i.X()+i.Y()+offset)%2; i.Z()<grid.Local().SizeTotal().Z(); i.Z() += 2)
            grid(i) = *iter++;

    }

    buffer_send_left.clear();
    buffer_send_right.clear();
    buffer_recv_left.clear();
    buffer_recv_right.clear();

    MPI_Cart_shift(comm, 1, 1, &neighbors.Left(), &neighbors.Right());

    if (halo_size_1.Y()) {

      buffer_send_left.reserve(grid.Local().Size().X() *
                               halo_size_1.Y() *
                               grid.Local().SizeTotal().Z() / 2);

      for (i.X()=grid.Local().Begin().X(); i.X()<grid.Local().End().X(); ++i.X())
        for (i.Y()=grid.Local().Begin().Y(); i.Y()<grid.Local().Begin().Y() + halo_size_1.Y(); ++i.Y())
          for (i.Z()=(i.X()+i.Y()+offset)%2; i.Z()<grid.Local().SizeTotal().Z(); i.Z() += 2)
            buffer_send_left.push_back(grid.GetVal(i));
    }

    if (halo_size_2.Y()) {

      buffer_send_right.reserve(grid.Local().Size().X() *
                                halo_size_2.Y() *
                                grid.Local().SizeTotal().Z() / 2);

      for (i.X()=grid.Local().Begin().X(); i.X()<grid.Local().End().X(); ++i.X())
        for (i.Y()=grid.Local().End().Y() - halo_size_2.Y(); i.Y()<grid.Local().End().Y(); ++i.Y())
          for (i.Z()=(i.X()+i.Y()+offset)%2; i.Z()<grid.Local().SizeTotal().Z(); i.Z() += 2)
            buffer_send_right.push_back(grid.GetVal(i));

    }

    buffer_recv_left.resize(buffer_send_left.size());
    buffer_recv_right.resize(buffer_send_right.size());

    MPI_Sendrecv(&buffer_send_right.front(), buffer_send_right.size(), MPI_DOUBLE, neighbors.Right(), 7,
                 &buffer_recv_left.front(), buffer_recv_left.size(), MPI_DOUBLE, neighbors.Left(), 7,
                 comm, MPI_STATUS_IGNORE);

    MPI_Sendrecv(&buffer_send_left.front(), buffer_send_left.size(), MPI_DOUBLE, neighbors.Left(), 8,
                 & buffer_recv_right.front(), buffer_recv_right.size(), MPI_DOUBLE, neighbors.Right(), 8,
                 comm, MPI_STATUS_IGNORE);

    if (halo_size_1.Y()) {

      iter = buffer_recv_left.begin();
      for (i.X()=grid.Local().Begin().X(); i.X()<grid.Local().End().X(); ++i.X())
        for (i.Y()=0; i.Y()<grid.Local().Begin().Y(); ++i.Y())
          for (i.Z()=(i.X()+i.Y()+offset)%2; i.Z()<grid.Local().SizeTotal().Z(); i.Z() += 2)
            grid(i) = *iter++;

    }

    if (halo_size_2.Y()) {

      iter = buffer_recv_right.begin();
      for (i.X()=grid.Local().Begin().X(); i.X()<grid.Local().End().X(); ++i.X())
        for (i.Y()=grid.Local().End().Y(); i.Y()<grid.Local().SizeTotal().Y(); ++i.Y())
          for (i.Z()=(i.X()+i.Y()+offset)%2; i.Z()<grid.Local().SizeTotal().Z(); i.Z() += 2)
            grid(i) = *iter++;

    }

    buffer_send_left.clear();
    buffer_send_right.clear();
    buffer_recv_left.clear();
    buffer_recv_right.clear();

    MPI_Cart_shift(comm, 0, 1, &neighbors.Left(), &neighbors.Right());

    if (halo_size_1.X()) {

      buffer_send_left.reserve(halo_size_1.X() *
                               grid.Local().SizeTotal().Y() *
                               grid.Local().SizeTotal().Z() / 2);

      for (i.X()=grid.Local().Begin().X(); i.X()<grid.Local().Begin().X() + halo_size_1.X(); ++i.X())
        for (i.Y()=0; i.Y()<grid.Local().SizeTotal().Y(); ++i.Y())
          for (i.Z()=(i.X()+i.Y()+offset)%2; i.Z()<grid.Local().SizeTotal().Z(); i.Z() += 2)
            buffer_send_left.push_back(grid.GetVal(i));
    }

    if (halo_size_2.X()) {

      buffer_send_right.reserve(halo_size_2.X() *
                                grid.Local().SizeTotal().Y() *
                                grid.Local().SizeTotal().Z() / 2);

      for (i.X()=grid.Local().End().X() - halo_size_2.X(); i.X()<grid.Local().End().X(); ++i.X())
        for (i.Y()=0; i.Y()<grid.Local().SizeTotal().Y(); ++i.Y())
          for (i.Z()=(i.X()+i.Y()+offset)%2; i.Z()<grid.Local().SizeTotal().Z(); i.Z() += 2)
            buffer_send_right.push_back(grid.GetVal(i));

    }

    buffer_recv_left.resize(buffer_send_left.size());
    buffer_recv_right.resize(buffer_send_right.size());

    MPI_Sendrecv(&buffer_send_right.front(), buffer_send_right.size(), MPI_DOUBLE, neighbors.Right(), 7,
                 &buffer_recv_left.front(), buffer_recv_left.size(), MPI_DOUBLE, neighbors.Left(), 7,
                 comm, MPI_STATUS_IGNORE);

    MPI_Sendrecv(&buffer_send_left.front(), buffer_send_left.size(), MPI_DOUBLE, neighbors.Left(), 8,
                 & buffer_recv_right.front(), buffer_recv_right.size(), MPI_DOUBLE, neighbors.Right(), 8,
                 comm, MPI_STATUS_IGNORE);

    if (halo_size_1.X()) {

      iter = buffer_recv_left.begin();
      for (i.X()=0; i.X()<grid.Local().Begin().X(); ++i.X())
        for (i.Y()=0; i.Y()<grid.Local().SizeTotal().Y(); ++i.Y())
          for (i.Z()=(i.X()+i.Y()+offset)%2; i.Z()<grid.Local().SizeTotal().Z(); i.Z() += 2)
            grid(i) = *iter++;
    }

    if (halo_size_2.X()) {

      iter = buffer_recv_right.begin();
      for (i.X()=grid.Local().End().X(); i.X()<grid.Local().SizeTotal().X(); ++i.X())
        for (i.Y()=0; i.Y()<grid.Local().SizeTotal().Y(); ++i.Y())
          for (i.Z()=(i.X()+i.Y()+offset)%2; i.Z()<grid.Local().SizeTotal().Z(); i.Z() += 2)
            grid(i) = *iter++;

    }

  }
}

Grid& CommMPI::GetCoarserGrid(Multigrid& multigrid)
{
  std::map<const Grid*, TempGrid*>::iterator iter = coarser_grids.find(&multigrid());

  if (iter == coarser_grids.end()) {
    TempGrid* temp_grid = new TempGrid();
    temp_grid->SetPropertiesToCoarser(multigrid(), BoundaryConditions());
    iter = coarser_grids.insert(std::make_pair(&multigrid(), temp_grid)).first;
  }

  TempGrid& grid_coarser = *iter->second;

  CommSubgrid(multigrid(multigrid.Level()-1), grid_coarser);

  return grid_coarser;
}

Grid& CommMPI::GetFinerGrid(Multigrid& multigrid)
{
  std::map<const Grid*, TempGrid*>::iterator iter = finer_grids.find(&multigrid());

  if (iter == finer_grids.end()) {
    TempGrid* temp_grid = new TempGrid();
    temp_grid->SetPropertiesToFiner(multigrid(), BoundaryConditions());
    iter = finer_grids.insert(std::make_pair(&multigrid(), temp_grid)).first;
  }

  TempGrid& grid_finer = *iter->second;

  CommSubgrid(multigrid(multigrid.Level()+1), grid_finer);

  return grid_finer;
}

Grid& CommMPI::GetGlobalCoarseGrid(Multigrid& multigrid)
{
  if (global_coarse_grid == NULL) {

    global_coarse_grid = new TempGrid();
    global_coarse_grid->SetPropertiesToGlobalCoarseGrid();

  }

  CommSubgrid(multigrid(multigrid.MinLevel()), *global_coarse_grid);

  return *global_coarse_grid;
}

//
// Decodes a derived datatype created by MPI_Type_create_subarray
// and adds values to grid.
// I chose this approach since MPI neither supports a point-to-point reduce nor
// a non-blocking global reduce.
// For more information see Chapter 4.1 "DERIVED DATATYPES" of the MPI standard
//
void CommMPI::ReceiveAndAddSubgrid(Grid& grid, const MPI_Comm& comm, MPI_Datatype& type,
                                   const int& rank, const int& tag)
{
  int counter = 0;
  int ibuffer[11];
  Index i;
  MPI_Datatype base_type;
  MPI_Aint address;

  MPI_Type_get_contents(type, 11, 0, 1, ibuffer, &address, &base_type);

  if (static_cast<int>(receive_buffer.size()) < ibuffer[4]*ibuffer[5]*ibuffer[6])
    receive_buffer.resize(ibuffer[4]*ibuffer[5]*ibuffer[6]);

  MPI_Recv(&(receive_buffer.front()), ibuffer[4]*ibuffer[5]*ibuffer[6], MPI_DOUBLE, rank, tag, comm, MPI_STATUS_IGNORE);

  for (i.X()=ibuffer[7]; i.X()<ibuffer[7]+ibuffer[4]; ++i.X())
    for (i.Y()=ibuffer[8]; i.Y()<ibuffer[8]+ibuffer[5]; ++i.Y())
      for (i.Z()=ibuffer[9]; i.Z()<ibuffer[9]+ibuffer[6]; ++i.Z())
        grid(i) += receive_buffer[counter++];
}

void CommMPI::CommParticles(const Grid& grid, std::list<Particle::Particle>& particles)
{
  Timer::Start("CommParticles");

  Factory& factory = MG::GetFactory();

  const vmg_int& num_particles_local = factory.GetObjectStorageVal<vmg_int>("PARTICLE_NUM_LOCAL");
  vmg_float* x = factory.GetObjectStorageArray<vmg_float>("PARTICLE_POS_ARRAY");
  vmg_float* q = factory.GetObjectStorageArray<vmg_float>("PARTICLE_CHARGE_ARRAY");

  int rank, size;
  MPI_Comm_rank(comm_global, &rank);
  MPI_Comm_size(comm_global, &size);

  Index index;
  int global_extent[6*size], send_sizes[size], recv_sizes[size];
  Index begin_remote[size], end_remote[size];
  std::vector<int> pos_indices[size], charge_indices[size];
  MPI_Datatype dt_temp;
  std::vector<vmg_float> buffer_x, buffer_q;
  std::vector<vmg_int> buffer_ind;

  std::memcpy(&global_extent[6*rank], grid.Global().BeginLocal().vec(), 3*sizeof(vmg_float));
  std::memcpy(&global_extent[6*rank+3], grid.Global().EndLocal().vec(), 3*sizeof(vmg_float));

  MPI_Allgather(MPI_IN_PLACE, 6, MPI_INT, global_extent, 6, MPI_INT, MPI_COMM_WORLD);

  for (int i=0; i<size; ++i) {
    begin_remote[i] = static_cast<Vector>(&global_extent[6*i]);
    end_remote[i] = static_cast<Vector>(&global_extent[6*i+3]);
  }

  for (int i=0; i<num_particles_local; ++i) {
    index = static_cast<Index>((Vector(&x[3*i]) - grid.Extent().Begin()) / grid.Extent().MeshWidth());
    for (int j=0; j<size; ++j)
      if (index.IsInBounds(begin_remote[j], end_remote[j])) {
        pos_indices[j].push_back(3*i);
        charge_indices[j].push_back(i);
        break;
      }
  }

  /*
   * Communicate which process gets how many particles
   */
  for (int i=0; i<size; ++i)
    send_sizes[i] = charge_indices[i].size();

  MPI_Alltoall(send_sizes, 1, MPI_INT, recv_sizes, 1, MPI_INT, comm_global);

  /*
   * Send particles
   */
  for (int i=0; i<size; ++i) {

    if (pos_indices[i].size() > 0) {

      MPI_Type_create_indexed_block(pos_indices[i].size(), 3, &pos_indices[i].front(), MPI_DOUBLE, &dt_temp);
      MPI_Type_commit(&dt_temp);
      MPI_Isend(x, 1, dt_temp, i, 0, MPI_COMM_WORLD, &Request());
      MPI_Type_free(&dt_temp);

      MPI_Type_create_indexed_block(charge_indices[i].size(), 1, &charge_indices[i].front(), MPI_DOUBLE, &dt_temp);
      MPI_Type_commit(&dt_temp);
      MPI_Isend(q, 1, dt_temp, i, 1, MPI_COMM_WORLD, &Request());
      MPI_Type_free(&dt_temp);

      MPI_Isend(&charge_indices[i].front(), charge_indices[i].size(), MPI_INT, i, 2, MPI_COMM_WORLD, &Request());

    }

  }

  /*
   * Receive particles
   */
  for (int i=0; i<size; ++i) {

    if (recv_sizes[i] > 0) {

      buffer_x.resize(3*recv_sizes[i]);
      buffer_q.resize(recv_sizes[i]);
      buffer_ind.resize(recv_sizes[i]);

      MPI_Recv(&buffer_x.front(), 3*recv_sizes[i], MPI_DOUBLE, i, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
      MPI_Recv(&buffer_q.front(), recv_sizes[i], MPI_DOUBLE, i, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
      MPI_Recv(&buffer_ind.front(), recv_sizes[i], MPI_INT, i, 2, MPI_COMM_WORLD, MPI_STATUS_IGNORE);

      particles.clear();

      for (int j=0; j<recv_sizes[i]; ++j)
        particles.push_back(Particle::Particle(&buffer_x[3*j], buffer_q[j], 0.0, 0.0, i, buffer_ind[j]));

    }

  }

  ClearRequests();

  Timer::Stop("CommParticles");
}

void CommMPI::CommLCListGhosts(const Grid& grid, Particle::LinkedCellList& lc)
{
  Timer::Start("CommLCListGhosts");

  std::vector<vmg_float> send_buffer_1, send_buffer_2, recv_buffer;
  vmg_int send_size_1, send_size_2, recv_size;
  std::list<Particle::Particle>::iterator p_iter;
  Grid::iterator g_iter;

  Tuple<int> neighbors;

  for (int i=2; i>=0; --i) {

    /*
     * Get ranks of neighbors
     */
    MPI_Cart_shift(comm_global, i, 1, &neighbors.Left(), &neighbors.Right());

    /*
     * Send to left
     */
    send_buffer_1.clear();

    for (g_iter=lc.Iterators().Boundary1()[i].Begin(); g_iter!=lc.Iterators().Boundary1()[i].End(); ++g_iter)
      lc(*g_iter).clear();

    for (g_iter=lc.Iterators().NearBoundary1()[i].Begin(); g_iter!=lc.Iterators().NearBoundary1()[i].End(); ++g_iter)
      for (p_iter=lc(*g_iter).begin(); p_iter!=lc(*g_iter).end(); ++p_iter) {

        if (grid.Global().BeginLocal()[i] == 0)
          for (int j=0; j<3; ++j)
            send_buffer_1.push_back(p_iter->Pos()[j] + grid.Extent().Size()[j]);
        else
          for (int j=0; j<3; ++j)
            send_buffer_1.push_back(p_iter->Pos()[j]);

        send_buffer_1.push_back(p_iter->Charge());

      }

    send_size_1 = send_buffer_1.size();

    MPI_Isend(&send_size_1, 1, MPI_INT, neighbors.Left(), 0, comm_global, &Request());

    if (send_size_1 > 0)
      MPI_Isend(&send_buffer_1.front(), send_size_1, MPI_DOUBLE, neighbors.Left(), 1, comm_global, &Request());

    /*
     * Send to right
     */
    send_buffer_2.clear();

    for (g_iter=lc.Iterators().Boundary2()[i].Begin(); g_iter!=lc.Iterators().Boundary2()[i].End(); ++g_iter)
      lc(*g_iter).clear();

    for (g_iter=lc.Iterators().NearBoundary2()[i].Begin(); g_iter!=lc.Iterators().NearBoundary2()[i].End(); ++g_iter)
      for (p_iter=lc(*g_iter).begin(); p_iter!=lc(*g_iter).end(); ++p_iter) {

        if (grid.Global().EndLocal()[i] == grid.Global().SizeGlobal()[i])
          for (int j=0; j<3; ++j)
            send_buffer_2.push_back(p_iter->Pos()[j] - grid.Extent().Size()[j]);
        else
          for (int j=0; j<3; ++j)
            send_buffer_2.push_back(p_iter->Pos()[j]);

        send_buffer_2.push_back(p_iter->Charge());

      }

    send_size_2 = send_buffer_2.size();

    MPI_Isend(&send_size_2, 1, MPI_INT, neighbors.Right(), 2, comm_global, &Request());

    if (send_size_2 > 0)
      MPI_Isend(&send_buffer_2.front(), send_size_2, MPI_DOUBLE, neighbors.Right(), 3, comm_global, &Request());

    /*
     * Receive from right
     */
    MPI_Recv(&recv_size, 1, MPI_INT, neighbors.Right(), 0, comm_global, MPI_STATUS_IGNORE);

    if (recv_size > 0) {

      recv_buffer.resize(recv_size);
      MPI_Recv(&recv_buffer.front(), recv_size, MPI_DOUBLE, neighbors.Right(), 1, comm_global, MPI_STATUS_IGNORE);

      for (vmg_int i=0; i<recv_size; i+=4)
        lc.AddParticleToComplete(&recv_buffer[i], recv_buffer[i+3]);

    }

    /*
     * Receive from left
     */
    MPI_Recv(&recv_size, 1, MPI_INT, neighbors.Left(), 2, comm_global, MPI_STATUS_IGNORE);

    if (recv_size > 0) {

      recv_buffer.resize(recv_size);
      MPI_Recv(&recv_buffer.front(), recv_size, MPI_DOUBLE, neighbors.Left(), 3, comm_global, MPI_STATUS_IGNORE);

      for (vmg_int i=0; i<recv_size; i+=4)
        lc.AddParticleToComplete(&recv_buffer[i], recv_buffer[i+3]);

    }


    ClearRequests();
  }

  Timer::Stop("CommLCListGhosts");
}

void CommMPI::CommAddPotential(std::list<Particle::Particle>& particles)
{
  Timer::Start("CommAddPotential");

  std::list<Particle::Particle>::iterator iter;

  if (!win_created) {
    vmg_float* p = MG::GetFactory().GetObjectStorageArray<vmg_float>("PARTICLE_POTENTIAL_ARRAY");
    const vmg_int& num_particles_local = MG::GetFactory().GetObjectStorageVal<vmg_int>("PARTICLE_NUM_LOCAL");
    MPI_Win_create(p, num_particles_local*sizeof(vmg_float), sizeof(vmg_float), info, comm_global, &win);
    win_created = true;
  }

  MPI_Win_fence(MPI_MODE_NOPRECEDE, win);

  for (iter=particles.begin(); iter!=particles.end(); ++iter)
    MPI_Accumulate(&iter->Pot(), 1, MPI_DOUBLE, iter->Rank(), iter->Index(), 1, MPI_DOUBLE, MPI_SUM, win);

  MPI_Win_fence(MPI_MODE_NOSTORE | MPI_MODE_NOSUCCEED, win);

  Timer::Stop("CommAddPotential");
}

void CommMPI::CommAddPotential(Particle::LinkedCellList& lc)
{
  Timer::Start("CommAddPotential");

  Grid::iterator g_iter;
  std::list<Particle::Particle>::iterator p_iter;

  if (!win_created) {
    vmg_float* p = MG::GetFactory().GetObjectStorageArray<vmg_float>("PARTICLE_POTENTIAL_ARRAY");
    const vmg_int& num_particles_local = MG::GetFactory().GetObjectStorageVal<vmg_int>("PARTICLE_NUM_LOCAL");
    MPI_Win_create(p, num_particles_local*sizeof(vmg_float), sizeof(vmg_float), info, comm_global, &win);
    win_created = true;
  }

  MPI_Win_fence(MPI_MODE_NOPRECEDE, win);

  for (g_iter=lc.Iterators().Local().Begin(); g_iter!=lc.Iterators().Local().End(); ++g_iter)
    for (p_iter=lc(*g_iter).begin(); p_iter!=lc(*g_iter).end(); ++p_iter)
      MPI_Accumulate(&p_iter->Pot(), 1, MPI_DOUBLE, p_iter->Rank(), p_iter->Index(), 1, MPI_DOUBLE, MPI_SUM, win);

  MPI_Win_fence(MPI_MODE_NOSTORE | MPI_MODE_NOSUCCEED, win);

  Timer::Stop("CommAddPotential");
}

vmg_float CommMPI::GlobalSum(const vmg_float& value)
{
  vmg_float result = value;

  MPI_Allreduce(MPI_IN_PLACE, &result, 1, MPI_DOUBLE, MPI_SUM, comm_global);

  return result;
}

vmg_float CommMPI::GlobalSumRoot(const vmg_float& value)
{
  vmg_float recv_buffer;
  vmg_float send_buffer = value;

  MPI_Reduce(&send_buffer, &recv_buffer, 1, MPI_DOUBLE, MPI_SUM, 0, comm_global);

  return recv_buffer;
}

void CommMPI::GlobalSumArray(vmg_float* array, const vmg_int& size)
{
  MPI_Allreduce(MPI_IN_PLACE, array, size, MPI_DOUBLE, MPI_SUM, comm_global);
}

void CommMPI::PrintString(const char* format, ...)
{
  va_list args;
  va_start(args, format);
  vsprintf(print_buffer, format, args);
  printf("VMG: Rank %d: %s\n", GlobalRank(), print_buffer);
  va_end(args);
}

void CommMPI::PrintStringOnce(const char* format, ...)
{
  if (GlobalRank() == 0) {
    va_list args;
    va_start(args, format);
    vsprintf(print_buffer, format, args);
    printf("VMG: Rank %d: %s\n", GlobalRank(), print_buffer);
    va_end(args);
  }
}

void CommMPI::PrintXML(const std::string& filename, const std::string& xml_data)
{
  MPI_File file;
  std::stringstream path, xml_header;

  pugi::xml_document doc;
  pugi::xml_node node_data = doc.append_child("Global").append_child("NumProcesses").append_child(pugi::node_pcdata);
  node_data.set_value(Helper::ToString(GlobalProcs().Product()).c_str());
  doc.save(xml_header);

  path << OutputPath() << filename;

  char* filename_array = Helper::GetCharArray(path.str());
  char* xml_header_array = Helper::GetCharArray(xml_header.str());
  char* str_array = Helper::GetCharArray(xml_data);

  MPI_File_open(MPI_COMM_SELF, filename_array, MPI_MODE_WRONLY|MPI_MODE_CREATE, MPI_INFO_NULL, &file);
  MPI_File_set_size(file, 0);
  MPI_File_write(file, xml_header_array, xml_header.str().size(), MPI_CHAR, MPI_STATUS_IGNORE);
  MPI_File_write(file, str_array, xml_data.size(), MPI_CHAR, MPI_STATUS_IGNORE);
  MPI_File_close(&file);

  delete [] filename_array;
  delete [] xml_header_array;
  delete [] str_array;
}

void CommMPI::PrintXMLAll(const std::string& filename, const std::string& xml_data)
{
  MPI_File file;
  std::stringstream path;

  path << OutputPath() << filename;

  char* filename_array = Helper::GetCharArray(path.str());
  char* str_array = Helper::GetCharArray(xml_data);

  MPI_File_open(comm_global, filename_array, MPI_MODE_WRONLY|MPI_MODE_CREATE, MPI_INFO_NULL, &file);
  MPI_File_set_size(file, 0);

  if (GlobalRank() == 0) {
    std::stringstream xml_header;
    pugi::xml_document doc;
    pugi::xml_node node_data = doc.append_child("Global").append_child("NumProcesses").append_child(pugi::node_pcdata);
    node_data.set_value(Helper::ToString(GlobalProcs().Product()).c_str());
    doc.save(xml_header);

    char* xml_header_array = Helper::GetCharArray(xml_header.str());

    MPI_File_write_shared(file, xml_header_array, xml_header.str().size(), MPI_CHAR, MPI_STATUS_IGNORE);

    delete [] xml_header_array;
  }

  MPI_File_write_ordered(file, str_array, xml_data.size(), MPI_CHAR, MPI_STATUS_IGNORE);
  MPI_File_close(&file);

  delete [] filename_array;
  delete [] str_array;
}

void CommMPI::PrintAllSettings()
{
  std::stringstream buf;
  MPI_File file;
  MPI_Status status;

  Multigrid* mg = MG::GetFactory().Get("SOL")->Cast<Multigrid>();

  buf << OutputPath() << "settings.txt";
  char *filename = Helper::GetCharArray(buf.str());

  MPI_File_open(comm_global, filename, MPI_MODE_WRONLY|MPI_MODE_CREATE, MPI_INFO_NULL, &file);
  MPI_File_set_size(file, 0);
  MPI_File_close(&file);

  for (int i=mg->MinLevel(); i<=mg->MaxLevel(); ++i) {

    MPI_Comm comm = comm_info.GetCommunicator((*mg)(i));

    if (comm != MPI_COMM_NULL) {

      MPI_File_open(comm, filename, MPI_MODE_WRONLY|MPI_MODE_CREATE|MPI_MODE_APPEND, MPI_INFO_NULL, &file);

      int comm_rank;
      MPI_Comm_rank(comm, &comm_rank);

      if (comm_rank == 0) {

        buf.str("");
        buf << "PROCESS INDEPENDENT INFORMATION:" << std::endl
            << "GRID LEVEL " << i << std::endl
            << "BEGINFINEST: " << (*mg)(i).Global().BeginFinest() << std::endl
            << "ENDFINEST: " << (*mg)(i).Global().EndFinest() << std::endl
            << "SIZEFINEST: " << (*mg)(i).Global().SizeFinest() << std::endl
            << "SIZEGLOBAL: " << (*mg)(i).Global().SizeGlobal() << std::endl
            << "BOUNDARYTYPE: " << (*mg)(i).Global().BoundaryType() << std::endl
            << "EXTENT" << std::endl
            << "SIZE: " << (*mg)(i).Extent().Size() << std::endl
            << "SIZEFACTOR: " << (*mg)(i).Extent().SizeFactor() << std::endl
            << "BEGIN: " << (*mg)(i).Extent().Begin() << std::endl
            << "END: " << (*mg)(i).Extent().End() << std::endl
            << "MESHWIDTH: " << (*mg)(i).Extent().MeshWidth() << std::endl
            << std::endl
            << "PROCESS DEPENDENT INFORMATION:" << std::endl;

        char* char_buf = Helper::GetCharArray(buf.str());

        MPI_File_write_shared(file, char_buf, buf.str().size(), MPI_CHAR, &status);

        delete [] char_buf;

      }

      buf.str("");
      buf << "PROCESS " << comm_rank << ":" << std::endl
          << "GLOBAL" << std::endl
          << "BEGINLOCAL: " << (*mg)(i).Global().BeginLocal() << std::endl
          << "ENDLOCAL: " << (*mg)(i).Global().EndLocal() << std::endl
          << "SIZELOCAL: " << (*mg)(i).Global().SizeLocal() << std::endl
          << "LOCAL" << std::endl
          << "SIZE: " << (*mg)(i).Local().Size() << std::endl
          << "SIZETOTAL: " << (*mg)(i).Local().SizeTotal() << std::endl
          << "BEGIN: " << (*mg)(i).Local().Begin() << std::endl
          << "END: " << (*mg)(i).Local().End() << std::endl
          << "HALOBEGIN1: " << (*mg)(i).Local().HaloBegin1() << std::endl
          << "HALOEND1: " << (*mg)(i).Local().HaloEnd1() << std::endl
          << "HALOBEGIN2: " << (*mg)(i).Local().HaloBegin2() << std::endl
          << "HALOEND2: " << (*mg)(i).Local().HaloEnd2() << std::endl
          << "BOUNDARYBEGIN1: " << (*mg)(i).Local().BoundaryBegin1() << std::endl
          << "BOUNDARYEND1: " << (*mg)(i).Local().BoundaryEnd1() << std::endl
          << "BOUNDARYBEGIN2: " << (*mg)(i).Local().BoundaryBegin2() << std::endl
          << "BOUNDARYEND2: " << (*mg)(i).Local().BoundaryEnd2() << std::endl
          << "ALIGNMENTBEGIN: " << (*mg)(i).Local().AlignmentBegin() << std::endl
          << "ALIGNMENTEND: " << (*mg)(i).Local().AlignmentEnd() << std::endl
          << std::endl;

      char* char_buf = Helper::GetCharArray(buf.str());

      MPI_File_write_ordered(file, char_buf, buf.str().size(), MPI_CHAR, &status);

      delete [] char_buf;

      MPI_File_close(&file);

    }
  }

  delete [] filename;

}

void CommMPI::PrintDefect(Grid& sol, Grid& rhs, const char* information)
{
  TempGrid *temp = MG::GetTempGrid();
  temp->SetProperties(sol);
  temp->ImportFromResidual(sol, rhs);

  PrintInnerGrid(*temp, information);
}

void CommMPI::PrintGrid(Grid& grid, const char* information)
{
  Index i;
  std::stringstream buf;

  Index begin = 0;
  Index end = grid.Local().End();
  Index begin_local, end_local, begin_global, end_global;

  CommToGhosts(grid);

  for (int j=0; j<3; ++j)
    if (grid.Global().EndLocal()[j] == grid.Global().SizeGlobal()[j])
      end[j] = grid.Local().SizeTotal()[j];

  for (i.Z()=begin.Z(); i.Z()<end.Z(); ++i.Z())
    for (i.Y()=begin.Y(); i.Y()<end.Y(); ++i.Y())
      for (i.X()=begin.X(); i.X()<end.X(); ++i.X())
        buf << std::scientific << grid.GetVal(i) << " ";

  begin_local = grid.Global().BeginLocal();
  end_local = grid.Global().EndLocal();
  begin_global = 0;

  for (int i=0; i<3; ++i) {

    if (BoundaryConditions()[i] == Dirichlet ||
        BoundaryConditions()[i] == Open) {

      if (begin_local[i] != 0)
        begin_local[i] -= grid.Local().HaloEnd1()[i] - grid.Local().HaloBegin1()[i];

      end_global[i] = grid.Global().SizeGlobal()[i];

    }else if (BoundaryConditions()[i] == Periodic) {

      if (end_local[i] == grid.Global().SizeGlobal()[i])
        end_local[i] += grid.Local().HaloEnd1()[i] - grid.Local().HaloBegin1()[i] + grid.Local().HaloEnd2()[i] - grid.Local().HaloBegin2()[i];
      else
        end_local[i] += grid.Local().HaloEnd1()[i] - grid.Local().HaloBegin1()[i];

      end_global[i] = grid.Global().SizeGlobal()[i] +
        grid.Local().HaloEnd1()[i] - grid.Local().HaloBegin1()[i] +
        grid.Local().HaloEnd2()[i] - grid.Local().HaloBegin2()[i];

    }

  }

  CreateOutputFiles(grid, buf, information,
                    begin_global, end_global,
                    begin_local, end_local);
}

void CommMPI::PrintCorrectedGrid(Grid& grid, const char* information)
{
  Index i;

  std::stringstream buf;

  Index begin = 0;
  Index end = grid.Local().End();
  Index begin_local, end_local, begin_global, end_global;

  CommToGhosts(grid);

  for (int j=0; j<3; ++j)
    if (grid.Global().EndLocal()[j] == grid.Global().SizeGlobal()[j])
      end[j] = grid.Local().SizeTotal()[j];

  for (i.Z()=begin.Z(); i.Z()<end.Z(); ++i.Z())
    for (i.Y()=begin.Y(); i.Y()<end.Y(); ++i.Y())
      for (i.X()=begin.X(); i.X()<end.X(); ++i.X())
        buf << std::scientific << grid.GetVal(i) - Grid::Correction() << " ";

  begin_local = grid.Global().BeginLocal();
  end_local = grid.Global().EndLocal();
  begin_global = 0;

  for (int i=0; i<3; ++i) {

    if (BoundaryConditions()[i] == Dirichlet ||
        BoundaryConditions()[i] == Open) {

      if (begin_local[i] != 0)
        --begin_local[i];

      end_global[i] = grid.Global().SizeGlobal()[i];

    }else if (BoundaryConditions()[i] == Periodic) {

      if (begin_local[i] != 0)
        --begin_local[i];

      if (end_local[i] == grid.Global().SizeGlobal()[i])
        end_local[i] += 2;

      end_global[i] = grid.Global().SizeGlobal()[i] + 2;

    }

  }

  CreateOutputFiles(grid, buf, information,
                    0, grid.Global().SizeGlobal()+2,
                    begin, end);
}

void CommMPI::PrintInnerGrid(Grid& grid, const char* information)
{
  Index i;
  Index begin_local, end_local, begin_global, end_global;
  std::stringstream buf;

  CommToGhosts(grid);

  Index begin = 0;
  Index end = grid.Local().End();

  /* VTK parallel ImageData files seem to require some overlap. */
  for (int j=0; j<3; ++j)
    if (grid.Global().BeginLocal()[j] == 0)
      ++begin[j];

  /* Write grid values into buffer */
  for (i.Z()=begin.Z(); i.Z()<end.Z(); ++i.Z())
    for (i.Y()=begin.Y(); i.Y()<end.Y(); ++i.Y())
      for (i.X()=begin.X(); i.X()<end.X(); ++i.X())
        buf << std::scientific << grid.GetVal(i) << " ";

  for (int j=0; j<3; ++j) {

    if (BoundaryConditions()[j] == Dirichlet ||
        BoundaryConditions()[j] == Open) {

      begin_global[j] = 1;
      end_global[j] = grid.Global().SizeGlobal()[j] - 1;

      if (grid.Global().BeginLocal()[j] == 0)
        begin_local[j] = 1;
      else
        begin_local[j] = grid.Global().BeginLocal()[j] - 1;

      if (grid.Global().EndLocal()[j] == grid.Global().SizeGlobal()[j])
        end_local[j] = grid.Global().EndLocal()[j] - 1;
      else
        end_local[j] = grid.Global().EndLocal()[j];

    }else if (BoundaryConditions()[j] == Periodic) {

      begin_global[j] = 0;
      end_global[j] = grid.Global().SizeGlobal()[j];

      if (grid.Global().BeginLocal()[j] == 0)
        begin_local[j] = grid.Global().BeginLocal()[j];
      else
        begin_local[j] = grid.Global().BeginLocal()[j] - 1;

      end_local[j] = grid.Global().EndLocal()[j];

    }else {

      assert(0 == "Boundary condition not handled properly.");

    }

  }

  CreateOutputFiles(grid, buf, information,
                    begin_global, end_global,
                    begin_local, end_local);
}

void CommMPI::PrintInnerCorrectedGrid(Grid& grid, const char* information)
{
  Index i;
  const Index begin = grid.Local().Begin();
  const Index end = grid.Local().End();
  std::stringstream buf;

  for (i.Z()=begin.Z(); i.Z()<end.Z(); ++i.Z())
    for (i.Y()=begin.Y(); i.Y()<end.Y(); ++i.Y())
      for (i.X()=begin.X(); i.X()<end.X(); ++i.X())
        buf << std::scientific << grid.GetVal(i) - Grid::Correction() << " ";

  CreateOutputFiles(grid, buf, information,
                    0, grid.Global().SizeGlobal(),
                    grid.Global().BeginLocal(), grid.Global().EndLocal());
}

void CommMPI::CreateOutputFiles(const Grid& grid, const std::stringstream& serial_data, const char* information,
                                const Index& begin_global, const Index& end_global,
                                const Index& begin_local, const Index& end_local)
{
  int output_count = OutputCount();

  MPI_Comm comm = comm_info.GetCommunicator(grid);

  if (comm != MPI_COMM_NULL) {

    MPI_File file;
    std::string conv_information = Helper::ReplaceWhitespaces(information, "_");

    CreateParallelOutputFile(grid, comm, output_count, conv_information.c_str(),
                             begin_global, end_global, begin_local, end_local);

    file = CreateSerialOutputFile(grid, comm, output_count, conv_information.c_str(),
                                  begin_global, end_global, begin_local, end_local);

    char *char_buf = Helper::GetCharArray(serial_data.str());
    MPI_File_write(file, char_buf, serial_data.str().size(), MPI_CHAR, MPI_STATUS_IGNORE);
    delete [] char_buf;

    FinalizeSerialOutputFile(file);

  }
}

void CommMPI::CreateParallelOutputFile(const Grid& grid, MPI_Comm& comm,
                                       const int& output_count, const char* information,
                                       const Index& begin_global, const Index& end_global,
                                       const Index& begin_local, const Index& end_local)
{
  int rank;
  MPI_File file;
  char parallel_filename[513], serial_filename[513];
  std::stringstream buf;

  MPI_Comm_rank(comm, &rank);

  sprintf(parallel_filename, "%s%04d.pvti", OutputPath().c_str(), output_count);
  sprintf(serial_filename, "%04d_%d.vti", output_count, rank);

  MPI_File_open(comm, parallel_filename, MPI_MODE_WRONLY|MPI_MODE_CREATE, MPI_INFO_NULL, &file);
  MPI_File_set_size(file, 0);

  if (rank == 0) {

    buf << "<?xml version=\"1.0\"?>" << std::endl
        << "<VTKFile type=\"PImageData\" version=\"0.1\" byte_order=\"LittleEndian\">" << std::endl
        << "  <PImageData WholeExtent=\"";

    for (int i=0; i<3; ++i)
      buf << begin_global[i] << " " << end_global[i] - 1 << " ";

    buf << "\"" << std::endl
        << "              GhostLevel=\"0\" Origin=\"0 0 0\" Spacing=\"";

    for (int i=0; i<3; ++i)
      buf << grid.Extent().MeshWidth()[i] << " ";

    buf << "\">" << std::endl
        << "    <PPointData Scalars=\"" << information << "\">" << std::endl
        << "      <PDataArray type=\"Float32\" Name=\"" << information << "\"/>" << std::endl
        << "    </PPointData>" << std::endl;

    char* char_buf = Helper::GetCharArray(buf.str());

    MPI_File_write_shared(file, char_buf, buf.str().size(), MPI_CHAR, MPI_STATUS_IGNORE);

    delete [] char_buf;
  }

  buf.str("");
  buf << "    <Piece Extent=\"";

  for (int i=0; i<3; ++i)
    buf << begin_local[i] << " " << end_local[i] - 1  << " ";

  buf << "\" Source=\"" << serial_filename << "\"/>" << std::endl;

  char* char_buf = Helper::GetCharArray(buf.str());

  MPI_File_write_ordered(file, char_buf, buf.str().size(), MPI_CHAR, MPI_STATUS_IGNORE);

  delete [] char_buf;

  if (rank == 0) {

    buf.str("");

    buf << "  </PImageData>" << std::endl
        << "</VTKFile>" << std::endl;

    char* char_buf = Helper::GetCharArray(buf.str());

    MPI_File_write_shared(file, char_buf, buf.str().size(), MPI_CHAR, MPI_STATUS_IGNORE);

    delete [] char_buf;
  }

  MPI_File_close(&file);
}

MPI_File CommMPI::CreateSerialOutputFile(const Grid& grid, MPI_Comm& comm,
                                         const int& output_count, const char* information,
                                         const Index& begin_global, const Index& end_global,
                                         const Index& begin_local, const Index& end_local)
{
  char serial_filename[513];
  int rank;
  MPI_File file;
  std::stringstream buf;

  MPI_Comm_rank(comm_global, &rank);

  sprintf(serial_filename, "%s%04d_%d.vti", OutputPath().c_str(), output_count, rank);

  MPI_File_open(MPI_COMM_SELF, serial_filename, MPI_MODE_WRONLY|MPI_MODE_CREATE, MPI_INFO_NULL, &file);

  buf << "<?xml version=\"1.0\"?>" << std::endl
      << "<VTKFile type=\"ImageData\" version=\"0.1\" byte_order=\"LittleEndian\">" << std::endl
      << "  <ImageData WholeExtent=\"";

  for (int i=0; i<3; ++i)
    buf << begin_global[i] << " " << end_global[i] - 1 << " ";

  buf << "\"" << std::endl
      << "             Origin=\"0 0 0\" Spacing=\"";

  for (int i=0; i<3; ++i)
    buf << grid.Extent().MeshWidth()[i] << " ";

  buf << "\">" << std::endl
      << "    <Piece Extent=\"";

  for (int i=0; i<3; ++i)
    buf << begin_local[i] << " " << end_local[i] - 1 << " ";

  buf << "\">" << std::endl
      << "      <PointData Scalars=\"" << information << "\">" << std::endl
      << "        <DataArray type=\"Float32\" Name=\"" << information << "\" format=\"ascii\">" << std::endl
      << "          ";

  char* char_buf = Helper::GetCharArray(buf.str());
  MPI_File_write(file, char_buf, buf.str().size(), MPI_CHAR, MPI_STATUS_IGNORE);
  delete [] char_buf;

  return file;
}

void CommMPI::FinalizeSerialOutputFile(MPI_File& file)
{
  std::stringstream buf;

  buf << std::endl
      << "        </DataArray>" << std::endl
      << "      </PointData>" << std::endl
      << "    </Piece>" << std::endl
      << "  </ImageData>" << std::endl
      << "</VTKFile>" << std::endl;

  char* char_buf = Helper::GetCharArray(buf.str());
  MPI_File_write(file, char_buf, buf.str().size(), MPI_CHAR, MPI_STATUS_IGNORE);
  delete [] char_buf;

  MPI_File_close(&file);
}

int CommMPI::GlobalRank() const
{
  int rank;

  MPI_Comm_rank(comm_global, &rank);

  return rank;
}

Index CommMPI::GlobalPos() const
{
  Index dims, periods, coords;

  MPI_Cart_get(comm_global, 3, dims.vec(), periods.vec(), coords.vec());

  return coords;
}

Index CommMPI::GlobalProcs() const
{
  Index dims, periods, coords;

  MPI_Cart_get(comm_global, 3, dims.vec(), periods.vec(), coords.vec());

  return dims;
}

void CommMPI::InitCommMPI(const MPI_Comm& comm)
{
  int status, size, rank;
  int dims[3] = {0, 0, 0};
  int periods[3];

  global_coarse_grid = NULL;

  for (int i=0; i<3; ++i)
    periods[i] = (BoundaryConditions()[i] == Periodic ? 1 : 0);

  MPI_Comm_size(comm, &size);
  MPI_Comm_rank(comm, &rank);

  MPI_Topo_test(comm, &status);

  if (status == MPI_CART) {

    comm_global = comm;

  }else {

    MPI_Dims_create(size, 3, dims);

    MPI_Cart_create(comm, 3, dims, periods, 1, &comm_global);

  }

  comm_info.SetCommGlobal(comm_global);

  MPI_Info_create(&info);
  char key[] = "no_locks";
  char val[] = "true";
  MPI_Info_set(info, key, val);

}

CommMPI::~CommMPI()
{
  std::map<const Grid*, TempGrid*>::iterator iter;

  for (iter=finer_grids.begin(); iter!=finer_grids.end(); ++iter)
    delete iter->second;

  for (iter=coarser_grids.begin(); iter!=coarser_grids.end(); ++iter)
    delete iter->second;

  MPI_Comm_free(&comm_global);

  if (win_created)
    MPI_Win_free(&win);

  MPI_Info_free(&info);

  delete global_coarse_grid;
}

std::string CommMPI::CreateOutputDirectory()
{
#ifdef HAVE_BOOST_FILESYSTEM
  std::string path, unique_path;
  std::stringstream unique_suffix;
  int suffix_counter = 0;
  char buffer[129];
  time_t rawtime;
  struct tm *timeinfo;
  int path_size;
  char* path_array;

  if (GlobalRank() == 0) {

    time(&rawtime);
    timeinfo = localtime(&rawtime);
    strftime(buffer, 128, "./output/%Y_%m_%d_%H_%M_%S/", timeinfo);
    path = buffer;

    if (!fs::exists("output"))
      fs::create_directory("output");

    unique_path = path;

    while (fs::exists(unique_path.c_str())) {

      unique_suffix.str("");
      unique_suffix << "_" << suffix_counter++ << "/";

      unique_path = path;
      unique_path.replace(unique_path.size()-1, 1, unique_suffix.str());

    }

    fs::create_directory(unique_path.c_str());

    path_size = unique_path.size() + 1;
    path_array = Helper::GetCharArray(unique_path);

    MPI_Bcast(&path_size, 1, MPI_INT, 0, comm_global);

  }else {

    MPI_Bcast(&path_size, 1, MPI_INT, 0, comm_global);
    path_array = new char[path_size];

  }

  MPI_Bcast(path_array, path_size, MPI_CHAR, 0, comm_global);

  unique_path = path_array;

  delete [] path_array;

  return unique_path;

#else

  return "./";

#endif
}


#endif /* HAVE_MPI */