CouplingIterationInterface created, similar to CouplingTimeInterface

AMDiS/src/AMDiS.h

@@ -39,6 +39,8 @@
 #include "CoarseningManager1d.h"
 #include "CoarseningManager2d.h"
 #include "CoarseningManager3d.h"
+#include "CouplingTimeInterface.h"
+#include "CouplingIterationInterface.h"
 #include "CreatorInterface.h"
 #include "CreatorMap.h"
 #include "Debug.h"
@@ -63,7 +65,7 @@
 #include "FixVec.h"
 #include "Flag.h"
 #include "Global.h"
-#include "Initfile.h" // replacement for Parameters.h
+#include "Initfile.h"
 #include "ITL_Preconditioner.h"
 #include "ITL_Solver.h"
 #include "Lagrange.h"

AMDiS/src/AMDiS_fwd.h

@@ -35,6 +35,7 @@ namespace AMDiS {
   class BoundaryManager;
   class CGSolver;
   class CoarseningManager;
+  class CouplingIterationInterface;
   class DiagonalPreconditioner; 
   class DOFAdmin; 
   class DOFContainer;

AMDiS/src/CouplingIterationInterface.cc

+// ============================================================================
+// ==                                                                        ==
+// == AMDiS - Adaptive multidimensional simulations                          ==
+// ==                                                                        ==
+// ==  http://www.amdis-fem.org                                              ==
+// ==                                                                        ==
+// ============================================================================
+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology 
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+
+/** \file CouplingIterationInterface.h */
+
+#include "Flag.h"
+#include "CouplingIterationInterface.h"
+
+namespace AMDiS {
+
+  /// add problem by number
+  void CouplingIterationInterface::addIterationInterface(ProblemIterationInterface* probIter, int number)
+  {
+    unsigned idx = std::min(problems.size(), static_cast<unsigned>(number)); 
+    if (number < 0)
+      idx = problems.size();
+    std::vector<ProblemIterationInterface*>::iterator pos = problems.begin() + idx;
+
+    problems.insert(pos, probIter);
+  };
+
+
+  /// Called before each adaption loop iteration.
+  void CouplingIterationInterface::beginIteration(AdaptInfo *adaptInfo)
+  { FUNCNAME("CouplingIterationInterface::beginIteration()");
+    MSG("\n");
+    MSG("begin of iteration number: %d/%d\n", 
+      adaptInfo->getTimestepNumber() + 1,
+      adaptInfo->getNumberOfTimesteps());
+    MSG("==================================================\n");
+  };
+
+
+  /** \brief
+    * Determines the execution order of the single adaption steps. If adapt is
+    * true, mesh adaption will be performed. This allows to avoid mesh adaption,
+    * e.g. in timestep adaption loops of timestep adaptive strategies.
+    */
+  Flag CouplingIterationInterface::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
+  { FUNCNAME("CouplingIterationInterface::oneIteration()");
+    Flag flag = 0;
+
+    for (unsigned i = 0; i < problems.size(); ++i) {
+      problems[i]->beginIteration(adaptInfo);
+      flag |= problems[i]->oneIteration(adaptInfo, toDo);
+      problems[i]->endIteration(adaptInfo);
+    }
+
+    return flag;
+  };
+
+
+  /// Called after each adaption loop iteration.
+  void CouplingIterationInterface::endIteration(AdaptInfo *adaptInfo)
+  { FUNCNAME("CouplingIterationInterface::endIteration()");
+    MSG("\n");
+    MSG("end of iteration number: %d\n", 
+      adaptInfo->getTimestepNumber() + 1);
+    MSG("==================================================\n");
+  };
+
+
+  /// Returns number of managed problems
+  int CouplingIterationInterface::getNumProblems()
+  {
+    int num= 0;
+    for (unsigned i = 0; i < problems.size(); ++i)
+      num += problems[i]->getNumProblems();
+    return num;
+  };
+
+
+  /** \brief
+    * Returns the problem with the given number. If only one problem
+    * is managed by this master problem, the number hasn't to be given.
+    */
+  ProblemStatBase *CouplingIterationInterface::getProblem(int number)
+  {
+    int maxNum = getNumProblems();
+    if (maxNum <= number || number <  0)
+      throw(std::runtime_error("Problem number out of range."));
+    
+    int sum = 0;
+    ProblemStatBase *probIter = NULL;
+    for (unsigned i = 0; i < problems.size(); ++i) {
+      if (sum + problems[i]->getNumProblems() <= number)
+        sum += problems[i]->getNumProblems();
+      else
+        probIter = problems[i]->getProblem(number - sum);
+    }
+    if (probIter == NULL)
+      throw(std::runtime_error("Problem not found. Should not happen, since number is in range."));
+    return probIter;
+  };
+
+
+  /// Returns the name of the problem.
+  std::string CouplingIterationInterface::getName(int number)
+  {
+    if (static_cast<int>(problems.size()) <= number || number <  0)
+      throw(std::runtime_error("Problem number out of range."));
+
+    return problems[number]->getName();
+  };
+
+} // namespace AMDiS

AMDiS/src/CouplingIterationInterface.h

+// ============================================================================
+// ==                                                                        ==
+// == AMDiS - Adaptive multidimensional simulations                          ==
+// ==                                                                        ==
+// ==  http://www.amdis-fem.org                                              ==
+// ==                                                                        ==
+// ============================================================================
+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology 
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+
+/** \file CouplingIterationInterface.h */
+
+#ifndef AMDIS_COUPLINGITERATIONINTERFACE_H
+#define AMDIS_COUPLINGITERATIONINTERFACE_H
+
+#include "Flag.h"
+#include "AdaptInfo.h"
+#include "ProblemIterationInterface.h"
+
+namespace AMDiS {
+
+  /** \brief
+   * Interface for master problems needed by the adaption loop. A master problem
+   * can handle one single or multiple coupled problems. In the latter case,
+   * the master problem can determine the execution order of the build, solve,
+   * estimate, and adapt steps of the single problems in \ref oneIteration().
+   * Standard execution order depends on the ordering of the problems in the
+   * problem-list, that is filled by addProblem. Alternatively one can a access
+   * each problem by an unique name.
+   */
+  class CouplingIterationInterface : public ProblemIterationInterface
+  {
+  public:
+    virtual ~CouplingIterationInterface() {}
+
+    /// add problem by number
+    virtual void addIterationInterface(ProblemIterationInterface* probIter, int number = -1);
+
+    /// Called before each adaption loop iteration.
+    virtual void beginIteration(AdaptInfo *adaptInfo);
+
+    /** \brief
+     * Determines the execution order of the single adaption steps. If adapt is
+     * true, mesh adaption will be performed. This allows to avoid mesh adaption,
+     * e.g. in timestep adaption loops of timestep adaptive strategies.
+     */
+    virtual Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION);
+
+    /// Called after each adaption loop iteration.
+    virtual void endIteration(AdaptInfo *adaptInfo);
+
+    /// Returns number of managed problems
+    virtual int getNumProblems();
+
+    /** \brief
+     * Returns the problem with the given number. If only one problem
+     * is managed by this master problem, the number hasn't to be given.
+     */
+    virtual ProblemStatBase *getProblem(int number = 0);
+
+    /// Returns the name of the problem with the given number.
+    virtual std::string getName(int number);
+    virtual std::string getName() { return getName(0); }
+
+    /// Function that serializes the problem plus information about the iteration.
+    virtual void serialize(std::ostream &out) {};
+
+    /// Function that deserializes the problem plus information about the iteration.
+    virtual void deserialize(std::istream &in) {};
+
+  protected:
+
+    /// vector/map of coupled stationary problems
+    std::vector<ProblemIterationInterface*> problems;
+  };
+
+}
+
+#endif // AMDIS_COUPLINGITERATIONINTERFACE_H
+

