source: src/comm/comm.cpp@ d24c2f

/**
 * @file   comm.hpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Wed Jun 16 13:21:06 2010
 *
 * @brief  Base class for communication.
 *
 */

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

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

#include "base/discretization.hpp"
#include "base/helper.hpp"
#include "base/stencil.hpp"
#include "comm/comm.hpp"
#include "comm/domain_decomposition.hpp"
#include "grid/grid.hpp"
#include "mg.hpp"

using namespace VMG;

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

  Grid::iterator grid_iter;
  Stencil::iterator stencil_iter;
  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 (grid_iter=rhs().Iterators().Local().Begin(); grid_iter!=rhs().Iterators().Local().End(); ++grid_iter) {
    val = rhs().GetVal(*grid_iter) - prefactor * A.GetDiag() * sol().GetVal(*grid_iter);
    for (stencil_iter=A.begin(); stencil_iter!=A.end(); ++stencil_iter)
      val -= prefactor * stencil_iter->Val() * sol().GetVal(grid_iter->X() + stencil_iter->Disp().X(),
                                                            grid_iter->Y() + stencil_iter->Disp().Y(),
                                                            grid_iter->Z() + stencil_iter->Disp().Z());
    norm += val*val;
  }

  norm = sqrt( Helper::pow_3(sol().MeshWidth()) * GlobalSum(norm) );

  return norm;
}

Grid& Comm::GetParticleGrid()
{
  if (particle_grid != NULL)
    return *particle_grid;

  const Multigrid& multigrid = *MG::GetRhs();
  const Grid& grid = multigrid(multigrid.GlobalMaxLevel());
  LocalIndices local = grid.Local();

  const int& near_field_cells = MG::GetFactory().GetObjectStorageVal<int>("PARTICLE_NEAR_FIELD_CELLS");

  // We don't need a boundary on this grid
  local.End() -= local.Begin();
  local.Begin() = 0;
  local.BoundaryBegin1() = local.BoundaryEnd1() = 0;
  local.BoundaryBegin2() = local.BoundaryEnd2() = 0;

  // Set grid size of intermediate temporary grid
  for (int i=0; i<3; ++i) {

    if (local.HasHalo1()[i]) {
      local.HaloBegin1()[i] = 0;
      local.HaloEnd1()[i] = local.Begin()[i] = near_field_cells+1;
      local.End()[i] = local.Begin()[i] + local.Size()[i];
    }

    if (local.HasHalo2()[i]) {
      local.HaloBegin2()[i] = local.End()[i];
      local.HaloEnd2()[i] = local.HaloBegin2()[i] + near_field_cells+1;
    }

  }

  local.SizeTotal() = local.Size() +
                      local.HaloEnd1() - local.HaloBegin1() +
                      local.HaloEnd2() - local.HaloBegin2();

  particle_grid = new Grid(grid.Global(), local, grid.Extent());

  return *particle_grid;
}

Comm::~Comm()
{
  delete dd;
  delete particle_grid;
}

const std::string& Comm::OutputPath()
{
  if (!output_directory_is_created) {
    output_path_str = CreateOutputDirectory();
    output_directory_is_created = true;
  }

  return output_path_str;
}