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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/level/level_operator_cs.cpp

    r66f24d rdfed1c  
    1414#include "base/index.hpp"
    1515#include "comm/comm.hpp"
     16#include "grid/grid_iterator.hpp"
    1617#include "grid/multigrid.hpp"
    1718#include "grid/tempgrid.hpp"
    … …  
    2122using namespace VMG;
    2223
    23 void LevelOperatorCS::Restrict(Multigrid& sol, Multigrid& rhs)
     24void LevelOperatorCS::Restrict(Grid& sol_f, Grid& rhs_f, Grid& sol_c, Grid& rhs_c)
    2425{
    25 #ifdef DEBUG
    26   sol().IsConsistent();
    27   rhs().IsConsistent();
    28   sol().IsCompatible(rhs());
     26#ifdef DEBUG_MATRIX_CHECKS
     27  sol_f.IsConsistent();
     28  rhs_f.IsConsistent();
    2929#endif
    3030
    31   Index i_c, i_f;
    32   Stencil::iterator iter;
     31#ifdef DEBUG
     32  sol_f.IsCompatible(rhs_f);
     33#endif
     34
     35  Grid::iterator iter_f, iter_c;
     36  Stencil::iterator stencil_iter;
    3337  vmg_float res;
    3438
    35   const int l_f = rhs.Level();
    36   const int l_c = l_f - 1;
    37 
    38   Index begin_f = rhs().Local().Begin();
    39   const Index end_f = rhs().Local().End() - 1;
    40   const Index begin_c = rhs(l_c).Local().AlignmentBegin();
    4139  const Stencil& op = OperatorRestrict();
    4240
    43   if (sol(l_c).Global().BoundaryType() == GlobalCoarsened &&
    44      (MG::GetComm()->BoundaryConditions() == Dirichlet || MG::GetComm()->BoundaryConditions() == Quasiperiodic))
    45     begin_f += 1;
     41  Index begin_f = rhs_f.Local().Begin() + 2*rhs_c.Global().BeginLocal() - rhs_f.Global().BeginLocal();
     42  Index end_f = rhs_f.Local().End();
    4643
    47   // Communication step
    48   MG::GetComm()->CommToGhosts(sol());
     44  /* Modify fine begin/end to align the points on both levels correctly */
     45  if (rhs_c.Global().BoundaryType() == GlobalCoarsened)
     46    for (int i=0; i<3; ++i) {
     47
     48      if (rhs_c.Local().HasBoundary1()[i])
     49        ++begin_f[i];
     50
     51      if (rhs_c.Local().HasBoundary2()[i])
     52        --end_f[i];
     53
     54    }
     55
     56  const GridIteratorSet bounds_f(begin_f, end_f);
     57  const GridIteratorSet bounds_c(rhs_c.Local().AlignmentBegin(), rhs_c.Local().AlignmentEnd());
    4958
    5059  // Compute residual
    5160  TempGrid *temp = MG::GetTempGrid();
    52   temp->SetProperties(rhs());
    53   temp->ImportFromResidual(sol(), rhs());
     61  temp->SetProperties(sol_f);
     62  temp->ImportFromResidual(sol_f, rhs_f);
    5463
    55   // Bring grids to the next coarser level
    56   rhs.ToCoarserLevel();
    57   sol.ToCoarserLevel();
     64  if (rhs_c.Global().BoundaryType() == GlobalCoarsened)
     65    rhs_c.ClearBoundary();
     66
     67  for (iter_f=bounds_f.Begin(), iter_c=bounds_c.Begin(); iter_c!=bounds_c.End(); iter_f+=2, ++iter_c) {
     68
     69    res = op.GetDiag() * temp->GetVal(*iter_f);
     70
     71    for (stencil_iter=op.begin(); stencil_iter!=op.end(); ++stencil_iter)
     72      res += stencil_iter->Val() * temp->GetVal(*iter_f + stencil_iter->Disp());
     73
     74    rhs_c(*iter_c) = res;
     75
     76  }
     77
     78#ifdef DEBUG_MATRIX_CHECKS
     79  rhs_c.IsConsistent();
     80#endif
     81}
     82
     83void LevelOperatorCS::Prolongate(Grid& sol_f, Grid& rhs_f, Grid& sol_c, Grid& rhs_c)
     84{
     85#ifdef DEBUG_MATRIX_CHECKS
     86  sol_c.IsConsistent();
     87  rhs_c.IsConsistent();
     88  sol_f.IsConsistent();
     89  rhs_f.IsConsistent();
     90#endif
    5891
    5992#ifdef DEBUG
    60   sol().IsConsistent();
    61   rhs().IsConsistent();
    62   sol().IsCompatible(rhs());
     93  sol_c.IsCompatible(rhs_c);
     94  sol_f.IsCompatible(rhs_f);
    6395#endif
    6496
    65   if (rhs().Global().BoundaryType() == GlobalCoarsened)
    66     rhs().ClearBoundary();
    67 
    68   // Set coarse grid values according to the given stencil
    69   for (i_c.X()=begin_c.X(), i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_c.X(), i_f.X()+=2)
    70     for (i_c.Y()=begin_c.Y(), i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_c.Y(), i_f.Y()+=2)
    71       for (i_c.Z()=begin_c.Z(), i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_c.Z(), i_f.Z()+=2) {
    72 
    73         res = op.GetDiag() * temp->GetVal(i_f);
    74 
    75         for (iter=op.begin(); iter!=op.end(); ++iter)
    76           res += iter->val * temp->GetVal(i_f + iter);
    77 
    78         rhs()(i_c) = res;
    79 
    80       }
    81 }
    82 
    83 void LevelOperatorCS::Prolongate(Multigrid& sol, Multigrid& rhs)
    84 {
    85 #ifdef DEBUG
    86   sol().IsConsistent();
    87   rhs().IsConsistent();
    88   sol().IsCompatible(rhs());
    89 #endif
    90 
    91   Stencil::iterator iter;
    92   Index i_c, i_f;
     97  Grid::iterator iter_f, iter_c;
     98  Stencil::iterator stencil_iter;
    9399  vmg_float val;
    94100
    95101  const Stencil& op = OperatorProlongate();
     102  Comm* comm = MG::GetComm();
    96103
    97   const int l_c = sol.Level();
    98   const int l_f = l_c + 1;
     104  Index begin_f = sol_f.Local().Begin() + 2*sol_c.Global().BeginLocal() - sol_f.Global().BeginLocal();
     105  Index end_f = sol_f.Local().End();
    99106
    100   const Index begin_c = sol().Local().AlignmentBegin();
    101   Index begin_f = sol(l_f).Local().Begin();
    102   const Index end_f = sol(l_f).Local().End();
     107  if (sol_c.Global().BoundaryType() == GlobalCoarsened)
     108    for (int i=0; i<3; ++i) {
    103109
    104   if (sol(l_c).Global().BoundaryType() == GlobalCoarsened &&
    105       (MG::GetComm()->BoundaryConditions() == Dirichlet || MG::GetComm()->BoundaryConditions() == Quasiperiodic))
    106     begin_f += 1;
     110      if (rhs_f.Local().HasBoundary1()[i])
     111        ++begin_f[i];
    107112
    108   sol.ToFinerLevel();
    109   rhs.ToFinerLevel();
     113      if (rhs_f.Local().HasBoundary2()[i])
     114        --end_f[i];
    110115
    111 #ifdef DEBUG
    112   sol().IsConsistent();
    113   rhs().IsConsistent();
    114   sol().IsCompatible(rhs());
    115 #endif
     116    }
    116117
    117   if (MG::GetComm()->BoundaryConditions() == Periodic) {
     118  const GridIteratorSet bounds_f(begin_f, end_f);
     119  const GridIteratorSet bounds_c(sol_c.Local().AlignmentBegin(), sol_c.Local().AlignmentEnd());
    118120
    119     sol().ClearHalo();
     121  Index apply_stencil_begin, apply_stencil_end;
    120122
    121     for (i_c.X()=begin_c.X(), i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_c.X()++, i_f.X()+=2)
    122       for (i_c.Y()=begin_c.Y(), i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_c.Y()++, i_f.Y()+=2)
    123         for (i_c.Z()=begin_c.Z(), i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_c.Z()++, i_f.Z()+=2) {
     123  for (int i=0; i<3; ++i) {
     124    apply_stencil_begin[i] = (sol_f.Local().HasBoundary1()[i] ? 1 : 0);
     125    apply_stencil_end[i] = (sol_f.Local().HasBoundary2()[i] ? sol_f.Local().End()[i] : sol_f.Local().SizeTotal()[i]);
     126  }
    124127
    125           val = sol(l_c).GetVal(i_c);
     128  sol_f.ClearHalo();
    126129
    127           sol()(i_f) += op.GetDiag() * val;
     130  for (iter_f=bounds_f.Begin(), iter_c=bounds_c.Begin(); iter_c!=bounds_c.End(); iter_f+=2, ++iter_c) {
    128131
    129           for (iter=op.begin(); iter!=op.end(); ++iter)
    130             sol()(i_f + iter) += iter->val * val;
     132    val = sol_c.GetVal(*iter_c);
    131133
    132         }
     134    sol_f(*iter_f) += op.GetDiag() * val;
    133135
    134     MG::GetComm()->CommFromGhosts(sol());
     136    for (stencil_iter = op.begin(); stencil_iter != op.end(); ++stencil_iter) {
    135137
    136   }else if (MG::GetComm()->BoundaryConditions() == Dirichlet ||
    137             MG::GetComm()->BoundaryConditions() == Quasiperiodic) {
     138      sol_f(*iter_f + stencil_iter->Disp()) += stencil_iter->Val() * val;
    138139
    139     for (i_c.X()=begin_c.X(), i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_c.X(), i_f.X()+=2)
    140       for (i_c.Y()=begin_c.Y(), i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_c.Y(), i_f.Y()+=2)
    141         for (i_c.Z()=begin_c.Z(), i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_c.Z(), i_f.Z()+=2) {
     140    }
    142141
    143           val = sol(l_c).GetVal(i_c);
     142  }
    144143
    145           sol()(i_f) += op.GetDiag() * val;
    146 
    147           for (iter=op.begin(); iter!=op.end(); ++iter)
    148             if (i_f.X()+iter->m >= sol().Local().Begin().X() && i_f.X()+iter->m < sol().Local().End().X() &&
    149                 i_f.Y()+iter->n >= sol().Local().Begin().Y() && i_f.Y()+iter->n < sol().Local().End().Y() &&
    150                 i_f.Z()+iter->o >= sol().Local().Begin().Z() && i_f.Z()+iter->o < sol().Local().End().Z()) {
    151               sol()(i_f + iter) += iter->val * val;
    152             }
    153         }
    154   }
     144  comm->CommFromGhosts(sol_f);
    155145}
Note: See TracChangeset for help on using the changeset viewer.