Changeset ac6d04 for src/level/level_operator_cs.cpp

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

File:

• src/level/level_operator_cs.cpp (modified) (2 diffs)

src/level/level_operator_cs.cpp

@@ -22 +22 @@
 using namespace VMG;
 
-void LevelOperatorCS::Restrict(Grid& sol_f, Grid& rhs_f, Grid& sol_c, Grid& rhs_c)
+void LevelOperatorCS::Restrict(Multigrid& sol, Multigrid& rhs)
 {
-#ifdef DEBUG_MATRIX_CHECKS
-  sol_f.IsConsistent();
-  rhs_f.IsConsistent();
-#endif
+  Grid::iterator iter_f, iter_c;
 
-#ifdef DEBUG
-  sol_f.IsCompatible(rhs_f);
-#endif
+  Comm& comm = *MG::GetComm();
 
-  Grid::iterator iter_f, iter_c;
-  Stencil::iterator stencil_iter;
-  vmg_float res;
+  Grid& sol_f = sol(sol.Level());
+  Grid& rhs_f = rhs(rhs.Level());
+
+  Grid& sol_c_dist = sol(sol.Level()-1);
+  Grid& rhs_c_dist = rhs(rhs.Level()-1);
+
+  Grid& sol_c_undist = comm.GetCoarserGrid(sol);
+  Grid& rhs_c_undist = comm.GetCoarserGrid(rhs);
 
   const Stencil& op = OperatorRestrict();
 
-  Index begin_f = rhs_f.Local().Begin() + 2*rhs_c.Global().BeginLocal() - rhs_f.Global().BeginLocal();
-  Index end_f = rhs_f.Local().End();
+  const Index begin_f_global = rhs_f.Global().LocalBegin() + rhs_f.Global().LocalBegin() % 2;
+  const Index end_f_global = rhs_f.Global().LocalEnd() - (rhs_f.Global().LocalEnd() - 1) % 2;
+
+  assert(begin_f_global % 2 == 0);
+  assert(end_f_global % 2 == 1);
+
+  Index begin_f = begin_f_global - rhs_f.Global().LocalBegin() + rhs_f.Local().Begin();
+  Index end_f = end_f_global - rhs_f.Global().LocalBegin() + rhs_f.Local().Begin();
 
   /* Modify fine begin/end to align the points on both levels correctly */
-  if (rhs_c.Global().BoundaryType() == GlobalCoarsened)
-    for (int i=0; i<3; ++i) {
+  if (rhs_c_undist.Global().BoundaryType() == GlobalCoarsened) {
+    begin_f += rhs_f.Local().BoundarySize1();
+    end_f -= rhs_f.Local().BoundarySize2();
+  }
 
-      if (rhs_c.Local().HasBoundary1()[i])
-        ++begin_f[i];
-
-      if (rhs_c.Local().HasBoundary2()[i])
-        --end_f[i];
-
-    }
+  const Index begin_c = rhs_c_undist.Local().Begin() + begin_f_global / 2 - rhs_c_undist.Global().LocalBegin();
+  const Index end_c = rhs_c_undist.Local().Begin() + end_f_global / 2 - rhs_c_undist.Global().LocalBegin() + 1;
 
   const GridIteratorSet bounds_f(begin_f, end_f);
-  const GridIteratorSet bounds_c(rhs_c.Local().AlignmentBegin(), rhs_c.Local().AlignmentEnd());
+  const GridIteratorSet bounds_c(begin_c, end_c);
 
   // Compute residual
@@ -61 +64 @@
   temp->SetProperties(sol_f);
   temp->ImportFromResidual(sol_f, rhs_f);
+  comm.CommToGhosts(*temp);
 
-  if (rhs_c.Global().BoundaryType() == GlobalCoarsened)
-    rhs_c.ClearBoundary();
+  for (iter_f=bounds_f.Begin(), iter_c=bounds_c.Begin(); iter_f!=bounds_f.End(); iter_f+=2, ++iter_c)
+    rhs_c_undist(*iter_c) = op.Apply(*temp, *iter_f);
 
-  for (iter_f=bounds_f.Begin(), iter_c=bounds_c.Begin(); iter_c!=bounds_c.End(); iter_f+=2, ++iter_c) {
+  comm.CommSubgrid(rhs_c_undist, rhs_c_dist, 1);
 
-    res = op.GetDiag() * temp->GetVal(*iter_f);
+  if (rhs_c_dist.Global().BoundaryType() == GlobalCoarsened)
+    rhs_c_dist.ClearBoundary();
 
-    for (stencil_iter=op.begin(); stencil_iter!=op.end(); ++stencil_iter)
-      res += stencil_iter->Val() * temp->GetVal(iter_f->X() + stencil_iter->Disp().X(),
-                                                iter_f->Y() + stencil_iter->Disp().Y(),
-                                                iter_f->Z() + stencil_iter->Disp().Z());
-
-    rhs_c(*iter_c) = res;
-
-  }
-
-#ifdef DEBUG_MATRIX_CHECKS
-  rhs_c.IsConsistent();
-#endif
+  sol.ToCoarserLevel();
+  rhs.ToCoarserLevel();
 }
 
