]> https://gitweb.dealii.org/ - dealii.git/commitdiff
update documentation, part II
authorMatthias Maier <tamiko@43-1.org>
Thu, 26 Dec 2019 15:31:40 +0000 (09:31 -0600)
committerMatthias Maier <tamiko@43-1.org>
Sat, 7 Mar 2020 18:28:54 +0000 (12:28 -0600)
examples/step-69/doc/intro.dox
examples/step-69/step-69.cc

index ec8b1890b690ff1dd8460c9e75d298c1d3d17fa3..1989403057705a5873dac72250baa18f94da8cf8 100644 (file)
@@ -3,8 +3,6 @@
   and Ignacio Tomas (Sandia National Laboratories, Albuquerque).
 </i>
 
-@dealiiTutorialDOI{10.5281/zenodo.TODO,https://zenodo.org/badge/DOI/10.5281/zenodo.TODO.svg}
-
 @note This tutorial step implements a first-order accurate <i>guaranteed
 maximum wavespeed method</i> based on a first-order <i>graph viscosity</i>
 for solving Euler's equations of gas dynamics @cite GuermondPopov2016. As
@@ -14,6 +12,8 @@ high-performance implementation of a second-order accurate scheme that uses
 <i>convex limiting</i> techniques, and strong stability-preserving (SSP)
 time integration, see @cite GuermondEtAl2018.
 
+@todo Add zenodo link
+
 
 <h1>Introduction</h1>
 
@@ -38,7 +38,9 @@ as elementary building blocks in higher-order schemes
 @cite GuermondEtAl2018. However, we hope that the reader still finds the
 tutorial step to be a good starting point (in particular with respect to
 the programming techniques) before jumping into full research codes such as
-the second-order scheme @cite GuermondEtAl2018 maintained TODO.
+the second-order scheme @cite GuermondEtAl2018.
+
+@todo Add link to repository and project of the second order code.
 
 
 <a name="eulerequations"></a>
@@ -173,8 +175,8 @@ instance @cite GuermondErn2004 Chapter 5 and references therein). Most
 time-dependent discretization approaches described in the deal.II tutorials
 are based on such a (semi-discrete) variational approach. Fundamentally,
 from an analysis perspective, variational discretizations are conceived in
-order to provide some notion of global (integral) stabiliy, meaning an estimate
-of the form
+order to provide some notion of global (integral) stabiliy, meaning an
+estimate of the form
 
 @f{align*}
   |\!|\!| \mathbf{u}_{h}(t) |\!|\!| \leq |\!|\!| \mathbf{u}_{h}(0) |\!|\!|
@@ -190,15 +192,17 @@ have proven to be some of the best approaches for simulations in the subsonic
 shockless regime and similarly benign regimes.
 
 However, in the transonic and supersonic regime, and shock-hydrodynamics
-applications the use of variational schemes might be questionable. In fact, at
-the time of this writing, most shock-hydrodynamics codes are still firmly
-grounded on finite volumes methods. The main reason for failure of variational
-schemes in such extreme regimes is the lack of pointwise stability. This stems
-from the fact that <i>a priori</i> bounds on integrated quantities (e.g.
-integrals of moments) have in general no implications on pointwise properties
-of the solution. While some of these problems might be alleviated by the
-(perpetual) chase of the right shock capturing scheme, finite difference-like
-and finite volume schemes still have an edge in many regards.
+
+applications the use of variational schemes might be questionable. In fact,
+at the time of this writing, most shock-hydrodynamics codes are still
+firmly grounded on finite volumes methods. The main reason for failure of
+variational schemes in such extreme regimes is the lack of pointwise
+stability. This stems from the fact that <i>a priori</i> bounds on
+integrated quantities (e.g. integrals of moments) have in general no
+implications on pointwise properties of the solution. While some of these
+problems might be alleviated by the (perpetual) chase of the right shock
+capturing scheme, finite difference-like and finite volume schemes still
+have an edge in many regards.
 
 In this tutorial step we therefore depart from variational schemes. We will
 present a completely algebraic formulation (with the flavor of a
@@ -209,13 +213,13 @@ collocation-type scheme) that preserves constraints pointwise, i.e.,
   \;\text{at every node}\;\mathbf{x}_i\;\text{of the mesh}.
 @f}
 
