Changeset dfed1c for src/interface


Ignore:
Timestamp:
Nov 22, 2011, 9:22:10 PM (14 years ago)
Author:
Julian Iseringhausen <isering@…>
Children:
facba0
Parents:
66f24d
Message:

Major vmg update.

git-svn-id: ​https://svn.version.fz-juelich.de/scafacos/trunk@1136 5161e1c8-67bf-11de-9fd5-51895aff932f

Location:
src/interface
Files:
2 added
2 deleted
6 edited

Legend:

Unmodified
Added
Removed

  • src/interface/interface.cpp

    r66f24d rdfed1c  
    1212#endif
    1313
     14#include <algorithm>
    1415#include <cmath>
    1516
    … …  
    1920using namespace VMG;
    2021
    21 Interface::Interface(BC boundary, int levelMin_, int levelMax_, vmg_float box_begin_, vmg_float box_end_,
    22                      int coarseningSteps, vmg_float alpha)
     22void Interface::InitInterface(const Vector& box_offset, const vmg_float& box_size,
     23                              const int& coarseningSteps, const vmg_float& alpha)
    2324{
    24   bc = boundary;
    25   levelMin = levelMin_;
    26   levelMax = levelMax_;
     25  Index size_factor_max = 1;
     26  int num_cells;
    2727
    28   bool restricted = false;
     28  const Vector box_center = box_offset + 0.5 * box_size;
    2929
    30   const vmg_float box_center = 0.5 * (box_begin_ + box_end_);
    31   const vmg_float box_size = box_end_ - box_begin_;
    32 
    33   for (int i=0; i<coarseningSteps; ++i) {
     30  for (int i=levelMax; i>=levelMin; --i) {
    3431
    3532    global.push_back(GlobalIndices());
    3633    extent.push_back(SpatialExtent());
    3734
    38     extent.back().SizeFactor() = Helper::intpow(2, static_cast<int>(log(pow(alpha, i+1)) / log(2.0) + 1.0));
     35    num_cells = Helper::intpow(2, i);
    3936
    40     extent.back().Size() = box_size * static_cast<Vector>(extent.back().SizeFactor());
    41     extent.back().Begin() = box_center - 0.5 * extent.back().Size();
    42     extent.back().End() = extent.back().Begin() + extent.back().Size();
    43     extent.back().MeshWidth() = pow(2.0, i-levelMax);
     37    extent.back().SizeFactor() = 1;
     38    extent.back().Size() = box_size;
     39    extent.back().Begin() = box_offset;
     40    extent.back().End() = box_offset + box_size;
     41    extent.back().MeshWidth() = box_size / num_cells;
    4442
    45     global.back().Size() = static_cast<Index>(static_cast<Vector>(extent.back().Size()) / extent.back().MeshWidth() + 0.5) + 1;
    46     global.back().Begin() = 0;
    47     global.back().End() = global.back().Size();
     43    global.back().SizeGlobal() = num_cells + 1;
     44    global.back().SizeLocal() = num_cells + 1;
     45    global.back().BeginFinest() = 0;
     46    global.back().BeginLocal() = 0;
     47    global.back().EndLocal() = num_cells + 1;
     48
     49    if (i == levelMax) {
     50      global.back().SizeFinest() = num_cells + 1;
     51      global.back().EndFinest() = num_cells + 1;
     52      global.back().BoundaryType() = GlobalMax;
     53    }else {
     54      global.back().SizeFinest() = (++global.rbegin())->SizeFinest();
     55      global.back().EndFinest() = (++global.rbegin())->SizeFinest();
     56      global.back().BoundaryType() = GlobalCoarsened;
     57    }
    4858
    4959  }
    5060
    51   for (int i=0; i<coarseningSteps; ++i) {
    52 
    53     if (restricted) {
    54       global[i].BoundaryType() = GlobalCoarsened;
    55     }else if (extent[i].SizeFactor() == extent.back().SizeFactor()) {
    56       global[i].BoundaryType() = GlobalMax;
    57       restricted = true;
    58     }else {
    59       global[i].BoundaryType() = LocallyRefined;
    60     }
    61 
    62     if (i==0) {
    63       global[0].AlignmentBegin() = (global[0].Size() - static_cast<int>( pow(2.0, levelMax) + 0.5) - 1) / 2;
    64     }else if (global[i].BoundaryType() == LocallyRefined)
    65       global[i].AlignmentBegin() = (extent[i].Size() - extent[i-1].Size()) * pow(2.0, levelMax-i-1) + 0.5;
    66     else if (global[i].BoundaryType() == GlobalMax)
    67       global[i].AlignmentBegin() = (extent[i].Size() - extent[i-1].Size()) * pow(2.0, levelMax-i-1) - 0.5;
    68     else
    69       global[i].AlignmentBegin() = 0;
    70 
    71     global[i].AlignmentEnd() = global[i].Size() - global[i].AlignmentBegin();
    72   }
    73 
    74   while (global.size() == 0 || global.back().Size().Max() > Helper::intpow(2,levelMin)+1) {
    75 
    76     global.push_back(GlobalIndices());
    77     extent.push_back(SpatialExtent());
    78 
    79     extent.back().SizeFactor() = (global.size() > 1 ? (++extent.rbegin())->Size() : 1);
    80     extent.back().Size() = (global.size() > 1 ?  (++extent.rbegin())->Size() : box_size);
    81     extent.back().Begin() = box_center - 0.5 * static_cast<Vector>(extent.back().Size());
    82     extent.back().End() = extent.back().Begin() + extent.back().Size();
    83     extent.back().MeshWidth() = pow(2.0, static_cast<int>(global.size())-levelMax-1);
    84 
    85     global.back().Size() = static_cast<Index>(static_cast<Vector>(extent.back().Size()) / extent.back().MeshWidth() + 0.5) + 1;
    86     global.back().Begin() = 0;
    87     global.back().End() = global.back().Size();
    88 
    89     if (global.size() > 1)
    90       global.back().BoundaryType() = GlobalCoarsened;
    91     else
    92       global.back().BoundaryType() = GlobalMax;
    93 
    94     global.back().AlignmentBegin() = 0;
    95     global.back().AlignmentEnd() = global.back().Size();
    96 
    97   }
    98 
    99   if (boundary == Periodic)
    100     for (std::vector<GlobalIndices>::iterator iter=global.begin(); iter!=global.end(); ++iter)
    101       iter->Size() -= 1;
     61  for (int i=0; i<3; ++i)
     62    if (bc[i] == Periodic)
     63      for (unsigned int j=0; j<global.size(); ++j) {
     64        global[j].SizeGlobal()[i] -= 1;
     65        global[j].SizeFinest()[i] -= Helper::intpow(2, j);
     66        global[j].SizeLocal()[i] -= 1;
     67        global[j].EndLocal()[i] -= 1;
     68        global[j].EndFinest()[i] -= Helper::intpow(2, j);
     69      }
    10270
    10371  levelMin = levelMax - global.size() + 1;

  • src/interface/interface.hpp

    r66f24d rdfed1c  
    1414
    1515#include "base/object.hpp"
    16 #include "grid/grid_indexing.hpp"
     16#include "grid/grid_properties.hpp"
    1717
    1818namespace VMG
    … …  
    2525{
    2626public:
    27   Interface(BC boundary, int levelMin, int levelMax, vmg_float box_begin, vmg_float box_end,
    28             int coarseningSteps=0, vmg_float alpha=1.6);
     27  Interface(const Boundary& boundary, const int& levelMin, const int& levelMax,
     28            const Vector& box_offset, const vmg_float box_size,
     29            int coarseningSteps=0, vmg_float alpha=1.6) :
     30    bc(boundary),
     31    levelMin(levelMin),
     32    levelMax(levelMax)
     33  {
     34    InitInterface(box_offset, box_size, coarseningSteps, alpha);
     35  }
    2936
    3037  virtual ~Interface() {}
    … …  
    3542  const std::vector<GlobalIndices>& Global() const {return global;}
    3643  const std::vector<SpatialExtent>& Extent() const {return extent;}
    37   const BC& BoundaryConditions() const {return bc;}
     44
     45  const GlobalIndices& Global(const int& level) const
     46  {
     47    return global[levelMax - level];
     48  }
     49
     50  const SpatialExtent& Extent(const int& level) const
     51  {
     52    return extent[levelMax - level];
     53  }
     54
     55  const Boundary& BoundaryConditions() const {return bc;}
    3856
    3957  const int& MinLevel() const {return levelMin;}
    … …  
    4159
    4260private:
     61  void InitInterface(const Vector& box_offset, const vmg_float& box_end,
     62                     const int& coarseningSteps, const vmg_float& alpha);
     63
    4364  std::vector<GlobalIndices> global;
    4465  std::vector<SpatialExtent> extent;
    4566
    46   BC bc;
     67  Boundary bc;
    4768
    4869  int levelMin, levelMax;

  • src/interface/interface_fcs.cpp

    r66f24d rdfed1c  
    1919
    2020#include "base/object.hpp"
    21 #include "comm/comm_serial.hpp"
     21#include "base/timer.hpp"
     22#include "comm/comm_mpi.hpp"
     23#include "comm/domain_decomposition_mpi.hpp"
    2224#include "interface/interface_fcs.h"
    2325#include "interface/interface_particles.hpp"
    … …  
    3335#include "solver/givens.hpp"
    3436#endif
    35 #include "solver/solver_dirichlet.hpp"
    36 #include "solver/solver_periodic.hpp"
     37#include "solver/solver_regular.hpp"
     38#include "solver/solver_singular.hpp"
    3739
    3840using namespace VMG;
    3941
    40 void VMGInterfaceFCSInit(vmg_int level, vmg_int periodic, vmg_int spline_degree,
    41                          vmg_int max_iter, vmg_int smoothing_steps, vmg_int gamma,
    42                          vmg_float precision, vmg_float box_begin, vmg_float box_end,
    43                          MPI_Comm mpi_comm)
    44 {
    45   const BC bc = (periodic ? Periodic : Dirichlet);
    46 
    47   /*
    48    * Register all the parameters for later use.
    49    */
    50   new ObjectStorage<int>("LEVEL_FCS", level);
    51   new ObjectStorage<int>("PERIODIC_FCS", periodic);
    52   new ObjectStorage<int>("SPLINE_DEGREE_FCS", spline_degree);
    53   new ObjectStorage<int>("MAX_ITER_FCS", max_iter);
    54   new ObjectStorage<int>("SMOOTHING_STEPS_FCS", smoothing_steps);
    55   new ObjectStorage<int>("GAMMA_FCS", gamma);
    56   new ObjectStorage<vmg_float>("PRECISION_FCS", precision);
    57   new ObjectStorage<vmg_float>("BOX_BEGIN_FCS", box_begin);
    58   new ObjectStorage<vmg_float>("BOX_END_FCS", box_end);
    59   new ObjectStorage<MPI_Comm>("MPI_COMM_FCS", mpi_comm);
     42namespace VMGBackupSettings
     43{
     44  static vmg_int level = -1;
     45  static vmg_int periodic[3] = {-1, -1, -1};
     46  static vmg_int max_iter = -1;
     47  static vmg_int smoothing_steps = -1;
     48  static vmg_int gamma = -1;
     49  static vmg_float precision = -1;
     50  static vmg_float box_offset[3];
     51  static vmg_float box_size = -1.0;
     52  static vmg_int near_field_cells = -1;
     53  static MPI_Comm mpi_comm;
     54}
     55
     56static void VMG_fcs_init(vmg_int level, vmg_int* periodic,vmg_int max_iter,
     57                         vmg_int smoothing_steps, vmg_int gamma, vmg_float precision,
     58                         vmg_float* box_offset, vmg_float box_size,
     59                         vmg_int near_field_cells, MPI_Comm mpi_comm)
     60{
     61  VMGBackupSettings::level = level;
     62  std::memcpy(VMGBackupSettings::periodic, periodic, 3*sizeof(vmg_int));
     63  VMGBackupSettings::max_iter = max_iter;
     64  VMGBackupSettings::smoothing_steps = smoothing_steps;
     65  VMGBackupSettings::gamma = gamma;
     66  VMGBackupSettings::precision = precision;
     67  std::memcpy(VMGBackupSettings::box_offset, box_offset, 3*sizeof(vmg_float));
     68  VMGBackupSettings::box_size = box_size;
     69  VMGBackupSettings::near_field_cells = near_field_cells;
     70  VMGBackupSettings::mpi_comm = mpi_comm;
     71
     72  const Boundary boundary(periodic[0] ? Periodic : Open,
     73                          periodic[1] ? Periodic : Open,
     74                          periodic[2] ? Periodic : Open);
     75
     76  const bool singular = periodic[0] * periodic[1] * periodic[2];
    6077
    6178  /*
    6279   * Register communication class.
    63    * For now, parallel version is not yet released.
    64    */
    65   Comm* comm = new CommSerial(bc);
     80   */
     81  Comm* comm = new CommMPI(boundary, new DomainDecompositionMPI());
    6682  comm->Register("COMM");
    6783
    … …  
    6985   * Register particle interface.
    7086   */
    71   Interface* interface = new VMG::InterfaceParticles(bc, 2, level, box_begin, box_end);
     87  Interface* interface = new VMG::InterfaceParticles(boundary, 2, level, Vector(box_offset), box_size, near_field_cells);
    7288  MG::SetInterface(interface, comm);
    7389
    … …  
    96112   */
    97113#ifdef HAVE_LAPACK
    98   Solver* solver = (periodic ?
    99                     static_cast<Solver*>(new DSYSV<SolverPeriodic>()) :
    100                     static_cast<Solver*>(new DGESV<SolverDirichlet>()));
     114  Solver* solver = (singular ?
     115                    static_cast<Solver*>(new DSYSV<SolverSingular>()) :
     116                    static_cast<Solver*>(new DGESV<SolverRegular>()));
    101117#else
    102   Solver* solver = (periodic ?
    103                     static_cast<Solver*>(new Givens<SolverPeriodic>()) :
    104                     static_cast<Solver*>(new Givens<SolverDirichlet>()));
     118  Solver* solver = (singular ?
     119                    static_cast<Solver*>(new Givens<SolverSingular>()) :
     120                    static_cast<Solver*>(new Givens<SolverRegular>()));
    105121#endif
    106122  solver->Register("SOLVER");
    … …  
    109125   * Set commands for the actual multigrid cycle
    110126   */
    111   if (periodic)
     127  if (singular)
    112128    Techniques::SetCorrectionSchemePeriodic(2, level, gamma);
    113129  else
    … …  
    122138  new ObjectStorage<int>("MAX_ITERATION", max_iter);
    123139
     140  new ObjectStorage<int>("PARTICLE_NEAR_FIELD_CELLS", near_field_cells);
     141
    124142  /*
    125143   * Check whether the library is correctly initialized now.
    … …  
    128146}
    129147
    130 void VMGInterfaceFCSRun(vmg_float* x, vmg_float* q, vmg_float* p, vmg_float* f, vmg_int num_particles_local)
     148void VMG_fcs_setup(vmg_int level, vmg_int* periodic, vmg_int max_iter,
     149                   vmg_int smoothing_steps, vmg_int gamma, vmg_float precision,
     150                   vmg_float* box_offset, vmg_float box_size,
     151                   vmg_int near_field_cells, MPI_Comm mpi_comm)
     152{
     153  if (VMGBackupSettings::level != level ||
     154      VMGBackupSettings::periodic[0] != periodic[0] ||
     155      VMGBackupSettings::periodic[1] != periodic[1] ||
     156      VMGBackupSettings::periodic[2] != periodic[2] ||
     157      VMGBackupSettings::max_iter != max_iter ||
     158      VMGBackupSettings::smoothing_steps != smoothing_steps ||
     159      VMGBackupSettings::gamma != gamma ||
     160      VMGBackupSettings::precision != precision ||
     161      VMGBackupSettings::box_offset[0] != box_offset[0] ||
     162      VMGBackupSettings::box_offset[1] != box_offset[1] ||
     163      VMGBackupSettings::box_offset[2] != box_offset[2] ||
     164      VMGBackupSettings::box_size != box_size ||
     165      VMGBackupSettings::near_field_cells != near_field_cells ||
     166      VMGBackupSettings::mpi_comm != mpi_comm) {
     167
     168    VMG_fcs_destroy();
     169    VMG_fcs_init(level, periodic, max_iter,
     170                 smoothing_steps, gamma, precision,
     171                 box_offset, box_size, near_field_cells,
     172                 mpi_comm);
     173
     174  }
     175}
     176
     177void VMG_fcs_run(vmg_float* x, vmg_float* q, vmg_float* p, vmg_float* f, vmg_int num_particles_local)
    131178{
    132179  /*
    133180   * Register parameters for later use.
    134181   */
    135   new ObjectStorage<vmg_float*>("X_FCS", x);
    136   new ObjectStorage<vmg_float*>("Q_FCS", q);
    137   new ObjectStorage<vmg_float*>("P_FCS", p);
    138   new ObjectStorage<vmg_float*>("F_FCS", f);
    139   new ObjectStorage<vmg_int>("NUM_PARTICLES_LOCAL_FCS", num_particles_local);
     182  new ObjectStorage<vmg_float*>("PARTICLE_POS_ARRAY", x);
     183  new ObjectStorage<vmg_float*>("PARTICLE_CHARGE_ARRAY", q);
     184  new ObjectStorage<vmg_float*>("PARTICLE_POTENTIAL_ARRAY", p);
     185  new ObjectStorage<vmg_float*>("PARTICLE_FORCE_ARRAY", f);
     186  new ObjectStorage<vmg_int>("PARTICLE_NUM_LOCAL", num_particles_local);
    140187
    141188  /*
    … …  
    145192}
    146193
    147 void VMGInterfaceFCSDestroy(void)
     194void VMG_fcs_print_timer()
     195{
     196  Timer::Print();
     197}
     198
     199void VMG_fcs_destroy(void)
    148200{
    149201  /*

  • src/interface/interface_fcs.h

    r66f24d rdfed1c  
    1313#endif
    1414
    15 void VMGInterfaceFCSInit(vmg_int level, vmg_int periodic, vmg_int spline_degree,
    16                          vmg_int max_iter, vmg_int smoothing_steps, vmg_int gamma,
    17                          vmg_float precision, vmg_float box_begin, vmg_float box_end,
    18                         MPI_Comm mpi_comm);
     15void VMG_fcs_setup(vmg_int level, vmg_int* periodic, vmg_int max_iter,
     16                   vmg_int smoothing_steps, vmg_int gamma, vmg_float precision,
     17                   vmg_float* box_offset, vmg_float box_size,
     18                   vmg_int near_field_cells, MPI_Comm mpi_comm);
    1919
    20 void VMGInterfaceFCSRun(vmg_float* x, vmg_float* q, vmg_float* p, vmg_float* f, vmg_int num_particles_local);
     20void VMG_fcs_run(vmg_float* pos, vmg_float* charge, vmg_float* potential, vmg_float* f, vmg_int num_particles_local);
    2121
    22 void VMGInterfaceFCSDestroy(void);
     22void VMG_fcs_print_timer(void);
     23
     24void VMG_fcs_destroy(void);
    2325
    2426#ifdef __cplusplus

  • src/interface/interface_particles.cpp

    r66f24d rdfed1c  
    88 */
    99
    10 
    1110#ifdef HAVE_CONFIG_H
    1211#include <config.h>
    1312#endif
    1413
     14#ifdef HAVE_MPI
     15#include <mpi.h>
     16#endif
     17
    1518#include <cmath>
     19#include <cstring>
    1620
    1721#include "base/helper.hpp"
    1822#include "base/index.hpp"
     23#include "base/linked_cell_list.hpp"
     24#include "base/timer.hpp"
    1925#include "base/vector.hpp"
    2026#include "comm/comm.hpp"
    … …  
    2935void InterfaceParticles::ImportRightHandSide(Multigrid& multigrid)
    3036{
     37  Timer::Start("Particle/Grid Interpolation");
     38
     39  Index index_global, index_local, index;
     40  Vector pos_rel, pos_abs, grid_val;
     41  vmg_float length;
     42
    3143  Factory& factory = MG::GetFactory();
    3244
    33   const int& deg = factory.GetObjectStorageVal<int>("SPLINE_DEGREE_FCS");
    34   const vmg_int& num_particles_local = factory.GetObjectStorageVal<vmg_int>("NUM_PARTICLES_LOCAL_FCS");
    35   const int deg_2 = deg / 2;
    36 
    37   vmg_float* x = factory.GetObjectStorageVal<vmg_float*>("X_FCS");
    38   vmg_float* q = factory.GetObjectStorageVal<vmg_float*>("Q_FCS");
    39 
    40   Index index_global, index_local, index;
    41   Vector pos_rel, pos_abs, index_float, grid_val;
    42   vmg_float foo;
    43   vmg_float MMtbl[3][deg];
     45  const int& near_field_cells = factory.GetObjectStorageVal<int>("PARTICLE_NEAR_FIELD_CELLS");
    4446
    4547  TempGrid* temp_grid = GetTempGrid(0);
    4648  Grid& grid = multigrid(multigrid.GlobalMaxLevel());
    4749  LocalIndices local = grid.Local();
     50
     51  //TODO This has to be replaced for a VMG::Vector for different mesh widths on the axis
     52  const vmg_float r_cut = (near_field_cells+0.5) * grid.Extent().MeshWidth().Max();
    4853
    4954  // We don't need a boundary on this grid
    … …  
    5358  local.BoundaryBegin2() = local.BoundaryEnd2() = 0;
    5459
    55   //
    56   // Check if we need to communicate and resize halo sizes
    57   //
    58   if ((local.HaloEnd1() - local.HaloBegin1()).Product() > 0) {
    59     local.HaloBegin1() = 0;
    60     local.HaloEnd1() = local.Begin() = deg_2;
    61     local.End() = local.Begin() + local.Size();
    62   }
    63 
    64   if ((local.HaloEnd2() - local.HaloBegin2()).Product() > 0) {
    65     local.HaloBegin2() = local.End();
    66     local.HaloEnd2() = local.HaloBegin2() + deg_2;
     60  // Set grid size of intermediate temporary grid
     61  for (int i=0; i<3; ++i) {
     62
     63    if (local.HasHalo1()[i]) {
     64      local.HaloBegin1()[i] = 0;
     65      local.HaloEnd1()[i] = local.Begin()[i] = near_field_cells+1;
     66      local.End()[i] = local.Begin()[i] + local.Size()[i];
     67    }
     68
     69    if (local.HasHalo2()[i]) {
     70      local.HaloBegin2()[i] = local.End()[i];
     71      local.HaloEnd2()[i] = local.HaloBegin2()[i] + near_field_cells+1;
     72    }
     73
    6774  }
    6875
    … …  
    7279
    7380  temp_grid->SetProperties(grid.Global(), local, grid.Extent());
    74 
    7581  temp_grid->Clear();
    7682
    77   for (vmg_int i=0; i<num_particles_local; ++i) {
    78 
    79     //TODO: Check for correct conversions here
    80     index_float = (static_cast<Vector>(&x[3*i]) - grid.Extent().Begin()) / grid.Extent().MeshWidth();
    81     index_global = static_cast<Index>(index_float);
    82     index_local = temp_grid->GlobalToLocalIndex(index_global);
    83 
    84     Helper::MM(deg, std::modf(index_float[0], &foo), MMtbl[0]);
    85     Helper::MM(deg, std::modf(index_float[1], &foo), MMtbl[1]);
    86     Helper::MM(deg, std::modf(index_float[2], &foo), MMtbl[2]);
    87 
    88     for (index.X() = 0; index.X() < deg; ++index.X()) {
    89 
    90       grid_val.X() = q[i] * MMtbl[0][index.X()];
    91 
    92       for (index.Y() = 0; index.Y() < deg; ++index.Y()) {
    93 
    94         grid_val.Y() = grid_val.X() * MMtbl[1][index.Y()];
    95 
    96         for (index.Z() = 0; index.Z() < deg; ++index.Z()) {
    97 
    98           (*temp_grid)(index_local + index - deg_2) -= grid_val.Y() * MMtbl[2][index.Z()];
    99 
    100         }
    101       }
    102     }
     83  /*
     84   * Distribute particles to their processes
     85   */
     86  particles.clear();
     87  MG::GetComm()->CommParticles(grid, particles);
     88
     89  /*
     90   * Charge assignment on the grid
     91   */
     92  std::list<Particle::Particle>::iterator iter;
     93
     94#ifdef DEBUG_OUTPUT
     95  vmg_float particle_charges = 0.0;
     96  for (iter=particles.begin(); iter!=particles.end(); ++iter)
     97    particle_charges += iter->Charge();
     98  particle_charges = MG::GetComm()->GlobalSumRoot(particle_charges);
     99  MG::GetComm()->PrintStringOnce("Particle list charge sum: %e", particle_charges);
     100#endif
     101
     102  for (iter=particles.begin(); iter!=particles.end(); ++iter) {
     103
     104    // Compute global and local grid index of the lower left corner of the grid cell containing the particle
     105    index_global = static_cast<Index>((iter->Pos() - temp_grid->Extent().Begin()) / temp_grid->Extent().MeshWidth());
     106    index_local = index_global - temp_grid->Global().BeginLocal() + temp_grid->Local().Begin();
     107
     108    // Iterate over all grid points which lie in the support of the interpolating B-Spline
     109    for (index.X()=-1*near_field_cells-1; index.X()<=near_field_cells+1; ++index.X())
     110      for (index.Y()=-1*near_field_cells-1; index.Y()<=near_field_cells+1; ++index.Y())
     111        for (index.Z()=-1*near_field_cells-1; index.Z()<=near_field_cells+1; ++index.Z()) {
     112
     113          // Compute distance from grid point to particle
     114          length = ( iter->Pos() - grid.Extent().Begin() - grid.Extent().MeshWidth() * (index_global+index) ).Length();
     115
     116          // If grid point lies in the support of the B-Spline, do charge interpolation
     117          if (length < r_cut)
     118            (*temp_grid)(index_local+index) += spl.EvaluateSpline(length) * iter->Charge();
     119
     120        }
     121
    103122  }
    104123
     124  // Communicate charges over halo
    105125  MG::GetComm()->CommFromGhosts(*temp_grid);
    106126
     127  // Assign charge values to the right hand side
    107128  for (index.X()=0; index.X()<grid.Local().Size().X(); ++index.X())
    108129    for (index.Y()=0; index.Y()<grid.Local().Size().Y(); ++index.Y())
    109130      for (index.Z()=0; index.Z()<grid.Local().Size().Z(); ++index.Z())
    110         grid(index + grid.Local().Begin()) = (*temp_grid)(index + temp_grid->Local().Begin());
     131        grid(index + grid.Local().Begin()) = temp_grid->GetVal(index + temp_grid->Local().Begin());
     132
     133  Timer::Stop("Particle/Grid Interpolation");
     134
     135#ifdef DEBUG_OUTPUT
     136  Grid::iterator grid_iter;
     137  vmg_float charge_sum = 0.0;
     138  for (grid_iter=grid.Iterators().Local().Begin(); grid_iter!=grid.Iterators().Local().End(); ++grid_iter)
     139    charge_sum += grid.GetVal(*grid_iter);
     140  charge_sum = MG::GetComm()->GlobalSum(charge_sum);
     141  MG::GetComm()->PrintStringOnce("Grid charge sum: %e", charge_sum);
     142#endif
    111143}
    112144
    113145void InterfaceParticles::ExportSolution(Grid& grid)
    114146{
     147  Timer::Start("Grid/Particle Interpolation");
     148
     149  Index index;
     150  vmg_float length;
     151  Vector dist_vec;
     152
     153#ifdef DEBUG_OUTPUT
     154  vmg_float e = 0.0;
     155  vmg_float e_long = 0.0;
     156  vmg_float e_self = 0.0;
     157#endif
     158
    115159  Factory& factory = MG::GetFactory();
    116 
    117   const int& deg = factory.GetObjectStorageVal<int>("SPLINE_DEGREE_FCS");
    118   const vmg_int& num_particles_local = factory.GetObjectStorageVal<vmg_int>("NUM_PARTICLES_LOCAL_FCS");
    119   const Vector h_inv = 1.0 / grid.Extent().MeshWidth();
    120 
    121   vmg_float* x = factory.GetObjectStorageVal<vmg_float*>("X_FCS");
    122   vmg_float* q = factory.GetObjectStorageVal<vmg_float*>("Q_FCS");
    123   vmg_float* f = factory.GetObjectStorageVal<vmg_float*>("F_FCS");
    124 
    125   Index index, index_global, index_local;
    126   Vector index_float, spl, der;
    127   vmg_float foo;
    128   vmg_float MMtbl[3][deg], MdMtbl[3][deg];
    129 
     160  Comm* comm = MG::GetComm();
     161
     162  /*
     163   * Get parameters and arrays
     164   */
     165  const vmg_int& num_particles_local = factory.GetObjectStorageVal<vmg_int>("PARTICLE_NUM_LOCAL");
     166  const int& near_field_cells = factory.GetObjectStorageVal<int>("PARTICLE_NEAR_FIELD_CELLS");
     167  vmg_float* p = factory.GetObjectStorageArray<vmg_float>("PARTICLE_POTENTIAL_ARRAY");
     168
     169  const vmg_float r_cut = (near_field_cells+0.5) * grid.Extent().MeshWidth().Max();
     170
     171  /*
     172   * Initialze potential
     173   */
     174  for (vmg_int i=0; i<num_particles_local; ++i)
     175    p[i] = 0.0;
     176
     177  /*
     178   * Copy potential values to a grid with sufficiently large halo size.
     179   * This may be optimized in future.
     180   * The parameters of this grid have been set in the import step.
     181   */
    130182  TempGrid* temp_grid = GetTempGrid(0);
    131183
    … …  
    133185    for (index.Y()=0; index.Y()<grid.Local().Size().Y(); ++index.Y())
    134186      for (index.Z()=0; index.Z()<grid.Local().Size().Z(); ++index.Z())
    135         (*temp_grid)(index + temp_grid->Local().Begin()) = grid(index + grid.Local().Begin());
    136 
    137   MG::GetComm()->CommToGhosts(*temp_grid);
    138 
    139   for (vmg_int i=0; i<num_particles_local; ++i) {
    140 
    141     index_float = (static_cast<Vector>(&x[3*i]) - grid.Extent().Begin()) / grid.Extent().MeshWidth();
    142     index_global = static_cast<Index>(index_float);
    143     index_local = temp_grid->GlobalToLocalIndex(index_global);
    144 
    145     Helper::MdM(deg, std::modf(index_float[0], &foo), MMtbl[0], MdMtbl[0]);
    146     Helper::MdM(deg, std::modf(index_float[1], &foo), MMtbl[1], MdMtbl[1]);
    147     Helper::MdM(deg, std::modf(index_float[2], &foo), MMtbl[2], MdMtbl[2]);
    148 
    149     for (int j=0; j<3; ++j)
    150       f[3*i+j] = 0.0;
    151 
    152     for (index.X() = 0; index.X() < deg; ++index.X()) {
    153 
    154       spl.X() = MMtbl[0][index.X()];
    155       der.X() = MdMtbl[0][index.X()];
    156 
    157       for (index.Y() = 0; index.Y() < deg; ++index.Y()) {
    158 
    159         spl.Y() = MMtbl[1][index.Y()];
    160         der.Y() = MdMtbl[1][index.Y()];
    161 
    162         for (index.Z() = 0; index.Z() < deg; ++index.Z()) {
    163 
    164           spl.Z() = MMtbl[2][index.Z()];
    165           der.Z() = MdMtbl[2][index.Z()];
    166 
    167           f[3*i  ] += temp_grid->GetVal(index + temp_grid->Local().Begin()) *
    168                       der[0]*spl[1]*spl[2] * h_inv[0];
    169 
    170           f[3*i+1] += temp_grid->GetVal(index + temp_grid->Local().Begin()) *
    171                       spl[0]*der[1]*spl[2] * h_inv[1];
    172 
    173           f[3*i+2] += temp_grid->GetVal(index + temp_grid->Local().Begin()) *
    174                       spl[0]*spl[1]*der[2] * h_inv[2];
    175 
    176         }
    177       }
    178     }
    179 
    180     for (int j=0; j<3; ++j)
    181       f[3*i+j] *= q[i];
     187        (*temp_grid)(index + temp_grid->Local().Begin()) = grid.GetVal(index + grid.Local().Begin());
     188
     189  comm->CommToGhosts(*temp_grid);
     190
     191  /*
     192   * Do potential back-interpolation
     193   */
     194  std::list<Particle::Particle>::iterator p_iter;
     195  for (p_iter=particles.begin(); p_iter!=particles.end(); ++p_iter) {
     196
     197    // Interpolate long range part of potential
     198    p_iter->Pot() = 0.5 * (Helper::InterpolateTrilinear(p_iter->Pos(), *temp_grid) - p_iter->Charge() * spl.GetAntiDerivativeZero());
     199
     200#ifdef DEBUG_OUTPUT
     201    e_long += 0.5 * p_iter->Charge() * Helper::InterpolateTrilinear(p_iter->Pos(), *temp_grid);
     202    e_self += 0.5 * p_iter->Charge() * p_iter->Charge() * spl.GetAntiDerivativeZero();
     203#endif
    182204
    183205  }
    184206
     207  comm->CommAddPotential(particles);
     208
     209  /*
     210   * Compute near field correction
     211   */
     212  Particle::LinkedCellList lc(particles, near_field_cells, grid);
     213  Particle::LinkedCellList::iterator p1, p2;
     214  Grid::iterator iter;
     215
     216  comm->CommLCListGhosts(grid, lc);
     217
     218  for (iter=lc.Iterators().Local().Begin(); iter!=lc.Iterators().Local().End(); ++iter)
     219    for (p1=lc(*iter).begin(); p1!=lc(*iter).end(); ++p1)
     220      for (index.X()=-1*near_field_cells-1; index.X()<=near_field_cells+1; ++index.X())
     221        for (index.Y()=-1*near_field_cells-1; index.Y()<=near_field_cells+1; ++index.Y())
     222          for (index.Z()=-1*near_field_cells-1; index.Z()<=near_field_cells+1; ++index.Z())
     223            for (p2=lc(*iter+index).begin(); p2!=lc(*iter+index).end(); ++p2)
     224              if (*p1 != *p2) {
     225
     226                length = (p2->Pos() - p1->Pos()).Length();
     227
     228                //TODO Rewrite this equation more efficiently
     229                if (length < r_cut)
     230                  p1->Pot() += 0.5 * p2->Charge() / length * (0.25*M_1_PI + spl.EvaluatePotential(length));
     231
     232              }
     233
     234  comm->CommAddPotential(lc);
     235
     236  Timer::Stop("Grid/Particle Interpolation");
     237
     238#ifdef DEBUG_OUTPUT
     239  vmg_float* q = factory.GetObjectStorageArray<vmg_float>("PARTICLE_CHARGE_ARRAY");
     240
     241  for (int i=0; i<num_particles_local; ++i)
     242    e += p[i] * q[i];
     243
     244  e = comm->GlobalSumRoot(e);
     245  e_long = comm->GlobalSumRoot(e_long);
     246  e_self = comm->GlobalSumRoot(e_self);
     247
     248  comm->PrintStringOnce("E_long:  %e", e_long);
     249  comm->PrintStringOnce("E_short: %e", e - e_long + e_self);
     250  comm->PrintStringOnce("E_self:  %e", e_self);
     251  comm->PrintStringOnce("E_total: %e", e);
     252
     253#endif /* DEBUG_OUTPUT */
     254
    185255}
    186 

  • src/interface/interface_particles.hpp

    r66f24d rdfed1c  
    1212#define INTERFACE_PARTIVLES_HPP
    1313
     14#include <list>
     15
    1416#include "base/has_tempgrids.hpp"
     17#include "base/linked_cell_list.hpp"
     18#include "base/particle.hpp"
     19#include "comm/mpi/has_request_vec.hpp"
    1520#include "interface/interface.hpp"
     21#include "samples/bspline.hpp"
    1622
    1723namespace VMG
    … …  
    2127class Multigrid;
    2228
    23 class InterfaceParticles : public Interface, public HasTempGrids
     29class InterfaceParticles : public Interface, public HasTempGrids, public HasRequestVec
    2430{
    2531public:
    26   InterfaceParticles(BC boundary, int levelMin, int levelMax, vmg_float box_begin, vmg_float box_end) :
    27     Interface(boundary, levelMin, levelMax, box_begin, box_end, 0, 1.0),
    28     HasTempGrids()
     32  InterfaceParticles(const Boundary& boundary, const int& levelMin, const int& levelMax,
     33                     const Vector& box_offset, const vmg_float& box_size,
     34                     const int& near_field_cells) :
     35    Interface(boundary, levelMin, levelMax, box_offset, box_size, 0, 1.0),
     36    HasTempGrids(),
     37    HasRequestVec(),
     38    spl((near_field_cells+0.5)*Extent(MaxLevel()).MeshWidth().Max())
    2939  {}
    3040
    3141  void ImportRightHandSide(Multigrid& multigrid);
    3242  void ExportSolution(Grid& grid);
     43
     44protected:
     45  BSpline spl;
     46
     47private:
     48  std::list<Particle::Particle> particles;
    3349};
    3450
Note: See TracChangeset for help on using the changeset viewer.