source: src/interface/interface_particles_cf.cpp@ 728149

/**
 * @file   interface_particles_cf.cpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Fri Sep 16 13:24:48 2011
 *
 */

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

#ifdef HAVE_MPI
#include <mpi.h>
#ifdef HAVE_MARMOT
#include <enhancempicalls.h>
#include <sourceinfompicalls.h>
#endif
#endif

#include <cassert>
#include <cstring>
#include <cstdio>
#include <fstream>
#include <sstream>

#include "base/factory.hpp"
#include "base/helper.hpp"
#include "interface/interface_particles.hpp"
#include "interface/interface_particles_cf.hpp"
#include "level/level_operator_cs.hpp"
#include "level/level_operator_fas.hpp"
#include "samples/discretization_poisson_fd.hpp"
#include "samples/stencils.hpp"
#include "samples/techniques.hpp"
#include "solver/solver_regular.hpp"
#include "solver/solver_singular.hpp"
#include "smoother/gsrb.hpp"
#include "thirdparty/pugixml/pugixml.hpp"

#ifdef HAVE_MPI
#include "comm/comm_mpi.hpp"
#include "comm/domain_decomposition_mpi.hpp"
#else
#include "comm/comm_serial.hpp"
#endif

#ifdef HAVE_LAPACK
#include "solver/dgesv.hpp"
#else
#include "solver/givens.hpp"
#endif

using namespace VMG;

