]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Write in-code comments
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 5 Jun 2019 16:28:57 +0000 (18:28 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 19 Jun 2019 15:45:57 +0000 (17:45 +0200)
examples/step-65/doc/builds-on
examples/step-65/step-65.cc

index 76df1d0b00d1037801b0bba21bce8debe876351b..2302bdb7e5e0f6731ae70c5368475cf48deb7623 100644 (file)
@@ -1 +1 @@
-step-49
\ No newline at end of file
+step-49
index 44bdb5abdba6ce9fb4682c4ba2b01badc4a713e8..ce3e1ae8ee05c0f9e1db5debbb63e11ee3b7a63b 100644 (file)
@@ -1,19 +1,24 @@
-// ---------------------------------------------------------------------
-//
-// 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.
-//
-// ---------------------------------------------------------------------
-
-// This tutorial program was contributed by Martin Kronbichler
+/* ---------------------------------------------------------------------
+ *
+ * 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.
+ *
+ * ---------------------------------------------------------------------
+
+ * This tutorial program was contributed by Martin Kronbichler
+ */
+
+// @sect3{Include files}
+
+// The include files for this tutorial are essentially the same as in step-6.
 
 #include <deal.II/base/timer.h>
 
@@ -30,7 +35,6 @@
 
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/mapping_q_cache.h>
 #include <deal.II/fe/mapping_q_generic.h>
 
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/numerics/error_estimator.h>
 #include <deal.II/numerics/vector_tools.h>
 
+#include <fstream>
+
+// The only new include file is the one for the MappingQCache class.
+#include <deal.II/fe/mapping_q_cache.h>
 
 
 namespace step65
 {
   using namespace dealii;
 
+
+  // @sect3{Analytical solution and coefficient}
+
+  // In this tutorial program, we want to solve the Poisson equation with a
+  // coefficient that jumps along a sphere of radius 0.5 and constant right
+  // hand side of value $f(\mathbf{x}) = -3$. Due to the jump in the
+  // coefficient, the analytical solution must have a kink where the
+  // coefficient switches from one value to the other. To keep things simple,
+  // we select an analytical solution that is quadratic in all components,
+  // i.e., $u(x,y,z) = x^2 + y^2 + z^2$ in the ball of radius 0.5 and
+  // $u(x,y,z) = 0.1(x^2 + y^2 + z^2) + 0.25-0.025$ in the outer part of the
+  // domain. This analytical solution is compatible with the right hand side
+  // in case the coefficient is 0.5 in the inner ball and 5 outside.
   template <int dim>
   class ExactSolution : public Function<dim>
   {
@@ -87,6 +108,14 @@ namespace step65
 
 
 
+  // @sect3{The PoissonProblem class}
+  //
+  // The implementation of the Poisson problem is very similar to what we used
+  // in the step-5 tutorial program. The two main differences are that we pass
+  // a mapping object to the various steps in the program in order to switch
+  // between two mapping representations as explained in the introduction, and
+  // the time object that will be used for measuring the run times in the
+  // various cases.
   template <int dim>
   class PoissonProblem
   {
@@ -100,7 +129,6 @@ namespace step65
     void assemble_system(const Mapping<dim> &mapping);
     void solve();
     void postprocess(const Mapping<dim> &mapping);
-    void output_results(const Mapping<dim> &mapping);
 
     Triangulation<dim> triangulation;
     FE_Q<dim>          fe;
@@ -117,6 +145,10 @@ namespace step65
 
 
 
+  // In the constructor, we set up the timer object to record wall times but
+  // be quiet during the normal execution. We will query it for timing details
+  // in the PoissonProblem::run() function. Furthermore, we select a
+  // relatively high polynomial degree of three.
   template <int dim>
   PoissonProblem<dim>::PoissonProblem()
     : fe(3)
@@ -126,13 +158,57 @@ namespace step65
 
 
 
+  // @sect3{Grid creation and initialization of the manifolds}
+  //
+  // This function presents the typical usage of
+  // TransfiniteInterpolationManifold. The first step is to create the desired
+  // grid, which can be done by composition of two grids from
+  // GridGenerator. The inner ball mesh is simple enough: We run
+  // GridGenerator::hyper_cube centered in the origin with radius 0.5 (third
+  // function argument). The second mesh is more interesting and constructed
+  // as follows: We want to have a mesh that is spherical in the interior but
+  // flat on the outer surface. Furthermore, the mesh topology of the inner
+  // ball should be compatible with the outer grid in the sense that their
+  // vertices coincide to allow the two grid to be merged. The grid coming out
+  // of GridGenerator::hyper_shell fulfills the requirements on the inner side
+  // in case it is created with $2d$ coarse cells (6 coarse cells in 3D which
+  // we are going to use) &endash; this is the same number of cells as there
+  // are boundary faces for the ball. For the outer surface, we use the fact
+  // that the 6 faces on the surface of the shell without a manifold attached
+  // would degenerate to the surface of a cube. What we are still missing is
+  // the radius of the outer shell boundary. Since we desire a cube of extent
+  // $[-1, 1]$ and the 6-cell shell puts its 8 outer vertices at the 8
+  // opposing diagonals, we must translate the points $(\pm 1, \pm 1, \pm 1)$
+  // into a radius: Clearly, the radius must be $\sqrt{d}$ in $d$ dimensions,
+  // i.e., $\sqrt{3}$ for the three-dimensional case we want to consider.
+  //
+  // Thus, we have a plan: After creating the inner triangulation for the ball
+  // and the one for the outer shell, we merge those two grids but remove all
+  // manifolds from the resulting triangulation to ensure that we have full
+  // control over manifolds. In particular, we want additional points added on
+  // the boundary during refinement to follow a flat manifold description. To
+  // start the process of adding more appropriate manifold ids, we assign the
+  // manifold id 0 to all mesh entities (cells, faces, lines), which will
+  // later be associated with the TransfiniteInterpolationManifold. Then, we
+  // must identify the faces and lines that are along the sphere of radius 0.5
+  // and mark them with a different manifold id to then assign a
+  // SphericalManifold to those. We will choose the manifold id of 1. Since we
+  // have thrown away all manifolds that pre-existed after calling
+  // GridGenerator::hyper_ball(), we manually go through the cells of the mesh
+  // and all their faces. We have found a face on the sphere if all four
+  // vertices have a radius of 0.5, or, as we write in the program, have
+  // $r^2=0.25$. Note that we call `cell->face(f)->set_all_manifold_ids(1)` to
+  // set the manifold id both on the faces and the surrounding
+  // lines. Furthermore, we want to distinguish the cells inside the ball and
+  // outside the ball by a material id for visualization according to the
+  // picture in the introduction.
   template <int dim>
   void PoissonProblem<dim>::create_grid()
   {
     Triangulation<dim> tria_outer, tria_inner;
+    GridGenerator::hyper_ball(tria_inner, Point<dim>(), 0.5);
     GridGenerator::hyper_shell(
       tria_outer, Point<dim>(), 0.5, std::sqrt(dim), 2 * dim);
-    GridGenerator::hyper_ball(tria_inner, Point<dim>(), 0.5);
     GridGenerator::merge_triangulations(tria_inner, tria_outer, triangulation);
     triangulation.reset_all_manifolds();
     triangulation.set_all_manifold_ids(0);
@@ -144,9 +220,11 @@ namespace step65
             for (unsigned int v = 0;
                  v < GeometryInfo<dim - 1>::vertices_per_cell;
                  ++v)
-              if (std::abs(cell->face(f)->vertex(v).norm_square() - 0.25) >
-                  1e-12)
-                face_at_sphere_boundary = false;
+              {
+                if (std::abs(cell->face(f)->vertex(v).norm_square() - 0.25) >
+                    1e-12)
+                  face_at_sphere_boundary = false;
+              }
             if (face_at_sphere_boundary)
               cell->face(f)->set_all_manifold_ids(1);
           }
@@ -155,6 +233,21 @@ namespace step65
         else
           cell->set_material_id(0);
       }
+
+    // With all cells, faces and lines marked appropriately, we can attach the
+    // Manifold objects to those numbers. The entities with manifold id 1 will
+    // get a spherical manifold, whereas the other entities, which have the
+    // manifold id 0, will be assigned the
+    // TransfiniteInterpolationManifold. As mentioned in the introduction, we
+    // must explicitly initialize the manifold to the mesh by a call to
+    // TransfiniteInterpolationManifold::initialize() in order to pick up the
+    // coarse mesh cells and the manifolds attached to the boundaries of those
+    // cells. Note that the manifolds are allowed to go out of scope, because
+    // the Triangulation object internally calls Manifold::clone() to have a
+    // valid object around.
+    //
+    // With all manifolds attached, we will finally go about and refine the
+    // mesh a few times to create a sufficiently large test case.
     triangulation.set_manifold(1, SphericalManifold<dim>());
     TransfiniteInterpolationManifold<dim> transfinite_manifold;
     transfinite_manifold.initialize(triangulation);
@@ -165,6 +258,16 @@ namespace step65
 
 
 
+  // @sect3{Setup of data structures}
+  //
+  // This function is well-known from other tutorials in that it enumerates
+  // the degrees of freedom, creates a constraint object and sets up a sparse
+  // matrix for the linear system. The only thing worth mentioning is the fact
+  // that we pass the mapping to the
+  // VectorTools::interpolate_boundary_values() function to ensure that our
+  // boundary values are evaluated on the high-order mesh used for
+  // assembly. In the present example, it does not really matter because the
+  // outer surfaces are flat, but for curved outer cells this is mandatory.
   template <int dim>
   void PoissonProblem<dim>::setup_system(const Mapping<dim> &mapping)
   {
@@ -194,12 +297,62 @@ namespace step65
   }
 
 
-
+  // @sect3{Assembly of the system matrix and right hand side}
+  //
+  // This function is also well-known from the previous tutorial programs. One
+  // thing to note is that we set the number of quadrature points to the
+  // polynomial degree plus two, not the degree plus one as in most other
+  // tutorials. This is because we expect some extra accuracy as the mapping
+  // also involves a degree one more than the polynomials for the
+  // solution.
+  //
+  // The only somewhat unusual code in the assembly is the way we compute the
+  // cell matrix. Rather than using three nested loop over the quadrature
+  // point index, the row and column of the matrix, we first collect the
+  // derivatives of the shape function, multiplied by the square root of the
+  // product of the coefficient and the integration factor `JxW` in a separate
+  // matrix `partial_matrix`. To compute the cell matrix, we then execute
+  // `cell_matrix = partial_matrix * transpose(partial_matrix)` in the line
+  // `partial_matrix.mTmult(cell_matrix, partial_matrix);`. To understand why
+  // this works, we realize that the matrix-matrix multiplication performs a
+  // summation over the columns of `partial_matrix`. If we denote the
+  // coefficient by $a(\mathbf{x}_q)$, the entries in the temporary matrix are
+  // $\sqrt{\text{det}(J) w_q a(x)} \frac{\partial \varphi_i(\boldsymbol
+  // \xi_q)}{\partial x_k}$. If we take the product of the <i>i</i>th row with
+  // the <i>j</i>th column of that matrix, we compute a nested sum involving
+  // $\sum_q \sum_{k=1}^d \sqrt{\text{det}(J) w_q a(x)} \frac{\partial
+  // \varphi_i(\boldsymbol \xi_q)}{\partial x_k} \sqrt{\text{det}(J) w_q a(x)}
+  // \frac{\partial \varphi_j(\boldsymbol \xi_q)}{\partial x_k} = \sum_q
+  // \sum_{k=1}^d\text{det}(J) w_q a(x)\frac{\partial \varphi_i(\boldsymbol
+  // \xi_q)}{\partial x_k} \frac{\partial \varphi_j(\boldsymbol
+  // \xi_q)}{\partial x_k}$, which is exactly the terms needed for the
+  // Laplacian.
+  //
+  // The reason for choosing this somewhat unusual scheme is due to the heavy
+  // work involved in computing the cell matrix for a relatively high
+  // polynomial degree in 3D. As we want to highlight the cost of the mapping
+  // in this tutorial program, we better do the assembly in an optimized way
+  // in order to not chase bottlenecks that have been solved by the community
+  // already. Matrix-matrix multiplication is one of the best optimized
+  // kernels in the HPC context, and the FullMatrix::mTmult() function will
+  // call into those optimized BLAS functions. If the user has provided a good
+  // BLAS library when configuring deal.II (like OpenBLAS or Intel's MKL), the
+  // computation of the cell matrix will execute close to the processor's peak
+  // arithmetic performance. As a side note, we mention that despite an
+  // optimized matrix-matrix multiplication, the current strategy is
+  // sub-optimal in terms of complexity as the work to be done is proportional
+  // to $(p+1)^9$ operations for degree $p$ (this also applies to the usual
+  // evaluation with FEValues). One could compute the cell matrix with
+  // $\mathcal O((p+1)^7)$ operations by utilizing the tensor product
+  // structure of the shape functions, as is done by the matrix-free framework
+  // in deal.II. We refer to step-37 and the documentation of the
+  // tensor-product-aware evaluators FEEvaluation for details on how an even
+  // more efficient cell matrix computation could be realized.
   template <int dim>
   void PoissonProblem<dim>::assemble_system(const Mapping<dim> &mapping)
   {
     TimerOutput::Scope scope(timer, "Assemble linear system");
-    QGauss<dim>        quadrature_formula(fe.degree + 1);
+    QGauss<dim>        quadrature_formula(fe.degree + 2);
     FEValues<dim>      fe_values(mapping,
                             fe,
                             quadrature_formula,
@@ -213,6 +366,7 @@ namespace step65
     Vector<double>     cell_rhs(dofs_per_cell);
     FullMatrix<double> partial_matrix(dofs_per_cell, dim * n_q_points);
     std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
     for (const auto &cell : dof_handler.active_cell_iterators())
       {
         cell_rhs = 0.;
@@ -243,6 +397,10 @@ namespace step65
 
 
 
+  // @sect3{Solution of the linear system}
+  //
+  // For solving the linear system, we pick a simple Jacobi-preconditioned
+  // variant, similar to the settings in the early tutorials.
   template <int dim>
   void PoissonProblem<dim>::solve()
   {
@@ -259,20 +417,40 @@ namespace step65
 
 
 
+  // @sect3{Output of the solution and computation of errors}
+  //
+  // In this function, we do various post-processing steps with the
+  // solution, all of which involve the mapping in one way or the other.
+  //
+  // The first operation we do is to write the solution as well as the
+  // material ids to a vtu file. This is similar to what was done in many
+  // other tutorial programs. The new ingredient presented in this tutorial
+  // program is the way we ensure that a high-order representation is written
+  // to the file. We need to consider two particular topics. Firstly, we tell
+  // the DataOut object via the DataOutBase::VtkFlags that we intend to
+  // interpret the subdivisions of the elements as a high-order Lagrange
+  // polynomial. Recent visualization programs, like ParaView of version 5.5
+  // or newer, can then render a high-order solution (one typically needs to
+  // adjust a parameter called "nonlinear subdivision level"). Secondly, we
+  // need to make sure that the mapping is passed to the
+  // DataOut::build_patches() method. Finally, the DataOut class only prints
+  // curved boundary cells by default, so we need to ensure that even inner
+  // cells are printed in a curved representation via the mapping.
   template <int dim>
   void PoissonProblem<dim>::postprocess(const Mapping<dim> &mapping)
   {
     {
       TimerOutput::Scope scope(timer, "Write output");
-      Timer              time;
-      DataOut<dim>       data_out;
-      data_out.attach_dof_handler(dof_handler);
-      data_out.add_data_vector(solution, "solution");
+
+      DataOut<dim> data_out;
 
       DataOutBase::VtkFlags flags;
       flags.write_higher_order_cells = true;
       data_out.set_flags(flags);
 
+      data_out.attach_dof_handler(dof_handler);
+      data_out.add_data_vector(solution, "solution");
+
       Vector<double> material_ids(triangulation.n_active_cells());
       for (const auto &cell : triangulation.active_cell_iterators())
         material_ids[cell->active_cell_index()] = cell->material_id();
@@ -290,6 +468,16 @@ namespace step65
       data_out.write_vtu(file);
     }
 
+    // The next operation in the postprocessing function is to compute the L2
+    // and H1 errors against the analytical solution. As the analytical
+    // solution is a quadratic polynomial, we expect a very accurate result at
+    // this point. If we were solving on a simple mesh with affine element
+    // shapes, we would expect the numerical result to coincide with the
+    // analyical solution up to roundoff accuracy. However, since we are used
+    // deformed cells following a sphere, which are only tracked by
+    // polynomials of degree 4 (one more than the degree for the finite
+    // elements), we will make an error around 1e-7. We could get more
+    // accuracy by increasing the polynomial degree or refining the mesh.
     {
       TimerOutput::Scope scope(timer, "Compute error norms");
       Vector<double>     norm_per_cell_p(triangulation.n_active_cells());
@@ -315,6 +503,14 @@ namespace step65
                 << norm_per_cell_p.l2_norm() << std::endl;
     }
 
+    // The final post-processing operation we do here is to compute an error
+    // estimation with the KellyErrorEstimator. We use the exact same settings
+    // as in the step-6 tutorial program, except for the fact that we also
+    // hand in the mapping to ensure that errors are evaluated along the
+    // curved element, consistent with the remainder of the program. However,
+    // we do not really use the result here to drive a mesh adaptation step
+    // (that would refine the mesh around the material interface along the
+    // sphere), as the focus here is on the cost of this operation.
     {
       TimerOutput::Scope scope(timer, "Compute error estimator");
       Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
@@ -325,46 +521,38 @@ namespace step65
         std::map<types::boundary_id, const Function<dim> *>(),
         solution,
         estimated_error_per_cell);
+      std::cout << "   Max cell-wise error estimate: "
+                << estimated_error_per_cell.linfty_norm() << std::endl;
     }
   }
 
 
 
-  template <int dim>
-  void PoissonProblem<dim>::output_results(const Mapping<dim> &mapping)
-  {
-    TimerOutput::Scope scope(timer, "Write output");
-    Timer              time;
-    DataOut<dim>       data_out;
-    data_out.attach_dof_handler(dof_handler);
-    data_out.add_data_vector(solution, "solution");
-
-    DataOutBase::VtkFlags flags;
-    flags.write_higher_order_cells = true;
-    data_out.set_flags(flags);
-
-    Vector<double> material_ids(triangulation.n_active_cells());
-    for (const auto &cell : triangulation.active_cell_iterators())
-      material_ids[cell->active_cell_index()] = cell->material_id();
-    data_out.add_data_vector(material_ids, "material_ids");
-
-    data_out.build_patches(mapping,
-                           fe.degree,
-                           DataOut<dim>::curved_inner_cells);
-
-    std::ofstream file(
-      ("solution-" +
-       std::to_string(triangulation.n_global_levels() - 10 + 2 * dim) + ".vtu")
-        .c_str());
-    data_out.write_vtu(file);
-  }
-
-
-
+  // @sect3{The PoissonProblem::run() function}
+  //
+  // Finally, we define the `run()` function that controls how we want to
+  // execute this program (which is called by the main() function in the usual
+  // way). We start by calling the `create_grid()` function that sets up our
+  // geometry with the appropriate manifolds. We then run two instances of a
+  // solver chain, starting from the setup of the equations, the assembly of
+  // the linear system, its solution with a simple iterative solver, and the
+  // postprocessing discussed above. The two instances differ in the way they
+  // use the mapping. The first uses a conventional MappingQGeneric mapping
+  // object which we initialize to a degree one more than we use for the
+  // finite element &endash; after all, we expect the geometry representation
+  // to be the bottleneck as the analytic solution is only a quadratic
+  // polynomial. (In reality, things are interlinked to quite some extent
+  // because the evaluation of the polynomials in real coordinates involves
+  // the mapping of a higher-degree polynomials, which represent some smooth
+  // rational functions. As a consequence, higher-degree polynomials still pay
+  // off, so it does not make sense to increase the degree of the mapping
+  // further.) Once the first pass is completed, we let the timer print a
+  // summary of the compute times of the individual stages.
   template <int dim>
   void PoissonProblem<dim>::run()
   {
     create_grid();
+
     {
       std::cout << std::endl
                 << "====== Running with the basic MappingQGeneric class ====== "
@@ -376,13 +564,21 @@ namespace step65
       assemble_system(mapping);
       solve();
       postprocess(mapping);
-      output_results(mapping);
 
       timer.print_summary();
+      timer.reset();
     }
 
-    timer.reset();
-
+    // For the second instance, we instead set up the MappingQCache class. Its
+    // use is very simple: After constructing it (with the degree, given that
+    // we want it to show the correct degree functionality in other contexts),
+    // we fill the cache via the MappingQCache::initialize() function. At this
+    // stage, we specify which mapping we want to use (obviously, the same
+    // MappingQGeneric as previously in order to repeat the same computations)
+    // for the cache, and then run through the same functions again, now
+    // handing in the modified mapping. In the end, we again print the
+    // accumulated wall times since the reset to see how the times compare to
+    // the original setting.
     {
       std::cout
         << "====== Running with the optimized MappingQCache class ====== "
@@ -401,15 +597,14 @@ namespace step65
       assemble_system(mapping);
       solve();
       postprocess(mapping);
-      output_results(mapping);
 
       timer.print_summary();
     }
   }
-
 } // namespace step65
 
 
+
 int main()
 {
   step65::PoissonProblem<3> test_program;

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.