]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Modernize step-51 8039/head
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 8 May 2019 15:08:01 +0000 (17:08 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 9 May 2019 07:42:39 +0000 (09:42 +0200)
1. Use range-based for loop for cell iterators
2. Move bodies of Function classes inline
3. Small changes to the documentation plus some cross-references

examples/step-51/doc/intro.dox
examples/step-51/doc/results.dox
examples/step-51/step-51.cc

index 26f510cb73f7c22c6fc72c28e5c95706f21bd4f6..1ab7b75c79e99a911ac5de340381957e9f5bd2bd 100644 (file)
@@ -20,7 +20,7 @@ each vertex <i>for each of the adjacent elements</i>, rather than just one,
 and similarly for edges and faces.  As an example of how fast the number of
 unknowns grows, consider the FE_DGPMonomial basis: each
 scalar solution component is represented by polynomials of degree $p$
-with $(1/dim!)*\prod_{i=1}^{dim}(p+i)$ degrees of freedom per
+with $(1/\text{dim}!) \prod_{i=1}^{\text{dim}}(p+i)$ degrees of freedom per
 element. Typically, all degrees of freedom in an element are coupled
 to all of the degrees of freedom in the adjacent elements.  The resulting
 discrete equations yield very large linear systems very quickly, especially
@@ -275,8 +275,9 @@ ingredients:
 <ol>
   <li> The computed solution gradient $\mathbf{q}_h$ converges at optimal rate,
    i.e., $\mathcal{O}(h^{p+1})$.
-  <li> The average of the scalar part of the solution, <i>u<sub>h</sub></i>,
-   on each cell $K$ super-converges at rate $\mathcal{O}(h^{p+2})$.
+  <li> The cell-wise average of the scalar part of the solution,
+   $\frac{(1,u_h)_K}{\text{vol}(K)}$, super-converges at rate
+   $\mathcal{O}(h^{p+2})$.
 </ol>
 
 We now introduce a new variable $u_h^* \in \mathcal{V}_h^{p+1}$, which we find
index d24090bfb870e4ccf2c4e172d9cd2d55ecaeec7e..6d88fd76a102172fce77d56a7b76faf0c9e4691f 100644 (file)
@@ -281,13 +281,15 @@ One final note on the efficiency comparison: We tried to use general-purpose
 sparse matrix structures and similar solvers (optimal AMG preconditioners for
 both without particular tuning of the AMG parameters on any of them) to give a
 fair picture of the cost versus accuracy of two methods, on a toy example. It
-should be noted however that geometric multigrid (GMG) for continuous finite elements is about a
-factor four to five faster for $p=3$ and $p=6$. The authors of this
-tutorial have not seen similarly advanced solvers for the HDG linear
-systems. Also, there are other implementation aspects for CG available such as
-fast matrix-free approaches as shown in step-37 that make higher order
-continuous elements more competitive. Again, it is not clear to the authors of
-the tutorial whether similar improvements could be made for HDG.
+should be noted however that geometric multigrid (GMG) for continuous finite
+elements is about a factor four to five faster for $p=3$ and $p=6$. As of
+2019, optimal-complexity iterative solvers for HDG are still under development
+in the research community. Also, there are other implementation aspects for CG
+available such as fast matrix-free approaches as shown in step-37 that make
+higher order continuous elements more competitive. Again, it is not clear to
+the authors of the tutorial whether similar improvements could be made for
+HDG. We refer to <a href="https://dx.doi.org/10.1137/16M110455X">Kronbichler
+and Wall (2018)</a> for a recent efficiency evaluation.
 
 
 <h3>Possibilities for improvements</h3>
@@ -370,13 +372,14 @@ improvements would make the most sense:
   problem size (the number of iterations for Q1 elements and global
   refinements starts at 35 for the small sizes but increase up to 701 for the
   largest size). To do better, one could for example use an algebraic
-  multigrid preconditioner from Trilinos. For diffusion-dominated
-  problems such as
-  the problem at hand with finer meshes, such a solver can be designed that
-  uses the matrix-vector products from the more efficient ChunkSparseMatrix on
-  the finest level, as long as we are not working in parallel with MPI. For
-  MPI-parallelized computation, a standard TrilinosWrappers::SparseMatrix can
-  be used.
+  multigrid preconditioner from Trilinos, or some more advanced variants as
+  the one discussed in <a
+  href="https://dx.doi.org/10.1137/16M110455X">Kronbichler and Wall
+  (2018)</a>. For diffusion-dominated problems such as the problem at hand
+  with finer meshes, such a solver can be designed that uses the matrix-vector
+  products from the more efficient ChunkSparseMatrix on the finest level, as
+  long as we are not working in parallel with MPI. For MPI-parallelized
+  computations, a standard TrilinosWrappers::SparseMatrix can be used.
 
   <li> Speed up assembly by pre-assembling parts that do not change from one
   cell to another (those that do neither contain variable coefficients nor
index 36a1d91145fa0c6cf969153191098fb582bf02e4..541247034b363f0d010e8d1ed11073d03df8fa81 100644 (file)
@@ -134,52 +134,40 @@ namespace Step51
       : Function<dim>()
     {}
 
-    virtual double value(const Point<dim> & p,
-                         const unsigned int component = 0) const override;
-
-    virtual Tensor<1, dim>
-    gradient(const Point<dim> & p,
-             const unsigned int component = 0) const override;
-  };
-
-
-
-  template <int dim>
-  double Solution<dim>::value(const Point<dim> &p, const unsigned int) const
-  {
-    double return_value = 0;
-    for (unsigned int i = 0; i < this->n_source_centers; ++i)
-      {
-        const Tensor<1, dim> x_minus_xi = p - this->source_centers[i];
-        return_value +=
-          std::exp(-x_minus_xi.norm_square() / (this->width * this->width));
-      }
-
-    return return_value / Utilities::fixed_power<dim>(
-                            std::sqrt(2. * numbers::PI) * this->width);
-  }
-
-
+    virtual double value(const Point<dim> &p,
+                         const unsigned int /*component*/ = 0) const override
+    {
+      double sum = 0;
+      for (unsigned int i = 0; i < this->n_source_centers; ++i)
+        {
+          const Tensor<1, dim> x_minus_xi = p - this->source_centers[i];
+          sum +=
+            std::exp(-x_minus_xi.norm_square() / (this->width * this->width));
+        }
 
-  template <int dim>
-  Tensor<1, dim> Solution<dim>::gradient(const Point<dim> &p,
-                                         const unsigned int) const
-  {
-    Tensor<1, dim> return_value;
+      return sum /
+             std::pow(2. * numbers::PI * this->width * this->width, dim / 2.);
+    }
 
-    for (unsigned int i = 0; i < this->n_source_centers; ++i)
-      {
-        const Tensor<1, dim> x_minus_xi = p - this->source_centers[i];
+    virtual Tensor<1, dim>
+    gradient(const Point<dim> &p,
+             const unsigned int /*component*/ = 0) const override
+    {
+      Tensor<1, dim> sum;
+      for (unsigned int i = 0; i < this->n_source_centers; ++i)
+        {
+          const Tensor<1, dim> x_minus_xi = p - this->source_centers[i];
 
-        return_value +=
-          (-2 / (this->width * this->width) *
-           std::exp(-x_minus_xi.norm_square() / (this->width * this->width)) *
-           x_minus_xi);
-      }
+          sum +=
+            (-2 / (this->width * this->width) *
+             std::exp(-x_minus_xi.norm_square() / (this->width * this->width)) *
+             x_minus_xi);
+        }
 
-    return return_value / Utilities::fixed_power<dim>(
-                            std::sqrt(2 * numbers::PI) * this->width);
-  }
+      return sum /
+             std::pow(2. * numbers::PI * this->width * this->width, dim / 2.);
+    }
+  };
 
 
 
