Skip to content
Snippets Groups Projects
Commit 51092637 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

extended rosenbrock stationary - dirichle BC corrections

parent 3d8ac690
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
extensions/time/ExtendedRosenbrockStationary.cc
+ 44
− 24
View file @ 51092637
......@@ -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);
......
extensions/time/ExtendedRosenbrockStationary.h
+ 4
− 41
View file @ 51092637
......@@ -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;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment