Skip to content

GitLab

Commit 36d40082 authored by Lukas Riedel's avatar Lukas Riedel
Browse files
Options

Merge branch 'feature/error-estimator-revamp' into development/dorie-v1.0

parents d5b3df7d 61efc64a
Hide whitespace changes
Inline Side-by-side
Showing with 178 additions and 75 deletions
dune/dorie/interface/adaptivity.hh
View file @ 36d40082
......@@ -5,12 +5,13 @@ namespace Dune{
namespace Dorie{
/// Adaptivity base class. Does nothing.
template<typename Traits, typename GFS, typename Param, typename U>
template<typename Traits, typename GFS, typename Param, typename Boundary, typename U>
class AdaptivityHandlerBase
{
private:
using Grid = typename Traits::Grid;
using GV = typename Traits::GV;
using RF = typename Traits::RF;
public:
AdaptivityHandlerBase () = default;
......@@ -19,12 +20,16 @@ public:
/// Do nothing.
/** \return Boolean if operators have to be updated
*/
virtual bool adapt_grid (Grid& grid, GV& gv, GFS& gfs, Param& param, U& uold, U& unew) { return false; }
virtual bool adapt_grid (Grid& grid, GV& gv, GFS& gfs, const Param& param, const Boundary& boundary, const RF time, U& uold, U& unew) { return false; }
};
/// Specialized SimulationBase class providing functions and members for adaptive grid refinement
template<typename Traits, typename GFS, typename Param, typename U>
class AdaptivityHandler : public AdaptivityHandlerBase<Traits,GFS,Param,U>
template<typename Traits,
typename GFS,
typename Param,
typename Boundary,
typename U>
class AdaptivityHandler : public AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U>
{
private:
static constexpr int order = Traits::fem_order;
......@@ -37,7 +42,7 @@ private:
/// Adaptivity GFS
using AGFS = typename AGFSHelper::Type;
/// Error estimator local operator
using ESTLOP = Dune::Dorie::FluxErrorEstimator<Traits,Param,typename AGFSHelper::FEM>;
using ESTLOP = Dune::Dorie::FluxErrorEstimator<Traits,Param,Boundary>;
/// Empty constraints container
using NoTrafo = Dune::PDELab::EmptyTransformation;
/// Solution vector type
......@@ -66,7 +71,7 @@ public:
* \param grid Grid to adapt (reference is not saved)
*/
AdaptivityHandler (Dune::ParameterTree& _inifile, Grid& grid)
: AdaptivityHandlerBase<Traits,GFS,Param,U>(),
: AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U>(),
inifile(_inifile),
maxLevel(inifile.get<int>("adaptivity.maxLevel")),
minLevel(inifile.get<int>("adaptivity.minLevel")),
......@@ -95,7 +100,15 @@ public:
* \param unew Solution vector
* \return Boolean if grid operators have to be set up again (true if grid changed)
*/
bool adapt_grid (Grid& grid, GV& gv, GFS& gfs, Param& param, U& uold, U& unew) override
bool adapt_grid (
Grid& grid,
GV& gv,
GFS& gfs,
const Param& param,
const Boundary& boundary,
const RF time,
U& uold,
U& unew) override
{
Dune::Timer timer2;
float t2;
......@@ -108,7 +121,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_sqrt, t_est, t_strtgy, t_mark, t_adapt;
if(verbose>1)
std::cout << " Refinement Step " << multiRefinementCounter << ": ";
......@@ -116,9 +129,10 @@ 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);
ESTLOP estlop(inifile, param, boundary);
estlop.setTime(time);
MBE mbe(0);
ESTGO estgo(gfs,p0gfs,estlop,mbe);
U0 eta(p0gfs,0.0);
......@@ -126,17 +140,21 @@ public:
t_setup = timer3.elapsed();
timer3.reset();
// Compute error estimate eta
estgo.residual(uold,eta);
t_est = timer3.elapsed();
timer3.reset();
// unsquare errors
using Dune::PDELab::Backend::native;
maxeta = native(eta)[0];
for (unsigned int i=0; i<eta.flatsize(); i++) {
native(eta)[i] = sqrt(native(eta)[i]); // eta contains squares
if (native(eta)[i] > maxeta) maxeta = native(eta)[i];
}
auto& eta_nat = native(eta);
std::transform(eta_nat.begin(), eta_nat.end(), eta_nat.begin(),
[](const auto val) { return std::sqrt(val); }
);
// compute largest error
maxeta = eta.infinity_norm();
if (verbose>1)
std::cout << "Largest Local Error: " << maxeta << " " << std::endl;
......@@ -186,7 +204,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;
......@@ -211,14 +228,14 @@ public:
* The constexpr boolean adapt_grid in Dune::Dorie::BaseTraits will
* define which create() function is enabled at compile time.
*/
template<typename Traits, typename GFS, typename Param, typename U>
template<typename Traits, typename GFS, typename Param, typename Boundary, typename U>
struct AdaptivityHandlerFactory
{
private:
static constexpr bool enabled = Traits::adapt_grid;
using Grid = typename Traits::Grid;
using AHB = AdaptivityHandlerBase<Traits,GFS,Param,U>;
using AHB = AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U>;
Dune::ParameterTree& inifile;
Grid& grid;
......@@ -252,7 +269,7 @@ public:
std::is_same<Grid,Dune::UGGrid<3>>::value,
"Adaptivity is only implemented for UGGrid!");
std::unique_ptr<AHB> p;
p = std::make_unique<AdaptivityHandler<Traits,GFS,Param,U>>(inifile,grid);
p = std::make_unique<AdaptivityHandler<Traits,GFS,Param,Boundary,U>>(inifile,grid);
return p;
}
};
......
dune/dorie/interface/simulation.cc
View file @ 36d40082
......@@ -165,12 +165,13 @@ bool Simulation<Traits>::compute_time_step ()
template<typename Traits>
void Simulation<Traits>::run ()
{
output->write_vtk_output(controller->getTime());
const auto t_start = controller->getTime();
output->write_vtk_output(t_start);
while(controller->doStep()) {
if(!compute_time_step()){
continue;
}
if(adaptivity->adapt_grid(*grid, gv, *gfs, *param, *uold, *unew)){ // reset operators if grid changes
if(adaptivity->adapt_grid(*grid, gv, *gfs, *param, *fboundary, t_start, *uold, *unew)){ // reset operators if grid changes
operator_setup();
}
output->write_vtk_output(controller->getTime());
......
dune/dorie/interface/simulation.hh
View file @ 36d40082
......@@ -83,9 +83,9 @@ protected:
/// Factory for creating the output writer
using OutputWriterFactory = Dune::Dorie::OutputWriterFactory<Traits,GFS,Parameters,U>;
/// Grid Adaptivity Handler base class
using AdaptivityHandler = Dune::Dorie::AdaptivityHandlerBase<Traits,GFS,Parameters,U>;
using AdaptivityHandler = Dune::Dorie::AdaptivityHandlerBase<Traits,GFS,Parameters,FlowBoundary,U>;
/// Factory for creating the AdaptivityHandler
using AdaptivityHandlerFactory = Dune::Dorie::AdaptivityHandlerFactory<Traits,GFS,Parameters,U>;
using AdaptivityHandlerFactory = Dune::Dorie::AdaptivityHandlerFactory<Traits,GFS,Parameters,FlowBoundary,U>;
Dune::MPIHelper& helper;
std::shared_ptr<Grid> grid;
......
dune/dorie/solver/operator_error_indicator.hh
View file @ 36d40082
......@@ -13,24 +13,13 @@
namespace Dune {
namespace Dorie {
//! Weighting at grid interface for error estimation
struct EstimatorWeights
{
enum Type { weightsOn, //!< Use saturated conductivities at interface for weighting (default)
weightsOff //!< No weighting
};
};
/// Local operator for residual-based error estimation.
/*
* A call to residual() of a grid operator space will assemble
* the quantity \eta_T^2 for each cell. Note that the squares
* of the cell indicator \eta_T is stored. To compute the global
* error estimate sum up all values and take the square root.
* the quantity \eta_T for each cell.
* Skeleton integral is mostly the same as in the DG operator.
*
*/
template<typename Traits, typename Parameter, typename FEM>
template<typename Traits, typename Parameter, typename Boundary>
class FluxErrorEstimator
: public Dune::PDELab::LocalOperatorDefaultFlags
{
......@@ -49,27 +38,27 @@ namespace Dune {
typedef typename Traits::Domain Domain;
typedef typename Traits::IntersectionDomain ID;
typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
typedef Dune::PDELab::LocalBasisCache<LocalBasisType> Cache;
// pattern assembly flags
enum { doPatternVolume = false };
enum { doPatternSkeleton = false };
// residual assembly flags
enum { doAlphaSkeleton = true };
enum { doAlphaBoundary = true };
Parameter& param;
EstimatorWeights::Type weights;
const Parameter& param;
const Boundary& boundary;
const RF penalty_factor;
int intorderadd;
int quadrature_factor;
const int intorderadd;
const int quadrature_factor;
RF time;
FluxErrorEstimator (const Dune::ParameterTree& config, Parameter& param_,
EstimatorWeights::Type weights_ = EstimatorWeights::weightsOn,
FluxErrorEstimator (const Dune::ParameterTree& config,
const Parameter& param_, const Boundary& boundary_,
int intorderadd_ = 2, int quadrature_factor_ = 2)
: param(param_), penalty_factor(config.get<RF>("dg.penaltyFactor")),
intorderadd(intorderadd_), quadrature_factor(quadrature_factor_)
: param(param_), boundary(boundary_),
penalty_factor(config.get<RF>("dg.penaltyFactor")),
intorderadd(intorderadd_), quadrature_factor(quadrature_factor_)
{ }
/// skeleton integral depending on test and ansatz functions
......@@ -173,14 +162,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();
......@@ -188,42 +169,139 @@ 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());
// jump in physical flux
const RF grad_jump = (satCond_s * condFactor_s * (n_F * gradu_s))
- (satCond_n * condFactor_n * (n_F * gradu_n));
// 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 );
// jump in solution
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;
// accumulate total jump and integrate
const RF total_jump = grad_jump + head_jump;
const RF total_jump = grad_jump/2.0 + head_jump;
sum += total_jump * total_jump * factor;
}
// geometric factor
DF h_T = std::max( diameter(ig.inside().geometry()),
diameter(ig.outside().geometry()) );
DF h_T_s = diameter(ig.inside().geometry());
DF h_T_n = diameter(ig.outside().geometry());
DF C_P_T = 1.0 / M_PI;
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 * sum );
r_n.accumulate(lfsv_n,0, h_T * sum );
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 );
}
template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
void alpha_boundary (const IG& ig,
const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
R& r_s) const
{
// get polynomial degree
const int degree = lfsu_s.finiteElement().localBasis().order();
const int intorder = intorderadd + quadrature_factor * degree;
// geometric factor of penalty
const RF h_F = ig.inside().geometry().volume() / ig.geometry().volume(); // Houston!
// get element geometry
auto gtface = ig.geometry();
// container for sum of errors
RF sum(0.0);
// loop over quadrature points and integrate normal flux
for (const auto& it : quadratureRule(gtface,intorder))
{
// check boundary condition type
const auto bc = boundary.bc(ig.intersection(), it.position(), time);
if (!BoundaryCondition::isDirichlet(bc)){
continue;
}
// position of quadrature point in local coordinates of elements
const Domain local_s = ig.geometryInInside().global(it.position());
// retrieve unit normal vector
const Dune::FieldVector<DF,dim> n_F = ig.unitOuterNormal(it.position());
// evaluate basis functions
std::vector<Scalar> phi_s(lfsu_s.size());
lfsu_s.finiteElement().localBasis().evaluateFunction(local_s,phi_s);
// evaluate u
RF u_s = 0.;
for (unsigned int i = 0; i<lfsu_s.size(); i++)
u_s += x_s(lfsu_s,i) * phi_s[i];
// evaluate gradient of basis functions (we assume Galerkin method lfsu = lfsv)
typedef typename LFSU::Traits::FiniteElementType::Traits
::LocalBasisType::Traits::JacobianType JacobianType;
std::vector<JacobianType> gradphi_s(lfsu_s.size());
lfsu_s.finiteElement().localBasis().evaluateJacobian(local_s,gradphi_s);
// transform gradients to real element
Dune::FieldMatrix<DF,dim,dim> jac;
std::vector<Vector> tgradphi_s(lfsu_s.size());
jac = ig.inside().geometry().jacobianInverseTransposed(local_s);
for (unsigned int i = 0; i<lfsu_s.size(); i++)
jac.mv(gradphi_s[i][0],tgradphi_s[i]);
// compute gradient of u
Vector gradu_s(0.);
for (unsigned int i = 0; i<lfsu_s.size(); i++)
gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
// offset for parameter queries
const RF eps = 1e-9;
// parameter queries: inside
Domain xGlobal_s = ig.inside().geometry().global(local_s);
xGlobal_s = xGlobal_s.axpy(-eps,n_F);
const RF head_s = param.headMillerToRef(u_s,xGlobal_s);
const RF satCond_s = param.condRefToMiller(param.K(xGlobal_s),xGlobal_s);
const RF saturation_s = param.saturation(head_s,xGlobal_s);
const RF condFactor_s = param.condFactor(saturation_s,xGlobal_s);
// boundary condition at position
const RF g = boundary.g(ig.intersection(), it.position(), time);
// compute jump in solution
const RF h_jump = u_s - g;
// apply gravity vector
gradu_s -= param.gravity();
// penalty term
const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);
// integration factor
const RF factor = it.weight() * ig.geometry().integrationElement(it.position());
const RF head_jump = satCond_s * condFactor_s * penalty * h_jump;
sum += head_jump * head_jump * factor;
}
DF h_T_s = diameter(ig.inside().geometry());
DF C_P_T = 1.0 / M_PI;
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);
// accumulate indicator
r_s.accumulate(lfsv_s,0, h_T_s * C_F_T_s * sum );
}
private:
// probably redundant
//! compute diameter of a grid cell
template<class GEO>
typename GEO::ctype diameter (const GEO& geo) const
{
......@@ -243,6 +321,13 @@ namespace Dune {
return hmax;
}
public:
//! set operator time
void setTime (RF t)
{
time = t;
}
};
}
......
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