AMDiS/src/CouplingTimeInterface.h

@@ -82,6 +82,12 @@ namespace AMDiS {
 	interfaces_[i]->transferInitialSolution(adaptInfo);
     }
 
+    /// Function that serializes the problem plus information about the iteration.
+    virtual void serialize(std::ostream &out) {};
+
+    /// Function that deserializes the problem plus information about the iteration.
+    virtual void deserialize(std::istream &in) {};
+
   protected:
     /// vector of coupled time interfaces
     std::vector<ProblemTimeInterface*> interfaces_;

AMDiS/src/CreatorMap.cc

@@ -149,5 +149,6 @@ namespace AMDiS {
   void CreatorMap<NonLinSolver>::addDefaultCreators()
   {
     addCreator("newton", new Newton::Creator);
+    addCreator("newtonArmijo", new NewtonArmijo::Creator);
   }
 }

AMDiS/src/DOFVector.cc

@@ -67,6 +67,99 @@ namespace AMDiS {
   }
 
 
+  template<>
+  const double& DOFVector<double>::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, double* values) const
+  {
+    FUNCNAME("DOFVector<double>::evalAtCoords()");
+  
+    Mesh *mesh = feSpace->getMesh();
+    const BasisFunction *basFcts = feSpace->getBasisFcts();
+  
+    int dim = mesh->getDim();
+    int nBasFcts = basFcts->getNumber();
+  
+    DegreeOfFreedom *localIndices = new DegreeOfFreedom[nBasFcts];
+    DimVec<double> lambda(dim, NO_INIT);
+  
+    ElInfo *elInfo = mesh->createNewElInfo();
+    double value = 0.0;
+    bool inside = false;
+  
+    if (oldElInfo && oldElInfo->getMacroElement()) {
+      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
+      delete oldElInfo;
+    } else
+      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);
+
+    if (oldElInfo)
+      oldElInfo = elInfo;
+
+    if (inside) {
+      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
+      ElementVector uh(lambda.size());
+      for (int i = 0; i < lambda.size(); i++)
+        uh[i] = operator[](localIndices[i]);
+      value = basFcts->evalUh(lambda, uh);
+    } else
+      throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));
+
+    delete [] localIndices;
+    if (oldElInfo == NULL)
+      delete elInfo;
+
+    if (values != NULL)
+      *values = value;
+    return value;
+  };
+
+
+  template<>
+  const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* values) const
+  {
+    FUNCNAME("DOFVector<double>::evalAtCoords()");
+  
+    Mesh *mesh = feSpace->getMesh();
+    const BasisFunction *basFcts = feSpace->getBasisFcts();
+  
+    int dim = mesh->getDim();
+    int nBasFcts = basFcts->getNumber();
+  
+    DegreeOfFreedom *localIndices = new DegreeOfFreedom[nBasFcts];
+    DimVec<double> lambda(dim, NO_INIT);
+  
+    ElInfo *elInfo = mesh->createNewElInfo();
+
+    static WorldVector<double> Values(DEFAULT_VALUE, 0.0);
+    WorldVector<double> *val = (NULL != values) ? values : &Values;
+
+    bool inside = false;
+  
+    if (oldElInfo && oldElInfo->getMacroElement()) {
+      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
+      delete oldElInfo;
+    } else
+      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);
+
+    if (oldElInfo)
+      oldElInfo = elInfo;
+
+    if (inside) {
+      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
+      mtl::dense_vector<WorldVector<double> > uh(lambda.size());
+      for (int i = 0; i < lambda.size(); i++)
+        uh[i] = operator[](localIndices[i]);
+      basFcts->evalUh(lambda, uh, val);
+    } else
+      throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));
+
+    delete [] localIndices;
+    if (oldElInfo == NULL)
+      delete elInfo;
+
+    return ((*val));
+  };
+
+
   template<>
   DOFVector<WorldVector<double> >*
   DOFVector<double>::getGradient(DOFVector<WorldVector<double> > *grad) const 

