step-69.h

Go to the documentation of this file.

 
 *             cell->get_dof_indices(dof_indices);
 *             std::transform(dof_indices.begin(),
 *                            dof_indices.end(),
 *                            dof_indices.begin(),
 *                            [&](types::global_dof_index index) {
 *                              return partitioner->global_to_local(index);
 *                            });
 *   
 *             /* And simply add, for each dof, a coupling to all other "local"
 *              * dofs on the cell: */
 *             for (const auto dof : dof_indices)
 *               dsp.add_entries(dof, dof_indices.begin(), dof_indices.end());
 *           }
 *   
 *         sparsity_pattern.copy_from(dsp);
 *   
 *         lumped_mass_matrix.reinit(sparsity_pattern);
 *         norm_matrix.reinit(sparsity_pattern);
 *         for (auto &matrix : cij_matrix)
 *           matrix.reinit(sparsity_pattern);
 *         for (auto &matrix : nij_matrix)
 *           matrix.reinit(sparsity_pattern);
 *       }
 *     }
 *   
 * @endcode
 * 
 * This concludes the setup of the DoFHandler and SparseMatrix objects.
 * Next, we have to assemble various matrices. We define a number of
 * helper functions and data structures in an anonymous namespace.
 * 
 
 * 
 * 
 * @code
 *     namespace
 *     {
 * @endcode
 * 
 * <code>CopyData</code> class that will be used to assemble the
 * offline data matrices using WorkStream. It acts as a container: it
 * is just a struct where WorkStream stores the local cell
 * contributions. Note that it also contains a class member
 * <code>local_boundary_normal_map</code> used to store the local
 * contributions required to compute the normals at the boundary.
 * 
 
 * 
 * 
 * @code
 *       template <int dim>
 *       struct CopyData
 *       {
 *         bool                                         is_artificial;
 *         std::vector<types::global_dof_index>         local_dof_indices;
 *         typename OfflineData<dim>::BoundaryNormalMap local_boundary_normal_map;
 *         FullMatrix<double>                           cell_lumped_mass_matrix;
 *         std::array<FullMatrix<double>, dim>          cell_cij_matrix;
 *       };
 *   
 * @endcode
 * 
 * Next we introduce a number of helper functions that are all
 * concerned about reading and writing matrix and vector entries. They
 * are mainly motivated by providing slightly more efficient code and
 * <a href="https://en.wikipedia.org/wiki/Syntactic_sugar"> syntactic
 * sugar</a> for otherwise somewhat tedious code.
 * 
 
 * 
 * The first function we introduce, <code>get_entry()</code>, will be
 * used to read the value stored at the entry pointed by a
 * SparsityPattern iterator <code>it</code> of <code>matrix</code>. The
 * function works around a small deficiency in the SparseMatrix
 * interface: The SparsityPattern is concerned with all index
 * operations of the sparse matrix stored in CRS format. As such the
 * iterator already knows the global index of the corresponding matrix
 * entry in the low-level vector stored in the SparseMatrix object. Due
 * to the lack of an interface in the SparseMatrix for accessing the
 * element directly with a SparsityPattern iterator, we unfortunately
 * have to create a temporary SparseMatrix iterator. We simply hide
 * this in the <code>get_entry()</code> function.
 * 
 
 * 
 * 
 * @code
 *       template <typename IteratorType>
 *       DEAL_II_ALWAYS_INLINE inline SparseMatrix<double>::value_type
 *       get_entry(const SparseMatrix<double> &matrix, const IteratorType &it)
 *       {
 *         const SparseMatrix<double>::const_iterator matrix_iterator(
 *           &matrix, it->global_index());
 *         return matrix_iterator->value();
 *       }
 *   
 * @endcode
 * 
 * The <code>set_entry()</code> helper is the inverse operation of
 * <code>get_value()</code>: Given an iterator and a value, it sets the
 * entry pointed to by the iterator in the matrix.
 * 
 
 * 
 * 
 * @code
 *       template <typename IteratorType>
 *       DEAL_II_ALWAYS_INLINE inline void
 *       set_entry(SparseMatrix<double> &           matrix,
 *                 const IteratorType &             it,
 *                 SparseMatrix<double>::value_type value)
 *       {
 *         SparseMatrix<double>::iterator matrix_iterator(&matrix,
 *                                                        it->global_index());
 *         matrix_iterator->value() = value;
 *       }
 *   
 * @endcode
 * 
 * <code>gather_get_entry()</code>: we note that @f$\mathbf{c}_{ij} \in
 * \mathbb{R}^d@f$. If @f$d=2@f$ then @f$\mathbf{c}_{ij} =
 * [\mathbf{c}_{ij}^1,\mathbf{c}_{ij}^2]^\top@f$. Which basically implies
 * that we need one matrix per space dimension to store the
 * @f$\mathbf{c}_{ij}@f$ vectors. Similar observation follows for the
 * matrix @f$\mathbf{n}_{ij}@f$. The purpose of
 * <code>gather_get_entry()</code> is to retrieve those entries and store
 * them into a <code>Tensor<1, dim></code> for our convenience.
 * 
 
 * 
 * 
 * @code
 *       template <std::size_t k, typename IteratorType>
 *       DEAL_II_ALWAYS_INLINE inline Tensor<1, k>
 *       gather_get_entry(const std::array<SparseMatrix<double>, k> &c_ij,
 *                        const IteratorType                         it)
 *       {
 *         Tensor<1, k> result;
 *         for (unsigned int j = 0; j < k; ++j)
 *           result[j] = get_entry(c_ij[j], it);
 *         return result;
 *       }
 *   
 * @endcode
 * 
 * <code>gather()</code> (first interface): this first function
 * signature, having three input arguments, will be used to retrieve
 * the individual components <code>(i,l)</code> of a matrix. The
 * functionality of <code>gather_get_entry()</code> and
 * <code>gather()</code> is very much the same, but their context is
 * different: the function <code>gather()</code> does not rely on an
 * iterator (that actually knows the value pointed to) but rather on the
 * indices <code>(i,l)</code> of the entry in order to retrieve its
 * actual value. We should expect <code>gather()</code> to be slightly
 * more expensive than <code>gather_get_entry()</code>. The use of
 * <code>gather()</code> will be limited to the task of computing the
 * algebraic viscosity @f$d_{ij}@f$ in the particular case that when
 * both @f$i@f$ and @f$j@f$ lie at the boundary.
 *     
 
 * 
 * @note The reader should be aware that accessing an arbitrary
 * <code>(i,l)</code> entry of a matrix (say for instance Trilinos or PETSc
 * matrices) is in general unacceptably expensive. Here is where we might
 * want to keep an eye on complexity: we want this operation to have
 * constant complexity, which is the case of the current implementation
 * using deal.II matrices.
 * 
 
 * 
 * 
 * @code
 *       template <std::size_t k>
 *       DEAL_II_ALWAYS_INLINE inline Tensor<1, k>
 *       gather(const std::array<SparseMatrix<double>, k> &n_ij,
 *              const unsigned int                         i,
 *              const unsigned int                         j)
 *       {
 *         Tensor<1, k> result;
 *         for (unsigned int l = 0; l < k; ++l)
 *           result[l] = n_ij[l](i, j);
 *         return result;
 *       }
 *   
 * @endcode
 * 
 * <code>gather()</code> (second interface): this second function
 * signature having two input arguments will be used to gather the
 * state at a node <code>i</code> and return it as a
 * <code>Tensor<1,n_solution_variables></code> for our convenience.
 * 
 
 * 
 * 
 * @code
 *       template <std::size_t k>
 *       DEAL_II_ALWAYS_INLINE inline Tensor<1, k>
 *       gather(const std::array<LinearAlgebra::distributed::Vector<double>, k> &U,
 *              const unsigned int                                               i)
 *       {
 *         Tensor<1, k> result;
 *         for (unsigned int j = 0; j < k; ++j)
 *           result[j] = U[j].local_element(i);
 *         return result;
 *       }
 *   
 * @endcode
 * 
 * <code>scatter()</code>: this function has three input arguments, the
 * first one is meant to be a "global object" (say a locally owned or
 * locally relevant vector), the second argument which could be a
 * <code>Tensor<1,n_solution_variables></code>, and the last argument
 * which represents a index of the global object. This function will be
 * primarily used to write the updated nodal values, stored as
 * <code>Tensor<1,n_solution_variables></code>, into the global objects.
 * 
 
 * 
 * 
 * @code
 *       template <std::size_t k, int k2>
 *       DEAL_II_ALWAYS_INLINE inline void
 *       scatter(std::array<LinearAlgebra::distributed::Vector<double>, k> &U,
 *               const Tensor<1, k2> &                                      tensor,
 *               const unsigned int                                         i)
 *       {
 *         static_assert(k == k2,
 *                       "The dimensions of the input arguments must agree");
 *         for (unsigned int j = 0; j < k; ++j)
 *           U[j].local_element(i) = tensor[j];
 *       }
 *     } // namespace
 *   
 * @endcode
 * 
 * We are now in a position to assemble all matrices stored in
 * <code>OfflineData</code>: the lumped mass entries @f$m_i@f$, the
 * vector-valued matrices @f$\mathbf{c}_{ij}@f$ and @f$\mathbf{n}_{ij} =
 * \frac{\mathbf{c}_{ij}}{|\mathbf{c}_{ij}|}@f$, and the boundary normals
 * @f$\boldsymbol{\nu}_i@f$.
 *   
 
 * 
 * In order to exploit thread parallelization we use the WorkStream approach
 * detailed in the @ref threads "Parallel computing with multiple processors"
 * accessing shared memory. As customary this requires
 * definition of
 * - Scratch data (i.e. input info required to carry out computations): in
 * this case it is <code>scratch_data</code>.
 * - The worker: in our case this is the <code>local_assemble_system()</code>
 * function that
 * actually computes the local (i.e. current cell) contributions from the
 * scratch data.
 * - A copy data: a struct that contains all the local assembly
 * contributions, in this case <code>CopyData<dim>()</code>.
 * - A copy data routine: in this case it is
 * <code>copy_local_to_global()</code> in charge of actually copying these
 * local contributions into the global objects (matrices and/or vectors)
 *   
 
 * 
 * Most of the following lines are spent in the definition of the worker
 * <code>local_assemble_system()</code> and the copy data routine
 * <code>copy_local_to_global()</code>. There is not much to say about the
 * WorkStream framework since the vast majority of ideas are reasonably
 * well-documented in @ref step_9 "step-9", @ref step_13 "step-13" and @ref step_32 "step-32" among others.
 *   
 
 * 
 * Finally, assuming that @f$\mathbf{x}_i@f$ is a support point at the boundary,
 * the (nodal) normals are defined as:
 *   
 
 * 
 * @f{align*}
 * \widehat{\boldsymbol{\nu}}_i \dealcoloneq
 * \frac{\int_{\partial\Omega} \phi_i \widehat{\boldsymbol{\nu}} \,
 * \, \mathrm{d}\mathbf{s}}{\big|\int_{\partial\Omega} \phi_i
 * \widehat{\boldsymbol{\nu}} \, \mathrm{d}\mathbf{s}\big|}
 * @f}
 *   
 
 * 
 * We will compute the numerator of this expression first and store it in
 * <code>OfflineData<dim>::BoundaryNormalMap</code>. We will normalize these
 * vectors in a posterior loop.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     void OfflineData<dim>::assemble()
 *     {
 *       lumped_mass_matrix = 0.;
 *       norm_matrix        = 0.;
 *       for (auto &matrix : cij_matrix)
 *         matrix = 0.;
 *       for (auto &matrix : nij_matrix)
 *         matrix = 0.;
 *   
 *       unsigned int dofs_per_cell =
 *         discretization->finite_element.n_dofs_per_cell();
 *       unsigned int n_q_points = discretization->quadrature.size();
 *   
 * @endcode
 * 
 * What follows is the initialization of the scratch data required by
 * WorkStream
 * 
 
 * 
 * 
 * @code
 *       MeshWorker::ScratchData<dim> scratch_data(
 *         discretization->mapping,
 *         discretization->finite_element,
 *         discretization->quadrature,
 *         update_values | update_gradients | update_quadrature_points |
 *           update_JxW_values,
 *         discretization->face_quadrature,
 *         update_normal_vectors | update_values | update_JxW_values);
 *   
 *       {
 *         TimerOutput::Scope scope(
 *           computing_timer,
 *           "offline_data - assemble lumped mass matrix, and c_ij");
 *   
 *         const auto local_assemble_system = 
 *           [&](const typename DoFHandler<dim>::cell_iterator &cell,
 *               MeshWorker::ScratchData<dim> &                 scratch,
 *               CopyData<dim> &                                copy) {
 *             copy.is_artificial = cell->is_artificial();
 *             if (copy.is_artificial)
 *               return;
 *   
 *             copy.local_boundary_normal_map.clear();
 *             copy.cell_lumped_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
 *             for (auto &matrix : copy.cell_cij_matrix)
 *               matrix.reinit(dofs_per_cell, dofs_per_cell);
 *   
 *             const auto &fe_values = scratch.reinit(cell);
 *   
 *             copy.local_dof_indices.resize(dofs_per_cell);
 *             cell->get_dof_indices(copy.local_dof_indices);
 *   
 *             std::transform(copy.local_dof_indices.begin(),
 *                            copy.local_dof_indices.end(),
 *                            copy.local_dof_indices.begin(),
 *                            [&](types::global_dof_index index) {
 *                              return partitioner->global_to_local(index);
 *                            });
 *   
 * @endcode
 * 
 * We compute the local contributions for the lumped mass matrix
 * entries @f$m_i@f$ and and vectors @f$c_{ij}@f$ in the usual fashion:
 * 
 * @code
 *             for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *               {
 *                 const auto JxW = fe_values.JxW(q_point);
 *   
 *                 for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                   {
 *                     const auto value_JxW =
 *                       fe_values.shape_value(j, q_point) * JxW;
 *                     const auto grad_JxW = fe_values.shape_grad(j, q_point) * JxW;
 *   
 *                     copy.cell_lumped_mass_matrix(j, j) += value_JxW;
 *   
 *                     for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *                       {
 *                         const auto value = fe_values.shape_value(i, q_point);
 *                         for (unsigned int d = 0; d < dim; ++d)
 *                           copy.cell_cij_matrix[d](i, j) += value * grad_JxW[d];
 *   
 *                       } /* i */
 *                   }     /* j */
 *               }         /* q */
 *   
 * @endcode
 * 
 * Now we have to compute the boundary normals. Note that the
 * following loop does not do much unless the element has faces on
 * the boundary of the domain.
 * 
 * @code
 *             for (const auto f : cell->face_indices())
 *               {
 *                 const auto face = cell->face(f);
 *                 const auto id   = face->boundary_id();
 *   
 *                 if (!face->at_boundary())
 *                   continue;
 *   
 *                 const auto &fe_face_values = scratch.reinit(cell, f);
 *   
 *                 const unsigned int n_face_q_points =
 *                   fe_face_values.get_quadrature().size();
 *   
 *                 for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                   {
 *                     if (!discretization->finite_element.has_support_on_face(j, f))
 *                       continue;
 *   
 * @endcode
 * 
 * Note that "normal" will only represent the contributions
 * from one of the faces in the support of the shape
 * function phi_j. So we cannot normalize this local
 * contribution right here, we have to take it "as is",
 * store it and pass it to the copy data routine. The
 * proper normalization requires an additional loop on
 * nodes. This is done in the copy function below.
 * 
 * @code
 *                     Tensor<1, dim> normal;
 *                     if (id == Boundaries::free_slip)
 *                       {
 *                         for (unsigned int q = 0; q < n_face_q_points; ++q)
 *                           normal += fe_face_values.normal_vector(q) *
 *                                     fe_face_values.shape_value(j, q);
 *                       }
 *   
 *                     const auto index = copy.local_dof_indices[j];
 *   
 *                     Point<dim> position;
 *                     for (const auto v : cell->vertex_indices())
 *                       if (cell->vertex_dof_index(v, 0) ==
 *                           partitioner->local_to_global(index))
 *                         {
 *                           position = cell->vertex(v);
 *                           break;
 *                         }
 *   
 *                     const auto old_id =
 *                       std::get<1>(copy.local_boundary_normal_map[index]);
 *                     copy.local_boundary_normal_map[index] =
 *                       std::make_tuple(normal, std::max(old_id, id), position);
 *                   }
 *               }
 *           };
 *   
 * @endcode
 * 
 * Last, we provide a copy_local_to_global function as required for
 * the WorkStream
 * 
 * @code
 *         const auto copy_local_to_global = [&](const CopyData<dim> &copy) {
 *           if (copy.is_artificial)
 *             return;
 *   
 *           for (const auto &it : copy.local_boundary_normal_map)
 *             {
 *               std::get<0>(boundary_normal_map[it.first]) +=
 *                 std::get<0>(it.second);
 *               std::get<1>(boundary_normal_map[it.first]) =
 *                 std::max(std::get<1>(boundary_normal_map[it.first]),
 *                          std::get<1>(it.second));
 *               std::get<2>(boundary_normal_map[it.first]) = std::get<2>(it.second);
 *             }
 *   
 *           lumped_mass_matrix.add(copy.local_dof_indices,
 *                                  copy.cell_lumped_mass_matrix);
 *   
 *           for (int k = 0; k < dim; ++k)
 *             {
 *               cij_matrix[k].add(copy.local_dof_indices, copy.cell_cij_matrix[k]);
 *               nij_matrix[k].add(copy.local_dof_indices, copy.cell_cij_matrix[k]);
 *             }
 *         };
 *   
 *         WorkStream::run(dof_handler.begin_active(),
 *                         dof_handler.end(),
 *                         local_assemble_system,
 *                         copy_local_to_global,
 *                         scratch_data,
 *                         CopyData<dim>());
 *       }
 *   
 * @endcode
 * 
 * At this point in time we are done with the computation of @f$m_i@f$ and
 * @f$\mathbf{c}_{ij}@f$, but so far the matrix <code>nij_matrix</code>
 * contains just a copy of the matrix <code>cij_matrix</code>.
 * That's not what we really want: we have to normalize its entries. In
 * addition, we have not filled the entries of the matrix
 * <code>norm_matrix</code>  and the vectors stored in the map
 * <code>OfflineData<dim>::BoundaryNormalMap</code> are not normalized.
 *     
 
 * 
 * In principle, this is just offline data, it doesn't make much sense
 * to over-optimize their computation, since their cost will get amortized
 * over the many time steps that we are going to use. However,
 * computing/storing the entries of the matrix
 * <code>norm_matrix</code> and the normalization of <code>nij_matrix</code>
 * are perfect to illustrate thread-parallel node-loops:
 * - we want to visit every node @f$i@f$ in the mesh/sparsity graph,
 * - and for every such node we want to visit to every @f$j@f$ such that
 * @f$\mathbf{c}_{ij} \not \equiv 0@f$.
 *     
 
 * 
 * From an algebraic point of view, this is equivalent to: visiting
 * every row in the matrix and for each one of these rows execute a loop on
 * the columns. Node-loops is a core theme of this tutorial step (see
 * the pseudo-code in the introduction) that will repeat over and over
 * again. That's why this is the right time to introduce them.
 *     
 
 * 
 * We have the thread parallelization capability
 * parallel::apply_to_subranges() that is somehow more general than the
 * WorkStream framework. In particular, parallel::apply_to_subranges() can
 * be used for our node-loops. This functionality requires four input
 * arguments which we explain in detail (for the specific case of our
 * thread-parallel node loops):
 * - The iterator <code>indices.begin()</code> points to a row index.
 * - The iterator <code>indices.end()</code> points to a numerically higher
 * row index.
 * - The function <code>on_subranges(index_begin, index_end)</code>
 * (where <code>index_begin</code> and <code>index_end</code>
 * define a sub-range within the range spanned by
 * the begin and end iterators defined in the two previous bullets)
 * applies an operation to every iterator in such subrange. We may as
 * well call <code>on_subranges</code> the "worker".
 * - Grainsize: minimum number of iterators (in this case representing
 * rows) processed by each thread. We decided for a minimum of 4096
 * rows.
 *     
 
 * 
 * A minor caveat here is that the iterators <code>indices.begin()</code>
 * and <code>indices.end()</code> supplied to
 * parallel::apply_to_subranges() have to be random access iterators:
 * internally, parallel::apply_to_subranges() will break the range
 * defined by the <code>indices.begin()</code> and
 * <code>indices.end()</code> iterators into subranges (we want to be
 * able to read any entry in those subranges with constant complexity).
 * In order to provide such iterators we resort to
 * std_cxx20::ranges::iota_view.
 *     
 
 * 
 * The bulk of the following piece of code is spent defining
 * the "worker" <code>on_subranges</code>: i.e. the  operation applied at
 * each row of the sub-range. Given a fixed <code>row_index</code>
 * we want to visit every column/entry in such row. In order to execute
 * such columns-loops we use
 * <a href="http://www.cplusplus.com/reference/algorithm/for_each/">
 * std::for_each</a>
 * from the standard library, where:
 * - <code>sparsity_pattern.begin(row_index)</code>
 * gives us an iterator starting at the first column of the row,
 * - <code>sparsity_pattern.end(row_index)</code> is an iterator pointing
 * at the last column of the row,
 * - the last argument required by `std::for_each` is the operation
 * applied at each nonzero entry (a lambda expression in this case)
 * of such row.
 *     
 
 * 
 * We note that, parallel::apply_to_subranges() will operate on
 * disjoint sets of rows (the subranges) and our goal is to write into
 * these rows. Because of the simple nature of the operations we want
 * to carry out (computation and storage of normals, and normalization
 * of the @f$\mathbf{c}_{ij}@f$ of entries) threads cannot conflict
 * attempting to write the same entry (we do not need a scheduler).
 * 
 * @code
 *       {
 *         TimerOutput::Scope scope(computing_timer,
 *                                  "offline_data - compute |c_ij|, and n_ij");
 *   
 *         const std_cxx20::ranges::iota_view<unsigned int, unsigned int> indices(
 *           0, n_locally_relevant);
 *   
 *         const auto on_subranges = 
 *           [&](const auto index_begin, const auto index_end) {
 *             for (const auto row_index :
 *                  std_cxx20::ranges::iota_view<unsigned int, unsigned int>(
 *                    *index_begin, *index_end))
 *               {
 * @endcode
 * 
 * First column-loop: we compute and store the entries of the
 * matrix norm_matrix and write normalized entries into the
 * matrix nij_matrix:
 * 
 * @code
 *                 std::for_each(sparsity_pattern.begin(row_index),
 *                               sparsity_pattern.end(row_index),
 *                               [&](const SparsityPatternIterators::Accessor &jt) {
 *                                 const auto c_ij =
 *                                   gather_get_entry(cij_matrix, &jt);
 *                                 const double norm = c_ij.norm();
 *   
 *                                 set_entry(norm_matrix, &jt, norm);
 *                                 for (unsigned int j = 0; j < dim; ++j)
 *                                   set_entry(nij_matrix[j], &jt, c_ij[j] / norm);
 *                               });
 *               }
 *           };
 *   
 *         parallel::apply_to_subranges(indices.begin(),
 *                                      indices.end(),
 *                                      on_subranges,
 *                                      4096);
 *   
 * @endcode
 * 
 * Finally, we normalize the vectors stored in
 * <code>OfflineData<dim>::BoundaryNormalMap</code>. This operation has
 * not been thread parallelized as it would neither illustrate any
 * important concept nor lead to any noticeable speed gain.
 * 
 * @code
 *         for (auto &it : boundary_normal_map)
 *           {
 *             auto &normal = std::get<0>(it.second);
 *             normal /= (normal.norm() + std::numeric_limits<double>::epsilon());
 *           }
 *       }
 *     }
 *   
 * @endcode
 * 
 * At this point we are very much done with anything related to offline data.
 * 
 
 * 
 * 
 * <a name="EquationofstateandapproximateRiemannsolver"></a> 
 * <h4>Equation of state and approximate Riemann solver</h4>
 * 
 
 * 
 * In this section we describe the implementation of the class members of
 * the <code>ProblemDescription</code> class. Most of the code here is
 * specific to the compressible Euler's equations with an ideal gas law.
 * If we wanted to re-purpose @ref step_69 "step-69" for a different conservation law
 * (say for: instance the shallow water equation) most of the
 * implementation of this class would have to change. But most of the
 * other classes (in particular those defining loop structures) would
 * remain unchanged.
 *   
 
 * 
 * We start by implementing a number of small member functions for
 * computing <code>momentum</code>, <code>internal_energy</code>,
 * <code>pressure</code>, <code>speed_of_sound</code>, and the flux
 * <code>f</code> of the system. The functionality of each one of these
 * functions is self-explanatory from their names.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     DEAL_II_ALWAYS_INLINE inline Tensor<1, dim>
 *     ProblemDescription<dim>::momentum(const state_type &U)
 *     {
 *       Tensor<1, dim> result;
 *       std::copy_n(&U[1], dim, &result[0]);
 *       return result;
 *     }
 *   
 *     template <int dim>
 *     DEAL_II_ALWAYS_INLINE inline double
 *     ProblemDescription<dim>::internal_energy(const state_type &U)
 *     {
 *       const double &rho = U[0];
 *       const auto    m   = momentum(U);
 *       const double &E   = U[dim + 1];
 *       return E - 0.5 * m.norm_square() / rho;
 *     }
 *   
 *     template <int dim>
 *     DEAL_II_ALWAYS_INLINE inline double
 *     ProblemDescription<dim>::pressure(const state_type &U)
 *     {
 *       return (gamma - 1.) * internal_energy(U);
 *     }
 *   
 *     template <int dim>
 *     DEAL_II_ALWAYS_INLINE inline double
 *     ProblemDescription<dim>::speed_of_sound(const state_type &U)
 *     {
 *       const double &rho = U[0];
 *       const double  p   = pressure(U);
 *   
 *       return std::sqrt(gamma * p / rho);
 *     }
 *   
 *     template <int dim>
 *     DEAL_II_ALWAYS_INLINE inline typename ProblemDescription<dim>::flux_type
 *     ProblemDescription<dim>::flux(const state_type &U)
 *     {
 *       const double &rho = U[0];
 *       const auto    m   = momentum(U);
 *       const auto    p   = pressure(U);
 *       const double &E   = U[dim + 1];
 *   
 *       flux_type result;
 *   
 *       result[0] = m;
 *       for (unsigned int i = 0; i < dim; ++i)
 *         {
 *           result[1 + i] = m * m[i] / rho;
 *           result[1 + i][i] += p;
 *         }
 *       result[dim + 1] = m / rho * (E + p);
 *   
 *       return result;
 *     }
 *   
 * @endcode
 * 
 * Now we discuss the computation of @f$\lambda_{\text{max}}
 * (\mathbf{U}_i^{n},\mathbf{U}_j^{n}, \textbf{n}_{ij})@f$. The analysis
 * and derivation of sharp upper-bounds of maximum wavespeeds of Riemann
 * problems is a very technical endeavor and we cannot include an
 * advanced discussion about it in this tutorial. In this portion of the
 * documentation we will limit ourselves to sketch the main functionality
 * of our implementation functions and point to specific academic
 * references in order to help the (interested) reader trace the
 * source (and proper mathematical justification) of these ideas.
 *   
 
 * 
 * In general, obtaining a sharp guaranteed upper-bound on the maximum
 * wavespeed requires solving a quite expensive scalar nonlinear problem.
 * This is typically done with an iterative solver. In order to simplify
 * the presentation in this example step we decided not to include such
 * an iterative scheme. Instead, we will just use an initial guess as a
 * guess for an upper bound on the maximum wavespeed. More precisely,
 * equations (2.11) (3.7), (3.8) and (4.3) of @cite GuermondPopov2016b
 * are enough to define a guaranteed upper bound on the maximum
 * wavespeed. This estimate is returned by a call to the function
 * <code>lambda_max_two_rarefaction()</code>. At its core the
 * construction of such an upper bound uses the so-called two-rarefaction
 * approximation for the intermediate pressure @f$p^*@f$, see for instance
 * Equation (4.46), page 128 in @cite Toro2009.
 *   
 
 * 
 * The estimate returned by <code>lambda_max_two_rarefaction()</code> is
 * guaranteed to be an upper bound, it is in general quite sharp, and
 * overall sufficient for our purposes. However, for some specific
 * situations (in particular when one of states is close to vacuum
 * conditions) such an estimate will be overly pessimistic. That's why we
 * used a second estimate to avoid this degeneracy that will be invoked
 * by a call to the function <code>lambda_max_expansion()</code>. The most
 * important function here is <code>compute_lambda_max()</code> which
 * takes the minimum between the estimates returned by
 * <code>lambda_max_two_rarefaction()</code> and
 * <code>lambda_max_expansion()</code>.
 *   
 
 * 
 * We start again by defining a couple of helper functions:
 *   
 
 * 
 * The first function takes a state <code>U</code> and a unit vector
 * <code>n_ij</code> and computes the <i>projected</i> 1d state in
 * direction of the unit vector.
 * 
 * @code
 *     namespace
 *     {
 *       template <int dim>
 *       DEAL_II_ALWAYS_INLINE inline std::array<double, 4> riemann_data_from_state(
 *         const typename ProblemDescription<dim>::state_type U,
 *         const Tensor<1, dim> &                             n_ij)
 *       {
 *         Tensor<1, 3> projected_U;
 *         projected_U[0] = U[0];
 *   
 * @endcode
 * 
 * For this, we have to change the momentum to @f$\textbf{m}\cdot
 * n_{ij}@f$ and have to subtract the kinetic energy of the
 * perpendicular part from the total energy:
 * 
 * @code
 *         const auto m   = ProblemDescription<dim>::momentum(U);
 *         projected_U[1] = n_ij * m;
 *   
 *         const auto perpendicular_m = m - projected_U[1] * n_ij;
 *         projected_U[2] = U[1 + dim] - 0.5 * perpendicular_m.norm_square() / U[0];
 *   
 * @endcode
 * 
 * We return the 1d state in <i>primitive</i> variables instead of
 * conserved quantities. The return array consists of density @f$\rho@f$,
 * velocity @f$u@f$, pressure @f$p@f$ and local speed of sound @f$a@f$:
 * 
 
 * 
 * 
 * @code
 *         return {{projected_U[0],
 *                  projected_U[1] / projected_U[0],
 *                  ProblemDescription<1>::pressure(projected_U),
 *                  ProblemDescription<1>::speed_of_sound(projected_U)}};
 *       }
 *   
 * @endcode
 * 
 * At this point we also define two small functions that return the
 * positive and negative part of a double.
 * 
 
 * 
 * 
 * @code
 *       DEAL_II_ALWAYS_INLINE inline double positive_part(const double number)
 *       {
 *         return std::max(number, 0.);
 *       }
 *   
 *   
 *       DEAL_II_ALWAYS_INLINE inline double negative_part(const double number)
 *       {
 *         return -std::min(number, 0.);
 *       }
 *   
 * @endcode
 * 
 * Next, we need two local wavenumbers that are defined in terms of a
 * primitive state @f$[\rho, u, p, a]@f$ and a given pressure @f$p^\ast@f$
 * @cite GuermondPopov2016  Eqn. (3.7):
 * @f{align*}
 * \lambda^- = u - a\,\sqrt{1 + \frac{\gamma+1}{2\gamma}
 * \left(\frac{p^\ast-p}{p}\right)_+}
 * @f}
 * Here, the @f$(\cdot)_{+}@f$ denotes the positive part of the given
 * argument.
 * 
 
 * 
 * 
 * @code
 *       DEAL_II_ALWAYS_INLINE inline double
 *       lambda1_minus(const std::array<double, 4> &riemann_data,
 *                     const double                 p_star)
 *       {
 *         /* Implements formula (3.7) in Guermond-Popov-2016 */
 *   
 *         constexpr double gamma = ProblemDescription<1>::gamma;
 *         const auto       u     = riemann_data[1];
 *         const auto       p     = riemann_data[2];
 *         const auto       a     = riemann_data[3];
 *   
 *         const double factor = (gamma + 1.0) / 2.0 / gamma;
 *         const double tmp    = positive_part((p_star - p) / p);
 *         return u - a * std::sqrt(1.0 + factor * tmp);
 *       }
 *   
 * @endcode
 * 
 * Analougously @cite GuermondPopov2016 Eqn. (3.8):
 * @f{align*}
 * \lambda^+ = u + a\,\sqrt{1 + \frac{\gamma+1}{2\gamma}
 * \left(\frac{p^\ast-p}{p}\right)_+}
 * @f}
 * 
 
 * 
 * 
 * @code
 *       DEAL_II_ALWAYS_INLINE inline double
 *       lambda3_plus(const std::array<double, 4> &riemann_data, const double p_star)
 *       {
 *         /* Implements formula (3.8) in Guermond-Popov-2016 */
 *   
 *         constexpr double gamma = ProblemDescription<1>::gamma;
 *         const auto       u     = riemann_data[1];
 *         const auto       p     = riemann_data[2];
 *         const auto       a     = riemann_data[3];
 *   
 *         const double factor = (gamma + 1.0) / 2.0 / gamma;
 *         const double tmp    = positive_part((p_star - p) / p);
 *         return u + a * std::sqrt(1.0 + factor * tmp);
 *       }
 *   
 * @endcode
 * 
 * All that is left to do is to compute the maximum of @f$\lambda^-@f$ and
 * @f$\lambda^+@f$ computed from the left and right primitive state
 * (@cite GuermondPopov2016 Eqn. (2.11)), where @f$p^\ast@f$ is given by
 * @cite GuermondPopov2016 Eqn (4.3):
 * 
 
 * 
 * 
 * @code
 *       DEAL_II_ALWAYS_INLINE inline double
 *       lambda_max_two_rarefaction(const std::array<double, 4> &riemann_data_i,
 *                                  const std::array<double, 4> &riemann_data_j)
 *       {
 *         constexpr double gamma = ProblemDescription<1>::gamma;
 *         const auto       u_i   = riemann_data_i[1];
 *         const auto       p_i   = riemann_data_i[2];
 *         const auto       a_i   = riemann_data_i[3];
 *         const auto       u_j   = riemann_data_j[1];
 *         const auto       p_j   = riemann_data_j[2];
 *         const auto       a_j   = riemann_data_j[3];
 *   
 *         const double numerator = a_i + a_j - (gamma - 1.) / 2. * (u_j - u_i);
 *   
 *         const double denominator =
 *           a_i * std::pow(p_i / p_j, -1. * (gamma - 1.) / 2. / gamma) + a_j * 1.;
 *   
 *         /* Formula (4.3) in Guermond-Popov-2016 */
 *   
 *         const double p_star =
 *           p_j * std::pow(numerator / denominator, 2. * gamma / (gamma - 1));
 *   
 *         const double lambda1 = lambda1_minus(riemann_data_i, p_star);
 *         const double lambda3 = lambda3_plus(riemann_data_j, p_star);
 *   
 *         /* Formula (2.11) in Guermond-Popov-2016 */
 *   
 *         return std::max(positive_part(lambda3), negative_part(lambda1));
 *       }
 *   
 * @endcode
 * 
 * We compute the second upper bound of the maximal wavespeed that is, in
 * general, not as sharp as the two-rarefaction estimate. But it will
 * save the day in the context of near vacuum conditions when the
 * two-rarefaction approximation might attain extreme values:
 * @f{align*}
 * \lambda_{\text{exp}} = \max(u_i,u_j) + 5. \max(a_i, a_j).
 * @f}
 * @note The constant 5.0 multiplying the maximum of the sound speeds
 * is <i>neither</i> an ad-hoc constant, <i>nor</i> a tuning parameter.
 * It defines an upper bound for any @f$\gamma \in [0,5/3]@f$. Do not play
 * with it!
 * 
 
 * 
 * 
 * @code
 *       DEAL_II_ALWAYS_INLINE inline double
 *       lambda_max_expansion(const std::array<double, 4> &riemann_data_i,
 *                            const std::array<double, 4> &riemann_data_j)
 *       {
 *         const auto u_i = riemann_data_i[1];
 *         const auto a_i = riemann_data_i[3];
 *         const auto u_j = riemann_data_j[1];
 *         const auto a_j = riemann_data_j[3];
 *   
 *         return std::max(std::abs(u_i), std::abs(u_j)) + 5. * std::max(a_i, a_j);
 *       }
 *     } // namespace
 *   
 * @endcode
 * 
 * The following is the main function that we are going to call in order to
 * compute @f$\lambda_{\text{max}} (\mathbf{U}_i^{n},\mathbf{U}_j^{n},
 * \textbf{n}_{ij})@f$. We simply compute both maximal wavespeed estimates
 * and return the minimum.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     DEAL_II_ALWAYS_INLINE inline double
 *     ProblemDescription<dim>::compute_lambda_max(const state_type &    U_i,
 *                                                 const state_type &    U_j,
 *                                                 const Tensor<1, dim> &n_ij)
 *     {
 *       const auto riemann_data_i = riemann_data_from_state(U_i, n_ij);
 *       const auto riemann_data_j = riemann_data_from_state(U_j, n_ij);
 *   
 *       const double lambda_1 =
 *         lambda_max_two_rarefaction(riemann_data_i, riemann_data_j);
 *   
 *       const double lambda_2 =
 *         lambda_max_expansion(riemann_data_i, riemann_data_j);
 *   
 *       return std::min(lambda_1, lambda_2);
 *     }
 *   
 * @endcode
 * 
 * We conclude this section by defining static arrays
 * <code>component_names</code> that contain strings describing the
 * component names of our state vector. We have template specializations
 * for dimensions one, two and three, that are used later in DataOut for
 * naming the corresponding components:
 * 
 
 * 
 * 
 * @code
 *     template <>
 *     const std::array<std::string, 3> ProblemDescription<1>::component_names{
 *       {"rho", "m", "E"}};
 *   
 *     template <>
 *     const std::array<std::string, 4> ProblemDescription<2>::component_names{
 *       {"rho", "m_1", "m_2", "E"}};
 *   
 *     template <>
 *     const std::array<std::string, 5> ProblemDescription<3>::component_names{
 *       {"rho", "m_1", "m_2", "m_3", "E"}};
 *   
 * @endcode
 * 
 * 
 * <a name="Initialvalues"></a> 
 * <h4>Initial values</h4>
 * 
 
 * 
 * The last preparatory step, before we discuss the implementation of the
 * forward Euler scheme, is to briefly implement the `InitialValues` class.
 *   
 
 * 
 * In the constructor we initialize all parameters with default values,
 * declare all parameters for the `ParameterAcceptor` class and connect the
 * <code>parse_parameters_call_back</code> slot to the respective signal.
 *   
 
 * 
 * The <code>parse_parameters_call_back</code> slot will be invoked from
 * ParameterAceptor after the call to ParameterAcceptor::initialize(). In
 * that regard, its use is appropriate for situations where the
 * parameters have to be postprocessed (in some sense) or some
 * consistency condition between the parameters has to be checked.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     InitialValues<dim>::InitialValues(const std::string &subsection)
 *       : ParameterAcceptor(subsection)
 *     {
 *       /* We wire up the slot InitialValues<dim>::parse_parameters_callback to
 *          the ParameterAcceptor::parse_parameters_call_back signal: */
 *       ParameterAcceptor::parse_parameters_call_back.connect(
 *         [&]() { this->parse_parameters_callback(); });
 *   
 *       initial_direction[0] = 1.;
 *       add_parameter("initial direction",
 *                     initial_direction,
 *                     "Initial direction of the uniform flow field");
 *   
 *       initial_1d_state[0] = ProblemDescription<dim>::gamma;
 *       initial_1d_state[1] = 3.;
 *       initial_1d_state[2] = 1.;
 *       add_parameter("initial 1d state",
 *                     initial_1d_state,
 *                     "Initial 1d state (rho, u, p) of the uniform flow field");
 *     }
 *   
 * @endcode
 * 
 * So far the constructor of <code>InitialValues</code> has defined
 * default values for the two private members
 * <code>initial_direction</code> and <code>initial_1d_state</code> and
 * added them to the parameter list. But we have not defined an
 * implementation of the only public member that we really care about,
 * which is <code>initial_state()</code> (the function that we are going to
 * call to actually evaluate the initial solution at the mesh nodes). At
 * the top of the function, we have to ensure that the provided initial
 * direction is not the zero vector.
 *   
 
 * 
 * @note As commented, we could have avoided using the method
 * <code>parse_parameters_call_back </code> and defined a class member
 * <code>setup()</code> in order to define the implementation of
 * <code>initial_state()</code>. But for illustrative purposes we want to
 * document a different way here and use the call back signal from
 * ParameterAcceptor.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     void InitialValues<dim>::parse_parameters_callback()
 *     {
 *       AssertThrow(initial_direction.norm() != 0.,
 *                   ExcMessage(
 *                     "Initial shock front direction is set to the zero vector."));
 *       initial_direction /= initial_direction.norm();
 *   
 * @endcode
 * 
 * Next, we implement the <code>initial_state</code> function object
 * with a lambda function computing a uniform flow field. For this we
 * have to translate a given primitive 1d state (density @f$\rho@f$,
 * velocity @f$u@f$, and pressure @f$p@f$) into a conserved n-dimensional state
 * (density @f$\rho@f$, momentum @f$\mathbf{m}@f$, and total energy @f$E@f$).
 * 
 
 * 
 * 
 * @code
 *       initial_state = [this](const Point<dim> & /*point*/, double /*t*/) {
 *         const double            rho   = initial_1d_state[0];
 *         const double            u     = initial_1d_state[1];
 *         const double            p     = initial_1d_state[2];
 *         static constexpr double gamma = ProblemDescription<dim>::gamma;
 *   
 *         state_type state;
 *   
 *         state[0] = rho;
 *         for (unsigned int i = 0; i < dim; ++i)
 *           state[1 + i] = rho * u * initial_direction[i];
 *   
 *         state[dim + 1] = p / (gamma - 1.) + 0.5 * rho * u * u;
 *   
 *         return state;
 *       };
 *     }
 *   
 * @endcode
 * 
 * 
 * <a name="TheForwardEulerstep"></a> 
 * <h4>The Forward Euler step</h4>
 * 
 
 * 
 * The constructor of the <code>%TimeStepping</code> class does not contain
 * any surprising code:
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     TimeStepping<dim>::TimeStepping(
 *       const MPI_Comm            mpi_communicator,
 *       TimerOutput &             computing_timer,
 *       const OfflineData<dim> &  offline_data,
 *       const InitialValues<dim> &initial_values,
 *       const std::string &       subsection /*= "TimeStepping"*/)
 *       : ParameterAcceptor(subsection)
 *       , mpi_communicator(mpi_communicator)
 *       , computing_timer(computing_timer)
 *       , offline_data(&offline_data)
 *       , initial_values(&initial_values)
 *     {
 *       cfl_update = 0.80;
 *       add_parameter("cfl update",
 *                     cfl_update,
 *                     "Relative CFL constant used for update");
 *     }
 *   
 * @endcode
 * 
 * In the class member <code>prepare()</code> we initialize the temporary
 * vector <code>temp</code> and the matrix <code>dij_matrix</code>. The
 * vector <code>temp</code> will be used to store the solution update
 * temporarily before its contents is swapped with the old vector.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     void TimeStepping<dim>::prepare()
 *     {
 *       TimerOutput::Scope scope(computing_timer,
 *                                "time_stepping - prepare scratch space");
 *   
 *       for (auto &it : temporary_vector)
 *         it.reinit(offline_data->partitioner);
 *   
 *       dij_matrix.reinit(offline_data->sparsity_pattern);
 *     }
 *   
 * @endcode
 * 
 * It is now time to implement the forward Euler step. Given a (writable
 * reference) to the old state <code>U</code> at time @f$t@f$ we update the
 * state <code>U</code> in place and return the chosen time-step size. We
 * first declare a number of read-only references to various different
 * variables and data structures. We do this is mainly to have shorter
 * variable names (e.g., <code>sparsity</code> instead of
 * <code>offline_data->sparsity_pattern</code>).
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     double TimeStepping<dim>::make_one_step(vector_type &U, const double t)
 *     {
 *       const auto &n_locally_owned    = offline_data->n_locally_owned;
 *       const auto &n_locally_relevant = offline_data->n_locally_relevant;
 *   
 *       const std_cxx20::ranges::iota_view<unsigned int, unsigned int>
 *         indices_owned(0, n_locally_owned);
 *       const std_cxx20::ranges::iota_view<unsigned int, unsigned int>
 *         indices_relevant(0, n_locally_relevant);
 *   
 *       const auto &sparsity = offline_data->sparsity_pattern;
 *   
 *       const auto &lumped_mass_matrix = offline_data->lumped_mass_matrix;
 *       const auto &norm_matrix        = offline_data->norm_matrix;
 *       const auto &nij_matrix         = offline_data->nij_matrix;
 *       const auto &cij_matrix         = offline_data->cij_matrix;
 *   
 *       const auto &boundary_normal_map = offline_data->boundary_normal_map;
 *   
 * @endcode
 * 
 * <b>Step 1</b>: Computing the @f$d_{ij}@f$ graph viscosity matrix.
 *     
 
 * 
 * It is important to highlight that the viscosity matrix has to be
 * symmetric, i.e., @f$d_{ij} = d_{ji}@f$. In this regard we note here that
 * @f$\int_{\Omega} \nabla \phi_j \phi_i \, \mathrm{d}\mathbf{x}= -
 * \int_{\Omega} \nabla \phi_i \phi_j \, \mathrm{d}\mathbf{x}@f$ (or
 * equivalently @f$\mathbf{c}_{ij} = - \mathbf{c}_{ji}@f$) provided either
 * @f$\mathbf{x}_i@f$ or @f$\mathbf{x}_j@f$ is a support point located away
 * from the boundary. In this case we can check that
 * @f$\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})@f$ by construction, which guarantees
 * the property @f$d_{ij} = d_{ji}@f$.
 *     
 
 * 
 * However, if both support points @f$\mathbf{x}_i@f$ or @f$\mathbf{x}_j@f$
 * happen to lie on the boundary, then, the equalities @f$\mathbf{c}_{ij} =
 * - \mathbf{c}_{ji}@f$ and @f$\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})@f$ do not
 * necessarily hold true. The only mathematically safe solution for this
 * dilemma is to compute both of them @f$d_{ij}@f$ and @f$d_{ji}@f$ and
 * take the maximum.
 *     
 
 * 
 * Overall, the computation of @f$d_{ij}@f$ is quite expensive. In
 * order to save some computing time we exploit the fact that the viscosity
 * matrix has to be symmetric (as mentioned above): we only compute
 * the upper-triangular entries of @f$d_{ij}@f$ and copy the
 * corresponding entries to the lower-triangular counterpart.
 *     
 
 * 
 * We use again parallel::apply_to_subranges() for thread-parallel for
 * loops. Pretty much all the ideas for parallel traversal that we
 * introduced when discussing the assembly of the matrix
 * <code>norm_matrix</code> and the normalization of
 * <code>nij_matrix</code> above are used here again.
 *     
 
 * 
 * We define again a "worker" function <code>on_subranges</code> that
 * computes the viscosity @f$d_{ij}@f$ for a subrange `[*index_begin,
 * *index_end)` of column indices:
 * 
 * @code
 *       {
 *         TimerOutput::Scope scope(computing_timer,
 *                                  "time_stepping - 1 compute d_ij");
 *   
 *         const auto on_subranges = 
 *           [&](const auto index_begin, const auto index_end) {
 *             for (const auto i :
 *                  std_cxx20::ranges::iota_view<unsigned int, unsigned int>(
 *                    *index_begin, *index_end))
 *               {
 *                 const auto U_i = gather(U, i);
 *   
 * @endcode
 * 
 * For a given column index i we iterate over the columns of the
 * sparsity pattern from <code>sparsity.begin(i)</code> to
 * <code>sparsity.end(i)</code>:
 * 
 * @code
 *                 for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
 *                   {
 *                     const auto j = jt->column();
 *   
 * @endcode
 * 
 * We only compute @f$d_{ij}@f$ if @f$j < i@f$ (upper triangular
 * entries) and later copy the values over to @f$d_{ji}@f$.
 * 
 * @code
 *                     if (j >= i)
 *                       continue;
 *   
 *                     const auto U_j = gather(U, j);
 *   
 *                     const auto   n_ij = gather_get_entry(nij_matrix, jt);
 *                     const double norm = get_entry(norm_matrix, jt);
 *   
 *                     const auto lambda_max =
 *                       ProblemDescription<dim>::compute_lambda_max(U_i, U_j, n_ij);
 *   
 *                     double d = norm * lambda_max;
 *   
 * @endcode
 * 
 * If both support points happen to be at the boundary we
 * have to compute @f$d_{ji}@f$ as well and then take
 * @f$\max(d_{ij},d_{ji})@f$. After this we can finally set the
 * upper triangular and lower triangular entries.
 * 
 * @code
 *                     if (boundary_normal_map.count(i) != 0 &&
 *                         boundary_normal_map.count(j) != 0)
 *                       {
 *                         const auto n_ji = gather(nij_matrix, j, i);
 *                         const auto lambda_max_2 =
 *                           ProblemDescription<dim>::compute_lambda_max(U_j,
 *                                                                       U_i,
 *                                                                       n_ji);
 *                         const double norm_2 = norm_matrix(j, i);
 *   
 *                         d = std::max(d, norm_2 * lambda_max_2);
 *                       }
 *   
 *                     set_entry(dij_matrix, jt, d);
 *                     dij_matrix(j, i) = d;
 *                   }
 *               }
 *           };
 *   
 *         parallel::apply_to_subranges(indices_relevant.begin(),
 *                                      indices_relevant.end(),
 *                                      on_subranges,
 *                                      4096);
 *       }
 *   
 * @endcode
 * 
 * <b>Step 2</b>: Compute diagonal entries @f$d_{ii}@f$ and
 * @f$\tau_{\text{max}}@f$.
 * 
 
 * 
 * So far we have computed all off-diagonal entries of the matrix
 * <code>dij_matrix</code>. We still have to fill its diagonal entries
 * defined as @f$d_{ii}^n = - \sum_{j \in \mathcal{I}(i)\backslash \{i\}}
 * d_{ij}^n@f$. We use again parallel::apply_to_subranges() for this
 * purpose. While computing the @f$d_{ii}@f$s we also determine the
 * largest admissible time-step, which is defined as
 * \f[
 * \tau_n \dealcoloneq c_{\text{cfl}}\,\min_{i\in\mathcal{V}}
 * \left(\frac{m_i}{-2\,d_{ii}^{n}}\right) \, .
 * \f]
 * Note that the operation @f$\min_{i \in \mathcal{V}}@f$ is intrinsically
 * global, it operates on all nodes: first we have to take the minimum
 * over all threads (of a given node) and then we have to take the
 * minimum over all MPI processes. In the current implementation:
 * - We store  <code>tau_max</code> (per node) as
 * <a
 * href="http://www.cplusplus.com/reference/atomic/atomic/"><code>std::atomic<double></code></a>.
 * The internal implementation of <code>std::atomic</code> will take
 * care of guarding any possible race condition when more than one
 * thread attempts to read and/or write <code>tau_max</code> at the
 * same time.
 * - In order to take the minimum over all MPI process we use the utility
 * function <code>Utilities::MPI::min</code>.
 * 
 
 * 
 * 
 * @code
 *       std::atomic<double> tau_max{std::numeric_limits<double>::infinity()};
 *   
 *       {
 *         TimerOutput::Scope scope(computing_timer,
 *                                  "time_stepping - 2 compute d_ii, and tau_max");
 *   
 * @endcode
 * 
 * on_subranges() will be executed on every thread individually. The
 * variable <code>tau_max_on_subrange</code> is thus stored thread
 * locally.
 * 
 
 * 
 * 
 * @code
 *         const auto on_subranges = 
 *           [&](const auto index_begin, const auto index_end) {
 *             double tau_max_on_subrange = std::numeric_limits<double>::infinity();
 *   
 *             for (const auto i :
 *                  std_cxx20::ranges::iota_view<unsigned int, unsigned int>(
 *                    *index_begin, *index_end))
 *               {
 *                 double d_sum = 0.;
 *   
 *                 for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
 *                   {
 *                     const auto j = jt->column();
 *   
 *                     if (j == i)
 *                       continue;
 *   
 *                     d_sum -= get_entry(dij_matrix, jt);
 *                   }
 *   
 * @endcode
 * 
 * We store the negative sum of the d_ij entries at the
 * diagonal position
 * 
 * @code
 *                 dij_matrix.diag_element(i) = d_sum;
 * @endcode
 * 
 * and compute the maximal local time-step size
 * <code>tau</code>:
 * 
 * @code
 *                 const double mass   = lumped_mass_matrix.diag_element(i);
 *                 const double tau    = cfl_update * mass / (-2. * d_sum);
 *                 tau_max_on_subrange = std::min(tau_max_on_subrange, tau);
 *               }
 *   
 * @endcode
 * 
 * <code>tau_max_on_subrange</code> contains the largest possible
 * time-step size computed for the (thread local) subrange. At this
 * point we have to synchronize the value over all threads. This is
 * were we use the <a
 * href="http://www.cplusplus.com/reference/atomic/atomic/"><code>std::atomic<double></code></a>
 * <i>compare exchange</i> update mechanism:
 * 
 * @code
 *             double current_tau_max = tau_max.load();
 *             while (current_tau_max > tau_max_on_subrange &&
 *                    !tau_max.compare_exchange_weak(current_tau_max,
 *                                                   tau_max_on_subrange))
 *               ;
 *           };
 *   
 *         parallel::apply_to_subranges(indices_relevant.begin(),
 *                                      indices_relevant.end(),
 *                                      on_subranges,
 *                                      4096);
 *   
 * @endcode
 * 
 * After all threads have finished we can simply synchronize the
 * value over all MPI processes:
 * 
 
 * 
 * 
 * @code
 *         tau_max.store(Utilities::MPI::min(tau_max.load(), mpi_communicator));
 *   
 * @endcode
 * 
 * This is a good point to verify that the computed
 * <code>tau_max</code> is indeed a valid floating point number.
 * 
 * @code
 *         AssertThrow(
 *           !std::isnan(tau_max.load()) && !std::isinf(tau_max.load()) &&
 *             tau_max.load() > 0.,
 *           ExcMessage(
 *             "I'm sorry, Dave. I'm afraid I can't do that. - We crashed."));
 *       }
 *   
 * @endcode
 * 
 * <b>Step 3</b>: Perform update.
 * 
 
 * 
 * At this point, we have computed all viscosity coefficients @f$d_{ij}@f$
 * and we know the maximal admissible time-step size
 * @f$\tau_{\text{max}}@f$. This means we can now compute the update:
 *     
 
 * 
 * \f[
 * \mathbf{U}_i^{n+1} = \mathbf{U}_i^{n} - \frac{\tau_{\text{max}} }{m_i}
 * \sum_{j \in \mathcal{I}(i)} (\mathbb{f}(\mathbf{U}_j^{n}) -
 * \mathbb{f}(\mathbf{U}_i^{n})) \cdot \mathbf{c}_{ij} - d_{ij}
 * (\mathbf{U}_j^{n} - \mathbf{U}_i^{n})
 * \f]
 *     
 
 * 
 * This update formula is slightly different from what was discussed in
 * the introduction (in the pseudo-code). However, it can be shown that
 * both equations are algebraically equivalent (they will produce the
 * same numerical values). We favor this second formula since it has
 * natural cancellation properties that might help avoid numerical
 * artifacts.
 * 
 
 * 
 * 
 * @code
 *       {
 *         TimerOutput::Scope scope(computing_timer,
 *                                  "time_stepping - 3 perform update");
 *   
 *         const auto on_subranges = 
 *           [&](const auto index_begin, const auto index_end) {
 *             for (const auto i :
 *                  boost::make_iterator_range(index_begin, index_end))
 *               {
 *                 Assert(i < n_locally_owned, ExcInternalError());
 *   
 *                 const auto U_i = gather(U, i);
 *   
 *                 const auto   f_i = ProblemDescription<dim>::flux(U_i);
 *                 const double m_i = lumped_mass_matrix.diag_element(i);
 *   
 *                 auto U_i_new = U_i;
 *   
 *                 for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
 *                   {
 *                     const auto j = jt->column();
 *   
 *                     const auto U_j = gather(U, j);
 *                     const auto f_j = ProblemDescription<dim>::flux(U_j);
 *   
 *                     const auto c_ij = gather_get_entry(cij_matrix, jt);
 *                     const auto d_ij = get_entry(dij_matrix, jt);
 *   
 *                     for (unsigned int k = 0; k < n_solution_variables; ++k)
 *                       {
 *                         U_i_new[k] +=
 *                           tau_max / m_i *
 *                           (-(f_j[k] - f_i[k]) * c_ij + d_ij * (U_j[k] - U_i[k]));
 *                       }
 *                   }
 *   
 *                 scatter(temporary_vector, U_i_new, i);
 *               }
 *           };
 *   
 *         parallel::apply_to_subranges(indices_owned.begin(),
 *                                      indices_owned.end(),
 *                                      on_subranges,
 *                                      4096);
 *       }
 *   
 * @endcode
 * 
 * <b>Step 4</b>: Fix up boundary states.
 * 
 
 * 
 * As a last step in the Forward Euler method, we have to fix up all
 * boundary states. As discussed in the intro we
 * - advance in time satisfying no boundary condition at all,
 * - at the end of the time step enforce boundary conditions strongly
 * in a post-processing step.
 *     
 
 * 
 * Here, we compute the correction
 * \f[
 * \mathbf{m}_i \dealcoloneq \mathbf{m}_i - (\boldsymbol{\nu}_i \cdot
 * \mathbf{m}_i) \boldsymbol{\nu}_i,
 * \f]
 * which removes the normal component of @f$\mathbf{m}@f$.
 * 
 
 * 
 * 
 * @code
 *       {
 *         TimerOutput::Scope scope(computing_timer,
 *                                  "time_stepping - 4 fix boundary states");
 *   
 *         for (const auto &it : boundary_normal_map)
 *           {
 *             const auto i = it.first;
 *   
 * @endcode
 * 
 * We only iterate over the locally owned subset:
 * 
 * @code
 *             if (i >= n_locally_owned)
 *               continue;
 *   
 *             const auto &normal   = std::get<0>(it.second);
 *             const auto &id       = std::get<1>(it.second);
 *             const auto &position = std::get<2>(it.second);
 *   
 *             auto U_i = gather(temporary_vector, i);
 *   
 * @endcode
 * 
 * On free slip boundaries we remove the normal component of the
 * momentum:
 * 
 * @code
 *             if (id == Boundaries::free_slip)
 *               {
 *                 auto m = ProblemDescription<dim>::momentum(U_i);
 *                 m -= (m * normal) * normal;
 *                 for (unsigned int k = 0; k < dim; ++k)
 *                   U_i[k + 1] = m[k];
 *               }
 *   
 * @endcode
 * 
 * On Dirichlet boundaries we enforce initial conditions
 * strongly:
 * 
 * @code
 *             else if (id == Boundaries::dirichlet)
 *               {
 *                 U_i = initial_values->initial_state(position, t + tau_max);
 *               }
 *   
 *             scatter(temporary_vector, U_i, i);
 *           }
 *       }
 *   
 * @endcode
 * 
 * <b>Step 5</b>: We now update the ghost layer over all MPI ranks,
 * swap the temporary vector with the solution vector <code>U</code>
 * (that will get returned by reference) and return the chosen
 * time-step size @f$\tau_{\text{max}}@f$:
 * 
 
 * 
 * 
 * @code
 *       for (auto &it : temporary_vector)
 *         it.update_ghost_values();
 *   
 *       U.swap(temporary_vector);
 *   
 *       return tau_max;
 *     }
 *   
 * @endcode
 * 
 * 
 * <a name="Schlierenpostprocessing"></a> 
 * <h4>Schlieren postprocessing</h4>
 *   
 
 * 
 * At various intervals we will output the current state <code>U</code>
 * of the solution together with a so-called Schlieren plot.
 * The constructor of the <code>SchlierenPostprocessor</code> class again
 * contains no surprises. We simply supply default values to and register
 * two parameters:
 * - schlieren_beta:
 * is an ad-hoc positive amplification factor in order to enhance the
 * contrast in the visualization. Its actual value is a matter of
 * taste.
 * - schlieren_index: is an integer indicating which component of the
 * state @f$[\rho, \mathbf{m},E]@f$ are we going to use in order to generate
 * the visualization.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     SchlierenPostprocessor<dim>::SchlierenPostprocessor(
 *       const MPI_Comm          mpi_communicator,
 *       TimerOutput &           computing_timer,
 *       const OfflineData<dim> &offline_data,
 *       const std::string &     subsection /*= "SchlierenPostprocessor"*/)
 *       : ParameterAcceptor(subsection)
 *       , mpi_communicator(mpi_communicator)
 *       , computing_timer(computing_timer)
 *       , offline_data(&offline_data)
 *     {
 *       schlieren_beta = 10.;
 *       add_parameter("schlieren beta",
 *                     schlieren_beta,
 *                     "Beta factor used in Schlieren-type postprocessor");
 *   
 *       schlieren_index = 0;
 *       add_parameter("schlieren index",
 *                     schlieren_index,
 *                     "Use the corresponding component of the state vector for the "
 *                     "schlieren plot");
 *     }
 *   
 * @endcode
 * 
 * Again, the <code>prepare()</code> function initializes two temporary
 * the vectors (<code>r</code> and <code>schlieren</code>).
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     void SchlierenPostprocessor<dim>::prepare()
 *     {
 *       TimerOutput::Scope scope(computing_timer,
 *                                "schlieren_postprocessor - prepare scratch space");
 *   
 *       r.reinit(offline_data->n_locally_relevant);
 *       schlieren.reinit(offline_data->partitioner);
 *     }
 *   
 * @endcode
 * 
 * We now discuss the implementation of the class member
 * <code>SchlierenPostprocessor<dim>::compute_schlieren()</code>, which
 * basically takes a component of the state vector <code>U</code> and
 * computes the Schlieren indicator for such component (the formula of the
 * Schlieren indicator can be found just before the declaration of the class
 * <code>SchlierenPostprocessor</code>). We start by noting
 * that this formula requires the "nodal gradients" @f$\nabla r_j@f$.
 * However, nodal values of gradients are not defined for @f$\mathcal{C}^0@f$
 * finite element functions. More generally, pointwise values of
 * gradients are not defined for @f$W^{1,p}(\Omega)@f$ functions. The
 * simplest technique we can use to recover gradients at nodes is
 * weighted-averaging i.e.
 *   
 
 * 
 * \f[ \nabla r_j \dealcoloneq \frac{1}{\int_{S_i} \omega_i(\mathbf{x}) \,
 * \mathrm{d}\mathbf{x}}
 * \int_{S_i} r_h(\mathbf{x}) \omega_i(\mathbf{x}) \, \mathrm{d}\mathbf{x}
 * \ \ \ \ \ \mathbf{(*)} \f]
 *   
 
 * 
 * where @f$S_i@f$ is the support of the shape function @f$\phi_i@f$, and
 * @f$\omega_i(\mathbf{x})@f$ is the weight. The weight could be any
 * positive function such as
 * @f$\omega_i(\mathbf{x}) \equiv 1@f$ (that would allow us to recover the usual
 * notion of mean value). But as usual, the goal is to reuse the off-line
 * data as much as possible. In this sense, the most natural
 * choice of weight is @f$\omega_i = \phi_i@f$. Inserting this choice of
 * weight and the expansion @f$r_h(\mathbf{x}) = \sum_{j \in \mathcal{V}}
 * r_j \phi_j(\mathbf{x})@f$ into @f$\mathbf{(*)}@f$ we get :
 *   
 
 * 
 * \f[ \nabla r_j \dealcoloneq \frac{1}{m_i} \sum_{j \in \mathcal{I}(i)} r_j
 * \mathbf{c}_{ij} \ \ \ \ \ \mathbf{(**)} \, . \f]
 *   
 
 * 
 * Using this last formula we can recover averaged nodal gradients without
 * resorting to any form of quadrature. This idea aligns quite well with
 * the whole spirit of edge-based schemes (or algebraic schemes) where
 * we want to operate on matrices and vectors as directly as
 * it could be possible avoiding by all means assembly of bilinear
 * forms, cell-loops, quadrature, or any other
 * intermediate construct/operation between the input arguments (the state
 * from the previous time-step) and the actual matrices and vectors
 * required to compute the update.
 *   
 
 * 
 * The second thing to note is that we have to compute global minimum and
 * maximum @f$\max_j |\nabla r_j|@f$ and @f$\min_j |\nabla r_j|@f$. Following the
 * same ideas used to compute the time step size in the class member
 * <code>%TimeStepping\<dim>::%step()</code> we define @f$\max_j |\nabla r_j|@f$
 * and @f$\min_j |\nabla r_j|@f$ as atomic doubles in order to resolve any
 * conflicts between threads. As usual, we use
 * <code>Utilities::MPI::max()</code> and
 * <code>Utilities::MPI::min()</code> to find the global maximum/minimum
 * among all MPI processes.
 *   
 
 * 
 * Finally, it is not possible to compute the Schlieren indicator in a single
 * loop over all nodes. The entire operation requires two loops over nodes:
 *   
 
 * 
 * - The first loop computes @f$|\nabla r_i|@f$ for all @f$i \in \mathcal{V}@f$ in
 * the mesh, and the bounds @f$\max_j |\nabla r_j|@f$ and
 * @f$\min_j |\nabla r_j|@f$.
 * - The second loop finally computes the Schlieren indicator using the
 * formula
 *   
 
 * 
 * \f[ \text{schlieren}[i] = e^{\beta \frac{ |\nabla r_i|
 * - \min_j |\nabla r_j| }{\max_j |\nabla r_j| - \min_j |\nabla r_j| } }
 * \, . \f]
 *   
 
 * 
 * This means that we will have to define two workers
 * <code>on_subranges</code> for each one of these stages.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     void SchlierenPostprocessor<dim>::compute_schlieren(const vector_type &U)
 *     {
 *       TimerOutput::Scope scope(
 *         computing_timer, "schlieren_postprocessor - compute schlieren plot");
 *   
 *       const auto &sparsity            = offline_data->sparsity_pattern;
 *       const auto &lumped_mass_matrix  = offline_data->lumped_mass_matrix;
 *       const auto &cij_matrix          = offline_data->cij_matrix;
 *       const auto &boundary_normal_map = offline_data->boundary_normal_map;
 *       const auto &n_locally_owned     = offline_data->n_locally_owned;
 *   
 *       const auto indices =
 *         std_cxx20::ranges::iota_view<unsigned int, unsigned int>(0,
 *                                                                  n_locally_owned);
 *   
 * @endcode
 * 
 * We define the r_i_max and r_i_min in the current MPI process as
 * atomic doubles in order to avoid race conditions between threads:
 * 
 * @code
 *       std::atomic<double> r_i_max{0.};
 *       std::atomic<double> r_i_min{std::numeric_limits<double>::infinity()};
 *   
 * @endcode
 * 
 * First loop: compute the averaged gradient at each node and the
 * global maxima and minima of the gradients.
 * 
 * @code
 *       {
 *         const auto on_subranges = 
 *           [&](const auto index_begin, const auto index_end) {
 *             double r_i_max_on_subrange = 0.;
 *             double r_i_min_on_subrange = std::numeric_limits<double>::infinity();
 *   
 *             for (const auto i :
 *                  boost::make_iterator_range(index_begin, index_end))
 *               {
 *                 Assert(i < n_locally_owned, ExcInternalError());
 *   
 *                 Tensor<1, dim> r_i;
 *   
 *                 for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
 *                   {
 *                     const auto j = jt->column();
 *   
 *                     if (i == j)
 *                       continue;
 *   
 *                     const auto U_js = U[schlieren_index].local_element(j);
 *                     const auto c_ij = gather_get_entry(cij_matrix, jt);
 *                     r_i += c_ij * U_js;
 *                   }
 *   
 * @endcode
 * 
 * We fix up the gradient r_i at free slip boundaries similarly to
 * how we fixed up boundary states in the forward Euler step.
 * This avoids sharp, artificial gradients in the Schlieren
 * plot at free slip boundaries and is a purely cosmetic choice.
 * 
 
 * 
 * 
 * @code
 *                 const auto bnm_it = boundary_normal_map.find(i);
 *                 if (bnm_it != boundary_normal_map.end())
 *                   {
 *                     const auto &normal = std::get<0>(bnm_it->second);
 *                     const auto &id     = std::get<1>(bnm_it->second);
 *   
 *                     if (id == Boundaries::free_slip)
 *                       r_i -= 1. * (r_i * normal) * normal;
 *                     else
 *                       r_i = 0.;
 *                   }
 *   
 * @endcode
 * 
 * We remind the reader that we are not interested in the nodal
 * gradients per se. We only want their norms in order to
 * compute the Schlieren indicator (weighted with the lumped
 * mass matrix @f$m_i@f$):
 * 
 * @code
 *                 const double m_i    = lumped_mass_matrix.diag_element(i);
 *                 r[i]                = r_i.norm() / m_i;
 *                 r_i_max_on_subrange = std::max(r_i_max_on_subrange, r[i]);
 *                 r_i_min_on_subrange = std::min(r_i_min_on_subrange, r[i]);
 *               }
 *   
 * @endcode
 * 
 * We compare the current_r_i_max and current_r_i_min (in the
 * current subrange) with r_i_max and r_i_min (for the current MPI
 * process) and update them if necessary:
 * 
 
 * 
 * 
 * @code
 *             double current_r_i_max = r_i_max.load();
 *             while (current_r_i_max < r_i_max_on_subrange &&
 *                    !r_i_max.compare_exchange_weak(current_r_i_max,
 *                                                   r_i_max_on_subrange))
 *               ;
 *   
 *             double current_r_i_min = r_i_min.load();
 *             while (current_r_i_min > r_i_min_on_subrange &&
 *                    !r_i_min.compare_exchange_weak(current_r_i_min,
 *                                                   r_i_min_on_subrange))
 *               ;
 *           };
 *   
 *         parallel::apply_to_subranges(indices.begin(),
 *                                      indices.end(),
 *                                      on_subranges,
 *                                      4096);
 *       }
 *   
 * @endcode
 * 
 * And synchronize <code>r_i_max</code> and <code>r_i_min</code> over
 * all MPI processes.
 * 
 
 * 
 * 
 * @code
 *       r_i_max.store(Utilities::MPI::max(r_i_max.load(), mpi_communicator));
 *       r_i_min.store(Utilities::MPI::min(r_i_min.load(), mpi_communicator));
 *   
 * @endcode
 * 
 * Second loop: we now have the vector <code>r</code> and the scalars
 * <code>r_i_max</code> and <code>r_i_min</code> at our disposal. We
 * are thus in a position to actually compute the Schlieren indicator.
 * 
 
 * 
 * 
 * @code
 *       {
 *         const auto on_subranges = 
 *           [&](const auto index_begin, const auto index_end) {
 *             for (const auto i :
 *                  boost::make_iterator_range(index_begin, index_end))
 *               {
 *                 Assert(i < n_locally_owned, ExcInternalError());
 *   
 *                 schlieren.local_element(i) =
 *                   1. - std::exp(-schlieren_beta * (r[i] - r_i_min) /
 *                                 (r_i_max - r_i_min));
 *               }
 *           };
 *   
 *         parallel::apply_to_subranges(indices.begin(),
 *                                      indices.end(),
 *                                      on_subranges,
 *                                      4096);
 *       }
 *   
 * @endcode
 * 
 * And finally, exchange ghost elements.
 * 
 * @code
 *       schlieren.update_ghost_values();
 *     }
 *   
 * @endcode
 * 
 * 
 * <a name="Themainloop"></a> 
 * <h4>The main loop</h4>
 *   
 
 * 
 * With all classes implemented it is time to create an instance of
 * <code>Discretization<dim></code>, <code>OfflineData<dim></code>,
 * <code>InitialValues<dim></code>, <code>%TimeStepping\<dim></code>, and
 * <code>SchlierenPostprocessor<dim></code>, and run the forward Euler
 * step in a loop.
 *   
 
 * 
 * In the constructor of <code>MainLoop<dim></code> we now initialize an
 * instance of all classes, and declare a number of parameters
 * controlling output. Most notable, we declare a boolean parameter
 * <code>resume</code> that will control whether the program attempts to
 * restart from an interrupted computation, or not.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     MainLoop<dim>::MainLoop(const MPI_Comm mpi_communicator)
 *       : ParameterAcceptor("A - MainLoop")
 *       , mpi_communicator(mpi_communicator)
 *       , computing_timer(mpi_communicator,
 *                         timer_output,
 *                         TimerOutput::never,
 *                         TimerOutput::cpu_and_wall_times)
 *       , pcout(std::cout, Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
 *       , discretization(mpi_communicator, computing_timer, "B - Discretization")
 *       , offline_data(mpi_communicator,
 *                      computing_timer,
 *                      discretization,
 *                      "C - OfflineData")
 *       , initial_values("D - InitialValues")
 *       , time_stepping(mpi_communicator,
 *                       computing_timer,
 *                       offline_data,
 *                       initial_values,
 *                       "E - TimeStepping")
 *       , schlieren_postprocessor(mpi_communicator,
 *                                 computing_timer,
 *                                 offline_data,
 *                                 "F - SchlierenPostprocessor")
 *     {
 *       base_name = "test";
 *       add_parameter("basename", base_name, "Base name for all output files");
 *   
 *       t_final = 4.;
 *       add_parameter("final time", t_final, "Final time");
 *   
 *       output_granularity = 0.02;
 *       add_parameter("output granularity",
 *                     output_granularity,
 *                     "time interval for output");
 *   
 *       asynchronous_writeback = true;
 *       add_parameter("asynchronous writeback",
 *                     asynchronous_writeback,
 *                     "Write out solution in a background thread performing IO");
 *   
 *       resume = false;
 *       add_parameter("resume", resume, "Resume an interrupted computation.");
 *     }
 *   
 * @endcode
 * 
 * We start by implementing a helper function <code>print_head()</code>
 * in an anonymous namespace that is used to output messages in the
 * terminal with some nice formatting.
 * 
 
 * 
 * 
 * @code
 *     namespace
 *     {
 *       void print_head(ConditionalOStream &pcout,
 *                       const std::string & header,
 *                       const std::string & secondary = "")
 *       {
 *         const auto header_size   = header.size();
 *         const auto padded_header = std::string((34 - header_size) / 2, ' ') +
 *                                    header +
 *                                    std::string((35 - header_size) / 2, ' ');
 *   
 *         const auto secondary_size = secondary.size();
 *         const auto padded_secondary =
 *           std::string((34 - secondary_size) / 2, ' ') + secondary +
 *           std::string((35 - secondary_size) / 2, ' ');
 *   
 *         /* clang-format off */
 *         pcout << std::endl;
 *         pcout << "    ####################################################" << std::endl;
 *         pcout << "    #########                                  #########" << std::endl;
 *         pcout << "    #########"     <<  padded_header   <<     "#########" << std::endl;
 *         pcout << "    #########"     << padded_secondary <<     "#########" << std::endl;
 *         pcout << "    #########                                  #########" << std::endl;
 *         pcout << "    ####################################################" << std::endl;
 *         pcout << std::endl;
 *         /* clang-format on */
 *       }
 *     } // namespace
 *   
 * @endcode
 * 
 * With <code>print_head</code> in place it is now time to implement the
 * <code>MainLoop<dim>::run()</code> that contains the main loop of our
 * program.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     void MainLoop<dim>::run()
 *     {
 * @endcode
 * 
 * We start by reading in parameters and initializing all objects. We
 * note here that the call to ParameterAcceptor::initialize reads in
 * all parameters from the parameter file (whose name is given as a
 * string argument). ParameterAcceptor handles a global
 * ParameterHandler that is initialized with subsections and parameter
 * declarations for all class instances that are derived from
 * ParameterAceptor. The call to initialize enters the subsection for
 * each each derived class, and sets all variables that were added
 * using ParameterAcceptor::add_parameter()
 * 
 
 * 
 * 
 * @code
 *       pcout << "Reading parameters and allocating objects... " << std::flush;
 *   
 *       ParameterAcceptor::initialize("step-69.prm");
 *       pcout << "done" << std::endl;
 *   
 * @endcode
 * 
 * Next we create the triangulation, assemble all matrices, set up
 * scratch space, and initialize the DataOut<dim> object. All of these
 * operations are pretty standard and discussed in detail in the
 * Discretization and OfflineData classes.
 * 
 
 * 
 * 
 * @code
 *       {
 *         print_head(pcout, "create triangulation");
 *   
 *         discretization.setup();
 *   
 *         if (resume)
 *           discretization.triangulation.load(base_name + "-checkpoint.mesh");
 *         else
 *           discretization.triangulation.refine_global(discretization.refinement);
 *   
 *         pcout << "Number of active cells:       "
 *               << discretization.triangulation.n_global_active_cells()
 *               << std::endl;
 *   
 *         print_head(pcout, "compute offline data");
 *         offline_data.setup();
 *         offline_data.assemble();
 *   
 *         pcout << "Number of degrees of freedom: "
 *               << offline_data.dof_handler.n_dofs() << std::endl;
 *   
 *         print_head(pcout, "set up time step");
 *         time_stepping.prepare();
 *         schlieren_postprocessor.prepare();
 *       }
 *   
 * @endcode
 * 
 * We will store the current time and state in the variable
 * <code>t</code> and vector <code>U</code>:
 * 
 
 * 
 * 
 * @code
 *       double       t            = 0.;
 *       unsigned int output_cycle = 0;
 *   
 *       vector_type U;
 *       for (auto &it : U)
 *         it.reinit(offline_data.partitioner);
 *   
 * @endcode
 * 
 * 
 * <a name="Resume"></a> 
 * <h5>Resume</h5>
 *     
 
 * 
 * By default the boolean <code>resume</code> is set to false, i.e. the
 * following code snippet is not run. However, if
 * <code>resume==true</code> we indicate that we have indeed an
 * interrupted computation and the program shall restart by reading in
 * an old state consisting of <code>t</code>,
 * <code>output_cycle</code>, and the state vector <code>U</code> from
 * checkpoint files.
 *     
 
 * 
 * A this point we have already read in the stored refinement history
 * of our parallel distributed mesh. What is missing are the actual
 * state vector <code>U</code>, the time and output cycle. We use the
 * SolutionTransfer class in combination with the
 * distributed::Triangulation::load() /
 * distributed::Triangulation::save() mechanism to read in the state
 * vector. A separate <code>boost::archive</code> is used to retrieve
 * <code>t</code> and <code>output_cycle</code>. The checkpoint files
 * will be created in the <code>output()</code> routine discussed
 * below.
 * 
 
 * 
 * 
 * @code
 *       if (resume)
 *         {
 *           print_head(pcout, "resume interrupted computation");
 *   
 *           parallel::distributed::
 *             SolutionTransfer<dim, LinearAlgebra::distributed::Vector<double>>
 *               solution_transfer(offline_data.dof_handler);
 *   
 *           std::vector<LinearAlgebra::distributed::Vector<double> *> vectors;
 *           std::transform(U.begin(),
 *                          U.end(),
 *                          std::back_inserter(vectors),
 *                          [](auto &it) { return &it; });
 *           solution_transfer.deserialize(vectors);
 *   
 *           for (auto &it : U)
 *             it.update_ghost_values();
 *   
 *           std::ifstream file(base_name + "-checkpoint.metadata",
 *                              std::ios::binary);
 *   
 *           boost::archive::binary_iarchive ia(file);
 *           ia >> t >> output_cycle;
 *         }
 *       else
 *         {
 *           print_head(pcout, "interpolate initial values");
 *           U = interpolate_initial_values();
 *         }
 *   
 * @endcode
 * 
 * With either the initial state set up, or an interrupted state
 * restored it is time to enter the main loop:
 * 
 
 * 
 * 
 * @code
 *       output(U, base_name, t, output_cycle++);
 *   
 *       print_head(pcout, "enter main loop");
 *   
 *       unsigned int timestep_number = 1;
 *       while (t < t_final)
 *         {
 * @endcode
 * 
 * We first print an informative status message
 * 
 
 * 
 * 
 * @code
 *           std::ostringstream head;
 *           std::ostringstream secondary;
 *   
 *           head << "Cycle  " << Utilities::int_to_string(timestep_number, 6)
 *                << "  ("                                                   
 *                << std::fixed << std::setprecision(1) << t / t_final * 100 
 *                << "%)";
 *           secondary << "at time t = " << std::setprecision(8) << std::fixed << t;
 *   
 *           print_head(pcout, head.str(), secondary.str());
 *   
 * @endcode
 * 
 * and then perform a single forward Euler step. Note that the
 * state vector <code>U</code> is updated in place and that
 * <code>time_stepping.make_one_step()</code> returns the chosen step
 * size.
 * 
 
 * 
 * 
 * @code
 *           t += time_stepping.make_one_step(U, t);
 *   
 * @endcode
 * 
 * Post processing, generating output and writing out the current
 * state is a CPU and IO intensive task that we cannot afford to do
 * every time step - in particular with explicit time stepping. We
 * thus only schedule output by calling the <code>output()</code>
 * function if we are past a threshold set by
 * <code>output_granularity</code>.
 * 
 
 * 
 * 
 * @code
 *           if (t > output_cycle * output_granularity)
 *             {
 *               checkpoint(U, base_name, t, output_cycle);
 *               output(U, base_name, t, output_cycle);
 *               ++output_cycle;
 *             }
 *   
 *           ++timestep_number;
 *         }
 *   
 * @endcode
 * 
 * We wait for any remaining background output thread to finish before
 * printing a summary and exiting.
 * 
 * @code
 *       if (background_thread_state.valid())
 *         background_thread_state.wait();
 *   
 *       computing_timer.print_summary();
 *       pcout << timer_output.str() << std::endl;
 *     }
 *   
 * @endcode
 * 
 * The <code>interpolate_initial_values</code> takes an initial time "t"
 * as input argument and populates a state vector <code>U</code> with the
 * help of the <code>InitialValues<dim>::initial_state</code> object.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     typename MainLoop<dim>::vector_type
 *     MainLoop<dim>::interpolate_initial_values(const double t)
 *     {
 *       pcout << "MainLoop<dim>::interpolate_initial_values(t = " << t << ')'
 *             << std::endl;
 *       TimerOutput::Scope scope(computing_timer,
 *                                "main_loop - setup scratch space");
 *   
 *       vector_type U;
 *   
 *       for (auto &it : U)
 *         it.reinit(offline_data.partitioner);
 *   
 *       constexpr auto n_solution_variables =
 *         ProblemDescription<dim>::n_solution_variables;
 *   
 * @endcode
 * 
 * The function signature of
 * <code>InitialValues<dim>::initial_state</code> is not quite right
 * for VectorTools::interpolate(). We work around this issue by, first,
 * creating a lambda function that for a given position <code>x</code>
 * returns just the value of the <code>i</code>th component. This
 * lambda in turn is converted to a Function<dim> object with the help of
 * the ScalarFunctionFromFunctionObject wrapper.
 * 
 * @code
 *       for (unsigned int i = 0; i < n_solution_variables; ++i)
 *         VectorTools::interpolate(offline_data.dof_handler,
 *                                  ScalarFunctionFromFunctionObject<dim, double>(
 *                                    [&](const Point<dim> &x) {
 *                                      return initial_values.initial_state(x, t)[i];
 *                                    }),
 *                                  U[i]);
 *   
 *       for (auto &it : U)
 *         it.update_ghost_values();
 *   
 *       return U;
 *     }
 *   
 * @endcode
 * 
 * 
 * <a name="Outputandcheckpointing"></a> 
 * <h5>Output and checkpointing</h5>
 * 
 
 * 
 * We checkpoint the current state by doing the precise inverse
 * operation to what we discussed for the <a href="Resume">resume
 * logic</a>:
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     void MainLoop<dim>::checkpoint(const typename MainLoop<dim>::vector_type &U,
 *                                    const std::string &name,
 *                                    const double       t,
 *                                    const unsigned int cycle)
 *     {
 *       print_head(pcout, "checkpoint computation");
 *   
 *       parallel::distributed::
 *         SolutionTransfer<dim, LinearAlgebra::distributed::Vector<double>>
 *           solution_transfer(offline_data.dof_handler);
 *   
 *       std::vector<const LinearAlgebra::distributed::Vector<double> *> vectors;
 *       std::transform(U.begin(),
 *                      U.end(),
 *                      std::back_inserter(vectors),
 *                      [](auto &it) { return &it; });
 *   
 *       solution_transfer.prepare_for_serialization(vectors);
 *   
 *       discretization.triangulation.save(name + "-checkpoint.mesh");
 *   
 *       if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
 *         {
 *           std::ofstream file(name + "-checkpoint.metadata", std::ios::binary);
 *           boost::archive::binary_oarchive oa(file);
 *           oa << t << cycle;
 *         }
 *     }
 *   
 * @endcode
 * 
 * Writing out the final vtk files is quite an IO intensive task that can
 * stall the main loop for a while. In order to avoid this we use an <a
 * href="https://en.wikipedia.org/wiki/Asynchronous_I/O">asynchronous
 * IO</a> strategy by creating a background thread that will perform IO
 * while the main loop is allowed to continue. In order for this to work
 * we have to be mindful of two things:
 * - Before running the <code>output_worker</code> thread, we have to create
 * a copy of the state vector <code>U</code>. We store it in the
 * vector <code>output_vector</code>.
 * - We have to avoid any MPI communication in the background thread,
 * otherwise the program might deadlock. This implies that we have to
 * run the postprocessing outside of the worker thread.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     void MainLoop<dim>::output(const typename MainLoop<dim>::vector_type &U,
 *                                const std::string &                        name,
 *                                const double                               t,
 *                                const unsigned int                         cycle)
 *     {
 *       pcout << "MainLoop<dim>::output(t = " << t << ')' << std::endl;
 *   
 * @endcode
 * 
 * If the asynchronous writeback option is set we launch a background
 * thread performing all the slow IO to disc. In that case we have to
 * make sure that the background thread actually finished running. If
 * not, we have to wait to for it to finish. We launch said background
 * thread with <a
 * href="https://en.cppreference.com/w/cpp/thread/async"><code>std::async()</code></a>
 * that returns a <a
 * href="https://en.cppreference.com/w/cpp/thread/future"><code>std::future</code></a>
 * object. This <code>std::future</code> object contains the return
 * value of the function, which is in our case simply
 * <code>void</code>.
 * 
 
 * 
 * 
 * @code
 *       if (background_thread_state.valid())
 *         {
 *           TimerOutput::Scope timer(computing_timer, "main_loop - stalled output");
 *           background_thread_state.wait();
 *         }
 *   
 *       constexpr auto n_solution_variables =
 *         ProblemDescription<dim>::n_solution_variables;
 *   
 * @endcode
 * 
 * At this point we make a copy of the state vector, run the schlieren
 * postprocessor, and run DataOut<dim>::build_patches() The actual
 * output code is standard: We create a DataOut instance, attach all
 * data vectors we want to output and call
 * DataOut<dim>::build_patches(). There is one twist, however. In order
 * to perform asynchronous IO on a background thread we create the
 * DataOut<dim> object as a shared pointer that we pass on to the
 * worker thread to ensure that once we exit this function and the
 * worker thread finishes the DataOut<dim> object gets destroyed again.
 * 
 
 * 
 * 
 * @code
 *       for (unsigned int i = 0; i < n_solution_variables; ++i)
 *         {
 *           output_vector[i] = U[i];
 *           output_vector[i].update_ghost_values();
 *         }
 *   
 *       schlieren_postprocessor.compute_schlieren(output_vector);
 *   
 *       std::unique_ptr<DataOut<dim>> data_out = std::make_unique<DataOut<dim>>();
 *       data_out->attach_dof_handler(offline_data.dof_handler);
 *   
 *       for (unsigned int i = 0; i < n_solution_variables; ++i)
 *         data_out->add_data_vector(output_vector[i],
 *                                   ProblemDescription<dim>::component_names[i]);
 *   
 *       data_out->add_data_vector(schlieren_postprocessor.schlieren,
 *                                 "schlieren_plot");
 *   
 *       data_out->build_patches(discretization.mapping,
 *                               discretization.finite_element.degree - 1);
 *   
 * @endcode
 * 
 * Next we create a lambda function for the background thread. We <a
 * href="https://en.cppreference.com/w/cpp/language/lambda">capture</a>
 * the <code>this</code> pointer as well as most of the arguments of
 * the output function by value so that we have access to them inside
 * the lambda function.
 *     
 
 * 
 * The first capture argument of the lambda function, `data_out_copy`
 * in essence creates a local variable inside the lambda function into
 * which we "move" the `data_out` variable from above. The way this works
 * is that we create a `std::unique_ptr` above that points to the DataOut
 * object. But we have no use for it any more after this point: We really
 * just want to move ownership from the current function to the lambda
 * function implemented in the following few lines. We could have done
 * this by using a `std::shared_ptr` above and giving the lambda function
 * a copy of that shared pointer; once the current function returns (but
 * maybe while the lambda function is still running), our local shared
 * pointer would go out of scope and stop pointing at the actual
 * object, at which point the lambda function simply becomes the sole
 * owner. But using the `std::unique_ptr` is conceptually cleaner as it
 * makes it clear that the current function's `data_out` variable isn't
 * even pointing to the object any more.
 * 
 * @code
 *       auto output_worker =
 *         [data_out_copy = std::move(data_out), this, name, t, cycle]() {
 *           const DataOutBase::VtkFlags flags(
 *             t, cycle, true, DataOutBase::CompressionLevel::best_speed);
 *           data_out_copy->set_flags(flags);
 *   
 *           data_out_copy->write_vtu_with_pvtu_record(
 *             "", name + "-solution", cycle, mpi_communicator, 6);
 *         };
 *   
 * @endcode
 * 
 * If the asynchronous writeback option is set we launch a new
 * background thread with the help of
 * <a
 * href="https://en.cppreference.com/w/cpp/thread/async"><code>std::async</code></a>
 * function. The function returns a <a
 * href="https://en.cppreference.com/w/cpp/thread/future"><code>std::future</code></a>
 * object that we can use to query the status of the background thread.
 * At this point we can return from the <code>output()</code> function
 * and resume with the time stepping in the main loop - the thread will
 * run in the background.
 * 
 * @code
 *       if (asynchronous_writeback)
 *         {
 *           background_thread_state =
 *             std::async(std::launch::async, std::move(output_worker));
 *         }
 *       else
 *         {
 *           output_worker();
 *         }
 *     }
 *   
 *   } // namespace Step69
 *   
 * @endcode
 * 
 * And finally, the main function.
 * 
 
 * 
 * 
 * @code
 *   int main(int argc, char *argv[])
 *   {
 *     try
 *       {
 *         constexpr int dim = 2;
 *   
 *         using namespace dealii;
 *         using namespace Step69;
 *   
 *         Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
 *   
 *         MPI_Comm      mpi_communicator(MPI_COMM_WORLD);
 *         MainLoop<dim> main_loop(mpi_communicator);
 *   
 *         main_loop.run();
 *       }
 *     catch (std::exception &exc)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Exception on processing: " << std::endl
 *                   << exc.what() << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *     catch (...)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Unknown exception!" << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       };
 *   }
 * @endcode
