* Parallel file i/o

AMDiS/src/AdaptInstationary.cc

@@ -204,7 +204,21 @@ namespace AMDiS {
       iterationTimestamp_ = time(NULL);
 
       problemTime_->initTimestep(adaptInfo_);
+
+#ifdef _OPENMP
+#pragma omp parallel sections
+      {
+#pragma omp section
+	problemTime_->startDelayedTimestepCalculation();
+
+#pragma omp section
+	oneTimestep();
+      }
+#else
+      problemTime_->startDelayedTimestepCalculation();
       oneTimestep();
+#endif
+
       problemTime_->closeTimestep(adaptInfo_);
 
       if (breakWhenStable && (adaptInfo_->getSolverIterations() == 0)) {
@@ -220,6 +234,8 @@ namespace AMDiS {
       }
     }
 
+    problemTime_->startDelayedTimestepCalculation();
+
     return errorCode;
   }
 

AMDiS/src/DataCollector.cc

@@ -56,6 +56,19 @@ namespace AMDiS {
     DELETE dofCoords_;
   }
 
+  void DataCollector::fillAllData()
+  {
+    if (!elementDataCollected_) {
+      startCollectingElementData();
+    }
+    if (!periodicDataCollected_) {
+      startCollectingPeriodicData();
+    }
+    if (!valueDataCollected_) {
+      startCollectingValueData();
+    }
+  }
+
   int DataCollector::startCollectingElementData()
   {
     FUNCNAME("DataCollector::startCollectingElementData()");
@@ -427,7 +440,7 @@ namespace AMDiS {
 
   ::std::list<PeriodicInfo>* DataCollector::getPeriodicInfos()
   {
-    if(!periodicDataCollected_) {
+    if (!periodicDataCollected_) {
       startCollectingPeriodicData();
     }
 
@@ -463,7 +476,7 @@ namespace AMDiS {
   
   int DataCollector::getNumberConnections()
   {
-    if(!periodicDataCollected_) {
+    if (!periodicDataCollected_) {
       startCollectingPeriodicData();
     }
 

AMDiS/src/DataCollector.h

@@ -61,6 +61,11 @@ namespace AMDiS {
 
       ~DataCollector();
 
+      /** \brief
+       * Fills the DataCollector with all possible datas.
+       */
+      void fillAllData();
+
       /** \brief
        * Returns list with element information.
        */

AMDiS/src/DuneSolver.cc

@@ -71,6 +71,13 @@ namespace AMDiS {
   {
     FUNCNAME("DuneSolver::solveSystem()");
 
+    // Calculate residual reduction which should be achieved by the dune solver.
+    *p = *x;
+    *p *= -1.0;
+    mv->matVec(NoTranspose, *p, *r);
+    *r += *b;
+    double reduction = this->tolerance / norm(r);
+
     typedef Dune::Matrix<DuneMatrix> SystemDuneMatrix;
     typedef Dune::BlockVector<DuneVector> SystemDuneVector;
 
@@ -117,7 +124,7 @@ namespace AMDiS {
     Dune::InverseOperator<SystemDuneVector, SystemDuneVector> *duneSolver =
       getSolver<SystemDuneMatrix, SystemDuneVector>(&op, dunePreconditioner);
 
-    duneSolver->apply(duneVecX, duneVecB, opRes);
+    duneSolver->apply(duneVecX, duneVecB, reduction, opRes);
 
     for (int i = 0; i < nComponents; i++) {
       DuneVector::Iterator duneVecIt;    
@@ -129,6 +136,16 @@ namespace AMDiS {
       }
     }
 
+    // Calculate and print the residual.
+    *p = *x;
+    *p *= -1.0;
+    mv->matVec(NoTranspose, *p, *r);
+    *r += *b;
+    
+    this->residual = norm(r);
+
+    MSG("Residual: %e\n", this->residual);
+
     delete dunePreconditioner;
     delete duneSolver;
    

AMDiS/src/DuneSolver.h

@@ -123,12 +123,20 @@ namespace AMDiS {
     /** \brief
      * Implements OEMSolver<VectorType>::init().
      */
-    void init() {};
+    void init() 
+    {
+      p = this->vectorCreator->create();
+      r = this->vectorCreator->create();
+    };
     
     /** \brief
      * Implements OEMSolver<VectorType>::exit().
      */
-    void exit() {};
+    void exit() 
+    {
+      this->vectorCreator->free(p);
+      this->vectorCreator->free(r);
+    };
 
     /** \brief
      * Creates a dune preconditioner using the informations in the used init file. 
@@ -168,6 +176,12 @@ namespace AMDiS {
      * to be allocated before.
      */
     void mapDOFVector(DOFVector<double> *dofVector, DuneVector *duneVector);
+
+  private:
+    /** \brief
+     * These vectors are justed to calculate the final residual of the solution.
+     */
+    VectorType *r, *p;
   };
 }
 

AMDiS/src/DuneSolver.hh

@@ -90,7 +90,7 @@ namespace AMDiS {
 								 Dune::Preconditioner<V, V>* precon)
   {
     if (solverType_ == "bicgstab") {      
-      return new Dune::BiCGSTABSolver<V>(*ma, *precon, this->tolerance, this->max_iter, verbose_);
+      return new Dune::BiCGSTABSolver<V>(*ma, *precon, this->tolerance, this->max_iter, verbose_);      
     } else if (solverType_ == "cg") {
       return new Dune::CGSolver<V>(*ma, *precon, this->tolerance, this->max_iter, verbose_);
     } else if (solverType_ == "gradient") {

AMDiS/src/ElementFileWriter.h

@@ -30,6 +30,8 @@ namespace AMDiS {
 		    Flag traverseFlag = Mesh::CALL_LEAF_EL,
 		    bool (*writeElem)(ElInfo*) = NULL);
 
+    void writeDelayedFiles() {};
+
   protected:
     /**
      * Writes element data in tecplot format.

AMDiS/src/FileWriter.cc

@@ -5,7 +5,6 @@
 #include "ValueWriter.h"
 #include "MacroWriter.h"
 #include "VtkWriter.h"
-#include "DataCollector.h"
 #include "FiniteElemSpace.h"
 #include "AdaptInfo.h"
 #include "Flag.h"
@@ -28,6 +27,8 @@ namespace AMDiS {
 
     solutionVecs_.resize(1);
     solutionVecs_[0] = vec;
+
+    dataCollectors_.resize(1);
   }
 
 
@@ -48,6 +49,8 @@ namespace AMDiS {
 
     feSpace = vecs[0]->getFESpace();
     solutionVecs_ = vecs;
+
+    dataCollectors_.resize(vecs.size());
   }
 
 
@@ -79,6 +82,8 @@ namespace AMDiS {
     }
 
     feSpace = vec->getFESpace();
+
+    dataCollectors_.resize(nTmpSolutions_);
   }
 
 
@@ -91,6 +96,12 @@ namespace AMDiS {
         DELETE solutionVecs_[i]; 
       }
     }
+
+    for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
+      if (dataCollectors_[i]) {
+	DELETE dataCollectors_[i];
+      }
+    }
   }
   
 
@@ -111,6 +122,9 @@ namespace AMDiS {
     indexDecimals = 3;
     tsModulo = 1;
     nTmpSolutions_ = 0;
+    delayWriting_ = 0;
+    writingIsDelayed_ = false;
+    delayedFilename_ = "";
     paraViewAnimationFrames_.resize(0);
 
     readParameters();
@@ -135,6 +149,7 @@ namespace AMDiS {
     GET_PARAMETER(0, name + "->index length", "%d", &indexLength);
     GET_PARAMETER(0, name + "->index decimals", "%d", &indexDecimals);
     GET_PARAMETER(0, name + "->write every i-th timestep", "%d", &tsModulo);
+    GET_PARAMETER(0, name + "->delay", "%d", &delayWriting_);
   }
 
   void FileWriter::writeFiles(AdaptInfo *adaptInfo,
@@ -144,10 +159,26 @@ namespace AMDiS {
 			      bool (*writeElem)(ElInfo*))
   {
     FUNCNAME("FileWriter::writeFiles()");
-
+    
     if ((adaptInfo->getTimestepNumber() % tsModulo != 0) && !force) 
       return;
 
+    if (writingIsDelayed_) {
+      ERROR_EXIT("This should not happen!\n");
+    }
+
+    if (writeElem) {
+      for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
+	dataCollectors_[i] = NEW DataCollector(feSpace, solutionVecs_[i], 
+				  level, flag, writeElem);
+      }
+    } else {
+      for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
+	dataCollectors_[i] = NEW DataCollector(feSpace, solutionVecs_[i], 
+				  traverseLevel, flag | traverseFlag, writeElement);
+      }
+    }
+
     ::std::string fn = filename;
 
     if (appendIndex) {
@@ -156,31 +187,29 @@ namespace AMDiS {
       TEST_EXIT(indexDecimals < indexLength)("index length <= index decimals\n");
     
       char formatStr[9];
-
-      sprintf(formatStr, "%%0%d.%df", indexLength, indexDecimals);
-
       char timeStr[20];
 
+      sprintf(formatStr, "%%0%d.%df", indexLength, indexDecimals);
       sprintf(timeStr, formatStr, adaptInfo ? adaptInfo->getTime() : 0.0);
 
       fn += timeStr;
     }
 
-    ::std::vector< DataCollector* > dc(solutionVecs_.size());
 
-    if (writeElem) {
-      for (int i = 0; i < static_cast<int>(dc.size()); i++) {
-	dc[i] = NEW DataCollector(feSpace, solutionVecs_[i], 
-				  level, flag, writeElem);
+    if (delayWriting_) {
+      if (writeTecPlotFormat || writeAMDiSFormat || writePeriodicFormat) {
+	ERROR_EXIT("Delay writing only supported for ParaView file format!\n");
       }
-    } else {
-      for (int i = 0; i < static_cast<int>(dc.size()); i++) {
-	dc[i] = NEW DataCollector(feSpace, solutionVecs_[i], 
-				  traverseLevel, flag | traverseFlag, writeElement);
+
+      for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
+	dataCollectors_[i]->fillAllData();
       }
-    }
 
-    
+      writingIsDelayed_ = true;
+      delayedFilename_ = fn;
+      return;
+    }
+ 
     if (writeTecPlotFormat) {
       TecPlotWriter<DOFVector<double> >::writeValues(solutionVecs_[0], 
 						     const_cast<char*>((fn + tecplotExt).c_str()), 
@@ -191,42 +220,67 @@ namespace AMDiS {
     if (writeAMDiSFormat) {
       TEST_EXIT(mesh)("no mesh\n");
 
-      MacroWriter::writeMacro(dc[0], 
-			     const_cast<char*>( (fn +  amdisMeshExt).c_str()), 
+      MacroWriter::writeMacro(dataCollectors_[0], 
+			     const_cast<char*>((fn +  amdisMeshExt).c_str()), 
 			     adaptInfo ? adaptInfo->getTime() : 0.0);
       MSG("macro file written to %s\n", (fn + amdisMeshExt).c_str());
 
 
-      ValueWriter::writeValues(dc[0],
+      ValueWriter::writeValues(dataCollectors_[0],
 			       (fn + amdisDataExt).c_str(),
 			       adaptInfo ? adaptInfo->getTime() : 0.0);
       MSG("value file written to %s\n", (fn + amdisDataExt).c_str());   
     }
 
     if (writePeriodicFormat) {
-      MacroWriter::writePeriodicFile(dc[0], 
+      MacroWriter::writePeriodicFile(dataCollectors_[0], 
 				     (fn + periodicFileExt).c_str());
       MSG("periodic file written to %s\n", (fn + periodicFileExt).c_str());
     }
     
     if (writeParaViewFormat) {
-      VtkWriter vtkWriter(&dc);
-      
-      vtkWriter.writeFile(const_cast<char*>( (fn + paraViewFileExt).c_str()));
+      VtkWriter vtkWriter(&dataCollectors_);     
+      vtkWriter.writeFile(const_cast<char*>((fn + paraViewFileExt).c_str()));
+
       MSG("ParaView file written to %s\n", (fn + paraViewFileExt).c_str());
     }
 
     if (writeParaViewAnimation) {
-      VtkWriter vtkWriter(&dc);
-
+      VtkWriter vtkWriter(&dataCollectors_);
       vtkWriter.updateAnimationFile(fn + paraViewFileExt, 
 				    &paraViewAnimationFrames_, 
-				    const_cast<char*>( (filename + ".pvd").c_str()));
+				    const_cast<char*>((filename + ".pvd").c_str()));
     }
     
 
-    for (int i = 0; i < static_cast<int>(dc.size()); i++) {
-      DELETE dc[i];
+    for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
+      DELETE dataCollectors_[i];
+    }
+  }
+
+  void FileWriter::writeDelayedFiles()
+  {
+    if (!writingIsDelayed_) {
+      return;
+    }
+
+    if (writeParaViewFormat) {
+      VtkWriter vtkWriter(&dataCollectors_);     
+      vtkWriter.writeFile(const_cast<char*>((delayedFilename_ + paraViewFileExt).c_str()));
+    }
+
+    if (writeParaViewAnimation) {
+      VtkWriter vtkWriter(&dataCollectors_);
+      vtkWriter.updateAnimationFile(delayedFilename_ + paraViewFileExt, 
+				    &paraViewAnimationFrames_, 
+				    const_cast<char*>((filename + ".pvd").c_str()));
     }
+
+    for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
+      DELETE dataCollectors_[i];
+    }
+
+    writingIsDelayed_ = false;
+    delayedFilename_ = "";
   }
 }

AMDiS/src/FileWriter.h

@@ -31,6 +31,7 @@
 #include "MemoryManager.h"
 #include "MatrixVector.h"
 #include "Mesh.h"
+#include "DataCollector.h"
 
 namespace AMDiS {
 
@@ -69,6 +70,8 @@ namespace AMDiS {
 			    Flag traverseFlag = Mesh::CALL_LEAF_EL,
 			    bool (*writeElem)(ElInfo*) = NULL) = 0;
 
+    virtual void writeDelayedFiles() = 0;
+
     void setTraverseProperties(int level, 
 			       Flag flag,
 			       bool (*writeElem)(ElInfo*))
@@ -150,6 +153,8 @@ namespace AMDiS {
 			    Flag traverseFlag = Mesh::CALL_LEAF_EL,
 			    bool (*writeElem)(ElInfo*) = NULL);
 
+    virtual void writeDelayedFiles();
+
   protected:
     /** \brief
      * Initialization of the filewriter.
@@ -272,6 +277,28 @@ namespace AMDiS {
      * stored in solutionVecs_ must be deleted in the destructor.
      */
     int nTmpSolutions_;
+
+    /** \brief
+     * Containers, which store the data to be written;
+     */
+    ::std::vector< DataCollector* > dataCollectors_;
+
+    /** \brief
+     * If set to 1, the FileWriter will delay the file writing to the future, where
+     * it can be executed in parallel with some other independent calculations.
+     */
+    int delayWriting_;
+
+    /** \brief
+     * If set to true, the FileWriter was filled with data, but the Files are currently
+     * not written to disk.
+     */
+    bool writingIsDelayed_;
+
+    /** \brief
+     * Here the filename for the file, which should be written to the next, is stored.
+     */
+    ::std::string delayedFilename_;
   };
 
 }

AMDiS/src/ProblemInstat.cc

@@ -95,7 +95,7 @@ namespace AMDiS {
       oldSolution = NEW DOFVector<double>(problemStat->getFESpace(), name + "->uOld");
       oldSolution->refineInterpol(true);
       oldSolution->setCoarsenOperation(COARSE_INTERPOL);
-      if(problemStat->getEstimator()) {
+      if (problemStat->getEstimator()) {
 	dynamic_cast<Estimator*>(problemStat->getEstimator())
 	  ->addUhOldToSystem(0, oldSolution);
       }
@@ -184,4 +184,14 @@ namespace AMDiS {
     oldSolution->copy(*(problemStat->getSolution())); 
   }
 
+  void ProblemInstatScal::startDelayedTimestepCalculation()
+  {
+    problemStat->writeDelayedFiles();
+  }
+
+  void ProblemInstatVec::startDelayedTimestepCalculation()
+  {
+    problemStat->writeDelayedFiles();
+  }
+
 }

AMDiS/src/ProblemInstat.h

@@ -121,6 +121,11 @@ namespace AMDiS {
      */
     virtual void solveInitialProblem(AdaptInfo *adaptInfo);  
 
+    /** \brief
+     * Implements the parallel file writing for instationary problems.
+     */
+    virtual void startDelayedTimestepCalculation() {};
+
   protected:
     /** \brief
      * Name of the problem.
@@ -199,6 +204,7 @@ namespace AMDiS {
      */
     virtual void transferInitialSolution(AdaptInfo *adaptInfo);  
 
+    virtual void startDelayedTimestepCalculation();
 
     virtual void serialize(::std::ostream &out) {};
 
@@ -287,6 +293,8 @@ namespace AMDiS {
      */
     virtual void transferInitialSolution(AdaptInfo *adaptInfo);  
 
+    virtual void startDelayedTimestepCalculation();
+
     virtual void serialize(::std::ostream &out) {};
 
     virtual void deserialize(::std::istream &in) {};

AMDiS/src/ProblemScal.cc

@@ -30,9 +30,15 @@ namespace AMDiS {
   ProblemScal *ProblemScal::traversePtr_ = NULL;
 
   void ProblemScal::writeFiles(AdaptInfo *adaptInfo, bool force) {
-    ::std::list<FileWriterInterface*>::iterator it;
-    for(it = fileWriters_.begin(); it != fileWriters_.end(); ++it) {
-      (*it)->writeFiles(adaptInfo, force);
+    for (int i = 0; i < static_cast<int>(fileWriters_.size()); i++) {
+      fileWriters_[i]->writeFiles(adaptInfo, force);
+    }
+  }
+
+  void ProblemScal::writeDelayedFiles()
+  {
+    for (int i = 0; i < static_cast<int>(fileWriters_.size()); i++) {
+      fileWriters_[i]->writeDelayedFiles();
     }
   }
 
@@ -59,11 +65,11 @@ namespace AMDiS {
 				   AbstractFunction<double, WorldVector<double> >* b) 
   {
     DirichletBC *dirichlet = new DirichletBC(type, b, feSpace_);
-    if(systemMatrix_) 
+    if (systemMatrix_) 
       systemMatrix_->getBoundaryManager()->addBoundaryCondition(dirichlet);
-    if(rhs_) 
+    if (rhs_) 
       rhs_->getBoundaryManager()->addBoundaryCondition(dirichlet);
-    if(solution_)
+    if (solution_)
       solution_->getBoundaryManager()->addBoundaryCondition(dirichlet);
   }
 
@@ -71,11 +77,11 @@ namespace AMDiS {
 				   DOFVector<double> *vec) 
   {
     DirichletBC *dirichlet = new DirichletBC(type, vec);
-    if(systemMatrix_) 
+    if (systemMatrix_) 
       systemMatrix_->getBoundaryManager()->addBoundaryCondition(dirichlet);
-    if(rhs_) 
+    if (rhs_) 
       rhs_->getBoundaryManager()->addBoundaryCondition(dirichlet);
-    if(solution_)
+    if (solution_)
       solution_->getBoundaryManager()->addBoundaryCondition(dirichlet);
   }
 
@@ -84,9 +90,9 @@ namespace AMDiS {
 			       AbstractFunction<double, WorldVector<double> > *r)
   {
     RobinBC *robin = new RobinBC(type, n, r, feSpace_);
-    if(rhs_) 
+    if (rhs_) 
       rhs_->getBoundaryManager()->addBoundaryCondition(robin);
-    if(systemMatrix_)
+    if (systemMatrix_)
       systemMatrix_->getBoundaryManager()->addBoundaryCondition(robin);
   }
 
@@ -103,9 +109,9 @@ namespace AMDiS {
 			       DOFVector<double> *r)
   {
     RobinBC *robin = new RobinBC(type, n, r, feSpace_);
-    if(rhs_) 
+    if (rhs_) 
       rhs_->getBoundaryManager()->addBoundaryCondition(robin);
-    if(systemMatrix_)
+    if (systemMatrix_)
       systemMatrix_->getBoundaryManager()->addBoundaryCondition(robin);
   }
 
@@ -113,13 +119,9 @@ namespace AMDiS {
   {
     PeriodicBC *periodic = new PeriodicBC(type, feSpace_);
 
-    //   TEST_EXIT(mesh_->getPeriodicBCMap()[type] == NULL)
-    //     ("periodic condition already set in mesh\n");
-    //   mesh_->getPeriodicBCMap()[type] = periodic;
-
-    if(systemMatrix_) 
+    if (systemMatrix_) 
       systemMatrix_->getBoundaryManager()->addBoundaryCondition(periodic);
-    if(rhs_) 
+    if (rhs_) 
       rhs_->getBoundaryManager()->addBoundaryCondition(periodic);
   }
 
@@ -131,7 +133,6 @@ namespace AMDiS {
     ::std::string meshName("");
 
     GET_PARAMETER(0, name_ + "->info", "%d", &info_);
-
     GET_PARAMETER(0, name_ + "->mesh", &meshName);
 
     TEST_EXIT(meshName != "")
@@ -145,7 +146,7 @@ namespace AMDiS {
     // create the mesh
     mesh_ = NEW Mesh(meshName, dim);
 
-    switch(dim) {
+    switch (dim) {
     case 1:
       coarseningManager_ = NEW CoarseningManager1d();
       refinementManager_ = NEW RefinementManager1d();
@@ -174,10 +175,11 @@ namespace AMDiS {
 
   Flag ProblemScal::markElements(AdaptInfo *adaptInfo) 
   { 
-    if(marker_)
+    if (marker_)
       return marker_->markMesh(adaptInfo, mesh_);
     else
       WARNING("no marker\n");
+
     return 0;
   }
 
@@ -190,9 +192,9 @@ namespace AMDiS {
 
   Flag ProblemScal::coarsenMesh(AdaptInfo *adaptInfo) 
   {
-    if(adaptInfo->isCoarseningAllowed(0))
+    if (adaptInfo->isCoarseningAllowed(0)) {
       return coarseningManager_->coarsenMesh(mesh_);
-    else {
+    } else {
       WARNING("coarsening not allowed\n");
       return 0;
     }
@@ -201,7 +203,7 @@ namespace AMDiS {
   void ProblemScal::solve(AdaptInfo *adaptInfo) 
   {
     FUNCNAME("Problem::solve()");
-    if(!solver_) {
+    if (!solver_) {
       WARNING("no solver\n");
       return;
     }
@@ -240,20 +242,19 @@ namespace AMDiS {
     } else {
       if (initFlag.isSet(CREATE_MESH) || 
 	 ((!adoptFlag.isSet(INIT_MESH))&&
-	  (initFlag.isSet(INIT_SYSTEM)||initFlag.isSet(INIT_FE_SPACE)))) 
-	{
-	  createMesh();
-	} 
+	  (initFlag.isSet(INIT_SYSTEM)||initFlag.isSet(INIT_FE_SPACE)))) {
+	createMesh();
+      } 
+      
       if (adoptProblem && 
 	 (adoptFlag.isSet(INIT_MESH) || 
 	  adoptFlag.isSet(INIT_SYSTEM) ||
-	  adoptFlag.isSet(INIT_FE_SPACE)))
-	{
-