From 638eebe9ae64b60e6659a762c6a8134f611e5328 Mon Sep 17 00:00:00 2001 From: "Ignacio Tomas (-EXP)" Date: Tue, 11 Feb 2020 19:08:39 -0700 Subject: [PATCH] Minor modifications --- examples/step-69/doc/intro.dox | 86 +++++++++++++--------------------- examples/step-69/step-69.cc | 62 +++++++++++++----------- 2 files changed, 68 insertions(+), 80 deletions(-) diff --git a/examples/step-69/doc/intro.dox b/examples/step-69/doc/intro.dox index 87307e1e92..6196338e54 100644 --- a/examples/step-69/doc/intro.dox +++ b/examples/step-69/doc/intro.dox @@ -34,7 +34,7 @@ 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- and process- (MPI) parallelization; efficient +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. @@ -83,16 +83,17 @@ $\mathbb{f}(\mathbf{u})$ is defined as 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. -For the tutorial we limit the discussion to the class of polytropic ideal gases -for which the pressure is given by +In this tutorial we limit the discussion to the class of polytropic +ideal gases for which the pressure is given by @f{align*} -p = p(\textbf{u}) := (\gamma -1) \Big(E - \frac{\|\textbf{m}\|^2}{2\,\rho} +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 $\|\,.\|$ denotes the Euclidian norm. +specific heats, and $|\cdot|_{\ell^2}$ denotes the Euclidian norm.

Solution theory

@@ -130,11 +131,11 @@ 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 | - \quad - \rho > 0, - \quad - \ E - \tfrac{|\textbf{m}|_{\ell^2}^2}{2 \rho} > 0, - \quad + \ + \rho > 0 \, , + \ + \ E - \tfrac{|\textbf{m}|_{\ell^2}^2}{2 \rho} > 0 \, , + \ s(\mathbf{u}) \geq \min_{x \in \Omega} s(\mathbf{u}_0(\mathbf{x})) \big\}. @f} @@ -180,8 +181,8 @@ $s_h(\cdot,\cdot)_{L^2(\Omega)}$ is some linear stabilization method instance @cite GuermondErn2004 Chapter 5 and references therein). Most time-dependent discretization approaches described in the deal.II tutorials are based on such a (semi-discrete) variational approach. Fundamentally, -from an analysis perspective, variational discretizations are conceived in -order to provide some notion of global (integral) stabiliy, meaning an +from an analysis perspective, variational discretizations are conceived +to provide some notion of global (integral) stabiliy, meaning an estimate of the form @f{align*} @@ -198,7 +199,6 @@ have proven to be some of the best approaches for simulations in the subsonic shockless regime and similarly benign regimes. However, in the transonic and supersonic regime, and shock-hydrodynamics - applications the use of variational schemes might be questionable. In fact, at the time of this writing, most shock-hydrodynamics codes are still firmly grounded on finite volumes methods. The main reason for failure of @@ -243,18 +243,18 @@ spaces $\pmb{\mathbb{V}}_h := \{\mathbb{V}_h\}^{d+2}$. Let $\mathbf{u}_h \mathbb{R}^{d+2}$ and $\phi_i$ is a scalar-valued shape function. @note For simplicity we will consider the usual Lagrange finite elements. -In such this context we let $\{\mathbf{x}_i\}_{i \in \mathcal{V}}$ denote -the set of all support points (see @ref GlossSupport "this glossary -entry"), where $\mathbf{x}_i \in \mathbb{R}^d$. Then each index $i \in +In such context, let $\{\mathbf{x}_i\}_{i \in \mathcal{V}}$ denote +the set of all support points (see @ref GlossSupport "this glossary entry"), +where $\mathbf{x}_i \in \mathbb{R}^d$. Then each index $i \in \mathcal{V}$ uniquely identifies a support point $\mathbf{x}_i$, as well as a scalar-valued shape function $\phi_i$. -With this notation at hand we can define the scheme as +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} + \mathbf{c}_{ij} - d_{ij} \mathbf{U}_j^{n} = \boldsymbol{0} \, , @f} Where @@ -280,9 +280,10 @@ The definition of $\lambda_{\text{max}} (\mathbf{U},\mathbf{V}, definition in order to focus first on some algorithmic and implementational questions. We note that - $m_i$ and $\mathbf{c}_{ij}$ do not evolve in time (provided we keep the - discretization fixed). It thus makes sense to assemble the matrices - once in a so called offline computation and reuse them in every - time step. They are part of what we are going to call off-line data. + discretization fixed). It thus makes sense to assemble these + matrices/vectors once in a so called offline computation and reuse + them in every time step. They are part of what we are going to call + off-line data. - At every time step we have to evaluate $\mathbb{f}(\mathbf{U}_j^{n})$ and $d_{ij} := \max \{ \lambda_{\text{max}} (\mathbf{U}_i^{n},\mathbf{U}_j^{n}, \textbf{n}_{ij}), @@ -297,8 +298,10 @@ $t_n$: @f{align*} &\textbf{for } i \in \mathcal{V} \\ -&\ \ \ \ \{\mathbf{c}_{ij}\}_{j \in \mathcal{I}(i)} \leftarrow \texttt{gather} (\textbf{c}, \mathcal{I}(i)) \\ -&\ \ \ \ \{\textbf{U}_j^n\}_{j \in \mathcal{I}(i)} \leftarrow \texttt{gather} (\textbf{U}^n, \mathcal{I}(i)) \\ +&\ \ \ \ \{\mathbf{c}_{ij}\}_{j \in \mathcal{I}(i)} \leftarrow +\texttt{gather_cij_vectors} (\textbf{c}, \mathcal{I}(i)) \\ +&\ \ \ \ \{\textbf{U}_j^n\}_{j \in \mathcal{I}(i)} \leftarrow +\texttt{gather_state_vectors} (\textbf{U}^n, \mathcal{I}(i)) \\ &\ \ \ \ \ \textbf{U}_i^{n+1} \leftarrow \mathbf{U}_i^{n} \\ &\ \ \ \ \textbf{for } j \in \mathcal{I}(i) \\ &\ \ \ \ \ \ \ \ \texttt{compute } d_{ij} \\ @@ -306,7 +309,7 @@ $t_n$: &\ \ \ \ \ \ \ \ \textbf{U}_i^{n+1} \leftarrow \textbf{U}_i^{n+1} - \frac{\tau_n}{m_i} \mathbb{f}(\mathbf{U}_j^{n})\cdot \mathbf{c}_{ij} + d_{ij} \mathbf{U}_j^{n} \\ &\ \ \ \ \textbf{end} \\ -&\ \ \ \ \texttt{scatter} (\textbf{U}^n, \mathcal{I}(i), \textbf{U}_i^n)) \\ +&\ \ \ \ \texttt{scatter_updated_state} (\textbf{U}_i^{n+1}) \\ &\textbf{end} @f} @@ -315,9 +318,9 @@ We note here that: - Here $\textbf{c}$ and $\textbf{U}^n$ are a global matrix and a global vector containing all the vectors $\mathbf{c}_{ij}$ and all the states $\mathbf{U}_j^n$ respectively. -- $\texttt{gather_cij_vectors}$ and $\texttt{gather_state_vectors}$ are -hypothetical implementations that collect (from global matrices and vectors) -only the quantities required to compute the update at the node $i$. +- $\texttt{gather_cij_vectors}$, $\texttt{gather_state_vectors}$, and +$\texttt{scatter_updated_state}$ are hypothetical implementations that +either collect (from) or write (into) global matrices and vectors. - Note that: if we assume a cartesian mesh in two space 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 @@ -325,34 +328,11 @@ 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() ) only collects values for the -local cell (just a subset of the patch). - -It is worth noting that, from a practitioner's point of view -fully-algebraic schemes (i.e. no bilinear forms, no cell-loops, and no -quadratures) are not unusual at all in the CFD community. There is rich history -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). - -We note that: -- This algorithm does not require any form of quadrature or cell-loops. -- Here, $\textbf{c}$ and $\textbf{U}^n$ are a global matrix and a global vector - containing all the vectors $\mathbf{c}_{ij}$ and all the states - $\mathbf{U}_j^n$, respectively. -- $\texttt{gather}$ and $\texttt{scatter}$ are helper functions - that collect from, or distribute values from global vectors and matrices. -- For an interior node $\mathbf{x}_i$ on a regular mesh in two space - dimensions (assuming a first-order polynomial space $\mathbb{Q}^1$) the - stencil $\mathcal{I}(i)$ contains nine entries. The update for a single - state $\textbf{U}_i^n$ thus depends on nine state-vectors - $\{\textbf{U}_j^n\}_{j \in \mathcal{I}(i)}$ (i.e., all the states in the - patch formed by the support of the shape function $\phi_i$). This is one - of the major differences compared to cell-based loop, where an update - typically only operates on states associated with a single cell. +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 +the time-step size $\tau$ has to be chosen subject to a CFL condition @f{align*} \tau_n = c_{\text{cfl}}\,\min_{ i\in\mathcal{V}}\left(\frac{m_i}{-2\,d_{ii}^{n}}\right), diff --git a/examples/step-69/step-69.cc b/examples/step-69/step-69.cc index 8512f40583..ea5f3e926e 100644 --- a/examples/step-69/step-69.cc +++ b/examples/step-69/step-69.cc @@ -37,7 +37,7 @@ // lac/la_parallel_vector.h. Instead of a Trilinos, or PETSc // specific matrix class, we will use a non-distributed // dealii::SparseMatrix (lac/sparse_matrix.h) to store the local -// part of the $c_{ij}$, $n_{ij}$ and $d_{ij}$ matrices. +// part of the $\mathbf{c}_{ij}$, $\mathbf{n}_{ij}$ and $d_{ij}$ matrices. #include #include #include @@ -95,11 +95,17 @@ // and scratch data object private and make methods and data structures // used by other classes public. // -// @note: A cleaner approach would be to guard access to all data +// @note A cleaner approach would be to guard access to all data // structures by getter/setter // functions. For the sake of brevity, we refrain from that approach, // though. +// +// We also note that the vast majority of classes is derived from +// ParameterAcceptor. This facilitates the population of all the global +// parameters into a single (global) ParameterHandler. More explanations +// about the use inheritance from ParameterAcceptor as a global subscription +// mechanism can be found in Step-59. namespace Step69 { @@ -120,7 +126,7 @@ namespace Step69 // // The class Discretization contains all data structures // concerning the mesh (triangulation) and discretization (mapping, - // finite element, quadrature) of the problem. We use the + // finite element, quadrature) of the problem. As mentioned, we use // ParameterAcceptor class to automatically populate problem-specific // parameters, such as the geometry information // (length, etc.) or the refinement level @@ -130,7 +136,7 @@ namespace Step69 // constructor, and defer the creation of the mesh to the // setup() method that can be called once all parameters are // read-in via ParameterAcceptor::initialize(). - // + template class Discretization : public ParameterAcceptor { @@ -165,10 +171,9 @@ namespace Step69 // // The class OfflineData contains pretty much all components // of the discretization that do not evolve in time, in particular, the - // DoFHandler, SparsityPattern, boundary maps, the lumped mass, $c_{ij}$, - // and $n_{ij}$ matrices. - // - // Here, the term offline refers to the fact that all the class + // DoFHandler, SparsityPattern, boundary maps, the lumped mass, + // $\mathbf{c}_{ij}$ and $\mathbf{n}_{ij}$ matrices. Here, the term + // offline refers to the fact that all the class // members of OfflineData have well-defined values // independent of the current time step. This means that they can be // initialized ahead of time (at time step zero) and are not meant @@ -347,7 +352,7 @@ namespace Step69 private: // We declare a private callback function that will be wired up to the - // ParameterAcceptor::parse_parameters_call_back signal + // ParameterAcceptor::parse_parameters_call_back signal. void parse_parameters_callback(); Tensor<1, dim> initial_direction; @@ -360,15 +365,16 @@ namespace Step69 // classes at hand we can now implement the explicit time-stepping scheme // that was introduced in the discussion above. The main method of the // TimeStep class is step(vector_type &U, double - // t). That takes a reference to a state vector U and - // a time point t as arguments, computes the updated solution, - // stores it in the vector temp, swaps its contents with the - // vector U, and returns the chosen step-size $\tau$. + // t) that takes a reference to a state vector U and + // a time point t (as input arguments) computes the updated + // solution, stores it in the vector temp, swaps its contents + // with the vector U, and returns the chosen step-size + // $\tau$. // // The other important method is prepare() which primarily // sets the proper partition and sparsity pattern for the temporary - // vector temp and the matrix dij_matrix. - // + // vector temp and the matrix dij_matrix + // respectively. template class TimeStep : public ParameterAcceptor @@ -518,7 +524,7 @@ namespace Step69 // The first major task at hand is the typical triplet of grid // generation, setup of data structures, and assembly. A notable novelty // in this example step is the use of the ParameterAcceptor class that we - // use to populate parameter values: We first initialize the + // use to populate parameter values: we first initialize the // ParameterAcceptor class by calling its constructor with a string // subsection denoting the correct subsection in the // parameter file. Then, in the constructor body every parameter value is @@ -566,7 +572,7 @@ namespace Step69 } // Note that in the previous constructor we only passed the MPI - // communicator to the triangulationbut we still have not + // communicator to the triangulation but we still have not // initialized the underlying geometry/mesh. As mentioned earlier, we // have to postpone this task to the setup() function that // gets called after the ParameterAcceptor::initialize() function has @@ -739,19 +745,21 @@ namespace Step69 // "holes" in the index range. The distributed matrices offered by // deal.II avoid this by translating from a global index range into a // contiguous local index range. But this is the precisely the type of - // index manipulation we want to avoid. + // index manipulation we want to avoid in our assembly loops. // - // Lucky enough, the Utilities::MPI::Partitioner used for distributed - // vectors provides exactly what we need: It manages a translation from + // The Utilities::MPI::Partitioner already implements the translation from // a global index range to a contiguous local (per MPI rank) index - // range. We therefore simply create a "local" sparsity pattern for the - // contiguous index range $[0,$n_locally_relevant$)$ and - // translate between global dof indices and the above local range with - // the help of the Utilities::MPI::Partitioner::global_to_local() - // function. All that is left to do is to ensure that we always access + // range (we don't have to reinvent the wheel). We just need to use that + // translation capability (once and only once) in order to create a + // "local" sparsity pattern for + // the contiguous index range $[0,$n_locally_relevant$)$. That + // capability can be invoked by + // Utilities::MPI::Partitioner::global_to_local() + // function. All that is left to do is to ensure that, when implementing + // our scatter and gather auxiliary functions, we always access // elements of a distributed vector by a call to - // LinearAlgebra::distributed::Vector::local_element(). That way we - // avoid index translations altogether. + // LinearAlgebra::distributed::Vector::local_element(). That way we avoid + // index translations altogether and operate exclusively with local indices. { TimerOutput::Scope t( -- 2.39.5