AMDiS/src/DOFVector.h

@@ -539,6 +539,11 @@ namespace AMDiS {
 
     void interpol(DOFVector<T> *v, double factor = 1.0);
 
+
+    /// eval DOFVector at given point p. If oldElInfo != NULL the search for the element, where p is inside,
+    /// starts from oldElInfo.
+    const T& evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo = NULL, T* value = NULL) const;
+
     /// Writes the data of the DOFVector to an output stream.
     void serialize(std::ostream &out) 
     {
@@ -577,6 +582,12 @@ namespace AMDiS {
   }; 
 
 
+  template<>
+  const double& DOFVector<double>::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, double* value) const;
+
+  template<>
+  const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* value) const;
+
   template<>
   void DOFVector<double>::refineInterpol(RCNeighbourList&, int); 
 

AMDiS/src/ITL_OEMSolver.hh

@@ -69,13 +69,16 @@ namespace AMDiS {
       FUNCNAME("ITL_OEMSolver_para_runner::solve()");
       TEST_EXIT(preconPair.l != NULL)("there is no left preconditioner");
       TEST_EXIT(preconPair.r != NULL)("there is no right preconditioner");
-      itl::cyclic_iteration<typename Vector::value_type> iter(b, oem.getMaxIterations(), oem.getRelative(), 
+
+      mtl::dense_vector<typename Matrix::value_type> r(A * x - b); 
+      itl::cyclic_iteration<typename Vector::value_type> iter(r, oem.getMaxIterations(), oem.getRelative(), 
 							      oem.getTolerance(), oem.getPrint_cycle());
 
       int error = solver(A, x, b, *(preconPair.l), *(preconPair.r), iter, ell);
       oem.setIterations(iter.iterations());
       oem.setResidual(iter.resid());
       oem.setErrorCode(error);
+
       return error;
     }
   private:
@@ -101,13 +104,16 @@ namespace AMDiS {
       FUNCNAME("ITL_OEMSolver_runner::solve()");
       TEST_EXIT(preconPair.l != NULL)("there is no left preconditioner");
       TEST_EXIT(preconPair.r != NULL)("there is no right preconditioner");
-      itl::cyclic_iteration<typename Matrix::value_type> iter(b, oem.getMaxIterations(), oem.getRelative(), 
+
+      mtl::dense_vector<typename Matrix::value_type> r(A * x - b); 
+      itl::cyclic_iteration<typename Matrix::value_type> iter(r, oem.getMaxIterations(), oem.getRelative(), 
 							      oem.getTolerance(), oem.getPrint_cycle());
 
       int error = solver(A, x, b, *(preconPair.l), *(preconPair.r), iter);
       oem.setErrorCode(error);
       oem.setIterations(iter.iterations());
       oem.setResidual(iter.resid());
+
       return error;
     }
     private:

AMDiS/src/MTL4Solver.h

@@ -55,6 +55,11 @@ namespace AMDiS {
 
       error = worker.solve(matrix ,mtl_x, mtl_b);
 
+      mtl::dense_vector<typename MTLMatrix::value_type> r(mtl_b); 
+      r -= matrix * mtl_x; 
+      double residual = two_norm(r);
+      MSG("MTL4Solver: ||b-Ax||= %e\n", residual);
+
       mtl_x >> xVecMap;
       return error;
     }

AMDiS/src/ProblemStat.cc

@@ -777,9 +777,12 @@ namespace AMDiS {
     }  
 
     if (asmMatrix) {
+      solver->setMultipleRhs(false);
       solverMatrix.setMatrix(*systemMatrix);
 
       INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
+    } else {
+      solver->setMultipleRhs(true);
     }
 
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
@@ -973,6 +976,7 @@ namespace AMDiS {
 		  TIME_USED(first, clock()));
   }
 
+
   bool ProblemStatSeq::dualMeshTraverseRequired()
   {
     FUNCNAME("ProblemStat::dualMeshTraverseRequired()");
@@ -1492,6 +1496,7 @@ namespace AMDiS {
       rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(robin);
   }
 
+
   void ProblemStatSeq::addPeriodicBC(BoundaryType type, int row, int col) 
   {
     FUNCNAME("ProblemStat::addPeriodicBC()");
@@ -1514,6 +1519,8 @@ namespace AMDiS {
   {
     FUNCNAME("ProblemStat::addBoundaryOperator()");
 
+    boundaryConditionSet = true;
+
     RobinBC *robin = 
       new RobinBC(type, NULL, op, componentSpaces[row], componentSpaces[col]);
 
@@ -1527,6 +1534,8 @@ namespace AMDiS {
   {
     FUNCNAME("ProblemStat::addBoundaryOperator()");
 
+    boundaryConditionSet = true;
+
     RobinBC *robin = 
       new RobinBC(type, op, NULL, componentSpaces[row]);