/**
 * @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/tuple.hpp"
#include "comm/comm_mpi.hpp"
#include "grid/grid.hpp"
#include "grid/multigrid.hpp"
#include "grid/tempgrid.hpp"
#include "mg.hpp"
#include "base/timer.hpp"
static char print_buffer[512];

using namespace VMG;

void CommMPI::Isend(Grid& grid, VMG::MPI::Datatype& type, const MPI_Comm& comm, const int& tag)
{
  if (type.Feasible())
    MPI_Isend(&grid(0), 1, type.Type(), type.Rank(), tag, comm, &Request());
}

void CommMPI::IsendAll(Grid& grid, std::vector<VMG::MPI::Datatype>& types, const MPI_Comm& comm, const int& tag_start)
{
  for (unsigned int i=0; i<types.size(); ++i)
    if (types[i].Feasible())
      MPI_Isend(&grid(0), 1, types[i].Type(), types[i].Rank(), tag_start+types[i].TagSend(), comm, &Request());
}

void CommMPI::Irecv(Grid& grid, VMG::MPI::Datatype& type, const MPI_Comm& comm, const int& tag)
{
  if (type.Feasible())
    MPI_Irecv(&grid(0), 1, type.Type(), type.Rank(), tag, comm, &Request());
}

void CommMPI::IrecvAll(Grid& grid, std::vector<VMG::MPI::Datatype>& types, const MPI_Comm& comm, const int& tag_start)
{
  for (unsigned int i=0; i<types.size(); ++i)
    if (types[i].Feasible())
      MPI_Irecv(&grid(0), 1, types[i].Type(), types[i].Rank(), tag_start+types[i].TagReceive(), comm, &Request());
}

void CommMPI::IrecvAllBuffered(std::vector<VMG::MPI::Datatype>& types, const MPI_Comm& comm, const int& tag_start)
{
  receive_buffer.resize(types.size());
  for (unsigned int i=0; i<types.size(); ++i)
    if (types[i].Feasible()) {
      receive_buffer[i].resize(types[i].Size());
      MPI_Irecv(&receive_buffer[i].front(), types[i].Size(), MPI_DOUBLE,
		types[i].Rank(), tag_start+types[i].TagReceive(), comm, &Request());
    }
}

void CommMPI::AddBufferAll(Grid& grid, std::vector<VMG::MPI::Datatype>& types)
{
  Index ind;

  for (unsigned int i=0; i<types.size(); ++i) {
    int c = 0;
    for (ind.X()=types[i].Starts().X(); ind.X()<types[i].Starts().X()+types[i].Subsizes().X(); ++ind.X())
      for (ind.Y()=types[i].Starts().Y(); ind.Y()<types[i].Starts().Y()+types[i].Subsizes().Y(); ++ind.Y())
  	for (ind.Z()=types[i].Starts().Z(); ind.Z()<types[i].Starts().Z()+types[i].Subsizes().Z(); ++ind.Z())
  	  grid(ind) += receive_buffer[i][c++];
  }

  receive_buffer.clear();
}

void CommMPI::CommSubgrid(Grid& grid_old, Grid& grid_new)
{
  MPI_Comm comm = settings.CommunicatorGlobal(grid_old);
  if (comm != MPI_COMM_NULL) {
    VMG::MPI::DatatypesGlobal& datatypes = settings.DatatypesGlobal(grid_old, grid_new);
    IsendAll(grid_old, datatypes.Send(), comm, 0411);
    IrecvAll(grid_new, datatypes.Receive(), comm, 0411);
    ClearRequests();
  }
}

void CommMPI::CommAddSubgrid(Grid& grid_old, Grid& grid_new)
{
  MPI_Comm comm = settings.CommunicatorGlobal(grid_old);
  if (comm != MPI_COMM_NULL) {
    VMG::MPI::DatatypesGlobal& datatypes = settings.DatatypesGlobal(grid_old, grid_new);
    IsendAll(grid_old, datatypes.Send(), comm, 1806);
    IrecvAllBuffered(datatypes.Receive(), comm, 1806);
    ClearRequests();
    AddBufferAll(grid_new, datatypes.Receive());
  }
}

void CommMPI::CommFromGhosts(Grid& grid)
{
  MPI_Comm comm = settings.CommunicatorLocal(grid);

  if (comm != MPI_COMM_NULL) {
    VMG::MPI::DatatypesLocal& types = settings.DatatypesLocal(grid);
    IsendAll(grid, types.Halo(), comm, 1337);
    IrecvAllBuffered(types.NB(), comm, 1337);
    ClearRequests();
    AddBufferAll(grid, types.NB());
  }
}

void CommMPI::CommFromGhostsAsyncStart(Grid& grid)
{
 MPI_Comm comm = settings.CommunicatorLocal(grid);

  if (comm != MPI_COMM_NULL) {
    VMG::MPI::DatatypesLocal& types = settings.DatatypesLocal(grid);
    IsendAll(grid, types.Halo(), comm, 0xc0ffee);
    IrecvAllBuffered(types.NB(), comm, 0xc0ffee);
  }
}

void CommMPI::CommFromGhostsAsyncFinish(Grid& grid)
{
  VMG::MPI::DatatypesLocal& types = settings.DatatypesLocal(grid);
  ClearRequests();
  AddBufferAll(grid, types.NB());
}

void CommMPI::CommToGhosts(Grid& grid)
{
  MPI_Comm comm = settings.CommunicatorLocal(grid);

  if (comm != MPI_COMM_NULL) {
    VMG::MPI::DatatypesLocal& types = settings.DatatypesLocal(grid);
    IsendAll(grid, types.NB(), comm, 2310);
    IrecvAll(grid, types.Halo(), comm, 2310);
    ClearRequests();
  }
}

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

  if (comm != MPI_COMM_NULL) {

    Index i;
    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++;

    }

  }
}

void CommMPI::CommToGhostsAsyncStart(Grid& grid)
{
  MPI_Comm comm = settings.CommunicatorLocal(grid);

  if (comm != MPI_COMM_NULL) {
    Timer::Start("CTGAsyncGetDatatypes");
    VMG::MPI::DatatypesLocal& types = settings.DatatypesLocal(grid);
    Timer::Stop("CTGAsyncGetDatatypes");
    Timer::Start("CTGAsyncIsend");
    IsendAll(grid, types.NB(), comm, 2412);
    Timer::Stop("CTGAsyncIsend");
    Timer::Start("CTGAsyncIrecv");
    IrecvAll(grid, types.Halo(), comm, 2412);
    Timer::Stop("CTGAsyncIrecv");
  }
}

void CommMPI::CommToGhostsAsyncFinish()
{
  ClearRequests();
}

void CommMPI::CommParticles(const Grid& grid, std::list<Particle::Particle>& particles)
{
  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;
  std::vector<int> global_extent, send_sizes, recv_sizes;
  std::vector<Index> begin_remote, end_remote;
  std::vector< std::vector<int> > pos_indices, charge_indices;
  std::vector<vmg_float> buffer_x, buffer_q;
  std::vector<vmg_int> buffer_ind;
  MPI_Datatype dt_temp;

  global_extent.resize(6*size);
  send_sizes.resize(size);
  recv_sizes.resize(size);
  begin_remote.resize(size);
  end_remote.resize(size);
  pos_indices.resize(size);
  charge_indices.resize(size);

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

  MPI_Allgather(MPI_IN_PLACE, 6, MPI_INT, &global_extent.front(), 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.front(), 1, MPI_INT, &recv_sizes.front(), 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();
}

void CommMPI::CommLCListGhosts(const Grid& grid, Particle::LinkedCellList& lc)
{
  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();
  }
}

void CommMPI::CommAddPotential(std::list<Particle::Particle>& particles)
{
  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);
}

void CommMPI::CommAddPotential(Particle::LinkedCellList& lc)
{
  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);
}

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 = settings.CommunicatorGlobal((*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 = settings.CommunicatorGlobal(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];

  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);

  }

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

}

CommMPI::~CommMPI()
{
  MPI_Comm_free(&comm_global);

  if (win_created)
    MPI_Win_free(&win);

  MPI_Info_free(&info);
}

Grid& CommMPI::GetCoarserGrid(Multigrid& multigrid)
{
  return settings.CoarserGrid(multigrid(multigrid.Level()));
}

Grid& CommMPI::GetFinerGrid(Multigrid& multigrid)
{
  return settings.FinerGrid(multigrid(multigrid.Level()));
}

Grid& CommMPI::GetGlobalCoarseGrid(Multigrid& multigrid)
{
  Grid& global_coarse_grid = settings.GlobalCoarseGrid();
  CommSubgrid(multigrid(multigrid.MinLevel()), global_coarse_grid);
  return 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 */