]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add a version of step-67 that uses PETSc TS time stepping. 15420/head
authorWolfgang Bangerth <bangerth@colostate.edu>
Wed, 21 Jun 2023 18:37:24 +0000 (12:37 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Wed, 21 Jun 2023 18:37:24 +0000 (12:37 -0600)
tests/petsc/step-67-with-petsc-ts.cc [new file with mode: 0644]
tests/petsc/step-67-with-petsc-ts.with_petsc=true.with_p4est=true.mpirun=1.output [new file with mode: 0644]

diff --git a/tests/petsc/step-67-with-petsc-ts.cc b/tests/petsc/step-67-with-petsc-ts.cc
new file mode 100644 (file)
index 0000000..82daec0
--- /dev/null
@@ -0,0 +1,2439 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2020 - 2022 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE.md at
+ * the top level directory of deal.II.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ * Author: Martin Kronbichler, 2020
+ */
+
+// A variation of step-67 that uses PETSc's TS library as a
+// time stepper
+
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/time_stepping.h>
+#include <deal.II/base/timer.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/vectorization.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/data_out.h>
+
+#ifdef DEAL_II_WITH_PETSC
+#  include <deal.II/lac/petsc_ts.h>
+#  include <deal.II/lac/petsc_vector.h>
+#endif
+
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+
+#include "../tests.h"
+
+
+// The following file includes the CellwiseInverseMassMatrix data structure
+// that we will use for the mass matrix inversion, the only new include
+// file for this tutorial program:
+#include <deal.II/matrix_free/operators.h>
+
+
+
+namespace Euler_DG
+{
+  using namespace dealii;
+
+  // Similarly to the other matrix-free tutorial programs, we collect all
+  // parameters that control the execution of the program at the top of the
+  // file. Besides the dimension and polynomial degree we want to run with, we
+  // also specify a number of points in the Gaussian quadrature formula we
+  // want to use for the nonlinear terms in the Euler equations. Furthermore,
+  // we specify the time interval for the time-dependent problem, and
+  // implement two different test cases. The first one is an analytical
+  // solution in 2d, whereas the second is a channel flow around a cylinder as
+  // described in the introduction. Depending on the test case, we also change
+  // the final time up to which we run the simulation, and a variable
+  // `output_tick` that specifies in which intervals we want to write output
+  // (assuming that the tick is larger than the time step size).
+  constexpr unsigned int testcase             = 0;
+  constexpr unsigned int dimension            = 2;
+  constexpr unsigned int n_global_refinements = 1;
+  constexpr unsigned int fe_degree            = 5;
+  constexpr unsigned int n_q_points_1d        = fe_degree + 2;
+
+  using Number = double;
+
+  constexpr double gamma       = 1.4;
+  constexpr double final_time  = 1.0;
+  constexpr double output_tick = 0.05;
+
+  // Next off are some details of the time integrator, namely a Courant number
+  // that scales the time step size in terms of the formula $\Delta t =
+  // \text{Cr} n_\text{stages} \frac{h}{(p+1)^{1.5} (\|\mathbf{u} +
+  // c)_\text{max}}$, as well as a selection of a few low-storage Runge--Kutta
+  // methods. We specify the Courant number per stage of the Runge--Kutta
+  // scheme, as this gives a more realistic expression of the numerical cost
+  // for schemes of various numbers of stages.
+  const double courant_number = 0.15 / std::pow(fe_degree, 1.5);
+  enum LowStorageRungeKuttaScheme
+  {
+    stage_3_order_3, /* Kennedy, Carpenter, Lewis, 2000 */
+    stage_5_order_4, /* Kennedy, Carpenter, Lewis, 2000 */
+    stage_7_order_4, /* Tselios, Simos, 2007 */
+    stage_9_order_5, /* Kennedy, Carpenter, Lewis, 2000 */
+  };
+  constexpr LowStorageRungeKuttaScheme lsrk_scheme = stage_5_order_4;
+
+  // Eventually, we select a detail of the spatial discretization, namely the
+  // numerical flux (Riemann solver) at the faces between cells. For this
+  // program, we have implemented a modified variant of the Lax--Friedrichs
+  // flux and the Harten--Lax--van Leer (HLL) flux.
+  enum EulerNumericalFlux
+  {
+    lax_friedrichs_modified,
+    harten_lax_vanleer,
+  };
+  constexpr EulerNumericalFlux numerical_flux_type = lax_friedrichs_modified;
+
+
+
+  // @sect3{Equation data}
+
+  // We now define a class with the exact solution for the test case 0 and one
+  // with a background flow field for test case 1 of the channel. Given that
+  // the Euler equations are a problem with $d+2$ equations in $d$ dimensions,
+  // we need to tell the Function base class about the correct number of
+  // components.
+  template <int dim>
+  class ExactSolution : public Function<dim>
+  {
+  public:
+    ExactSolution(const double time)
+      : Function<dim>(dim + 2, time)
+    {}
+
+    virtual double
+    value(const Point<dim> &p, const unsigned int component = 0) const override;
+  };
+
+
+
+  // As far as the actual function implemented is concerned, the analytical
+  // test case is an isentropic vortex case (see e.g. the book by Hesthaven
+  // and Warburton, Example 6.1 in Section 6.6 on page 209) which fulfills the
+  // Euler equations with zero force term on the right hand side. Given that
+  // definition, we return either the density, the momentum, or the energy
+  // depending on which component is requested. Note that the original
+  // definition of the density involves the $\frac{1}{\gamma -1}$-th power of
+  // some expression. Since `std::pow()` has pretty slow implementations on
+  // some systems, we replace it by logarithm followed by exponentiation (of
+  // base 2), which is mathematically equivalent but usually much better
+  // optimized. This formula might lose accuracy in the last digits
+  // for very small numbers compared to `std::pow()`, but we are happy with
+  // it anyway, since small numbers map to data close to 1.
+  //
+  // For the channel test case, we simply select a density of 1, a velocity of
+  // 0.4 in $x$ direction and zero in the other directions, and an energy that
+  // corresponds to a speed of sound of 1.3 measured against the background
+  // velocity field, computed from the relation $E = \frac{c^2}{\gamma (\gamma
+  // -1)} + \frac 12 \rho \|u\|^2$.
+  template <int dim>
+  double
+  ExactSolution<dim>::value(const Point<dim> & x,
+                            const unsigned int component) const
+  {
+    const double t = this->get_time();
+
+    switch (testcase)
+      {
+        case 0:
+          {
+            Assert(dim == 2, ExcNotImplemented());
+            const double beta = 5;
+
+            Point<dim> x0;
+            x0[0] = 5.;
+            const double radius_sqr =
+              (x - x0).norm_square() - 2. * (x[0] - x0[0]) * t + t * t;
+            const double factor =
+              beta / (numbers::PI * 2) * std::exp(1. - radius_sqr);
+            const double density_log = std::log2(
+              std::abs(1. - (gamma - 1.) / gamma * 0.25 * factor * factor));
+            const double density = std::exp2(density_log * (1. / (gamma - 1.)));
+            const double u       = 1. - factor * (x[1] - x0[1]);
+            const double v       = factor * (x[0] - t - x0[0]);
+
+            if (component == 0)
+              return density;
+            else if (component == 1)
+              return density * u;
+            else if (component == 2)
+              return density * v;
+            else
+              {
+                const double pressure =
+                  std::exp2(density_log * (gamma / (gamma - 1.)));
+                return pressure / (gamma - 1.) +
+                       0.5 * (density * u * u + density * v * v);
+              }
+          }
+
+        case 1:
+          {
+            if (component == 0)
+              return 1.;
+            else if (component == 1)
+              return 0.4;
+            else if (component == dim + 1)
+              return 3.097857142857143;
+            else
+              return 0.;
+          }
+
+        default:
+          Assert(false, ExcNotImplemented());
+          return 0.;
+      }
+  }
+
+
+
+  // @sect3{Low-storage explicit Runge--Kutta time integrators}
+
+  // The next few lines implement a few low-storage variants of Runge--Kutta
+  // methods. These methods have specific Butcher tableaux with coefficients
+  // $b_i$ and $a_i$ as shown in the introduction. As usual in Runge--Kutta
+  // method, we can deduce time steps, $c_i = \sum_{j=1}^{i-2} b_i + a_{i-1}$
+  // from those coefficients. The main advantage of this kind of scheme is the
+  // fact that only two vectors are needed per stage, namely the accumulated
+  // part of the solution $\mathbf{w}$ (that will hold the solution
+  // $\mathbf{w}^{n+1}$ at the new time $t^{n+1}$ after the last stage), the
+  // update vector $\mathbf{r}_i$ that gets evaluated during the stages, plus
+  // one vector $\mathbf{k}_i$ to hold the evaluation of the operator. Such a
+  // Runge--Kutta setup reduces the memory storage and memory access. As the
+  // memory bandwidth is often the performance-limiting factor on modern
+  // hardware when the evaluation of the differential operator is
+  // well-optimized, performance can be improved over standard time
+  // integrators. This is true also when taking into account that a
+  // conventional Runge--Kutta scheme might allow for slightly larger time
+  // steps as more free parameters allow for better stability properties.
+  //
+  // In this tutorial programs, we concentrate on a few variants of
+  // low-storage schemes defined in the article by Kennedy, Carpenter, and
+  // Lewis (2000), as well as one variant described by Tselios and Simos
+  // (2007). There is a large series of other schemes available, which could
+  // be addressed by additional sets of coefficients or slightly different
+  // update formulas.
+  //
+  // We define a single class for the four integrators, distinguished by the
+  // enum described above. To each scheme, we then fill the vectors for the
+  // $b_i$ and $a_i$ to the given variables in the class.
+  class LowStorageRungeKuttaIntegrator
+  {
+  public:
+    LowStorageRungeKuttaIntegrator(const LowStorageRungeKuttaScheme scheme)
+    {
+      TimeStepping::runge_kutta_method lsrk;
+      // First comes the three-stage scheme of order three by Kennedy et al.
+      // (2000). While its stability region is significantly smaller than for
+      // the other schemes, it only involves three stages, so it is very
+      // competitive in terms of the work per stage.
+      switch (scheme)
+        {
+          case stage_3_order_3:
+            {
+              lsrk = TimeStepping::LOW_STORAGE_RK_STAGE3_ORDER3;
+              break;
+            }
+
+            // The next scheme is a five-stage scheme of order four, again
+            // defined in the paper by Kennedy et al. (2000).
+          case stage_5_order_4:
+            {
+              lsrk = TimeStepping::LOW_STORAGE_RK_STAGE5_ORDER4;
+              break;
+            }
+
+            // The following scheme of seven stages and order four has been
+            // explicitly derived for acoustics problems. It is a balance of
+            // accuracy for imaginary eigenvalues among fourth order schemes,
+            // combined with a large stability region. Since DG schemes are
+            // dissipative among the highest frequencies, this does not
+            // necessarily translate to the highest possible time step per
+            // stage. In the context of the present tutorial program, the
+            // numerical flux plays a crucial role in the dissipation and thus
+            // also the maximal stable time step size. For the modified
+            // Lax--Friedrichs flux, this scheme is similar to the
+            // `stage_5_order_4` scheme in terms of step size per stage if only
+            // stability is considered, but somewhat less efficient for the HLL
+            // flux.
+          case stage_7_order_4:
+            {
+              lsrk = TimeStepping::LOW_STORAGE_RK_STAGE7_ORDER4;
+              break;
+            }
+
+            // The last scheme included here is the nine-stage scheme of order
+            // five from Kennedy et al. (2000). It is the most accurate among
+            // the schemes used here, but the higher order of accuracy
+            // sacrifices some stability, so the step length normalized per
+            // stage is less than for the fourth order schemes.
+          case stage_9_order_5:
+            {
+              lsrk = TimeStepping::LOW_STORAGE_RK_STAGE9_ORDER5;
+              break;
+            }
+
+          default:
+            AssertThrow(false, ExcNotImplemented());
+        }
+      TimeStepping::LowStorageRungeKutta<
+        LinearAlgebra::distributed::Vector<Number>>
+        rk_integrator(lsrk);
+      rk_integrator.get_coefficients(ai, bi, ci);
+    }
+
+    unsigned int
+    n_stages() const
+    {
+      return bi.size();
+    }
+
+    // The main function of the time integrator is to go through the stages,
+    // evaluate the operator, prepare the $\mathbf{r}_i$ vector for the next
+    // evaluation, and update the solution vector $\mathbf{w}$. We hand off
+    // the work to the `pde_operator` involved in order to be able to merge
+    // the vector operations of the Runge--Kutta setup with the evaluation of
+    // the differential operator for better performance, so all we do here is
+    // to delegate the vectors and coefficients.
+    //
+    // We separately call the operator for the first stage because we need
+    // slightly modified arguments there: We evaluate the solution from
+    // the old solution $\mathbf{w}^n$ rather than a $\mathbf r_i$ vector, so
+    // the first argument is `solution`. We here let the stage vector
+    // $\mathbf{r}_i$ also hold the temporary result of the evaluation, as it
+    // is not used otherwise. For all subsequent stages, we use the vector
+    // `vec_ki` as the second vector argument to store the result of the
+    // operator evaluation. Finally, when we are at the last stage, we must
+    // skip the computation of the vector $\mathbf{r}_{s+1}$ as there is no
+    // coefficient $a_s$ available (nor will it be used).
+    template <typename VectorType, typename Operator>
+    void
+    perform_time_step(const Operator &pde_operator,
+                      const double    current_time,
+                      const double    time_step,
+                      VectorType &    solution,
+                      VectorType &    vec_ri,
+                      VectorType &    vec_ki) const
+    {
+      AssertDimension(ai.size() + 1, bi.size());
+
+      pde_operator.perform_stage(current_time,
+                                 bi[0] * time_step,
+                                 ai[0] * time_step,
+                                 solution,
+                                 vec_ri,
+                                 solution,
+                                 vec_ri);
+
+      for (unsigned int stage = 1; stage < bi.size(); ++stage)
+        {
+          const double c_i = ci[stage];
+          pde_operator.perform_stage(current_time + c_i * time_step,
+                                     bi[stage] * time_step,
+                                     (stage == bi.size() - 1 ?
+                                        0 :
+                                        ai[stage] * time_step),
+                                     vec_ri,
+                                     vec_ki,
+                                     solution,
+                                     vec_ri);
+        }
+    }
+
+  private:
+    std::vector<double> bi;
+    std::vector<double> ai;
+    std::vector<double> ci;
+  };
+
+
+  // @sect3{Implementation of point-wise operations of the Euler equations}
+
+  // In the following functions, we implement the various problem-specific
+  // operators pertaining to the Euler equations. Each function acts on the
+  // vector of conserved variables $[\rho, \rho\mathbf{u}, E]$ that we hold in
+  // the solution vectors, and computes various derived quantities.
+  //
+  // First out is the computation of the velocity, that we derive from the
+  // momentum variable $\rho \mathbf{u}$ by division by $\rho$. One thing to
+  // note here is that we decorate all those functions with the keyword
+  // `DEAL_II_ALWAYS_INLINE`. This is a special macro that maps to a
+  // compiler-specific keyword that tells the compiler to never create a
+  // function call for any of those functions, and instead move the
+  // implementation <a
+  // href="https://en.wikipedia.org/wiki/Inline_function">inline</a> to where
+  // they are called. This is critical for performance because we call into some
+  // of those functions millions or billions of times: For example, we both use
+  // the velocity for the computation of the flux further down, but also for the
+  // computation of the pressure, and both of these places are evaluated at
+  // every quadrature point of every cell. Making sure these functions are
+  // inlined ensures not only that the processor does not have to execute a jump
+  // instruction into the function (and the corresponding return jump), but also
+  // that the compiler can re-use intermediate information from one function's
+  // context in code that comes after the place where the function was called.
+  // (We note that compilers are generally quite good at figuring out which
+  // functions to inline by themselves. Here is a place where compilers may or
+  // may not have figured it out by themselves but where we know for sure that
+  // inlining is a win.)
+  //
+  // Another trick we apply is a separate variable for the inverse density
+  // $\frac{1}{\rho}$. This enables the compiler to only perform a single
+  // division for the flux, despite the division being used at several
+  // places. As divisions are around ten to twenty times as expensive as
+  // multiplications or additions, avoiding redundant divisions is crucial for
+  // performance. We note that taking the inverse first and later multiplying
+  // with it is not equivalent to a division in floating point arithmetic due
+  // to roundoff effects, so the compiler is not allowed to exchange one way by
+  // the other with standard optimization flags. However, it is also not
+  // particularly difficult to write the code in the right way.
+  //
+  // To summarize, the chosen strategy of always inlining and careful
+  // definition of expensive arithmetic operations allows us to write compact
+  // code without passing all intermediate results around, despite making sure
+  // that the code maps to excellent machine code.
+  template <int dim, typename Number>
+  inline DEAL_II_ALWAYS_INLINE //
+    Tensor<1, dim, Number>
+    euler_velocity(const Tensor<1, dim + 2, Number> &conserved_variables)
+  {
+    const Number inverse_density = Number(1.) / conserved_variables[0];
+
+    Tensor<1, dim, Number> velocity;
+    for (unsigned int d = 0; d < dim; ++d)
+      velocity[d] = conserved_variables[1 + d] * inverse_density;
+
+    return velocity;
+  }
+
+  // The next function computes the pressure from the vector of conserved
+  // variables, using the formula $p = (\gamma - 1) \left(E - \frac 12 \rho
+  // \mathbf{u}\cdot \mathbf{u}\right)$. As explained above, we use the
+  // velocity from the `euler_velocity()` function. Note that we need to
+  // specify the first template argument `dim` here because the compiler is
+  // not able to deduce it from the arguments of the tensor, whereas the
+  // second argument (number type) can be automatically deduced.
+  template <int dim, typename Number>
+  inline DEAL_II_ALWAYS_INLINE //
+    Number
+    euler_pressure(const Tensor<1, dim + 2, Number> &conserved_variables)
+  {
+    const Tensor<1, dim, Number> velocity =
+      euler_velocity<dim>(conserved_variables);
+
+    Number rho_u_dot_u = conserved_variables[1] * velocity[0];
+    for (unsigned int d = 1; d < dim; ++d)
+      rho_u_dot_u += conserved_variables[1 + d] * velocity[d];
+
+    return (gamma - 1.) * (conserved_variables[dim + 1] - 0.5 * rho_u_dot_u);
+  }
+
+  // Here is the definition of the Euler flux function, i.e., the definition
+  // of the actual equation. Given the velocity and pressure (that the
+  // compiler optimization will make sure are done only once), this is
+  // straight-forward given the equation stated in the introduction.
+  template <int dim, typename Number>
+  inline DEAL_II_ALWAYS_INLINE //
+    Tensor<1, dim + 2, Tensor<1, dim, Number>>
+    euler_flux(const Tensor<1, dim + 2, Number> &conserved_variables)
+  {
+    const Tensor<1, dim, Number> velocity =
+      euler_velocity<dim>(conserved_variables);
+    const Number pressure = euler_pressure<dim>(conserved_variables);
+
+    Tensor<1, dim + 2, Tensor<1, dim, Number>> flux;
+    for (unsigned int d = 0; d < dim; ++d)
+      {
+        flux[0][d] = conserved_variables[1 + d];
+        for (unsigned int e = 0; e < dim; ++e)
+          flux[e + 1][d] = conserved_variables[e + 1] * velocity[d];
+        flux[d + 1][d] += pressure;
+        flux[dim + 1][d] =
+          velocity[d] * (conserved_variables[dim + 1] + pressure);
+      }
+
+    return flux;
+  }
+
+  // This next function is a helper to simplify the implementation of the
+  // numerical flux, implementing the action of a tensor of tensors (with
+  // non-standard outer dimension of size `dim + 2`, so the standard overloads
+  // provided by deal.II's tensor classes do not apply here) with another
+  // tensor of the same inner dimension, i.e., a matrix-vector product.
+  template <int n_components, int dim, typename Number>
+  inline DEAL_II_ALWAYS_INLINE //
+    Tensor<1, n_components, Number>
+    operator*(const Tensor<1, n_components, Tensor<1, dim, Number>> &matrix,
+              const Tensor<1, dim, Number> &                         vector)
+  {
+    Tensor<1, n_components, Number> result;
+    for (unsigned int d = 0; d < n_components; ++d)
+      result[d] = matrix[d] * vector;
+    return result;
+  }
+
+  // This function implements the numerical flux (Riemann solver). It gets the
+  // state from the two sides of an interface and the normal vector, oriented
+  // from the side of the solution $\mathbf{w}^-$ towards the solution
+  // $\mathbf{w}^+$. In finite volume methods which rely on piece-wise
+  // constant data, the numerical flux is the central ingredient as it is the
+  // only place where the physical information is entered. In DG methods, the
+  // numerical flux is less central due to the polynomials within the elements
+  // and the physical flux used there. As a result of higher-degree
+  // interpolation with consistent values from both sides in the limit of a
+  // continuous solution, the numerical flux can be seen as a control of the
+  // jump of the solution from both sides to weakly impose continuity. It is
+  // important to realize that a numerical flux alone cannot stabilize a
+  // high-order DG method in the presence of shocks, and thus any DG method
+  // must be combined with further shock-capturing techniques to handle those
+  // cases. In this tutorial, we focus on wave-like solutions of the Euler
+  // equations in the subsonic regime without strong discontinuities where our
+  // basic scheme is sufficient.
+  //
+  // Nonetheless, the numerical flux is decisive in terms of the numerical
+  // dissipation of the overall scheme and influences the admissible time step
+  // size with explicit Runge--Kutta methods. We consider two choices, a
+  // modified Lax--Friedrichs scheme and the widely used Harten--Lax--van Leer
+  // (HLL) flux. For both variants, we first need to get the velocities and
+  // pressures from both sides of the interface and evaluate the physical
+  // Euler flux.
+  //
+  // For the local Lax--Friedrichs flux, the definition is $\hat{\mathbf{F}}
+  // =\frac{\mathbf{F}(\mathbf{w}^-)+\mathbf{F}(\mathbf{w}^+)}{2} +
+  // \frac{\lambda}{2}\left[\mathbf{w}^--\mathbf{w}^+\right]\otimes
+  // \mathbf{n^-}$, where the factor $\lambda =
+  // \max\left(\|\mathbf{u}^-\|+c^-, \|\mathbf{u}^+\|+c^+\right)$ gives the
+  // maximal wave speed and $c = \sqrt{\gamma p / \rho}$ is the speed of
+  // sound. Here, we choose two modifications of that expression for reasons
+  // of computational efficiency, given the small impact of the flux on the
+  // solution. For the above definition of the factor $\lambda$, we would need
+  // to take four square roots, two for the two velocity norms and two for the
+  // speed of sound on either side. The first modification is hence to rather
+  // use $\sqrt{\|\mathbf{u}\|^2+c^2}$ as an estimate of the maximal speed
+  // (which is at most a factor of 2 away from the actual maximum, as shown in
+  // the introduction). This allows us to pull the square root out of the
+  // maximum and get away with a single square root computation. The second
+  // modification is to further relax on the parameter $\lambda$---the smaller
+  // it is, the smaller the dissipation factor (which is multiplied by the
+  // jump in $\mathbf{w}$, which might result in a smaller or bigger
+  // dissipation in the end). This allows us to fit the spectrum into the
+  // stability region of the explicit Runge--Kutta integrator with bigger time
+  // steps. However, we cannot make dissipation too small because otherwise
+  // imaginary eigenvalues grow larger. Finally, the current conservative
+  // formulation is not energy-stable in the limit of $\lambda\to 0$ as it is
+  // not skew-symmetric, and would need additional measures such as split-form
+  // DG schemes in that case.
+  //
+  // For the HLL flux, we follow the formula from literature, introducing an
+  // additional weighting of the two states from Lax--Friedrichs by a
+  // parameter $s$. It is derived from the physical transport directions of
+  // the Euler equations in terms of the current direction of velocity and
+  // sound speed. For the velocity, we here choose a simple arithmetic average
+  // which is sufficient for DG scenarios and moderate jumps in material
+  // parameters.
+  //
+  // Since the numerical flux is multiplied by the normal vector in the weak
+  // form, we multiply by the result by the normal vector for all terms in the
+  // equation. In these multiplications, the `operator*` defined above enables
+  // a compact notation similar to the mathematical definition.
+  //
+  // In this and the following functions, we use variable suffixes `_m` and
+  // `_p` to indicate quantities derived from $\mathbf{w}^-$ and $\mathbf{w}^+$,
+  // i.e., values "here" and "there" relative to the current cell when looking
+  // at a neighbor cell.
+  template <int dim, typename Number>
+  inline DEAL_II_ALWAYS_INLINE //
+    Tensor<1, dim + 2, Number>
+    euler_numerical_flux(const Tensor<1, dim + 2, Number> &u_m,
+                         const Tensor<1, dim + 2, Number> &u_p,
+                         const Tensor<1, dim, Number> &    normal)
+  {
+    const auto velocity_m = euler_velocity<dim>(u_m);
+    const auto velocity_p = euler_velocity<dim>(u_p);
+
+    const auto pressure_m = euler_pressure<dim>(u_m);
+    const auto pressure_p = euler_pressure<dim>(u_p);
+
+    const auto flux_m = euler_flux<dim>(u_m);
+    const auto flux_p = euler_flux<dim>(u_p);
+
+    switch (numerical_flux_type)
+      {
+        case lax_friedrichs_modified:
+          {
+            const auto lambda =
+              0.5 * std::sqrt(std::max(velocity_p.norm_square() +
+                                         gamma * pressure_p * (1. / u_p[0]),
+                                       velocity_m.norm_square() +
+                                         gamma * pressure_m * (1. / u_m[0])));
+
+            return 0.5 * (flux_m * normal + flux_p * normal) +
+                   0.5 * lambda * (u_m - u_p);
+          }
+
+        case harten_lax_vanleer:
+          {
+            const auto avg_velocity_normal =
+              0.5 * ((velocity_m + velocity_p) * normal);
+            const auto   avg_c = std::sqrt(std::abs(
+              0.5 * gamma *
+              (pressure_p * (1. / u_p[0]) + pressure_m * (1. / u_m[0]))));
+            const Number s_pos =
+              std::max(Number(), avg_velocity_normal + avg_c);
+            const Number s_neg =
+              std::min(Number(), avg_velocity_normal - avg_c);
+            const Number inverse_s = Number(1.) / (s_pos - s_neg);
+
+            return inverse_s *
+                   ((s_pos * (flux_m * normal) - s_neg * (flux_p * normal)) -
+                    s_pos * s_neg * (u_m - u_p));
+          }
+
+        default:
+          {
+            Assert(false, ExcNotImplemented());
+            return {};
+          }
+      }
+  }
+
+
+
+  // This and the next function are helper functions to provide compact
+  // evaluation calls as multiple points get batched together via a
+  // VectorizedArray argument (see the step-37 tutorial for details). This
+  // function is used for the subsonic outflow boundary conditions where we
+  // need to set the energy component to a prescribed value. The next one
+  // requests the solution on all components and is used for inflow boundaries
+  // where all components of the solution are set.
+  template <int dim, typename Number>
+  VectorizedArray<Number>
+  evaluate_function(const Function<dim> &                      function,
+                    const Point<dim, VectorizedArray<Number>> &p_vectorized,
+                    const unsigned int                         component)
+  {
+    VectorizedArray<Number> result;
+    for (unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v)
+      {
+        Point<dim> p;
+        for (unsigned int d = 0; d < dim; ++d)
+          p[d] = p_vectorized[d][v];
+        result[v] = function.value(p, component);
+      }
+    return result;
+  }
+
+
+  template <int dim, typename Number, int n_components = dim + 2>
+  Tensor<1, n_components, VectorizedArray<Number>>
+  evaluate_function(const Function<dim> &                      function,
+                    const Point<dim, VectorizedArray<Number>> &p_vectorized)
+  {
+    AssertDimension(function.n_components, n_components);
+    Tensor<1, n_components, VectorizedArray<Number>> result;
+    for (unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v)
+      {
+        Point<dim> p;
+        for (unsigned int d = 0; d < dim; ++d)
+          p[d] = p_vectorized[d][v];
+        for (unsigned int d = 0; d < n_components; ++d)
+          result[d][v] = function.value(p, d);
+      }
+    return result;
+  }
+
+
+
+  // @sect3{The EulerOperation class}
+
+  // This class implements the evaluators for the Euler problem, in analogy to
+  // the `LaplaceOperator` class of step-37 or step-59. Since the present
+  // operator is non-linear and does not require a matrix interface (to be
+  // handed over to preconditioners), we skip the various `vmult` functions
+  // otherwise present in matrix-free operators and only implement an `apply`
+  // function as well as the combination of `apply` with the required vector
+  // updates for the low-storage Runge--Kutta time integrator mentioned above
+  // (called `perform_stage`). Furthermore, we have added three additional
+  // functions involving matrix-free routines, namely one to compute an
+  // estimate of the time step scaling (that is combined with the Courant
+  // number for the actual time step size) based on the velocity and speed of
+  // sound in the elements, one for the projection of solutions (specializing
+  // VectorTools::project() for the DG case), and one to compute the errors
+  // against a possible analytical solution or norms against some background
+  // state.
+  //
+  // The rest of the class is similar to other matrix-free tutorials. As
+  // discussed in the introduction, we provide a few functions to allow a user
+  // to pass in various forms of boundary conditions on different parts of the
+  // domain boundary marked by types::boundary_id variables, as well as
+  // possible body forces.
+  template <int dim, int degree, int n_points_1d>
+  class EulerOperator
+  {
+  public:
+    static constexpr unsigned int n_quadrature_points_1d = n_points_1d;
+
+    EulerOperator(TimerOutput &timer_output);
+
+    void
+    reinit(const Mapping<dim> &mapping, const DoFHandler<dim> &dof_handler);
+
+    void
+    set_inflow_boundary(const types::boundary_id       boundary_id,
+                        std::unique_ptr<Function<dim>> inflow_function);
+
+    void
+    set_subsonic_outflow_boundary(
+      const types::boundary_id       boundary_id,
+      std::unique_ptr<Function<dim>> outflow_energy);
+
+    void
+    set_wall_boundary(const types::boundary_id boundary_id);
+
+    void
+    set_body_force(std::unique_ptr<Function<dim>> body_force);
+
+    void
+    apply(const double                                      current_time,
+          const LinearAlgebra::distributed::Vector<Number> &src,
+          LinearAlgebra::distributed::Vector<Number> &      dst) const;
+
+    void
+    perform_stage(const Number cur_time,
+                  const Number factor_solution,
+                  const Number factor_ai,
+                  const LinearAlgebra::distributed::Vector<Number> &current_ri,
+                  LinearAlgebra::distributed::Vector<Number> &      vec_ki,
+                  LinearAlgebra::distributed::Vector<Number> &      solution,
+                  LinearAlgebra::distributed::Vector<Number> &next_ri) const;
+
+    void
+    project(const Function<dim> &                       function,
+            LinearAlgebra::distributed::Vector<Number> &solution) const;
+
+    std::array<double, 3>
+    compute_errors(
+      const Function<dim> &                             function,
+      const LinearAlgebra::distributed::Vector<Number> &solution) const;
+
+    double
+    compute_cell_transport_speed(
+      const LinearAlgebra::distributed::Vector<Number> &solution) const;
+
+    void
+    initialize_vector(LinearAlgebra::distributed::Vector<Number> &vector) const;
+
+  private:
+    MatrixFree<dim, Number> data;
+
+    TimerOutput &timer;
+
+    std::map<types::boundary_id, std::unique_ptr<Function<dim>>>
+      inflow_boundaries;
+    std::map<types::boundary_id, std::unique_ptr<Function<dim>>>
+                                   subsonic_outflow_boundaries;
+    std::set<types::boundary_id>   wall_boundaries;
+    std::unique_ptr<Function<dim>> body_force;
+
+    void
+    local_apply_inverse_mass_matrix(
+      const MatrixFree<dim, Number> &                   data,
+      LinearAlgebra::distributed::Vector<Number> &      dst,
+      const LinearAlgebra::distributed::Vector<Number> &src,
+      const std::pair<unsigned int, unsigned int> &     cell_range) const;
+
+    void
+    local_apply_cell(
+      const MatrixFree<dim, Number> &                   data,
+      LinearAlgebra::distributed::Vector<Number> &      dst,
+      const LinearAlgebra::distributed::Vector<Number> &src,
+      const std::pair<unsigned int, unsigned int> &     cell_range) const;
+
+    void
+    local_apply_face(
+      const MatrixFree<dim, Number> &                   data,
+      LinearAlgebra::distributed::Vector<Number> &      dst,
+      const LinearAlgebra::distributed::Vector<Number> &src,
+      const std::pair<unsigned int, unsigned int> &     face_range) const;
+
+    void
+    local_apply_boundary_face(
+      const MatrixFree<dim, Number> &                   data,
+      LinearAlgebra::distributed::Vector<Number> &      dst,
+      const LinearAlgebra::distributed::Vector<Number> &src,
+      const std::pair<unsigned int, unsigned int> &     face_range) const;
+  };
+
+
+
+  template <int dim, int degree, int n_points_1d>
+  EulerOperator<dim, degree, n_points_1d>::EulerOperator(TimerOutput &timer)
+    : timer(timer)
+  {}
+
+
+
+  // For the initialization of the Euler operator, we set up the MatrixFree
+  // variable contained in the class. This can be done given a mapping to
+  // describe possible curved boundaries as well as a DoFHandler object
+  // describing the degrees of freedom. Since we use a discontinuous Galerkin
+  // discretization in this tutorial program where no constraints are imposed
+  // strongly on the solution field, we do not need to pass in an
+  // AffineConstraints object and rather use a dummy for the
+  // construction. With respect to quadrature, we want to select two different
+  // ways of computing the underlying integrals: The first is a flexible one,
+  // based on a template parameter `n_points_1d` (that will be assigned the
+  // `n_q_points_1d` value specified at the top of this file). More accurate
+  // integration is necessary to avoid the aliasing problem due to the
+  // variable coefficients in the Euler operator. The second less accurate
+  // quadrature formula is a tight one based on `fe_degree+1` and needed for
+  // the inverse mass matrix. While that formula provides an exact inverse
+  // only on affine element shapes and not on deformed elements, it enables
+  // the fast inversion of the mass matrix by tensor product techniques,
+  // necessary to ensure optimal computational efficiency overall.
+  template <int dim, int degree, int n_points_1d>
+  void
+  EulerOperator<dim, degree, n_points_1d>::reinit(
+    const Mapping<dim> &   mapping,
+    const DoFHandler<dim> &dof_handler)
+  {
+    const std::vector<const DoFHandler<dim> *> dof_handlers = {&dof_handler};
+    const AffineConstraints<double>            dummy;
+    const std::vector<const AffineConstraints<double> *> constraints = {&dummy};
+    const std::vector<Quadrature<1>> quadratures = {QGauss<1>(n_q_points_1d),
+                                                    QGauss<1>(fe_degree + 1)};
+
+    typename MatrixFree<dim, Number>::AdditionalData additional_data;
+    additional_data.mapping_update_flags =
+      (update_gradients | update_JxW_values | update_quadrature_points |
+       update_values);
+    additional_data.mapping_update_flags_inner_faces =
+      (update_JxW_values | update_quadrature_points | update_normal_vectors |
+       update_values);
+    additional_data.mapping_update_flags_boundary_faces =
+      (update_JxW_values | update_quadrature_points | update_normal_vectors |
+       update_values);
+    additional_data.tasks_parallel_scheme =
+      MatrixFree<dim, Number>::AdditionalData::none;
+
+    data.reinit(
+      mapping, dof_handlers, constraints, quadratures, additional_data);
+  }
+
+
+
+  template <int dim, int degree, int n_points_1d>
+  void
+  EulerOperator<dim, degree, n_points_1d>::initialize_vector(
+    LinearAlgebra::distributed::Vector<Number> &vector) const
+  {
+    data.initialize_dof_vector(vector);
+  }
+
+
+
+  // The subsequent four member functions are the ones that must be called from
+  // outside to specify the various types of boundaries. For an inflow boundary,
+  // we must specify all components in terms of density $\rho$, momentum $\rho
+  // \mathbf{u}$ and energy $E$. Given this information, we then store the
+  // function alongside the respective boundary id in a map member variable of
+  // this class. Likewise, we proceed for the subsonic outflow boundaries (where
+  // we request a function as well, which we use to retrieve the energy) and for
+  // wall (no-penetration) boundaries where we impose zero normal velocity (no
+  // function necessary, so we only request the boundary id). For the present
+  // DG code where boundary conditions are solely applied as part of the weak
+  // form (during time integration), the call to set the boundary conditions
+  // can appear both before or after the `reinit()` call to this class. This
+  // is different from continuous finite element codes where the boundary
+  // conditions determine the content of the AffineConstraints object that is
+  // sent into MatrixFree for initialization, thus requiring to be set before
+  // the initialization of the matrix-free data structures.
+  //
+  // The checks added in each of the four function are used to
+  // ensure that boundary conditions are mutually exclusive on the various
+  // parts of the boundary, i.e., that a user does not accidentally designate a
+  // boundary as both an inflow and say a subsonic outflow boundary.
+  template <int dim, int degree, int n_points_1d>
+  void
+  EulerOperator<dim, degree, n_points_1d>::set_inflow_boundary(
+    const types::boundary_id       boundary_id,
+    std::unique_ptr<Function<dim>> inflow_function)
+  {
+    AssertThrow(subsonic_outflow_boundaries.find(boundary_id) ==
+                    subsonic_outflow_boundaries.end() &&
+                  wall_boundaries.find(boundary_id) == wall_boundaries.end(),
+                ExcMessage("You already set the boundary with id " +
+                           std::to_string(static_cast<int>(boundary_id)) +
+                           " to another type of boundary before now setting " +
+                           "it as inflow"));
+    AssertThrow(inflow_function->n_components == dim + 2,
+                ExcMessage("Expected function with dim+2 components"));
+
+    inflow_boundaries[boundary_id] = std::move(inflow_function);
+  }
+
+
+  template <int dim, int degree, int n_points_1d>
+  void
+  EulerOperator<dim, degree, n_points_1d>::set_subsonic_outflow_boundary(
+    const types::boundary_id       boundary_id,
+    std::unique_ptr<Function<dim>> outflow_function)
+  {
+    AssertThrow(inflow_boundaries.find(boundary_id) ==
+                    inflow_boundaries.end() &&
+                  wall_boundaries.find(boundary_id) == wall_boundaries.end(),
+                ExcMessage("You already set the boundary with id " +
+                           std::to_string(static_cast<int>(boundary_id)) +
+                           " to another type of boundary before now setting " +
+                           "it as subsonic outflow"));
+    AssertThrow(outflow_function->n_components == dim + 2,
+                ExcMessage("Expected function with dim+2 components"));
+
+    subsonic_outflow_boundaries[boundary_id] = std::move(outflow_function);
+  }
+
+
+  template <int dim, int degree, int n_points_1d>
+  void
+  EulerOperator<dim, degree, n_points_1d>::set_wall_boundary(
+    const types::boundary_id boundary_id)
+  {
+    AssertThrow(inflow_boundaries.find(boundary_id) ==
+                    inflow_boundaries.end() &&
+                  subsonic_outflow_boundaries.find(boundary_id) ==
+                    subsonic_outflow_boundaries.end(),
+                ExcMessage("You already set the boundary with id " +
+                           std::to_string(static_cast<int>(boundary_id)) +
+                           " to another type of boundary before now setting " +
+                           "it as wall boundary"));
+
+    wall_boundaries.insert(boundary_id);
+  }
+
+
+  template <int dim, int degree, int n_points_1d>
+  void
+  EulerOperator<dim, degree, n_points_1d>::set_body_force(
+    std::unique_ptr<Function<dim>> body_force)
+  {
+    AssertDimension(body_force->n_components, dim);
+
+    this->body_force = std::move(body_force);
+  }
+
+
+
+  // @sect4{Local evaluators}
+
+  // Now we proceed to the local evaluators for the Euler problem. The
+  // evaluators are relatively simple and follow what has been presented in
+  // step-37, step-48, or step-59. The first notable difference is the fact
+  // that we use an FEEvaluation with a non-standard number of quadrature
+  // points. Whereas we previously always set the number of quadrature points
+  // to equal the polynomial degree plus one (ensuring exact integration on
+  // affine element shapes), we now set the number quadrature points as a
+  // separate variable (e.g. the polynomial degree plus two or three halves of
+  // the polynomial degree) to more accurately handle nonlinear terms. Since
+  // the evaluator is fed with the appropriate loop lengths via the template
+  // argument and keeps the number of quadrature points in the whole cell in
+  // the variable FEEvaluation::n_q_points, we now automatically operate on
+  // the more accurate formula without further changes.
+  //
+  // The second difference is due to the fact that we are now evaluating a
+  // multi-component system, as opposed to the scalar systems considered
+  // previously. The matrix-free framework provides several ways to handle the
+  // multi-component case. The variant shown here utilizes an FEEvaluation
+  // object with multiple components embedded into it, specified by the fourth
+  // template argument `dim + 2` for the components in the Euler system. As a
+  // consequence, the return type of FEEvaluation::get_value() is not a scalar
+  // any more (that would return a VectorizedArray type, collecting data from
+  // several elements), but a Tensor of `dim+2` components. The functionality
+  // is otherwise similar to the scalar case; it is handled by a template
+  // specialization of a base class, called FEEvaluationAccess. An alternative
+  // variant would have been to use several FEEvaluation objects, a scalar one
+  // for the density, a vector-valued one with `dim` components for the
+  // momentum, and another scalar evaluator for the energy. To ensure that
+  // those components point to the correct part of the solution, the
+  // constructor of FEEvaluation takes three optional integer arguments after
+  // the required MatrixFree field, namely the number of the DoFHandler for
+  // multi-DoFHandler systems (taking the first by default), the number of the
+  // quadrature point in case there are multiple Quadrature objects (see more
+  // below), and as a third argument the component within a vector system. As
+  // we have a single vector for all components, we would go with the third
+  // argument, and set it to `0` for the density, `1` for the vector-valued
+  // momentum, and `dim+1` for the energy slot. FEEvaluation then picks the
+  // appropriate subrange of the solution vector during
+  // FEEvaluationBase::read_dof_values() and
+  // FEEvaluation::distributed_local_to_global() or the more compact
+  // FEEvaluation::gather_evaluate() and FEEvaluation::integrate_scatter()
+  // calls.
+  //
+  // When it comes to the evaluation of the body force vector, we distinguish
+  // between two cases for efficiency reasons: In case we have a constant
+  // function (derived from Functions::ConstantFunction), we can precompute
+  // the value outside the loop over quadrature points and simply use the
+  // value everywhere. For a more general function, we instead need to call
+  // the `evaluate_function()` method we provided above; this path is more
+  // expensive because we need to access the memory associated with the
+  // quadrature point data.
+  //
+  // The rest follows the other tutorial programs. Since we have implemented
+  // all physics for the Euler equations in the separate `euler_flux()`
+  // function, all we have to do here is to call this function
+  // given the current solution evaluated at quadrature points, returned by
+  // `phi.get_value(q)`, and tell the FEEvaluation object to queue the flux
+  // for testing it by the gradients of the shape functions (which is a Tensor
+  // of outer `dim+2` components, each holding a tensor of `dim` components
+  // for the $x,y,z$ component of the Euler flux). One final thing worth
+  // mentioning is the order in which we queue the data for testing by the
+  // value of the test function, `phi.submit_value()`, in case we are given an
+  // external function: We must do this after calling `phi.get_value(q)`,
+  // because `get_value()` (reading the solution) and `submit_value()`
+  // (queuing the value for multiplication by the test function and summation
+  // over quadrature points) access the same underlying data field. Here it
+  // would be easy to achieve also without temporary variable `w_q` since
+  // there is no mixing between values and gradients. For more complicated
+  // setups, one has to first copy out e.g. both the value and gradient at a
+  // quadrature point and then queue results again by
+  // FEEvaluationBase::submit_value() and FEEvaluationBase::submit_gradient().
+  //
+  // As a final note, we mention that we do not use the first MatrixFree
+  // argument of this function, which is a call-back from MatrixFree::loop().
+  // The interfaces imposes the present list of arguments, but since we are in
+  // a member function where the MatrixFree object is already available as the
+  // `data` variable, we stick with that to avoid confusion.
+  template <int dim, int degree, int n_points_1d>
+  void
+  EulerOperator<dim, degree, n_points_1d>::local_apply_cell(
+    const MatrixFree<dim, Number> &,
+    LinearAlgebra::distributed::Vector<Number> &      dst,
+    const LinearAlgebra::distributed::Vector<Number> &src,
+    const std::pair<unsigned int, unsigned int> &     cell_range) const
+  {
+    FEEvaluation<dim, degree, n_points_1d, dim + 2, Number> phi(data);
+
+    Tensor<1, dim, VectorizedArray<Number>> constant_body_force;
+    const Functions::ConstantFunction<dim> *constant_function =
+      dynamic_cast<Functions::ConstantFunction<dim> *>(body_force.get());
+
+    if (constant_function)
+      constant_body_force = evaluate_function<dim, Number, dim>(
+        *constant_function, Point<dim, VectorizedArray<Number>>());
+
+    for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
+      {
+        phi.reinit(cell);
+        phi.gather_evaluate(src, EvaluationFlags::values);
+
+        for (unsigned int q = 0; q < phi.n_q_points; ++q)
+          {
+            const auto w_q = phi.get_value(q);
+            phi.submit_gradient(euler_flux<dim>(w_q), q);
+            if (body_force.get() != nullptr)
+              {
+                const Tensor<1, dim, VectorizedArray<Number>> force =
+                  constant_function ? constant_body_force :
+                                      evaluate_function<dim, Number, dim>(
+                                        *body_force, phi.quadrature_point(q));
+
+                Tensor<1, dim + 2, VectorizedArray<Number>> forcing;
+                for (unsigned int d = 0; d < dim; ++d)
+                  forcing[d + 1] = w_q[0] * force[d];
+                for (unsigned int d = 0; d < dim; ++d)
+                  forcing[dim + 1] += force[d] * w_q[d + 1];
+
+                phi.submit_value(forcing, q);
+              }
+          }
+
+        phi.integrate_scatter(((body_force.get() != nullptr) ?
+                                 EvaluationFlags::values :
+                                 EvaluationFlags::nothing) |
+                                EvaluationFlags::gradients,
+                              dst);
+      }
+  }
+
+
+
+  // The next function concerns the computation of integrals on interior
+  // faces, where we need evaluators from both cells adjacent to the face. We
+  // associate the variable `phi_m` with the solution component $\mathbf{w}^-$
+  // and the variable `phi_p` with the solution component $\mathbf{w}^+$. We
+  // distinguish the two sides in the constructor of FEFaceEvaluation by the
+  // second argument, with `true` for the interior side and `false` for the
+  // exterior side, with interior and exterior denoting the orientation with
+  // respect to the normal vector.
+  //
+  // Note that the calls FEFaceEvaluation::gather_evaluate() and
+  // FEFaceEvaluation::integrate_scatter() combine the access to the vectors
+  // and the sum factorization parts. This combined operation not only saves a
+  // line of code, but also contains an important optimization: Given that we
+  // use a nodal basis in terms of the Lagrange polynomials in the points of
+  // the Gauss-Lobatto quadrature formula, only $(p+1)^{d-1}$ out of the
+  // $(p+1)^d$ basis functions evaluate to non-zero on each face. Thus, the
+  // evaluator only accesses the necessary data in the vector and skips the
+  // parts which are multiplied by zero. If we had first read the vector, we
+  // would have needed to load all data from the vector, as the call in
+  // isolation would not know what data is required in subsequent
+  // operations. If the subsequent FEFaceEvaluation::evaluate() call requests
+  // values and derivatives, indeed all $(p+1)^d$ vector entries for each
+  // component are needed, as the normal derivative is nonzero for all basis
+  // functions.
+  //
+  // The arguments to the evaluators as well as the procedure is similar to
+  // the cell evaluation. We again use the more accurate (over-)integration
+  // scheme due to the nonlinear terms, specified as the third template
+  // argument in the list. At the quadrature points, we then go to our
+  // free-standing function for the numerical flux. It receives the solution
+  // evaluated at quadrature points from both sides (i.e., $\mathbf{w}^-$ and
+  // $\mathbf{w}^+$), as well as the normal vector onto the minus side. As
+  // explained above, the numerical flux is already multiplied by the normal
+  // vector from the minus side. We need to switch the sign because the
+  // boundary term comes with a minus sign in the weak form derived in the
+  // introduction. The flux is then queued for testing both on the minus sign
+  // and on the plus sign, with switched sign as the normal vector from the
+  // plus side is exactly opposed to the one from the minus side.
+  template <int dim, int degree, int n_points_1d>
+  void
+  EulerOperator<dim, degree, n_points_1d>::local_apply_face(
+    const MatrixFree<dim, Number> &,
+    LinearAlgebra::distributed::Vector<Number> &      dst,
+    const LinearAlgebra::distributed::Vector<Number> &src,
+    const std::pair<unsigned int, unsigned int> &     face_range) const
+  {
+    FEFaceEvaluation<dim, degree, n_points_1d, dim + 2, Number> phi_m(data,
+                                                                      true);
+    FEFaceEvaluation<dim, degree, n_points_1d, dim + 2, Number> phi_p(data,
+                                                                      false);
+
+    for (unsigned int face = face_range.first; face < face_range.second; ++face)
+      {
+        phi_p.reinit(face);
+        phi_p.gather_evaluate(src, EvaluationFlags::values);
+
+        phi_m.reinit(face);
+        phi_m.gather_evaluate(src, EvaluationFlags::values);
+
+        for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+          {
+            const auto numerical_flux =
+              euler_numerical_flux<dim>(phi_m.get_value(q),
+                                        phi_p.get_value(q),
+                                        phi_m.get_normal_vector(q));
+            phi_m.submit_value(-numerical_flux, q);
+            phi_p.submit_value(numerical_flux, q);
+          }
+
+        phi_p.integrate_scatter(EvaluationFlags::values, dst);
+        phi_m.integrate_scatter(EvaluationFlags::values, dst);
+      }
+  }
+
+
+
+  // For faces located at the boundary, we need to impose the appropriate
+  // boundary conditions. In this tutorial program, we implement four cases as
+  // mentioned above. (A fifth case, for supersonic outflow conditions is
+  // discussed in the "Results" section below.) The discontinuous Galerkin
+  // method imposes boundary conditions not as constraints, but only
+  // weakly. Thus, the various conditions are imposed by finding an appropriate
+  // <i>exterior</i> quantity $\mathbf{w}^+$ that is then handed to the
+  // numerical flux function also used for the interior faces. In essence,
+  // we "pretend" a state on the outside of the domain in such a way that
+  // if that were reality, the solution of the PDE would satisfy the boundary
+  // conditions we want.
+  //
+  // For wall boundaries, we need to impose a no-normal-flux condition on the
+  // momentum variable, whereas we use a Neumann condition for the density and
+  // energy with $\rho^+ = \rho^-$ and $E^+ = E^-$. To achieve the no-normal
+  // flux condition, we set the exterior values to the interior values and
+  // subtract two times the velocity in wall-normal direction, i.e., in the
+  // direction of the normal vector.
+  //
+  // For inflow boundaries, we simply set the given Dirichlet data
+  // $\mathbf{w}_\mathrm{D}$ as a boundary value. An alternative would have been
+  // to use $\mathbf{w}^+ = -\mathbf{w}^- + 2 \mathbf{w}_\mathrm{D}$, the
+  // so-called mirror principle.
+  //
+  // The imposition of outflow is essentially a Neumann condition, i.e.,
+  // setting $\mathbf{w}^+ = \mathbf{w}^-$. For the case of subsonic outflow,
+  // we still need to impose a value for the energy, which we derive from the
+  // respective function. A special step is needed for the case of
+  // <i>backflow</i>, i.e., the case where there is a momentum flux into the
+  // domain on the Neumann portion. According to the literature (a fact that can
+  // be derived by appropriate energy arguments), we must switch to another
+  // variant of the flux on inflow parts, see Gravemeier, Comerford,
+  // Yoshihara, Ismail, Wall, "A novel formulation for Neumann inflow
+  // conditions in biomechanics", Int. J. Numer. Meth. Biomed. Eng., vol. 28
+  // (2012). Here, the momentum term needs to be added once again, which
+  // corresponds to removing the flux contribution on the momentum
+  // variables. We do this in a post-processing step, and only for the case
+  // when we both are at an outflow boundary and the dot product between the
+  // normal vector and the momentum (or, equivalently, velocity) is
+  // negative. As we work on data of several quadrature points at once for
+  // SIMD vectorizations, we here need to explicitly loop over the array
+  // entries of the SIMD array.
+  //
+  // In the implementation below, we check for the various types
+  // of boundaries at the level of quadrature points. Of course, we could also
+  // have moved the decision out of the quadrature point loop and treat entire
+  // faces as of the same kind, which avoids some map/set lookups in the inner
+  // loop over quadrature points. However, the loss of efficiency is hardly
+  // noticeable, so we opt for the simpler code here. Also note that the final
+  // `else` clause will catch the case when some part of the boundary was not
+  // assigned any boundary condition via `EulerOperator::set_..._boundary(...)`.
+  template <int dim, int degree, int n_points_1d>
+  void
+  EulerOperator<dim, degree, n_points_1d>::local_apply_boundary_face(
+    const MatrixFree<dim, Number> &,
+    LinearAlgebra::distributed::Vector<Number> &      dst,
+    const LinearAlgebra::distributed::Vector<Number> &src,
+    const std::pair<unsigned int, unsigned int> &     face_range) const
+  {
+    FEFaceEvaluation<dim, degree, n_points_1d, dim + 2, Number> phi(data, true);
+
+    for (unsigned int face = face_range.first; face < face_range.second; ++face)
+      {
+        phi.reinit(face);
+        phi.gather_evaluate(src, EvaluationFlags::values);
+
+        for (unsigned int q = 0; q < phi.n_q_points; ++q)
+          {
+            const auto w_m    = phi.get_value(q);
+            const auto normal = phi.get_normal_vector(q);
+
+            auto rho_u_dot_n = w_m[1] * normal[0];
+            for (unsigned int d = 1; d < dim; ++d)
+              rho_u_dot_n += w_m[1 + d] * normal[d];
+
+            bool at_outflow = false;
+
+            Tensor<1, dim + 2, VectorizedArray<Number>> w_p;
+            const auto boundary_id = data.get_boundary_id(face);
+            if (wall_boundaries.find(boundary_id) != wall_boundaries.end())
+              {
+                w_p[0] = w_m[0];
+                for (unsigned int d = 0; d < dim; ++d)
+                  w_p[d + 1] = w_m[d + 1] - 2. * rho_u_dot_n * normal[d];
+                w_p[dim + 1] = w_m[dim + 1];
+              }
+            else if (inflow_boundaries.find(boundary_id) !=
+                     inflow_boundaries.end())
+              w_p =
+                evaluate_function(*inflow_boundaries.find(boundary_id)->second,
+                                  phi.quadrature_point(q));
+            else if (subsonic_outflow_boundaries.find(boundary_id) !=
+                     subsonic_outflow_boundaries.end())
+              {
+                w_p          = w_m;
+                w_p[dim + 1] = evaluate_function(
+                  *subsonic_outflow_boundaries.find(boundary_id)->second,
+                  phi.quadrature_point(q),
+                  dim + 1);
+                at_outflow = true;
+              }
+            else
+              AssertThrow(false,
+                          ExcMessage("Unknown boundary id, did "
+                                     "you set a boundary condition for "
+                                     "this part of the domain boundary?"));
+
+            auto flux = euler_numerical_flux<dim>(w_m, w_p, normal);
+
+            if (at_outflow)
+              for (unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v)
+                {
+                  if (rho_u_dot_n[v] < -1e-12)
+                    for (unsigned int d = 0; d < dim; ++d)
+                      flux[d + 1][v] = 0.;
+                }
+
+            phi.submit_value(-flux, q);
+          }
+
+        phi.integrate_scatter(EvaluationFlags::values, dst);
+      }
+  }
+
+
+
+  // The next function implements the inverse mass matrix operation. The
+  // algorithms and rationale have been discussed extensively in the
+  // introduction, so we here limit ourselves to the technicalities of the
+  // MatrixFreeOperators::CellwiseInverseMassMatrix class. It does similar
+  // operations as the forward evaluation of the mass matrix, except with a
+  // different interpolation matrix, representing the inverse $S^{-1}$
+  // factors. These represent a change of basis from the specified basis (in
+  // this case, the Lagrange basis in the points of the Gauss--Lobatto
+  // quadrature formula) to the Lagrange basis in the points of the Gauss
+  // quadrature formula. In the latter basis, we can apply the inverse of the
+  // point-wise `JxW` factor, i.e., the quadrature weight times the
+  // determinant of the Jacobian of the mapping from reference to real
+  // coordinates. Once this is done, the basis is changed back to the nodal
+  // Gauss-Lobatto basis again. All of these operations are done by the
+  // `apply()` function below. What we need to provide is the local fields to
+  // operate on (which we extract from the global vector by an FEEvaluation
+  // object) and write the results back to the destination vector of the mass
+  // matrix operation.
+  //
+  // One thing to note is that we added two integer arguments (that are
+  // optional) to the constructor of FEEvaluation, the first being 0
+  // (selecting among the DoFHandler in multi-DoFHandler systems; here, we
+  // only have one) and the second being 1 to make the quadrature formula
+  // selection. As we use the quadrature formula 0 for the over-integration of
+  // nonlinear terms, we use the formula 1 with the default $p+1$ (or
+  // `fe_degree+1` in terms of the variable name) points for the mass
+  // matrix. This leads to square contributions to the mass matrix and ensures
+  // exact integration, as explained in the introduction.
+  template <int dim, int degree, int n_points_1d>
+  void
+  EulerOperator<dim, degree, n_points_1d>::local_apply_inverse_mass_matrix(
+    const MatrixFree<dim, Number> &,
+    LinearAlgebra::distributed::Vector<Number> &      dst,
+    const LinearAlgebra::distributed::Vector<Number> &src,
+    const std::pair<unsigned int, unsigned int> &     cell_range) const
+  {
+    FEEvaluation<dim, degree, degree + 1, dim + 2, Number> phi(data, 0, 1);
+    MatrixFreeOperators::CellwiseInverseMassMatrix<dim, degree, dim + 2, Number>
+      inverse(phi);
+
+    for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
+      {
+        phi.reinit(cell);
+        phi.read_dof_values(src);
+
+        inverse.apply(phi.begin_dof_values(), phi.begin_dof_values());
+
+        phi.set_dof_values(dst);
+      }
+  }
+
+
+
+  // @sect4{The apply() and related functions}
+
+  // We now come to the function which implements the evaluation of the Euler
+  // operator as a whole, i.e., $\mathcal M^{-1} \mathcal L(t, \mathbf{w})$,
+  // calling into the local evaluators presented above. The steps should be
+  // clear from the previous code. One thing to note is that we need to adjust
+  // the time in the functions we have associated with the various parts of
+  // the boundary, in order to be consistent with the equation in case the
+  // boundary data is time-dependent. Then, we call MatrixFree::loop() to
+  // perform the cell and face integrals, including the necessary ghost data
+  // exchange in the `src` vector. The seventh argument to the function,
+  // `true`, specifies that we want to zero the `dst` vector as part of the
+  // loop, before we start accumulating integrals into it. This variant is
+  // preferred over explicitly calling `dst = 0.;` before the loop as the
+  // zeroing operation is done on a subrange of the vector in parts that are
+  // written by the integrals nearby. This enhances data locality and allows
+  // for caching, saving one roundtrip of vector data to main memory and
+  // enhancing performance. The last two arguments to the loop determine which
+  // data is exchanged: Since we only access the values of the shape functions
+  // one faces, typical of first-order hyperbolic problems, and since we have
+  // a nodal basis with nodes at the reference element surface, we only need
+  // to exchange those parts. This again saves precious memory bandwidth.
+  //
+  // Once the spatial operator $\mathcal L$ is applied, we need to make a
+  // second round and apply the inverse mass matrix. Here, we call
+  // MatrixFree::cell_loop() since only cell integrals appear. The cell loop
+  // is cheaper than the full loop as access only goes to the degrees of
+  // freedom associated with the locally owned cells, which is simply the
+  // locally owned degrees of freedom for DG discretizations. Thus, no ghost
+  // exchange is needed here.
+  //
+  // Around all these functions, we put timer scopes to record the
+  // computational time for statistics about the contributions of the various
+  // parts.
+  template <int dim, int degree, int n_points_1d>
+  void
+  EulerOperator<dim, degree, n_points_1d>::apply(
+    const double                                      current_time,
+    const LinearAlgebra::distributed::Vector<Number> &src,
+    LinearAlgebra::distributed::Vector<Number> &      dst) const
+  {
+    {
+      TimerOutput::Scope t(timer, "apply - integrals");
+
+      for (auto &i : inflow_boundaries)
+        i.second->set_time(current_time);
+      for (auto &i : subsonic_outflow_boundaries)
+        i.second->set_time(current_time);
+
+      data.loop(&EulerOperator::local_apply_cell,
+                &EulerOperator::local_apply_face,
+                &EulerOperator::local_apply_boundary_face,
+                this,
+                dst,
+                src,
+                true,
+                MatrixFree<dim, Number>::DataAccessOnFaces::values,
+                MatrixFree<dim, Number>::DataAccessOnFaces::values);
+    }
+
+    {
+      TimerOutput::Scope t(timer, "apply - inverse mass");
+
+      data.cell_loop(&EulerOperator::local_apply_inverse_mass_matrix,
+                     this,
+                     dst,
+                     dst);
+    }
+  }
+
+
+
+  // Let us move to the function that does an entire stage of a Runge--Kutta
+  // update. It calls EulerOperator::apply() followed by some updates
+  // to the vectors, namely `next_ri = solution + factor_ai * k_i` and
+  // `solution += factor_solution * k_i`. Rather than performing these
+  // steps through the vector interfaces, we here present an alternative
+  // strategy that is faster on cache-based architectures. As the memory
+  // consumed by the vectors is often much larger than what fits into caches,
+  // the data has to effectively come from the slow RAM memory. The situation
+  // can be improved by loop fusion, i.e., performing both the updates to
+  // `next_ki` and `solution` within a single sweep. In that case, we would
+  // read the two vectors `rhs` and `solution` and write into `next_ki` and
+  // `solution`, compared to at least 4 reads and two writes in the baseline
+  // case. Here, we go one step further and perform the loop immediately when
+  // the mass matrix inversion has finished on a part of the
+  // vector. MatrixFree::cell_loop() provides a mechanism to attach an
+  // `std::function` both before the loop over cells first touches a vector
+  // entry (which we do not use here, but is e.g. used for zeroing the vector)
+  // and a second `std::function` to be called after the loop last touches
+  // an entry. The callback is in form of a range over the given vector (in
+  // terms of the local index numbering in the MPI universe) that can be
+  // addressed by `local_element()` functions.
+  //
+  // For this second callback, we create a lambda that works on a range and
+  // write the respective update on this range. Ideally, we would add the
+  // `DEAL_II_OPENMP_SIMD_PRAGMA` before the local loop to suggest to the
+  // compiler to SIMD parallelize this loop (which means in practice that we
+  // ensure that there is no overlap, also called aliasing, between the index
+  // ranges of the pointers we use inside the loops). It turns out that at the
+  // time of this writing, GCC 7.2 fails to compile an OpenMP pragma inside a
+  // lambda function, so we comment this pragma out below. If your compiler is
+  // newer, you should be able to uncomment these lines again.
+  //
+  // Note that we select a different code path for the last
+  // Runge--Kutta stage when we do not need to update the `next_ri`
+  // vector. This strategy gives a considerable speedup. Whereas the inverse
+  // mass matrix and vector updates take more than 60% of the computational
+  // time with default vector updates on a 40-core machine, the percentage is
+  // around 35% with the more optimized variant. In other words, this is a
+  // speedup of around a third.
+  template <int dim, int degree, int n_points_1d>
+  void
+  EulerOperator<dim, degree, n_points_1d>::perform_stage(
+    const Number                                      current_time,
+    const Number                                      factor_solution,
+    const Number                                      factor_ai,
+    const LinearAlgebra::distributed::Vector<Number> &current_ri,
+    LinearAlgebra::distributed::Vector<Number> &      vec_ki,
+    LinearAlgebra::distributed::Vector<Number> &      solution,
+    LinearAlgebra::distributed::Vector<Number> &      next_ri) const
+  {
+    {
+      TimerOutput::Scope t(timer, "rk_stage - integrals L_h");
+
+      for (auto &i : inflow_boundaries)
+        i.second->set_time(current_time);
+      for (auto &i : subsonic_outflow_boundaries)
+        i.second->set_time(current_time);
+
+      data.loop(&EulerOperator::local_apply_cell,
+                &EulerOperator::local_apply_face,
+                &EulerOperator::local_apply_boundary_face,
+                this,
+                vec_ki,
+                current_ri,
+                true,
+                MatrixFree<dim, Number>::DataAccessOnFaces::values,
+                MatrixFree<dim, Number>::DataAccessOnFaces::values);
+    }
+
+
+    {
+      TimerOutput::Scope t(timer, "rk_stage - inv mass + vec upd");
+      data.cell_loop(
+        &EulerOperator::local_apply_inverse_mass_matrix,
+        this,
+        next_ri,
+        vec_ki,
+        std::function<void(const unsigned int, const unsigned int)>(),
+        [&](const unsigned int start_range, const unsigned int end_range) {
+          const Number ai = factor_ai;
+          const Number bi = factor_solution;
+          if (ai == Number())
+            {
+              /* DEAL_II_OPENMP_SIMD_PRAGMA */
+              for (unsigned int i = start_range; i < end_range; ++i)
+                {
+                  const Number k_i          = next_ri.local_element(i);
+                  const Number sol_i        = solution.local_element(i);
+                  solution.local_element(i) = sol_i + bi * k_i;
+                }
+            }
+          else
+            {
+              /* DEAL_II_OPENMP_SIMD_PRAGMA */
+              for (unsigned int i = start_range; i < end_range; ++i)
+                {
+                  const Number k_i          = next_ri.local_element(i);
+                  const Number sol_i        = solution.local_element(i);
+                  solution.local_element(i) = sol_i + bi * k_i;
+                  next_ri.local_element(i)  = sol_i + ai * k_i;
+                }
+            }
+        });
+    }
+  }
+
+
+
+  // Having discussed the implementation of the functions that deal with
+  // advancing the solution by one time step, let us now move to functions
+  // that implement other, ancillary operations. Specifically, these are
+  // functions that compute projections, evaluate errors, and compute the speed
+  // of information transport on a cell.
+  //
+  // The first of these functions is essentially equivalent to
+  // VectorTools::project(), just much faster because it is specialized for DG
+  // elements where there is no need to set up and solve a linear system, as
+  // each element has independent basis functions. The reason why we show the
+  // code here, besides a small speedup of this non-critical operation, is that
+  // it shows additional functionality provided by
+  // MatrixFreeOperators::CellwiseInverseMassMatrix.
+  //
+  // The projection operation works as follows: If we denote the matrix of
+  // shape functions evaluated at quadrature points by $S$, the projection on
+  // cell $K$ is an operation of the form $\underbrace{S J^K S^\mathrm
+  // T}_{\mathcal M^K} \mathbf{w}^K = S J^K
+  // \tilde{\mathbf{w}}(\mathbf{x}_q)_{q=1:n_q}$, where $J^K$ is the diagonal
+  // matrix containing the determinant of the Jacobian times the quadrature
+  // weight (JxW), $\mathcal M^K$ is the cell-wise mass matrix, and
+  // $\tilde{\mathbf{w}}(\mathbf{x}_q)_{q=1:n_q}$ is the evaluation of the
+  // field to be projected onto quadrature points. (In reality the matrix $S$
+  // has additional structure through the tensor product, as explained in the
+  // introduction.) This system can now equivalently be written as
+  // $\mathbf{w}^K = \left(S J^K S^\mathrm T\right)^{-1} S J^K
+  // \tilde{\mathbf{w}}(\mathbf{x}_q)_{q=1:n_q} = S^{-\mathrm T}
+  // \left(J^K\right)^{-1} S^{-1} S J^K
+  // \tilde{\mathbf{w}}(\mathbf{x}_q)_{q=1:n_q}$. Now, the term $S^{-1} S$ and
+  // then $\left(J^K\right)^{-1} J^K$ cancel, resulting in the final
+  // expression $\mathbf{w}^K = S^{-\mathrm T}
+  // \tilde{\mathbf{w}}(\mathbf{x}_q)_{q=1:n_q}$. This operation is
+  // implemented by
+  // MatrixFreeOperators::CellwiseInverseMassMatrix::transform_from_q_points_to_basis().
+  // The name is derived from the fact that this projection is simply
+  // the multiplication by $S^{-\mathrm T}$, a basis change from the
+  // nodal basis in the points of the Gaussian quadrature to the given finite
+  // element basis. Note that we call FEEvaluation::set_dof_values() to write
+  // the result into the vector, overwriting previous content, rather than
+  // accumulating the results as typical in integration tasks -- we can do
+  // this because every vector entry has contributions from only a single
+  // cell for discontinuous Galerkin discretizations.
+  template <int dim, int degree, int n_points_1d>
+  void
+  EulerOperator<dim, degree, n_points_1d>::project(
+    const Function<dim> &                       function,
+    LinearAlgebra::distributed::Vector<Number> &solution) const
+  {
+    FEEvaluation<dim, degree, degree + 1, dim + 2, Number> phi(data, 0, 1);
+    MatrixFreeOperators::CellwiseInverseMassMatrix<dim, degree, dim + 2, Number>
+      inverse(phi);
+    solution.zero_out_ghost_values();
+    for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell)
+      {
+        phi.reinit(cell);
+        for (unsigned int q = 0; q < phi.n_q_points; ++q)
+          phi.submit_dof_value(evaluate_function(function,
+                                                 phi.quadrature_point(q)),
+                               q);
+        inverse.transform_from_q_points_to_basis(dim + 2,
+                                                 phi.begin_dof_values(),
+                                                 phi.begin_dof_values());
+        phi.set_dof_values(solution);
+      }
+  }
+
+
+
+  // The next function again repeats functionality also provided by the
+  // deal.II library, namely VectorTools::integrate_difference(). We here show
+  // the explicit code to highlight how the vectorization across several cells
+  // works and how to accumulate results via that interface: Recall that each
+  // <i>lane</i> of the vectorized array holds data from a different cell. By
+  // the loop over all cell batches that are owned by the current MPI process,
+  // we could then fill a VectorizedArray of results; to obtain a global sum,
+  // we would need to further go on and sum across the entries in the SIMD
+  // array. However, such a procedure is not stable as the SIMD array could in
+  // fact not hold valid data for all its lanes. This happens when the number
+  // of locally owned cells is not a multiple of the SIMD width. To avoid
+  // invalid data, we must explicitly skip those invalid lanes when accessing
+  // the data. While one could imagine that we could make it work by simply
+  // setting the empty lanes to zero (and thus, not contribute to a sum), the
+  // situation is more complicated than that: What if we were to compute a
+  // velocity out of the momentum? Then, we would need to divide by the
+  // density, which is zero -- the result would consequently be NaN and
+  // contaminate the result. This trap is avoided by accumulating the results
+  // from the valid SIMD range as we loop through the cell batches, using the
+  // function MatrixFree::n_active_entries_per_cell_batch() to give us the
+  // number of lanes with valid data. It equals VectorizedArray::size() on
+  // most cells, but can be less on the last cell batch if the number of cells
+  // has a remainder compared to the SIMD width.
+  template <int dim, int degree, int n_points_1d>
+  std::array<double, 3>
+  EulerOperator<dim, degree, n_points_1d>::compute_errors(
+    const Function<dim> &                             function,
+    const LinearAlgebra::distributed::Vector<Number> &solution) const
+  {
+    TimerOutput::Scope t(timer, "compute errors");
+    double             errors_squared[3] = {};
+    FEEvaluation<dim, degree, n_points_1d, dim + 2, Number> phi(data, 0, 0);
+
+    for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell)
+      {
+        phi.reinit(cell);
+        phi.gather_evaluate(solution, EvaluationFlags::values);
+        VectorizedArray<Number> local_errors_squared[3] = {};
+        for (unsigned int q = 0; q < phi.n_q_points; ++q)
+          {
+            const auto error =
+              evaluate_function(function, phi.quadrature_point(q)) -
+              phi.get_value(q);
+            const auto JxW = phi.JxW(q);
+
+            local_errors_squared[0] += error[0] * error[0] * JxW;
+            for (unsigned int d = 0; d < dim; ++d)
+              local_errors_squared[1] += (error[d + 1] * error[d + 1]) * JxW;
+            local_errors_squared[2] += (error[dim + 1] * error[dim + 1]) * JxW;
+          }
+        for (unsigned int v = 0; v < data.n_active_entries_per_cell_batch(cell);
+             ++v)
+          for (unsigned int d = 0; d < 3; ++d)
+            errors_squared[d] += local_errors_squared[d][v];
+      }
+
+    Utilities::MPI::sum(errors_squared, MPI_COMM_WORLD, errors_squared);
+
+    std::array<double, 3> errors;
+    for (unsigned int d = 0; d < 3; ++d)
+      errors[d] = std::sqrt(errors_squared[d]);
+
+    return errors;
+  }
+
+
+
+  // This final function of the EulerOperator class is used to estimate the
+  // transport speed, scaled by the mesh size, that is relevant for setting
+  // the time step size in the explicit time integrator. In the Euler
+  // equations, there are two speeds of transport, namely the convective
+  // velocity $\mathbf{u}$ and the propagation of sound waves with sound
+  // speed $c = \sqrt{\gamma p/\rho}$ relative to the medium moving at
+  // velocity $\mathbf u$.
+  //
+  // In the formula for the time step size, we are interested not by
+  // these absolute speeds, but by the amount of time it takes for
+  // information to cross a single cell. For information transported along with
+  // the medium, $\mathbf u$ is scaled by the mesh size,
+  // so an estimate of the maximal velocity can be obtained by computing
+  // $\|J^{-\mathrm T} \mathbf{u}\|_\infty$, where $J$ is the Jacobian of the
+  // transformation from real to the reference domain. Note that
+  // FEEvaluationBase::inverse_jacobian() returns the inverse and transpose
+  // Jacobian, representing the metric term from real to reference
+  // coordinates, so we do not need to transpose it again. We store this limit
+  // in the variable `convective_limit` in the code below.
+  //
+  // The sound propagation is isotropic, so we need to take mesh sizes in any
+  // direction into account. The appropriate mesh size scaling is then given
+  // by the minimal singular value of $J$ or, equivalently, the maximal
+  // singular value of $J^{-1}$. Note that one could approximate this quantity
+  // by the minimal distance between vertices of a cell when ignoring curved
+  // cells. To get the maximal singular value of the Jacobian, the general
+  // strategy would be some LAPACK function. Since all we need here is an
+  // estimate, we can avoid the hassle of decomposing a tensor of
+  // VectorizedArray numbers into several matrices and go into an (expensive)
+  // eigenvalue function without vectorization, and instead use a few
+  // iterations (five in the code below) of the power method applied to
+  // $J^{-1}J^{-\mathrm T}$. The speed of convergence of this method depends
+  // on the ratio of the largest to the next largest eigenvalue and the
+  // initial guess, which is the vector of all ones. This might suggest that
+  // we get slow convergence on cells close to a cube shape where all
+  // lengths are almost the same. However, this slow convergence means that
+  // the result will sit between the two largest singular values, which both
+  // are close to the maximal value anyway. In all other cases, convergence
+  // will be quick. Thus, we can merely hardcode 5 iterations here and be
+  // confident that the result is good.
+  template <int dim, int degree, int n_points_1d>
+  double
+  EulerOperator<dim, degree, n_points_1d>::compute_cell_transport_speed(
+    const LinearAlgebra::distributed::Vector<Number> &solution) const
+  {
+    TimerOutput::Scope t(timer, "compute transport speed");
+    Number             max_transport = 0;
+    FEEvaluation<dim, degree, degree + 1, dim + 2, Number> phi(data, 0, 1);
+
+    for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell)
+      {
+        phi.reinit(cell);
+        phi.gather_evaluate(solution, EvaluationFlags::values);
+        VectorizedArray<Number> local_max = 0.;
+        for (unsigned int q = 0; q < phi.n_q_points; ++q)
+          {
+            const auto solution = phi.get_value(q);
+            const auto velocity = euler_velocity<dim>(solution);
+            const auto pressure = euler_pressure<dim>(solution);
+
+            const auto inverse_jacobian = phi.inverse_jacobian(q);
+            const auto convective_speed = inverse_jacobian * velocity;
+            VectorizedArray<Number> convective_limit = 0.;
+            for (unsigned int d = 0; d < dim; ++d)
+              convective_limit =
+                std::max(convective_limit, std::abs(convective_speed[d]));
+
+            const auto speed_of_sound =
+              std::sqrt(gamma * pressure * (1. / solution[0]));
+
+            Tensor<1, dim, VectorizedArray<Number>> eigenvector;
+            for (unsigned int d = 0; d < dim; ++d)
+              eigenvector[d] = 1.;
+            for (unsigned int i = 0; i < 5; ++i)
+              {
+                eigenvector = transpose(inverse_jacobian) *
+                              (inverse_jacobian * eigenvector);
+                VectorizedArray<Number> eigenvector_norm = 0.;
+                for (unsigned int d = 0; d < dim; ++d)
+                  eigenvector_norm =
+                    std::max(eigenvector_norm, std::abs(eigenvector[d]));
+                eigenvector /= eigenvector_norm;
+              }
+            const auto jac_times_ev   = inverse_jacobian * eigenvector;
+            const auto max_eigenvalue = std::sqrt(
+              (jac_times_ev * jac_times_ev) / (eigenvector * eigenvector));
+            local_max =
+              std::max(local_max,
+                       max_eigenvalue * speed_of_sound + convective_limit);
+          }
+
+        // Similarly to the previous function, we must make sure to accumulate
+        // speed only on the valid cells of a cell batch.
+        for (unsigned int v = 0; v < data.n_active_entries_per_cell_batch(cell);
+             ++v)
+          for (unsigned int d = 0; d < 3; ++d)
+            max_transport = std::max(max_transport, local_max[v]);
+      }
+
+    max_transport = Utilities::MPI::max(max_transport, MPI_COMM_WORLD);
+
+    return max_transport;
+  }
+
+
+
+  // @sect3{The EulerProblem class}
+
+  // This class combines the EulerOperator class with the time integrator and
+  // the usual global data structures such as FiniteElement and DoFHandler, to
+  // actually run the simulations of the Euler problem.
+  //
+  // The member variables are a triangulation, a finite element, a mapping (to
+  // create high-order curved surfaces, see e.g. step-10), and a DoFHandler to
+  // describe the degrees of freedom. In addition, we keep an instance of the
+  // EulerOperator described above around, which will do all heavy lifting in
+  // terms of integrals, and some parameters for time integration like the
+  // current time or the time step size.
+  //
+  // Furthermore, we use a PostProcessor instance to write some additional
+  // information to the output file, in similarity to what was done in
+  // step-33. The interface of the DataPostprocessor class is intuitive,
+  // requiring us to provide information about what needs to be evaluated
+  // (typically only the values of the solution, except for the Schlieren plot
+  // that we only enable in 2d where it makes sense), and the names of what
+  // gets evaluated. Note that it would also be possible to extract most
+  // information by calculator tools within visualization programs such as
+  // ParaView, but it is so much more convenient to do it already when writing
+  // the output.
+  template <int dim>
+  class EulerProblem
+  {
+  public:
+    EulerProblem();
+
+    void
+    run();
+
+  private:
+    void
+    make_grid_and_dofs();
+
+    void
+    output_results(const unsigned int result_number);
+
+    LinearAlgebra::distributed::Vector<Number> solution;
+
+    ConditionalOStream pcout;
+
+#ifdef DEAL_II_WITH_P4EST
+    parallel::distributed::Triangulation<dim> triangulation;
+#else
+    Triangulation<dim> triangulation;
+#endif
+
+    FESystem<dim>   fe;
+    MappingQ<dim>   mapping;
+    DoFHandler<dim> dof_handler;
+
+    TimerOutput timer;
+
+    EulerOperator<dim, fe_degree, n_q_points_1d> euler_operator;
+
+    double time, time_step;
+
+    class Postprocessor : public DataPostprocessor<dim>
+    {
+    public:
+      Postprocessor();
+
+      virtual void
+      evaluate_vector_field(
+        const DataPostprocessorInputs::Vector<dim> &inputs,
+        std::vector<Vector<double>> &computed_quantities) const override;
+
+      virtual std::vector<std::string>
+      get_names() const override;
+
+      virtual std::vector<
+        DataComponentInterpretation::DataComponentInterpretation>
+      get_data_component_interpretation() const override;
+
+      virtual UpdateFlags
+      get_needed_update_flags() const override;
+
+    private:
+      const bool do_schlieren_plot;
+    };
+  };
+
+
+
+  template <int dim>
+  EulerProblem<dim>::Postprocessor::Postprocessor()
+    : do_schlieren_plot(dim == 2)
+  {}
+
+
+
+  // For the main evaluation of the field variables, we first check that the
+  // lengths of the arrays equal the expected values (the lengths `2*dim+4` or
+  // `2*dim+5` are derived from the sizes of the names we specify in the
+  // get_names() function below). Then we loop over all evaluation points and
+  // fill the respective information: First we fill the primal solution
+  // variables of density $\rho$, momentum $\rho \mathbf{u}$ and energy $E$,
+  // then we compute the derived velocity $\mathbf u$, the pressure $p$, the
+  // speed of sound $c=\sqrt{\gamma p / \rho}$, as well as the Schlieren plot
+  // showing $s = |\nabla \rho|^2$ in case it is enabled. (See step-69 for
+  // another example where we create a Schlieren plot.)
+  template <int dim>
+  void
+  EulerProblem<dim>::Postprocessor::evaluate_vector_field(
+    const DataPostprocessorInputs::Vector<dim> &inputs,
+    std::vector<Vector<double>> &               computed_quantities) const
+  {
+    const unsigned int n_evaluation_points = inputs.solution_values.size();
+
+    if (do_schlieren_plot == true)
+      Assert(inputs.solution_gradients.size() == n_evaluation_points,
+             ExcInternalError());
+
+    Assert(computed_quantities.size() == n_evaluation_points,
+           ExcInternalError());
+    Assert(inputs.solution_values[0].size() == dim + 2, ExcInternalError());
+    Assert(computed_quantities[0].size() ==
+             dim + 2 + (do_schlieren_plot == true ? 1 : 0),
+           ExcInternalError());
+
+    for (unsigned int p = 0; p < n_evaluation_points; ++p)
+      {
+        Tensor<1, dim + 2> solution;
+        for (unsigned int d = 0; d < dim + 2; ++d)
+          solution[d] = inputs.solution_values[p](d);
+
+        const double         density  = solution[0];
+        const Tensor<1, dim> velocity = euler_velocity<dim>(solution);
+        const double         pressure = euler_pressure<dim>(solution);
+
+        for (unsigned int d = 0; d < dim; ++d)
+          computed_quantities[p](d) = velocity[d];
+        computed_quantities[p](dim)     = pressure;
+        computed_quantities[p](dim + 1) = std::sqrt(gamma * pressure / density);
+
+        if (do_schlieren_plot == true)
+          computed_quantities[p](dim + 2) =
+            inputs.solution_gradients[p][0] * inputs.solution_gradients[p][0];
+      }
+  }
+
+
+
+  template <int dim>
+  std::vector<std::string>
+  EulerProblem<dim>::Postprocessor::get_names() const
+  {
+    std::vector<std::string> names;
+    for (unsigned int d = 0; d < dim; ++d)
+      names.emplace_back("velocity");
+    names.emplace_back("pressure");
+    names.emplace_back("speed_of_sound");
+
+    if (do_schlieren_plot == true)
+      names.emplace_back("schlieren_plot");
+
+    return names;
+  }
+
+
+
+  // For the interpretation of quantities, we have scalar density, energy,
+  // pressure, speed of sound, and the Schlieren plot, and vectors for the
+  // momentum and the velocity.
+  template <int dim>
+  std::vector<DataComponentInterpretation::DataComponentInterpretation>
+  EulerProblem<dim>::Postprocessor::get_data_component_interpretation() const
+  {
+    std::vector<DataComponentInterpretation::DataComponentInterpretation>
+      interpretation;
+    for (unsigned int d = 0; d < dim; ++d)
+      interpretation.push_back(
+        DataComponentInterpretation::component_is_part_of_vector);
+    interpretation.push_back(DataComponentInterpretation::component_is_scalar);
+    interpretation.push_back(DataComponentInterpretation::component_is_scalar);
+
+    if (do_schlieren_plot == true)
+      interpretation.push_back(
+        DataComponentInterpretation::component_is_scalar);
+
+    return interpretation;
+  }
+
+
+
+  // With respect to the necessary update flags, we only need the values for
+  // all quantities but the Schlieren plot, which is based on the density
+  // gradient.
+  template <int dim>
+  UpdateFlags
+  EulerProblem<dim>::Postprocessor::get_needed_update_flags() const
+  {
+    if (do_schlieren_plot == true)
+      return update_values | update_gradients;
+    else
+      return update_values;
+  }
+
+
+
+  // The constructor for this class is unsurprising: We set up a parallel
+  // triangulation based on the `MPI_COMM_WORLD` communicator, a vector finite
+  // element with `dim+2` components for density, momentum, and energy, a
+  // high-order mapping of the same degree as the underlying finite element,
+  // and initialize the time and time step to zero.
+  template <int dim>
+  EulerProblem<dim>::EulerProblem()
+    : pcout(deallog.get_file_stream(),
+            Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+#ifdef DEAL_II_WITH_P4EST
+    , triangulation(MPI_COMM_WORLD)
+#endif
+    , fe(FE_DGQ<dim>(fe_degree), dim + 2)
+    , mapping(fe_degree)
+    , dof_handler(triangulation)
+    , timer(pcout, TimerOutput::never, TimerOutput::wall_times)
+    , euler_operator(timer)
+    , time(0)
+    , time_step(0)
+  {}
+
+
+
+  // As a mesh, this tutorial program implements two options, depending on the
+  // global variable `testcase`: For the analytical variant (`testcase==0`),
+  // the domain is $(0, 10) \times (-5, 5)$, with Dirichlet boundary
+  // conditions (inflow) all around the domain. For `testcase==1`, we set the
+  // domain to a cylinder in a rectangular box, derived from the flow past
+  // cylinder testcase for incompressible viscous flow by Sch&auml;fer and
+  // Turek (1996). Here, we have a larger variety of boundaries. The inflow
+  // part at the left of the channel is given the inflow type, for which we
+  // choose a constant inflow profile, whereas we set a subsonic outflow at
+  // the right. For the boundary around the cylinder (boundary id equal to 2)
+  // as well as the channel walls (boundary id equal to 3) we use the wall
+  // boundary type, which is no-normal flow. Furthermore, for the 3d cylinder
+  // we also add a gravity force in vertical direction. Having the base mesh
+  // in place (including the manifolds set by
+  // GridGenerator::channel_with_cylinder()), we can then perform the
+  // specified number of global refinements, create the unknown numbering from
+  // the DoFHandler, and hand the DoFHandler and Mapping objects to the
+  // initialization of the EulerOperator.
+  template <int dim>
+  void
+  EulerProblem<dim>::make_grid_and_dofs()
+  {
+    switch (testcase)
+      {
+        case 0:
+          {
+            Point<dim> lower_left;
+            for (unsigned int d = 1; d < dim; ++d)
+              lower_left[d] = -5;
+
+            Point<dim> upper_right;
+            upper_right[0] = 10;
+            for (unsigned int d = 1; d < dim; ++d)
+              upper_right[d] = 5;
+
+            GridGenerator::hyper_rectangle(triangulation,
+                                           lower_left,
+                                           upper_right);
+            triangulation.refine_global(2);
+
+            euler_operator.set_inflow_boundary(
+              0, std::make_unique<ExactSolution<dim>>(0));
+
+            break;
+          }
+
+        case 1:
+          {
+            GridGenerator::channel_with_cylinder(
+              triangulation, 0.03, 1, 0, true);
+
+            euler_operator.set_inflow_boundary(
+              0, std::make_unique<ExactSolution<dim>>(0));
+            euler_operator.set_subsonic_outflow_boundary(
+              1, std::make_unique<ExactSolution<dim>>(0));
+
+            euler_operator.set_wall_boundary(2);
+            euler_operator.set_wall_boundary(3);
+
+            if (dim == 3)
+              euler_operator.set_body_force(
+                std::make_unique<Functions::ConstantFunction<dim>>(
+                  std::vector<double>({0., 0., -0.2})));
+
+            break;
+          }
+
+        default:
+          Assert(false, ExcNotImplemented());
+      }
+
+    triangulation.refine_global(n_global_refinements);
+
+    dof_handler.distribute_dofs(fe);
+
+    euler_operator.reinit(mapping, dof_handler);
+    euler_operator.initialize_vector(solution);
+
+    // In the following, we output some statistics about the problem. Because we
+    // often end up with quite large numbers of cells or degrees of freedom, we
+    // would like to print them with a comma to separate each set of three
+    // digits. This can be done via "locales", although the way this works is
+    // not particularly intuitive. step-32 explains this in slightly more
+    // detail.
+    std::locale s = pcout.get_stream().getloc();
+    pcout.get_stream().imbue(std::locale(""));
+    pcout << "Number of degrees of freedom: " << dof_handler.n_dofs()
+          << " ( = " << (dim + 2) << " [vars] x "
+          << triangulation.n_global_active_cells() << " [cells] x "
+          << Utilities::pow(fe_degree + 1, dim) << " [dofs/cell/var] )"
+          << std::endl;
+    pcout.get_stream().imbue(s);
+  }
+
+
+
+  // For output, we first let the Euler operator compute the errors of the
+  // numerical results. More precisely, we compute the error against the
+  // analytical result for the analytical solution case, whereas we compute
+  // the deviation against the background field with constant density and
+  // energy and constant velocity in $x$ direction for the second test case.
+  //
+  // The next step is to create output. This is similar to what is done in
+  // step-33: We let the postprocessor defined above control most of the
+  // output, except for the primal field that we write directly. For the
+  // analytical solution test case, we also perform another projection of the
+  // analytical solution and print the difference between that field and the
+  // numerical solution. Once we have defined all quantities to be written, we
+  // build the patches for output. Similarly to step-65, we create a
+  // high-order VTK output by setting the appropriate flag, which enables us
+  // to visualize fields of high polynomial degrees. Finally, we call the
+  // `DataOutInterface::write_vtu_in_parallel()` function to write the result
+  // to the given file name. This function uses special MPI parallel write
+  // facilities, which are typically more optimized for parallel file systems
+  // than the standard library's `std::ofstream` variants used in most other
+  // tutorial programs. A particularly nice feature of the
+  // `write_vtu_in_parallel()` function is the fact that it can combine output
+  // from all MPI ranks into a single file, making it unnecessary to have a
+  // central record of all such files (namely, the "pvtu" file).
+  //
+  // For parallel programs, it is often instructive to look at the partitioning
+  // of cells among processors. To this end, one can pass a vector of numbers
+  // to DataOut::add_data_vector() that contains as many entries as the
+  // current processor has active cells; these numbers should then be the
+  // rank of the processor that owns each of these cells. Such a vector
+  // could, for example, be obtained from
+  // GridTools::get_subdomain_association(). On the other hand, on each MPI
+  // process, DataOut will only read those entries that correspond to locally
+  // owned cells, and these of course all have the same value: namely, the rank
+  // of the current process. What is in the remaining entries of the vector
+  // doesn't actually matter, and so we can just get away with a cheap trick: We
+  // just fill *all* values of the vector we give to DataOut::add_data_vector()
+  // with the rank of the current MPI process. The key is that on each process,
+  // only the entries corresponding to the locally owned cells will be read,
+  // ignoring the (wrong) values in other entries. The fact that every process
+  // submits a vector in which the correct subset of entries is correct is all
+  // that is necessary.
+  template <int dim>
+  void
+  EulerProblem<dim>::output_results(const unsigned int result_number)
+  {
+    return;
+
+    const std::array<double, 3> errors =
+      euler_operator.compute_errors(ExactSolution<dim>(time), solution);
+    const std::string quantity_name = testcase == 0 ? "error" : "norm";
+
+    pcout << "Time:" << std::setw(8) << std::setprecision(3) << time
+          << ", dt: " << std::setw(8) << std::setprecision(2) << time_step
+          << ", " << quantity_name << " rho: " << std::setprecision(4)
+          << std::setw(10) << errors[0] << ", rho * u: " << std::setprecision(4)
+          << std::setw(10) << errors[1] << ", energy:" << std::setprecision(4)
+          << std::setw(10) << errors[2] << std::endl;
+
+    {
+      TimerOutput::Scope t(timer, "output");
+
+      Postprocessor postprocessor;
+      DataOut<dim>  data_out;
+
+      DataOutBase::VtkFlags flags;
+      flags.write_higher_order_cells = true;
+      data_out.set_flags(flags);
+
+      data_out.attach_dof_handler(dof_handler);
+      {
+        std::vector<std::string> names;
+        names.emplace_back("density");
+        for (unsigned int d = 0; d < dim; ++d)
+          names.emplace_back("momentum");
+        names.emplace_back("energy");
+
+        std::vector<DataComponentInterpretation::DataComponentInterpretation>
+          interpretation;
+        interpretation.push_back(
+          DataComponentInterpretation::component_is_scalar);
+        for (unsigned int d = 0; d < dim; ++d)
+          interpretation.push_back(
+            DataComponentInterpretation::component_is_part_of_vector);
+        interpretation.push_back(
+          DataComponentInterpretation::component_is_scalar);
+
+        data_out.add_data_vector(dof_handler, solution, names, interpretation);
+      }
+      data_out.add_data_vector(solution, postprocessor);
+
+      LinearAlgebra::distributed::Vector<Number> reference;
+      if (testcase == 0 && dim == 2)
+        {
+          reference.reinit(solution);
+          euler_operator.project(ExactSolution<dim>(time), reference);
+          reference.sadd(-1., 1, solution);
+          std::vector<std::string> names;
+          names.emplace_back("error_density");
+          for (unsigned int d = 0; d < dim; ++d)
+            names.emplace_back("error_momentum");
+          names.emplace_back("error_energy");
+
+          std::vector<DataComponentInterpretation::DataComponentInterpretation>
+            interpretation;
+          interpretation.push_back(
+            DataComponentInterpretation::component_is_scalar);
+          for (unsigned int d = 0; d < dim; ++d)
+            interpretation.push_back(
+              DataComponentInterpretation::component_is_part_of_vector);
+          interpretation.push_back(
+            DataComponentInterpretation::component_is_scalar);
+
+          data_out.add_data_vector(dof_handler,
+                                   reference,
+                                   names,
+                                   interpretation);
+        }
+
+      Vector<double> mpi_owner(triangulation.n_active_cells());
+      mpi_owner = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+      data_out.add_data_vector(mpi_owner, "owner");
+
+      data_out.build_patches(mapping,
+                             fe.degree,
+                             DataOut<dim>::curved_inner_cells);
+
+      const std::string filename =
+        "solution_" + Utilities::int_to_string(result_number, 3) + ".vtu";
+      data_out.write_vtu_in_parallel(filename, MPI_COMM_WORLD);
+    }
+  }
+
+
+
+  // The EulerProblem::run() function puts all pieces together. It starts off
+  // by calling the function that creates the mesh and sets up data structures,
+  // and then initializing the time integrator and the two temporary vectors of
+  // the low-storage integrator. We call these vectors `rk_register_1` and
+  // `rk_register_2`, and use the first vector to represent the quantity
+  // $\mathbf{r}_i$ and the second one for $\mathbf{k}_i$ in the formulas for
+  // the Runge--Kutta scheme outlined in the introduction. Before we start the
+  // time loop, we compute the time step size by the
+  // `EulerOperator::compute_cell_transport_speed()` function. For reasons of
+  // comparison, we compare the result obtained there with the minimal mesh
+  // size and print them to screen. For velocities and speeds of sound close
+  // to unity as in this tutorial program, the predicted effective mesh size
+  // will be close, but they could vary if scaling were different.
+  template <int dim>
+  void
+  EulerProblem<dim>::run()
+  {
+    {
+      const unsigned int n_vect_number = VectorizedArray<Number>::size();
+      const unsigned int n_vect_bits   = 8 * sizeof(Number) * n_vect_number;
+
+      pcout << "Running with "
+            << Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD)
+            << " MPI processes" << std::endl;
+      // pcout << "Vectorization over " << n_vect_number << ' '
+      //       << (std::is_same<Number, double>::value ? "doubles" : "floats")
+      //       << " = " << n_vect_bits << " bits ("
+      //       << Utilities::System::get_current_vectorization_level() << ')'
+      //       << std::endl;
+    }
+
+    make_grid_and_dofs();
+
+    const LowStorageRungeKuttaIntegrator integrator(lsrk_scheme);
+
+    LinearAlgebra::distributed::Vector<Number> rk_register_1;
+    LinearAlgebra::distributed::Vector<Number> rk_register_2;
+    rk_register_1.reinit(solution);
+    rk_register_2.reinit(solution);
+
+    euler_operator.project(ExactSolution<dim>(time), solution);
+
+    double min_vertex_distance = std::numeric_limits<double>::max();
+    for (const auto &cell : triangulation.active_cell_iterators())
+      if (cell->is_locally_owned())
+        min_vertex_distance =
+          std::min(min_vertex_distance, cell->minimum_vertex_distance());
+    min_vertex_distance =
+      Utilities::MPI::min(min_vertex_distance, MPI_COMM_WORLD);
+
+    time_step = courant_number * integrator.n_stages() /
+                euler_operator.compute_cell_transport_speed(solution);
+    pcout << "Time step size: " << time_step
+          << ", minimal h: " << min_vertex_distance
+          << ", initial transport scaling: "
+          << 1. / euler_operator.compute_cell_transport_speed(solution)
+          << std::endl
+          << std::endl;
+
+    output_results(0);
+
+    // Now we are ready to start the time loop, which we run until the time
+    // has reached the desired end time. Every 5 time steps, we compute a new
+    // estimate for the time step -- since the solution is nonlinear, it is
+    // most effective to adapt the value during the course of the
+    // simulation. In case the Courant number was chosen too aggressively, the
+    // simulation will typically blow up with time step NaN, so that is easy
+    // to detect here. One thing to note is that roundoff errors might
+    // propagate to the leading digits due to an interaction of slightly
+    // different time step selections that in turn lead to slightly different
+    // solutions. To decrease this sensitivity, it is common practice to round
+    // or truncate the time step size to a few digits, e.g. 3 in this case. In
+    // case the current time is near the prescribed 'tick' value for output
+    // (e.g. 0.02), we also write the output. After the end of the time loop,
+    // we summarize the computation by printing some statistics, which is
+    // mostly done by the TimerOutput::print_wall_time_statistics() function.
+    unsigned int timestep_number = 0;
+
+    // Here we detail how to use PETSc TS explicit solvers of Runge Kutta type.
+    // We use two small helper functions to copy to and from
+    // instances of `LinearAlgebra::distributed::Vector`, since the
+    // current implementation does not directly support this vector class.
+    using VectorType = PETScWrappers::MPI::Vector;
+
+    auto copy_from_petsc = [](const VectorType &                          src,
+                              LinearAlgebra::distributed::Vector<Number> &dst) {
+      auto lr = src.local_range();
+      for (unsigned int i = lr.first; i < lr.second; i++)
+        dst[i] = src[i];
+    };
+
+    auto copy_to_petsc =
+      [](const LinearAlgebra::distributed::Vector<Number> &src,
+         VectorType &                                      dst) {
+        auto lr = dst.local_range();
+        for (unsigned int i = lr.first; i < lr.second; i++)
+          dst[i] = src[i];
+        dst.compress(VectorOperation::insert);
+      };
+
+    // We construct our solver with parameters like initial time step and final
+    // time, together with the solver type for Runge-Kutta. We also decide to
+    // use adaptive time stepping based on <a
+    // href="https://dl.acm.org/doi/10.1145/641876.641877"> digital signal
+    // processing</a>, and set adaptive tolerances to -1 to use default
+    // tolerances from PETSc
+    PETScWrappers::TimeStepperData tsdata;
+    tsdata.final_time         = final_time;
+    tsdata.initial_step_size  = time_step;
+    tsdata.ts_type            = "rk";
+    tsdata.ts_adapt_type      = "dsp";
+    tsdata.absolute_tolerance = -1.0;
+    tsdata.relative_tolerance = -1.0;
+
+    PETScWrappers::TimeStepper<VectorType> petsc_time_stepper(tsdata);
+
+    // PETScWrappers::TimeStepper can solve ODE/DAE problems in explicit form
+    // i.e., $\dot{u} = G(u,t)$, or in fully implicit form $F(\dot{u},u,t) = 0$.
+    // This function implements the explicit interface for $G$.
+    petsc_time_stepper.explicit_function =
+      [&](const double t, const VectorType &y, VectorType &f) {
+        time = t;
+        copy_from_petsc(y, rk_register_1);
+        euler_operator.apply(t, rk_register_1, rk_register_2);
+        copy_to_petsc(rk_register_2, f);
+        return 0;
+      };
+
+    // We can also attach a monitoring callback to inspect the solution.
+    petsc_time_stepper.monitor =
+      [&](const double t, const VectorType &y, const unsigned int) {
+        time      = t;
+        time_step = petsc_time_stepper.get_time_step();
+        copy_from_petsc(y, solution);
+
+        static int prev_tick = 0;
+        int        curr_tick = static_cast<int>(t / output_tick);
+        if (curr_tick != prev_tick)
+          {
+            pcout << "Monitor is at t=" << t << " with ||y||=" << y.l2_norm()
+                  << std::endl;
+
+            output_results(
+              static_cast<unsigned int>(std::round(t / output_tick)));
+          }
+
+        prev_tick = curr_tick;
+        return 0;
+      };
+
+    // We finally set the initial condition and integrate the ODE.
+    VectorType psolution(solution.locally_owned_elements(),
+                         solution.get_mpi_communicator());
+    copy_to_petsc(solution, psolution);
+    petsc_time_stepper.solve(psolution);
+    copy_from_petsc(psolution, solution);
+
+    pcout << std::endl;
+  }
+
+} // namespace Euler_DG
+
+
+
+// The main() function is not surprising and follows what was done in all
+// previous MPI programs: As we run an MPI program, we need to call `MPI_Init()`
+// and `MPI_Finalize()`, which we do through the
+// Utilities::MPI::MPI_InitFinalize data structure. Note that we run the program
+// only with MPI, and set the thread count to 1.
+int
+main(int argc, char **argv)
+{
+  using namespace Euler_DG;
+  using namespace dealii;
+
+  initlog();
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  try
+    {
+      deallog.depth_console(0);
+
+      EulerProblem<dimension> euler_problem;
+      euler_problem.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;
+    }
+
+  return 0;
+}
diff --git a/tests/petsc/step-67-with-petsc-ts.with_petsc=true.with_p4est=true.mpirun=1.output b/tests/petsc/step-67-with-petsc-ts.with_petsc=true.with_p4est=true.mpirun=1.output
new file mode 100644 (file)
index 0000000..f45824c
--- /dev/null
@@ -0,0 +1,26 @@
+
+Running with 1 MPI processes
+Number of degrees of freedom: 9216 ( = 4 [vars] x 64 [cells] x 36 [dofs/cell/var] )
+Time step size: 0.0275808, minimal h: 1.25000, initial transport scaling: 0.411150
+
+Monitor is at t=0.0693617 with ||y||=158.219
+Monitor is at t=0.116085 with ||y||=158.220
+Monitor is at t=0.170333 with ||y||=158.221
+Monitor is at t=0.223603 with ||y||=158.221
+Monitor is at t=0.259770 with ||y||=158.221
+Monitor is at t=0.332301 with ||y||=158.220
+Monitor is at t=0.373992 with ||y||=158.219
+Monitor is at t=0.417647 with ||y||=158.218
+Monitor is at t=0.462060 with ||y||=158.217
+Monitor is at t=0.503602 with ||y||=158.216
+Monitor is at t=0.578443 with ||y||=158.214
+Monitor is at t=0.618551 with ||y||=158.214
+Monitor is at t=0.660433 with ||y||=158.214
+Monitor is at t=0.703591 with ||y||=158.214
+Monitor is at t=0.787191 with ||y||=158.215
+Monitor is at t=0.825503 with ||y||=158.216
+Monitor is at t=0.864442 with ||y||=158.217
+Monitor is at t=0.905591 with ||y||=158.217
+Monitor is at t=0.989302 with ||y||=158.219
+Monitor is at t=1.03147 with ||y||=158.219
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.