Changeset f57182 for src/level/level_operator_fas.cpp

Timestamp:

Mar 30, 2013, 2:44:52 PM (13 years ago)

Author:

Julian Iseringhausen <isering@…>

Children:

8180d8

Parents:

d13e27

git-author:

Julian Iseringhausen <isering@…> (06/11/12 14:02:16)

git-committer:

Julian Iseringhausen <isering@…> (03/30/13 14:44:52)

Message:

Open boundary conditions.

Conflicts:

lib/vmg/src/Makefile.am
lib/vmg/src/base/factory.cpp
lib/vmg/test/unit_test/library/dirichlet_fas_lr_mpi.cpp

File:

• src/level/level_operator_fas.cpp (modified) (6 diffs)

src/level/level_operator_fas.cpp

@@ -30 +30 @@
 #endif
 
+#include <limits>
+
 #include "base/discretization.hpp"
+#include "base/helper.hpp"
 #include "base/index.hpp"
 #include "comm/comm.hpp"
@@ -56 +59 @@
   Grid& rhs_c_undist = comm.GetCoarserGrid(rhs);
 
+  Index begin_c = rhs_c_undist.Local().FinerBegin();
+  Index end_c = rhs_c_undist.Local().FinerEnd();
+
   Index begin_f = rhs_f.Local().Begin();
   Index end_f = rhs_f.Local().End();
 
   if (rhs_c_undist.Global().BoundaryType() == GlobalCoarsened) {
+    begin_c += rhs_c_undist.Local().BoundarySize1();
+    end_c -= rhs_c_undist.Local().BoundarySize2();
     begin_f += rhs_f.Local().BoundarySize1();
     end_f -= rhs_f.Local().BoundarySize2();
   }
 
+  Helper::AssertVectorsEqual(rhs_c_undist.GetSpatialPos(begin_c), rhs_f.GetSpatialPos(begin_f));
+
   const GridIteratorSet bounds_f(begin_f, end_f);
-  const GridIteratorSet bounds_c(rhs_c_undist.Local().FinerBegin(), rhs_c_undist.Local().FinerEnd());
+  const GridIteratorSet bounds_c(begin_c, end_c);
+
+  assert(bounds_c.Begin().GetEnd()-bounds_c.Begin().GetBegin() == ((bounds_f.Begin().GetEnd()-bounds_f.Begin().GetBegin())-1)/2+1);
 
   const Stencil& op_res = OperatorRestrict();
@@ -75 +87 @@
   comm.CommToGhosts(sol_f);
   comm.CommSubgrid(sol_c_dist, sol_c_undist, 0);
+  comm.CommSubgrid(rhs_c_dist, rhs_c_undist, 0);
 
   MG::GetDiscretization()->SetInnerBoundary(sol_f, rhs_f, sol_c_undist);
@@ -83 +96 @@
   res_grid->ImportFromResidual(sol_f, rhs_f);
 
-  for (iter_f=bounds_f.Begin(), iter_c=bounds_c.Begin(); iter_c!=bounds_c.End(); iter_f+=2, ++iter_c)
+  for (iter_f=bounds_f.Begin(), iter_c=bounds_c.Begin(); iter_c!=bounds_c.End(); iter_f+=2, ++iter_c) {
+    Helper::AssertVectorsEqual(sol_c_undist.GetSpatialPos(*iter_c), sol_f.GetSpatialPos(*iter_f));
     sol_c_undist(*iter_c) = op_sol.Apply(sol_f, *iter_f);
+  }
 
-  for (iter_f=bounds_f.Begin(), iter_c=bounds_c.Begin(); iter_c!=bounds_c.End(); iter_f+=2, ++iter_c)
-    rhs_c_undist(*iter_c) = op_res.Apply(*res_grid, *iter_f);
+  comm.CommToGhosts(sol_c_undist);
+
+  for (iter_f=bounds_f.Begin(), iter_c=bounds_c.Begin(); iter_c!=bounds_c.End(); iter_f+=2, ++iter_c) {
+    Helper::AssertVectorsEqual(rhs_c_undist.GetSpatialPos(*iter_c), (*res_grid).GetSpatialPos(*iter_f));
+    rhs_c_undist(*iter_c) = op_res.Apply(*res_grid, *iter_f) + prefactor * op_pde.Apply(sol_c_undist, *iter_c);
+  }
+
+  /*
+  if (rhs_c_undist.Global().BoundaryType() == GlobalCoarsened)
+    rhs_c_undist.SetBoundary(sol_c_undist);
+  */
 
   comm.CommSubgrid(sol_c_undist, sol_c_dist, 1);
-  comm.CommToGhosts(sol_c_dist);
-
   comm.CommSubgrid(rhs_c_undist, rhs_c_dist, 1);
 
@@ -98 +120 @@
   sol_old->SetGrid(sol_c_dist);
 
-  for (iter_c=rhs_c_dist.Iterators().Local().Begin(); iter_c!=rhs_c_dist.Iterators().Local().End(); ++iter_c)
-    rhs_c_dist(*iter_c) += prefactor * op_pde.Apply(sol_c_dist, *iter_c);
-
-  comm.CommToGhosts(rhs_c_dist);
+  //  comm.CommToGhosts(rhs_c_dist);
 
   sol.ToCoarserLevel();
@@ -125 +144 @@
   TempGrid *sol_old = this->GetTempGrid(sol_c.Level());
 
+  Index begin_c = sol_c.Local().FinerBegin();
+  Index end_c = sol_c.Local().FinerEnd();
+
   Index begin_f = sol_f_undist.Local().Begin();
   Index end_f = sol_f_undist.Local().End();
 
   if (sol_c.Global().BoundaryType() == GlobalCoarsened) {
-    begin_f += rhs_f_undist.Local().BoundarySize1();
-    end_f -= rhs_f_undist.Local().BoundarySize2();
+    begin_c += sol_c.Local().BoundarySize1();
+    end_c -= sol_c.Local().BoundarySize2();
+    begin_f += sol_f_undist.Local().BoundarySize1();
+    end_f -= sol_f_undist.Local().BoundarySize2();
   }
 
+  Helper::AssertVectorsEqual(sol_c.GetSpatialPos(begin_c), sol_f_undist.GetSpatialPos(begin_f));
+
   const GridIteratorSet bounds_f(begin_f, end_f);
-  const GridIteratorSet bounds_c(sol_c.Local().FinerBegin(), sol_c.Local().FinerEnd());
+  const GridIteratorSet bounds_c(begin_c, end_c);
 
+  assert(bounds_c.Begin().GetEnd()-bounds_c.Begin().GetBegin() == ((bounds_f.Begin().GetEnd()-bounds_f.Begin().GetBegin())-1)/2+1);
+
+
+  comm.CommSubgrid(sol_f_dist, sol_f_undist, 0);
   sol_f_undist.ClearHalo();
 
   for (iter_f=bounds_f.Begin(), iter_c=bounds_c.Begin(); iter_f!=bounds_f.End(); iter_f+=2, ++iter_c) {
+    Helper::AssertVectorsEqual(sol_c.GetSpatialPos(*iter_c), sol_f_undist.GetSpatialPos(*iter_f));
     val = sol_c.GetVal(*iter_c) - sol_old->GetVal(*iter_c);
     sol_f_undist(*iter_f) += op.GetDiag() * val;