InterfaceParticlesCF::InterfaceParticlesCF(const char* filename)
{
  Factory& factory = MG::GetFactory();

  int max_level, max_iterations;
  int pre_smoothing_steps, post_smoothing_steps;
  int cycle_type, near_field_cells, discretization_order;
  vmg_float precision;
  std::string lop_str, datafile_str;

  pugi::xml_document doc;

#ifdef HAVE_MPI
  int rank;
  long file_size;
  char* file_content = NULL;

  MPI_Comm_rank(MPI_COMM_WORLD, &rank);

  if (rank == 0) {
    std::ifstream file(filename);

    if (!file.is_open() || !file.good()) {
      MPI_Abort(MPI_COMM_WORLD, MPI_ERR_OTHER);
    }

    assert(file.is_open() && file.good());

    std::filebuf* buf = file.rdbuf();
    file_size = buf->pubseekoff(0, std::ios::end, std::ios::in);
    buf->pubseekpos(0, std::ios::in);
    file_content = static_cast<char*>(pugi::get_memory_allocation_function()(file_size*sizeof(char)));
    buf->sgetn(file_content, file_size);

    file.close();
  }

  //Can't probe on collective communication
  MPI_Bcast(&file_size, 1, MPI_LONG, 0, MPI_COMM_WORLD);
  if (rank != 0)
    file_content = static_cast<char*>(pugi::get_memory_allocation_function()(file_size*sizeof(char)));
  MPI_Bcast(file_content, file_size, MPI_CHAR, 0, MPI_COMM_WORLD);

  doc.load_buffer_inplace_own(file_content, file_size);
#else
  doc.load_file(filename);
#endif

  pugi::xml_node xml_conf = doc.document_element().child("conf");

  max_level = Helper::ToValWithDefault<int>(xml_conf.child_value("max_level"), 7);
  max_iterations = Helper::ToValWithDefault<int>(xml_conf.child_value("max_iterations"), 20);
  pre_smoothing_steps = Helper::ToValWithDefault<int>(xml_conf.child_value("pre_smoothing_steps"), 3);
  post_smoothing_steps = Helper::ToValWithDefault<int>(xml_conf.child_value("post_smoothing_steps"), 3);
  cycle_type = Helper::ToValWithDefault<int>(xml_conf.child_value("cycle_type"), 1);
  near_field_cells = Helper::ToValWithDefault<int>(xml_conf.child_value("near_field_cells"), 3);
  discretization_order = Helper::ToValWithDefault<int>(xml_conf.child_value("discretization_order"), 2);
  precision = Helper::ToValWithDefault<vmg_float>(xml_conf.child_value("precision"), 1.0e-7);
  lop_str = xml_conf.child_value("level_operator");
  datafile_str = xml_conf.child_value("datafile");

  // TODO: Hardcoded periodic boundary conditions are just for testing
  Boundary boundary(Periodic, Periodic, Periodic);

#ifdef HAVE_MPI
  Comm* comm = new CommMPI(boundary, new DomainDecompositionMPI());
#else
  Comm* comm = new CommSerial(boundary);
#endif
  comm->Register("COMM");

  Interface* interface = new InterfaceParticles(boundary, 2, max_level, 0.0, 1.0, near_field_cells, 0, 1.0);
  MG::SetInterface(interface, comm);

  Discretization* discretization = new DiscretizationPoissonFD(discretization_order);
  discretization->Register("DISCRETIZATION");

  LevelOperator* lop;
  if (!lop_str.compare("cs")) {
    lop = new LevelOperatorCS(Stencils::RestrictionFullWeight, Stencils::InterpolationTrilinear);
    Techniques::SetCorrectionSchemePeriodic(interface->MinLevel(), interface->MaxLevel(), cycle_type);
  } else if (!lop_str.compare("fas")) {
    lop = new LevelOperatorFAS(Stencils::RestrictionFullWeight, Stencils::Injection, Stencils::InterpolationTrilinear);
    Techniques::SetFullApproximationSchemePeriodic(interface->MinLevel(), interface->MaxLevel(), cycle_type);
  } else if (!lop_str.compare("cs_debug")) {
    lop = new LevelOperatorCS(Stencils::RestrictionFullWeight, Stencils::InterpolationTrilinear);
    Techniques::SetCorrectionSchemePeriodicDebug(interface->MinLevel(), interface->MaxLevel(), cycle_type);
  } else if (!lop_str.compare("fas_debug")) {
    lop = new LevelOperatorFAS(Stencils::RestrictionFullWeight, Stencils::Injection, Stencils::InterpolationTrilinear);
    Techniques::SetFullApproximationSchemePeriodicDebug(interface->MinLevel(), interface->MaxLevel(), cycle_type);
  } else {
    lop = new LevelOperatorCS(Stencils::RestrictionFullWeight, Stencils::InterpolationTrilinear);
    Techniques::SetCorrectionSchemePeriodic(interface->MinLevel(), interface->MaxLevel(), cycle_type);
  }
  lop->Register("LEVEL_OPERATOR");

  Smoother* smoother = new GaussSeidelRB();
  smoother->Register("SMOOTHER");

#ifdef HAVE_LAPACK
  Solver* solver =  new DGESV<SolverSingular>();
#else
  Solver* solver = new Givens<SolverSingular>();
#endif
  solver->Register("SOLVER");

  factory.RegisterObjectStorage("MAX_ITERATION", max_iterations);
  factory.RegisterObjectStorage("PRESMOOTHSTEPS", pre_smoothing_steps);
  factory.RegisterObjectStorage("POSTSMOOTHSTEPS", post_smoothing_steps);
  factory.RegisterObjectStorage("PRECISION", precision);

  factory.RegisterObjectStorage("PARTICLE_NEAR_FIELD_CELLS", near_field_cells);

  if (rank == 0) {
    std::printf("VMG settings:\n");
    std::printf("  Maximum level:                %d\n", max_level);
    std::printf("  Maximum number of iterations: %d\n", max_iterations);
    std::printf("  Gamma:                        %d\n", cycle_type);
    std::printf("  Pre smoothing steps:          %d\n", pre_smoothing_steps);
    std::printf("  Post smoothing steps:         %d\n", post_smoothing_steps);
    std::printf("  Precision:                    %e\n", precision);
    std::printf("  Near field cells:             %d\n", near_field_cells);
    std::printf("  Discretization order:         %d\n", discretization_order);
  }

  LoadDatafile(filename, datafile_str);
}


