From 8e4f16acba7035befdf99b18fd902274ab544b4b Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Fri, 12 Aug 2011 08:12:40 +0000
Subject: [PATCH] CouplingIterationInterface created, similar to
CouplingTimeInterface
---
AMDiS/src/AMDiS.h | 4 +-
AMDiS/src/AMDiS_fwd.h | 1 +
AMDiS/src/CouplingIterationInterface.cc | 123 ++++++++++++++++++++++++
AMDiS/src/CouplingIterationInterface.h | 90 +++++++++++++++++
AMDiS/src/CouplingTimeInterface.h | 6 ++
AMDiS/src/CreatorMap.cc | 1 +
AMDiS/src/DOFVector.cc | 93 ++++++++++++++++++
AMDiS/src/DOFVector.h | 11 +++
AMDiS/src/ITL_OEMSolver.hh | 10 +-
AMDiS/src/MTL4Solver.h | 5 +
AMDiS/src/ProblemStat.cc | 9 ++
11 files changed, 350 insertions(+), 3 deletions(-)
create mode 100644 AMDiS/src/CouplingIterationInterface.cc
create mode 100644 AMDiS/src/CouplingIterationInterface.h
diff --git a/AMDiS/src/AMDiS.h b/AMDiS/src/AMDiS.h
index fae091ce..82b62be9 100644
--- a/AMDiS/src/AMDiS.h
+++ b/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"
diff --git a/AMDiS/src/AMDiS_fwd.h b/AMDiS/src/AMDiS_fwd.h
index a1070295..e3b96b80 100644
--- a/AMDiS/src/AMDiS_fwd.h
+++ b/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;
diff --git a/AMDiS/src/CouplingIterationInterface.cc b/AMDiS/src/CouplingIterationInterface.cc
new file mode 100644
index 00000000..1f609e47
--- /dev/null
+++ b/AMDiS/src/CouplingIterationInterface.cc
@@ -0,0 +1,123 @@
+// ============================================================================
+// == ==
+// == 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
diff --git a/AMDiS/src/CouplingIterationInterface.h b/AMDiS/src/CouplingIterationInterface.h
new file mode 100644
index 00000000..39c4e9df
--- /dev/null
+++ b/AMDiS/src/CouplingIterationInterface.h
@@ -0,0 +1,90 @@
+// ============================================================================
+// == ==
+// == 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
+
diff --git a/AMDiS/src/CouplingTimeInterface.h b/AMDiS/src/CouplingTimeInterface.h
index 4eeb3374..ec89be4e 100644
--- a/AMDiS/src/CouplingTimeInterface.h
+++ b/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_;
diff --git a/AMDiS/src/CreatorMap.cc b/AMDiS/src/CreatorMap.cc
index d21ba5c1..830a1321 100644
--- a/AMDiS/src/CreatorMap.cc
+++ b/AMDiS/src/CreatorMap.cc
@@ -149,5 +149,6 @@ namespace AMDiS {
void CreatorMap<NonLinSolver>::addDefaultCreators()
{
addCreator("newton", new Newton::Creator);
+ addCreator("newtonArmijo", new NewtonArmijo::Creator);
}
}
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index 001cf56c..0d390f63 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/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
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index 9488d148..e8e28eb3 100644
--- a/AMDiS/src/DOFVector.h
+++ b/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);
diff --git a/AMDiS/src/ITL_OEMSolver.hh b/AMDiS/src/ITL_OEMSolver.hh
index abd5dc72..a8514b0a 100644
--- a/AMDiS/src/ITL_OEMSolver.hh
+++ b/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:
diff --git a/AMDiS/src/MTL4Solver.h b/AMDiS/src/MTL4Solver.h
index 525306d4..c48b9808 100644
--- a/AMDiS/src/MTL4Solver.h
+++ b/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;
}
diff --git a/AMDiS/src/ProblemStat.cc b/AMDiS/src/ProblemStat.cc
index 3183c724..f437f85c 100644
--- a/AMDiS/src/ProblemStat.cc
+++ b/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]);
--
GitLab