]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make text presentation more uniform
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 24 Mar 2020 15:23:02 +0000 (16:23 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 24 Mar 2020 15:23:02 +0000 (16:23 +0100)
examples/step-67/doc/intro.dox
examples/step-67/doc/results.dox
examples/step-67/step-67.cc

index 4470ed139a578fab8fd8ef66ef0c942027414bce..02a599009e35d17618196b64392cbf9a54b8631d 100644 (file)
@@ -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 <i>effective</i> 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 <i>barely</i> 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
 <a href="#Results">results section</a>.) In fact, it
 is so cheap that it is limited by the bandwidth of reading the source vector,
index 2c31e0e9dd820b2ae89a945f90f8dd5576759a31..8e1aca2b07370257cba986674e012fc8f0754261 100644 (file)
@@ -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
index 4f413b2c871553fd501d13ee3d97d145b1d29633..d680b405295da38b50a941221da1c24cd3acb4e8 100644 (file)
@@ -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 <typename VectorType, typename Operator>
     void perform_time_step(const Operator &pde_operator,
                            const double    current_time,
                            const double    time_step,
                            VectorType &    solution,
-                           VectorType &    vec_Ti,
-                           VectorType &    vec_Ki)
+                           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<Number> &current_Ti,
-                  LinearAlgebra::distributed::Vector<Number> &      vec_Ki,
+                  const LinearAlgebra::distributed::Vector<Number> &current_ri,
+                  LinearAlgebra::distributed::Vector<Number> &      vec_ki,
                   LinearAlgebra::distributed::Vector<Number> &      solution,
-                  LinearAlgebra::distributed::Vector<Number> &next_Ti) const;
+                  LinearAlgebra::distributed::Vector<Number> &next_ri) const;
 
     void project(const Function<dim> &                       function,
                  LinearAlgebra::distributed::Vector<Number> &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<Number> &current_Ti,
-    LinearAlgebra::distributed::Vector<Number> &      vec_Ki,
+    const LinearAlgebra::distributed::Vector<Number> &current_ri,
+    LinearAlgebra::distributed::Vector<Number> &      vec_ki,
     LinearAlgebra::distributed::Vector<Number> &      solution,
-    LinearAlgebra::distributed::Vector<Number> &      next_Ti) const
+    LinearAlgebra::distributed::Vector<Number> &      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<dim, Number>::DataAccessOnFaces::values,
                 MatrixFree<dim, Number>::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<void(const unsigned int, const unsigned int)>(),
         [&](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 <int dim>
   void EulerProblem<dim>::run()
   {

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.