/* ---------------------------------------------------------------------
*
- * Copyright (C) 2019 by the deal.II authors
+ * Copyright (C) 2020 by the deal.II authors
*
* This file is part of the deal.II library.
*
*
* ---------------------------------------------------------------------
+ *
+ * Author: Martin Kronbichler, 2020
*/
-// @sect3{Include files}
+// The include files are similar to the previous matrix-free tutorial programs
+// step-37, step-48, and step-59
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.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/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.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>
-#include <deal.II/base/mpi.h>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+// This includes the CellwiseInverseMassMatrix data structure for the
+// definition of the interface to mass matrix inversion, the only new include
+// file for this tutorial program.
+#include <deal.II/matrix_free/operators.h>
-namespace step67
+
+
+namespace Euler_DG
{
using namespace dealii;
-}
-// @sect3{The main function}
+ // 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 = 3;
+ 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 = testcase == 0 ? 10 : 2.0;
+ constexpr double output_tick = testcase == 0 ? 1 : 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 select the function with 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. We note that this formula might lose accuracy in the last
+ // digits for very small numbers compared to `std::pow()`, 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
+ {
+ double t = this->get_time();
+ if (testcase == 0)
+ {
+ Assert(dim == 2, ExcNotImplemented());
+ double beta = 5;
+ Point<dim> x0;
+ x0[0] = 5.;
+ double radius_sqr =
+ (x - x0).norm_square() - 2. * (x[0] - x0[0]) * t + t * t;
+ 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);
+ }
+ }
+ else
+ {
+ if (component == 0)
+ return 1.;
+ else if (component == 1)
+ return 0.4;
+ else if (component == dim + 1)
+ return 3.097857142857143;
+ else
+ 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 tableaus 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{T}_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)
+ {
+ // 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.
+ if (scheme == stage_3_order_3)
+ {
+ bi = {{0.245170287303492, 0.184896052186740, 0.569933660509768}};
+ ai = {{0.755726351946097, 0.386954477304099}};
+ }
+ // The next scheme is a five-stage scheme of order four, again defined
+ // in the paper by Kennedy et al. (2000).
+ else if (scheme == stage_5_order_4)
+ {
+ bi = {{1153189308089. / 22510343858157.,
+ 1772645290293. / 4653164025191.,
+ -1672844663538. / 4480602732383.,
+ 2114624349019. / 3568978502595.,
+ 5198255086312. / 14908931495163.}};
+ ai = {{970286171893. / 4311952581923.,
+ 6584761158862. / 12103376702013.,
+ 2251764453980. / 15575788980749.,
+ 26877169314380. / 34165994151039.}};
+ }
+ // This 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
+ // disspiation 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.
+ else if (scheme == stage_7_order_4)
+ {
+ bi = {{0.0941840925477795334,
+ 0.149683694803496998,
+ 0.285204742060440058,
+ -0.122201846148053668,
+ 0.0605151571191401122,
+ 0.345986987898399296,
+ 0.186627171718797670}};
+ ai = {{0.241566650129646868 + bi[0],
+ 0.0423866513027719953 + bi[1],
+ 0.215602732678803776 + bi[2],
+ 0.232328007537583987 + bi[3],
+ 0.256223412574146438 + bi[4],
+ 0.0978694102142697230 + bi[5]}};
+ }
+ // 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.
+ else if (scheme == stage_9_order_5)
+ {
+ bi = {{2274579626619. / 23610510767302.,
+ 693987741272. / 12394497460941.,
+ -347131529483. / 15096185902911.,
+ 1144057200723. / 32081666971178.,
+ 1562491064753. / 11797114684756.,
+ 13113619727965. / 44346030145118.,
+ 393957816125. / 7825732611452.,
+ 720647959663. / 6565743875477.,
+ 3559252274877. / 14424734981077.}};
+ ai = {{1107026461565. / 5417078080134.,
+ 38141181049399. / 41724347789894.,
+ 493273079041. / 11940823631197.,
+ 1851571280403. / 6147804934346.,
+ 11782306865191. / 62590030070788.,
+ 9452544825720. / 13648368537481.,
+ 4435885630781. / 26285702406235.,
+ 2357909744247. / 11371140753790.}};
+ }
+ else
+ AssertThrow(false, ExcNotImplemented());
+ }
+
+ 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{T}_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: Here, we evaluate the solution from
+ // the old solution $\mathbf{w}^n$ rather than a $\mathbf T_i$ vector, so
+ // the first argument is `solution`. We here let the stage vector
+ // $\mathbf{T}_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{T}_{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_Ti,
+ VectorType & vec_Ki)
+ {
+ AssertDimension(ai.size() + 1, bi.size());
+ pde_operator.perform_stage(current_time,
+ bi[0] * time_step,
+ ai[0] * time_step,
+ solution,
+ vec_Ti,
+ solution,
+ vec_Ti);
+ double sum_previous_bi = 0;
+ for (unsigned int stage = 1; stage < bi.size(); ++stage)
+ {
+ const double c_i = sum_previous_bi + ai[stage - 1];
+ pde_operator.perform_stage(current_time + c_i * time_step,
+ bi[stage] * time_step,
+ (stage == bi.size() - 1 ?
+ 0 :
+ ai[stage] * time_step),
+ vec_Ti,
+ vec_Ki,
+ solution,
+ vec_Ti);
+ sum_previous_bi += bi[stage - 1];
+ }
+ }
+
+ private:
+ std::vector<double> bi;
+ std::vector<double> ai;
+ };
+
+
+
+ // @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 inline to where they are called. This is critical for
+ // performance because we repeatedly call into some of those functions: For
+ // example, we both use the velocity for the computation of the flux further
+ // down, but also for the computation of the pressure. Keeping these
+ // functions inline means that the repeated use is seen by the compiler
+ // during the optimization passes, and it eventually only keeps a single
+ // call around. If it were a separate function, it gets more complicated or
+ // impossible because already computed temporary information cannot be
+ // passed around.
+ //
+ // 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 do it 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_u = conserved_variables[1] * velocity[0];
+ for (unsigned int d = 1; d < dim; ++d)
+ rho_u_u += conserved_variables[1 + d] * velocity[d];
+
+ return (gamma - 1.) * (conserved_variables[dim + 1] - 0.5 * rho_u_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 the
+ // basic scheme is already very powerful.
+ //
+ // 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{\lambda 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.
+ 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);
+
+ if (numerical_flux_type == 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);
+ }
+ else if (numerical_flux_type == 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));
+ }
+ }
+
+
+
+ // 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,
+ const std::shared_ptr<Function<dim>> &inflow_function);
+
+ void set_subsonic_outflow_boundary(
+ const types::boundary_id boundary_id,
+ const std::shared_ptr<Function<dim>> &outflow_energy);
+
+ void set_wall_boundary(const types::boundary_id boundary_id);
+
+ void set_body_force(const std::shared_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> ¤t_Ti,
+ LinearAlgebra::distributed::Vector<Number> & vec_Ki,
+ LinearAlgebra::distributed::Vector<Number> & solution,
+ LinearAlgebra::distributed::Vector<Number> &next_Ti) const;
+
+ void project(const Function<dim> & function,
+ LinearAlgebra::distributed::Vector<Number> &solution) const;
+
+ Tensor<1, 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::shared_ptr<Function<dim>>>
+ inflow_boundaries;
+ std::map<types::boundary_id, std::shared_ptr<Function<dim>>>
+ subsonic_outflow_boundaries;
+ std::set<types::boundary_id> wall_boundaries;
+ std::shared_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> & cell_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> & cell_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)
+ {
+ std::vector<const DoFHandler<dim> *> dof_handlers({&dof_handler});
+ AffineConstraints<double> dummy;
+ std::vector<const AffineConstraints<double> *> constraints({&dummy});
+ 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 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 assign a
+ // boundary to both an inflow and say a subsonic outflow.
+ 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,
+ const std::shared_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] = 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,
+ const std::shared_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] = 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(
+ const std::shared_ptr<Function<dim>> &body_force)
+ {
+ AssertDimension(body_force->n_components, dim);
+ this->body_force = 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 halfs 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 the `euler_flux` function
+ // given the current solution interpolated 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().
+ template <int dim, int degree, int n_points_1d>
+ void EulerOperator<dim, degree, n_points_1d>::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
+ {
+ 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, true, false);
+
+ 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, true, 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, true, false);
+
+ phi_m.reinit(face);
+ phi_m.gather_evaluate(src, true, false);
+
+ 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(true, false, dst);
+ phi_m.integrate_scatter(true, false, 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. The discontinuous Galerkin method sets these values
+ // weakly, so 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.
+ //
+ // 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 value to the interior value 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 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. 28
+ // (2012). Here, the momentum term needs to be added once again, which
+ // translates 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 implement the 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, 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, true, false);
+
+ 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?"));
+
+ 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(true, false, dst);
+ }
+ }
+
+
+
+ // This 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 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 vecor 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> & data,
+ 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);
+ }
+ }
+
+
-// Finally, the main function. There isn't much to do here, only to call the
-// two subfunctions, which produce the two grids.
+ // @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 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);
+ }
+ }
+
+
+
+ // This function implements EulerOperator::apply() followed by some updates
+ // to the vectors, namely `next_Ti = 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 performed 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. We add the `DEAL_II_OPENMP_SIMD_PRAGMA` before the local loop
+ // to suggest the compiler to SIMD parallelize this loop (which means in
+ // practice that we ensure that there is no overlapping, also called
+ // aliasing, between the index ranges of the pointers we use inside the
+ // loops). Note that we select a different code path for the last
+ // Runge--Kutta stage when we do not need to update the `next_Ti`
+ // 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> ¤t_Ti,
+ LinearAlgebra::distributed::Vector<Number> & vec_Ki,
+ LinearAlgebra::distributed::Vector<Number> & solution,
+ LinearAlgebra::distributed::Vector<Number> & next_Ti) 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_Ti,
+ 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_Ti,
+ 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_Ti.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_Ti.local_element(i);
+ const Number sol_i = solution.local_element(i);
+ solution.local_element(i) = sol_i + bi * K_i;
+ next_Ti.local_element(i) = sol_i + ai * K_i;
+ }
+ }
+ });
+ }
+ }
+
+
+
+ // This function 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 $\Omega_e$ is an operation of the form $\underbrace{S J^e S^\mathrm
+ // T}_{\mathcal M^e} \mathbf{w}^e = S J^e
+ // \tilde{\mathbf{w}}(\mathbf{x}_q)_{q=1:n_q}$, where $J^e$ is the diagonal
+ // matrix containing the determinant of the Jacobian times the quadrature
+ // weight (JxW), $\mathcal M^e$ 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}^e = \left(S J^e S^\mathrm T\right)^{-1} S J^e
+ // \tilde{\mathbf{w}}(\mathbf{x}_q)_{q=1:n_q} = S^{-\mathrm T}
+ // \left(J^e\right)^{-1} S^{-1} S J^e
+ // \tilde{\mathbf{w}}(\mathbf{x}_q)_{q=1:n_q}$. Now, the term $S^{-1} S$ and
+ // then $\left(J^e\right)^{-1} J^e$ cancel, resulting in the final
+ // expression $\mathbf{w}^e = 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 nothing else
+ // than 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.
+ 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_ghosts();
+ for (unsigned int cell = 0; cell < data.n_macro_cells(); ++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>
+ Tensor<1, 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");
+ Tensor<1, 3> errors;
+ 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, true, false);
+ Tensor<1, 3, VectorizedArray<Number>> local_errors;
+ 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[0] += error[0] * error[0] * JxW;
+ for (unsigned int d = 0; d < dim; ++d)
+ local_errors[1] += error[d + 1] * error[d + 1] * JxW;
+ local_errors[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[d] += local_errors[d][v];
+ }
+ errors = Utilities::MPI::sum(errors, MPI_COMM_WORLD);
+ for (unsigned int d = 0; d < 3; ++d)
+ errors[d] = std::sqrt(errors[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 via $\mathbf{u}$ and the propagation of sound waves with sound
+ // speed $c = \sqrt{\gamma p/\rho}$. The former is scaled by the mesh size,
+ // so an estimate of the maximal velocity can be obtained by computing
+ // $\|J^{-\mathrm T} \mathbf{u}\|_\inf$, 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 are 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, true, false);
+ 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 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;
+ MappingQGeneric<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
+ // in case it is enabled.
+ 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());
+
+ if (do_schlieren_plot == true)
+ {
+ Assert(computed_quantities[0].size() == 2 * dim + 5,
+ ExcInternalError());
+ }
+ else
+ {
+ Assert(computed_quantities[0].size() == 2 * dim + 4,
+ ExcInternalError());
+ }
+
+ for (unsigned int q = 0; q < n_evaluation_points; ++q)
+ {
+ Tensor<1, dim + 2> solution;
+ for (unsigned int d = 0; d < dim + 2; ++d)
+ solution[d] = inputs.solution_values[q](d);
+ for (unsigned int d = 0; d < dim + 2; ++d)
+ computed_quantities[q](d) = solution[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[q](dim + 2 + d) = velocity[d];
+ computed_quantities[q](2 * dim + 2) = pressure;
+ computed_quantities[q](2 * dim + 3) =
+ std::sqrt(gamma * pressure / density);
+ if (do_schlieren_plot == true)
+ computed_quantities[q](2 * dim + 4) =
+ inputs.solution_gradients[q][0] * inputs.solution_gradients[q][0];
+ }
+ }
+
+
+
+ template <int dim>
+ std::vector<std::string> EulerProblem<dim>::Postprocessor::get_names() const
+ {
+ 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");
+ 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;
+ 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);
+ 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(std::cout, 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 according to 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ä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 cyclinder
+ // 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()
+ {
+ if (testcase == 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;
+
+ std::vector<unsigned int> refinements(dim, 1);
+ GridGenerator::subdivided_hyper_rectangle(triangulation,
+ refinements,
+ lower_left,
+ upper_right);
+ triangulation.refine_global(2);
+
+ euler_operator.set_inflow_boundary(
+ 0, std::make_shared<ExactSolution<dim>>(0));
+ }
+ else if (testcase == 1)
+ {
+ GridGenerator::channel_with_cylinder(triangulation, 0.03, 1, 0, true);
+ euler_operator.set_inflow_boundary(
+ 0, std::make_shared<ExactSolution<dim>>(0));
+ euler_operator.set_subsonic_outflow_boundary(
+ 1, std::make_shared<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_shared<Functions::ConstantFunction<dim>>(
+ std::vector<double>({0., 0., -0.2})));
+ }
+
+ triangulation.refine_global(n_global_refinements);
+
+ dof_handler.distribute_dofs(fe);
+
+ std::locale s = pcout.get_stream().getloc();
+ pcout.get_stream().imbue(std::locale("en_US.UTF-8"));
+ 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);
+
+ euler_operator.reinit(mapping, dof_handler);
+ euler_operator.initialize_vector(solution);
+ }
+
+
+
+ // 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. 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, obviating a VTU
+ // master file.
+ template <int dim>
+ void EulerProblem<dim>::output_results(const unsigned int result_number)
+ {
+ Tensor<1, 3> errors =
+ euler_operator.compute_errors(ExactSolution<dim>(time), solution);
+ 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");
+
+ 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);
+ Postprocessor postprocessor;
+ 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 of
+ // by calling the function that creates the mesh and sets up data structures
+ // and initializing the time integrator and the two temporary vectors of the
+ // low-storage integrator. 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()
+ << "), VECTORIZATION_LEVEL=" << DEAL_II_COMPILER_VECTORIZATION_LEVEL
+ << std::endl;
+ }
+
+ make_grid_and_dofs();
+
+ 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);
+
+ typename Triangulation<dim>::active_cell_iterator cell = triangulation
+ .begin_active(),
+ endc =
+ triangulation.end();
+ double min_vertex_distance = std::numeric_limits<double>::max();
+ for (; cell != endc; ++cell)
+ 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;
+
+ while (time < FINAL_TIME - 1e-12)
+ {
+ ++timestep_number;
+ if (timestep_number % 5 == 0)
+ time_step =
+ courant_number * integrator.n_stages() /
+ Utilities::truncate_to_n_digits(
+ euler_operator.compute_cell_transport_speed(solution), 3);
+
+ {
+ TimerOutput::Scope t(timer, "rk time step total");
+ integrator.perform_time_step(euler_operator,
+ time,
+ time_step,
+ solution,
+ rk_register_1,
+ rk_register_2);
+ }
+
+ time += time_step;
+
+ if (static_cast<int>(time / output_tick) !=
+ static_cast<int>((time - time_step) / output_tick) ||
+ time >= FINAL_TIME - 1e-12)
+ output_results(
+ static_cast<unsigned int>(std::round(time / output_tick)));
+ }
+
+ timer.print_wall_time_statistics(MPI_COMM_WORLD);
+ pcout << std::endl;
+ }
+
+} // namespace Euler_DG
+
+
+
+// The main() function is not surprising and follows what was done in step-59:
+// 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)
{
- dealii::Utilities::MPI::MPI_InitFinalize mpi_init_finalize(argc, argv, 1);
+ using namespace Euler_DG;
+ using namespace dealii;
+
+ 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;
}
projects / dealii.git / commitdiff