]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Minor modifications
authorIgnacio Tomas (-EXP) <itomas@sandia.gov>
Wed, 12 Feb 2020 02:08:39 +0000 (19:08 -0700)
committerMatthias Maier <tamiko@43-1.org>
Sat, 7 Mar 2020 18:28:55 +0000 (12:28 -0600)
examples/step-69/doc/intro.dox
examples/step-69/step-69.cc

index 87307e1e929b053b082dce84cdbe92744d4b7e45..6196338e54090d2961ca78c74eba4b42a31f8a62 100644 (file)
@@ -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
 <a href="https://en.wikipedia.org/wiki/Heat_capacity_ratio">ratio of
-specific heats</a>, and $\|\,.\|$ denotes the Euclidian norm.
+specific heats</a>, and $|\cdot|_{\ell^2}$ denotes the Euclidian norm.
 
 
 <h4>Solution theory</h4>
@@ -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 <i>offline computation</i> 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 <i>offline computation</i> 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<dim, spacedim>.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<dim, spacedim>.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),
index 8512f40583250ec3cce0d895c59735ac0640a698..ea5f3e926e8abd3fabf906e94a5466fb97b03bdf 100644 (file)
@@ -37,7 +37,7 @@
 // <code>lac/la_parallel_vector.h</code>. Instead of a Trilinos, or PETSc
 // specific matrix class, we will use a non-distributed
 // dealii::SparseMatrix (<code>lac/sparse_matrix.h</code>) 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 <deal.II/base/conditional_ostream.h>
 #include <deal.II/base/parallel.h>
 #include <deal.II/base/parameter_acceptor.h>
 // 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 <a
 // href="https://en.wikipedia.org/wiki/Mutator_method">getter/setter
 // functions</a>. 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 <code>Discretization</code> 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
   // (<code>length</code>, etc.) or the refinement level
@@ -130,7 +136,7 @@ namespace Step69
   // constructor, and defer the creation of the mesh to the
   // <code>setup()</code> method that can be called once all parameters are
   // read-in via ParameterAcceptor::initialize().
-  //
+
   template <int dim>
   class Discretization : public ParameterAcceptor
   {
@@ -165,10 +171,9 @@ namespace Step69
   //
   // The class <code>OfflineData</code> 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 <i>offline</i> 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 
+  // <i>offline</i> refers to the fact that all the class
   // members of <code>OfflineData</code> have well-defined values
   // independent of the current time step. This means that they can be
   // initialized ahead of time (at <i>time step zero</i>) 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
   // <code>TimeStep</code> class is <code>step(vector_type &U, double
-  // t)</code>. That takes a reference to a state vector <code>U</code> and
-  // a time point <code>t</code> as arguments, computes the updated solution,
-  // stores it in the vector <code>temp</code>, swaps its contents with the
-  // vector <code>U</code>, and returns the chosen step-size $\tau$.
+  // t)</code> that takes a reference to a state vector <code>U</code> and
+  // a time point <code>t</code> (as input arguments) computes the updated 
+  // solution, stores it in the vector <code>temp</code>, swaps its contents 
+  // with the vector <code>U</code>, and returns the chosen step-size 
+  // $\tau$.
   //
   // The other important method is <code>prepare()</code> which primarily
   // sets the proper partition and sparsity pattern for the temporary
-  // vector <code>temp</code> and the matrix <code>dij_matrix</code>.
-  //
+  // vector <code>temp</code> and the matrix <code>dij_matrix</code> 
+  // respectively.
 
   template <int dim>
   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
   // <code>subsection</code> 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 <code>triangulation</code>but we still have not
+  // communicator to the <code>triangulation</code> but we still have not
   // initialized the underlying geometry/mesh. As mentioned earlier, we
   // have to postpone this task to the <code>setup()</code> 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,$<code>n_locally_relevant</code>$)$ 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,$<code>n_locally_relevant</code>$)$. 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(

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.