Work on parallel jump residual estimator.

AMDiS/src/est/ResidualEstimator.cc

@@ -20,6 +20,7 @@
 
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
 #include <mpi.h>
+#include "parallel/MeshDistributor.h"
 #endif
 
 namespace AMDiS {
@@ -152,7 +153,128 @@ namespace AMDiS {
 	  secondOrderTerms[system] = secondOrderTerms[system] || (*it)->secondOrderTerms();
       }
     }
+
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+    initParallel();
+#endif
+  }
+
+
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+  void ResidualEstimator::initParallel()
+  {
+    FUNCNAME("ResidualEstimator::initParallel()");
+
+    if (C1 == 0.0)
+      return;
+
+    DOFVector<WorldVector<double> > coords(uh[0]->getFeSpace(), "tmp");
+    mesh->getDofIndexCoords(uh[0]->getFeSpace(), coords);
+
+    InteriorBoundary &intBoundary = 
+      MeshDistributor::globalMeshDistributor->getIntBoundary();
+    RankToBoundMap& otherBound = intBoundary.getOther();
+
+    ElInfo *elInfo = mesh->createNewElInfo();
+    elInfo->setFillFlag(Mesh::FILL_COORDS);
+
+    StdMpi<vector<double> > stdMpiDet(MeshDistributor::globalMeshDistributor->getMpiComm());
+    StdMpi<vector<vector<WorldVector<double> > > > stdMpiGrdUh(MeshDistributor::globalMeshDistributor->getMpiComm());
+
+    for (RankToBoundMap::iterator it = otherBound.begin();
+	 it != otherBound.end(); ++it) {
+      for (unsigned int i = 0; i < it->second.size(); i++) {
+	BoundaryObject &bObj = it->second[i].rankObj;
+	if (bObj.subObj == VERTEX)
+	  continue;
+
+	vector<BoundaryObject> subBound;
+	bObj.el->getSubBoundary(bObj, subBound);
+
+	WorldVector<int> faceIndices;
+
+	for (unsigned int j = 0; j < subBound.size(); j++) {
+	  Element *el = subBound[j].el;	  
+	  int oppV = subBound[j].ithObj;
+
+	  elInfo->setElement(el);
+	  el->sortFaceIndices(oppV, faceIndices);
+
+	  elInfo->getCoord(oppV) = coords[el->getDof(oppV, 0)];
+	  for (int k = 0; k < dim; k++)
+	    elInfo->getCoord(faceIndices[k]) = coords[el->getDof(k, 0)];
+
+	  double detNeigh = abs(elInfo->calcGrdLambda(*lambdaNeigh));
+	  stdMpiDet.getSendData(it->first).push_back(detNeigh);
+
+	  for (int system = 0; system < nSystems; system++) {
+	    if (matrix[system] == NULL || secondOrderTerms[system] == false) 
+	      continue;
+
+	    uh[system]->getLocalVector(el, uhNeigh[system]);
+
+	    for (int iq = 0; iq < nPointsSurface; iq++) {
+
+	      (*lambda)[oppV] = 0.0;
+	      for (int k = 0; k < dim; k++)
+		(*lambda)[faceIndices[k]] = surfaceQuad->getLambda(iq, k);
+
+	      basFcts[system]->evalGrdUh(*lambda, *lambdaNeigh, uhNeigh[system], grdUhNeigh[iq]);
+	    }
+
+	    stdMpiGrdUh.getSendData(it->first).push_back(grdUhNeigh);
+	  }
+	}
+      }
+    }
+
+    stdMpiDet.updateSendDataSize();
+    stdMpiGrdUh.updateSendDataSize();
+
+    RankToBoundMap& ownBound = intBoundary.getOwn();
+
+    for (RankToBoundMap::iterator it = ownBound.begin();
+	 it != ownBound.end(); ++it) {
+
+      for (unsigned int i = 0; i < it->second.size(); i++) {
+	BoundaryObject &bObj = it->second[i].rankObj;
+	if (bObj.subObj != VERTEX) {
+	  stdMpiDet.recv(it->first);
+	  stdMpiGrdUh.recv(it->first);
+	  break;
+	}
+      }
+    }
+
+    stdMpiDet.startCommunication();
+    stdMpiGrdUh.startCommunication();
+
+
+    for (RankToBoundMap::iterator it = ownBound.begin();
+	 it != ownBound.end(); ++it) {
+      vector<BoundaryObject> subBound;
+	
+      for (unsigned int i = 0; i < it->second.size(); i++) {
+	BoundaryObject &bObj = it->second[i].rankObj;
+	if (bObj.subObj == VERTEX)
+	  continue;
+
+	bObj.el->getSubBoundary(bObj, subBound);
+      }
+      
+      if (subBound.size() == 0)
+	continue;
+
+      TEST_EXIT_DBG(subBound.size() == stdMpiDet.getRecvData(it->first).size())
+	("Should not happen: %d %d from rank %d\n", subBound.size(), stdMpiDet.getRecvData(it->first).size(), it->first);
+      
+      TEST_EXIT_DBG(subBound.size() == stdMpiGrdUh.getRecvData(it->first).size())
+	("Should not happen!\n");      
+    }    
+
+    MSG("INIT PARALLEL FINISCHED!\n");
   }