<a name="Results"></a>
<a name="Results"></a><h1>Results</h1>
 
 
Running the program with default parameters in release mode takes about 1
minute on a 4 core machine (with hyperthreading):
@verbatim
# mpirun -np 4 ./step-69 | tee output
Reading parameters and allocating objects... done
 
    ####################################################
    #########                                  #########
    #########       create triangulation       #########
    #########                                  #########
    ####################################################
 
Number of active cells:       36864
 
    ####################################################
    #########                                  #########
    #########       compute offline data       #########
    #########                                  #########
    ####################################################
 
Number of degrees of freedom: 37376
 
    ####################################################
    #########                                  #########
    #########         set up time step         #########
    #########                                  #########
    ####################################################
 
    ####################################################
    #########                                  #########
    #########    interpolate initial values    #########
    #########                                  #########
    #########                                  #########
    ####################################################
 
TimeLoop<dim>::interpolate_initial_values(t = 0)
TimeLoop<dim>::output(t = 0, checkpoint = 0)
 
    ####################################################
    #########                                  #########
    #########         enter main loop          #########
    #########                                  #########
    #########                                  #########
    ####################################################
 
    ####################################################
    #########                                  #########
    #########      Cycle  000001  (0.0%)       #########
    #########      at time t = 0.00000000      #########
    #########                                  #########
    ####################################################
 
