Set dirichlet boundaries to correct values.

[[Imported from SVN: r517]]

istl/paamg/test/anisotropic.hh

@@ -125,6 +125,19 @@ void fillValues(int N, M& mat, int overlapStart, int overlapEnd, int start, int
   }
 }
 
+template<int BS, class G, class L, int s>
+void setBoundary(Dune::BlockVector<Dune::FieldVector<double,BS> >& vec, const G& n, Dune::ParallelIndexSet<G,L,s>& indices)
+{
+  typedef typename Dune::ParallelIndexSet<G,L,s>::const_iterator Iter;
+  for(Iter i=indices.begin(); i != indices.end(); ++i) {
+    G x = i->global()/n;
+    G y = i->global()%n;
+
+    if(x==0 || y ==0 || x==n-1 || y==n-1 || i->local().attribute()==GridAttributes::copy)
+      vec[i->local()]=0;
+  }
+}
+
 template<int BS, class G, class L, int s>
 Dune::BCRSMatrix<Dune::FieldMatrix<double,BS,BS> > setupAnisotropic2d(int N, Dune::ParallelIndexSet<G,L,s>& indices, int *nout, double eps)
 {

istl/paamg/test/parallelamgtest.cc

@@ -9,10 +9,68 @@
 #include <dune/istl/schwarz.hh>
 #include <dune/istl/owneroverlapcopy.hh>
 #include <dune/common/mpicollectivecommunication.hh>
+#include <string>
+
+template<class T>
+class DoubleStepPreconditioner
+  : public Dune::Preconditioner<typename T::domain_type, typename T::range_type>
+{
+public:
+  typedef typename T::domain_type X;
+  typedef typename T::range_type Y;
+
+  enum {category = T::category};
+
+  DoubleStepPreconditioner(T& preconditioner_)
+    : preconditioner(&preconditioner_)
+  {}
+
+  virtual void pre (X& x, Y& b)
+  {
+    preconditioner->pre(x,b);
+  }
+
+  virtual void apply(X& v, const Y& d)
+  {
+    preconditioner->apply(v,d);
+    preconditioner->apply(v,d);
+  }
+
+  virtual void post (X& x)
+  {
+    preconditioner->post(x);
+  }
+private:
+  T* preconditioner;
+};
+
+
+class MPIError {
+public:
+  /** @brief Constructor. */
+  MPIError(std::string s, int e) : errorstring(s), errorcode(e){}
+  /** @brief The error string. */
+  std::string errorstring;
+  /** @brief The mpi error code. */
+  int errorcode;
+};
+
+void MPI_err_handler(MPI_Comm *comm, int *err_code, ...){
+  char *err_string=new char[MPI_MAX_ERROR_STRING];
+  int err_length;
+  MPI_Error_string(*err_code, err_string, &err_length);
+  std::string s(err_string, err_length);
+  std::cerr << "An MPI Error ocurred:"<<std::endl<<s<<std::endl;
+  delete[] err_string;
+  throw MPIError(s, *err_code);
+}
 
 int main(int argc, char** argv)
 {
   MPI_Init(&argc, &argv);
+  MPI_Errhandler handler;
+  MPI_Errhandler_create(MPI_err_handler, &handler);
+  MPI_Errhandler_set(MPI_COMM_WORLD, handler);
 
   const int BS=1;
   int N=250;
@@ -44,31 +102,43 @@ int main(int argc, char** argv)
 
   BCRSMat mat = setupAnisotropic2d<BS>(N, comm.indexSet(), &n, 1);
 
+  const BCRSMat& cmat = mat;
+
   comm.remoteIndices().rebuild<false>();
 
-  Vector b(mat.N()), x(mat.M());
+  Vector b(cmat.N()), x(cmat.M());
 
   b=0;
   x=100;
 
-  if(N<6)
-    Dune::printmatrix(std::cout, mat, "A", "row");
+  setBoundary(x, N, comm.indexSet());
+
+  Vector b1=b, x1=x;
 
+  if(N<6 && rank==0) {
+    Dune::printmatrix(std::cout, cmat, "A", "row");
+    Dune::printvector(std::cout, x, "x", "row");
+    Dune::printvector(std::cout, b, "b", "row");
+    Dune::printvector(std::cout, b1, "b1", "row");
+    Dune::printvector(std::cout, x1, "x1", "row");
+  }
 
   Dune::Timer watch;
 
   watch.reset();
-  Operator fop(mat, comm);
+  Operator fop(cmat, comm);
 
   typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<BCRSMat,Dune::Amg::FirstDiagonal> >
   Criterion;
-  typedef Dune::SeqSSOR<BCRSMat,Vector,Vector> Smoother;
-  typedef Dune::ParSSOR<BCRSMat,Vector,Vector,Communication> ParSmoother;
+  typedef Dune::SeqJac<BCRSMat,Vector,Vector> Smoother;
+  typedef Dune::BlockPreconditioner<Vector,Vector,Communication,Smoother> ParSmoother;
+  //typedef Dune::ParSSOR<BCRSMat,Vector,Vector,Communication> ParSmoother;
   typedef Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
+  //typedef Dune::Amg::BlockPreconditionerConstructionArgs<Smoother,Communication> SmootherArgs;
 
   SmootherArgs smootherArgs;
 
-  smootherArgs.iterations = 2;
+  smootherArgs.iterations = 1;
 
 
   Criterion criterion(15,coarsenTarget);
@@ -84,10 +154,17 @@ int main(int argc, char** argv)
   std::cout<<"Building hierarchy took "<<buildtime<<" seconds"<<std::endl;
 
   Dune::OverlappingSchwarzScalarProduct<Vector,Communication> sp(comm);
-  Dune::CGSolver<Vector> amgCG(fop, sp, amg, 10e-8, 80, 2);
+
+  Dune::CGSolver<Vector> amgCG(fop, sp, amg, 10e-8, 800, (rank==0) ? 2 : 0);
   watch.reset();
-  Dune::InverseOperatorResult r;
+  Dune::InverseOperatorResult r, r1;
+
   amgCG.apply(x,b,r);
+  MPI_Barrier(MPI_COMM_WORLD);
+  Smoother ssm(fop.getmat(),1,.8);
+  ParSmoother sm(ssm,comm);
+  DoubleStepPreconditioner<ParSmoother> dsp(sm);
+  Dune::CGSolver<Vector> cg(fop, sp, sm, 10e-08, 800, (rank==0) ? 2 : 0);
 
   double solvetime = watch.elapsed();
 
@@ -96,11 +173,11 @@ int main(int argc, char** argv)
   std::cout<<"AMG building took "<<(buildtime/r.elapsed*r.iterations)<<" iterations"<<std::endl;
   std::cout<<"AMG building together with slving took "<<buildtime+solvetime<<std::endl;
 
-  /*
-     watch.reset();
-     cg.apply(x,b,r);
+  watch.reset();
+
+  cg.apply(x1,b1,r1);
+
+  //std::cout<<"CG solving took "<<watch.elapsed()<<" seconds"<<std::endl;
 
-     std::cout<<"CG solving took "<<watch.elapsed()<<" seconds"<<std::endl;
-   */
   MPI_Finalize();
 }