+#endif
 
 
   void ResidualEstimator::exit(bool output)
@@ -356,7 +478,6 @@ namespace AMDiS {
     FUNCNAME("ResidualEstimator::computeJumpResidual()");
 
     double result = 0.0;
-    int dow = Global::getGeo(WORLD);
     Element *el = elInfo->getElement();
     const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
     double det = elInfo->getDet();
@@ -373,10 +494,8 @@ namespace AMDiS {
       el->sortFaceIndices(face, faceIndEl);
       neigh->sortFaceIndices(oppV, faceIndNeigh);	
       neighInfo->setElement(const_cast<Element*>(neigh));
-      neighInfo->setFillFlag(Mesh::FILL_COORDS);
-      
-      for (int i = 0; i < dow; i++)
-	neighInfo->getCoord(oppV)[i] = elInfo->getOppCoord(face)[i];
+      neighInfo->setFillFlag(Mesh::FILL_COORDS);     
+      neighInfo->getCoord(oppV) = elInfo->getOppCoord(face);
       
       // periodic leaf data ?
       ElementData *ldp = el->getElementData()->getElementData(PERIODIC);	
@@ -411,8 +530,7 @@ namespace AMDiS {
 	for (int i = 0; i < dim; i++) {
 	  int i1 = faceIndEl[i];
 	  int i2 = faceIndNeigh[i];
-	  for (int j = 0; j < dow; j++)
-	    neighInfo->getCoord(i2)[j] = elInfo->getCoord(i1)[j];
+	  neighInfo->getCoord(i2) = elInfo->getCoord(i1);
 	}
       }
 	
@@ -441,7 +559,7 @@ namespace AMDiS {
 	  
 	  (*lambda)[oppV] = 0.0;
 	  for (int i = 0; i < dim; i++)
-	    (*lambda)[faceIndNeigh[i]] = surfaceQuad->getLambda(iq, i);		  
+	    (*lambda)[faceIndNeigh[i]] = surfaceQuad->getLambda(iq, i);
 	  
 	  basFcts[system]->evalGrdUh(*lambda, *lambdaNeigh, uhNeigh[system], grdUhNeigh[iq]);
 	  

AMDiS/src/est/ResidualEstimator.h

@@ -85,13 +85,17 @@ namespace AMDiS {
     /// Constructor.
     ResidualEstimator(std::string name, int r);
 
-    virtual void init(double timestep);
+    void init(double timestep);
+
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+    void initParallel();
+#endif
 
     /// Estimates the error on an element. For more information about the
     /// parameter, see the description \ref Estimator::estimateElement.
-    virtual void estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo = NULL);
+    void estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo = NULL);
 
-    virtual void exit(bool output = true);
+    void exit(bool output = true);
 
   protected:
     /// Computes the element residual for a given element.

AMDiS/src/est/SimpleResidualEstimator.cc

@@ -230,7 +230,6 @@ namespace AMDiS {
 
     // === Init temporary variables. ===
     double result = 0.0;
-    int dow = Global::getGeo(WORLD);
     Element *el = elInfo->getElement();
     const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
     double det = elInfo->getDet();
@@ -251,16 +250,13 @@ namespace AMDiS {
       el->sortFaceIndices(face, faceIndEl);
       neigh->sortFaceIndices(oppV, faceIndNeigh);	
       neighInfo->setElement(const_cast<Element*>(neigh));
-      neighInfo->setFillFlag(Mesh::FILL_COORDS);
-      
-      for (int i = 0; i < dow; i++)
-	neighInfo->getCoord(oppV)[i] = elInfo->getOppCoord(face)[i];
+      neighInfo->setFillFlag(Mesh::FILL_COORDS);     
+      neighInfo->getCoord(oppV) = elInfo->getOppCoord(face);
       
       for (int i = 0; i < dim; i++) {
 	int i1 = faceIndEl[i];
 	int i2 = faceIndNeigh[i];
-	for (int j = 0; j < dow; j++)
-	  neighInfo->getCoord(i2)[j] = elInfo->getCoord(i1)[j];
+	neighInfo->getCoord(i2) = elInfo->getCoord(i1);
       }      
 	
       // Compute determinant of the neighbouring element.

AMDiS/src/parallel/StdMpi.cc

@@ -26,7 +26,7 @@ namespace AMDiS {
   MPI_Datatype StdMpiHelper<vector<std::pair<int, int> > >::mpiDataType = MPI_INT;
   MPI_Datatype StdMpiHelper<vector<WorldVector<double> > >::mpiDataType = MPI_DOUBLE;
   MPI_Datatype StdMpiHelper<map<WorldVector<double>, int> >::mpiDataType = MPI_DOUBLE;
-
+  MPI_Datatype StdMpiHelper<vector<vector<WorldVector<double> > > >::mpiDataType = MPI_DOUBLE;
 
   // T = int
 
@@ -449,4 +449,48 @@ namespace AMDiS {
     }
   }
 
+
+  // T = vector<vector<WorldVector<double> > >
+
+  int StdMpiHelper<vector<vector<WorldVector<double> > > >::getBufferSize(vector<vector<WorldVector<double> > > &data)
+  {
+    int result = 1;
+
+    for (unsigned int i = 0; i < data.size(); i++)
+      result += 1 + (data[i].size() * Global::getGeo(WORLD));
+
+    return result;
+  }
+
+  void StdMpiHelper<vector<vector<WorldVector<double> > > >::createBuffer(vector<vector<WorldVector<double> > > &data, 
+									  double* buf)
+  {
+    int counter = 0;
+
+    buf[counter++] = static_cast<double>(data.size());
+
+    for (unsigned int i = 0; i < data.size(); i++) {
+      buf[counter++] = data[i].size();
+      for (unsigned int j = 0; j < data[i].size(); j++)
+	for (unsigned int k = 0; k < Global::getGeo(WORLD); k++) 
+	  buf[counter++] = data[i][j][k];
+    }
+  }
+
+  void StdMpiHelper<vector<vector<WorldVector<double> > > >::makeFromBuffer(vector<vector<WorldVector<double> > > &data, 
+									    double* buf, int bufSize)
+  {
+    int counter = 0;
+    data.resize(static_cast<int>(buf[counter++]));
+
+    for (unsigned int i = 0; i < data.size(); i++) {
+      data[i].resize(static_cast<int>(buf[counter++]));
+      for (unsigned int j = 0; j < data[i].size(); j++) 
+	for (unsigned int k = 0; k < Global::getGeo(WORLD); k++) 
+	  data[i][j][k] = buf[counter++];
+    }
+
+    TEST_EXIT_DBG(counter == bufSize)("There is something very wrong!\n");
+  }
+
 }

AMDiS/src/parallel/StdMpi.h

@@ -229,6 +229,22 @@ namespace AMDiS {
 
 
 
+  template<>
+  struct StdMpiHelper<vector<vector<WorldVector<double> > > > {
+    static MPI_Datatype mpiDataType;
+
+    typedef double cppDataType;
+
+    typedef vector<vector<WorldVector<double> > >  DataType;
+
+    static int getBufferSize(DataType &data);
+
+    static void createBuffer(DataType &data, double* buf);
+
+    static void makeFromBuffer(DataType &data, double* buf, int bufSize);
+  };
+
+
   /** \brief
    * This class is used to easily send and receive STL containers using MPI.
    */
@@ -302,14 +318,14 @@ namespace AMDiS {
 
 
     /// Returns sending data, see \ref sendData.
-    map<int, SendT>& getSendData()
+    inline map<int, SendT>& getSendData()
     {
       return sendData;
     }
 
 
     /// Returns the data that should be send to a specific rank, see \ref sendData.
-    SendT& getSendData(int rank)
+    inline SendT& getSendData(int rank)
     {
       return sendData[rank];
     }