-void LevelOperatorCS::Prolongate(Grid& sol_f, Grid& rhs_f, Grid& sol_c, Grid& rhs_c)
+void LevelOperatorCS::Prolongate(Multigrid& sol, Multigrid& rhs)
 {
-#ifdef DEBUG_MATRIX_CHECKS
-  sol_c.IsConsistent();
-  rhs_c.IsConsistent();
-  sol_f.IsConsistent();
-  rhs_f.IsConsistent();
-#endif
-
-#ifdef DEBUG
-  sol_c.IsCompatible(rhs_c);
-  sol_f.IsCompatible(rhs_f);
-#endif
-
   Grid::iterator iter_f, iter_c;
   Stencil::iterator stencil_iter;
   vmg_float val;
 
+  Comm& comm = *MG::GetComm();
+
+  Grid& sol_c = sol(sol.Level());
+  Grid& rhs_c = rhs(rhs.Level());
+
+  Grid& sol_f_dist = sol(sol.Level()+1);
+  Grid& rhs_f_dist = rhs(rhs.Level()+1);
+
+  Grid& sol_f_undist = comm.GetFinerGrid(sol);
+  Grid& rhs_f_undist = comm.GetFinerGrid(rhs);
+
   const Stencil& op = OperatorProlongate();
-  Comm* comm = MG::GetComm();
 
-  Index begin_f = sol_f.Local().Begin() + 2*sol_c.Global().BeginLocal() - sol_f.Global().BeginLocal();
-  Index end_f = sol_f.Local().End();
+  Index begin_f = sol_f_undist.Local().Begin() + 2*sol_c.Global().LocalBegin() - sol_f_undist.Global().LocalBegin();
+  Index end_f = sol_f_undist.Local().End();
 
-  if (sol_c.Global().BoundaryType() == GlobalCoarsened)
-    for (int i=0; i<3; ++i) {
-
-      if (rhs_f.Local().HasBoundary1()[i])
-        ++begin_f[i];
-
-      if (rhs_f.Local().HasBoundary2()[i])
-        --end_f[i];
-
-    }
+  if (sol_c.Global().BoundaryType() == GlobalCoarsened) {
+    begin_f += rhs_f_undist.Local().BoundarySize1();
+    end_f -= rhs_f_undist.Local().BoundarySize2();
+  }
 
   const GridIteratorSet bounds_f(begin_f, end_f);
-  const GridIteratorSet bounds_c(sol_c.Local().AlignmentBegin(), sol_c.Local().AlignmentEnd());
+  const GridIteratorSet bounds_c(sol_c.Local().FinerBeginFoo(), sol_c.Local().FinerEndFoo());
 
-  Index apply_stencil_begin, apply_stencil_end;
+  sol_f_undist.ClearHalo();
 
-  for (int i=0; i<3; ++i) {
-    apply_stencil_begin[i] = (sol_f.Local().HasBoundary1()[i] ? 1 : 0);
-    apply_stencil_end[i] = (sol_f.Local().HasBoundary2()[i] ? sol_f.Local().End()[i] : sol_f.Local().SizeTotal()[i]);
-  }
-
-  sol_f.ClearHalo();
+  comm.CommSubgrid(sol_f_dist, sol_f_undist, 0);
 
   for (iter_f=bounds_f.Begin(), iter_c=bounds_c.Begin(); iter_c!=bounds_c.End(); iter_f+=2, ++iter_c) {
     val = sol_c.GetVal(*iter_c);
-    sol_f(*iter_f) += op.GetDiag() * val;
+    sol_f_undist(*iter_f) += op.GetDiag() * val;
     for (stencil_iter = op.begin(); stencil_iter != op.end(); ++stencil_iter)
-      sol_f(iter_f->X() + stencil_iter->Disp().X(),
-            iter_f->Y() + stencil_iter->Disp().Y(),
-            iter_f->Z() + stencil_iter->Disp().Z()) += stencil_iter->Val() * val;
+      sol_f_undist(iter_f->X() + stencil_iter->Disp().X(),
+                   iter_f->Y() + stencil_iter->Disp().Y(),
+                   iter_f->Z() + stencil_iter->Disp().Z()) += stencil_iter->Val() * val;
   }
 
-  comm->CommFromGhosts(sol_f);
+  comm.CommFromGhosts(sol_f_undist);
+  comm.CommSubgrid(sol_f_undist, sol_f_dist, 1);
+
+  sol.ToFinerLevel();
+  rhs.ToFinerLevel();
 }