-However, contrary to finite difference/volume schemes, the
-scheme implemented in this step maximizes the use of finite element
-software infrastructure, works in any mesh, in any space dimension, and is
-theoretically guaranteed to always work, all the time, no exception. This
-illustrates that deal.ii can be used far beyond the context of variational
-schemes in Hilbert spaces and that a large number of classes, modules and
-namespaces from deal.ii can be adapted for such purpose.
+Contrary to finite difference/volume schemes, the scheme implemented in
+this step maximizes the use of finite element software infrastructure,
+works in any mesh, in any space dimension, and is theoretically guaranteed
+to always work, all the time, no exception. This illustrates that deal.ii
+can be used far beyond the context of variational schemes in Hilbert spaces
+and that a large number of classes, modules and namespaces from deal.ii can
+be adapted for such purpose.
 
 
 <h3>Description of the scheme </h3>
@@ -232,15 +236,14 @@ spaces $\pmb{\mathbb{V}}_h := \{\mathbb{V}_h\}^{d+2}$. Let $\mathbf{u}_h
 \in \mathcal{V}} \mathbf{U}_i \phi_i$ where $\mathbf{U}_i \in
 \mathbb{R}^{d+2}$ and $\phi_i$ is a scalar-valued shape function.
 
-<b>Note.</b>
-For simplicity we will consider the usual Lagrange finite elements. In such
-context $\{\mathbf{x}_i\}_{i \in \mathcal{V}}$ be the set of all "support
-points" (see @ref GlossSupport "this glossary entry") where $\mathbf{x}_i \in
-\mathbb{R}^d$. Then each integer index $i \in \mathcal{V}$
-uniquely identifies a support point $\mathbf{x}_i$ and/or scalar-valued shape
-function $\phi_i$.
+@note For simplicity we will consider the usual Lagrange finite elements.
+In such this context we let $\{\mathbf{x}_i\}_{i \in \mathcal{V}}$ denote
+the set of all support points (see @ref GlossSupport "this glossary
+entry"), where $\mathbf{x}_i \in \mathbb{R}^d$. Then each index $i \in
+\mathcal{V}$ uniquely identifies a support point $\mathbf{x}_i$, as well as
+a scalar-valued shape function $\phi_i$.
 
-With this notation we can define the scheme as
+With this notation at hand we can define the scheme as
 
 @f{align*}
   m_i \frac{\mathbf{U}_i^{n+1} - \mathbf{U}_i^{n}}{\tau}
@@ -266,44 +269,40 @@ Where
   - $\textbf{n}_{ij} = \frac{\mathbf{c}_{ij}}{ \|\mathbf{c}_{ij}\|_{\ell^2} }$
 
 The definition of $\lambda_{\text{max}} (\mathbf{U},\mathbf{V},
-\textbf{n})$ is far from trivial and we will postpone their definition in
-order to focus on the computational/coding issues of this tutorial Step.
-For the time being let's note that
-  - $m_i$ and $\mathbf{c}_{ij}$ do not evolve in time. It makes sense to
-    compute and store them once, and later recall them at very time step.
-    They are part of what we are going to call off-line data.
+\textbf{n})$ is far from trivial and we will postpone the precise
+definition in order to focus first on some algorithmic and implementational
+questions. We note that
+  - $m_i$ and $\mathbf{c}_{ij}$ do not evolve in time (provided we keep the
+    discretization fixed). It thus makes sense to assemble the matrices
+    once in a so called <i>offline computation</i> and reuse them in every
+    time step. They are part of what we are going to call off-line data.
   - At every time step we have to evaluate $\mathbb{f}(\mathbf{U}_j^{n})$ and
     $d_{ij} := \max \{ \lambda_{\text{max}}
     (\mathbf{U}_i^{n},\mathbf{U}_j^{n}, \textbf{n}_{ij}),
     \lambda_{\text{max}} (\mathbf{U}_j^{n}, \mathbf{U}_i^{n},
-    \textbf{n}_{ji}) \} \|\mathbf{c}_{ij}\|_{\ell^2} $
+    \textbf{n}_{ji}) \} \|\mathbf{c}_{ij}\|_{\ell^2} $, which will
+    constitute the bulk of the computational cost.
 
