diff --git a/extensions/time/ExtendedRosenbrockStationary.cc b/extensions/time/ExtendedRosenbrockStationary.cc
index 1fb320792197f1313cbe9289a97e3b81159990e8..b260cead555027c82206eb5beba9b54b8cf2b02e 100644
--- a/extensions/time/ExtendedRosenbrockStationary.cc
+++ b/extensions/time/ExtendedRosenbrockStationary.cc
@@ -38,7 +38,8 @@ using namespace AMDiS;
timeRhsVec = new SystemVector(*solution);
newUn = new SystemVector(*solution);
tmp = new SystemVector(*solution);
- lowSol = new SystemVector(*solution);
+ if (!fixedTimestep)
+ lowSol = new SystemVector(*solution);
stageSolution->set(0.0);
unVec->set(0.0);
@@ -64,13 +65,15 @@ using namespace AMDiS;
}
*newUn = *unVec;
- *lowSol = *unVec;
+ if (!fixedTimestep)
+ *lowSol = *unVec;
for (int i = 0; i < rm->getStages(); i++) {
stageTime = (*oldTime) + rm->getAlphaI(i) * (*tauPtr);
tauGammaI = rm->getGammaI(i) * (*tauPtr);
+ // stage-solution: u_s(i) = u_old + sum_{j=0}^{i-1} a_ij*U_j
*stageSolution = *unVec;
for (int j = 0; j < i; j++) {
*tmp = *(stageSolutions[j]);
@@ -78,14 +81,19 @@ using namespace AMDiS;
*stageSolution += *tmp;
}
+ // Dirichlet-BC implemented as additional algebraic equation u = g(x,t) on boundary
+ // => U_i = -u_s(i) + g(x,t_s(i)) + tau*gamma_i* d_t(g)(t_old) on boundary
+ // where u_s(i) = ith stage-solution, t_s(i) = ith stage-time
for (unsigned int j = 0; j < boundaries.size(); j++) {
boundaries[j].vec->interpol(boundaries[j].fct);
- *(boundaries[j].vec) -= *(stageSolution->getDOFVector(boundaries[j].row));
+ *(boundaries[j].vec) -= *(stageSolution->getDOFVector(boundaries[j].col));
+
if (boundaries[j].fctDt != NULL) {
- DOFVector<double>* tmpDt = new DOFVector<double>(getFeSpace(boundaries[j].row), "tmp");
+ // time derivative of dirichlet bc is given
+ DOFVector<double>* tmpDt = new DOFVector<double>(getFeSpace(boundaries[j].col), "tmp");
tmpDt->interpol(boundaries[j].fctDt);
*tmpDt *= tauGammaI;
- *(boundaries[j].vec) -= *tmpDt;
+ *(boundaries[j].vec) += *tmpDt; // alternativ transformDOF(...) benutzen, für x = x + a*y
delete tmpDt;
}
}
@@ -97,27 +105,39 @@ using namespace AMDiS;
*timeRhsVec += *tmp;
}
- super::buildAfterCoarsen(adaptInfo, flag | UPDATE_DIRICHLET_BC, (i == 0), asmVector);
- super::solve(adaptInfo, i == 0, i + 1 < rm->getStages());
+ super::buildAfterCoarsen(adaptInfo, flag | UPDATE_DIRICHLET_BC, (i == 0), true);
+
+// mtl::io::matrix_market_ostream out("matrix.mat");
+// SolverMatrix<Matrix<DOFMatrix*> > solverMatrix;
+// solverMatrix.setMatrix(*getSystemMatrix());
+// out << solverMatrix.getMatrix();
+// out.close();
+
+ super::solve(adaptInfo, i == 0, true);
*(stageSolutions[i]) = *solution;
*tmp = *solution;
*tmp *= rm->getM1(i);
-
*newUn += *tmp;
- *tmp = *solution;
- *tmp *= rm->getM2(i);
- *lowSol += *tmp;
+ if (!fixedTimestep) {
+ // lower order approximation, with coefficients of embedded method
+ // used for local error estimator
+ *tmp = *solution;
+ *tmp *= rm->getM2(i);
+ *lowSol += *tmp;
+ }
}
stageTime = (*oldTime) + (*tauPtr);
- for (int i = 0; i < nComponents; i++) {
- (*(lowSol->getDOFVector(i))) -= (*(newUn->getDOFVector(i)));
- adaptInfo->setTimeEstSum(lowSol->getDOFVector(i)->l2norm(), i+componentShift);
- }
+ if (!fixedTimestep) {
+ for (int i = 0; i < nComponents; i++) {
+ (*(lowSol->getDOFVector(i))) -= (*(newUn->getDOFVector(i)));
+ adaptInfo->setTimeEstSum(lowSol->getDOFVector(i)->l2norm(), i+componentShift);
+ }
+ }
}
@@ -156,11 +176,11 @@ using namespace AMDiS;
TEST_EXIT(invTauGamma)("This should not happen!\n");
Operator *op = new Operator(componentSpaces[row], componentSpaces[col]);
- op->addZeroOrderTerm(new Simple_ZOT(factor));
+ op->addTerm(new Simple_ZOT(factor));
super::addMatrixOperator(op, row, col, invTauGamma, invTauGamma);
Operator *opRhs = new Operator(componentSpaces[row]);
- opRhs->addZeroOrderTerm(new Phase_ZOT(timeRhsVec->getDOFVector(col), factor));
+ opRhs->addTerm(new Phase_ZOT(timeRhsVec->getDOFVector(col), factor));
super::addVectorOperator(opRhs, row);
}
@@ -171,7 +191,7 @@ using namespace AMDiS;
{
FUNCNAME("ExtendedRosenbrockStationary::addDirichletBC()");
- DOFVector<double>* vec = new DOFVector<double>(componentSpaces[row], "vec");
+ DOFVector<double>* vec = new DOFVector<double>(componentSpaces[col], "vec");
MyRosenbrockBoundary bound(fct, fctDt, vec, row, col);
boundaries.push_back(bound);
@@ -182,7 +202,7 @@ using namespace AMDiS;
void ExtendedRosenbrockStationary::addSingularDirichletBC(WorldVector<double> &pos, int row, int col,
AbstractFunction<double, WorldVector<double> > &fct)
{
- DOFVector<double>* vec = new DOFVector<double>(componentSpaces[row], "vec");
+ DOFVector<double>* vec = new DOFVector<double>(componentSpaces[col], "vec");
MyRosenbrockBoundary bound(&fct, NULL, vec, row, col);
boundaries.push_back(bound);
@@ -193,7 +213,7 @@ using namespace AMDiS;
void ExtendedRosenbrockStationary::addSingularDirichletBC(DegreeOfFreedom idx, int row, int col,
AbstractFunction<double, WorldVector<double> > &fct)
{
- DOFVector<double>* vec = new DOFVector<double>(componentSpaces[row], "vec");
+ DOFVector<double>* vec = new DOFVector<double>(componentSpaces[col], "vec");
MyRosenbrockBoundary bound(&fct, NULL, vec, row, col);
boundaries.push_back(bound);
@@ -204,7 +224,7 @@ using namespace AMDiS;
void ExtendedRosenbrockStationary::addImplicitDirichletBC(AbstractFunction<double, WorldVector<double> > &signedDist, int row, int col,
AbstractFunction<double, WorldVector<double> > &fct)
{
- DOFVector<double>* vec = new DOFVector<double>(componentSpaces[row], "vec");
+ DOFVector<double>* vec = new DOFVector<double>(componentSpaces[col], "vec");
MyRosenbrockBoundary bound(&fct, NULL, vec, row, col);
boundaries.push_back(bound);
@@ -215,7 +235,7 @@ using namespace AMDiS;
void ExtendedRosenbrockStationary::addImplicitDirichletBC(DOFVector<double> &signedDist, int row, int col,
AbstractFunction<double, WorldVector<double> > &fct)
{
- DOFVector<double>* vec = new DOFVector<double>(componentSpaces[row], "vec");
+ DOFVector<double>* vec = new DOFVector<double>(componentSpaces[col], "vec");
MyRosenbrockBoundary bound(&fct, NULL, vec, row, col);
boundaries.push_back(bound);
@@ -227,7 +247,7 @@ using namespace AMDiS;
int row, int col,
AbstractFunction<double, WorldVector<double> > &fct)
{
- DOFVector<double>* vec = new DOFVector<double>(componentSpaces[row], "vec");
+ DOFVector<double>* vec = new DOFVector<double>(componentSpaces[col], "vec");
MyRosenbrockBoundary bound(&fct, NULL, vec, row, col);
boundaries.push_back(bound);
@@ -239,7 +259,7 @@ using namespace AMDiS;
int row, int col,
AbstractFunction<double, WorldVector<double> > &fct)
{
- DOFVector<double>* vec = new DOFVector<double>(componentSpaces[row], "vec");
+ DOFVector<double>* vec = new DOFVector<double>(componentSpaces[col], "vec");
MyRosenbrockBoundary bound(&fct, NULL, vec, row, col);
boundaries.push_back(bound);
diff --git a/extensions/time/ExtendedRosenbrockStationary.h b/extensions/time/ExtendedRosenbrockStationary.h
index 50c1994d689965bf78eec0c38b89f4a44271ddc4..604482056a85150fb9409e41d00a8a731a82ec44 100644
--- a/extensions/time/ExtendedRosenbrockStationary.h
+++ b/extensions/time/ExtendedRosenbrockStationary.h
@@ -68,6 +68,7 @@ using namespace AMDiS;
: super(name),
componentShift(componentShift_),
first(true),
+ fixedTimestep(false),
minusOne(-1.0),
tauPtr(NULL),
tauGamma(NULL),
@@ -104,6 +105,8 @@ using namespace AMDiS;
rm = method;
init();
}
+
+ void setFixedTimestep(bool fixed) { fixedTimestep = fixed; }
void reset()
{
@@ -171,7 +174,6 @@ using namespace AMDiS;
return &stageTime;
}
-
double* getTauGammaI()
{
return &tauGammaI;
@@ -193,36 +195,6 @@ using namespace AMDiS;
ERROR_EXIT("Not yet supported!\n");
}
- /// Adds a Neumann boundary condition, where the rhs is given by an
- /// abstract function.
-// void addNeumannBC(BoundaryType type, int row, int col,
-// AbstractFunction<double, WorldVector<double> > *n)
-// {
-// FUNCNAME("ExtendedRosenbrockStationary::addNeumannBC()");
-//
-// ERROR_EXIT("Not yet supported!\n");
-// }
-
- /// Adds a Neumann boundary condition, where the rhs is given by an DOF
- /// vector.
-// void addNeumannBC(BoundaryType type, int row, int col,
-// DOFVector<double> *n)
-// {
-// FUNCNAME("ExtendedRosenbrockStationary::addNeumannBC()");
-//
-// ERROR_EXIT("Not yet supported!\n");
-// }
-
- /// Adds Robin boundary condition.
-// void addRobinBC(BoundaryType type, int row, int col,
-// AbstractFunction<double, WorldVector<double> > *n,
-// AbstractFunction<double, WorldVector<double> > *r)
-// {
-// FUNCNAME("ExtendedRosenbrockStationary::addRobinBC()");
-//
-// ERROR_EXIT("Not yet supported!\n");
-// }
-
/// special boundary conditions defined in ExtendedProblemStat
void addSingularDirichletBC(WorldVector<double> &pos, int row, int col,
@@ -245,16 +217,6 @@ using namespace AMDiS;
int row, int col,
AbstractFunction<double, WorldVector<double> > &fct);
-#if 0
- /// Adds a periodic boundary condition.
- void addPeriodicBC(BoundaryType type, int row, int col)
- {
- FUNCNAME("ExtendedRosenbrockStationary::addPeriodicBC()");
-
- ERROR_EXIT("Not yet supported!\n");
- }
-#endif
-
protected:
void init();
@@ -268,6 +230,7 @@ using namespace AMDiS;
int componentShift;
bool first;
+ bool fixedTimestep;
double minusOne;