Changeset ac6d04 for src/interface


Ignore:
Timestamp:
Apr 10, 2012, 1:55:49 PM (14 years ago)
Author:
Julian Iseringhausen <isering@…>
Children:
a40eea
Parents:
d24c2f
Message:

Merge recent changes of the vmg library into ScaFaCos.

Includes a fix for the communication problems on Jugene.

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

Location:
src/interface
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • src/interface/interface.cpp

    rd24c2f rac6d04  
    2323                              const int& coarseningSteps, const vmg_float& alpha)
    2424{
    25   Index size_factor_max = 1;
    26   int num_cells;
     25  int i;
     26  Index num_cells = Helper::intpow(2, levelMax);
     27  Index size_factor;
    2728
    2829  const Vector box_center = box_offset + 0.5 * box_size;
    2930
    30   for (int i=levelMax; i>=levelMin; --i) {
     31  for (i=0; i<coarseningSteps; ++i) {
    3132
    3233    global.push_back(GlobalIndices());
    3334    extent.push_back(SpatialExtent());
    3435
    35     num_cells = Helper::intpow(2, i);
     36    for (int j=0; j<3; ++j)
     37      size_factor[j] = (bc[j] == Periodic ? 1 : Helper::intpow(2, static_cast<int>(log(pow(alpha, i+1)) / log(2.0) + 1.0)));
    3638
    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;
     39    extent.back().Size() = box_size * static_cast<Vector>(size_factor);
     40    extent.back().Begin() = box_center - 0.5 * extent.back().Size();
     41    extent.back().End() = extent.back().Begin() + extent.back().Size();
     42    extent.back().MeshWidth() = pow(2.0, i-levelMax);
    4243
    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;
     44    num_cells = Helper::intpow(2,levelMax-i) * size_factor;
    4845
    49     if (i == levelMax) {
    50       global.back().SizeFinest() = num_cells + 1;
    51       global.back().EndFinest() = num_cells + 1;
    52       global.back().BoundaryType() = GlobalMax;
     46    global.back().GlobalSize() = num_cells + 1;
     47    global.back().LocalSize() = num_cells + 1;
     48    global.back().LocalBegin() = 0;
     49    global.back().LocalEnd() = num_cells + 1;
     50
     51    global.back().FinestAbsSize() = Helper::intpow(2,i) * num_cells + 1;
     52    global.back().FinestAbsBegin() = ((global.back().FinestAbsSize()-1) * (1-size_factor)) / (2*size_factor);
     53    global.back().FinestAbsEnd() = global.back().FinestAbsBegin() + global.back().FinestAbsSize();
     54
     55    if (i==0)
     56      global.back().GlobalFinerBegin() = (num_cells - Helper::intpow(2, levelMax-i))/2;
     57    else
     58      global.back().GlobalFinerBegin() = (global.back().FinestAbsSize() - (++global.rbegin())->FinestAbsSize()) / Helper::intpow(2,i+1);
     59
     60    global.back().GlobalFinerEnd() = global.back().GlobalSize() - global.back().GlobalFinerBegin();
     61    global.back().GlobalFinerSize() = global.back().GlobalFinerEnd() - global.back().GlobalFinerBegin();
     62
     63    global.back().LocalFinerBegin() = global.back().GlobalFinerBegin();
     64    global.back().LocalFinerEnd() = global.back().GlobalFinerEnd();
     65    global.back().LocalFinerSize() = global.back().GlobalFinerSize();
     66  }
     67
     68  while (global.size() == 0 || global.back().GlobalSize().Min() > Helper::intpow(2, levelMin)+1) {
     69
     70    if (global.size() > 0)
     71      num_cells /= 2;
     72
     73    global.push_back(GlobalIndices());
     74    extent.push_back(SpatialExtent());
     75
     76    if (global.size() == 1) {
     77      extent.back().Size() = box_size;
     78      extent.back().Begin() = box_offset;
     79      extent.back().End() = box_offset + box_size;
     80      extent.back().MeshWidth() = box_size / static_cast<Vector>(num_cells);
    5381    }else {
    54       global.back().SizeFinest() = (++global.rbegin())->SizeFinest();
    55       global.back().EndFinest() = (++global.rbegin())->SizeFinest();
    56       global.back().BoundaryType() = GlobalCoarsened;
     82      extent.back().Size() = (++extent.rbegin())->Size();
     83      extent.back().Begin() = (++extent.rbegin())->Begin();
     84      extent.back().End() = (++extent.rbegin())->End();
     85      extent.back().MeshWidth() = 2.0 * (++extent.rbegin())->MeshWidth();
    5786    }
     87
     88    global.back().GlobalSize() = num_cells + 1;
     89    global.back().LocalSize() = num_cells + 1;
     90    global.back().LocalBegin() = 0;
     91    global.back().LocalEnd() = num_cells + 1;
     92
     93    if (global.size() == 1) {
     94      global.back().FinestAbsBegin() = 0;
     95      global.back().FinestAbsEnd() = global.back().GlobalSize();
     96      global.back().FinestAbsSize() = global.back().GlobalSize();
     97    }else {
     98      global.back().FinestAbsBegin() = (++global.rbegin())->FinestAbsBegin();
     99      global.back().FinestAbsEnd() = (++global.rbegin())->FinestAbsEnd();
     100      global.back().FinestAbsSize() = (++global.rbegin())->FinestAbsSize();
     101    }
     102
     103    global.back().GlobalFinerBegin() = 0;
     104    global.back().GlobalFinerEnd() = global.back().GlobalSize();
     105    global.back().GlobalFinerSize() = global.back().GlobalSize();
     106
     107    global.back().LocalFinerBegin() = global.back().GlobalFinerBegin();
     108    global.back().LocalFinerEnd() = global.back().GlobalFinerEnd();
     109    global.back().LocalFinerSize() = global.back().GlobalFinerSize();
    58110
    59111  }
    60112
    61   for (int i=0; i<3; ++i)
     113  for (i=0; i<3; ++i)
    62114    if (bc[i] == Periodic)
    63115      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);
     116        global[j].GlobalSize()[i] -= 1;
     117        global[j].FinestAbsSize()[i] -= Helper::intpow(2, j);
     118        global[j].FinestAbsEnd()[i] -= Helper::intpow(2, j);
     119        global[j].LocalSize()[i] -= 1;
     120        global[j].LocalEnd()[i] -= 1;
     121        global[j].GlobalFinerSize()[i] -= 1;
     122        global[j].GlobalFinerEnd()[i] -= 1;
     123        global[j].LocalFinerSize()[i] -= 1;
     124        global[j].LocalFinerEnd()[i] -= 1;
    69125      }
    70126
    71127  levelMin = levelMax - global.size() + 1;
     128
     129  global.back().BoundaryType() = GlobalCoarsened;
     130
     131  for (i=global.size()-2; i>=0; --i) {
     132    if (global[i].FinestAbsSize().Product() >= global[i+1].FinestAbsSize().Product()) {
     133      global[i].BoundaryType() = GlobalCoarsened;
     134    }else {
     135      global[i].BoundaryType() = LocallyRefined;
     136      global[i+1].BoundaryType() = GlobalMax;
     137      break;
     138    }
     139  }
     140
     141  for (; i>=0; --i)
     142    global[i].BoundaryType() = LocallyRefined;
     143
     144  if (global.front().BoundaryType() != LocallyRefined &&
     145      global.front().BoundaryType() != GlobalMax)
     146    global.front().BoundaryType() = GlobalMax;
     147
    72148}

  • src/interface/interface_fcs.cpp

    rd24c2f rac6d04  
    1717
    1818#include <mpi.h>
     19#ifdef HAVE_MARMOT
     20#include <enhancempicalls.h>
     21#include <sourceinfompicalls.h>
     22#endif
    1923
    2024#include "base/object.hpp"
    … …  
    2226#include "comm/comm_mpi.hpp"
    2327#include "comm/domain_decomposition_mpi.hpp"
     28#include "comm/mpi/error_handler.hpp"
    2429#include "interface/interface_fcs.h"
    2530#include "interface/interface_particles.hpp"
    2631#include "level/level_operator_cs.hpp"
     32#include "level/level_operator_fas.hpp"
    2733#include "samples/discretization_poisson_fd.hpp"
     34#include "samples/discretization_poisson_fv.hpp"
    2835#include "samples/stencils.hpp"
    2936#include "samples/techniques.hpp"
    … …  
    7077  VMGBackupSettings::mpi_comm = mpi_comm;
    7178
     79  MPI_Errhandler mpiErrorHandler;
     80  MPI_Comm_create_errhandler(VMG::MPI::ConvertToException, &mpiErrorHandler);
     81  MPI_Comm_set_errhandler(MPI_COMM_WORLD, mpiErrorHandler);
     82
    7283  const Boundary boundary(periodic[0] ? Periodic : Open,
    7384                          periodic[1] ? Periodic : Open,
    … …  
    8596   * Register particle interface.
    8697   */
    87   Interface* interface = new VMG::InterfaceParticles(boundary, 2, level, Vector(box_offset), box_size, near_field_cells);
     98  Interface* interface;
     99  if (singular)
     100    interface = new InterfaceParticles(boundary, 2, level, Vector(box_offset), box_size, near_field_cells, 0, 1.0);
     101  else
     102    interface = new InterfaceParticles(boundary, 2, level, Vector(box_offset), box_size, near_field_cells, 2, 1.6);
    88103  MG::SetInterface(interface, comm);
    89104
    90105  /*
    91    * Use a finite difference discretization of the Poisson equation
    92    */
    93   Discretization* discretization = new DiscretizationPoissonFD();
     106   * Define discretization of the Poisson equation.
     107   * Finite volumes for locally refined grids and
     108   * finite differences otherwise.
     109   */
     110  Discretization* discretization;
     111  if (singular)
     112    discretization = new DiscretizationPoissonFD();
     113  else
     114    discretization = new DiscretizationPoissonFV();
    94115  discretization->Register("DISCRETIZATION");
    95116
    … …  
    99120   * with a half-weight version once debugging is finished.
    100121   */
    101   LevelOperator* lop = new LevelOperatorCS(Stencils::RestrictionFullWeight, Stencils::InterpolationTrilinear);
     122  LevelOperator* lop;
     123  if (singular)
     124    lop = new LevelOperatorCS(Stencils::RestrictionFullWeight, Stencils::InterpolationTrilinear);
     125  else
     126    lop = new LevelOperatorFAS(Stencils::RestrictionFullWeight, Stencils::Injection, Stencils::InterpolationTrilinear);
    102127  lop->Register("LEVEL_OPERATOR");
    103128
    … …  
    126151   */
    127152  if (singular)
    128     Techniques::SetCorrectionSchemePeriodic(2, level, gamma);
    129   else
    130     Techniques::SetCorrectionSchemeDirichlet(2, level, gamma);
     153    Techniques::SetCorrectionSchemePeriodic(interface->MinLevel(), interface->MaxLevel(), gamma);
     154  else
     155    Techniques::SetFullApproximationSchemeDirichlet(interface->MinLevel(), interface->MaxLevel(), gamma);
    131156
    132157  /*
    … …  
    137162  new ObjectStorage<vmg_float>("PRECISION", precision);
    138163  new ObjectStorage<int>("MAX_ITERATION", max_iter);
    139 
    140164  new ObjectStorage<int>("PARTICLE_NEAR_FIELD_CELLS", near_field_cells);
    141165
    … …  
    180204}
    181205
     206int VMG_fcs_check()
     207{
     208  const int& near_field_cells = MG::GetFactory().GetObjectStorageVal<int>("PARTICLE_NEAR_FIELD_CELLS");
     209  const Multigrid& multigrid = *MG::GetRhs();
     210  const Grid& grid = multigrid(multigrid.MaxLevel());
     211
     212  if (!grid.Global().LocalSize().IsComponentwiseGreater(near_field_cells))
     213    return 1;
     214
     215  return 0;
     216}
     217
    182218void VMG_fcs_run(vmg_float* x, vmg_float* q, vmg_float* p, vmg_float* f, vmg_int num_particles_local)
    183219{

  • src/interface/interface_fcs.h

    rd24c2f rac6d04  
    1313#endif
    1414
    15 void 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);
     15void VMG_fcs_setup(fcs_int level, fcs_int* periodic, fcs_int max_iter,
     16                   fcs_int smoothing_steps, fcs_int gamma, fcs_float precision,
     17                   fcs_float* box_offset, fcs_float box_size,
     18                   fcs_int near_field_cells, MPI_Comm mpi_comm);
    1919
    20 void VMG_fcs_run(vmg_float* pos, vmg_float* charge, vmg_float* potential, vmg_float* f, vmg_int num_particles_local);
     20int VMG_fcs_check();
     21
     22void VMG_fcs_run(fcs_float* pos, fcs_float* charge, fcs_float* potential, fcs_float* f, fcs_int num_particles_local);
    2123
    2224void VMG_fcs_print_timer(void);

  • src/interface/interface_particles.cpp

    rd24c2f rac6d04  
    1414#ifdef HAVE_MPI
    1515#include <mpi.h>
     16#ifdef HAVE_MARMOT
     17#include <enhancempicalls.h>
     18#include <sourceinfompicalls.h>
     19#endif
    1620#endif
    1721
    … …  
    2226#include "base/index.hpp"
    2327#include "base/linked_cell_list.hpp"
     28#include "base/math.hpp"
    2429#include "base/vector.hpp"
    2530#include "comm/comm.hpp"
    … …  
    4247  const int& near_field_cells = factory.GetObjectStorageVal<int>("PARTICLE_NEAR_FIELD_CELLS");
    4348
    44   Grid& grid = multigrid(multigrid.GlobalMaxLevel());
     49  Grid& grid = multigrid(multigrid.MaxLevel());
    4550  Grid& particle_grid = comm.GetParticleGrid();
    4651
    4752  particle_grid.Clear();
     53
     54  assert(particle_grid.Global().LocalSize().IsComponentwiseGreater(near_field_cells));
    4855
    4956  /*
    … …  
    6471  particle_charges = MG::GetComm()->GlobalSumRoot(particle_charges);
    6572  comm.PrintStringOnce("Particle list charge sum: %e", particle_charges);
     73  comm.PrintString("Local number of particles: %d", particles.size());
    6674#endif
    6775
    … …  
    7684    for (index.Y()=0; index.Y()<grid.Local().Size().Y(); ++index.Y())
    7785      for (index.Z()=0; index.Z()<grid.Local().Size().Z(); ++index.Z())
    78         grid(index + grid.Local().Begin()) = particle_grid.GetVal(index + particle_grid.Local().Begin());
     86        grid(index + grid.Local().Begin()) = 4.0 * Math::pi * particle_grid.GetVal(index + particle_grid.Local().Begin());
    7987
    8088#ifdef DEBUG_OUTPUT
    … …  
    98106  vmg_float e_long = 0.0;
    99107  vmg_float e_self = 0.0;
     108  vmg_float e_short = 0.0;
     109  vmg_float e_shortcorr = 0.0;
    100110#endif
    101111
    … …  
    142152  for (p_iter=particles.begin(); p_iter!=particles.end(); ++p_iter) {
    143153
    144     // Interpolate long range part of potential
    145     p_iter->Pot() = 0.5 * (Helper::InterpolateTrilinear(p_iter->Pos(), particle_grid) - p_iter->Charge() * spl.GetAntiDerivativeZero());
    146 
    147 #ifdef DEBUG_OUTPUT
     154    // Interpolate long range part of potential minus self-interaction
     155    p_iter->Pot() = Helper::InterpolateTrilinear(p_iter->Pos(), particle_grid)
     156      - p_iter->Charge() * spl.GetAntiDerivativeAtZero();
     157
     158#ifdef DEBUG_OUTPUT
     159    // TODO The scaling with 1/2 may not be correct. Check that.
    148160    e_long += 0.5 * p_iter->Charge() * Helper::InterpolateTrilinear(p_iter->Pos(), particle_grid);
    149     e_self += 0.5 * p_iter->Charge() * p_iter->Charge() * spl.GetAntiDerivativeZero();
     161    e_self += 0.5 * p_iter->Charge() * p_iter->Charge() * spl.GetAntiDerivativeAtZero();
    150162#endif
    151163
    152164  }
    153 
    154   comm.CommAddPotential(particles);
    155165
    156166  /*
    … …  
    161171  Grid::iterator iter;
    162172
    163   comm.CommLCListGhosts(grid, lc);
    164 
    165   for (iter=lc.Iterators().Local().Begin(); iter!=lc.Iterators().Local().End(); ++iter)
    166     for (p1=lc(*iter).begin(); p1!=lc(*iter).end(); ++p1)
     173  comm.CommLCListToGhosts(grid, lc);
     174
     175  for (iter=lc.Iterators().Local().Begin(); iter!=lc.Iterators().Local().End(); ++iter) {
     176    std::list<Particle::Particle*>& lc_1 = lc(*iter);
     177    for (p1=lc_1.begin(); p1!=lc_1.end(); ++p1)
    167178      for (index.X()=-1*near_field_cells-1; index.X()<=near_field_cells+1; ++index.X())
    168179        for (index.Y()=-1*near_field_cells-1; index.Y()<=near_field_cells+1; ++index.Y())
    169           for (index.Z()=-1*near_field_cells-1; index.Z()<=near_field_cells+1; ++index.Z())
    170             for (p2=lc(*iter+index).begin(); p2!=lc(*iter+index).end(); ++p2)
     180          for (index.Z()=-1*near_field_cells-1; index.Z()<=near_field_cells+1; ++index.Z()) {
     181            std::list<Particle::Particle*>& lc_2 = lc(*iter+index);
     182            for (p2=lc_2.begin(); p2!=lc_2.end(); ++p2)
    171183              if (*p1 != *p2) {
    172184
    173                 length = (p2->Pos() - p1->Pos()).Length();
    174 
    175                 //TODO Rewrite this equation more efficiently
    176                 if (length < r_cut)
    177                   p1->Pot() += 0.5 * p2->Charge() / length * (0.25*M_1_PI + spl.EvaluatePotential(length));
    178 
     185                length = ((*p2)->Pos() - (*p1)->Pos()).Length();
     186
     187                //TODO Rewrite this equation more efficiently
     188                if (length < r_cut) {
     189                  (*p1)->Pot() += (*p2)->Charge() / length * (1.0 - spl.EvaluatePotential(length));
     190
     191#ifdef DEBUG_OUTPUT
     192                  e_short += 0.5 * (*p1)->Charge() * (*p2)->Charge() / length;
     193                  e_shortcorr -= 0.5 * (*p1)->Charge() * (*p2)->Charge() / length * spl.EvaluatePotential(length);
     194#endif
     195
     196                }
    179197              }
    180 
    181   comm.CommAddPotential(lc);
     198          }
     199  }
     200
     201  comm.CommParticlesBack(particles);
    182202
    183203#ifdef DEBUG_OUTPUT
    … …  
    185205
    186206  for (int i=0; i<num_particles_local; ++i)
    187     e += p[i] * q[i];
     207    e += 0.5 * p[i] * q[i];
    188208
    189209  e = comm.GlobalSumRoot(e);
    190210  e_long = comm.GlobalSumRoot(e_long);
    191211  e_self = comm.GlobalSumRoot(e_self);
     212  e_short = comm.GlobalSumRoot(e_short);
     213  e_shortcorr = comm.GlobalSumRoot(e_shortcorr);
    192214
    193215  comm.PrintStringOnce("E_long:  %e", e_long);
    194   comm.PrintStringOnce("E_short: %e", e - e_long + e_self);
     216  comm.PrintStringOnce("E_short: %e", e_short);
     217  comm.PrintStringOnce("E_shortcorr: %e", e_shortcorr);
     218  comm.PrintStringOnce("E_short*: %e", e - e_long + e_self);
    195219  comm.PrintStringOnce("E_self:  %e", e_self);
    196220  comm.PrintStringOnce("E_total: %e", e);

  • src/interface/interface_particles.hpp

    rd24c2f rac6d04  
    3131  InterfaceParticles(const Boundary& boundary, const int& levelMin, const int& levelMax,
    3232                     const Vector& box_offset, const vmg_float& box_size,
    33                      const int& near_field_cells) :
    34     Interface(boundary, levelMin, levelMax, box_offset, box_size, 0, 1.0),
     33                     const int& near_field_cells,
     34                     const int& coarsening_steps, const vmg_float& alpha) :
     35    Interface(boundary, levelMin, levelMax, box_offset, box_size, coarsening_steps, alpha),
    3536    HasRequestVec(),
    3637    spl((near_field_cells+0.5)*Extent(MaxLevel()).MeshWidth().Max())

  • src/interface/interface_particles_cf.cpp

    rd24c2f rac6d04  
    1212#ifdef HAVE_MPI
    1313#include <mpi.h>
     14#ifdef HAVE_MARMOT
     15#include <enhancempicalls.h>
     16#include <sourceinfompicalls.h>
     17#endif
    1418#endif
    1519
    … …  
    119123  comm->Register("COMM");
    120124
    121   Interface* interface = new InterfaceParticles(boundary, 2, max_level, 0.0, 1.0, near_field_cells);
     125  Interface* interface = new InterfaceParticles(boundary, 2, max_level, 0.0, 1.0, near_field_cells, 0, 1.0);
    122126  MG::SetInterface(interface, comm);
    123127
Note: See TracChangeset for help on using the changeset viewer.