-<<<<<<< HEAD
-Before we start with the description of the implementation of this scheme, it
-is worth saying a thing or two about the "assembly" of this system. Consider
-for instance a hypothetical pseudo-code, illustrating
-a possible strategy to compute the solution $\textbf{U}^{n+1}$:
-=======
-Before we start with the description of the implementation of this scheme,
-it is worth saying a thing or two about the "assembly" of this system.
-Consider for instance a hypothetical pseudo-code, illustrating a possible
-strategy to compute the solution $\textbf{U}^{n+1}$:
->>>>>>> ee60914d6e... reindent, layout changes and address some review comments
+Consider the following pseudo-code, illustrating a possible straight
+forward strategy for computing the solution $\textbf{U}^{n+1}$ at a new
+time $t_{n+1} = t_n + \tau_n$ given a known state $\textbf{U}^{n}$ at time
+$t_n$:
 
 @f{align*}
-&\textbf{For } i \in \mathcal{V} \\
-&\ \ \ \  \{\mathbf{c}_{ij}\}_{j \in \mathcal{I}(i)} :=
-\texttt{gather_cij_vectors}(\textbf{c}, \mathcal{I}(i)) \\
-&\ \ \ \ \{\textbf{U}_j^n\}_{j \in \mathcal{I}(i)} :=
-\texttt{gather_state_vectors}(\textbf{U}^n, \mathcal{I}(i)) \\
-&\ \ \ \ \ \textbf{U}_i^{n+1} := \mathbf{U}_i^{n} \\
-&\ \ \ \ \textbf{For } j \in \mathcal{I}(i) \\
-&\ \ \ \ \ \ \ \  \textbf{U}_i^{n+1} := \textbf{U}_i^{n+1} - \frac{\tau}{m_i}
- \mathbb{f}(\mathbf{U}_j^{n})\cdot
-   \mathbf{c}_{ij} + d_{ij} \mathbf{U}_j^{n} \\
-&\ \ \ \ \textbf{EndFor} \\
-&\textbf{EndFor}
+&\textbf{for } i \in \mathcal{V} \\
+&\ \ \ \  \{\mathbf{c}_{ij}\}_{j \in \mathcal{I}(i)} \leftarrow \texttt{gather} (\textbf{c}, \mathcal{I}(i)) \\
+&\ \ \ \ \{\textbf{U}_j^n\}_{j \in \mathcal{I}(i)} \leftarrow \texttt{gather} (\textbf{U}^n, \mathcal{I}(i)) \\
+&\ \ \ \ \ \textbf{U}_i^{n+1} \leftarrow \mathbf{U}_i^{n} \\
+&\ \ \ \ \textbf{for } j \in \mathcal{I}(i) \\
+&\ \ \ \ \ \ \ \  \texttt{compute } d_{ij} \\
+&\ \ \ \ \ \ \ \  \texttt{compute } \mathbb{f}(\mathbf{U}_j^{n}) \\
+&\ \ \ \ \ \ \ \  \textbf{U}_i^{n+1} \leftarrow \textbf{U}_i^{n+1} - \frac{\tau_n}{m_i}
+ \mathbb{f}(\mathbf{U}_j^{n})\cdot \mathbf{c}_{ij} + d_{ij} \mathbf{U}_j^{n} \\
+&\ \ \ \ \textbf{end} \\
+&\ \ \ \ \texttt{scatter} (\textbf{U}^n, \mathcal{I}(i), \textbf{U}_i^n)) \\
+&\textbf{end}
 @f}
+<<<<<<< HEAD
 We note here that:
 - This "assembly" does not require any form of quadrature or cell-loops.
 - Here $\textbf{c}$ and $\textbf{U}^n$ are a global matrix and a global vector
@@ -329,22 +328,43 @@ of application of this kind of schemes, also called "edge-based" or
 "graph-based" finite element schemes (see for instance @cite Rainald2008 for
 more historical references).
 
