From: Matthias Maier Date: Thu, 26 Dec 2019 14:41:19 +0000 (-0600) Subject: reindent, layout changes, and address some review comments X-Git-Tag: v9.2.0-rc1~443^2~35 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=470d470ebbb4b846255c4a13f14f0efab5566e25;p=dealii.git reindent, layout changes, and address some review comments --- diff --git a/examples/step-69/doc/builds-on b/examples/step-69/doc/builds-on index e69de29bb2..48f6065a7c 100644 --- a/examples/step-69/doc/builds-on +++ b/examples/step-69/doc/builds-on @@ -0,0 +1 @@ +step-33 step-40 diff --git a/examples/step-69/doc/intro.dox b/examples/step-69/doc/intro.dox index 4ec7cc963d..ec8b1890b6 100644 --- a/examples/step-69/doc/intro.dox +++ b/examples/step-69/doc/intro.dox @@ -12,13 +12,7 @@ such it is presented primarily for educational purposes. For actual research computations you might want to consider exploring a corresponding high-performance implementation of a second-order accurate scheme that uses convex limiting techniques, and strong stability-preserving (SSP) -time integration @cite GuermondEtAl2018. The repository can be found at -TODO. - - - - - +time integration, see @cite GuermondEtAl2018.

Introduction

@@ -47,11 +41,6 @@ the programming techniques) before jumping into full research codes such as the second-order scheme @cite GuermondEtAl2018 maintained TODO. - - - - -

Euler's equations of gas dynamics

@@ -98,10 +87,6 @@ where the factor $\gamma \in (1,5/3]$ denotes the specific heats, and $\|\,.\|$ denotes the Euclidian norm. - - - -

Solution theory

Hyperbolic conservation laws, such as @@ -168,10 +153,6 @@ formulate a scheme that ensures that the discrete approximation of $\mathbf{u}(\mathbf{x},t)$ remains in $\mathcal{B}$. - - - -

Variational versus collocation-type discretizations

Following Step-9, Step-12, and Step-33, at this point it might look tempting @@ -236,9 +217,6 @@ illustrates that deal.ii can be used far beyond the context of variational schemes in Hilbert spaces and that a large number of classes, modules and namespaces from deal.ii can be adapted for such purpose. - - -

Description of the scheme

@@ -300,10 +278,17 @@ For the time being let's note that \lambda_{\text{max}} (\mathbf{U}_j^{n}, \mathbf{U}_i^{n}, \textbf{n}_{ji}) \} \|\mathbf{c}_{ij}\|_{\ell^2} $ +<<<<<<< HEAD Before we start with the description of the implementation of this scheme, it is worth saying a thing or two about the "assembly" of this system. Consider for instance a hypothetical pseudo-code, illustrating a possible strategy to compute the solution $\textbf{U}^{n+1}$: +======= +Before we start with the description of the implementation of this scheme, +it is worth saying a thing or two about the "assembly" of this system. +Consider for instance a hypothetical pseudo-code, illustrating a possible +strategy to compute the solution $\textbf{U}^{n+1}$: +>>>>>>> ee60914d6e... reindent, layout changes and address some review comments @f{align*} &\textbf{For } i \in \mathcal{V} \\ @@ -344,24 +329,22 @@ of application of this kind of schemes, also called "edge-based" or "graph-based" finite element schemes (see for instance @cite Rainald2008 for more historical references). +<<<<<<< HEAD This pseudo-code was introduced only to prepare the mindset of the reader for what is going to be presented in the in the next section. The actual implementation described in the next section is somewhat different from what is described in the pseudo-code but shares the same core mentality: we do not loop on cells but rather we loop on the edges of the sparsity graph (hence the name "edge-based" code) in order to assemble the system. +======= +This pseudo-code was introduced only to prepare the mindset of the reader +for what is going to be presented in the in the next section. The actual +implementation described in the next section is somewhat different from +what is described in the pseudo-code but shares the same core mentality: we +do not loop on cells but rather we loop on the edges of the sparsity graph +(hence the name "edge-based" code) in order to assemble the system. +>>>>>>> ee60914d6e... reindent, layout changes and address some review comments - - -

Implementation of the scheme

