Commit 80f5f335 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* Fixed bug in ResidualEstimator

parent 06e1522f
......@@ -205,7 +205,7 @@ namespace AMDiS {
problemTime_->initTimestep(adaptInfo_);
#ifdef _OPENMP
/*#ifdef _OPENMP
#pragma omp parallel sections
{
#pragma omp section
......@@ -215,9 +215,10 @@ namespace AMDiS {
oneTimestep();
}
#else
problemTime_->startDelayedTimestepCalculation();
*/
// problemTime_->startDelayedTimestepCalculation();
oneTimestep();
#endif
//#endif
problemTime_->closeTimestep(adaptInfo_);
......
......@@ -58,15 +58,21 @@ namespace AMDiS {
void DataCollector::fillAllData()
{
::std::cout << "START1" << ::std::endl;
if (!elementDataCollected_) {
startCollectingElementData();
}
::std::cout << "START2" << ::std::endl;
if (!periodicDataCollected_) {
startCollectingPeriodicData();
}
::std::cout << "START3" << ::std::endl;
if (!valueDataCollected_) {
startCollectingValueData();
}
::std::cout << "START4" << ::std::endl;
}
int DataCollector::startCollectingElementData()
......@@ -123,13 +129,12 @@ namespace AMDiS {
ElInfo *elInfo = stack.traverseFirst(mesh_,
level_,
traverseFlag_ | Mesh::FILL_COORDS);
while(elInfo) {
while (elInfo) {
if (!writeElem_ || writeElem_(elInfo))
addValueData(elInfo);
addValueData(elInfo);
elInfo = stack.traverseNext(elInfo);
}
// Remove all interpolation marks and, instead, set to each
// interpolation point its continous index starting from 0.
int i = 0;
......@@ -146,7 +151,7 @@ namespace AMDiS {
addInterpData(elInfo);
elInfo = stack.traverseNext(elInfo);
}
valueDataCollected_ = true;
return(0);
......
......@@ -197,14 +197,19 @@ namespace AMDiS {
if (delayWriting_) {
::std::cout << "REIN!" << ::std::endl; ::std::cout.flush();
if (writeTecPlotFormat || writeAMDiSFormat || writePeriodicFormat) {
ERROR_EXIT("Delay writing only supported for ParaView file format!\n");
}
::std::cout << "DC=" << dataCollectors_.size()<< ::std::endl;
for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
dataCollectors_[i]->fillAllData();
}
::std::cout << "BIS HIERHIN!" << ::std::endl; ::std::cout.flush();
writingIsDelayed_ = true;
delayedFilename_ = fn;
return;
......
......@@ -526,7 +526,7 @@ namespace AMDiS {
clock_t first = clock();
int iter = solver_->solve(matVec_, solution_, rhs_, leftPrecon_, rightPrecon_);
#ifdef _OPENMP
INFO(info_, 8)("solution of discrete system needed %.5f seconds system time / %.5f seconds wallclock time\n",
TIME_USED(first, clock()),
......
......@@ -77,6 +77,9 @@ namespace AMDiS {
riq = GET_MEMORY(double, numPoints);
grdUh_qp = NULL;
D2uhqp = NULL;
TraverseStack stack;
ElInfo *elInfo = NULL;
......@@ -117,12 +120,15 @@ namespace AMDiS {
}
FREE_MEMORY(uhEl, double*, numSystems);
if (timestep)
FREE_MEMORY(uhOldEl, double*, numSystems);
if (timestep) {
FREE_MEMORY(uhOldEl, double*, numSystems);
FREE_MEMORY(uhQP, double, numPoints);
FREE_MEMORY(uhOldQP, double, numPoints);
} else {
if (uhQP != NULL) {
FREE_MEMORY(uhQP, double, numPoints);
}
}
if (output) {
......@@ -135,6 +141,13 @@ namespace AMDiS {
FREE_MEMORY(riq, double, numPoints);
FREE_MEMORY(basFcts, const BasisFunction*, numSystems);
FREE_MEMORY(quadFast, FastQuadrature*, numSystems);
if (grdUh_qp != NULL) {
FREE_MEMORY(grdUh_qp, WorldVector<double>, numPoints);
}
if (D2uhqp != NULL) {
FREE_MEMORY(D2uhqp, WorldMatrix<double>, numPoints);
}
}
void ResidualEstimator::estimateElement(ElInfo *elInfo)
......@@ -150,9 +163,6 @@ namespace AMDiS {
::std::vector<Operator*>::iterator it;
WorldVector<double> *grdUh_qp = NULL;
WorldMatrix<double> *D2uhqp = NULL;
el = elInfo->getElement();
double det = elInfo->getDet();
......@@ -212,20 +222,17 @@ namespace AMDiS {
for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin();
it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd();
++it) {
if ((*it)->zeroOrderTerms() && !uhQP) {
if ((uhQP == NULL) && (*it)->zeroOrderTerms()) {
uhQP = GET_MEMORY(double, numPoints);
uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
}
if ((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi()
&& !grdUh_qp) {
if ((grdUh_qp == NULL) && ((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi())) {
grdUh_qp = new WorldVector<double>[numPoints];
uh[system]->getGrdAtQPs(elInfo, NULL, quadFast[system], grdUh_qp);
}
if ((*it)->secondOrderTerms() && !D2uhqp) {
if (degree > 2) {
D2uhqp = new WorldMatrix<double>[numPoints];
uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);
}
if ((D2uhqp == NULL) && (degree > 2) && (*it)->secondOrderTerms()) {
D2uhqp = new WorldMatrix<double>[numPoints];
uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);
}
}
......@@ -240,7 +247,7 @@ namespace AMDiS {
matrix[system],
fh[system],
quad,
riq);
riq);
}
}
......
......@@ -148,6 +148,10 @@ namespace AMDiS {
double *uhOldQP;
double *riq;
WorldVector<double> *grdUh_qp;
WorldMatrix<double> *D2uhqp;
};
}
......
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