-<<<<<<< HEAD
-This pseudo-code was introduced only to prepare the mindset of the reader for
-what is going to be presented in the in the next section. The
-actual implementation described in the next section is somewhat different from
-what is described in the pseudo-code but shares the same core mentality: we do
-not loop on cells but rather we loop on the edges of the sparsity graph (hence
-the name "edge-based" code) in order to assemble the system.
-=======
-This pseudo-code was introduced only to prepare the mindset of the reader
-for what is going to be presented in the in the next section. The actual
-implementation described in the next section is somewhat different from
-what is described in the pseudo-code but shares the same core mentality: we
-do not loop on cells but rather we loop on the edges of the sparsity graph
-(hence the name "edge-based" code) in order to assemble the system.
->>>>>>> ee60914d6e... reindent, layout changes and address some review comments
+We note that:
+- This algorithm does not require any form of quadrature or cell-loops.
+- Here, $\textbf{c}$ and $\textbf{U}^n$ are a global matrix and a global vector
+  containing all the vectors $\mathbf{c}_{ij}$ and all the states
+  $\mathbf{U}_j^n$, respectively.
+- $\texttt{gather}$ and $\texttt{scatter}$ are helper functions
+  that collect from, or distribute values from global vectors and matrices.
+- For an interior node $\mathbf{x}_i$ on a regular mesh in two space
+  dimensions (assuming a first-order polynomial space $\mathbb{Q}^1$) the
+  stencil $\mathcal{I}(i)$ contains nine entries. The update for a single
+  state $\textbf{U}_i^n$ thus depends on nine state-vectors
+  $\{\textbf{U}_j^n\}_{j \in \mathcal{I}(i)}$ (i.e., all the states in the
+  patch formed by the support of the shape function $\phi_i$). This is one
+  of the major differences compared to cell-based loop, where an update
+  typically only operates on states associated with a single cell.
+
+The actual implementation will deviate from above code in one key aspect:
+The time-step size $\tau$ has to be chosen subject to a CFL condition
+@f{align*}
+  \tau_n = c_{\text{cfl}}\,\min_{
+  i\in\mathcal{V}}\left(\frac{m_i}{-2\,d_{ii}^{n}}\right),
+@f}
+where $0<c_{\text{cfl}}\le1$ is a chosen constant. This will require to
+compute all $d_{ij}$ in a separte step prior to actually performing above
+update. The core principle remains unchanged, though: we do not loop over
+cells but rather over all edges of the sparsity graph.
+
+@note It is not uncommon to encounter such fully algebraic schemes (i.e.,
+no bilinear forms, no cell loops, and no quadrature) outside of the finite
+element community in the wieder CFD community. There is a rich history of
+application of this kind of schemes, also called <i>edge-based</i> or
+<i>graph-based</i> finite element schemes (see for instance @cite
+Rainald2008 for a historical overview).
 
+@todo Explain what to do for slip, dirichlet and do-nothing boundary
+conditions.
 
 <h3>Implementation of the scheme </h3>
 
+@todo Maybe comment on some key features of the implementation?!
index 8823c0f2cbda5974e704bc5d28696f9446d398ec..948fb2dc44c53bed76e904682e5f86ab0f0ff4be 100644 (file)
@@ -1,6 +1,6 @@
 /* ---------------------------------------------------------------------
  *
- * Copyright (C) 2012 - 2019 by the deal.II authors
+ * Copyright (C) 2019 - 2020 by the deal.II authors
  *
  * This file is part of the deal.II library.
  *
  */
 
 // @sect3{Include files}
-// The set of include files is quite standard. The most intriguing part
-// is that: either though this code is a "thread and mpi parallel"
-// we are using neither Trilinos nor PETSC vectors. Actually we are using dealii
-// distributed vectors <code>la_parallel_vector.h</code> and the regular dealii
-// sparse matrices <code>sparse_matrix.h</code>
+
+// The set of include files is quite standard. The most intriguing part is
+// the fact that we will rely solely on deal.II data structures for MPI
+// parallelization, in particular distributed::Triangulation and
+// LinearAlgebra::distributed::Vector included through
+// <code>distributed/tria.h</code> and
+// <code>lac/la_parallel_vector.h</code>. Instead of a Trilinos, or PETSc
+// specific matrix class, we will use a non-distributed
+// dealii::SparseMatrix (<code>lac/sparse_matrix.h</code>) to store the local
+// part of the $c_{ij}$, $n_{ij}$ and $d_{ij}$ matrices.
 #include <deal.II/base/conditional_ostream.h>
