Algorithms for the solution of nonlinear problems often rely on adaptive grid/mesh refinement. Triangulation of the geometry is usually chosen a priori, but by adaptive mesh refinement, based on the solution behavior the triangulation will be refined until a certain tolerance is met. A coarse grid is refined by splitting elements into smaller ones and in the same way, fine grids can be coarsened. Multiple refinement steps can be performed, if necessary. Refinement is performed directly by Kaskade 7.
In numerical analysis, adaptive mesh refinement (AMR) is a method of adapting the accuracy of a solution within certain sensitive or turbulent regions of simulation, dynamically and during the time the solution is being calculated. – Wikipedia entry
Hence, error estimation and adaptive grid refinement is intertwined in Kaskade 7. First, the error in the given finite element approximation is estimated a posteriori. Afterwards, an adaptive mesh refinement method refines elements based on these error estimates. The error estimator takes care of grid refinement as long as the estimated error exceeds a certain threshold. Beware of singularities in the geometry (e.g. due to sharp corners) since this can lead to unnecessary refinements.
In general, these steps are repeated until the desired tolerance is met:
Adaptive mesh refinement can improve the convergence of the error or decrease the number of degrees of freedom required to achieve a certain accuracy.
Have a look at tutorial xyz to get an impression of how to use the (embedded) error estimator in a practical problem with adaptive grid refinement.
The general setup to illustrate the used Kaskade 7 methods in the subsequent two sections is as follows:
/* define necessary parameters */
int verbosity = 1; // for debugging output
int refinements = 1; // initial grid refinement
int order = 2; // finite element ansatz order
int dim = 3;
/* create grid */
using Grid = Dune::UGGrid<dim>;
GridManager<Grid> gridManager(createUnitCube<Grid>(0.5));
gridManager.globalRefine(refinements);
/* define spaces and variables */
using H1Space = FEFunctionSpace<ContinuousLagrangeMapper<double,Grid::LeafGridView> >;
using Spaces = boost::fusion::vector<H1Space const*>;
using VariableSet = VariableSetDescription<Spaces,VariableDescriptions>;
H1Space temperatureSpace(gridManager,gridManager.grid().leafGridView(),order);
Spaces spaces(&temperatureSpace);
std::string varNames[1] = {"u"};
VariableSet ansatVars(spaces,varNames);
/* set tolerances */
double aTolX = 1e-3; // absolute tolerance
double rTolX = 1e-3; // relative tolerance
bool accurate = true; // marks if desired accuracy is reached
std::vector<std::pair<double,double> > tolX(ansatzVars.noOfVariables);
for (int i=0; i<tolX.size(); ++i) {
tolX[i] = std::make_pair(aTolX,rTolX);
}
VariableSet::VariableSet dx; // approximate solution, to be calculated
/* solution */
do {
[... SOLVE EQUATION ...]
[... ESTIMATE ERROR ...]
[... REFINE MESH ...]
} while (!accurate);
To control the discretization error, the calculation of the solution of a system is usually followed by an estimation of the discretization error and by adaptive mesh refinement as long as the estimated error exceeds a threshold. For this, first the cells have to be identified where mesh refinement is necessary. Second, these cells are marked. In the end, the mesh is refined at the selected cells.
This snippet depicts the use of the projectHierarchically method which yields a representation of the estimated error:
err = dx;
projectHierarchically(ansatzVars,err);
err -= dx; // this is a representation of the estimated error
The solution dx is projected onto the ansatz space with polynomials of order one less.
The embedded error estimator needs an approximation of order p. Then the difference between the order p-1 projection from projectHierarchically and the order p solution yields a measure for the error. For instance, providing a solution with P2 elements, the embedded part of order 1 (P1) is substracted. Therefore, the embedded error estimator is not applicable for linear finite elements - use the hierarchic error estimator in the next section instead.
This snippet depicts the use of the embeddedErrorEstimator method which refines the grid if the tolerances are not met, having a representation err of the estimated error:
accurate = embeddedErrorEstimator(ansatzVars,err,dx,IdentityScaling(),tolX,gridManager,verbosity);
This method returns true if the error is below the tolerance and the mesh hat not been refined. Depending on the underlying equation, other scaling methods than IdentityScaling() can be implemented.
The embedded error estimator marks cells for refinement before adapting the grid by use of the GridManager class.
Two important parameters for controlling the adaptive process are the absolute tolerance atolX and the relative tolerance rTolX. Sufficient accuracy is assumed if the local error \( \epsilon \) of the approximate solution \( u \) satisfies \( |\epsilon|^2 \leq \mathrm{aTolX}^2 + \mathrm{rTolX}^2 |u|^2 \) (and the normalized error respectively) and the adaptive refinement is stopped.
The use of the hierarchic error estimator is more intricate. In contrary to the embedded error estimator, it can be applied to linear (P1) finite elements. We use error estimates from an extended space, as shown in BornemannErdmannKornhuber1996. For a continuous problem with \( u \in H \) \[ a(u,v) = (f,v) \quad \forall v \in H=H_0^1(\Omega) \] the resulting finite element approximation in space \( S \) is given by \[ a(u_S,v) = (f,v) \quad \forall v\in S \] with \( u_S \in S \). We are interested in the a posteriori energy norm \( || e || \) of the error \( e = u-u_S \). This error solves the defect equation for \( e \in H \) \[ a(e,v) = (f,v)-a(u_S,v) \quad \forall v\in H \] and is approximated by a larger finite element space \( Q = S+V \). The solution \( e_Q\in Q \) of the discretized defect equation \[ a(e_Q,v)=(f,v)-a(u_S,v) \quad \forall v\in Q \] satisfies \( e_Q = u_Q-u_S \) with \[ a(u_Q,v) = (f,v) \quad \in Q \] for \( u_Q \in Q \). The quantity \( || e_Q || \) is the first candidate for an a posteriori error estimate. If we select the space \( Q \subset H \) of higher order, e.g. piecewise quadratic, finite elements, \( Q \) can be regarded as hierarchical extension of \( S \).
Hence, for the usage of the hierarchic error estimator, an (probably higher order) extension space and the according assembler, variable sets, coefficient vectors and so on have to be created together with the error estimator itself:
using H1ExSpaces = boost::fusion::vector<H1Space<Grid> const*,H1HierarchicExtensionSpace<Grid> const*>;
using ExVariableDescriptions = boost::fusion::vector<VariableDescription<1,1,0> >;
using ExVSD = VariableSetDescription<H1ExSpaces,ExVariableDescriptions>;
using ErrorEstimator = HierarchicErrorEstimator<LinearizationAt<Functional>,ExVSD>;
using EstAssembler = VariationalFunctionalAssembler<ErrorEstimator>;
using ExCoefficientVectors = ExVSD::CoefficientVectorRepresentation<0,1>::type;
using AssEstOperator = AssembledGalerkinOperator<EstAssembler>;
During the error estimation step, the defect problem has to be solved and the error indicators have to be transferred to the cells. If the norm of the error is above a certain tolerance, the mesh will be refined. The procedure is as follows and can be viewed in detail in the code of example tutorials/HB_errorEstimation:
do{
[... SOLVE EQUATION ...]
dx.data = solution.data;
x += dx;
// hierarchical error estimation:
if (!tolX.empty()){
[... solve defect equation ...]
[... transfer error indicators to cells ...]
}
errNorm = sqrt(errNorm);
// evaluate error estimator information:
if (!tolX.empty())
{
if (errNorm<requested)
{
accurate = true ;
}
else
{
[... refine mesh ...]
}
}
} while (!accurate)
Mesh refinement is performed by marking the affected cells ci by gridManager.mark(1,ci) and calling gridManager.adaptAtOnce() to refine each marked leaf cell of the grid and prolongate FE functions to the new leaf grid view.
auto cend = variableSetDescription.gridView.end<0>();
for (auto ci=variableSetDescription.gridView.begin<0>(); ci!=cend; ++ci)
if (fabs(errorDistribution[is.index(*ci)]) >= errLevel) gridManager.mark(1,ci); // mark
accurate = !gridManager.adaptAtOnce(); // refine