[...]
 
    ####################################################
    #########                                  #########
    #########     Cycle  007553  (100.0%)      #########
    #########      at time t = 3.99984036      #########
    #########                                  #########
    ####################################################
 
TimeLoop<dim>::output(t = 4.00038, checkpoint = 1)
 
+------------------------------------------------------------------------+------------+------------+
| Total CPU time elapsed since start                                     |       357s |            |
|                                                                        |            |            |
| Section                                                    | no. calls |  CPU time  | % of total |
+------------------------------------------------------------+-----------+------------+------------+
| discretization - setup                                     |         1 |     0.113s |         0% |
| offline_data - assemble lumped mass matrix, and c_ij       |         1 |     0.167s |         0% |
| offline_data - compute |c_ij|, and n_ij                    |         1 |   0.00255s |         0% |
| offline_data - create sparsity pattern and set up matrices |         1 |    0.0224s |         0% |
| offline_data - distribute dofs                             |         1 |    0.0617s |         0% |
| offline_data - fix slip boundary c_ij                      |         1 |    0.0329s |         0% |
| schlieren_postprocessor - compute schlieren plot           |       201 |     0.811s |      0.23% |
| schlieren_postprocessor - prepare scratch space            |         1 |   7.6e-05s |         0% |
| time_loop - setup scratch space                            |         1 |     0.127s |         0% |
| time_loop - stalled output                                 |       200 |  0.000685s |         0% |
| time_step - 1 compute d_ij                                 |      7553 |       240s |        67% |
| time_step - 2 compute d_ii, and tau_max                    |      7553 |      11.5s |       3.2% |
| time_step - 3 perform update                               |      7553 |       101s |        28% |
| time_step - 4 fix boundary states                          |      7553 |     0.724s |       0.2% |
| time_step - prepare scratch space                          |         1 |   0.00245s |         0% |
+------------------------------------------------------------+-----------+------------+------------+
@endverbatim
 
