]> https://gitweb.dealii.org/ - dealii.git/commitdiff
hp::Refinement: Error prediction. 8341/head
authorMarc Fehling <marc.fehling@gmx.net>
Thu, 20 Jun 2019 17:45:07 +0000 (19:45 +0200)
committerMarc Fehling <marc.fehling@gmx.net>
Mon, 8 Jul 2019 13:15:26 +0000 (15:15 +0200)
include/deal.II/hp/refinement.h
source/hp/refinement.cc
source/hp/refinement.inst.in
tests/hp/error_prediction.cc [new file with mode: 0644]
tests/hp/error_prediction.output [new file with mode: 0644]

index f526d3527d68e3fe0cce239a0de90e7cdbee4521..c125dd8e3ec7fbf282c86cca5696b7bd3b5b1ec3 100644 (file)
@@ -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.
+     * <table>
+     *   <tr><th>Adaptation type <th colspan="2">Prediction formula
+     *   <tr><td>no adaptation
+     *       <td>$\eta_{K,\text{pred}} = \eta_{K} \, \gamma_\text{n}$
+     *       <td>$\gamma_\text{n} \in (0,\infty)$
+     *   <tr><td>p adaptation
+     *       <td>$\eta_{K,\text{pred}} = \eta_{K} \,
+     *            \gamma_\text{p}^{(p_{K,\text{future}} - p_K)}$
+     *       <td>$\gamma_\text{p} \in (0,1)$
+     *   <tr><td>h refinement
+     *       <td>$\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$
+     *       <td rowspan="2">$\gamma_\text{h} \in (0,\infty)$
+     *   <tr><td>h coarsening
+     *       <td>$\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$
+     * </table>
+     *
+     * 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 <int dim, typename Number, int spacedim>
+    void
+    predict_error(const hp::DoFHandler<dim, spacedim> &dof_handler,
+                  const Vector<Number> &               error_indicators,
+                  Vector<Number> &                     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
      * @{
index 82eb819177c040721955d13aa4462919abb1e8f6..97d1f44dfb2392463376a5ca37097a5659b47ad2 100644 (file)
@@ -260,6 +260,75 @@ namespace hp
 
 
 
+    /**
+     * Error prediction
+     */
+    template <int dim, typename Number, int spacedim>
+    void
+    predict_error(const hp::DoFHandler<dim, spacedim> &dof_handler,
+                  const Vector<Number> &               error_indicators,
+                  Vector<Number> &                     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<dim>::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
      */
index 480c873bd6119ae0bf26ec605f0fd5a5be289ee7..2980563f27e04aa0a7359e5f006b429f82b40f99 100644 (file)
@@ -74,6 +74,15 @@ for (deal_II_dimension : DIMENSIONS; S : REAL_SCALARS;
           const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
           const Vector<S> &,
           const Vector<S> &);
+
+        template void
+        predict_error<deal_II_dimension, S, deal_II_space_dimension>(
+          const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
+          const Vector<S> &,
+          Vector<S> &,
+          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 (file)
index 0000000..f8837ee
--- /dev/null
@@ -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 <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/refinement.h>
+
+#include <deal.II/lac/vector.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+test()
+{
+  const unsigned int n_cells = 4;
+
+  // ----- setup -----
+  Triangulation<dim>        tria;
+  std::vector<unsigned int> rep(dim, 1);
+  rep[0] = n_cells;
+  Point<dim> 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<dim> fes;
+  for (unsigned int d = 1; d <= 3; ++d)
+    fes.push_back(FE_Q<dim>(d));
+
+  hp::DoFHandler<dim> 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<float> 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 (file)
index 0000000..1d11f6b
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.