Commit 211268c7 authored by Lukas Riedel's avatar Lukas Riedel 📝

Completely switch to new error estimator

Error indicator now does not contain squares anymore.
Removed VTK output of error estimators.
Removed all unused variables inside the estimator.
parent c21549d7
......@@ -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");
......@@ -111,7 +108,7 @@ public:
do { // Loop for multiple refinements
Dune::Timer timer3;
float t_setup,t_est,t_sqrt,t_strtgy,t_mark,t_adapt;
float t_setup,t_est,t_strtgy,t_mark,t_adapt;
if(verbose>1)
std::cout << " Refinement Step " << multiRefinementCounter << ": ";
......@@ -119,64 +116,26 @@ public:
double eta_alpha(0.0);
double eta_beta(0.0);
// Compute error estimate eta
// set up helper GFS
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();
// Compute error estimate eta
estgo.residual(uold,eta);
t_est = timer3.elapsed();
timer3.reset();
using Dune::PDELab::Backend::native;
maxeta = native(eta_newer)[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_newer)[i] > maxeta)
maxeta = native(eta_newer)[i];
}
// compute largest error
maxeta = eta.infinity_norm();
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();
// Stop adaptivity if target error is reached
......@@ -222,7 +181,6 @@ public:
std::cout << std::setprecision(4) << std::scientific;
std::cout << "::: setup time" << std::setw(12) << t_setup << std::endl;
std::cout << "::: error estimation time" << std::setw(12) << t_est << std::endl;
std::cout << "::: unsquaring errors time" << std::setw(12) << t_sqrt << std::endl;
std::cout << "::: strategy application time" << std::setw(12) << t_strtgy << std::endl;
std::cout << "::: grid marking time" << std::setw(12) << t_mark << std::endl;
std::cout << "::: grid adaptation time" << std::setw(12) << t_adapt << std::endl;
......
......@@ -62,21 +62,14 @@ namespace Dune {
enum { doAlphaSkeleton = true };
Parameter& param;
EstimatorWeights::Type weights;
const RF penalty_factor;
int intorderadd;
int quadrature_factor;
const Impl impl;
const int intorderadd;
const int quadrature_factor;
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
......@@ -181,14 +174,6 @@ namespace Dune {
// compute jump in solution
const RF h_jump = u_s - u_n;
// compute weights
RF omega_s = 0.5, omega_n = 0.5;
if(weights == EstimatorWeights::weightsOn)
{
omega_s = satCond_n / (satCond_s + satCond_n);
omega_n = satCond_s / (satCond_s + satCond_n);
}
// apply gravity vector
gradu_s -= param.gravity();
gradu_n -= param.gravity();
......@@ -196,62 +181,33 @@ namespace Dune {
// penalty term
const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);
// compute numerical flux
const RF numFlux_s = satCond_s * ( - (gradu_s * n_F) + penalty * h_jump);
const RF numFlux_n = satCond_n * ( - (gradu_n * n_F) + penalty * h_jump);
// compute jump in physical flux
const RF jump = omega_s * condFactor_s * numFlux_s - omega_n * condFactor_n * numFlux_n;
// 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));
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;
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;
const RF total_jump = grad_jump/2.0 + head_jump;
const RF total_jump = grad_jump/2.0 + head_jump;
sum += total_jump * total_jump * factor;
}
sum += total_jump * total_jump * factor;
}
// correct for the number of neighbors
if (impl == Impl::old) {
sum *= 1. / quadrature_rule.size();
}
DF h_T_s = diameter(ig.inside().geometry());
DF h_T_n = diameter(ig.outside().geometry());
// what is this variable for?
DF h_T_s = 1.0;
DF h_T_n = 1.0;
if (impl == Impl::corrected || impl == Impl::renew) {
h_T_s = diameter(ig.inside().geometry());
h_T_n = diameter(ig.outside().geometry());
}
DF C_F_T_s = 1.0;
DF C_F_T_n = 1.0;
DF C_P_T = 1.0 / M_PI;
if (impl == Impl::corrected || impl == Impl::renew) {
C_F_T_s = h_T_s * ig.geometry().volume() / ig.inside().geometry().volume();
C_F_T_s *= C_P_T * (2.0 / dim + C_P_T);
C_F_T_n = h_T_n * ig.geometry().volume() / ig.outside().geometry().volume();
C_F_T_n *= C_P_T * (2.0 / dim + C_P_T);
}
DF C_F_T_s = h_T_s * ig.geometry().volume() / ig.inside().geometry().volume();
C_F_T_s *= C_P_T * (2.0 / dim + C_P_T);
DF C_F_T_n = h_T_n * ig.geometry().volume() / ig.outside().geometry().volume();
C_F_T_n *= C_P_T * (2.0 / dim + C_P_T);
// accumulate indicator
r_s.accumulate(lfsv_s,0, h_T_s * C_F_T_s * sum );
r_n.accumulate(lfsv_n,0, h_T_n * C_F_T_n * sum );
r_s.accumulate(lfsv_s,0, std::sqrt( h_T_s * C_F_T_s * sum ));
r_n.accumulate(lfsv_n,0, std::sqrt( h_T_n * C_F_T_n * sum ));
}
private:
......
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