source: src/interface/interface_particles.cpp@ dfed1c

/**
 * @file   interface_particles.cpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Mon Apr 18 12:56:48 2011
 *
 * @brief  VMG::InterfaceParticles
 *
 */

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

#ifdef HAVE_MPI
#include <mpi.h>
#endif

#include <cmath>
#include <cstring>

#include "base/helper.hpp"
#include "base/index.hpp"
#include "base/linked_cell_list.hpp"
#include "base/timer.hpp"
#include "base/vector.hpp"
#include "comm/comm.hpp"
#include "grid/grid.hpp"
#include "grid/multigrid.hpp"
#include "grid/tempgrid.hpp"
#include "interface/interface_particles.hpp"
#include "mg.hpp"

using namespace VMG;

void InterfaceParticles::ImportRightHandSide(Multigrid& multigrid)
{
  Timer::Start("Particle/Grid Interpolation");

  Index index_global, index_local, index;
  Vector pos_rel, pos_abs, grid_val;
  vmg_float length;

  Factory& factory = MG::GetFactory();

  const int& near_field_cells = factory.GetObjectStorageVal<int>("PARTICLE_NEAR_FIELD_CELLS");

  TempGrid* temp_grid = GetTempGrid(0);
  Grid& grid = multigrid(multigrid.GlobalMaxLevel());
  LocalIndices local = grid.Local();

  //TODO This has to be replaced for a VMG::Vector for different mesh widths on the axis
  const vmg_float r_cut = (near_field_cells+0.5) * grid.Extent().MeshWidth().Max();

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

  temp_grid->SetProperties(grid.Global(), local, grid.Extent());
  temp_grid->Clear();

  /*
   * Distribute particles to their processes
   */
  particles.clear();
  MG::GetComm()->CommParticles(grid, particles);

  /*
   * Charge assignment on the grid
   */
  std::list<Particle::Particle>::iterator iter;

#ifdef DEBUG_OUTPUT
  vmg_float particle_charges = 0.0;
  for (iter=particles.begin(); iter!=particles.end(); ++iter)
    particle_charges += iter->Charge();
  particle_charges = MG::GetComm()->GlobalSumRoot(particle_charges);
  MG::GetComm()->PrintStringOnce("Particle list charge sum: %e", particle_charges);
#endif

  for (iter=particles.begin(); iter!=particles.end(); ++iter) {

    // Compute global and local grid index of the lower left corner of the grid cell containing the particle
    index_global = static_cast<Index>((iter->Pos() - temp_grid->Extent().Begin()) / temp_grid->Extent().MeshWidth());
    index_local = index_global - temp_grid->Global().BeginLocal() + temp_grid->Local().Begin();

    // Iterate over all grid points which lie in the support of the interpolating B-Spline
    for (index.X()=-1*near_field_cells-1; index.X()<=near_field_cells+1; ++index.X())
      for (index.Y()=-1*near_field_cells-1; index.Y()<=near_field_cells+1; ++index.Y())
        for (index.Z()=-1*near_field_cells-1; index.Z()<=near_field_cells+1; ++index.Z()) {

          // Compute distance from grid point to particle
          length = ( iter->Pos() - grid.Extent().Begin() - grid.Extent().MeshWidth() * (index_global+index) ).Length();

          // If grid point lies in the support of the B-Spline, do charge interpolation
          if (length < r_cut)
            (*temp_grid)(index_local+index) += spl.EvaluateSpline(length) * iter->Charge();

        }

  }

  // Communicate charges over halo
  MG::GetComm()->CommFromGhosts(*temp_grid);

  // Assign charge values to the right hand side
  for (index.X()=0; index.X()<grid.Local().Size().X(); ++index.X())
    for (index.Y()=0; index.Y()<grid.Local().Size().Y(); ++index.Y())
      for (index.Z()=0; index.Z()<grid.Local().Size().Z(); ++index.Z())
        grid(index + grid.Local().Begin()) = temp_grid->GetVal(index + temp_grid->Local().Begin());

  Timer::Stop("Particle/Grid Interpolation");

#ifdef DEBUG_OUTPUT
  Grid::iterator grid_iter;
  vmg_float charge_sum = 0.0;
  for (grid_iter=grid.Iterators().Local().Begin(); grid_iter!=grid.Iterators().Local().End(); ++grid_iter)
    charge_sum += grid.GetVal(*grid_iter);
  charge_sum = MG::GetComm()->GlobalSum(charge_sum);
  MG::GetComm()->PrintStringOnce("Grid charge sum: %e", charge_sum);
#endif
}

