Skip to content
GitLab
Commit 5de011fb authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

RecoveryEstimator modified to have a structure similar to ResidualEstimator

parent 04249ebf
Hide whitespace changes
Inline Side-by-side
AMDiS/src/ProblemVec.cc
View file @ 5de011fb
......@@ -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",
......
AMDiS/src/RecoveryEstimator.cc
View file @ 5de011fb
......@@ -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;
}
}
AMDiS/src/RecoveryEstimator.h
View file @ 5de011fb
......@@ -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;
};
}
......
AMDiS/src/UmfPackSolver.h
View file @ 5de011fb
......@@ -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;
}
......
Supports Markdown
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