Commit 2b3163c9 authored by Lukas Riedel's avatar Lukas Riedel 📝

completely switched to new error estimator implementation. removed vtk output of error estimates

parent f809a0b1
......@@ -58,8 +58,6 @@ private:
const float adaptivityThreshold; //!< Global error threshold below which no refinement is applied
const int verbose; //!< output verbosity of this object
int count;
public:
/// Initialize members from config file parameters.
......@@ -76,8 +74,7 @@ public:
alpha(inifile.get<float>("adaptivity.refinementFraction")),
beta(inifile.get<float>("adaptivity.coarseningFraction")),
adaptivityThreshold(inifile.get<float>("adaptivity.threshold")),
verbose(inifile.get<int>("output.verbose")),
count(0)
verbose(inifile.get<int>("output.verbose"))
{
if(strategy!="elementFraction" && strategy!="errorFraction" && strategy!="threshold" && strategy!="targetTolerance")
DUNE_THROW(Dune::IOError,"Valid values of adaptivity.markingStrategy: elementFraction, errorFraction, threshold, targetTolerance");
......@@ -121,22 +118,11 @@ public:
// Compute error estimate eta
AGFS p0gfs = AGFSHelper::create(gv);
ESTLOP estlop(inifile,param,Impl::old);
ESTLOP estlop(inifile,param);
MBE mbe(0);
ESTGO estgo(gfs,p0gfs,estlop,mbe);
U0 eta(p0gfs,0.0);
// compute new estimates
ESTLOP estlop_new(inifile,param,Impl::corrected);
ESTGO estgo_new(gfs,p0gfs,estlop_new,mbe);
U0 eta_new(p0gfs,0.0);
estgo_new.residual(uold,eta_new);
ESTLOP estlop_newer(inifile,param,Impl::renew);
ESTGO estgo_newer(gfs,p0gfs,estlop_newer,mbe);
U0 eta_newer(p0gfs,0.0);
estgo_newer.residual(uold,eta_newer);
t_setup = timer3.elapsed();
timer3.reset();
......@@ -149,32 +135,11 @@ public:
maxeta = native(eta)[0];
for (unsigned int i=0; i<eta.flatsize(); i++) {
native(eta)[i] = sqrt(native(eta)[i]); // eta contains squares
native(eta_new)[i] = sqrt(native(eta_new)[i]);
native(eta_newer)[i] = sqrt(native(eta_newer)[i]);
if (native(eta)[i] > maxeta) maxeta = native(eta)[i];
}
if (verbose>1)
std::cout << "Largest Local Error: " << maxeta << " " << std::endl;
// print both estimates into one file
Dune::VTKWriter<GV> vtkwriter(gv);
using ETADGF = typename Dune::PDELab::DiscreteGridFunction<AGFS,U0>;
ETADGF etadgf(p0gfs,eta);
ETADGF etadgf_new(p0gfs,eta_new);
ETADGF etadgf_newer(p0gfs,eta_newer);
auto etadgf_vtk = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<ETADGF>>(etadgf,"err_est_old");
auto etadgf_new_vtk = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<ETADGF>>(etadgf_new,"err_est_corrected");
auto etadgf_newer_vtk = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<ETADGF>>(etadgf_newer,"err_est_new");
vtkwriter.addCellData(etadgf_vtk);
vtkwriter.addCellData(etadgf_new_vtk);
vtkwriter.addCellData(etadgf_newer_vtk);
vtkwriter.write("err_est" + std::to_string(count));
++count;
// use new scheme
eta = eta_newer;
t_sqrt = timer3.elapsed();
timer3.reset();
......
......@@ -21,8 +21,6 @@ namespace Dune {
};
};
enum Impl { old, corrected, renew };
/// Local operator for residual-based error estimation.
/*
* A call to residual() of a grid operator space will assemble
......@@ -67,16 +65,11 @@ namespace Dune {
int intorderadd;
int quadrature_factor;
const Impl impl;
FluxErrorEstimator (const Dune::ParameterTree& config, Parameter& param_,
const Impl impl_,
EstimatorWeights::Type weights_ = EstimatorWeights::weightsOn,
int intorderadd_ = 2, int quadrature_factor_ = 2)
: param(param_), penalty_factor(config.get<RF>("dg.penaltyFactor")),
intorderadd(intorderadd_), quadrature_factor(quadrature_factor_),
impl(impl_)
intorderadd(intorderadd_), quadrature_factor(quadrature_factor_)
{ }
/// skeleton integral depending on test and ansatz functions
......@@ -108,8 +101,7 @@ namespace Dune {
RF sum(0.0);
// loop over quadrature points and integrate normal flux
const auto quadrature_rule = quadratureRule(gtface,intorder);
for (const auto& it : quadrature_rule)
for (const auto& it : quadratureRule(gtface,intorder))
{
// retrieve unit normal vector
const Dune::FieldVector<DF,dim> n_F = ig.unitOuterNormal(it.position());
......@@ -206,36 +198,24 @@ namespace Dune {
// integration factor
const RF factor = it.weight() * ig.geometry().integrationElement(it.position());
if (impl == Impl::old) {
sum += jump * jump * factor * factor;
}
else if (impl == Impl::corrected) {
sum += jump * jump * factor;
}
else if (impl == Impl::renew) {
const RF grad_jump = (satCond_s * condFactor_s * (n_F * gradu_s))
- (satCond_n * condFactor_n * (n_F * gradu_n));
// jump in physical flux
const RF grad_jump = (satCond_s * condFactor_s * (n_F * gradu_s))
- (satCond_n * condFactor_n * (n_F * gradu_n));
const RF gamma = 2.0 * satCond_s * condFactor_s * satCond_n * condFactor_n / ( satCond_s * condFactor_s + satCond_n * condFactor_n );
const RF head_jump = gamma * penalty * h_jump;
// harmonic average of conductivities
const RF gamma = 2.0 * satCond_s * condFactor_s * satCond_n * condFactor_n / ( satCond_s * condFactor_s + satCond_n * condFactor_n );
const RF total_jump = grad_jump + head_jump;
// jump in solution
const RF head_jump = gamma * penalty * h_jump;
sum += total_jump * total_jump * factor;
}
}
// correct for the number of neighbors
if (impl == Impl::old) {
sum *= 1. / quadrature_rule.size();
// accumulate total jump and integrate
const RF total_jump = grad_jump + head_jump;
sum += total_jump * total_jump * factor;
}
// what is this variable for?
DF h_T = 1.0;
if (impl == Impl::corrected || impl == Impl::renew) {
h_T = std::max( diameter(ig.inside().geometry()),
// geometric factor
DF h_T = std::max( diameter(ig.inside().geometry()),
diameter(ig.outside().geometry()) );
}
// accumulate indicator
r_s.accumulate(lfsv_s,0, h_T * sum );
......
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