void InterfaceParticles::ExportSolution(Grid& grid)
{
  Timer::Start("Grid/Particle Interpolation");

  Index index;
  vmg_float length;
  Vector dist_vec;

#ifdef DEBUG_OUTPUT
  vmg_float e = 0.0;
  vmg_float e_long = 0.0;
  vmg_float e_self = 0.0;
#endif

  Factory& factory = MG::GetFactory();
  Comm* comm = MG::GetComm();

  /*
   * Get parameters and arrays
   */
  const vmg_int& num_particles_local = factory.GetObjectStorageVal<vmg_int>("PARTICLE_NUM_LOCAL");
  const int& near_field_cells = factory.GetObjectStorageVal<int>("PARTICLE_NEAR_FIELD_CELLS");
  vmg_float* p = factory.GetObjectStorageArray<vmg_float>("PARTICLE_POTENTIAL_ARRAY");

  const vmg_float r_cut = (near_field_cells+0.5) * grid.Extent().MeshWidth().Max();

  /*
   * Initialze potential
   */
  for (vmg_int i=0; i<num_particles_local; ++i)
    p[i] = 0.0;

  /*
   * Copy potential values to a grid with sufficiently large halo size.
   * This may be optimized in future.
   * The parameters of this grid have been set in the import step.
   */
  TempGrid* temp_grid = GetTempGrid(0);

  for (index.X()=0; index.X()<grid.Local().Size().X(); ++index.X())
    for (index.Y()=0; index.Y()<grid.Local().Size().Y(); ++index.Y())
      for (index.Z()=0; index.Z()<grid.Local().Size().Z(); ++index.Z())
        (*temp_grid)(index + temp_grid->Local().Begin()) = grid.GetVal(index + grid.Local().Begin());

  comm->CommToGhosts(*temp_grid);

  /*
   * Do potential back-interpolation
   */
  std::list<Particle::Particle>::iterator p_iter;
  for (p_iter=particles.begin(); p_iter!=particles.end(); ++p_iter) {

    // Interpolate long range part of potential
    p_iter->Pot() = 0.5 * (Helper::InterpolateTrilinear(p_iter->Pos(), *temp_grid) - p_iter->Charge() * spl.GetAntiDerivativeZero());

#ifdef DEBUG_OUTPUT
    e_long += 0.5 * p_iter->Charge() * Helper::InterpolateTrilinear(p_iter->Pos(), *temp_grid);
    e_self += 0.5 * p_iter->Charge() * p_iter->Charge() * spl.GetAntiDerivativeZero();
#endif

  }

  comm->CommAddPotential(particles);

  /*
   * Compute near field correction
   */
  Particle::LinkedCellList lc(particles, near_field_cells, grid);
  Particle::LinkedCellList::iterator p1, p2;
  Grid::iterator iter;

  comm->CommLCListGhosts(grid, lc);

  for (iter=lc.Iterators().Local().Begin(); iter!=lc.Iterators().Local().End(); ++iter)
    for (p1=lc(*iter).begin(); p1!=lc(*iter).end(); ++p1)
      for (index.X()=-1*near_field_cells-1; index.X()<=near_field_cells+1; ++index.X())
        for (index.Y()=-1*near_field_cells-1; index.Y()<=near_field_cells+1; ++index.Y())
          for (index.Z()=-1*near_field_cells-1; index.Z()<=near_field_cells+1; ++index.Z())
            for (p2=lc(*iter+index).begin(); p2!=lc(*iter+index).end(); ++p2)
              if (*p1 != *p2) {

                length = (p2->Pos() - p1->Pos()).Length();

                //TODO Rewrite this equation more efficiently
                if (length < r_cut)
                  p1->Pot() += 0.5 * p2->Charge() / length * (0.25*M_1_PI + spl.EvaluatePotential(length));

              }

  comm->CommAddPotential(lc);

  Timer::Stop("Grid/Particle Interpolation");

#ifdef DEBUG_OUTPUT
  vmg_float* q = factory.GetObjectStorageArray<vmg_float>("PARTICLE_CHARGE_ARRAY");

  for (int i=0; i<num_particles_local; ++i)
    e += p[i] * q[i];

  e = comm->GlobalSumRoot(e);
  e_long = comm->GlobalSumRoot(e_long);
  e_self = comm->GlobalSumRoot(e_self);

  comm->PrintStringOnce("E_long:  %e", e_long);
  comm->PrintStringOnce("E_short: %e", e - e_long + e_self);
  comm->PrintStringOnce("E_self:  %e", e_self);
  comm->PrintStringOnce("E_total: %e", e);

#endif /* DEBUG_OUTPUT */

}