@@ -196,21 +184,17 @@ namespace Step51
     {}
 
     virtual void vector_value(const Point<dim> &p,
-                              Vector<double> &  v) const override;
+                              Vector<double> &  v) const override
+    {
+      AssertDimension(v.size(), dim + 1);
+      Solution<dim>  solution;
+      Tensor<1, dim> grad = solution.gradient(p);
+      for (unsigned int d = 0; d < dim; ++d)
+        v[d] = -grad[d];
+      v[dim] = solution.value(p);
+    }
   };
 
-  template <int dim>
-  void SolutionAndGradient<dim>::vector_value(const Point<dim> &p,
-                                              Vector<double> &  v) const
-  {
-    AssertDimension(v.size(), dim + 1);
-    Solution<dim>  solution;
-    Tensor<1, dim> grad = solution.gradient(p);
-    for (unsigned int d = 0; d < dim; ++d)
-      v[d] = -grad[d];
-    v[dim] = solution.value(p);
-  }
-
 
 
   // Next comes the implementation of the convection velocity. As described in
@@ -224,42 +208,37 @@ namespace Step51
       : TensorFunction<1, dim>()
     {}
 
-    virtual Tensor<1, dim> value(const Point<dim> &p) const override;
+    virtual Tensor<1, dim> value(const Point<dim> &p) const override
+    {
+      Tensor<1, dim> convection;
+      switch (dim)
+        {
+          case 1:
+            convection[0] = 1;
+            break;
+          case 2:
+            convection[0] = p[1];
+            convection[1] = -p[0];
+            break;
+          case 3:
+            convection[0] = p[1];
+            convection[1] = -p[0];
+            convection[2] = 1;
+            break;
+          default:
+            Assert(false, ExcNotImplemented());
+        }
+      return convection;
+    }
   };
 
 
 
