Changeset 4571da


Ignore:
Timestamp:
Apr 27, 2012, 11:34:57 PM (14 years ago)
Author:
Julian Iseringhausen <isering@…>
Children:
1a92cf
Parents:
b2154a3
Message:

vmg: Implement fourth-order discretization of the Poisson equation.

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

Files:
4 added
1 deleted
32 edited

Legend:

Unmodified
Added
Removed

  • src/Makefile.am

    rb2154a3 r4571da  
    2929        base/proxy.cpp \
    3030        base/proxy.hpp \
     31        base/stencil.cpp \
    3132        base/stencil.hpp \
    3233        base/timer.cpp \
    … …  
    120121        samples/bspline.cpp \
    121122        samples/bspline.hpp \
    122         samples/discretization_poisson_fd_collatz.hpp \
     123        samples/discretization_poisson_fd.cpp \
    123124        samples/discretization_poisson_fd.hpp \
    124125        samples/discretization_poisson_fv.cpp \
    … …  
    130131        smoother/gsrb.cpp \
    131132        smoother/gsrb.hpp \
     133        smoother/jacobi.cpp \
     134        smoother/jacobi.hpp \
    132135        smoother/smoother.cpp \
    133136        smoother/smoother.hpp \

  • src/base/discretization.hpp

    rb2154a3 r4571da  
    2828{
    2929public:
    30   Discretization() : stencil(1.0) {}
    31   Discretization(const Stencil& stencil_) : stencil(stencil_) {}
     30  Discretization() :
     31    stencil(1.0),
     32    order(2)
     33  {}
     34
     35  Discretization(const int& order) :
     36    stencil(1.0),
     37    order(order)
     38  {}
     39
     40  Discretization(const Stencil& stencil_, const int& order) :
     41    stencil(stencil_),
     42    order(order)
     43  {}
    3244
    3345  Discretization(std::string id) :
    3446    Object(id),
    35     stencil(1.0)
     47    stencil(1.0),
     48    order(2)
    3649  {}
    3750
    38   Discretization(std::string id, const Stencil& stencil_) :
     51  Discretization(std::string id, const int& order) :
    3952    Object(id),
    40     stencil(stencil_)
     53    stencil(1.0),
     54    order(order)
     55  {}
     56
     57  Discretization(std::string id, const Stencil& stencil_, const int& order) :
     58    Object(id),
     59    stencil(stencil_),
     60    order(order)
    4161  {}
    4262
    … …  
    4565  virtual vmg_float OperatorPrefactor(const Grid& grid) const = 0; ///< Returns the prefactor of the operator.
    4666
    47   /**
     67  virtual void ModifyRightHandSide() {}
     68
     69  /**
    4870   * This function gets called whenever boundary points at inner boundaries are needed.
    4971   * Inner boundaries occur when using adaptive grid refinement.
    50    * 
     72   *
    5173   * @param sol_fine Solution vector / fine level
    5274   * @param rhs_fine Right handside vector / fine level
    … …  
    6082
    6183private:
    62   virtual void SetInnerBoundaryCompute(Grid& sol_fine, Grid& rhs_fine, Grid& sol_coarse) const = 0;
     84  virtual void SetInnerBoundaryCompute(Grid& sol_fine, Grid& rhs_fine, Grid& sol_coarse) const {}
    6385
    6486protected:
    6587  VMG::Stencil stencil;
     88  int order;
    6689};
    6790

  • src/base/stencil.hpp

    rb2154a3 r4571da  
    103103  }
    104104
     105  void Apply(Grid& grid) const;
     106
    105107private:
    106108  std::vector<Displacement> disp;

  • src/commands/com_import_rhs.cpp

    rb2154a3 r4571da  
    1414
    1515#include "base/command.hpp"
     16#include "base/discretization.hpp"
    1617#include "comm/comm.hpp"
    1718#include "grid/multigrid.hpp"
    … …  
    2930
    3031    MG::GetInterface()->ImportRightHandSide(*MG::GetRhs());
     32    MG::GetDiscretization()->ModifyRightHandSide();
    3133
    3234    MPE_EVENT_END()

  • src/grid/grid.cpp

    rb2154a3 r4571da  
    2525using namespace VMG;
    2626
    27 vmg_float Grid::correction;
    28 
    2927void Grid::InitGrid()
    3028{
    … …  
    112110    avg += GetVal(*iter);
    113111
    114   avg /= Local().Size().Product();
     112  avg = MG::GetComm()->GlobalSum(avg);
     113  avg /= Global().GlobalSize().Product();
     114
     115#ifdef DEBUG_OUTPUT
     116  MG::GetComm()->PrintStringOnce("Global constraint enforcement: %e", avg);
     117#endif
    115118
    116119  for (iter = Iterators().Local().Begin(); iter != Iterators().Local().End(); ++iter)
    … …  
    126129    val += GetVal(*iter);
    127130
    128   val = MG::GetComm()->GlobalSum(val);
    129 
    130   if (fabs(val) > Global().GlobalSize().Product() * std::numeric_limits<vmg_float>::epsilon()) {
     131  val = MG::GetComm()->GlobalSum(val) / Global().GlobalSize().Product();
     132
     133  if (std::abs(val) > std::numeric_limits<vmg_float>::epsilon()) {
    131134
    132135#ifdef DEBUG_OUTPUT
    … …  
    134137#endif
    135138
    136     val *= Helper::pow_2(this->Extent().MeshWidth().Max());
    137139    for (iter = Iterators().Local().Begin(); iter != Iterators().Local().End(); ++iter)
    138140      (*this)(*iter) -= val;
    … …  
    155157#endif
    156158
    157   for (Grid::iterator iter = Iterators().CompleteGrid().Begin(); iter != Iterators().CompleteGrid().End(); ++iter)
    158     (*this)(*iter) = rhs.GetVal(*iter);
     159  std::memcpy(grid, rhs.grid, local.SizeTotal().Product()*sizeof(vmg_float));
    159160}
    160161

  • src/grid/grid.hpp

    rb2154a3 r4571da  
    8181  const GridIteratorSuite& Iterators() const {return iterators;}
    8282
    83   static vmg_float& Correction() {return Grid::correction;}
    84 
    8583  void Clear();         ///< Overwrites all grid points on current level with zeros
    8684  void ClearInner();
    … …  
    9391  const vmg_float& GetVal(int x, int y, int z) const; ///< Returns the value of a requested gridpoint.
    9492  const vmg_float& GetVal(const Index& index) const;
    95 
    96   vmg_float GetCorrectedVal(int x, int y, int z) const {return this->GetVal(x, y, z) - correction;} ///< Return grid
    97   vmg_float GetCorrectedVal(const Index& index) const {return this->GetVal(index) - correction;}
    9893
    9994  void ForceDiscreteCompatibilityCondition();
    … …  
    142137
    143138  Multigrid* father;
    144 
    145   static vmg_float correction;
    146139};
    147140

  • src/grid/grid_iterator_suite.cpp

    rb2154a3 r4571da  
    2121  local.SetBounds(local_.Begin(), local_.End());
    2222  complete_grid.SetBounds(0, local_.SizeTotal());
    23   inner_local_grid.SetBounds(local_.Begin()+1, local_.End()-1);
     23  inner_local_grid.SetBounds(local_.Begin()+local_.HaloSize1(), local_.End()-local_.HaloSize2());
    2424
    2525  halo_1.X().SetBounds(Index(local_.HaloBegin1().X(), 0, 0),
    … …  
    6161
    6262  near_boundary_1.X().SetBounds(Index(local_.Begin().X(), 0, 0),
    63                                 Index(local_.Begin().X()+local_.HaloEnd1().X()-local_.HaloBegin1().X(),
     63                                Index(local_.Begin().X()+local_.HaloSize1().X(),
    6464                                      local_.SizeTotal().Y(),
    6565                                      local_.SizeTotal().Z()));
    … …  
    6767  near_boundary_1.Y().SetBounds(Index(local_.Begin().X(), local_.Begin().Y(), 0),
    6868                                Index(local_.End().X(),
    69                                       local_.Begin().Y()+local_.HaloEnd1().Y()-local_.HaloBegin1().Y(),
     69                                      local_.Begin().Y()+local_.HaloSize1().Y(),
    7070                                      local_.SizeTotal().Z()));
    7171
    … …  
    7373                                Index(local_.End().X(),
    7474                                      local_.End().Y(),
    75                                       local_.Begin().Z()+local_.HaloEnd1().Z()-local_.HaloBegin1().Z()));
     75                                      local_.Begin().Z()+local_.HaloSize1().Z()));
    7676
    77   near_boundary_2.X().SetBounds(Index(local_.End().X()-local_.HaloEnd2().X()+local_.HaloBegin2().X(),
     77  near_boundary_2.X().SetBounds(Index(local_.End().X()-local_.HaloSize2().X(),
    7878                                      0,
    7979                                      0),
    … …  
    8181
    8282  near_boundary_2.Y().SetBounds(Index(local_.Begin().X(),
    83                                       local_.End().Y()-local_.HaloEnd2().Y()+local_.HaloBegin2().Y(),
     83                                      local_.End().Y()-local_.HaloSize2().Y(),
    8484                                      0),
    8585                                Index(local_.End().X(), local_.End().Y(), local_.SizeTotal().Z()));
    … …  
    8787  near_boundary_2.Z().SetBounds(Index(local_.Begin().X(),
    8888                                      local_.Begin().Y(),
    89                                       local_.End().Z()-local_.HaloEnd2().Z()+local_.HaloBegin2().Z()),
     89                                      local_.End().Z()-local_.HaloSize2().Z()),
    9090                                local_.End());
    9191}

  • src/interface/interface_fcs.cpp

    rb2154a3 r4571da  
    3535#include "samples/stencils.hpp"
    3636#include "samples/techniques.hpp"
    37 #include "smoother/gsrb.cpp"
     37#include "smoother/gsrb.hpp"
    3838#ifdef HAVE_LAPACK
    3939#include "solver/dgesv.hpp"
    … …  
    5959  static vmg_int near_field_cells = -1;
    6060  static vmg_int interpolation_degree = -1;
     61  static vmg_int discretization_order = -1;
    6162  static MPI_Comm mpi_comm;
    6263}
    … …  
    6667                         vmg_float* box_offset, vmg_float box_size,
    6768                         vmg_int near_field_cells, vmg_int interpolation_degree,
    68                          MPI_Comm mpi_comm)
     69                         vmg_int discretization_order, MPI_Comm mpi_comm)
    6970{
    7071  VMGBackupSettings::level = level;
    … …  
    7879  VMGBackupSettings::near_field_cells = near_field_cells;
    7980  VMGBackupSettings::interpolation_degree = interpolation_degree;
     81  VMGBackupSettings::discretization_order = discretization_order;
    8082  VMGBackupSettings::mpi_comm = mpi_comm;
    8183
    … …  
    113115  Discretization* discretization;
    114116  if (singular)
    115     discretization = new DiscretizationPoissonFD();
    116   else
    117     discretization = new DiscretizationPoissonFV();
     117    discretization = new DiscretizationPoissonFD(discretization_order);
     118  else
     119    discretization = new DiscretizationPoissonFV(discretization_order);
    118120  discretization->Register("DISCRETIZATION");
    119121
    … …  
    154156   */
    155157  if (singular)
    156     Techniques::SetCorrectionSchemePeriodicParticle(interface->MinLevel(), interface->MaxLevel(), gamma);
     158    Techniques::SetCorrectionSchemePeriodic(interface->MinLevel(), interface->MaxLevel(), gamma);
    157159  else
    158160    Techniques::SetFullApproximationSchemeDirichlet(interface->MinLevel(), interface->MaxLevel(), gamma);
    … …  
    183185                   vmg_float* box_offset, vmg_float box_size,
    184186                   vmg_int near_field_cells, vmg_int interpolation_degree,
    185                    MPI_Comm mpi_comm)
     187                   vmg_int discretization_order, MPI_Comm mpi_comm)
    186188{
    187189  if (VMGBackupSettings::level != level ||
    … …  
    199201      VMGBackupSettings::near_field_cells != near_field_cells ||
    200202      VMGBackupSettings::interpolation_degree != interpolation_degree ||
     203      VMGBackupSettings::discretization_order != discretization_order ||
    201204      VMGBackupSettings::mpi_comm != mpi_comm) {
    202205
    … …  
    205208                 smoothing_steps, gamma, precision,
    206209                 box_offset, box_size, near_field_cells,
    207                  interpolation_degree, mpi_comm);
     210                 interpolation_degree, discretization_order,
     211                 mpi_comm);
    208212
    209213  }

  • src/interface/interface_fcs.h

    rb2154a3 r4571da  
    1717                   fcs_float* box_offset, fcs_float box_size,
    1818                   fcs_int near_field_cells, fcs_int interpolation_degree,
    19                    MPI_Comm mpi_comm);
     19                   fcs_int discretization_order, MPI_Comm mpi_comm);
    2020
    2121int VMG_fcs_check();

  • src/interface/interface_particles.cpp

    rb2154a3 r4571da  
    2020#endif
    2121
     22#include <algorithm>
    2223#include <cmath>
    2324#include <cstring>
    … …  
    8586    for (index.Y()=0; index.Y()<grid.Local().Size().Y(); ++index.Y())
    8687      for (index.Z()=0; index.Z()<grid.Local().Size().Z(); ++index.Z())
    87         grid(index + grid.Local().Begin()) = 4.0 * Math::pi * particle_grid.GetVal(index + particle_grid.Local().Begin());
     88        grid(index + grid.Local().Begin()) = 4.0 * Math::pi * particle_grid.GetVal(index + particle_grid.Local().Begin());
    8889
    8990#ifdef DEBUG_OUTPUT
    … …  
    99100void InterfaceParticles::ExportSolution(Grid& grid)
    100101{
    101   Index index;
     102  Index i;
    102103  vmg_float length;
    103   Vector dist_vec;
    104104
    105105#ifdef DEBUG_OUTPUT
    … …  
    107107  vmg_float e_long = 0.0;
    108108  vmg_float e_self = 0.0;
    109   vmg_float e_short = 0.0;
    110   vmg_float e_shortcorr = 0.0;
     109  vmg_float e_short_peak = 0.0;
     110  vmg_float e_short_spline = 0.0;
    111111#endif
    112112
    … …  
    125125  InterpolatePolynomial ip(interpolation_degree);
    126126
    127   const vmg_float r_cut = (near_field_cells+0.5) * grid.Extent().MeshWidth().Max();
    128 
    129   /*
    130    * Initialize potential
    131    */
    132   for (vmg_int i=0; i<num_particles_local; ++i) {
    133     p[i] = 0.0;
    134     for (vmg_int j=0;j<3;++j)
    135       f[3*i+j] = 0.;
    136   }
     127  const vmg_float r_cut = near_field_cells * grid.Extent().MeshWidth().Max();
     128
     129  // Initialize field
     130  std::fill(f, f+3*num_particles_local, 0.0);
    137131
    138132  /*
    … …  
    143137  Grid& particle_grid = comm.GetParticleGrid();
    144138
    145   for (index.X()=0; index.X()<grid.Local().Size().X(); ++index.X())
    146     for (index.Y()=0; index.Y()<grid.Local().Size().Y(); ++index.Y())
    147       for (index.Z()=0; index.Z()<grid.Local().Size().Z(); ++index.Z())
    148         particle_grid(index + particle_grid.Local().Begin()) = grid.GetVal(index + grid.Local().Begin());
     139  for (i.X()=0; i.X()<grid.Local().Size().X(); ++i.X())
     140    for (i.Y()=0; i.Y()<grid.Local().Size().Y(); ++i.Y())
     141      for (i.Z()=0; i.Z()<grid.Local().Size().Z(); ++i.Z())
     142        particle_grid(i + particle_grid.Local().Begin()) = grid.GetVal(i + grid.Local().Begin());
    149143
    150144  comm.CommToGhosts(particle_grid);
    … …  
    170164
    171165#ifdef DEBUG_OUTPUT
    172     // TODO The scaling with 1/2 may not be correct. Check that.
    173       e_long += 0.5 * (*p1)->Charge() * Helper::InterpolateTrilinear((*p1)->Pos(), particle_grid);
     166      e_long += 0.5 * (*p1)->Charge() * ip.Evaluate((*p1)->Pos());
    174167      e_self += 0.5 * (*p1)->Charge() * (*p1)->Charge() * spl.GetAntiDerivativeAtZero();
    175168#endif
    176169
    177       for (index.X()=-1*near_field_cells-1; index.X()<=near_field_cells+1; ++index.X())
    178         for (index.Y()=-1*near_field_cells-1; index.Y()<=near_field_cells+1; ++index.Y())
    179           for (index.Z()=-1*near_field_cells-1; index.Z()<=near_field_cells+1; ++index.Z()) {
    180 
    181             std::list<Particle::Particle*>& lc_2 = lc(*iter+index);
     170      for (i.X()=-1*near_field_cells; i.X()<=near_field_cells; ++i.X())
     171        for (i.Y()=-1*near_field_cells; i.Y()<=near_field_cells; ++i.Y())
     172          for (i.Z()=-1*near_field_cells; i.Z()<=near_field_cells; ++i.Z()) {
     173
     174            std::list<Particle::Particle*>& lc_2 = lc(*iter+i);
    182175
    183176            for (p2=lc_2.begin(); p2!=lc_2.end(); ++p2)
    … …  
    186179                length = ((*p2)->Pos() - (*p1)->Pos()).Length();
    187180
    188                 //TODO Rewrite this equation more efficiently
    189181                if (length < r_cut) {
    190182                  (*p1)->Pot() += (*p2)->Charge() / length * (1.0 + spl.EvaluatePotential(length));
    191183
    192184#ifdef DEBUG_OUTPUT
    193                   e_short += 0.5 * (*p1)->Charge() * (*p2)->Charge() / length;
    194                   e_shortcorr -= 0.5 * (*p1)->Charge() * (*p2)->Charge() / length * spl.EvaluatePotential(length);
     185                  e_short_peak += 0.5 * (*p1)->Charge() * (*p2)->Charge() / length;
     186                  e_short_spline += 0.5 * (*p1)->Charge() * (*p2)->Charge() / length * spl.EvaluatePotential(length);
    195187#endif
    196188                }
    … …  
    205197  vmg_float* q = factory.GetObjectStorageArray<vmg_float>("PARTICLE_CHARGE_ARRAY");
    206198
    207   for (int i=0; i<num_particles_local; ++i)
    208     e += 0.5 * p[i] * q[i];
    209 
     199  e_long = comm.GlobalSumRoot(e_long);
     200  e_short_peak = comm.GlobalSumRoot(e_short_peak);
     201  e_short_spline = comm.GlobalSumRoot(e_short_spline);
     202  e_self = comm.GlobalSumRoot(e_self);
     203
     204  for (int j=0; j<num_particles_local; ++j)
     205    e += 0.5 * p[j] * q[j];
    210206  e = comm.GlobalSumRoot(e);
    211   e_long = comm.GlobalSumRoot(e_long);
    212   e_self = comm.GlobalSumRoot(e_self);
    213   e_short = comm.GlobalSumRoot(e_short);
    214   e_shortcorr = comm.GlobalSumRoot(e_shortcorr);
    215 
    216   comm.PrintStringOnce("E_long:  %e", e_long);
    217   comm.PrintStringOnce("E_short: %e", e_short);
    218   comm.PrintStringOnce("E_shortcorr: %e", e_shortcorr);
    219   comm.PrintStringOnce("E_short*: %e", e - e_long + e_self);
    220   comm.PrintStringOnce("E_self:  %e", e_self);
    221   comm.PrintStringOnce("E_total: %e", e);
     207
     208  comm.PrintStringOnce("E_long:         %e", e_long);
     209  comm.PrintStringOnce("E_short_peak:   %e", e_short_peak);
     210  comm.PrintStringOnce("E_short_spline: %e", e_short_spline);
     211  comm.PrintStringOnce("E_self:         %e", e_self);
     212  comm.PrintStringOnce("E_total:        %e", e);
     213  comm.PrintStringOnce("E_total*:       %e", e_long + e_short_peak + e_short_spline - e_self);
    222214
    223215#endif /* DEBUG_OUTPUT */

  • src/interface/interface_particles.hpp

    rb2154a3 r4571da  
    3535    Interface(boundary, levelMin, levelMax, box_offset, box_size, coarsening_steps, alpha),
    3636    HasRequestVec(),
    37     spl((near_field_cells+0.5)*Extent(MaxLevel()).MeshWidth().Max())
     37    spl(near_field_cells*Extent(MaxLevel()).MeshWidth().Max())
    3838  {}
    3939

  • src/interface/interface_particles_cf.cpp

    rb2154a3 r4571da  
    5959  int max_level, max_iterations;
    6060  int pre_smoothing_steps, post_smoothing_steps;
    61   int gamma, near_field_cells;
     61  int gamma, near_field_cells, discretization_order;
    6262  vmg_float precision;
    6363  std::string lop_str, datafile_str;
    … …  
    109109  gamma = Helper::ToValWithDefault<int>(xml_conf.child_value("gamma"), 1);
    110110  near_field_cells = Helper::ToValWithDefault<int>(xml_conf.child_value("near_field_cells"), 3);
     111  discretization_order = Helper::ToValWithDefault<int>(xml_conf.child_value("discretization_order"), 2);
    111112  precision = Helper::ToValWithDefault<vmg_float>(xml_conf.child_value("precision"), 1.0e-7);
    112113  lop_str = xml_conf.child_value("level_operator");
    … …  
    126127  MG::SetInterface(interface, comm);
    127128
    128   Discretization* discretization = new DiscretizationPoissonFD();
     129  Discretization* discretization = new DiscretizationPoissonFD(discretization_order);
    129130  discretization->Register("DISCRETIZATION");
    130131
    … …  
    174175    std::printf("  Precision:                    %e\n", precision);
    175176    std::printf("  Near field cells:             %d\n", near_field_cells);
     177    std::printf("  Discretization order:         %d\n", discretization_order);
    176178  }
    177179

  • src/mg.cpp

    rb2154a3 r4571da  
    223223}
    224224
     225VMG::Grid& MG::GetRhsMaxLevel()
     226{
     227  return (*MG::GetRhs())(MG::GetRhs()->MaxLevel());
     228}
     229
     230VMG::Grid& MG::GetSolMaxLevel()
     231{
     232  return (*MG::GetSol())(MG::GetSol()->MaxLevel());
     233}
     234
    225235Smoother* MG::GetSmoother()
    226236{

  • src/mg.hpp

    rb2154a3 r4571da  
    4747  static VMG::Multigrid* GetRhs();
    4848  static VMG::Multigrid* GetSol();
     49  static VMG::Grid& GetRhsMaxLevel();
     50  static VMG::Grid& GetSolMaxLevel();
    4951  static VMG::Smoother* GetSmoother();
    5052  static VMG::Solver* GetSolver();

  • src/samples/bspline.hpp

    rb2154a3 r4571da  
    2626  BSpline(const vmg_float& width);
    2727
    28   vmg_float EvaluateSpline(const vmg_float& val)
     28  vmg_float EvaluateSpline(const vmg_float& val) const
    2929  {
    3030    for (unsigned int i=0; i<intervals.size(); ++i)
    … …  
    3434  }
    3535
    36   void SetSpline(Grid& grid, const Particle::Particle p, const int& near_field_cells)
     36  void SetSpline(Grid& grid, const Particle::Particle p, const int& near_field_cells) const
    3737  {
    3838    assert(p.Pos().X() >= grid.Extent().Begin().X() && p.Pos().X() < grid.Extent().End().X());
    … …  
    8383  }
    8484
    85   vmg_float EvaluatePotential(const vmg_float& val)
     85  vmg_float EvaluatePotential(const vmg_float& val) const
    8686  {
    8787    for (unsigned int i=0; i<intervals.size(); ++i)

  • src/samples/discretization_poisson_fd.hpp

    rb2154a3 r4571da  
    2121{
    2222public:
    23   DiscretizationPoissonFD()
    24   {
    25     stencil.SetDiag(6.0);
    26     stencil.push_back(-1,  0,  0, -1.0);
    27     stencil.push_back( 1,  0,  0, -1.0);
    28     stencil.push_back( 0, -1,  0, -1.0);
    29     stencil.push_back( 0,  1,  0, -1.0);
    30     stencil.push_back( 0,  0, -1, -1.0);
    31     stencil.push_back( 0,  0,  1, -1.0);
    32   }
     23  DiscretizationPoissonFD(const int& order);
    3324
    3425  vmg_float OperatorPrefactor(const Grid& grid) const
    … …  
    3728  }
    3829
    39 private:
    40   void SetInnerBoundaryCompute(Grid& sol_fine, Grid& rhs_fine, Grid& sol_coarse) const {}
     30  void ModifyRightHandSide();
    4131};
    4232

  • src/samples/discretization_poisson_fv.hpp

    rb2154a3 r4571da  
    2323{
    2424public:
    25   DiscretizationPoissonFV()
     25  DiscretizationPoissonFV(const int& order) :
     26    Discretization(order)
    2627  {
    2728    stencil.SetDiag(6.0);

  • src/samples/techniques.hpp

    rb2154a3 r4571da  
    108108    AddCycleGamma(*loop, maxLevel-minLevel+1, gamma);
    109109
     110    loop->AddCommand("SetAverageToZero", "SOL");
    110111    loop->AddCommand("ComputeResidualNorm", "RESIDUAL");
    111112    loop->AddCommand("CheckResidual", "RESIDUAL");
    … …  
    140141    AddCycleGammaDebug(*loop, maxLevel-minLevel+1, gamma);
    141142
     143    loop->AddCommand("SetAverageToZero", "SOL");
    142144    loop->AddCommand("ComputeResidualNorm", "RESIDUAL");
    143145    loop->AddCommand("CheckResidual", "RESIDUAL");
    … …  
    170172    AddCycleGamma(*loop, maxLevel-minLevel+1, gamma);
    171173
     174    loop->AddCommand("SetAverageToZero", "SOL");
    172175    loop->AddCommand("ComputeResidualNorm", "RESIDUAL");
    173176    loop->AddCommand("CheckResidual", "RESIDUAL");
    … …  
    201204    AddCycleGammaDebug(*loop, maxLevel-minLevel+1, gamma);
    202205
     206    loop->AddCommand("SetAverageToZero", "SOL");
    203207    loop->AddCommand("ComputeResidualNorm", "RESIDUAL");
    204208    loop->AddCommand("CheckResidual", "RESIDUAL");

  • src/smoother/gs.cpp

    rb2154a3 r4571da  
    4343  for (grid_iter = rhs.Iterators().Local().Begin(); grid_iter != rhs.Iterators().Local().End(); ++grid_iter) {
    4444
    45     temp = prefactor_inv * rhs.GetCorrectedVal(*grid_iter);
     45    temp = prefactor_inv * rhs.GetVal(*grid_iter);
    4646
    4747    for (stencil_iter=A.begin(); stencil_iter!=A.end(); ++stencil_iter)

  • src/smoother/gsrb.cpp

    rb2154a3 r4571da  
    1212#endif
    1313
     14#ifdef HAVE_MPI
    1415#include <mpi.h>
    15 
    16 #include <sstream>
     16#endif
    1717
    1818#include "base/discretization.hpp"

  • src/solver/solver_singular.cpp

    rb2154a3 r4571da  
    123123
    124124      }
    125 
    126   // Set correction
    127   Grid::Correction() = this->Sol(this->Size()-1);
    128125}

  • test/unit_test/library/dirichlet_cs.cpp

    rb2154a3 r4571da  
    4848    MG::SetInterface(interface, comm);
    4949
    50     Discretization* discretization = new DiscretizationPoissonFD();
     50    Discretization* discretization = new DiscretizationPoissonFD(2);
    5151    discretization->Register("DISCRETIZATION");
    5252

  • test/unit_test/library/dirichlet_cs_mpi.cpp

    rb2154a3 r4571da  
    5656    MG::SetInterface(interface, comm);
    5757
    58     Discretization* discretization = new DiscretizationPoissonFD();
     58    Discretization* discretization = new DiscretizationPoissonFD(2);
    5959    discretization->Register("DISCRETIZATION");
    6060

  • test/unit_test/library/dirichlet_fas.cpp

    rb2154a3 r4571da  
    4848    MG::SetInterface(interface, comm);
    4949
    50     Discretization* discretization = new DiscretizationPoissonFD();
     50    Discretization* discretization = new DiscretizationPoissonFD(2);
    5151    discretization->Register("DISCRETIZATION");
    5252

  • test/unit_test/library/dirichlet_fas_lr.cpp

    rb2154a3 r4571da  
    4848    MG::SetInterface(interface, comm);
    4949
    50     Discretization* discretization = new DiscretizationPoissonFV();
     50    Discretization* discretization = new DiscretizationPoissonFV(2);
    5151    discretization->Register("DISCRETIZATION");
    5252

  • test/unit_test/library/dirichlet_fas_lr_mpi.cpp

    rb2154a3 r4571da  
    5757    MG::SetInterface(interface, comm);
    5858
    59     Discretization* discretization = new DiscretizationPoissonFV();
     59    Discretization* discretization = new DiscretizationPoissonFV(2);
    6060    discretization->Register("DISCRETIZATION");
    6161

  • test/unit_test/library/dirichlet_fas_mpi.cpp

    rb2154a3 r4571da  
    5757    MG::SetInterface(interface, comm);
    5858
    59     Discretization* discretization = new DiscretizationPoissonFD();
     59    Discretization* discretization = new DiscretizationPoissonFD(2);
    6060    discretization->Register("DISCRETIZATION");
    6161

  • test/unit_test/library/periodic_cs.cpp

    rb2154a3 r4571da  
    4848    MG::SetInterface(interface, comm);
    4949
    50     Discretization* discretization = new DiscretizationPoissonFD();
     50    Discretization* discretization = new DiscretizationPoissonFD(2);
    5151    discretization->Register("DISCRETIZATION");
    5252

  • test/unit_test/library/periodic_cs_mpi.cpp

    rb2154a3 r4571da  
    5656    MG::SetInterface(interface, comm);
    5757
    58     Discretization* discretization = new DiscretizationPoissonFD();
     58    Discretization* discretization = new DiscretizationPoissonFD(2);
    5959    discretization->Register("DISCRETIZATION");
    6060

  • test/unit_test/library/periodic_fas.cpp

    rb2154a3 r4571da  
    4747    MG::SetInterface(interface, comm);
    4848
    49     Discretization* discretization = new DiscretizationPoissonFD();
     49    Discretization* discretization = new DiscretizationPoissonFD(2);
    5050    discretization->Register("DISCRETIZATION");
    5151

  • test/unit_test/library/periodic_fas_mpi.cpp

    rb2154a3 r4571da  
    5656  MG::SetInterface(interface, comm);
    5757
    58   Discretization* discretization = new DiscretizationPoissonFD();
     58  Discretization* discretization = new DiscretizationPoissonFD(2);
    5959  discretization->Register("DISCRETIZATION");
    6060

  • test/unit_test/unit_test/smoother_test.cpp

    rb2154a3 r4571da  
    3434    comm->Register("COMM");
    3535
    36     Discretization* discretization = new DiscretizationPoissonFD();
     36    Discretization* discretization = new DiscretizationPoissonFD(2);
    3737    discretization->Register("DISCRETIZATION");
    3838
Note: See TracChangeset for help on using the changeset viewer.