source: src/comm/comm_serial.cpp@ 66f24d

/**
 * @file   comm_serial.cpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Mon Apr 18 12:28:12 2011
 *
 * @brief  VMG::CommSerial
 *
 */

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

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

#include <iostream>
#include <fstream>
#include <cassert>
#include <cmath>

#include "base/discretization.hpp"
#include "base/stencil.hpp"
#include "base/vector.hpp"
#include "comm/comm_serial.hpp"
#include "comm/comm.hpp"
#include "grid/multigrid.hpp"
#include "grid/tempgrid.hpp"
#include "mg.hpp"

using namespace VMG;

inline vmg_float pow_3(vmg_float val)
{
  return val*val*val;
}

void CommSerial::Init()
{
  output_count = 0;
  error_norm_count = 0;
  residual_count = 0;

  /* Create directory to output defects */
#ifdef HAVE_BOOST_FILESYSTEM
  char buffer[129];
  time_t rawtime;
  struct tm *timeinfo;

  time(&rawtime);
  timeinfo = localtime(&rawtime);
  strftime(buffer, 128, "output/%Y_%m_%d_%H_%M_%S/", timeinfo);
  defects_path_str = buffer;
#else
  defects_path_str = "";
#endif
}

void CommSerial::CreateOutputDirectory()
{
#ifdef HAVE_BOOST_FILESYSTEM
  if (!fs::exists("output"))
    fs::create_directory("output");

  if (!fs::exists(defects_path_str.c_str()))
    fs::create_directory(defects_path_str.c_str());
#endif
}

Grid& CommSerial::GetFinerGrid(Multigrid& multigrid)
{
  return multigrid(multigrid.Level()+1);
}

Grid& CommSerial::GetCoarserGrid(Multigrid& multigrid)
{
  return multigrid(multigrid.Level()-1);
}

vmg_float CommSerial::ComputeResidualNorm(Multigrid& sol, Multigrid& rhs)
{
#ifdef DEBUG
  sol().IsCompatible(rhs());
  sol().IsConsistent();
  rhs().IsConsistent();
#endif

  Stencil::iterator iter;
  Index i;
  vmg_float val;
  vmg_float norm = 0.0;

  const vmg_float prefactor = MG::GetDiscretization()->OperatorPrefactor(sol());
  const Stencil& A = MG::GetDiscretization()->GetStencil();

  this->CommToGhosts(sol());

  if (sol().Global().BoundaryType() == LocallyRefined)
    MG::GetDiscretization()->SetInnerBoundary(sol(), rhs(), sol(sol.Level()-1));

  for (i.X()=rhs().Local().Begin().X(); i.X()<rhs().Local().End().X(); ++i.X())
    for (i.Y()=rhs().Local().Begin().Y(); i.Y()<rhs().Local().End().Y(); ++i.Y())
      for (i.Z()=rhs().Local().Begin().Z(); i.Z()<rhs().Local().End().Z(); ++i.Z()) {

        val = rhs().GetVal(i) - prefactor * A.GetDiag() * sol().GetVal(i);

        for (iter = A.begin(); iter != A.end(); iter++)
          val -= prefactor * iter->val * sol().GetVal(i+iter);

        norm += val*val;

      }

  norm = sqrt( pow_3(sol().MeshWidth()) * norm );

#ifdef DEBUG
  std::ofstream out;
  char path_str[129];

  CreateOutputDirectory();

  sprintf(path_str, "%sresidual.txt", defects_path_str.c_str());

  out.open(path_str, std::ios::app);

  assert(out.good());

  out << residual_count++ << " " << norm << std::endl;

  out.close();
#endif

  return norm;
}

