From 143c33239e921036f9681ed8736e637df08b1dab Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 16 Dec 2019 15:49:35 +0100 Subject: [PATCH] Add code for step-67 tutorial program --- examples/step-67/doc/builds-on | 2 +- examples/step-67/step-67.cc | 2188 +++++++++++++++++++++++++++++++- 2 files changed, 2180 insertions(+), 10 deletions(-) diff --git a/examples/step-67/doc/builds-on b/examples/step-67/doc/builds-on index 8943c29acb..d191181421 100644 --- a/examples/step-67/doc/builds-on +++ b/examples/step-67/doc/builds-on @@ -1 +1 @@ -step-48, step-59 +step-33, step-48, step-59 diff --git a/examples/step-67/step-67.cc b/examples/step-67/step-67.cc index d3e5722851..f3c0c533fe 100644 --- a/examples/step-67/step-67.cc +++ b/examples/step-67/step-67.cc @@ -1,6 +1,6 @@ /* --------------------------------------------------------------------- * - * 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. * @@ -13,23 +13,2193 @@ * * --------------------------------------------------------------------- + * + * 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 +#include +#include +#include +#include +#include + +#include + +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include -#include +#include +#include +#include +// 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 -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 + class ExactSolution : public Function + { + public: + ExactSolution(const double time) + : Function(dim + 2, time) + {} + + virtual double value(const Point & 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 + double ExactSolution::value(const Point & x, + const unsigned int component) const + { + double t = this->get_time(); + if (testcase == 0) + { + Assert(dim == 2, ExcNotImplemented()); + double beta = 5; + Point 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 + 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 bi; + std::vector 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 + 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 + inline DEAL_II_ALWAYS_INLINE // + Number + euler_pressure(const Tensor<1, dim + 2, Number> &conserved_variables) + { + const Tensor<1, dim, Number> velocity = + euler_velocity(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 + 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(conserved_variables); + const Number pressure = euler_pressure(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 + 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 + 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(u_m); + const auto velocity_p = euler_velocity(u_p); + + const auto pressure_m = euler_pressure(u_m); + const auto pressure_p = euler_pressure(u_p); + + const auto flux_m = euler_flux(u_m); + const auto flux_p = euler_flux(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 + VectorizedArray + evaluate_function(const Function & function, + const Point> &p_vectorized, + const unsigned int component) + { + VectorizedArray result; + for (unsigned int v = 0; v < VectorizedArray::size(); ++v) + { + Point p; + for (unsigned int d = 0; d < dim; ++d) + p[d] = p_vectorized[d][v]; + result[v] = function.value(p, component); + } + return result; + } + + template + Tensor<1, n_components, VectorizedArray> + evaluate_function(const Function & function, + const Point> &p_vectorized) + { + AssertDimension(function.n_components, n_components); + Tensor<1, n_components, VectorizedArray> result; + for (unsigned int v = 0; v < VectorizedArray::size(); ++v) + { + Point 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 + class EulerOperator + { + public: + static constexpr unsigned int n_quadrature_points_1d = n_points_1d; + + EulerOperator(TimerOutput &timer_output); + + void reinit(const Mapping & mapping, + const DoFHandler &dof_handler); + + void + set_inflow_boundary(const types::boundary_id boundary_id, + const std::shared_ptr> &inflow_function); + + void set_subsonic_outflow_boundary( + const types::boundary_id boundary_id, + const std::shared_ptr> &outflow_energy); + + void set_wall_boundary(const types::boundary_id boundary_id); + + void set_body_force(const std::shared_ptr> &body_force); + + void apply(const double current_time, + const LinearAlgebra::distributed::Vector &src, + LinearAlgebra::distributed::Vector & dst) const; + + void + perform_stage(const Number cur_time, + const Number factor_solution, + const Number factor_ai, + const LinearAlgebra::distributed::Vector ¤t_Ti, + LinearAlgebra::distributed::Vector & vec_Ki, + LinearAlgebra::distributed::Vector & solution, + LinearAlgebra::distributed::Vector &next_Ti) const; + + void project(const Function & function, + LinearAlgebra::distributed::Vector &solution) const; + + Tensor<1, 3> compute_errors( + const Function & function, + const LinearAlgebra::distributed::Vector &solution) const; + + double compute_cell_transport_speed( + const LinearAlgebra::distributed::Vector &solution) const; + + void + initialize_vector(LinearAlgebra::distributed::Vector &vector) const; + + private: + MatrixFree data; + + TimerOutput &timer; + + std::map>> + inflow_boundaries; + std::map>> + subsonic_outflow_boundaries; + std::set wall_boundaries; + std::shared_ptr> body_force; + + void local_apply_inverse_mass_matrix( + const MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & cell_range) const; + + void local_apply_cell( + const MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & cell_range) const; + + void local_apply_face( + const MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & cell_range) const; + + void local_apply_boundary_face( + const MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & cell_range) const; + }; + + + + template + EulerOperator::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 + void EulerOperator::reinit( + const Mapping & mapping, + const DoFHandler &dof_handler) + { + std::vector *> dof_handlers({&dof_handler}); + AffineConstraints dummy; + std::vector *> constraints({&dummy}); + std::vector> quadratures( + {QGauss<1>(n_q_points_1d), QGauss<1>(fe_degree + 1)}); + typename MatrixFree::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::AdditionalData::none; + + data.reinit( + mapping, dof_handlers, constraints, quadratures, additional_data); + } + + + + template + void EulerOperator::initialize_vector( + LinearAlgebra::distributed::Vector &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 + void EulerOperator::set_inflow_boundary( + const types::boundary_id boundary_id, + const std::shared_ptr> &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(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 + void EulerOperator::set_subsonic_outflow_boundary( + const types::boundary_id boundary_id, + const std::shared_ptr> &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(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 + void EulerOperator::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(boundary_id)) + + " to another type of boundary before now setting " + + "it as wall boundary")); + wall_boundaries.insert(boundary_id); + } + + template + void EulerOperator::set_body_force( + const std::shared_ptr> &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 + void EulerOperator::local_apply_cell( + const MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & cell_range) const + { + FEEvaluation phi(data); + Tensor<1, dim, VectorizedArray> constant_body_force; + const Functions::ConstantFunction * constant_function = + dynamic_cast *>(body_force.get()); + if (constant_function) + constant_body_force = evaluate_function( + *constant_function, Point>()); + + 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(w_q), q); + if (body_force.get() != nullptr) + { + const Tensor<1, dim, VectorizedArray> force = + constant_function ? constant_body_force : + evaluate_function( + *body_force, phi.quadrature_point(q)); + Tensor<1, dim + 2, VectorizedArray> 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 + void EulerOperator::local_apply_face( + const MatrixFree &, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & face_range) const + { + FEFaceEvaluation phi_m(data, + true); + FEFaceEvaluation 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(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 + // exterior 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 + // backflow, 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 + void EulerOperator::local_apply_boundary_face( + const MatrixFree &, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & face_range) const + { + FEFaceEvaluation 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> 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(w_m, w_p, normal); + + if (at_outflow) + for (unsigned int v = 0; v < VectorizedArray::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 + void EulerOperator::local_apply_inverse_mass_matrix( + const MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & cell_range) const + { + FEEvaluation phi(data, 0, 1); + MatrixFreeOperators::CellwiseInverseMassMatrix + 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 + void EulerOperator::apply( + const double current_time, + const LinearAlgebra::distributed::Vector &src, + LinearAlgebra::distributed::Vector & 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::DataAccessOnFaces::values, + MatrixFree::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 + void EulerOperator::perform_stage( + const Number current_time, + const Number factor_solution, + const Number factor_ai, + const LinearAlgebra::distributed::Vector ¤t_Ti, + LinearAlgebra::distributed::Vector & vec_Ki, + LinearAlgebra::distributed::Vector & solution, + LinearAlgebra::distributed::Vector & 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::DataAccessOnFaces::values, + MatrixFree::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(), + [&](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 + void EulerOperator::project( + const Function & function, + LinearAlgebra::distributed::Vector &solution) const + { + FEEvaluation phi(data, 0, 1); + MatrixFreeOperators::CellwiseInverseMassMatrix + 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 + // lane 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 + Tensor<1, 3> EulerOperator::compute_errors( + const Function & function, + const LinearAlgebra::distributed::Vector &solution) const + { + TimerOutput::Scope t(timer, "compute errors"); + Tensor<1, 3> errors; + FEEvaluation 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> 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 + double EulerOperator::compute_cell_transport_speed( + const LinearAlgebra::distributed::Vector &solution) const + { + TimerOutput::Scope t(timer, "compute transport speed"); + Number max_transport = 0; + FEEvaluation 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 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(solution); + const auto pressure = euler_pressure(solution); + + const auto inverse_jacobian = phi.inverse_jacobian(q); + const auto convective_speed = inverse_jacobian * velocity; + VectorizedArray 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> 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 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 + class EulerProblem + { + public: + EulerProblem(); + + void run(); + + private: + void make_grid_and_dofs(); + + void output_results(const unsigned int result_number); + + LinearAlgebra::distributed::Vector solution; + + ConditionalOStream pcout; + +#ifdef DEAL_II_WITH_P4EST + parallel::distributed::Triangulation triangulation; +#else + Triangulation triangulation; +#endif + + FESystem fe; + MappingQGeneric mapping; + DoFHandler dof_handler; + + TimerOutput timer; + + EulerOperator euler_operator; + + double time, time_step; + + class Postprocessor : public DataPostprocessor + { + public: + Postprocessor(); + + virtual void evaluate_vector_field( + const DataPostprocessorInputs::Vector &inputs, + std::vector> &computed_quantities) const override; + + virtual std::vector 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 + EulerProblem::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 + void EulerProblem::Postprocessor::evaluate_vector_field( + const DataPostprocessorInputs::Vector &inputs, + std::vector> & 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(solution); + const double pressure = euler_pressure(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 + std::vector EulerProblem::Postprocessor::get_names() const + { + std::vector 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 + std::vector + EulerProblem::Postprocessor::get_data_component_interpretation() const + { + std::vector + 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 + UpdateFlags EulerProblem::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 + EulerProblem::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(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 + void EulerProblem::make_grid_and_dofs() + { + if (testcase == 0) + { + Point lower_left; + for (unsigned int d = 1; d < dim; ++d) + lower_left[d] = -5; + Point upper_right; + upper_right[0] = 10; + for (unsigned int d = 1; d < dim; ++d) + upper_right[d] = 5; + + std::vector 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>(0)); + } + else if (testcase == 1) + { + GridGenerator::channel_with_cylinder(triangulation, 0.03, 1, 0, true); + euler_operator.set_inflow_boundary( + 0, std::make_shared>(0)); + euler_operator.set_subsonic_outflow_boundary( + 1, std::make_shared>(0)); + euler_operator.set_wall_boundary(2); + euler_operator.set_wall_boundary(3); + if (dim == 3) + euler_operator.set_body_force( + std::make_shared>( + std::vector({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 + void EulerProblem::output_results(const unsigned int result_number) + { + Tensor<1, 3> errors = + euler_operator.compute_errors(ExactSolution(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 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 reference; + if (testcase == 0 && dim == 2) + { + reference.reinit(solution); + euler_operator.project(ExactSolution(time), reference); + reference.sadd(-1., 1, solution); + std::vector 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 + 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 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::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 + void EulerProblem::run() + { + { + const unsigned int n_vect_number = VectorizedArray::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::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 rk_register_1; + LinearAlgebra::distributed::Vector rk_register_2; + rk_register_1.reinit(solution); + rk_register_2.reinit(solution); + + euler_operator.project(ExactSolution(time), solution); + + typename Triangulation::active_cell_iterator cell = triangulation + .begin_active(), + endc = + triangulation.end(); + double min_vertex_distance = std::numeric_limits::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(time / output_tick) != + static_cast((time - time_step) / output_tick) || + time >= FINAL_TIME - 1e-12) + output_results( + static_cast(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 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; } -- 2.39.5