From 78c2ce51e16d86b9782f190b76a253d61fe76777 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 24 Mar 2020 16:23:02 +0100 Subject: [PATCH] Make text presentation more uniform --- examples/step-67/doc/intro.dox | 57 ++++++++++++---------- examples/step-67/doc/results.dox | 8 ++- examples/step-67/step-67.cc | 84 ++++++++++++++++---------------- 3 files changed, 76 insertions(+), 73 deletions(-) diff --git a/examples/step-67/doc/intro.dox b/examples/step-67/doc/intro.dox index 4470ed139a..02a599009e 100644 --- a/examples/step-67/doc/intro.dox +++ b/examples/step-67/doc/intro.dox @@ -79,7 +79,7 @@ For spatial discretization, we use a high-order discontinuous Galerkin (DG) discretization, using a solution expansion of the form @f[ \mathbf{w}_h(\mathbf{x}, t) = -\sum_{j=0}^{n_\mathbf{dofs}} \boldsymbol{\varphi}_j(\mathbf{x}) {w}_j(t). +\sum_{j=1}^{n_\mathbf{dofs}} \boldsymbol{\varphi}_j(\mathbf{x}) {w}_j(t). @f] Here, $\boldsymbol{\varphi}_j$ denotes the $j$th basis function, written in vector form with separate shape functions for the different components and @@ -89,7 +89,7 @@ functions and the time dependence in the unknown coefficients $w_j$. As opposed to the continuous finite element method where some shape functions span across element boundaries, the shape functions are local to a single element in DG methods, with a discontinuity from one element to the next. The -connectio of the solution from one cell to its neighbors is instead +connection of the solution from one cell to its neighbors is instead imposed by the numerical fluxes specified below. This allows for some additional flexibility, for example to introduce directionality in the numerical method by, e.g., upwinding. @@ -252,7 +252,7 @@ simply use arithmetic means, $\bar{\mathbf{u}} = \frac{\mathbf{u}^- + \sqrt{\gamma p^{\pm} / \rho^{\pm}}$, in this tutorial program, and leave other variants to a possible extension. We also note that the HLL flux has been extended in the literature to the so-called HLLC flux, where C stands for the -ability to represent so-called contact discontinuities. +ability to represent contact discontinuities. At the boundaries with no neighboring state $\mathbf{w}^+$ available, it is common practice to deduce suitable exterior values from the boundary @@ -379,15 +379,18 @@ polynomial degree $p$ and the physics of the Euler operator, i.e., the eigenvalues of the linearization of $\mathbf F(\mathbf w)$ with respect to $\mathbf{w}$. In this program, we set the time step as follows: @f[ -\Delta t = \frac{Cr}{p^{1.5}}\left(\frac{1}{\max\left[\frac{\|\mathbf{u}\|}{h_u} + - \frac{c}{h_c}\right]}\right), +\Delta t = \frac{\mathrm{Cr}}{p^{1.5}}\left(\frac{1} + {\max\left[\frac{\|\mathbf{u}\|}{h_u} + \frac{c}{h_c}\right]}\right), @f] -with the maximum taken over all quadrature points and all cells. The power -$p^{1.5}$ used for the polynomial scaling is heuristic and represents the -closest fit for polynomial degrees between 1 and 8, see e.g. -@cite SchoederKormann2018. In the limit of higher degrees, $p>10$, a scaling -of $p^2$ is more accurate, related to the inverse estimates typically used for +with the maximum taken over all quadrature points and all cells. The +dimensionless number $\mathrm{Cr}$ denotes the Courant number and can be +chosen up to a maximally stable number $\mathrm{Cr}_\text{max}$, whose value +depends on the selected time stepping method and its stability properties. The +power $p^{1.5}$ used for the polynomial scaling is heuristic and represents +the closest fit for polynomial degrees between 1 and 8, see e.g. @cite +SchoederKormann2018. In the limit of higher degrees, $p>10$, a scaling of +$p^2$ is more accurate, related to the inverse estimates typically used for interior penalty methods. Regarding the effective mesh sizes $h_u$ and $h_c$ used in the formula, we note that the convective transport is directional. Thus an appropriate scaling is to use the element length in the @@ -460,21 +463,21 @@ c_5 & b_1 & b_2 & b_3 & a_4 \\ @f] With such a definition, the update to $\mathbf{w}_h^n$ shares the storage with the information for the intermediate values $\mathbf{k}_i$. Starting with -$\mathbf{w}^{n+1}=\mathbf{w}^n$ and $\mathbf{t}_1 = \mathbf{w}^n$, the update +$\mathbf{w}^{n+1}=\mathbf{w}^n$ and $\mathbf{r}_1 = \mathbf{w}^n$, the update in each of the $s$ stages simplifies to @f[ \begin{aligned} \mathbf{k}_i &= -\mathcal M^{-1} \mathcal L_h\left(t^n+c_i\Delta t, \mathbf{t}_{i} \right),\\ -\mathbf{t}_{i+1} &= \mathbf{w}_h^{n+1} + \Delta t \, a_i \mathbf{k}_i,\\ +\mathcal M^{-1} \mathcal L_h\left(t^n+c_i\Delta t, \mathbf{r}_{i} \right),\\ +\mathbf{r}_{i+1} &= \mathbf{w}_h^{n+1} + \Delta t \, a_i \mathbf{k}_i,\\ \mathbf{w}_h^{n+1} &= \mathbf{w}_h^{n+1} + \Delta t \, b_i \mathbf{k}_i. \end{aligned} @f] -Besides the vector $\mathbf w_h$ that is successively updated, this scheme +Besides the vector $\mathbf w_h^{n+1}$ that is successively updated, this scheme only needs two auxiliary vectors, namely the vector $\mathbf{k}_i$ to hold the -evaluation of the differential operator, and the vector $\mathbf{t}_i$ that +evaluation of the differential operator, and the vector $\mathbf{r}_i$ that holds the right-hand side for the differential operator application. In -subsequent stages $i$, the values $\mathbf{k}_i$ and $\mathbf{t}_i$ can use +subsequent stages $i$, the values $\mathbf{k}_i$ and $\mathbf{r}_i$ can use the same storage. The main advantages of low-storage variants are the reduced memory consumption @@ -587,12 +590,12 @@ variable-coefficient integrand. It has been shown in the DG literature that aliasing errors can introduce unphysical oscillations in the numerical solution for barely resolved simulations. The fact that aliasing mostly affects coarse resolutions -- whereas finer meshes with the same scheme -otherwise work fine -- is not surprising because well-resolved simulations -have tend to be smooth on length-scales of a cell (i.e., they have +work fine -- is not surprising because well-resolved simulations +tend to be smooth on length-scales of a cell (i.e., they have small coefficients in the higher polynomial degrees that are missed by too few quadrature points, whereas the main solution contribution in the lower polynomial degrees is still well-captured -- this is simply a consequence of Taylor's -theorem). To address this topic, various approaches have been proposed in the +theorem. To address this topic, various approaches have been proposed in the DG literature. One technique is filtering which damps the solution components pertaining to higher polynomial degrees. As the chosen nodal basis is not hierarchical, this would mean to transform from the nodal basis into a @@ -608,11 +611,11 @@ over-integration or consistent integration. The latter name is most common in the context of the incompressible Navier-Stokes equations, where the $\mathbf{u}\otimes \mathbf{u}$ nonlinearity results in polynomial integrands of degree $3p$ (when also considering the test function), which can be -integrated exactly with $\textbf{floor}\left(\frac{3p}{2}\right)+1$ quadrature +integrated exactly with $\textrm{floor}\left(\frac{3p}{2}\right)+1$ quadrature points per direction as long as the element geometry is affine. In the context of the Euler equations with non-polynomial integrands, the choice is less clear. Depending on the variation in the various variables both -$\textbf{floor}\left(\frac{3p}{2}\right)+1$ or $2p+1$ points (integrating +$\textrm{floor}\left(\frac{3p}{2}\right)+1$ or $2p+1$ points (integrating exactly polynomials of degree $3p$ or $4p$, respectively) are common. To reflect this variability in the choice of quadrature in the program, we @@ -646,7 +649,7 @@ the error made by taking a diagonal mass matrix and ignoring the rest (a variant of mass lumping, though not the one with an additional integration error as utilized in step-48) has been shown to not alter discretization accuracy. The Lagrange basis in the points of Gaussian quadrature is sometimes -also referred to as a collocation setup , as the nodal points of the +also referred to as a collocation setup, as the nodal points of the polynomials coincide (= are "co-located") with the points of quadrature, obviating some interpolation operations. Given the fact that we want to use more quadrature points for nonlinear terms in $\mathcal L_h$, however, the collocation @@ -702,15 +705,15 @@ $S_{\mathrm{1D}}^{-\mathrm{T}}$ rather than $S_{\mathrm{1D}}$. The class MatrixFreeOperators::CellwiseInverseMassMatrix implements this operation: It changes from the basis contained in the finite element (in this case, FE_DGQ) to the Lagrange basis in Gaussian quadrature points. Here, the -diagonal mass matrix can be evaluate, which is simply the inverse of -the `JxW` factors (i.e., the quadrature weight times the determinant of the +inverse of a diagonal mass matrix can be evaluated, which is simply the inverse +of the `JxW` factors (i.e., the quadrature weight times the determinant of the Jacobian from reference to real coordinates). Once this is done, we can change back to the standard nodal Gauss-Lobatto basis. The advantage of this particular way of applying the inverse mass matrix is -the fact that it is of similar cost as the forward application of a mass -matrix, and cheaper than the evaluation of the spatial operator $\mathcal L_h$ -which is more costly due to over-integration and face integrals. (We +a cost similar to the forward application of a mass matrix, which is cheaper +than the evaluation of the spatial operator $\mathcal L_h$ +with over-integration and face integrals. (We will demonstrate this with detailed timing information in the results section.) In fact, it is so cheap that it is limited by the bandwidth of reading the source vector, diff --git a/examples/step-67/doc/results.dox b/examples/step-67/doc/results.dox index 2c31e0e9dd..8e1aca2b07 100644 --- a/examples/step-67/doc/results.dox +++ b/examples/step-67/doc/results.dox @@ -45,9 +45,7 @@ speed of sound (and thus faster propagation) leaves the computational domain between times 5 and 6.5. After that point, the flow is simply uniform in the same direction, and the maximum velocity of the gas is reduced compared to the previous state where the uniform velocity was overlaid -by the vortex. Our time step formula recognizes this and only -uses the acoustic limit in the last part of the simulation when -determining the time step size. +by the vortex. Our time step formula recognizes this effect. The final block of output shows detailed information about the timing of individual parts of the programs; it breaks this down by showing @@ -124,8 +122,8 @@ vector updates (797 vs 77.5 seconds). The reason for this non-optimal increase in run time can be traced back to cache effects on the given hardware (with 40 MB of L2 cache and 55 MB of L3 cache): While not all of the relevant data fits into caches for 9.4 million DoFs (one vector takes 75 MB and we have three -vectors plus some additional data in MatrixFree), there is capacity for almost -half of one vector nonetheless. Given that modern caches are more sophisticated than +vectors plus some additional data in MatrixFree), there is capacity for one and +a half vector nonetheless. Given that modern caches are more sophisticated than the naive least-recently-used strategy (where we would have little re-use as the data is used in a streaming-like fashion), we can assume that a sizeable fraction of data can indeed be delivered from caches for the 9.4 million DoFs diff --git a/examples/step-67/step-67.cc b/examples/step-67/step-67.cc index 4f413b2c87..d680b40529 100644 --- a/examples/step-67/step-67.cc +++ b/examples/step-67/step-67.cc @@ -211,8 +211,8 @@ namespace Euler_DG // 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 + // update vector $\mathbf{r}_i$ that gets evaluated during the stages, plus + // one vector $\mathbf{k}_i$ to hold the evaluation of the operator. Such a // Runge--Kutta setup reduces the memory storage and memory access. As the // memory bandwidth is often the performance-limiting factor on modern // hardware when the evaluation of the differential operator is @@ -322,7 +322,7 @@ namespace Euler_DG } // 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 + // evaluate the operator, prepare the $\mathbf{r}_i$ vector for the next // evaluation, and update the solution vector $\mathbf{w}$. We hand off // the work to the `pde_operator` involved in order to be able to merge // the vector operations of the Runge--Kutta setup with the evaluation of @@ -331,30 +331,30 @@ namespace Euler_DG // // 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 old solution $\mathbf{w}^n$ rather than a $\mathbf r_i$ vector, so // the first argument is `solution`. We here let the stage vector - // $\mathbf{T}_i$ also hold the temporary result of the evaluation, as it + // $\mathbf{r}_i$ also hold the temporary result of the evaluation, as it // is not used otherwise. For all subsequent stages, we use the vector - // `vec_Ki` as the second vector argument to store the result of the + // `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 + // skip the computation of the vector $\mathbf{r}_{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) + VectorType & vec_ri, + 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, + vec_ri, solution, - vec_Ti); + vec_ri); double sum_previous_bi = 0; for (unsigned int stage = 1; stage < bi.size(); ++stage) { @@ -364,10 +364,10 @@ namespace Euler_DG (stage == bi.size() - 1 ? 0 : ai[stage] * time_step), - vec_Ti, - vec_Ki, + vec_ri, + vec_ki, solution, - vec_Ti); + vec_ri); sum_previous_bi += bi[stage - 1]; } } @@ -700,10 +700,10 @@ namespace Euler_DG 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, + const LinearAlgebra::distributed::Vector ¤t_ri, + LinearAlgebra::distributed::Vector & vec_ki, LinearAlgebra::distributed::Vector & solution, - LinearAlgebra::distributed::Vector &next_Ti) const; + LinearAlgebra::distributed::Vector &next_ri) const; void project(const Function & function, LinearAlgebra::distributed::Vector &solution) const; @@ -1335,15 +1335,15 @@ namespace Euler_DG // 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 + // to the vectors, namely `next_ri = solution + factor_ai * k_i` and + // `solution += factor_solution * k_i`. Rather than performing these // steps through the vector interfaces, we here present an alternative // strategy that is faster on cache-based architectures. As the memory // consumed by the vectors is often much larger than what fits into caches, // the data has to effectively come from the slow RAM memory. The situation // can be improved by loop fusion, i.e., performing both the updates to - // `next_Ki` and `solution` within a single sweep. In that case, we would - // read the two vectors `rhs` and `solution` and write into `next_Ki` and + // `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 @@ -1360,7 +1360,7 @@ namespace Euler_DG // 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` + // Runge--Kutta stage when we do not need to update the `next_ri` // vector. This strategy gives a considerable speedup. Whereas the inverse // mass matrix and vector updates take more than 60% of the computational // time with default vector updates on a 40-core machine, the percentage is @@ -1371,10 +1371,10 @@ namespace Euler_DG const Number current_time, const Number factor_solution, const Number factor_ai, - const LinearAlgebra::distributed::Vector ¤t_Ti, - LinearAlgebra::distributed::Vector & vec_Ki, + const LinearAlgebra::distributed::Vector ¤t_ri, + LinearAlgebra::distributed::Vector & vec_ki, LinearAlgebra::distributed::Vector & solution, - LinearAlgebra::distributed::Vector & next_Ti) const + LinearAlgebra::distributed::Vector & next_ri) const { { TimerOutput::Scope t(timer, "rk_stage - integrals L_h"); @@ -1388,8 +1388,8 @@ namespace Euler_DG &EulerOperator::local_apply_face, &EulerOperator::local_apply_boundary_face, this, - vec_Ki, - current_Ti, + vec_ki, + current_ri, true, MatrixFree::DataAccessOnFaces::values, MatrixFree::DataAccessOnFaces::values); @@ -1401,8 +1401,8 @@ namespace Euler_DG data.cell_loop( &EulerOperator::local_apply_inverse_mass_matrix, this, - next_Ti, - vec_Ki, + next_ri, + vec_ki, std::function(), [&](const unsigned int start_range, const unsigned int end_range) { const Number ai = factor_ai; @@ -1412,9 +1412,9 @@ namespace Euler_DG 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 k_i = next_ri.local_element(i); const Number sol_i = solution.local_element(i); - solution.local_element(i) = sol_i + bi * K_i; + solution.local_element(i) = sol_i + bi * k_i; } } else @@ -1422,10 +1422,10 @@ namespace Euler_DG 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 k_i = next_ri.local_element(i); const Number sol_i = solution.local_element(i); - solution.local_element(i) = sol_i + bi * K_i; - next_Ti.local_element(i) = sol_i + ai * K_i; + solution.local_element(i) = sol_i + bi * k_i; + next_ri.local_element(i) = sol_i + ai * k_i; } } }); @@ -2044,13 +2044,15 @@ namespace Euler_DG // 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. + // low-storage integrator. We call these vectors `rk_register_1` and + // `rk_register_2`, and use the first vector to represent the quantity + // $\mathbf{r}_i$ and the second one for $\mathbf{k}_i$. 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() { -- 2.39.5