Commit e9313900 by Praetorius, Simon

### some new baseproblems

parent 7d720fce
 ... ... @@ -72,9 +72,6 @@ namespace AMDiS { void reset(); protected: /// Pointer to the Rosenbrock method that should be used. RosenbrockMethod *rosenbrockMethod; /// Pointer to the stationary problem; ProblemStationary *rosenbrockStat; ... ... @@ -84,7 +81,6 @@ namespace AMDiS { /// If true, the last timestep was rejected. bool lastTimestepRejected; /// If true, more than one of the last timesteps were rejected. bool succRejection; ... ... @@ -95,24 +91,12 @@ namespace AMDiS { /// Timestep of the last accepted time iteration. double tauAcc; /// Error estimation of the last accepted time iteration. double estAcc; /// Timestep of the last rejected time iteration. double tauRej; /// Error estimation of the last rejected time iteration. double estRej; /// Current timestep. double tau; /// The value tau * gamma, where gamma is a value of the used /// Rosenbrock method. double tauGamma, minusTauGamma, invTauGamma, minusInvTauGamma; double oldTime; /// If true, the first timestep is calculated with different timesteps. /// This is usually used to make a study how the time error estimator /// behavous for different timesteps. ... ... @@ -121,8 +105,6 @@ namespace AMDiS { /// If \ref dbgTimestepStudy is set to true, then this array contains the /// timesteps for which the first timestep should be calculated. vector dbgTimesteps; vector errorEstWeights; }; ... ...
 ... ... @@ -31,15 +31,11 @@ namespace AMDiS { succRejection(false), fixFirstTimesteps(0), tau(1.0), tauGamma(1.0), minusTauGamma(-1.0), invTauGamma(1.0), minusInvTauGamma(-1.0), oldTime(adaptInfo->getStartTime()), dbgTimestepStudy(false) { FUNCNAME("ExtendedRosenbrockAdaptInstationary::ExtendedRosenbrockAdaptInstationary()"); initConstructor(problemStat); rosenbrockStat->setOldTime(adaptInfo->getStartTime()); } ... ... @@ -57,14 +53,11 @@ namespace AMDiS { succRejection(false), fixFirstTimesteps(0), tau(1.0), tauGamma(1.0), minusTauGamma(-1.0), invTauGamma(1.0), minusInvTauGamma(-1.0), dbgTimestepStudy(false) { FUNCNAME("ExtendedRosenbrockAdaptInstationary::ExtendedRosenbrockAdaptInstationary()"); initConstructor(&problemStat); rosenbrockStat->setOldTime(adaptInfo->getStartTime()); } ... ... @@ -73,43 +66,12 @@ namespace AMDiS { { FUNCNAME("ExtendedRosenbrockAdaptInstationary::initConstructor()"); std::string str(""); std::string initFileStr(name + "->rosenbrock method"); Parameters::get(initFileStr, str); RosenbrockMethodCreator *creator = dynamic_cast(CreatorMap::getCreator(str, initFileStr)); rosenbrockMethod = creator->create(); TEST_EXIT_DBG(rosenbrockMethod)("This should not happen!\n"); Parameters::get(name + "->fix first timesteps", fixFirstTimesteps); rosenbrockStat->setRosenbrockMethod(rosenbrockMethod); adaptInfo->setRosenbrockMode(true); rosenbrockStat->setTauGamma(&tauGamma, &minusTauGamma, &invTauGamma, &minusInvTauGamma, &oldTime); rosenbrockStat->setTau(&tau); Parameters::get(name + "->fix first timesteps", fixFirstTimesteps); Parameters::get(name + "->rosenbrock->timestep study", dbgTimestepStudy); Parameters::get(name + "->rosenbrock->timestep study steps", dbgTimesteps); Parameters::get(name + "->rosenbrock->error weights", errorEstWeights); if (errorEstWeights.size() == 0) { WARNING("No weights for time-error estimation given.\n"); } errorEstWeights.resize(adaptInfo->getSize(), 0.0); double SumErrorEstWeights = vector_operations::sum(errorEstWeights); if (SumErrorEstWeights(adaptInfo->getSize()); } else { for (size_t i = 0; i < errorEstWeights.size(); i++) errorEstWeights[i] *= 1.0/SumErrorEstWeights; } } ... ... @@ -122,15 +84,9 @@ namespace AMDiS { lastTimestepRejected = false; succRejection = false; fixFirstTimesteps = 0; tau = 1.0; tauGamma = 1.0; minusTauGamma = -1.0; invTauGamma = 1.0; minusInvTauGamma = -1.0; oldTime = 0.0; } /* template void ExtendedRosenbrockAdaptInstationary::oneTimestep() { ... ... @@ -160,11 +116,6 @@ namespace AMDiS { adaptInfo->setTime(adaptInfo->getTime() + adaptInfo->getTimestep()); problemTime->setTime(adaptInfo); tau = adaptInfo->getTimestep(); tauGamma = tau * rosenbrockMethod->getGamma(); minusTauGamma = -tauGamma; invTauGamma = 1.0 / (tau * rosenbrockMethod->getGamma()); minusInvTauGamma = -1.0; INFO(info, 6)("time = %e timestep = %e\n", adaptInfo->getTime(), adaptInfo->getTimestep()); ... ... @@ -176,13 +127,10 @@ namespace AMDiS { problemIteration->oneIteration(adaptInfo, FULL_ITERATION); problemIteration->endIteration(adaptInfo); double errorEst = 0.0; for (size_t i = 0; i < adaptInfo->getSize(); i++) { errorEst += errorEstWeights[i] * adaptInfo->getTimeEstSum(i); } double errorEst = adaptInfo->getTimeEst(); INFO(info, 6)("overall time estimate = %e\n", errorEst); double newTimestep = 0.0; double order = rosenbrockMethod->getOrder(); if (errorEst < timeTol || fixFirstTimesteps > 0 || firstTimestep) { double fac = 1.0; ... ... @@ -190,39 +138,23 @@ namespace AMDiS { if (fixFirstTimesteps > 0) { newTimestep = adaptInfo->getTimestep(); } else { if (firstTimestep || succRejection) { fac = pow((timeTol / errorEst), 1.0 / order); } else { fac = adaptInfo->getTimestep() / tauAcc * pow((timeTol * estAcc / (errorEst * errorEst)), 1.0 / order); } fac = std::min(fac, 3.0); newTimestep = fac * adaptInfo->getTimestep(); newTimestep *= 0.95; if (firstTimestep || succRejection) newTimestep = rosenbrockStat->getNewTimestep(timeTol, adaptInfo->getTimestep()); else newTimestep = rosenbrockStat->getNewTimestep(timeTol, tauAcc, adaptInfo->getTimestep()); } tauAcc = adaptInfo->getTimestep(); estAcc = errorEst; lastTimestepRejected = false; succRejection = false; } else { if (lastTimestepRejected) { succRejection = true; double reducedP = log(errorEst / estRej) / log(adaptInfo->getTimestep() / tauRej); if (reducedP < order && reducedP > 0.0) order = reducedP; succRejection = true; rosenbrockStat->reduceOrder(tauRej, adaptInfo->getTimestep()); } newTimestep = pow((timeTol / errorEst), 1.0 / order) * adaptInfo->getTimestep(); newTimestep *= 0.95; newTimestep = rosenbrockStat->getNewTimestep(timeTol, adaptInfo->getTimestep(), false); tauRej = adaptInfo->getTimestep(); estRej = errorEst; tauRej = adaptInfo->getTimestep(); lastTimestepRejected = true; } ... ... @@ -267,12 +199,77 @@ namespace AMDiS { } while (rejected || (dbgTimestepStudy && (studyTimestep + 1 < dbgTimesteps.size()))); rosenbrockStat->acceptTimestep(); oldTime = adaptInfo->getTime(); rosenbrockStat->acceptTimestep(adaptInfo); adaptInfo->setLastProcessedTimestep(adaptInfo->getTimestep()); adaptInfo->incTimestepNumber(); }*/ template void ExtendedRosenbrockAdaptInstationary::oneTimestep() { FUNCNAME("ExtendedRosenbrockAdaptInstationary::oneTimestep()"); // estimate before first adaption if (adaptInfo->getTime() <= adaptInfo->getStartTime()) problemIteration->oneIteration(adaptInfo, ESTIMATE); double timeTol = adaptInfo->getTimeTolerance(0); // increment time adaptInfo->setTime(adaptInfo->getTime() + adaptInfo->getTimestep()); problemTime->setTime(adaptInfo); tau = adaptInfo->getTimestep(); INFO(info, 6)("time = %e timestep = %e\n", adaptInfo->getTime(), adaptInfo->getTimestep()); adaptInfo->setSpaceIteration(0); // do the iteration problemIteration->beginIteration(adaptInfo); problemIteration->oneIteration(adaptInfo, FULL_ITERATION); problemIteration->endIteration(adaptInfo); double errorEst = adaptInfo->getTimeEst(); INFO(info, 6)("overall time estimate = %e\n", errorEst); double newTimestep = 0.0; double fac = 1.0; if (fixFirstTimesteps > 0) { newTimestep = adaptInfo->getTimestep(); } else { if (firstTimestep) newTimestep = rosenbrockStat->getNewTimestep(timeTol, adaptInfo->getTimestep()); else newTimestep = rosenbrockStat->getNewTimestep(timeTol, tauAcc, adaptInfo->getTimestep()); } tauAcc = adaptInfo->getTimestep(); for (int i = 0; i < adaptInfo->getSize(); i++) INFO(info, 6)("time estimate for component %d = %e\n", i, adaptInfo->getTimeEstSum(i)); INFO(info, 6)("Accepted timestep at time = %e with timestep = %e\n", adaptInfo->getTime() - adaptInfo->getTimestep(), adaptInfo->getTimestep()); adaptInfo->setTimestep(newTimestep); if (errorEst > timeTol) MSG("Accepted timestep but tolerance not reached \n"); if (firstTimestep) firstTimestep = false; if (fixFirstTimesteps > 0) fixFirstTimesteps--; rosenbrockStat->acceptTimestep(adaptInfo); adaptInfo->setLastProcessedTimestep(adaptInfo->getTimestep()); adaptInfo->incTimestepNumber(); } }
 ... ... @@ -12,11 +12,11 @@ #include "time/ExtendedRosenbrockStationary.h" #include "io/VtkWriter.h" #include "ExtendedRosenbrockStationary.h" #include "SystemVector.h" #include "OEMSolver.h" #include "Debug.h" #include "POperators_ZOT.h" #include "VectorOperations.h" #ifdef HAVE_PARALLEL_DOMAIN_AMDIS #include "parallel/MeshDistributor.h" ... ... @@ -24,12 +24,43 @@ using namespace AMDiS; void ExtendedRosenbrockStationary::acceptTimestep() { *solution = *newUn; *unVec = *newUn; } void ExtendedRosenbrockStationary::initialize(Flag initFlag, ProblemStat *adoptProblem, Flag adoptFlag) { FUNCNAME_DBG("ExtendedRosenbrockStationary::initialize()"); super::initialize(initFlag, adoptProblem, adoptFlag); // create rosenbrock method std::string str(""); std::string initFileStr(name + "->rosenbrock->method"); Parameters::get(initFileStr, str); RosenbrockMethodCreator *creator = dynamic_cast(CreatorMap::getCreator(str, initFileStr)); rm = creator->create(); TEST_EXIT_DBG(rm)("No Rosenbrock method created!\n"); init(); // read error weights Parameters::get(name + "->rosenbrock->error weights", errorEstWeights); if (errorEstWeights.size() == 0) { WARNING("No weights for time-error estimation given.\n"); } errorEstWeights.resize(super::getNumComponents(), 0.0); double SumErrorEstWeights = vector_operations::sum(errorEstWeights); if (SumErrorEstWeights(super::getNumComponents()); } else { for (size_t i = 0; i < errorEstWeights.size(); i++) errorEstWeights[i] *= 1.0/SumErrorEstWeights; } std::cout << "errorEstWeights = " << errorEstWeights << "\n"; } void ExtendedRosenbrockStationary::init() { ... ... @@ -49,13 +80,43 @@ using namespace AMDiS; stageSolutions[i] = new SystemVector(*solution); stageSolutions[i]->set(0.0); } order = rm->getOrder(); } Flag ExtendedRosenbrockStationary::oneIteration(AdaptInfo *adaptInfo, Flag toDo) { order = rm->getOrder(); Flag flag = 0, markFlag = 0; if (toDo.isSet(MARK)) markFlag = problem->markElements(adaptInfo); else markFlag = 3; // refine if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_REFINED)) flag = problem->refineMesh(adaptInfo); // coarsen if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_COARSENED)) flag |= problem->coarsenMesh(adaptInfo); if (toDo.isSet(BUILD) || toDo.isSet(SOLVE)) { flag = stageIteration(adaptInfo, toDo, true, true); estimateTimeError(adaptInfo); } if (toDo.isSet(ESTIMATE)) problem->estimate(adaptInfo); void ExtendedRosenbrockStationary::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag, return flag; } Flag ExtendedRosenbrockStationary::stageIteration(AdaptInfo *adaptInfo, Flag flag, bool asmMatrix, bool asmVector) { FUNCNAME("ExtendedRosenbrockStationary::buildAfterCoarsen()"); { FUNCNAME("ExtendedRosenbrockStationary::stageIteration()"); TEST_EXIT(tauPtr)("No tau pointer defined in stationary problem!\n"); ... ... @@ -67,10 +128,12 @@ using namespace AMDiS; *newUn = *unVec; if (!fixedTimestep) *lowSol = *unVec; invTauGamma = 1.0 / ((*tauPtr) * rm->getGamma()); for (int i = 0; i < rm->getStages(); i++) { stageTime = (*oldTime) + rm->getAlphaI(i) * (*tauPtr); 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 ... ... @@ -107,12 +170,6 @@ using namespace AMDiS; super::buildAfterCoarsen(adaptInfo, flag | UPDATE_DIRICHLET_BC, (i == 0), true); // mtl::io::matrix_market_ostream out("matrix.mat"); // SolverMatrix > solverMatrix; // solverMatrix.setMatrix(*getSystemMatrix()); // out << solverMatrix.getMatrix(); // out.close(); super::solve(adaptInfo, i == 0, true); *(stageSolutions[i]) = *solution; ... ... @@ -130,21 +187,53 @@ using namespace AMDiS; } } stageTime = (*oldTime) + (*tauPtr); stageTime = oldTime + (*tauPtr); Flag flag_; return flag_; } double ExtendedRosenbrockStationary::estimateTimeError(AdaptInfo* adaptInfo) { oldErrorEst = errorEst; errorEst = 0.0; if (!fixedTimestep) { for (int i = 0; i < nComponents; i++) { for (int i = 0; i < lowSol->getSize(); i++) { (*(lowSol->getDOFVector(i))) -= (*(newUn->getDOFVector(i))); adaptInfo->setTimeEstSum(lowSol->getDOFVector(i)->l2norm(), i+componentShift); errorEst += errorEstWeights[i] * adaptInfo->getTimeEstSum(i+componentShift); } } adaptInfo->setTimeEst(errorEst); return errorEst; } void ExtendedRosenbrockStationary::acceptTimestep(AdaptInfo* adaptInfo) { *solution = *newUn; *unVec = *newUn; oldTime = adaptInfo->getTime(); } void ExtendedRosenbrockStationary::solve(AdaptInfo *adaptInfo, bool, bool) {} double ExtendedRosenbrockStationary::getNewTimestep(double tol, double tau, bool restrict) { double fac = pow((tol / errorEst), 1.0 / order); if (restrict) fac = std::min(3.0, fac); return 0.95 * (fac * tau); } double ExtendedRosenbrockStationary::getNewTimestep(double tol, double oldTau, double tau, bool restrict) { double fac = tau / oldTau * pow((tol * oldErrorEst / sqr(errorEst)), 1.0 / order); if (restrict) fac = std::min(3.0, fac); return 0.95 * (fac * tau); } void ExtendedRosenbrockStationary::addOperator(Operator &op, int row, int col, double *factor, double *estFactor) { ... ... @@ -173,11 +262,9 @@ using namespace AMDiS; { FUNCNAME("ExtendedRosenbrockStationary::addTimeOperator()"); TEST_EXIT(invTauGamma)("This should not happen!\n"); Operator *op = new Operator(componentSpaces[row], componentSpaces[col]); op->addTerm(new Simple_ZOT(factor)); super::addMatrixOperator(op, row, col, invTauGamma, invTauGamma); super::addMatrixOperator(op, row, col, &invTauGamma, &invTauGamma); Operator *opRhs = new Operator(componentSpaces[row]); opRhs->addTerm(new Phase_ZOT(timeRhsVec->getDOFVector(col), factor)); ... ... @@ -200,10 +287,10 @@ using namespace AMDiS; void ExtendedRosenbrockStationary::addSingularDirichletBC(WorldVector &pos, int row, int col, AbstractFunction > &fct) AbstractFunction > *fct) { DOFVector* vec = new DOFVector(componentSpaces[col], "vec"); MyRosenbrockBoundary bound(&fct, NULL, vec, row, col); MyRosenbrockBoundary bound(fct, NULL, vec, row, col); boundaries.push_back(bound); super::addSingularDirichletBC(pos, row, col, *vec); ... ... @@ -211,10 +298,10 @@ using namespace AMDiS; void ExtendedRosenbrockStationary::addSingularDirichletBC(DegreeOfFreedom idx, int row, int col, AbstractFunction > &fct) AbstractFunction > *fct) { DOFVector* vec = new DOFVector(componentSpaces[col], "vec"); MyRosenbrockBoundary bound(&fct, NULL, vec, row, col); MyRosenbrockBoundary bound(fct, NULL, vec, row, col); boundaries.push_back(bound); super::addSingularDirichletBC(idx, row, col, *vec); ... ...
 ... ... @@ -69,53 +69,30 @@ using namespace AMDiS; componentShift(componentShift_), first(true), fixedTimestep(false), order(0), minusOne(-1.0), tauPtr(NULL), tauGamma(NULL), minusTauGamma(NULL), invTauGamma(NULL), minusInvTauGamma(NULL), stageTime(0.0) invTauGamma(1.0), stageTime(0.0), errorEst(0.0), oldErrorEst(0.0) {} DOFVector* getUnVec(int i) { return unVec->getDOFVector(i); } DOFVector* getStageSolution(int i) { return stageSolution->getDOFVector(i); } DOFVector* getTimeRhsVec(int i) { return timeRhsVec->getDOFVector(i); } void acceptTimestep(); void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag, bool asmMatrix, bool asmVector); virtual void initialize(Flag initFlag, ProblemStat *adoptProblem = NULL, Flag adoptFlag = INIT_NOTHING); void solve(AdaptInfo *adaptInfo, bool a = true, bool b = false); Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo); void setRosenbrockMethod(RosenbrockMethod *method) { rm = method; init(); } Flag stageIteration(AdaptInfo *adaptInfo, Flag flag, bool asmMatrix, bool asmVector); void setFixedTimestep(bool fixed) { fixedTimestep = fixed; } virtual double estimateTimeError(AdaptInfo* adaptInfo); void reset() { first = true; stageSolution->set(0.0); unVec->set(0.0); for (int i = 0; i < rm->getStages(); i++) stageSolutions[i]->set(0.0); } void acceptTimestep(AdaptInfo* adaptInfo); double getNewTimestep(double tol, double tau, bool restrict = true); double getNewTimestep(double tol, double oldTau, double tau, bool restrict = true); void addOperator(Operator &op, int row, int col,