One thing that becomes evident is the fact that the program spends two
thirds of the execution time computing the graph viscosity d_ij and about a
third of the execution time in performing the update, where computing the
flux @f$f(U)@f$ is the expensive operation. The preset default resolution is
about 37k gridpoints, which amounts to about 148k spatial degrees of
freedom in 2D. An animated schlieren plot of the solution looks as follows:
 
<img src="https://www.dealii.org/images/steps/developer/step-69.coarse.gif" alt="" height="300">
 
It is evident that 37k gridpoints for the first-order method is nowhere
near the resolution needed to resolve any flow features. For comparison,
here is a "reference" computation with a second-order method and about 9.5M
gridpoints (<a href="https://github.com/conservation-laws/ryujin">github
project page</a>):
 
<img src="https://www.dealii.org/images/steps/developer/step-69.2nd-order.t400.jpg" alt="" height="300">
 
So, we give the first-order method a second chance and run it with about
2.4M gridpoints on a small compute server:
 
@verbatim
# mpirun -np 16 ./step-69 | tee output
 
[...]
 
    ####################################################
    #########                                  #########
    #########     Cycle  070216  (100.0%)      #########
    #########      at time t = 3.99999231      #########
    #########                                  #########
    ####################################################
 
TimeLoop<dim>::output(t = 4.00006, checkpoint = 1)
 