- - - - - - - - diff --git a/examples/step-69/step-69.cc b/examples/step-69/step-69.cc index 296217c8c7..95a1ec813a 100644 --- a/examples/step-69/step-69.cc +++ b/examples/step-69/step-69.cc @@ -19,7 +19,7 @@ */ // @sect3{Include files} -// The set of include files is quite standard. The most intriguing part +// The set of include files is quite standard. The most intriguing part // is that: either though this code is a "thread and mpi parallel" // we are using neither Trilinos nor PETSC vectors. Actually we are using dealii // distributed vectors la_parallel_vector.h and the regular dealii @@ -828,10 +828,10 @@ namespace Step69 // \frac{\boldsymbol{\nu}_i}{|\boldsymbol{\nu}_i|}$ where // $\boldsymbol{\nu}_i = \sum_{F \subset \text{supp}(\phi_i)} // \sum_{\mathbf{x}_{q,F}} \nu(\mathbf{x}_{q,F}) - // \phi_i(\mathbf{x}_{q,F})$, here: $F \subset \partial \Omega$ denotes - // faces of elements at the boundary of the domain, and $\mathbf{x}_{q,F}$ + // \phi_i(\mathbf{x}_{q,F})$, here: $F \subset \partial \Omega$ denotes + // faces of elements at the boundary of the domain, and $\mathbf{x}_{q,F}$ // are quadrature points on such face. - // Other more sophisticated definitions for $\nu_i$ are + // Other more sophisticated definitions for $\nu_i$ are // possible but none of them have much influence in theory or practice. // We remind the reader that CopyData includes the class member // local_boundary_normal_map in order to store these local @@ -941,10 +941,10 @@ namespace Step69 if (!discretization->finite_element.has_support_on_face(j, f)) continue; - /* Note that "normal" will only represent the contributions - from one of the faces in the support of the shape - function \phi_j. So we cannot normalize this local - contribution right here, we have to take it "as is" and pass + /* Note that "normal" will only represent the contributions + from one of the faces in the support of the shape + function \phi_j. So we cannot normalize this local + contribution right here, we have to take it "as is" and pass it to the copy data routine. */ Tensor<1, dim> normal; if (id == Boundary::slip) @@ -985,8 +985,13 @@ namespace Step69 for (const auto &it : local_boundary_normal_map) { - auto &[normal, id, position] = boundary_normal_map[it.first]; - auto &[new_normal, new_id, new_position] = it.second; + auto &normal = std::get<0>(boundary_normal_map[it.first]); + auto &id = std::get<1>(boundary_normal_map[it.first]); + auto &position = std::get<2>(boundary_normal_map[it.first]); + + const auto &new_normal = std::get<0>(it.second); + const auto &new_id = std::get<1>(it.second); + const auto &new_position = std::get<2>(it.second); normal += new_normal; id = std::max(id, new_id); @@ -1015,30 +1020,30 @@ namespace Step69 // contains a just copy of the matrix cij_matrix. // That's not what we really // want: we have to normalize its entries. In addition, we have not even - // touched the entries of the matrix norm_matrix yet, and the + // touched the entries of the matrix norm_matrix yet, and the // vectors stored in the map // OfflineData::BoundaryNormalMap are not normalized. // // In principle, this is just offline data, it doesn't make much sense // to over-optimize their computation, since their cost will get amortized - // over the many time steps that we are going to use. However, + // over the many time steps that we are going to use. However, // computing/storing the entries of the matrix // norm_matrix and the normalization of nij_matrix - // are perfect to illustrate thread-parallel node-loops: - // - We want to visit every node $i$ in the mesh/sparsity graph, + // are perfect to illustrate thread-parallel node-loops: + // - We want to visit every node $i$ in the mesh/sparsity graph, // - and for every such node we want to visit to every $j$ such that // $\mathbf{c}_{ij} \not \equiv 0$. // - // From an algebraic point of view, this is equivalent to: visiting + // From an algebraic point of view, this is equivalent to: visiting // every row in the matrix (equivalently sparsity // pattern) and for each one of these rows execute a loop on the columns. // Node-loops is a core theme of this tutorial step (see the pseudo-code - // in the introduction) that will repeat over and over again. That's why + // in the introduction) that will repeat over and over again. That's why // this is the right time to introduce them. // // We have the thread paralellization capability // parallel::apply_to_subranges that is somehow more general than the - // WorkStream framework. In particular, it can be used for our + // WorkStream framework. In particular, it can be used for our // node-loops. // This functionality requires four input arguments: // - A begin iterator: indices.begin() @@ -1073,8 +1078,8 @@ namespace Step69 // element schemes when they are properly implemented. // // Finally, we normalize the vector stored in - // OfflineData::BoundaryNormalMap. This operation has - // not been thread paralellized as it would not illustrate any important + // OfflineData::BoundaryNormalMap. This operation has + // not been thread paralellized as it would not illustrate any important // concept. { @@ -1125,22 +1130,22 @@ namespace Step69 /* This is not thread parallelized, too bad! */ for (auto &it : boundary_normal_map) { - auto &[normal, id, _] = it.second; + auto &normal = std::get<0>(it.second); normal /= (normal.norm() + std::numeric_limits::epsilon()); } } // In order to implement reflecting boundary conditions // $\mathbf{m} \cdot \boldsymbol{\nu}_i =0$ (or equivalently $\mathbf{v} - // \cdot \boldsymbol{\nu}_i =0$ ) the vectors $\mathbf{c}_{ij}$ at the + // \cdot \boldsymbol{\nu}_i =0$ ) the vectors $\mathbf{c}_{ij}$ at the // boundary have to be modified as: // // $\mathbf{c}_{ij} += \int_{\partial \Omega} // (\boldsymbol{\nu}_j - \boldsymbol{\nu}(s)) \phi_j \, \mathrm{d}s$ // // Otherwise we will not be able to claim conservation. The ideas repeat - // themselves: we use Workstream in order to compute this correction, most - // of the following code is about the definition of the worker + // themselves: we use Workstream in order to compute this correction, most + // of the following code is about the definition of the worker // local_assemble_system. { @@ -1247,7 +1252,7 @@ namespace Step69 // Now we define the implementation of momentum, // internal_energy, pressure, // speed_of_sound, and f (the flux of the system). - // The functionality of each one of these functions is self-explanatory from + // The functionality of each one of these functions is self-explanatory from // their names. template @@ -1308,15 +1313,15 @@ namespace Step69 return result; } - // The following function, riemann_data_from_state, takes the - // full state $\mathbf{u} = [\rho,\mathbf{m},E]^\top$ defines a new + // The following function, riemann_data_from_state, takes the + // full state $\mathbf{u} = [\rho,\mathbf{m},E]^\top$ defines a new // "projected state" defined as // // $\widetilde{\mathbf{u}} = [\rho, // \mathbf{m} - (\mathbf{m}\cdot \mathbf{n}_{ij})\mathbf{n}_{ij}, // E - \tfrac{(\mathbf{m}\cdot \mathbf{n}_{ij})^2}{2\rho} ]^\top$ // - // Projected states appear naturally when attempting to compute a maximum + // Projected states appear naturally when attempting to compute a maximum // wavespeed appearing in Riemann problems. namespace @@ -1414,7 +1419,7 @@ namespace Step69 return std::max(std::abs(u_i), std::abs(u_j)) + 5. * std::max(a_i, a_j); } - } /* End of namespace dedicated to the computation of the maximum wavespeed */ + } // namespace // Placeholder here. @@ -1721,7 +1726,9 @@ namespace Step69 if (i >= n_locally_owned) continue; - const auto &[normal, id, position] = it->second; + const auto &normal = std::get<0>(it->second); + const auto &id = std::get<1>(it->second); + const auto &position = std::get<2>(it->second); /* Skip constrained degrees of freedom */ if (++sparsity.begin(i) == sparsity.end(i)) @@ -1741,7 +1748,7 @@ namespace Step69 /* On boundary 2 enforce initial conditions: */ - if (id == Boundary::dirichlet) + else if (id == Boundary::dirichlet) { U_i = initial_values->initial_state(position, t + tau_max); } @@ -1845,15 +1852,13 @@ namespace Step69 const auto bnm_it = boundary_normal_map.find(i); if (bnm_it != boundary_normal_map.end()) { - const auto [normal, id, _] = bnm_it->second; + const auto &normal = std::get<0>(bnm_it->second); + const auto &id = std::get<1>(bnm_it->second); + if (id == Boundary::slip) - { - r_i -= 1. * (r_i * normal) * normal; - } + r_i -= 1. * (r_i * normal) * normal; else - { - r_i = 0.; - } + r_i = 0.; } const double m_i = lumped_mass_matrix.diag_element(i);