From 5de011fba6b61d425ccae66cc567354e6cf46e2e Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Tue, 21 Dec 2010 22:46:47 +0000
Subject: [PATCH] RecoveryEstimator modified to have a structure similar to
ResidualEstimator
---
AMDiS/src/ProblemVec.cc | 4 +-
AMDiS/src/RecoveryEstimator.cc | 154 ++++++++++++++++++---------------
AMDiS/src/RecoveryEstimator.h | 42 +++++++--
AMDiS/src/UmfPackSolver.h | 19 ++--
4 files changed, 135 insertions(+), 84 deletions(-)
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index dca63831..41403824 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -772,12 +772,12 @@ namespace AMDiS {
assembleFlag);
}
- if (asmMatrix) {
+// if (asmMatrix) {
solverMatrix.setMatrix(*systemMatrix);
createPrecon();
INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
- }
+// }
#ifdef _OPENMP
INFO(info, 8)("buildAfterCoarsen needed %.5f seconds system time / %.5f seconds wallclock time\n",
diff --git a/AMDiS/src/RecoveryEstimator.cc b/AMDiS/src/RecoveryEstimator.cc
index 501ee0a2..67cde50d 100644
--- a/AMDiS/src/RecoveryEstimator.cc
+++ b/AMDiS/src/RecoveryEstimator.cc
@@ -15,12 +15,13 @@
namespace AMDiS {
- RecoveryEstimator::RecoveryEstimator(std::string name, DOFVector<double> *uh, int r)
+ RecoveryEstimator::RecoveryEstimator(std::string name, DOFVector<double> *uh_, int r)
: Estimator(name, r),
- uh_(uh),
- relative_(0),
+ uh(uh_),
+ relative(0),
C(1.0),
method(0),
+ addEstimationToOld(false),
feSpace(NULL),
f_vec(NULL),
f_scal(NULL),
@@ -29,14 +30,14 @@ namespace AMDiS {
{
FUNCNAME("RecoveryEstimator::constructor()");
- GET_PARAMETER(0, name + "->rec method", "%d", &method);
- GET_PARAMETER(0, name + "->rel error", "%d", &relative_);
+ GET_PARAMETER(0, name + "->rec method", "%d", &method); // 0, 1, or 2 (see Recovery.h)
+ GET_PARAMETER(0, name + "->rel error", "%d", &relative); // 0 or 1
GET_PARAMETER(0, name + "->C", "%f", &C);
C = C > 1e-25 ? sqr(C) : 1.0;
if (norm == H1_NORM) {
- feSpace = uh->getFeSpace();
+ feSpace = uh_->getFeSpace();
degree = feSpace->getBasisFcts()->getDegree();
if (degree <= 2 && C != 1.0) {
@@ -44,10 +45,10 @@ namespace AMDiS {
WAIT;
}
} else {
- degree = uh->getFeSpace()->getBasisFcts()->getDegree() + 1;
+ degree = uh_->getFeSpace()->getBasisFcts()->getDegree() + 1;
feSpace = FiniteElemSpace::provideFeSpace(NULL,
- Lagrange::getLagrange(uh->getFeSpace()->getMesh()->getDim(), degree),
- uh->getFeSpace()->getMesh(),
+ Lagrange::getLagrange(uh_->getFeSpace()->getMesh()->getDim(), degree),
+ uh_->getFeSpace()->getMesh(),
name + "->feSpace");
if (method == 2) {
@@ -66,72 +67,114 @@ namespace AMDiS {
rec_struct = new Recovery(norm, method);
}
-
- double RecoveryEstimator::estimate(double timestep)
+ void RecoveryEstimator::init(double ts)
{
- FUNCNAME("RecoveryEstimator::estimate()");
-
- const BasisFunction *basFcts = uh_->getFeSpace()->getBasisFcts();
- Mesh *mesh = uh_->getFeSpace()->getMesh();
+ FUNCNAME("RecoveryEstimator::init()");
+
+ basFcts = uh->getFeSpace()->getBasisFcts();
int dim = mesh->getDim();
- int dow = Global::getGeo(WORLD);
- double h1Norm2 = 0.0;
-
+ h1Norm2 = 0.0;
+
if (norm == H1_NORM) { // sets recovery gradient.
if (method == 2)
- rec_grd = rec_struct->recovery(uh_, f_vec, f_scal, aux_vec);
+ rec_grd = rec_struct->recovery(uh, f_vec, f_scal, aux_vec);
else
- rec_grd = rec_struct->recovery(uh_, feSpace, f_vec, f_scal, aux_vec);
+ rec_grd = rec_struct->recovery(uh, feSpace, f_vec, f_scal, aux_vec);
rec_basFcts = rec_grd->getFeSpace()->getBasisFcts();
} else { // sets higher-order recovery solution.
- rec_uh = rec_struct->recoveryUh(uh_, feSpace);
+ rec_uh = rec_struct->recoveryUh(uh, feSpace);
rec_basFcts = rec_uh->getFeSpace()->getBasisFcts();
}
int deg = 2 * std::max(basFcts->getDegree(), rec_basFcts->getDegree());
- Quadrature *quad = Quadrature::provideQuadrature(dim, deg);
- int nPoints = quad->getNumPoints();
+ quad = Quadrature::provideQuadrature(dim, deg);
+ nPoints = quad->getNumPoints();
- WorldVector<double> quad_pt;
- WorldVector<double> *grdAtQP = new WorldVector<double>[nPoints];
+// WorldVector<double> quad_pt;
+ grdAtQP = new WorldVector<double>[nPoints];
- mtl::dense_vector<WorldVector<double> > recoveryGrdAtQP(nPoints);
- mtl::dense_vector<double> uhAtQP(nPoints), recoveryUhAtQP(nPoints);
-
- FastQuadrature *quadFast =
- FastQuadrature::provideFastQuadrature(basFcts, *quad, INIT_PHI | INIT_GRD_PHI);
- FastQuadrature *rec_quadFast =
- FastQuadrature::provideFastQuadrature(rec_basFcts, *quad, INIT_PHI | INIT_GRD_PHI);
+ recoveryGrdAtQP = mtl::dense_vector<WorldVector<double> >(nPoints);
+ uhAtQP = mtl::dense_vector<double>(nPoints);
+ recoveryUhAtQP = mtl::dense_vector<double>(nPoints);
+ quadFast = FastQuadrature::provideFastQuadrature(basFcts, *quad, INIT_PHI | INIT_GRD_PHI);
+ rec_quadFast = FastQuadrature::provideFastQuadrature(rec_basFcts, *quad, INIT_PHI | INIT_GRD_PHI);
+
// clear error indicators
+ if(!addEstimationToOld) {
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
elInfo->getElement()->setEstimation(0.0, row);
elInfo = stack.traverseNext(elInfo);
}
+ }
+
+ est_sum = 0.0;
+ est_max = 0.0;
+ est_t_sum = 0.0;
+ est_t_max = 0.0;
+
+ traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA;
+
+ }
+
+ void RecoveryEstimator::exit(bool output)
+ {
+ FUNCNAME("RecoveryEstimator::exit()");
+
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ double send_est_sum = est_sum;
+ double send_est_max = est_max;
- // traverse mesh
- est_sum = est_max = est_t_sum = 0.0;
+ MPI::COMM_WORLD.Allreduce(&send_est_sum, &est_sum, 1, MPI_DOUBLE, MPI_SUM);
+ MPI::COMM_WORLD.Allreduce(&send_est_max, &est_max, 1, MPI_DOUBLE, MPI_MAX);
+#endif
+
+ // Computing relative errors
+ if (relative) {
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+ while (elInfo) {
+ double estEl = elInfo->getElement()->getEstimation(row);
+ estEl /= h1Norm2;
+ elInfo->getElement()->setEstimation(estEl, row);
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ est_max /= h1Norm2;
+ est_sum /= h1Norm2;
+ }
+
+ est_sum = sqrt(est_sum);
+
+ delete [] grdAtQP;
- elInfo = stack.traverseFirst(mesh, -1,
- Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS |
- Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA);
- while (elInfo) {
+ if (output) {
+ MSG("estimate for component %d = %.8e\n", row, est_sum);
+ }
+ }
+
+ void RecoveryEstimator::estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo)
+ {
+ FUNCNAME("RecoveryEstimator::estimateElement()");
+
Element *el = elInfo->getElement();
double det = elInfo->getDet();
double errEl = 0.0;
+ double estEl = el->getEstimation(row);
+ int dow = Global::getGeo(WORLD);
if (norm == H1_NORM) {
// get gradient and recovery gradient at quadrature points
- uh_->getGrdAtQPs(elInfo, NULL, quadFast, grdAtQP);
+ uh->getGrdAtQPs(elInfo, NULL, quadFast, grdAtQP);
rec_grd->getVecAtQPs(elInfo, NULL, rec_quadFast, recoveryGrdAtQP);
if (f_scal) {
if (aux_vec)
aux_vec->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
else
- uh_->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
+ uh->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
}
// calc h1 error
@@ -151,7 +194,7 @@ namespace AMDiS {
}
} else {
// get vector and recovery vector at quadrature points
- uh_->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
+ uh->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
rec_uh->getVecAtQPs(elInfo, NULL, rec_quadFast, recoveryUhAtQP);
// calc l2 error
@@ -159,12 +202,12 @@ namespace AMDiS {
errEl += quad->getWeight(i) * sqr(recoveryUhAtQP[i] - uhAtQP[i]);
}
- double estEl = C * det * errEl;
+ estEl += C * det * errEl;
el->setEstimation(estEl, row);
est_sum += estEl;
est_max = std::max(est_max, estEl);
- if (relative_) {
+ if (relative) {
double normEl = 0.0;
if (norm == H1_NORM) {
@@ -181,31 +224,6 @@ namespace AMDiS {
h1Norm2 += det * normEl;
}
-
- elInfo = stack.traverseNext(elInfo);
- }
-
- // Computing relative errors
- if (relative_) {
- elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
- while (elInfo) {
- double estEl = elInfo->getElement()->getEstimation(row);
- estEl /= h1Norm2;
- elInfo->getElement()->setEstimation(estEl, row);
- elInfo = stack.traverseNext(elInfo);
- }
-
- est_max /= h1Norm2;
- est_sum /= h1Norm2;
- }
-
- est_sum = sqrt(est_sum);
-
- delete [] grdAtQP;
-
- MSG("estimate = %.8e\n", est_sum);
-
- return est_sum;
}
}
diff --git a/AMDiS/src/RecoveryEstimator.h b/AMDiS/src/RecoveryEstimator.h
index 636c37b8..218a484a 100644
--- a/AMDiS/src/RecoveryEstimator.h
+++ b/AMDiS/src/RecoveryEstimator.h
@@ -59,19 +59,24 @@ namespace AMDiS {
};
/// constructor
- RecoveryEstimator(std::string name, DOFVector<double> *uh, int r = -1);
+ RecoveryEstimator(std::string name, DOFVector<double> *uh_, int r = -1);
/// destructor.
virtual ~RecoveryEstimator() {}
- /// implements \ref Estimator::estimate().
- virtual double estimate(double timestep = 0.0);
+ /// implements \ref Estimator::init(double).
+ void init(double ts);
+ /// implements \ref Estimator::estimateElement(ElInfo*, DualElInfo*).
+ void estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo);
+
+ /// implements \ref Estimator::exit(bool).
+ void exit(bool output);
/// Sets uh.
- inline void setUh(DOFVector<double> *uh)
+ inline void setUh(DOFVector<double> *uh_)
{
- uh_ = uh;
+ uh = uh_;
}
/// Sets f.
@@ -92,6 +97,10 @@ namespace AMDiS {
aux_vec = uh;
}
+ inline void setAddEstimationToOld(bool value)
+ {
+ addEstimationToOld = value;
+ }
/// Gets recovery gradient.
inline DOFVector<WorldVector<double> >* getRecGrd()
@@ -108,10 +117,10 @@ namespace AMDiS {
protected:
/// finite element solution
- DOFVector<double> *uh_;
+ DOFVector<double> *uh;
/// absolute or relative error?
- int relative_;
+ int relative;
/// constant for scaling the estimator
double C;
@@ -119,6 +128,8 @@ namespace AMDiS {
/// recovery method
int method;
+ bool addEstimationToOld;
+
/// Working finite element space
const FiniteElemSpace *feSpace;
@@ -143,6 +154,23 @@ namespace AMDiS {
/// Recovery structure.
Recovery *rec_struct;
+
+ /// Number of quadrature points.
+ int nPoints;
+
+ Quadrature *quad;
+ FastQuadrature *quadFast, *rec_quadFast;
+
+ /// Basis functions
+ const BasisFunction *basFcts;
+
+ double h1Norm2;
+
+ WorldVector<double> quad_pt;
+ WorldVector<double> *grdAtQP;
+ mtl::dense_vector<WorldVector<double> > recoveryGrdAtQP;
+ mtl::dense_vector<double> uhAtQP;
+ mtl::dense_vector<double> recoveryUhAtQP;
};
}
diff --git a/AMDiS/src/UmfPackSolver.h b/AMDiS/src/UmfPackSolver.h
index c8f6e5ed..25906f4a 100644
--- a/AMDiS/src/UmfPackSolver.h
+++ b/AMDiS/src/UmfPackSolver.h
@@ -87,19 +87,24 @@ namespace AMDiS {
else
solver = new mtl::matrix::umfpack::solver<matrix_type>(A, UMFPACK_STRATEGY_SYMMETRIC);
} else {
- if (!multipleRhs)
+ if (!multipleRhs) {
if (store_symbolic)
solver->update_numeric();
else
- solver->update();
+ solver->update();
+ }
}
int code = (*solver)(x, b);
- mtl::dense_vector<value_type> r(b);
- r -= A * x;
- residual = two_norm(r);
- std::cout << "UmfPackSolver: ||b-Ax|| = " << residual << "\n";
-
+ if (info > 0) {
+ mtl::dense_vector<value_type> r(b);
+ r -= A * x;
+ residual = two_norm(r);
+ std::cout << "UmfPackSolver: ||b-Ax|| = " << residual << "\n";
+ if (residual > tolerance) {
+ WARNING("Tolerance tol=%e could not be reached!\n",tolerance);
+ }
+ }
return code;
}
--
GitLab