[...]
 
+------------------------------------------------------------------------+------------+------------+
| Total wallclock time elapsed since start                               |  6.75e+03s |            |
|                                                                        |            |            |
| Section                                                    | no. calls |  wall time | % of total |
+------------------------------------------------------------+-----------+------------+------------+
| discretization - setup                                     |         1 |      1.97s |         0% |
| offline_data - assemble lumped mass matrix, and c_ij       |         1 |      1.19s |         0% |
| offline_data - compute |c_ij|, and n_ij                    |         1 |    0.0172s |         0% |
| offline_data - create sparsity pattern and set up matrices |         1 |     0.413s |         0% |
| offline_data - distribute dofs                             |         1 |      1.05s |         0% |
| offline_data - fix slip boundary c_ij                      |         1 |     0.252s |         0% |
| schlieren_postprocessor - compute schlieren plot           |       201 |      1.82s |         0% |
| schlieren_postprocessor - prepare scratch space            |         1 |  0.000497s |         0% |
| time_loop - setup scratch space                            |         1 |      1.45s |         0% |
| time_loop - stalled output                                 |       200 |   0.00342s |         0% |
| time_step - 1 compute d_ij                                 |     70216 |  4.38e+03s |        65% |
| time_step - 2 compute d_ii, and tau_max                    |     70216 |       419s |       6.2% |
| time_step - 3 perform update                               |     70216 |  1.87e+03s |        28% |
| time_step - 4 fix boundary states                          |     70216 |        24s |      0.36% |
| time_step - prepare scratch space                          |         1 |    0.0227s |         0% |
+------------------------------------------------------------+-----------+------------+------------+
@endverbatim
 