-#include <deal.II/base/graph_coloring.h>
 #include <deal.II/base/parallel.h>
 #include <deal.II/base/parameter_acceptor.h>
 #include <deal.II/base/partitioner.h>
 #include <deal.II/base/quadrature.h>
 #include <deal.II/base/timer.h>
 #include <deal.II/base/work_stream.h>
+
 #include <deal.II/distributed/tria.h>
+
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_renumbering.h>
 #include <deal.II/dofs/dof_tools.h>
+
 #include <deal.II/fe/fe.h>
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/mapping.h>
 #include <deal.II/fe/mapping_q.h>
+
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/manifold_lib.h>
+
 #include <deal.II/lac/dynamic_sparsity_pattern.h>
 #include <deal.II/lac/la_parallel_vector.h>
 #include <deal.II/lac/sparse_matrix.h>
 #include <deal.II/lac/sparse_matrix.templates.h>
 #include <deal.II/lac/vector.h>
+
 #include <deal.II/meshworker/scratch_data.h>
+
 #include <deal.II/numerics/data_out.h>
 #include <deal.II/numerics/vector_tools.h>
 
+// In addition to above deal.II specific includes, we also include four
+// boost headers. The first two are for binary archives that we will use
+// for implementing a checkpointing and restart mechanism.
 #include <boost/archive/binary_iarchive.hpp>
 #include <boost/archive/binary_oarchive.hpp>
-#include <boost/core/demangle.hpp>
+
+// The last two boost header files are for creating custom iterator ranges
+// over integer intervals.
 #include <boost/range/irange.hpp>
 #include <boost/range/iterator_range.hpp>
 
-#include <filesystem>
-
-
-// @sect3{Declaration/s of the namespace Step69}
+// @sect3{Class template declarations}
 