void CommSerial::CommFromGhosts(Grid& mesh)
{
  if (this->BoundaryConditions() == Periodic) {

    // Copy values in x-direction
    for (int i=mesh.Local().HaloBegin1().X(); i<mesh.Local().HaloEnd1().X(); i++)
      for (int j=0; j<mesh.Local().SizeTotal().Y(); j++)
        for (int k=0; k<mesh.Local().SizeTotal().Z(); k++) {
          assert(i+mesh.Local().Size().X() >= mesh.Local().Begin().X() && i+mesh.Local().Size().X() < mesh.Local().End().X());
          mesh(i+mesh.Local().Size().X(),j,k) += mesh.GetVal(i,j,k);
        }

    for (int i=mesh.Local().HaloBegin2().X(); i<mesh.Local().HaloEnd2().X(); i++)
      for (int j=0; j<mesh.Local().SizeTotal().Y(); j++)
        for (int k=0; k<mesh.Local().SizeTotal().Z(); k++) {
          assert(i-mesh.Local().Size().X() >= mesh.Local().Begin().X() && i-mesh.Local().Size().X() < mesh.Local().End().X());
          mesh(i-mesh.Local().Size().X(),j,k) += mesh.GetVal(i,j,k);
        }

    // Copy values in y-direction
    for (int i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
      for (int j=mesh.Local().HaloBegin1().Y(); j<mesh.Local().HaloEnd1().Y(); j++)
        for (int k=0; k<mesh.Local().SizeTotal().Z(); k++) {

          mesh(i,j+mesh.Local().Size().Y(),k) += mesh.GetVal(i,j,k);
        }
    for (int i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
      for (int j=mesh.Local().HaloBegin2().Y(); j<mesh.Local().HaloEnd2().Y(); j++)
        for (int k=0; k<mesh.Local().SizeTotal().Z(); k++)
          mesh(i,j-mesh.Local().Size().Y(),k) += mesh.GetVal(i,j,k);

    // Copy values in z-direction
    for (int i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
      for (int j=mesh.Local().Begin().Y(); j<mesh.Local().End().Y(); j++)
        for (int k=mesh.Local().HaloBegin1().Z(); k<mesh.Local().HaloEnd1().Z(); k++)
          mesh(i,j,k+mesh.Local().Size().Z()) += mesh.GetVal(i,j,k);

    for (int i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
      for (int j=mesh.Local().Begin().Y(); j<mesh.Local().End().Y(); j++)
        for (int k=mesh.Local().HaloBegin2().Z(); k<mesh.Local().HaloEnd2().Z(); k++)
          mesh(i,j,k-mesh.Local().Size().Z()) += mesh.GetVal(i,j,k);
  }

#ifdef DEBUG
  mesh.IsConsistent();
#endif
}

void CommSerial::CommToGhosts(Grid& mesh)
{
#ifdef DEBUG
  mesh.IsConsistent();
#endif

  int i,j,k;

  if (this->BoundaryConditions() == Periodic) {

    // Copy values in z-direction
    for (i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
      for (j=mesh.Local().Begin().Y(); j<mesh.Local().End().Y(); j++)
        for (k=mesh.Local().HaloBegin1().Z(); k<mesh.Local().HaloEnd1().Z(); k++)
          mesh(i,j,k) = mesh.GetVal(i,j,k+mesh.Local().Size().Z());

    for (i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
      for (j=mesh.Local().Begin().Y(); j<mesh.Local().End().Y(); j++)
        for (k=mesh.Local().HaloBegin2().Z(); k<mesh.Local().HaloEnd2().Z(); k++)
          mesh(i,j,k) = mesh.GetVal(i,j,k-mesh.Local().Size().Z());

    // Copy values in y-direction
    for (i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
      for (j=mesh.Local().HaloBegin1().Y(); j<mesh.Local().HaloEnd1().Y(); j++)
        for (k=0; k<mesh.Local().SizeTotal().Z(); k++)
          mesh(i,j,k) = mesh.GetVal(i,j+mesh.Local().Size().Y(),k);

    for (i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
      for (j=mesh.Local().HaloBegin2().Y(); j<mesh.Local().HaloEnd2().Y(); j++)
        for (k=0; k<mesh.Local().SizeTotal().Z(); k++)
          mesh(i,j,k) = mesh.GetVal(i,j-mesh.Local().Size().Y(),k);

    // Copy values in x-direction
    for (i=mesh.Local().HaloBegin1().X(); i<mesh.Local().HaloEnd1().X(); i++)
      for (j=0; j<mesh.Local().SizeTotal().Y(); j++)
        for (k=0; k<mesh.Local().SizeTotal().Z(); k++)
          mesh(i,j,k) = mesh.GetVal(i+mesh.Local().Size().X(),j,k);

    for (i=mesh.Local().HaloBegin2().X(); i<mesh.Local().HaloEnd2().X(); i++)
      for (j=0; j<mesh.Local().SizeTotal().Y(); j++)
        for (k=0; k<mesh.Local().SizeTotal().Z(); k++)
          mesh(i,j,k) = mesh.GetVal(i-mesh.Local().Size().X(),j,k);
  }

#ifdef DEBUG
  mesh.IsConsistent();
#endif
}

void CommSerial::OpenFileAndPrintHeader(std::ofstream& out, const Grid& mesh, const char* information, bool inner)
{
  char file_str[129], path_str[129];
  struct tm *timeinfo;
  time_t rawtime;

  CreateOutputDirectory();

  output_count++;

  time(&rawtime);
  timeinfo = localtime(&rawtime);

  strftime(file_str, 128, "%H_%M_%S", timeinfo);

  sprintf(path_str, "%s%04d_%s.dat", defects_path_str.c_str(), output_count, file_str);

  out.open(path_str, std::ios::trunc);

  if (!out.good()) {
    printf("Multigrid: File %s not accessible.\n", path_str);
    return;
  }

  out << "# vtk DataFile Version 2.0" << std::endl
      << output_count << ": " << information << std::endl
      << "ASCII" << std::endl
      << "DATASET STRUCTURED_POINTS" << std::endl;

  if (inner) {
    out << "DIMENSIONS "
        << mesh.Local().Size().X() << " "
        << mesh.Local().Size().Y() << " "
        << mesh.Local().Size().Z() << std::endl;
  }else {
    out << "DIMENSIONS "
        << mesh.Local().SizeTotal().X() << " "
        << mesh.Local().SizeTotal().Y() << " "
        << mesh.Local().SizeTotal().Z() << std::endl;
  }

  out   << "ORIGIN 0 0 0" << std::endl
        << "SPACING " << mesh.MeshWidth() << " "
        << mesh.MeshWidth() << " "
        << mesh.MeshWidth() << std::endl;

  if (inner)
    out << "POINT_DATA " << mesh.Local().Size().Product() << std::endl;
  else
    out << "POINT_DATA " << mesh.Local().SizeTotal().Product() << std::endl;

  out << "SCALARS residual double 1" << std::endl
      << "LOOKUP_TABLE default" << std::endl;
}

void CommSerial::PrintDefect(const Grid& sol, const Grid& rhs, const char* information)
{
  std::ofstream out;

  TempGrid *temp = MG::GetTempGrid();
  temp->SetProperties(sol);
  temp->ImportFromResidual(sol, rhs);

#ifdef DEBUG
  temp->IsConsistent();
#endif

  OpenFileAndPrintHeader(out, *temp, information, true);

  assert(out.good());

  if (!out.good())
    return;

  for (int k=temp->Local().Begin().Z(); k<temp->Local().End().Z(); k++)
    for (int j=temp->Local().Begin().Y(); j<temp->Local().End().Y(); j++)
      for (int i=temp->Local().Begin().X(); i<temp->Local().End().X(); i++)
        out << std::scientific << temp->GetVal(i,j,k) << std::endl;

  out.close();
}

void CommSerial::PrintGrid(const Grid& mesh, const char* information)
{
  std::ofstream out;

  OpenFileAndPrintHeader(out, mesh, information, false);

  for (int k=0; k<mesh.Local().SizeTotal().Z(); k++)
    for (int j=0; j<mesh.Local().SizeTotal().Y(); j++)
      for (int i=0; i<mesh.Local().SizeTotal().X(); i++)
        out << std::scientific << mesh.GetVal(i,j,k) << std::endl;

  out.close();
}

void CommSerial::PrintCorrectedGrid(const Grid& mesh, const char* information)
{
  std::ofstream out;

  OpenFileAndPrintHeader(out, mesh, information, false);

  for (int k=0; k<mesh.Local().SizeTotal().Z(); k++)
    for (int j=0; j<mesh.Local().SizeTotal().Y(); j++)
      for (int i=0; i<mesh.Local().SizeTotal().X(); i++)
        out << std::scientific << mesh.GetCorrectedVal(i,j,k) << std::endl;

  out.close();
}

void CommSerial::PrintInnerGrid(const Grid& mesh, const char* information)
{
  std::ofstream out;

  OpenFileAndPrintHeader(out, mesh, information, true);

  for (int k=mesh.Local().Begin().Z(); k<mesh.Local().End().Z(); k++)
    for (int j=mesh.Local().Begin().Y(); j<mesh.Local().End().Y(); j++)
      for (int i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
        out << std::scientific << mesh.GetVal(i,j,k) << std::endl;

  out.close();
}

void CommSerial::PrintInnerCorrectedGrid(const Grid& mesh, const char* information)
{
  std::ofstream out;

  OpenFileAndPrintHeader(out, mesh, information, true);

  for (int k=mesh.Local().Begin().Z(); k<mesh.Local().End().Z(); k++)
    for (int j=mesh.Local().Begin().Y(); j<mesh.Local().End().Y(); j++)
      for (int i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
        out << std::scientific << mesh.GetCorrectedVal(i,j,k) << std::endl;

  out.close();
}

void CommSerial::DebugPrintError(const Grid& sol, const char* information)
{
  Index i;
  std::ofstream out;

  OpenFileAndPrintHeader(out, sol, information, false);

  for (i.Z()=0; i.Z()<sol.Local().SizeTotal().Z(); ++i.Z())
    for (i.Y()=0; i.Y()<sol.Local().SizeTotal().Y(); ++i.Y())
      for (i.X()=0; i.X()<sol.Local().SizeTotal().X(); ++i.X()) {
        out << std::scientific << sol.DebugKnownSolution(i) - sol.GetVal(i) << std::endl;

      }
}

inline int GetIndex(const Grid& grid, int i, int j, int k)
{
  if (grid.Global().BoundaryType() == LocallyRefined)
    return k + grid.Local().Size().Z() * (j + grid.Local().Size().Y() * i);
  else
    return k + grid.Local().SizeTotal().Z() * (j + grid.Local().SizeTotal().Y() * i);
}

void CommSerial::PrintGridStructureLevel(Grid& grid, std::ofstream& out)
{
  const vmg_float sp = grid.MeshWidth();
  int numLines;

  if (grid.Global().BoundaryType() == LocallyRefined)
    numLines = grid.Local().Size().X() * grid.Local().Size().Y() +
      grid.Local().Size().Y() * grid.Local().Size().Z() +
      grid.Local().Size().X() * grid.Local().Size().Z();
  else
    numLines = grid.Local().SizeTotal().X() * grid.Local().SizeTotal().Y() +
      grid.Local().SizeTotal().Y() * grid.Local().SizeTotal().Z() +
      grid.Local().SizeTotal().X() * grid.Local().SizeTotal().Z();

  out << "    <Piece";

  if (grid.Global().BoundaryType() == LocallyRefined) {
    out << " NumberOfPoints=\"" << grid.Local().Size().Product()  << "\""
        << " NumberOfVerts=\"" << grid.Local().Size().Product()  << "\"";
  }else {
    out << " NumberOfPoints=\"" << grid.Local().SizeTotal().Product()  << "\""
        << " NumberOfVerts=\"" << grid.Local().SizeTotal().Product()  << "\"";
  }

  out << " NumberOfLines=\"" << numLines << "\""
      << " NumberOfStrips=\"0\""
      << " NumberOfPolys=\"0\"" << ">" << std::endl
      << "      <Points>" << std::endl
      << "        <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl
      << "          ";

  if (grid.Global().BoundaryType() == LocallyRefined) {
    for (int i=0; i<grid.Local().Size().X(); i++)
      for (int j=0; j<grid.Local().Size().Y(); j++)
        for (int k=0; k<grid.Local().Size().Z(); k++)
          out << grid.Extent().Begin().X() + i * sp << " "
              << grid.Extent().Begin().Y() + j * sp << " "
              << grid.Extent().Begin().Z() + k * sp << " ";
  }else {
    for (int i=0; i<grid.Local().SizeTotal().X(); i++)
      for (int j=0; j<grid.Local().SizeTotal().Y(); j++)
        for (int k=0; k<grid.Local().SizeTotal().Z(); k++)
          out << grid.Extent().Begin().X() + i * sp << " "
              << grid.Extent().Begin().Y() + j * sp << " "
              << grid.Extent().Begin().Z() + k * sp << " ";
  }

  out << std::endl
      << "        </DataArray>" << std::endl
      << "      </Points>" << std::endl
      << "      <Verts>" << std::endl
      << "        <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << std::endl
      << "          ";

  if (grid.Global().BoundaryType() == LocallyRefined) {
    for (int i=0; i<grid.Local().Size().Product(); i++)
      out << i << " ";
  }else {
    for (int i=0; i<grid.Local().SizeTotal().Product(); i++)
      out << i << " ";
  }

  out << std::endl
      << "        </DataArray>" << std::endl
      << "        <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << std::endl
      << "          ";

  if (grid.Global().BoundaryType() == LocallyRefined) {
    for (int i=1; i<=grid.Local().Size().Product(); i++)
      out  << i << " ";
  }else {
    for (int i=1; i<=grid.Local().SizeTotal().Product(); i++)
      out  << i << " ";
  }

  out << std::endl
      << "        </DataArray>" << std::endl
      << "      </Verts>" << std::endl
      << "      <Lines>" << std::endl
      << "        <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << std::endl
      << "          ";

  if (grid.Global().BoundaryType() == LocallyRefined) {
    for (int i=0; i<grid.Local().Size().X(); i++)
      for (int j=0; j<grid.Local().Size().Y(); j++)
        out << GetIndex(grid,i,j,0) << " "
            << GetIndex(grid,i,j,grid.Local().Size().Z()-1) << " ";

    for (int j=0; j<grid.Local().Size().Y(); j++)
      for (int k=0; k<grid.Local().Size().Z(); k++)
        out << GetIndex(grid,0,j,k) << " "
            << GetIndex(grid,grid.Local().Size().X()-1,j,k) << " ";

    for (int i=0; i<grid.Local().Size().X(); i++)
      for (int k=0; k<grid.Local().Size().Z(); k++)
        out << GetIndex(grid,i,0,k) << " "
            << GetIndex(grid,i,grid.Local().Size().Y()-1,k) << " ";
  }else {
    for (int i=0; i<grid.Local().SizeTotal().X(); i++)
      for (int j=0; j<grid.Local().SizeTotal().Y(); j++)
        out << GetIndex(grid,i,j,0) << " "
            << GetIndex(grid,i,j,grid.Local().SizeTotal().Z()-1) << " ";

    for (int j=0; j<grid.Local().SizeTotal().Y(); j++)
      for (int k=0; k<grid.Local().SizeTotal().Z(); k++)
        out << GetIndex(grid,0,j,k) << " "
            << GetIndex(grid,grid.Local().SizeTotal().X()-1,j,k) << " ";

    for (int i=0; i<grid.Local().SizeTotal().X(); i++)
      for (int k=0; k<grid.Local().SizeTotal().Z(); k++)
        out << GetIndex(grid,i,0,k) << " "
            << GetIndex(grid,i,grid.Local().SizeTotal().Y()-1,k) << " ";
  }

  out << std::endl
      << "        </DataArray>" << std::endl
      << "        <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << std::endl
      << "          ";

  for (int i=1; i<=numLines; i++)
    out << 2*i << " ";

  out << std::endl
      << "        </DataArray>" << std::endl
      << "      </Lines>" << std::endl
      << "    </Piece>" << std::endl;
}

void CommSerial::DebugPrintGridStructure(Multigrid& grid)
{
  std::ofstream out;
  char path_str[129];

  CreateOutputDirectory();

  sprintf(path_str, "%sgrid.vtp", defects_path_str.c_str());

  out.open(path_str, std::ios::trunc);

  if (!out.good()) {
    printf("Multigrid: File %s not accessible.\n", path_str);
    return;
  }

  out << "<?xml version=\"1.0\"?>" << std::endl
      << "<VTKFile type=\"PolyData\">" << std::endl
      << "  <PolyData>" << std::endl;

  for (int i=grid.MinLevel(); i<=grid.MaxLevel(); i++)
    PrintGridStructureLevel(grid(i), out);

  out << "  </PolyData>" << std::endl
      << "</VTKFile>" << std::endl;

  out.close();
}