From 802ccc5d889f0c6e86620574fc912ec39a828059 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Mon, 18 Aug 2008 14:17:50 +0000
Subject: [PATCH] * Operator VecGrad_FOT added * AdaptInfo::reset() does not
set timestep to 0.0, instead it reads the parameter from init file again *
ParaReal-stuff changed
---
AMDiS/src/AdaptInfo.h | 2 +
AMDiS/src/Operator.cc | 32 ++++++++++++++-
AMDiS/src/Operator.h | 61 +++++++++++++++++++++++++++++
AMDiS/src/parareal/AdaptParaReal.cc | 1 +
AMDiS/src/parareal/AdaptParaReal.h | 2 +-
AMDiS/src/parareal/ProblemBase.h | 12 +++++-
6 files changed, 105 insertions(+), 5 deletions(-)
diff --git a/AMDiS/src/AdaptInfo.h b/AMDiS/src/AdaptInfo.h
index d4be7b40..0bf7beca 100644
--- a/AMDiS/src/AdaptInfo.h
+++ b/AMDiS/src/AdaptInfo.h
@@ -213,6 +213,8 @@ namespace AMDiS {
timestepNumber = 0;
solverIterations = 0;
solverResidual = 0.0;
+
+ GET_PARAMETER(0, name + "->timestep", "%f", ×tep);
};
/** \brief
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index ab73a1f8..af480064 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -407,8 +407,6 @@ namespace AMDiS {
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
}
- // bis hierhin erstmal
-
void VecAndGradAtQP_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
@@ -2157,4 +2155,34 @@ namespace AMDiS {
}
}
+ void VecGrad_FOT::initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+ {
+ vecAtQPs = subAssembler->getVectorAtQPs(vec1, elInfo, quad);
+ gradAtQPs = subAssembler->getGradientsAtQPs(vec2, elInfo, quad);
+ }
+
+ void VecGrad_FOT::getLb(const ElInfo *elInfo, int numPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
+ const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ for (int iq = 0; iq < numPoints; iq++) {
+ lb(Lambda, (*vecFct)(vecAtQPs[iq], gradAtQPs[iq]), Lb[iq], 1.0);
+ }
+ }
+
+ void VecGrad_FOT::eval(int numPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double factor) const
+ {
+ if (grdUhAtQP) {
+ for (int iq = 0; iq < numPoints; iq++) {
+ WorldVector<double> b = (*vecFct)(vecAtQPs[iq], gradAtQPs[iq]);
+ result[iq] += b * grdUhAtQP[iq] * factor;
+ }
+ }
+ }
+
}
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 24c70d9c..c99b5a08 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -1900,6 +1900,67 @@ namespace AMDiS {
AbstractFunction<WorldVector<double>, WorldVector<double> > *g;
};
+
+ /**
+ * \ingroup Assembler
+ *
+ * \brief
+ * First order term: \f$ b(v(\vec{x}), \nabla w(\vec{x})) \cdot \nabla u(\vec{x})\f$.
+ */
+ class VecGrad_FOT : public FirstOrderTerm
+ {
+ public:
+ /** \brief
+ * Constructor.
+ */
+ VecGrad_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
+ BinaryAbstractFunction<WorldVector<double>, double,
+ WorldVector<double> > *vecFct_)
+ : FirstOrderTerm(vecFct_->getDegree()), vec1(dv1), vec2(dv2), vecFct(vecFct_)
+ {};
+
+ /** \brief
+ * Implementation of \ref OperatorTerm::initElement().
+ */
+ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
+ Quadrature *quad = NULL);
+
+ /** \brief
+ * Implements FirstOrderTerm::getLb().
+ */
+ void getLb(const ElInfo *elInfo, int qPoint,
+ VectorOfFixVecs<DimVec<double> >& Lb) const;
+
+ /** \brief
+ * Implements FirstOrderTerm::eval().
+ */
+ void eval(int numPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double factor) const;
+
+ protected:
+ /** \brief
+ * DOFVector to be evaluated at quadrature points.
+ */
+ DOFVectorBase<double>* vec1;
+ DOFVectorBase<double>* vec2;
+
+ /** \brief
+ * Vector v at quadrature points.
+ */
+ double *vecAtQPs;
+
+ WorldVector<double> *gradAtQPs;
+
+ /** \brief
+ * Function for b.
+ */
+ BinaryAbstractFunction<WorldVector<double>, double, WorldVector<double> > *vecFct;
+ };
+
// ============================================================================
class General_FOT : public FirstOrderTerm
diff --git a/AMDiS/src/parareal/AdaptParaReal.cc b/AMDiS/src/parareal/AdaptParaReal.cc
index f02e5754..9ddf1176 100644
--- a/AMDiS/src/parareal/AdaptParaReal.cc
+++ b/AMDiS/src/parareal/AdaptParaReal.cc
@@ -12,6 +12,7 @@ namespace AMDiS {
adaptInfo_->setTimestep(coarseTimestep);
pararealProb->setStoreSolutions(false, true);
pararealProb->setStoreInitSolution(true);
+ pararealProb->setInitSolutionVec(NULL);
AdaptInstationary::adapt();
// ParaReal iterations
diff --git a/AMDiS/src/parareal/AdaptParaReal.h b/AMDiS/src/parareal/AdaptParaReal.h
index c6a66416..943ddb5f 100644
--- a/AMDiS/src/parareal/AdaptParaReal.h
+++ b/AMDiS/src/parareal/AdaptParaReal.h
@@ -52,7 +52,7 @@ namespace AMDiS {
TEST_EXIT(fineTimestep > 0.0)("fineTimestep must be greater than zero!\n");
TEST_EXIT(coarseTimestep > 0.0)("coarseTimestep must be greater than zero!\n");
TEST_EXIT(fineTimestep < coarseTimestep)("fineTimestep must be smaller than coarseTimestep!\n");
- TEST_EXIT(pararealIter > 0)("ParaReal Iterations must be greater than zero!\n");
+ TEST_EXIT(pararealIter >= 0)("ParaReal Iterations must be >= 0!\n");
}
/** \brief
diff --git a/AMDiS/src/parareal/ProblemBase.h b/AMDiS/src/parareal/ProblemBase.h
index d129a5d0..a47b2707 100644
--- a/AMDiS/src/parareal/ProblemBase.h
+++ b/AMDiS/src/parareal/ProblemBase.h
@@ -165,7 +165,11 @@ namespace AMDiS {
}
virtual void solveInitialProblem(AdaptInfo *adaptInfo) {
- ProblemInstatScal::solveInitialProblem(adaptInfo);
+ if (initSolutionVec == NULL) {
+ ProblemInstatScal::solveInitialProblem(adaptInfo);
+ } else {
+ *(problemStat->getSolution()) = *initSolutionVec;
+ }
if (storeInitSolution) {
storeSolution(problemStat->getSolution());
@@ -194,7 +198,11 @@ namespace AMDiS {
}
virtual void solveInitialProblem(AdaptInfo *adaptInfo) {
- ProblemInstatVec::solveInitialProblem(adaptInfo);
+ if (initSolutionVec == NULL) {
+ ProblemInstatVec::solveInitialProblem(adaptInfo);
+ } else {
+ *(problemStat->getSolution()) = *initSolutionVec;
+ }
if (storeInitSolution) {
storeSolution(problemStat->getSolution());
--
GitLab