And with the following result:
 
<img src="https://www.dealii.org/images/steps/developer/step-69.fine.gif" alt="" height="300">
 
That's substantially better, although of course at the price of having run
the code for roughly 2 hours on 16 cores.
 
 
 
<a name="extensions"></a>
<a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
 
 
The program showcased here is really only first-order accurate, as
discussed above. The pictures above illustrate how much diffusion that
introduces and how far the solution is from one that actually resolves
the features we care about.
 
This can be fixed, but it would exceed what a *tutorial* is about.
Nevertheless, it is worth showing what one can achieve by adding a
second-order scheme. For example, here is a video computed with <a
href=https://conservation-laws.org/>the following research code</a>
that shows (with a different color scheme) a 2d simulation that corresponds
to the cases shown above:
 
@htmlonly
<p align="center">
  <iframe width="560" height="315" src="https://www.youtube.com/embed/xIwJZlsXpZ4"
   frameborder="0"
   allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture"
   allowfullscreen></iframe>
 </p>
@endhtmlonly
 
This simulation was done with 38 million degrees of freedom
(continuous @f$Q_1@f$ finite elements) per component of the solution
vector. The exquisite detail of the solution is remarkable for these
kinds of simulations, including in the sub-sonic region behind the
obstacle.
 
One can also with relative ease further extend this to the 3d case:
 
@htmlonly
<p align="center">
  <iframe width="560" height="315" src="https://www.youtube.com/embed/vBCRAF_c8m8"
   frameborder="0"
   allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture"
   allowfullscreen></iframe>
 </p>
@endhtmlonly
 
Solving this becomes expensive, however: The simulation was done with
1,817 million degrees of freedom (continuous @f$Q_1@f$ finite elements)
per component (for a total of 9.09 billion spatial degrees of freedom)
and ran on 30,720 MPI ranks. The code achieved an average throughput of
969M grid points per second (0.04M gridpoints per second per CPU). The
front and back wall show a "Schlieren plot": the magnitude of the
gradient of the density on an exponential scale from white (low) to
black (high). All other cutplanes and the surface of the obstacle show
the magnitude of the vorticity on a white (low) - yellow (medium) -
red (high) scale. (The scales of the individual cutplanes have been
adjusted for a nicer visualization.)
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-69.cc"
*/