void InterfaceParticlesCF::LoadDatafile(const std::string& conf_filename, const std::string& data_filename)
{
  Factory& factory = MG::GetFactory();

#ifdef HAVE_MPI
  int rank;
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  if (rank == 0) {
#endif

    std::ifstream file(data_filename.c_str());

    if (!file.is_open() || !file.good()) {
      std::stringstream search_file;
      search_file << "./" << data_filename;
      file.clear();
      file.open(search_file.str().c_str());
    }

    if (!file.is_open() || !file.good()) {
      std::stringstream search_file;
      search_file << conf_filename.substr(0, conf_filename.find_last_of('/')+1) << data_filename;
      file.clear();
      file.open(search_file.str().c_str());
    }

    assert(file.is_open() && file.good());

    std::stringstream data;
    vmg_float x, y, z, q;
    data << file.rdbuf();

    while (data.good()) {
      data >> x >> y >> z >> q;
      if (data.good()) {
        data_vec.push_back(x);
        data_vec.push_back(y);
        data_vec.push_back(z);
        data_vec.push_back(q);
      }
    }
#ifdef HAVE_MPI
  }
#endif

  factory.RegisterObjectStorage("PARTICLE_NUM_LOCAL", static_cast<int>(data_vec.size()/4));

  vmg_float* x_arr = factory.RegisterObjectStorageArray<vmg_float>("PARTICLE_POS_ARRAY", 3*(data_vec.size()/4));
  vmg_float* q_arr = factory.RegisterObjectStorageArray<vmg_float>("PARTICLE_CHARGE_ARRAY", data_vec.size()/4);
  vmg_float* p_arr = factory.RegisterObjectStorageArray<vmg_float>("PARTICLE_POTENTIAL_ARRAY", data_vec.size()/4);
  vmg_float* f_arr = factory.RegisterObjectStorageArray<vmg_float>("PARTICLE_FIELD_ARRAY", 3*(data_vec.size()/4));

  for (unsigned int i=0; i<data_vec.size()/4; ++i) {
    std::memcpy(&x_arr[3*i], &data_vec[4*i], 3*sizeof(vmg_float));
    q_arr[i] = data_vec[4*i+3];
    p_arr[i] = 0.0;
    f_arr[3*i+0] = 0.0;
    f_arr[3*i+1] = 0.0;
    f_arr[3*i+2] = 0.0;
  }

  data_vec.clear();
}

// void InterfaceParticlesCF::CommParticles(const Grid& grid)
// {
// #ifdef HAVE_MPI
//   int rank, size;
//   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
//   MPI_Comm_size(MPI_COMM_WORLD, &size);

//   std::vector<double> send_buffer[size];

//   int global_extent[6];

//   for (int i=0; i<3; ++i) {
//     global_extent[i] = grid.Global().BeginLocal()[i];
//     global_extent[i+3] = grid.Global().EndLocal()[i];
//   }

//   if (rank == 0) {

//     int global_extent_all[6*size];
//     Index index;

//     MPI_Gather(global_extent, 6, MPI_INT, global_extent_all, 6, MPI_INT, 0, MPI_COMM_WORLD);

//     for (unsigned int i=0; i<data_vec.size()/4; ++i) {
//       index = (static_cast<Vector>(&(data_vec[4*i])) - grid.Extent().Begin()) / grid.Extent().MeshWidth();

//       for (int j=0; j<size; ++j) {
//      if (index[0] >= global_extent_all[6*j+0] && index[0] < global_extent_all[6*j+3] &&
//          index[1] >= global_extent_all[6*j+1] && index[1] < global_extent_all[6*j+4] &&
//          index[2] >= global_extent_all[6*j+2] && index[2] < global_extent_all[6*j+5]) {
//        for (int k=0; k<4; ++k)
//          send_buffer[j].push_back(data_vec[4*i+k]);
//        break;
//      }
//       }

//     }

//     data_vec.clear();

//     for (int i=0; i<size; ++i)
//       MPI_Isend(&send_buffer[i].front(), send_buffer[i].size(), MPI_DOUBLE, i, i, MPI_COMM_WORLD, &Request());

//   }else {
//     MPI_Gather(global_extent, 6, MPI_INT, NULL, 0, MPI_INT, 0, MPI_COMM_WORLD);
//   }

//   MPI_Status status;
//   int count;
//   MPI_Probe(0, rank, MPI_COMM_WORLD, &status);
//   MPI_Get_count(&status, MPI_DOUBLE, &count);

//   data_vec.resize(count);

//   MPI_Recv(&data_vec.front(), count, MPI_DOUBLE, 0, rank, MPI_COMM_WORLD, MPI_STATUS_IGNORE);

// #endif /* HAVE_MPI */
// }