+// Next we declare all data structures and class templates of the example
+// step.
 namespace Step69
 {
-  enum Boundary : dealii::types::boundary_id
+  using namespace dealii;
+
+  // We start with an enum describing all possible boundary conditions
+  // encountered in this tutorial step. Such an enum allows us to refer to
+  // boundary types by a mnemonic (such as
+  // <code>Boundary::do_nothing</code>) rather than a numerical value.
+  enum Boundary : types::boundary_id
   {
     do_nothing = 0,
     slip       = 2,
     dirichlet  = 3,
   };
 
-
-  // @sect4{Declaration of <code>Discretization</code> class template}
-
-  // The main goal of this class is to digest the input file and act as a
-  // "container" of members that may be changed/decided at run time (through the
-  // input file). It was natural to derive this class from
-  // <code>dealii::ParameterAcceptor</code>. This class is in charge of
-  // initializing mpi comunicator, geometry dimensions, triangulation, mapping,
-  // finite element, mapping, and quadratures. If we think of the class
-  // <code>Discretization</code> as a "container": it doesn't contain any
-  // memmory demanding class member such a dof_handlers, vectors or matrices.
-  // The most memmory thirsty class member is the <code>
-  // dealii::parallel::distributed::Triangulation<dim></code>.
-
+  // @sect4{class <code>Discretization</code>}
+  //
+  // The class <code>Discretization</code> contains all data structures
+  // concerning the mesh (triangulation) and discretization (mapping,
+  // finite element, quadrature) of the problem. We use the
+  // ParameterAcceptor class to automatically populate problem-specific
+  // parameters, such as the geometry information
+  // (<code>immersed_disc_length</code>, etc.) or the refinement level
+  // (<code>refinement</code>) from a parameter file. This requires us to
+  // split the initialization of data structures into two functions: We
+  // initialize everything that does not depend on parameters in the
+  // constructor, and defer the creation of the mesh to the
+  // <code>setup()</code> method that can be called once all parameters are
+  // read-in via ParameterAcceptor::initialize().
+  //
   template <int dim>
-  class Discretization : public dealii::ParameterAcceptor
+  class Discretization : public ParameterAcceptor
   {
   public:
-    Discretization(const MPI_Comm &     mpi_communicator,
-                   dealii::TimerOutput &computing_timer,
-                   const std::string &  subsection = "Discretization");
+    Discretization(const MPI_Comm &   mpi_communicator,
+                   TimerOutput &      computing_timer,
+                   const std::string &subsection = "Discretization");
 
     void setup();
 
     const MPI_Comm &mpi_communicator;
 
-    dealii::parallel::distributed::Triangulation<dim> triangulation;
+    parallel::distributed::Triangulation<dim> triangulation;
 
-    const dealii::MappingQ<dim>   mapping;
-    const dealii::FE_Q<dim>       finite_element;
-    const dealii::QGauss<dim>     quadrature;
-    const dealii::QGauss<dim - 1> face_quadrature;
+    const MappingQ<dim>   mapping;
+    const FE_Q<dim>       finite_element;
+    const QGauss<dim>     quadrature;
+    const QGauss<dim - 1> face_quadrature;
 
   private:
-    dealii::TimerOutput &computing_timer;
+    TimerOutput &computing_timer;
 
     double immersed_disc_length;
     double immersed_disc_height;
@@ -115,87 +137,110 @@ namespace Step69
     unsigned int refinement;
   };
 
-  // @sect4{Declaration of <code>OfflineData</code> class template}
-
-  // The class OfflineData is initializes (and "owns")
-  // pretty much all the components of the discretization that
-  // do not evolve in time. In particular: dof_handler, sparsity
-  // patterns, boundary maps, lumped mass matrix, and other matrices
-  // and vectors that do not change in time are members of this class.
-  // The term "offline" here refers to the idea that all the class members
-  // of <code>OfflineData</code> are initialized and assigned values
-  // a "time step zero" and are not meant to be modified at any other later
-  // time step. For instance, the sparsity pattern should not
-  // change as we advance in time (we are not doing any form of adaptivity in
-  // space). Similarly, the entries of the vector
-  // <code>lumped_mass_matrix</code> should not be modified as we advance in
-  // time either.
+  // @sect4{class <code>OfflineData</code>}
   //
-  // Placeholder: Say something about BoundaryNormalMap.
+  // The class <code>OfflineData</code> contains pretty much all components
+  // of the discretization that do not evolve in time, in particular, the
+  // DoFHandler, SparsityPattern, boundary maps, the lumped mass, $c_ij$,
+  // and $n_ij$ matrices.
+  //
+  // Here, the term <i>offline</i> refers to the fact that all the class
+  // members of <code>OfflineData</code> have well-defined values
+  // independent of the current time step. This means that they can be
+  // initialized ahead of time (at <i>time step zero</i>) and are not meant
+  // to be modified at any other later time step. For instance, the
+  // sparsity pattern should not change as we advance in time (we are not
+  // doing any form of adaptivity in space). Similarly, the entries of the
+  // lumped mass matrix should not be modified as we advance in time
+  // either.
+  //
+  // We also compute and store a <code>boundary_normal_map</code> that
+  // contains a map from a global index of type `types:global_dof_index` of
+  // a boundary degree of freedom to a tuple consisting of a normal vector,
+  // the boundary id, and the position associated with the degree of
+  // freedom. We actually have to compuate and store this geometric
+  // information in this class because we won't have access to geometric
+  // (or cell-based) information later on in the algebraic loops over the
+  // sparsity pattern.
+  //
+  // @note Even though this class currently does not have any parameters
+  // that could be read in from a parameter file we nevertheless dervie
+  // from ParameterAcceptor and follow the same idiom of providing a
+  // <code>setup()</code> (and <code>assemble()</code>) method as for the
+  // class Discretization.
 
   template <int dim>
-  class OfflineData : public dealii::ParameterAcceptor
+  class OfflineData : public ParameterAcceptor
   {
   public:
     using BoundaryNormalMap =
-      std::map<dealii::types::global_dof_index,
-               std::tuple<dealii::Tensor<1, dim> /* normal */,
-                          dealii::types::boundary_id /* id */,
-                          dealii::Point<dim>> /* position */>;
+      std::map<types::global_dof_index,
+               std::tuple<Tensor<1, dim>, types::boundary_id, Point<dim>>>;
 
     OfflineData(const MPI_Comm &           mpi_communicator,
-                dealii::TimerOutput &      computing_timer,
+                TimerOutput &              computing_timer,
                 const Discretization<dim> &discretization,
                 const std::string &        subsection = "OfflineData");
 
     void setup();
     void assemble();
 
-    dealii::DoFHandler<dim> dof_handler;
+    DoFHandler<dim> dof_handler;
 
-    std::shared_ptr<const dealii::Utilities::MPI::Partitioner> partitioner;
+    std::shared_ptr<const Utilities::MPI::Partitioner> partitioner;
 
     unsigned int n_locally_owned;
     unsigned int n_locally_relevant;
 
-    dealii::SparsityPattern sparsity_pattern;
+    SparsityPattern sparsity_pattern;
 
     BoundaryNormalMap boundary_normal_map;
 
-    dealii::SparseMatrix<double>                  lumped_mass_matrix;
-    std::array<dealii::SparseMatrix<double>, dim> cij_matrix;
-    std::array<dealii::SparseMatrix<double>, dim> nij_matrix;
-    dealii::SparseMatrix<double>                  norm_matrix;
+    SparseMatrix<double>                  lumped_mass_matrix;
+    std::array<SparseMatrix<double>, dim> cij_matrix;
+    std::array<SparseMatrix<double>, dim> nij_matrix;
+    SparseMatrix<double>                  norm_matrix;
 
   private:
-    const MPI_Comm &     mpi_communicator;
-    dealii::TimerOutput &computing_timer;
+    const MPI_Comm &mpi_communicator;
+    TimerOutput &   computing_timer;
 
-    dealii::SmartPointer<const Discretization<dim>> discretization;
+    SmartPointer<const Discretization<dim>> discretization;
   };
 
-  // @sect4{Declaration of <code>ProblemDescription</code> class template}
-
-  // Most of the implementations of the members of this class will be utility
-  // classes/functions specific for Euler's equations:
-  // - The type alias <code>rank1_type</code> will be used for the states
-  // $\mathbf{U}_i^n$
-  // - The type alias <code>rank2_type</code> will be used for the fluxes
-  // $\mathbb{f}(\mathbf{U}_j^n)$.
-  // - The implementation of <code>momentum</code> will extract $\textbf{m}$
-  //   (out of the state vector $[\rho,\textbf{m},E]$) and store it in a
-  //   <code>Tensor<1, dim></code> for our convenience.
-  // - The implementation of <code>internal_energy</code> will compute
-  //   $E - \frac{|\textbf{m}|^2}{2\rho}$ from the state vector
+  // @sect4{class <code>ProblemDescription</code>}
+  //
+  // The member functions of this class are utility functions specific to
+  // Euler's equations:
+  // - The type alias <code>rank1_type</code> is used for the states
+  //   $\mathbf{U}_i^n$
+  // - The type alias <code>rank2_type</code> is used for the fluxes
+  //   $\mathbb{f}(\mathbf{U}_j^n)$.
+  // - The <code>momentum</code> function extracts $\textbf{m}$
+  //   out of the state vector $[\rho,\textbf{m},E]$) and stores it in a
+  //   <code>Tensor<1, dim></code>.
+  // - The <code>internal_energy</code> function computes $E -
+  //   \frac{|\textbf{m}|^2}{2\rho}$ from a given state vector
   //   $[\rho,\textbf{m},E]$.
   //
-  // The purpose of the remaining class members <code>component_names</code>,
-  // <code>pressure</code>, and <code>speed_of_sound</code>,
-  // is evident from their names. Most notably, the last
-  // one <code>compute_lambda_max</code> is in charge of computing
-  // $\lambda_{max}(\mathbf{U},\mathbf{V},\mathbf{n})$ which is required
-  // to compute the first order viscosity $d_{ij}$ as detailed in the section
-  // <b>Description of the scheme</b>.
+  // The purpose of the class members <code>component_names</code>,
+  // <code>pressure</code>, and <code>speed_of_sound</code>, is evident
+  // from their names. We also provide a function
+  // <code>compute_lambda_max</code>, that computes the wave speed estimate
+  // mentioned above, $\lambda_{max}(\mathbf{U},\mathbf{V},\mathbf{n})$,
+  // which is used in the computation of the $d_{ij}$ matrix.
+  //
+  // @note The <code>DEAL_II_ALWAYS_INLINE</code> macro expands to a
+  // (compiler specific) pragma that ensures that the corresponding
+  // function defined in this class is always inlined, i.e., the function
+  // body is put in place for every invocation of the function, and no call
+  // (and code indirection) is generated. This is stronger than the
+  // <code>inline</code> keyword, which is more or less a (mild) suggestion
+  // to the compiler that the programmer things it would be beneficial to
+  // inline the function. <code>DEAL_II_ALWAYS_INLINE</code> should only be
+  // used rarely and with caution in situations such as this one, where we
+  // actually know (due to benchmarking) that inlining the function in
+  // question actually improves performance.
 
   template <int dim>
   class ProblemDescription
@@ -203,15 +248,14 @@ namespace Step69
   public:
     static constexpr unsigned int problem_dimension = 2 + dim;
 
-    using rank1_type = dealii::Tensor<1, problem_dimension>;
-    using rank2_type =
-      dealii::Tensor<1, problem_dimension, dealii::Tensor<1, dim>>;
+    using rank1_type = Tensor<1, problem_dimension>;
+    using rank2_type = Tensor<1, problem_dimension, Tensor<1, dim>>;
 
     const static std::array<std::string, dim + 2> component_names;
 
     static constexpr double gamma = 7. / 5.;
 
-    static DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, dim>
+    static DEAL_II_ALWAYS_INLINE inline Tensor<1, dim>
     momentum(const rank1_type U);
 
     static DEAL_II_ALWAYS_INLINE inline double
@@ -225,34 +269,43 @@ namespace Step69
     static DEAL_II_ALWAYS_INLINE inline rank2_type f(const rank1_type U);
 
     static DEAL_II_ALWAYS_INLINE inline double
-    compute_lambda_max(const rank1_type              U_i,
-                       const rank1_type              U_j,
-                       const dealii::Tensor<1, dim> &n_ij);
+    compute_lambda_max(const rank1_type      U_i,
+                       const rank1_type      U_j,
+                       const Tensor<1, dim> &n_ij);
   };
 