-  template <int dim>
-  Tensor<1, dim> ConvectionVelocity<dim>::value(const Point<dim> &p) const
-  {
-    Tensor<1, dim> convection;
-    switch (dim)
-      {
-        case 1:
-          convection[0] = 1;
-          break;
-        case 2:
-          convection[0] = p[1];
-          convection[1] = -p[0];
-          break;
-        case 3:
-          convection[0] = p[1];
-          convection[1] = -p[0];
-          convection[2] = 1;
-          break;
-        default:
-          Assert(false, ExcNotImplemented());
-      }
-    return convection;
-  }
-
-
-
-  // The last function we implement is the right hand side for the manufactured
-  // solution. It is very similar to step-7, with the exception that we now have
-  // a convection term instead of the reaction term. Since the velocity field is
-  // incompressible, i.e. $\nabla \cdot \mathbf{c} = 0$, this term simply reads
-  // $\mathbf{c} \nabla u$.
+  // The last function we implement is the right hand side for the
+  // manufactured solution. It is very similar to step-7, with the exception
+  // that we now have a convection term instead of the reaction term. Since
+  // the velocity field is incompressible, i.e., $\nabla \cdot \mathbf{c} =
+  // 0$, the advection term simply reads $\mathbf{c} \nabla u$.
   template <int dim>
   class RightHandSide : public Function<dim>, protected SolutionBase<dim>
   {
@@ -268,48 +247,43 @@ namespace Step51
       : Function<dim>()
     {}
 
-    virtual double value(const Point<dim> & p,
-                         const unsigned int component = 0) const override;
-
-  private:
-    const ConvectionVelocity<dim> convection_velocity;
-  };
+    virtual double value(const Point<dim> &p,
+                         const unsigned int /*component*/ = 0) const override
+    {
+      ConvectionVelocity<dim> convection_velocity;
+      Tensor<1, dim>          convection = convection_velocity.value(p);
+      double                  sum        = 0;
+      for (unsigned int i = 0; i < this->n_source_centers; ++i)
+        {
+          const Tensor<1, dim> x_minus_xi = p - this->source_centers[i];
 
+          sum +=
+            ((2 * dim - 2 * convection * x_minus_xi -
+              4 * x_minus_xi.norm_square() / (this->width * this->width)) /
+             (this->width * this->width) *
+             std::exp(-x_minus_xi.norm_square() / (this->width * this->width)));
+        }
 
-  template <int dim>
-  double RightHandSide<dim>::value(const Point<dim> &p,
-                                   const unsigned int) const
-  {
-    Tensor<1, dim> convection   = convection_velocity.value(p);
-    double         return_value = 0;
-    for (unsigned int i = 0; i < this->n_source_centers; ++i)
-      {
-        const Tensor<1, dim> x_minus_xi = p - this->source_centers[i];
+      return sum /
+             std::pow(2. * numbers::PI * this->width * this->width, dim / 2.);
+    }
+  };
 
-        return_value +=
-          ((2 * dim - 2 * convection * x_minus_xi -
-            4 * x_minus_xi.norm_square() / (this->width * this->width)) /
-           (this->width * this->width) *
-           std::exp(-x_minus_xi.norm_square() / (this->width * this->width)));
-      }
 
