From: Matthias Maier Date: Thu, 13 Feb 2020 19:38:44 +0000 (-0600) Subject: address @bangerth's review comments X-Git-Tag: v9.2.0-rc1~443^2~16 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6d552f9d5bc7ea412868a16f6faf49fb2a8efc08;p=dealii.git address @bangerth's review comments --- diff --git a/examples/step-69/doc/intro.dox b/examples/step-69/doc/intro.dox index ba3a649081..2cbc8fb574 100644 --- a/examples/step-69/doc/intro.dox +++ b/examples/step-69/doc/intro.dox @@ -27,17 +27,16 @@ time integration, see @cite GuermondEtAl2018. This tutorial presents a first-order scheme for solving compressible Euler's equations that is based on three ingredients: a -collocation-type discretization of Euler's equations in context of -finite elements; a graph-viscosity stabilization based on a +collocation-type discretization of Euler's equations in the context +of finite elements; a graph-viscosity stabilization based on a guaranteed upper bound of the local wave speed; and explicit time-stepping. As such, the ideas and techniques presented in this tutorial step are drastically different from those used in Step-33, which focuses on the use of automatic differentiation. From a programming perspective this tutorial will focus on a number of techniques found in large-scale -computations: hybrid thread-MPI parallelization; efficient -local numbering of degrees of freedom; concurrent post-processing and -write-out of results using worker threads; as well as checkpointing and -restart. +computations: hybrid thread-MPI parallelization; efficient local numbering +of degrees of freedom; concurrent post-processing and write-out of results +using worker threads; as well as checkpointing and restart. It should be noted that first-order schemes in the context of hyperbolic conservation laws require prohibitively many degrees of freedom to resolve @@ -54,11 +53,9 @@ the second-order scheme @cite GuermondEtAl2018. The compressible Euler's equations of gas dynamics are written in conservative form as follows: - @f{align} \mathbf{u}_t + \text{div} \, \mathbb{f}(\mathbf{u}) = \boldsymbol{0} , @f} - where $\mathbf{u}(\textbf{x},t):\mathbb{R}^{d} \times \mathbb{R} \rightarrow \mathbb{R}^{d+2}$, and $\mathbb{f}(\mathbf{u}):\mathbb{R}^{d+2} \rightarrow \mathbb{R}^{(d+2) \times d}$, and $d \geq 1$ is the space @@ -69,7 +66,6 @@ $\textbf{u} = [\rho, \textbf{m},E]^{\top}$: where $\rho \in \mathbb{R}^+$ denotes the density, $\textbf{m} \in \mathbb{R}^d$ is the momentum, and $E \in \mathbb{R}^+$ is the total energy of the system. The flux of the system $\mathbb{f}(\mathbf{u})$ is defined as - @f{align*} \mathbb{f}(\textbf{u}) = @@ -79,7 +75,6 @@ $\mathbb{f}(\mathbf{u})$ is defined as \tfrac{\textbf{m}^\top}{\rho} (E + p) \end{bmatrix}, @f} - where $\mathbb{I} \in \mathbb{R}^{d \times d}$ is the identity matrix and $\otimes$ denotes the tensor product. Here, we have introduced the pressure $p$ that, in general, is defined by an closed-form equation of state. @@ -90,7 +85,6 @@ p = p(\textbf{u}) := (\gamma -1) \Big(E - \tfrac{|\textbf{m}|_{\ell^2}^2}{2\,\rho} \Big), @f} - where the factor $\gamma \in (1,5/3]$ denotes the ratio of specific heats, and $|\cdot|_{\ell^2}$ denotes the Euclidian norm. @@ -115,19 +109,16 @@ understanding of hyperbolic conservation laws was to assume that the solution is formally defined as $\mathbf{u} := \lim_{\epsilon \rightarrow 0^+} \mathbf{u}^{\epsilon}$ where $\mathbf{u}^{\epsilon}$ is the solution of the parabolic regularization - @f{align} \mathbf{u}_t^{\epsilon} + \text{div} \, \mathbb{f}(\mathbf{u}^{\epsilon}) = {\epsilon} \Delta \mathbf{u}^{\epsilon}. @f} - Such solutions, which are understood as the solution recovered in the zero-viscosity limit, are often refered to as viscosity solutions. Global existence and uniqueness of such solutions is a widely open issue. However, we know at least that if such viscosity solutions exists they have to satisfy the constraint $\textbf{u}(\mathbf{x},t) \in \mathcal{B}$ for all $\mathbf{x} \in \Omega$ and $t \geq 0$ where - @f{align} \mathcal{B} = \big\{ \textbf{u} = [\rho, \textbf{m},E]^{\top} \in \mathbb{R}^{d+2} \, \big | @@ -139,9 +130,7 @@ all $\mathbf{x} \in \Omega$ and $t \geq 0$ where s(\mathbf{u}) \geq \min_{x \in \Omega} s(\mathbf{u}_0(\mathbf{x})) \big\}. @f} - Here, $s(\mathbf{u})$ denotes the specific entropy - @f{align} s(\mathbf{u}) = \ln \Big(\frac{p(\mathbf{u})}{\rho^{\gamma}}\Big). @f} @@ -167,14 +156,12 @@ $\mathbf{u}(\mathbf{x},t)$ remains in $\mathcal{B}$. Following Step-9, Step-12, and Step-33, at this point it might look tempting to base a discretization of Euler's equations on a (semi-discrete) variational formulation: - @f{align*} (\partial_t\mathbf{u}_{h},\textbf{v}_h)_{L^2(\Omega)} - ( \mathbb{f}(\mathbf{u}_{h}) ,\text{grad} \, \textbf{v}_{h})_{L^2(\Omega)} + s_h(\mathbf{u}_{h},\textbf{v}_h)_{L^2(\Omega)} = \boldsymbol{0} \quad\forall \textbf{v}_h \in \mathbb{V}_h. @f} - Here, $\mathbb{V}_h$ is an appropriate finite element space, and $s_h(\cdot,\cdot)_{L^2(\Omega)}$ is some linear stabilization method (possibly complemented with some ad-hoc shock-capturing technique, see for @@ -184,11 +171,9 @@ are based on such a (semi-discrete) variational approach. Fundamentally, from an analysis perspective, variational discretizations are conceived to provide some notion of global (integral) stabiliy, meaning an estimate of the form - @f{align*} |\!|\!| \mathbf{u}_{h}(t) |\!|\!| \leq |\!|\!| \mathbf{u}_{h}(0) |\!|\!| @f} - holds true, where $|\!|\!| \cdot |\!|\!| $ could represent the $L^2(\Omega)$-norm or, more generally, some discrete (possibly mesh dependent) energy-norm. Variational discretizations of hyperbolic @@ -213,12 +198,10 @@ have an edge in many regards. In this tutorial step we therefore depart from variational schemes. We will present a completely algebraic formulation (with the flavor of a collocation-type scheme) that preserves constraints pointwise, i.e., - @f{align*} \textbf{u}_h(\mathbf{x}_i,t) \in \mathcal{B} \;\text{at every node}\;\mathbf{x}_i\;\text{of the mesh}. @f} - Contrary to finite difference/volume schemes, the scheme implemented in this step maximizes the use of finite element software infrastructure, works in any mesh, in any space dimension, and is theoretically guaranteed @@ -250,14 +233,12 @@ where $\mathbf{x}_i \in \mathbb{R}^d$. Then each index $i \in a scalar-valued shape function $\phi_i$. With this notation at hand we can define the scheme as: - @f{align*} m_i \frac{\mathbf{U}_i^{n+1} - \mathbf{U}_i^{n}}{\tau} + \sum_{j \in \mathcal{I}(i)} \mathbb{f}(\mathbf{U}_j^{n})\cdot \mathbf{c}_{ij} - d_{ij} \mathbf{U}_j^{n} = \boldsymbol{0} \, , @f} - -Where +where - $m_i := \int_{\Omega} \phi_i \, \mathrm{d}\mathbf{x}$ - $\tau$ is the time step size - $\mathbf{c}_{ij} := \int_{\Omega} \nabla\phi_j\phi_i \, @@ -295,7 +276,6 @@ Consider the following pseudo-code, illustrating a possible straight forward strategy for computing the solution $\textbf{U}^{n+1}$ at a new time $t_{n+1} = t_n + \tau_n$ given a known state $\textbf{U}^{n}$ at time $t_n$: - @f{align*} &\textbf{for } i \in \mathcal{V} \\ &\ \ \ \ \{\mathbf{c}_{ij}\}_{j \in \mathcal{I}(i)} \leftarrow @@ -325,11 +305,12 @@ either collect (from) or write (into) global matrices and vectors. dimensions, first-order polynomial space $\mathbb{Q}^1$, and that $\mathbf{x}_i$ is an interior node (i.e. $\mathbf{x}_i$ is not on the boundary of the domain ) then: $\{\textbf{U}_j^n\}_{j \in \mathcal{I}(i)}$ should contain -nine state-vectors (i.e. all the states in the patch/macro element associated to -the shape function $\phi_i$). This is one of the major differences with the -usual cell-based loop where the gather functionality (encoded in -FEValuesBase.get_function_values() in the case of deal.ii) only -collects values for the local cell (just a subset of the patch). +nine state vector elements (i.e. all the states in the patch/macro element +associated to the shape function $\phi_i$). This is one of the major +differences with the usual cell-based loop where the gather functionality +(encoded in FEValuesBase.get_function_values() in the case +of deal.ii) only collects values for the local cell (just a subset of the +patch). The actual implementation will deviate from above code in one key aspect: the time-step size $\tau$ has to be chosen subject to a CFL condition diff --git a/examples/step-69/doc/kind b/examples/step-69/doc/kind index 86a44aa1ef..e62f4e7222 100644 --- a/examples/step-69/doc/kind +++ b/examples/step-69/doc/kind @@ -1 +1 @@ -time dependent +fluids diff --git a/examples/step-69/step-69.cc b/examples/step-69/step-69.cc index daf47dfb17..3f860ae092 100644 --- a/examples/step-69/step-69.cc +++ b/examples/step-69/step-69.cc @@ -909,15 +909,15 @@ namespace Step69 // indices (i,l) of the entry in order to retrieve its // actual value. We should expect gather to be slightly // more expensive than gather_get_entry. The use of - // gather will be limited to the task of computing the - // algebraic viscosity $d_{ij}$ in the particular case that when + // gather will be limited to the task of computing the + // algebraic viscosity $d_{ij}$ in the particular case that when // both $i$ and $j$ lie at the boundary. // - // @note The reader should be aware that accessing an arbitrary - // (i,l) entry of a matrix (say for instance Trilinos or PETSc - // matrices) is in general unacceptably expensive. Here is where we might - // want to keep an eye on complexity: we want this operation to have - // constant complexity, which is the case of the current implementation + // @note The reader should be aware that accessing an arbitrary + // (i,l) entry of a matrix (say for instance Trilinos or PETSc + // matrices) is in general unacceptably expensive. Here is where we might + // want to keep an eye on complexity: we want this operation to have + // constant complexity, which is the case of the current implementation // using deal.ii matrices. template @@ -932,7 +932,7 @@ namespace Step69 // gather (second interface): this second function // signature having two input arguments will be used to gather the - // state at a node i and return it as a + // state at a node i and return it as a // Tensor<1,problem_dimension> for our convenience. template @@ -946,7 +946,7 @@ namespace Step69 } // scatter: this function has three input arguments, the - // first one is meant to be a "global object" (say a locally owned or + // first one is meant to be a "global object" (say a locally owned or // locally relevant vector), the second argument which could be a // Tensor<1,problem_dimension>, and the last argument // which represents a index of the global object. This function will be @@ -1002,9 +1002,9 @@ namespace Step69 // here $T$ denotes elements, // $\text{supp}(\phi_i)$ the support of the shape function $\phi_i$, // $F$ are faces of the element $T$, and $\mathbf{x}_{q,F}$ - // are quadrature points on such face. Note that this formula for - // $\widehat{\boldsymbol{\nu}}_i$ is nothing else than some form of - // weighted averaging. Other more sophisticated definitions for $\nu_i$ + // are quadrature points on such face. Note that this formula for + // $\widehat{\boldsymbol{\nu}}_i$ is nothing else than some form of + // weighted averaging. Other more sophisticated definitions for $\nu_i$ // are possible but none of them have much influence in theory or practice. template @@ -1434,7 +1434,7 @@ namespace Step69 // In this section we describe the implementation of the class members of // the ProblemDescription class. Most of the code here is - // specific for compressible Euler's equations with an ideal gas law. + // specific for compressible Euler's equations with an ideal gas law. // If we wanted to re-purpose Step-69 for a different conservation law // (say for: instance the shallow water equation) most of the // implementation of this class would have to change. But most of the other @@ -1529,14 +1529,14 @@ namespace Step69 // approximation for the intermediate pressure $p^*$, see for instance // Equation (4.46), page 128 in @cite Toro2009. // - // The estimate returned by lambda_max_two_rarefaction - // is guaranteed to be an upper bound, it is in general quite sharp, and - // overall sufficient for our purposes. However, for some specific situations - // (in particular when one of states is close to vacuum conditions) such - // an estimate will be overly pessimistic. That's why we used a second - // estimate to avoid this degeneracy that will be invoked by a call to the - // function lambda_max_expansion. The most important function - // here is compute_lambda_max which takes the minimum between + // The estimate returned by lambda_max_two_rarefaction + // is guaranteed to be an upper bound, it is in general quite sharp, and + // overall sufficient for our purposes. However, for some specific situations + // (in particular when one of states is close to vacuum conditions) such + // an estimate will be overly pessimistic. That's why we used a second + // estimate to avoid this degeneracy that will be invoked by a call to the + // function lambda_max_expansion. The most important function + // here is compute_lambda_max which takes the minimum between // the estimates returned by lambda_max_two_rarefaction and // lambda_max_expansion. // @@ -1596,10 +1596,10 @@ namespace Step69 // primitive state $[\rho, u, p, a]$ and a given pressure $p^\ast$ // @cite GuermondPopov2016 Eqn. (3.7): // @f{align*} - // \lambda^- = u - a\,\sqrt{1 + \frac{\gamma+1}{2\gamma} + // \lambda^- = u - a\,\sqrt{1 + \frac{\gamma+1}{2\gamma} // \left(\frac{p^\ast-p}{p}\right)_+} // @f} - // Here, the $(\cdot)_{+}$ denotes the positive part of the given + // Here, the $(\cdot)_{+}$ denotes the positive part of the given // argument. DEAL_II_ALWAYS_INLINE inline double @@ -1620,7 +1620,7 @@ namespace Step69 // Analougously @cite GuermondPopov2016 Eqn. (3.8): // @f{align*} - // \lambda^+ = u + a\,\sqrt{1 + \frac{\gamma+1}{2\gamma} + // \lambda^+ = u + a\,\sqrt{1 + \frac{\gamma+1}{2\gamma} // \left(\frac{p^\ast-p}{p}\right)_+} // @f} @@ -1699,7 +1699,7 @@ namespace Step69 } } // namespace - // The following is the main function that we are going to call in order to + // The following is the main function that we are going to call in order to // compute $\lambda_{\text{max}} (\mathbf{U}_i^{n},\mathbf{U}_j^{n}, // \textbf{n}_{ij})$. We simply compute both maximal wavespeed estimates // and return the minimum. @@ -1811,8 +1811,8 @@ namespace Step69 static constexpr auto gamma = ProblemDescription::gamma; // The following lambda function translates a given primitive 1d state - // (density $\rho$, velocity $u$, and pressure $p$) into a - // conserved n-dimensional state (density $\rho$, momentum + // (density $\rho$, velocity $u$, and pressure $p$) into a + // conserved n-dimensional state (density $\rho$, momentum // $\mathbf{m}$, and total energy $E$). Note that we // capture // the this pointer and thus access to @@ -1935,12 +1935,12 @@ namespace Step69 // \mathbf{U}_j^{n}, \textbf{n}_{ij}) = \lambda_{\text{max}} // (\mathbf{U}_j^{n}, \mathbf{U}_i^{n}, \textbf{n}_{ji})$ do not // necessarily hold true. The only mathematically safe solution for this - // dilemma is to compute both of them $d_{ij}$ and $d_{ji}$ and + // dilemma is to compute both of them $d_{ij}$ and $d_{ji}$ and // take the maximum. // // Overall, the computation of $d_{ij}$ is quite expensive. In // order to save some computing time we exploit the fact that the viscosity - // matrix has to be symmetric (as mentioned above): we only compute + // matrix has to be symmetric (as mentioned above): we only compute // the upper-triangular entries of $d_{ij}$ and copy the // corresponding entries to the lower-triangular counterpart. // @@ -1968,7 +1968,7 @@ namespace Step69 { const auto j = jt->column(); - // We only compute $d_{ij}$ if $j < i$ (upper triangular + // We only compute $d_{ij}$ if $j < i$ (upper triangular // entries) and later copy the values over to $d_{ji}$. if (j >= i) continue; @@ -1984,7 +1984,7 @@ namespace Step69 double d = norm * lambda_max; // If both support points happen to be at the boundary we - // have to compute $d_{ji}$ as well and then take + // have to compute $d_{ji}$ as well and then take // $max(d_{ij},d_{ji})$: if (boundary_normal_map.count(i) != 0 && boundary_normal_map.count(j) != 0)