-  // @sect4{Declaration of <code>InitialValues</code> class template}
-
-  // Placeholder here
+  // @sect4{class <code>InitialValues</code>}
+  //
+  // The class <code>InitialValues</code> only public data type is a
+  // function <code>initial_state</code> that computes the initial state of
+  // a given point and time. For the purpose of this example step we simply
+  // implement a homogeneous uniform flow field for which the direction and
+  // a 1D primitive state (density, velocity, pressure) are read from the
+  // parameter file.
+  //
+  // Instead of implementing yet another <code>setup()</code> function we
+  // use a callback function <code>parse_parameters_callback</code> that
+  // can be hooked up to corresponding signal of the ParameterAcceptor,
+  // <code>ParameterAcceptor::parse_parameters_call_back.connect(...)</code>.
 
   template <int dim>
-  class InitialValues : public dealii::ParameterAcceptor
+  class InitialValues : public ParameterAcceptor
   {
   public:
     using rank1_type = typename ProblemDescription<dim>::rank1_type;
 
     InitialValues(const std::string &subsection = "InitialValues");
 
-    std::function<rank1_type(const dealii::Point<dim> &point, double t)>
-      initial_state;
+    std::function<rank1_type(const Point<dim> &point, double t)> initial_state;
 
   private:
     void parse_parameters_callback();
 
-    dealii::Tensor<1, dim> initial_direction;
-    dealii::Tensor<1, 3>   initial_1d_state;
+    Tensor<1, dim> initial_direction;
+    Tensor<1, 3>   initial_1d_state;
   };
 
-  // @sect4{Declaration of <code>TimeStep</code> class template}
+  // @sect4{class <code>TimeStep</code>}
 
   // Placeholder here
 
@@ -393,14 +446,7 @@ namespace Step69
     vector_type output_vector;
   };
 
-} // namespace Step69
-
-
-// @sect3{Implementation of the classes in namespace <code>Step69</code>}
-
-namespace Step69
-{
-  using namespace dealii;
+  // @sect3{Implementation of the classes in namespace <code>Step69</code>}
 
   // @sect4{Implementation of the members of the class <code>Discretization</code>}
 

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.