From 012cfaa443573594a529445fc7deb54b4d066a52 Mon Sep 17 00:00:00 2001 From: Marc Fehling Date: Thu, 20 Jun 2019 19:45:07 +0200 Subject: [PATCH] hp::Refinement: Error prediction. --- include/deal.II/hp/refinement.h | 99 ++++++++++++++++++++++++ source/hp/refinement.cc | 69 +++++++++++++++++ source/hp/refinement.inst.in | 9 +++ tests/hp/error_prediction.cc | 124 +++++++++++++++++++++++++++++++ tests/hp/error_prediction.output | 19 +++++ 5 files changed, 320 insertions(+) create mode 100644 tests/hp/error_prediction.cc create mode 100644 tests/hp/error_prediction.output diff --git a/include/deal.II/hp/refinement.h b/include/deal.II/hp/refinement.h index f526d3527d..c125dd8e3e 100644 --- a/include/deal.II/hp/refinement.h +++ b/include/deal.II/hp/refinement.h @@ -280,6 +280,105 @@ namespace hp * @} */ + /** + * @name Error prediction + * @{ + */ + + /** + * Predict how the current @p error_indicators will adapt after refinement + * and coarsening has happened on the provided @p dof_handler, and write its + * results to @p predicted_errors. Each entry of @p error_indicators and + * @p predicted_errors corresponds to an active cell on the underlying + * Triangulation, thus each container has to be of size + * Triangulation::n_active_cells(). + * + * For h adaptation, we expect the local error $\eta_K$ on cell $K$ to be + * proportional to $(h_K)^{p_K}$ in the energy norm, where $h_K$ denotes the + * cell diameter and $p_K$ the polynomial degree of the currently assigned + * finite element. Here, we assume that the finite element will not change + * in the adaptation process so that $p_K = \text{const}$. However during + * coarsening, the finite elements on siblings may be different, and their + * parent cell will be assigned to their least dominating finite element + * that belongs to its most general child. Thus, we will always interpolate + * on an enclosing finite element space. Additionaly assuming that the + * finite elements on the cells to be coarsened are sufficient to represent + * the solution correct (e.g. at least quadratic basis functions for a + * quadratic solution), we are confident to say that the error will not + * change by sole interpolation on the larger finite element space. + * + * Further, we expect that the local error will be divided equally on + * all $2^{dim}$ children during refinement, whereas local errors on + * siblings will be summed up on the parent cell in case of coarsening. When + * transferring the predicted error to the coarsened mesh, make sure to + * configure your CellDataTransfer object with + * GridTools::CoarseningStrategies::sum() as a coarsening strategy. + * + * For p adaptation, the local error is expected to converge exponentially + * with the polynomial degree of the assigned finite element. Each increase + * or decrease of the degree will thus change its value by a user-defined + * control parameter @p gamma_p. The assumption of exponential convergence + * is only valid if both h and p adaptive methods are combined. An exception + * is thrown if a cell is flagged for both h and p adaptation at once. + * + * The prediction algorithm is formulated as follows with control parameters + * @p gamma_p, @p gamma_h and @p gamma_n that may be used to influence + * prediction for each adaptation type individually. + * + *
Adaptation type Prediction formula + *
no adaptation + * $\eta_{K,\text{pred}} = \eta_{K} \, \gamma_\text{n}$ + * $\gamma_\text{n} \in (0,\infty)$ + *
p adaptation + * $\eta_{K,\text{pred}} = \eta_{K} \, + * \gamma_\text{p}^{(p_{K,\text{future}} - p_K)}$ + * $\gamma_\text{p} \in (0,1)$ + *
h refinement + * $\eta_{K_c,\text{pred}} = \eta_{K} \, + * \gamma_\text{h} \, 0.5^{p_K} \, 0.5^{\text{dim}} + * \quad \forall K_c \text{ children of } K$ + * $\gamma_\text{h} \in (0,\infty)$ + *
h coarsening + * $\eta_{K,\text{pred}} = \sum\limits_{K_c} \eta_{K_c} / + * (\gamma_\text{h} \, 0.5^{p_{K_c}}) + * \quad \forall K_c \text{ children of } K$ + *
+ * + * For more theoretical details see + * @code{.bib} + * @article{Melenk2001, + * author = {Melenk, Jens Markus and Wohlmuth, Barbara I.}, + * title = {{On residual-based a posteriori error estimation + * in hp-FEM}}, + * journal = {{Advances in Computational Mathematics}}, + * volume = {15}, + * number = {1}, + * pages = {311--331}, + * publisher = {Springer US}, + * year = {2001}, + * doi = {10.1023/A:1014268310921} + * } + * @endcode + * + * @note This feature is currently only implemented for isotropic refinement. + * + * @note We want to predict the error by how adaptation will actually happen. + * Thus, this function needs to be called after + * Triangulation::prepare_for_coarsening_and_refinement(). + */ + template + void + predict_error(const hp::DoFHandler &dof_handler, + const Vector & error_indicators, + Vector & predicted_errors, + const double gamma_p = std::sqrt(0.1), + const double gamma_h = 1., + const double gamma_n = 1.); + + /** + * @} + */ + /** * @name Decide between h and p adaptivity * @{ diff --git a/source/hp/refinement.cc b/source/hp/refinement.cc index 82eb819177..97d1f44dfb 100644 --- a/source/hp/refinement.cc +++ b/source/hp/refinement.cc @@ -260,6 +260,75 @@ namespace hp + /** + * Error prediction + */ + template + void + predict_error(const hp::DoFHandler &dof_handler, + const Vector & error_indicators, + Vector & predicted_errors, + const double gamma_p, + const double gamma_h, + const double gamma_n) + { + AssertDimension(dof_handler.get_triangulation().n_active_cells(), + error_indicators.size()); + AssertDimension(dof_handler.get_triangulation().n_active_cells(), + predicted_errors.size()); + Assert(0 < gamma_p && gamma_p < 1, + dealii::GridRefinement::ExcInvalidParameterValue()); + Assert(0 < gamma_h, dealii::GridRefinement::ExcInvalidParameterValue()); + Assert(0 < gamma_n, dealii::GridRefinement::ExcInvalidParameterValue()); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + const unsigned int active_cell_index = cell->active_cell_index(); + + if (cell->future_fe_index_set()) // p adaptation + { + Assert(!cell->refine_flag_set() && !cell->coarsen_flag_set(), + ExcMessage("Cell has to be either flagged for h or p " + "adaptation, and not for both!")); + + const int degree_difference = + dof_handler.get_fe_collection()[cell->future_fe_index()] + .degree - + cell->get_fe().degree; + + predicted_errors[active_cell_index] = + error_indicators[active_cell_index] * + std::pow(gamma_p, degree_difference); + } + else if (cell->refine_flag_set()) // h refinement + { + Assert( + cell->refine_flag_set() == + RefinementCase::isotropic_refinement, + ExcMessage( + "Error prediction is only valid for isotropic refinement!")); + + predicted_errors[active_cell_index] = + error_indicators[active_cell_index] * + (gamma_h * std::pow(.5, dim + cell->get_fe().degree)); + } + else if (cell->coarsen_flag_set()) // h coarsening + { + predicted_errors[active_cell_index] = + error_indicators[active_cell_index] / + (gamma_h * std::pow(.5, cell->get_fe().degree)); + } + else // no changes + { + predicted_errors[active_cell_index] = + error_indicators[active_cell_index] * gamma_n; + } + } + } + + + /** * Decide between h and p adaptivity */ diff --git a/source/hp/refinement.inst.in b/source/hp/refinement.inst.in index 480c873bd6..2980563f27 100644 --- a/source/hp/refinement.inst.in +++ b/source/hp/refinement.inst.in @@ -74,6 +74,15 @@ for (deal_II_dimension : DIMENSIONS; S : REAL_SCALARS; const hp::DoFHandler &, const Vector &, const Vector &); + + template void + predict_error( + const hp::DoFHandler &, + const Vector &, + Vector &, + const double, + const double, + const double); \} \} #endif diff --git a/tests/hp/error_prediction.cc b/tests/hp/error_prediction.cc new file mode 100644 index 0000000000..f8837eecbc --- /dev/null +++ b/tests/hp/error_prediction.cc @@ -0,0 +1,124 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// validate error prediction algorithm for hp adaptive methods + + +#include + +#include +#include + +#include +#include +#include + +#include + +#include "../tests.h" + + + +template +void +test() +{ + const unsigned int n_cells = 4; + + // ----- setup ----- + Triangulation tria; + std::vector rep(dim, 1); + rep[0] = n_cells; + Point p1, p2; + for (unsigned int d = 0; d < dim; ++d) + { + p1[d] = 0; + p2[d] = (d == 0) ? n_cells : 1; + } + GridGenerator::subdivided_hyper_rectangle(tria, rep, p1, p2); + + hp::FECollection fes; + for (unsigned int d = 1; d <= 3; ++d) + fes.push_back(FE_Q(d)); + + hp::DoFHandler dh(tria); + dh.set_fe(fes); + for (const auto &cell : dh.active_cell_iterators()) + { + if (cell->id().to_string() == "0_0:") + cell->set_refine_flag(); + else if (cell->id().to_string() == "1_0:") + cell->set_coarsen_flag(); + else if (cell->id().to_string() == "2_0:") + cell->set_future_fe_index(2); + } + + // ----- predict ----- + Vector error_indicators, predicted_errors; + error_indicators.reinit(tria.n_active_cells()); + predicted_errors.reinit(tria.n_active_cells()); + for (unsigned int i = 0; i < tria.n_active_cells(); ++i) + error_indicators[i] = 10; + + hp::Refinement::predict_error(dh, + error_indicators, + predicted_errors, + /*gamma_p=*/0.5, + /*gamma_h=*/1., + /*gamma_n=*/1.); + + // ----- verify ------ + deallog << "ncells:" << tria.n_active_cells() << std::endl; + for (const auto &cell : dh.active_cell_iterators()) + { + deallog << " cell:" << cell->id().to_string() + << " fe_deg:" << cell->get_fe().degree + << " error:" << error_indicators[cell->active_cell_index()] + << " predict:" << predicted_errors[cell->active_cell_index()]; + + if (cell->refine_flag_set()) + deallog << " refining"; + else if (cell->coarsen_flag_set()) + deallog << " coarsening"; + else if (cell->future_fe_index_set()) + deallog << " future_fe_deg:" + << dh.get_fe_collection()[cell->future_fe_index()].degree; + + deallog << std::endl; + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char *argv[]) +{ + initlog(); + deallog << std::setprecision(3); + + deallog.push("1d"); + test<1>(); + deallog.pop(); + deallog.push("2d"); + test<2>(); + deallog.pop(); + deallog.push("3d"); + test<3>(); + deallog.pop(); +} diff --git a/tests/hp/error_prediction.output b/tests/hp/error_prediction.output new file mode 100644 index 0000000000..1d11f6b76b --- /dev/null +++ b/tests/hp/error_prediction.output @@ -0,0 +1,19 @@ + +DEAL:1d::ncells:4 +DEAL:1d:: cell:0_0: fe_deg:1 error:10.0 predict:2.50 refining +DEAL:1d:: cell:1_0: fe_deg:1 error:10.0 predict:20.0 coarsening +DEAL:1d:: cell:2_0: fe_deg:1 error:10.0 predict:2.50 future_fe_deg:3 +DEAL:1d:: cell:3_0: fe_deg:1 error:10.0 predict:10.0 +DEAL:1d::OK +DEAL:2d::ncells:4 +DEAL:2d:: cell:0_0: fe_deg:1 error:10.0 predict:1.25 refining +DEAL:2d:: cell:1_0: fe_deg:1 error:10.0 predict:20.0 coarsening +DEAL:2d:: cell:2_0: fe_deg:1 error:10.0 predict:2.50 future_fe_deg:3 +DEAL:2d:: cell:3_0: fe_deg:1 error:10.0 predict:10.0 +DEAL:2d::OK +DEAL:3d::ncells:4 +DEAL:3d:: cell:0_0: fe_deg:1 error:10.0 predict:0.625 refining +DEAL:3d:: cell:1_0: fe_deg:1 error:10.0 predict:20.0 coarsening +DEAL:3d:: cell:2_0: fe_deg:1 error:10.0 predict:2.50 future_fe_deg:3 +DEAL:3d:: cell:3_0: fe_deg:1 error:10.0 predict:10.0 +DEAL:3d::OK -- 2.39.5