-    return return_value / Utilities::fixed_power<dim>(
-                            std::sqrt(2 * numbers::PI) * this->width);
-  }
 
   // @sect3{The HDG solver class}
 
   // The HDG solution procedure follows closely that of step-7. The major
-  // difference is the use of three different sets of <code>DoFHandler</code>
-  // and FE objects, along with the <code>ChunkSparseMatrix</code> and the
-  // corresponding solutions vectors. We also use WorkStream to enable a
-  // multithreaded local solution process which exploits the embarrassingly
-  // parallel nature of the local solver. For WorkStream, we define the local
-  // operations on a cell and a copy function into the global matrix and
-  // vector. We do this both for the assembly (which is run twice, once when we
-  // generate the system matrix and once when we compute the element-interior
-  // solutions from the skeleton values) and for the postprocessing where
-  // we extract a solution that converges at higher order.
+  // difference is the use of three different sets of DoFHandler and FE
+  // objects, along with the ChunkSparseMatrix and the corresponding solutions
+  // vectors. We also use WorkStream to enable a multithreaded local solution
+  // process which exploits the embarrassingly parallel nature of the local
+  // solver. For WorkStream, we define the local operations on a cell and a
+  // copy function into the global matrix and vector. We do this both for the
+  // assembly (which is run twice, once when we generate the system matrix and
+  // once when we compute the element-interior solutions from the skeleton
+  // values) and for the postprocessing where we extract a solution that
+  // converges at higher order.
   template <int dim>
   class HDG
   {
@@ -384,12 +358,19 @@ namespace Step51
 
     // The degrees of freedom corresponding to the skeleton strongly enforce
     // Dirichlet boundary conditions, just as in a continuous Galerkin finite
-    // element method.  We can enforce the boundary conditions in an analogous
-    // manner through the use of AffineConstraints constructs. In
-    // addition, hanging nodes are handled in the same way as for
-    // continuous finite elements: For the face elements which
-    // only define degrees of freedom on the face, this process sets the
-    // solution on the refined to be the one from the coarse side.
+    // element method. We can enforce the boundary conditions in an analogous
+    // manner via an AffineConstraints object. In addition, hanging nodes are
+    // handled in the same way as for continuous finite elements: For the face
+    // elements which only define degrees of freedom on the face, this process
+    // sets the solution on the refined side to coincide with the
+    // representation on the coarse side.
+    //
+    // Note that for HDG, the elimination of hanging nodes is not the only
+    // possibility &mdash; in terms of the HDG theory, one could also use the
+    // unknowns from the refined side and express the local solution on the
+    // coarse side through the trace values on the refined side. However, such
+    // a setup is not as easily implemented in terms of deal.II loops and not
+    // further analyzed.
     AffineConstraints<double> constraints;
 
     // The usage of the ChunkSparseMatrix class is similar to the usual sparse
@@ -407,11 +388,10 @@ namespace Step51
   // @sect3{The HDG class implementation}
 
   // @sect4{Constructor}
-  // The constructor is similar to those in other examples,
-  // with the exception of handling multiple <code>DoFHandler</code> and
-  // <code>FiniteElement</code> objects. Note that we create a system of finite
-  // elements for the local DG part, including the gradient/flux part and the
-  // scalar part.
+  // The constructor is similar to those in other examples, with the exception
+  // of handling multiple DoFHandler and FiniteElement objects. Note that we
+  // create a system of finite elements for the local DG part, including the
+  // gradient/flux part and the scalar part.
   template <int dim>
   HDG<dim>::HDG(const unsigned int degree, const RefinementMode refinement_mode)
     : fe_local(FE_DGQ<dim>(degree), dim, FE_DGQ<dim>(degree), 1)
@@ -428,7 +408,7 @@ namespace Step51
   // @sect4{HDG::setup_system}
   // The system for an HDG solution is setup in an analogous manner to most
   // of the other tutorial programs.  We are careful to distribute dofs with
-  // all of our <code>DoFHandler</code> objects.  The @p solution and @p system_matrix
+  // all of our DoFHandler objects.  The @p solution and @p system_matrix
   // objects go with the global skeleton solution.
   template <int dim>
   void HDG<dim>::setup_system()
@@ -458,7 +438,7 @@ namespace Step51
     constraints.close();
 
     // When creating the chunk sparsity pattern, we first create the usual
-    // compressed sparsity pattern and then set the chunk size, which is equal
+    // dynamic sparsity pattern and then set the chunk size, which is equal
     // to the number of dofs on a face, when copying this into the final
     // sparsity pattern.
     {
@@ -508,11 +488,11 @@ namespace Step51
 
   // @sect4{HDG::ScratchData}
   // @p ScratchData contains persistent data for each
-  // thread within <code>WorkStream</code>.  The <code>FEValues</code>, matrix,
+  // thread within WorkStream.  The FEValues, matrix,
   // and vector objects should be familiar by now.  There are two objects that
-  // need to be discussed: @p std::vector<std::vector<unsigned int> >
-  // fe_local_support_on_face and @p std::vector<std::vector<unsigned int> >
-  // fe_support_on_face.  These are used to indicate whether or not the finite
+  // need to be discussed: `std::vector<std::vector<unsigned int> >
+  // fe_local_support_on_face` and `std::vector<std::vector<unsigned int> >
+  // fe_support_on_face`.  These are used to indicate whether or not the finite
   // elements chosen have support (non-zero values) on a given face of the
   // reference cell for the local part associated to @p fe_local and the
   // skeleton part @p fe. We extract this information in the
@@ -621,7 +601,7 @@ namespace Step51
 
 
   // @sect4{HDG::PostProcessScratchData}
-  // @p PostProcessScratchData contains the data used by <code>WorkStream</code>
+  // @p PostProcessScratchData contains the data used by WorkStream
   // when post-processing the local solution $u^*$.  It is similar, but much
   // simpler, than @p ScratchData.
   template <int dim>
@@ -669,7 +649,7 @@ namespace Step51
 
 
   // @sect4{HDG::assemble_system}
-  // The @p assemble_system function is similar to <code>Step-32</code>, where
+  // The @p assemble_system function is similar to the one on Step-32, where
   // the quadrature formula and the update flags are set up, and then
   // <code>WorkStream</code> is used to do the work in a multi-threaded
   // manner.  The @p trace_reconstruct input parameter is used to decide
@@ -1211,8 +1191,8 @@ namespace Step51
   // codimension-1 surfaces
   // of the triangulation.  Our @p output_results function writes all local solutions
   // to the same vtk file, even though they correspond to different
-  // <code>DoFHandler</code> objects.  The graphical output for the skeleton
-  // variable is done through use of the <code>DataOutFaces</code> class.
+  // DoFHandler objects.  The graphical output for the skeleton
+  // variable is done through use of the DataOutFaces class.
   template <int dim>
   void HDG<dim>::output_results(const unsigned int cycle)
   {
@@ -1359,9 +1339,7 @@ namespace Step51
     // conditions. Since we re-create the triangulation every time for global
     // refinement, the flags are set in every refinement step, not just at the
     // beginning.
-    typename Triangulation<dim>::cell_iterator cell = triangulation.begin(),
-                                               endc = triangulation.end();
-    for (; cell != endc; ++cell)
+    for (const auto &cell : triangulation.cell_iterators())
       for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
            ++face)
         if (cell->face(face)->at_boundary())

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.