Commit e3fbdd85 authored by Praetorius, Simon's avatar Praetorius, Simon

fixed issue with RosenBrockStationary and changed ProblemStat::buildAfterCoarsen to name assemble

parent a877b6da
...@@ -835,8 +835,8 @@ namespace AMDiS { ...@@ -835,8 +835,8 @@ namespace AMDiS {
} }
void ProblemStatSeq::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag, void ProblemStatSeq::assemble(AdaptInfo *adaptInfo, Flag flag,
bool asmMatrix, bool asmVector) bool asmMatrix, bool asmVector)
{ {
FUNCNAME("ProblemStat::buildAfterCoarsen()"); FUNCNAME("ProblemStat::buildAfterCoarsen()");
......
...@@ -139,14 +139,22 @@ namespace AMDiS { ...@@ -139,14 +139,22 @@ namespace AMDiS {
/// Implementation of ProblemStatBase::buildAfterCoarsen(). /// Implementation of ProblemStatBase::buildAfterCoarsen().
/// Assembles \ref A and \ref rhs. With the last two parameters, assembling /// Assembles \ref A and \ref rhs. With the last two parameters, assembling
/// can be restricted to matrices or vectors only. /// can be restricted to matrices or vectors only.
virtual void assemble(AdaptInfo *adaptInfo, Flag flag,
bool assembleMatrix = true,
bool assembleVector = true);
/// Alias for \ref assemble
void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag, void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
bool assembleMatrix = true, bool assembleMatrix = true,
bool assembleVector = true) override; bool assembleVector = true) override
{
assemble(adaptInfo, flag, assembleMatrix, assembleVector);
}
/// assemble all operators of matrix and vector side /// assemble all operators of matrix and vector side
void assemble(AdaptInfo* adaptInfo) void assemble(AdaptInfo* adaptInfo)
{ {
buildAfterCoarsen(adaptInfo, 0, true, true); assemble(adaptInfo, 0, true, true);
} }
bool dualMeshTraverseRequired(); bool dualMeshTraverseRequired();
......
...@@ -70,8 +70,12 @@ namespace AMDiS { ...@@ -70,8 +70,12 @@ namespace AMDiS {
if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_COARSENED)) if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_COARSENED))
flag |= problem->coarsenMesh(adaptInfo); flag |= problem->coarsenMesh(adaptInfo);
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Parallel::MeshDistributor::globalMeshDistributor->checkMeshChange();
#endif
if (toDo.isSet(BUILD) || toDo.isSet(SOLVE)) { if (toDo.isSet(BUILD) || toDo.isSet(SOLVE)) {
flag = stageIteration(adaptInfo, toDo, true, true); flag = stageIteration(adaptInfo);
estimateTimeError(adaptInfo); estimateTimeError(adaptInfo);
} }
...@@ -82,17 +86,22 @@ namespace AMDiS { ...@@ -82,17 +86,22 @@ namespace AMDiS {
} }
Flag RosenbrockStationary::stageIteration(AdaptInfo *adaptInfo, Flag flag, Flag RosenbrockStationary::stageIteration(AdaptInfo *adaptInfo)
bool asmMatrix, bool asmVector)
{ {
FUNCNAME("RosenbrockStationary::stageIteration()"); FUNCNAME("RosenbrockStationary::stageIteration()");
TEST_EXIT(tauPtr)("No tau pointer defined in stationary problem!\n"); TEST_EXIT(tauPtr)("No tau pointer defined in stationary problem!\n");
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
if (first) { if (first) {
first = false; first = false;
*unVec = *solution; *unVec = *solution;
} }
#else
// NOTE: To guarantee that unVec is synchronized with solution
// after possible mesh change.
*unVec = *solution;
#endif
*newUn = *unVec; *newUn = *unVec;
*lowSol = *unVec; *lowSol = *unVec;
...@@ -103,40 +112,42 @@ namespace AMDiS { ...@@ -103,40 +112,42 @@ namespace AMDiS {
// stage-solution: u_s(i) = u_old + sum_{j=0}^{i-1} a_ij*U_j // stage-solution: u_s(i) = u_old + sum_{j=0}^{i-1} a_ij*U_j
*stageSolution = *unVec; *stageSolution = *unVec;
for (int j = 0; j < i; j++) { for (int j = 0; j < i; j++) {
*tmp = *(stageSolutions[j]); *tmp = *(stageSolutions[j]);
*tmp *= rm->getA(i, j); *tmp *= rm->getA(i, j);
*stageSolution += *tmp; *stageSolution += *tmp;
} }
// Dirichlet-BC implemented as additional algebraic equation u = g(x,t) on boundary // 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 // => 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 // where u_s(i) = ith stage-solution, t_s(i) = ith stage-time
for (unsigned int j = 0; j < boundaries.size(); j++) { for (unsigned int j = 0; j < boundaries.size(); j++) {
boundaries[j].vec->interpol(boundaries[j].fct); boundaries[j].vec->interpol(boundaries[j].fct);
*(boundaries[j].vec) -= *(stageSolution->getDOFVector(boundaries[j].col)); *(boundaries[j].vec) -= *(stageSolution->getDOFVector(boundaries[j].col));
if (boundaries[j].fctDt != NULL) { if (boundaries[j].fctDt != NULL) {
// time derivative of dirichlet bc is given // time derivative of dirichlet bc is given
DOFVector<double> tmpDt(getFeSpace(boundaries[j].col), "tmp"); DOFVector<double> tmpDt(getFeSpace(boundaries[j].col), "tmp");
tmpDt.interpol(boundaries[j].fctDt); tmpDt.interpol(boundaries[j].fctDt);
tmpDt *= tauGammaI; tmpDt *= tauGammaI;
*(boundaries[j].vec) += tmpDt; *(boundaries[j].vec) += tmpDt;
} }
} }
// timeRhs: sum_{j=0}^{i-1} c_ij / tau * U_j // timeRhs: sum_{j=0}^{i-1} c_ij / tau * U_j
timeRhsVec->set(0.0); timeRhsVec->set(0.0);
for (int j = 0; j < i; j++) { for (int j = 0; j < i; j++) {
*tmp = *(stageSolutions[j]); *tmp = *(stageSolutions[j]);
*tmp *= (rm->getC(i, j) / *tauPtr); *tmp *= (rm->getC(i, j) / *tauPtr);
*timeRhsVec += *tmp; *timeRhsVec += *tmp;
} }
// assemble and solve stage equation // assemble and solve stage equation
ProblemStat::buildAfterCoarsen(adaptInfo, flag, (i == 0), true); Flag flag = BUILD | SOLVE;
ProblemStat::assemble(adaptInfo, flag, (i == 0), true);
#if defined HAVE_PARALLEL_PETSC #if defined HAVE_PARALLEL_PETSC
// TODO: Problems with reuse of Matrix with parallel PETSC-Solvers // NOTE: Problems with reuse of Matrix with parallel PETSC-Solvers
// Thus, Rosenbrock not efficient but working (Rainer) // Thus, Rosenbrock not efficient but working (Rainer)
// TODO: Change implementation
ProblemStat::solve(adaptInfo, true , false); ProblemStat::solve(adaptInfo, true , false);
#else #else
ProblemStat::solve(adaptInfo, i == 0, i + 1 < rm->getStages()); ProblemStat::solve(adaptInfo, i == 0, i + 1 < rm->getStages());
......
...@@ -162,8 +162,7 @@ namespace AMDiS { ...@@ -162,8 +162,7 @@ namespace AMDiS {
Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo) override; Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo) override;
virtual Flag stageIteration(AdaptInfo *adaptInfo, Flag flag, virtual Flag stageIteration(AdaptInfo *adaptInfo);
bool asmMatrix, bool asmVector);
virtual void estimateTimeError(AdaptInfo* adaptInfo); virtual void estimateTimeError(AdaptInfo* adaptInfo);
......
Markdown is supported
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