From b8c597b3189f9c4c2c2e062ba43035e9ca4cb8e2 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Wed, 5 Jun 2019 18:28:57 +0200 Subject: [PATCH] Write in-code comments --- examples/step-65/doc/builds-on | 2 +- examples/step-65/step-65.cc | 323 ++++++++++++++++++++++++++------- 2 files changed, 260 insertions(+), 65 deletions(-) diff --git a/examples/step-65/doc/builds-on b/examples/step-65/doc/builds-on index 76df1d0b00..2302bdb7e5 100644 --- a/examples/step-65/doc/builds-on +++ b/examples/step-65/doc/builds-on @@ -1 +1 @@ -step-49 \ No newline at end of file +step-49 diff --git a/examples/step-65/step-65.cc b/examples/step-65/step-65.cc index 44bdb5abdb..ce3e1ae8ee 100644 --- a/examples/step-65/step-65.cc +++ b/examples/step-65/step-65.cc @@ -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 @@ -30,7 +35,6 @@ #include #include -#include #include #include @@ -40,12 +44,29 @@ #include #include +#include + +// The only new include file is the one for the MappingQCache class. +#include 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 class ExactSolution : public Function { @@ -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 class PoissonProblem { @@ -100,7 +129,6 @@ namespace step65 void assemble_system(const Mapping &mapping); void solve(); void postprocess(const Mapping &mapping); - void output_results(const Mapping &mapping); Triangulation triangulation; FE_Q 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 PoissonProblem::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 void PoissonProblem::create_grid() { Triangulation tria_outer, tria_inner; + GridGenerator::hyper_ball(tria_inner, Point(), 0.5); GridGenerator::hyper_shell( tria_outer, Point(), 0.5, std::sqrt(dim), 2 * dim); - GridGenerator::hyper_ball(tria_inner, Point(), 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::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()); TransfiniteInterpolationManifold 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 void PoissonProblem::setup_system(const Mapping &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 ith row with + // the jth 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 void PoissonProblem::assemble_system(const Mapping &mapping) { TimerOutput::Scope scope(timer, "Assemble linear system"); - QGauss quadrature_formula(fe.degree + 1); + QGauss quadrature_formula(fe.degree + 2); FEValues fe_values(mapping, fe, quadrature_formula, @@ -213,6 +366,7 @@ namespace step65 Vector cell_rhs(dofs_per_cell); FullMatrix partial_matrix(dofs_per_cell, dim * n_q_points); std::vector 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 void PoissonProblem::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 void PoissonProblem::postprocess(const Mapping &mapping) { { TimerOutput::Scope scope(timer, "Write output"); - Timer time; - DataOut data_out; - data_out.attach_dof_handler(dof_handler); - data_out.add_data_vector(solution, "solution"); + + DataOut 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 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 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 estimated_error_per_cell(triangulation.n_active_cells()); @@ -325,46 +521,38 @@ namespace step65 std::map *>(), solution, estimated_error_per_cell); + std::cout << " Max cell-wise error estimate: " + << estimated_error_per_cell.linfty_norm() << std::endl; } } - template - void PoissonProblem::output_results(const Mapping &mapping) - { - TimerOutput::Scope scope(timer, "Write output"); - Timer time; - DataOut 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 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::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 void PoissonProblem::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; -- 2.39.5