From: Matthias Maier Date: Wed, 19 Feb 2020 22:18:23 +0000 (-0600) Subject: address some more review comments X-Git-Tag: v9.2.0-rc1~443^2~10 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ee4c2c6a46c16ce0908bccbc95306574fca6f8e1;p=dealii.git address some more review comments --- diff --git a/examples/step-69/doc/intro.dox b/examples/step-69/doc/intro.dox index 6fc90ea6e0..239c97cc26 100644 --- a/examples/step-69/doc/intro.dox +++ b/examples/step-69/doc/intro.dox @@ -23,6 +23,7 @@ time integration, see @cite GuermondEtAl2018. @dealiiTutorialDOI{10.5281/zenodo.3643899,https://zenodo.org/badge/DOI/10.5281/zenodo.3643899.svg} +ratio of -specific heats, and $|\cdot|_{\ell^2}$ denotes the Euclidian norm. +where the factor $\gamma \in (1,5/3]$ denotes the ratio of specific +heats.

Solution theory

@@ -111,10 +112,15 @@ solution is formally defined as $\mathbf{u} := \lim_{\epsilon \rightarrow of the parabolic regularization @f{align} \mathbf{u}_t^{\epsilon} + \text{div} \, \mathbb{f}(\mathbf{u}^{\epsilon}) -= {\epsilon} \Delta \mathbf{u}^{\epsilon}. +- {\epsilon} \Delta \mathbf{u}^{\epsilon} = 0. @f} Such solutions, which are understood as the solution recovered in the zero-viscosity limit, are often refered to as viscosity solutions. +(This is, because physically $\epsilon$ can be understood as related to the viscosity of the +fluid, i.e., a quantity that indicates the amount of friction neighboring gas particles moving at +different speeds exert on each other. The Euler equations themselves are derived under +the assumption of no friction, but can physically be expected to describe the limiting +case of vanishing friction or viscosity.) Global existence and uniqueness of such solutions is an 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 @@ -125,7 +131,7 @@ all $\mathbf{x} \in \Omega$ and $t \geq 0$ where \ \rho > 0 \, , \ - \ E - \tfrac{|\textbf{m}|_{\ell^2}^2}{2 \rho} > 0 \, , + \ E - \tfrac{|\textbf{m}|^2}{2 \rho} > 0 \, , \ s(\mathbf{u}) \geq \min_{x \in \Omega} s(\mathbf{u}_0(\mathbf{x})) \big\}. @@ -145,13 +151,15 @@ be satisfied at every point of the domain, not just in an averaged In context of a numerical approximation, a violation of such a constraint has dire consequences: it almost surely leads to catrastrophic failure of -the numerical scheme; loss of hyperbolicity, and overall, loss of -well-posedness of the (discrete) problem. In the following we will +the numerical scheme, loss of hyperbolicity, and overall, loss of +well-posedness of the (discrete) problem. It would also mean that we have computed +something that can not be interpreted physically. (For example, what are we to make +of a computed solution with a negative density?) In the following we will formulate a scheme that ensures that the discrete approximation of $\mathbf{u}(\mathbf{x},t)$ remains in $\mathcal{B}$. -

Variational versus collocation-type discretizations

+

Variational versus collocation-type discretizations

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 @@ -232,29 +240,31 @@ 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: -@f{align*} +With this notation at hand we can define the (explicit time stepping) +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} - \sum_{j \in \mathcal{I}(i)} + d_{ij} \mathbf{U}_j^{n} = \boldsymbol{0} \, , @f} where - - $m_i := \int_{\Omega} \phi_i \, \mathrm{d}\mathbf{x}$ + - $m_i \dealcoloneq \int_{\Omega} \phi_i \, \mathrm{d}\mathbf{x}$ - $\tau$ is the time step size - - $\mathbf{c}_{ij} := \int_{\Omega} \nabla\phi_j\phi_i \, + - $\mathbf{c}_{ij} \dealcoloneq \int_{\Omega} \nabla\phi_j\phi_i \, \mathrm{d}\mathbf{x}$ (note that $\mathbf{c}_{ij}\in \mathbb{R}^d$) - - $\mathcal{I}(i) := \{j \in \mathcal{V} \ | \ \mathbf{c}_{ij} \not \equiv + - $\mathcal{I}(i) \dealcoloneq \{j \in \mathcal{V} \ | \ \mathbf{c}_{ij} \not \equiv \boldsymbol{0}\} \cup \{i\}$. We will refer to $\mathcal{I}(i)$ as the "stencil" (or adjacency list) at the support point $i$. - $\mathbb{f}(\mathbf{U}_j^{n})$ is just the flux $\mathbb{f}$ of the hyperbolic system evaluated at the state $\mathbf{U}_j^{n}$ stored at the support point $j$. - - $d_{ij} := \max \{ \lambda_{\text{max}} + - $d_{ij} \dealcoloneq \max \{ \lambda_{\text{max}} (\mathbf{U}_i^{n},\mathbf{U}_j^{n}, \textbf{n}_{ij}), \lambda_{\text{max}} (\mathbf{U}_j^{n}, \mathbf{U}_i^{n}, - \textbf{n}_{ji}) \} \|\mathbf{c}_{ij}\|_{\ell^2}$ if $i \not = j$ + \textbf{n}_{ji}) \} \|\mathbf{c}_{ij}\|$ if $i \not = j$ - $d_{ii} = - \sum_{j \in \mathcal{I}(i)\backslash \{i\}} d_{ij}$ - - $\textbf{n}_{ij} = \frac{\mathbf{c}_{ij}}{ \|\mathbf{c}_{ij}\|_{\ell^2} }$ + - $\textbf{n}_{ij} = \frac{\mathbf{c}_{ij}}{ \|\mathbf{c}_{ij}\| }$ The definition of $\lambda_{\text{max}} (\mathbf{U},\mathbf{V}, \textbf{n})$ is far from trivial and we will postpone the precise @@ -266,10 +276,10 @@ questions. We note that 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}} + $d_{ij} \dealcoloneq \max \{ \lambda_{\text{max}} (\mathbf{U}_i^{n},\mathbf{U}_j^{n}, \textbf{n}_{ij}), \lambda_{\text{max}} (\mathbf{U}_j^{n}, \mathbf{U}_i^{n}, - \textbf{n}_{ji}) \} \|\mathbf{c}_{ij}\|_{\ell^2} $, which will + \textbf{n}_{ji}) \} \|\mathbf{c}_{ij}\| $, which will constitute the bulk of the computational cost. Consider the following pseudo-code, illustrating a possible straight @@ -301,7 +311,7 @@ $\mathbf{U}_j^n$ respectively. - $\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 +- 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 of the domain) then: $\{\textbf{U}_j^n\}_{j \in \mathcal{I}(i)}$ should contain @@ -333,16 +343,17 @@ application of this kind of schemes, also called edge-based or

Stable boundary conditions and conservation properties.

In the example considered in this tutorial step we use three different types of -boundary conditions: essential-like boundary conditions (we prescribe a state in -the left portion of our domain), outflow boundary conditions (also called -"do-nothing" boundary conditions) at the right boundary of the domain, and -"reflecting" boundary conditions $\mathbf{m} \cdot \boldsymbol{\nu} = 0$ (also -called "slip" boundary conditions) at the top, bottom, and surface of the -obstacle. We will not discuss much about essential and "do-nothing" boundary -conditions since their implementation is relatively easy and the reader will be -able to pick-up the implementation directly from the (documented) source code. -In this portion of the introduction we will focus only on the "reflecting" -boundary conditions which are somewhat more tricky. +boundary conditions: essential-like boundary conditions (we prescribe a +state at the left boundary of our domain), outflow boundary conditions +(also called "do-nothing" boundary conditions) at the right boundary of the +domain, and "reflecting" boundary conditions $\mathbf{m} \cdot +\boldsymbol{\nu} = 0$ (also called "slip" boundary conditions) at the top, +bottom, and surface of the obstacle. We will not discuss much about +essential and "do-nothing" boundary conditions since their implementation +is relatively easy and the reader will be able to pick-up the +implementation directly from the (documented) source code. In this portion +of the introduction we will focus only on the "reflecting" boundary +conditions which are somewhat more tricky. @note At the time of this writing (early 2020) it is not unreasonable to say that both analysis and implementation of stable boundary conditions for @@ -362,7 +373,7 @@ For the case of the reflecting boundary conditions we will proceed as follows: reflecting boundary conditions strongly in a post-processing step where we execute the projection @f{align*} - \mathbf{m}_i := \mathbf{m}_i - (\boldsymbol{\nu}_i \cdot \mathbf{m}_i) + \mathbf{m}_i \dealcoloneq \mathbf{m}_i - (\boldsymbol{\nu}_i \cdot \mathbf{m}_i) \boldsymbol{\nu}_i \ \ \text{for all }\mathbf{x}_i \in \partial\Omega^r @f} that removes the normal component of $\mathbf{m}$. Here the definition of @@ -396,7 +407,7 @@ Dirichlet-like boundary conditions is stable under CFL conditions and does not introduce any loss in accuracy. If $\mathbf{u}_t + \text{div} \, \mathbb{f}(\mathbf{u}) = \boldsymbol{0}$ -represents Euler's equation with reflecting boundary conditions in the entirety +represents Euler's equation with reflecting boundary conditions on the entirety of the boundary (i.e. $\partial\Omega^r \equiv \partial\Omega$) and we integrate in space and time $\int_{\Omega}\int_{t_1}^{t_2}$ we would obtain @f{align*} @@ -425,5 +436,3 @@ boundary} This consistent modification of the $\mathbf{c}_{ij}$ is a direct consequence of simple integration by parts arguments, see page 12 of @cite GuermondEtAl2018 for more details. - - diff --git a/examples/step-69/doc/results.dox b/examples/step-69/doc/results.dox index 9b07843d81..b72a09cd3a 100644 --- a/examples/step-69/doc/results.dox +++ b/examples/step-69/doc/results.dox @@ -146,3 +146,6 @@ TimeLoop::output(t = 4.00006, checkpoint = 1) And with the following result: + +That's substantially better, although of course at the price of having run +the code for roughly 6 hours on 16 cores. diff --git a/examples/step-69/step-69.cc b/examples/step-69/step-69.cc index 6cd9cb4bd8..af53f0cbd4 100644 --- a/examples/step-69/step-69.cc +++ b/examples/step-69/step-69.cc @@ -90,10 +90,10 @@ // structures, and parameters into individual classes. A single class thus // usually centers around either a single data structure (such as the // Triangulation) in the Discretization class, or a single -// method (such as the step() function of the +// method (such as the make_one_step() function of the // TimeStep class). We typically declare parameter variables // and scratch data object `private` and make methods and data structures -// used by other classes public. +// used by other classes `public`. // // @note A cleaner approach would be to guard access to all data // structures by Boundary::do_nothing) rather than a numerical value. - enum Boundary : types::boundary_id - { - do_nothing = 0, - slip = 1, - dirichlet = 2, - }; + // We start with defining a number of types::boundary_id constants used + // throughout the tutorial step. This allows us to refer to boundary + // types by a mnemonic (such as do_nothing) rather than a + // numerical value. + + constexpr types::boundary_id do_nothing = 0; + constexpr types::boundary_id slip = 1; + constexpr types::boundary_id dirichlet = 2; // @sect4{The Discretization class} // @@ -139,13 +137,13 @@ namespace Step69 class Discretization : public ParameterAcceptor { public: - Discretization(const MPI_Comm & mpi_communicator, + Discretization(const MPI_Comm mpi_communicator, TimerOutput & computing_timer, const std::string &subsection = "Discretization"); void setup(); - const MPI_Comm &mpi_communicator; + const MPI_Comm mpi_communicator; parallel::distributed::Triangulation triangulation; @@ -203,7 +201,7 @@ namespace Step69 std::map, types::boundary_id, Point>>; - OfflineData(const MPI_Comm & mpi_communicator, + OfflineData(const MPI_Comm mpi_communicator, TimerOutput & computing_timer, const Discretization &discretization, const std::string & subsection = "OfflineData"); @@ -228,16 +226,16 @@ namespace Step69 SparseMatrix norm_matrix; private: - const MPI_Comm &mpi_communicator; - TimerOutput & computing_timer; + const MPI_Comm mpi_communicator; + TimerOutput & computing_timer; SmartPointer> discretization; }; // @sect4{The ProblemDescription class} // - // The member functions of this class are utility functions specific to - // Euler's equations: + // The member functions of this class are utility functions and data + // structures specific to Euler's equations: // - The type alias rank1_type is used for the states // $\mathbf{U}_i^n$ // - The type alias rank2_type is used for the fluxes @@ -268,24 +266,15 @@ namespace Step69 // actually know (due to benchmarking) that inlining the function in // question improves performance. // - // Finally, we observe that: - // - This is the only class in this tutorial step that is tied to a - // particular "physics" or "hyperbolic conservation law" (in this - // case Euler's equations). All the other classes are primarily - // "discretization" classes, very much agnostic of the particular physics - // being solved. - // - This is a "pure static" class (the antithesis of a - // "pure virtual" class). It's just a convenient way to wrap-up a - // collection of related methods into a single object. Note that we will - // be able to invoke such methods without without creating an instance of - // the class. Similarly, we will not have to provide a constructor - // for this class. + // Finally, we observe that this is the only class in this tutorial step + // that is tied to a particular "physics" or "hyperbolic conservation + // law" (in this case Euler's equations). All the other classes are + // primarily "discretization" classes, very much agnostic of the + // particular physics being solved. template class ProblemDescription { public: - /* constexpr tells the compiler to evaluate "2 + dim" just once at compile - time rather than everytime problem_dimension is invoked. */ static constexpr unsigned int problem_dimension = 2 + dim; using rank1_type = Tensor<1, problem_dimension>; @@ -332,8 +321,8 @@ namespace Step69 // another setup() method to be called (by-hand) after the // call to ParameterAcceptor::initialize() we provide an // "implementation" for the class member - // parse_parameters_call_back which is automatically called when - // invoking ParameterAcceptor::initialize() for every class + // parse_parameters_call_back() which is automatically + // called when invoking ParameterAcceptor::initialize() for every class // that inherits from ParameterAceptor. template class InitialValues : public ParameterAcceptor @@ -359,12 +348,12 @@ namespace Step69 // With the OfflineData and ProblemDescription // 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 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$. + // TimeStep class is make_one_step(vector_type &U, + // double 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 @@ -380,11 +369,10 @@ namespace Step69 using rank1_type = typename ProblemDescription::rank1_type; using rank2_type = typename ProblemDescription::rank2_type; - typedef std::array, - problem_dimension> - vector_type; + using vector_type = + std::array, problem_dimension>; - TimeStep(const MPI_Comm & mpi_communicator, + TimeStep(const MPI_Comm mpi_communicator, TimerOutput & computing_timer, const OfflineData & offline_data, const InitialValues &initial_values, @@ -392,11 +380,11 @@ namespace Step69 void prepare(); - double step(vector_type &U, double t); + double make_one_step(vector_type &U, double t); private: - const MPI_Comm &mpi_communicator; - TimerOutput & computing_timer; + const MPI_Comm mpi_communicator; + TimerOutput & computing_timer; SmartPointer> offline_data; SmartPointer> initial_values; @@ -411,12 +399,12 @@ namespace Step69 // @sect4{The SchlierenPostprocessor class} // // At its core, the Schlieren class implements the class member - // compute_schlieren. The main purpose of this class member + // compute_schlieren(). The main purpose of this class member // is to compute an auxiliary finite element field // schlieren, that is defined at each node by // \f[ \text{schlieren}[i] = e^{\beta \frac{ |\nabla r_i| // - \min_j |\nabla r_j| }{\max_j |\nabla r_j| - \min_j |\nabla r_j| } }, \f] - // where $r$ can in principle be any scalar quantitiy. In practice + // where $r$ can in principle be any scalar quantity. In practice // though, the density is a natural candidate, viz. $r := \rho$. // Schlieren // postprocessing is a standard method for enhancing the contrast of a @@ -435,7 +423,7 @@ namespace Step69 std::array, problem_dimension>; SchlierenPostprocessor( - const MPI_Comm & mpi_communicator, + const MPI_Comm mpi_communicator, TimerOutput & computing_timer, const OfflineData &offline_data, const std::string & subsection = "SchlierenPostprocessor"); @@ -447,8 +435,8 @@ namespace Step69 LinearAlgebra::distributed::Vector schlieren; private: - const MPI_Comm &mpi_communicator; - TimerOutput & computing_timer; + const MPI_Comm mpi_communicator; + TimerOutput & computing_timer; SmartPointer> offline_data; @@ -472,12 +460,12 @@ namespace Step69 public: using vector_type = typename TimeStep::vector_type; - TimeLoop(const MPI_Comm &mpi_comm); + TimeLoop(const MPI_Comm mpi_communnicator); void run(); private: - vector_type interpolate_initial_values(double t = 0); + vector_type interpolate_initial_values(const double t = 0); void output(const vector_type &U, const std::string &name, @@ -485,7 +473,7 @@ namespace Step69 unsigned int cycle, bool checkpoint = false); - const MPI_Comm & mpi_communicator; + const MPI_Comm mpi_communicator; std::ostringstream timer_output; TimerOutput computing_timer; @@ -504,8 +492,6 @@ namespace Step69 TimeStep time_step; SchlierenPostprocessor schlieren_postprocessor; - std::unique_ptr filestream; - std::thread output_thread; vector_type output_vector; }; @@ -525,7 +511,7 @@ namespace Step69 // ParameterAcceptor class with a call to // ParameterAcceptor::add_parameter(). template - Discretization::Discretization(const MPI_Comm & mpi_communicator, + Discretization::Discretization(const MPI_Comm mpi_communicator, TimerOutput & computing_timer, const std::string &subsection) : ParameterAcceptor(subsection) @@ -538,29 +524,25 @@ namespace Step69 , computing_timer(computing_timer) { length = 4.; - add_parameter("immersed disc - length", - length, - "Immersed disc: length of computational domain"); + add_parameter("length", length, "length of computational domain"); height = 2.; - add_parameter("immersed disc - height", - height, - "Immersed disc: height of computational domain"); + add_parameter("height", height, "height of computational domain"); disc_position = 0.6; - add_parameter("immersed disc - object position", + add_parameter("object position", disc_position, - "Immersed disc: x position of immersed disc center point"); + "x position of immersed disc center point"); disc_diameter = 0.5; - add_parameter("immersed disc - object diameter", + add_parameter("object diameter", disc_diameter, - "Immersed disc: diameter of immersed disc"); + "diameter of immersed disc"); refinement = 5; - add_parameter("initial refinement", + add_parameter("refinement", refinement, - "Initial refinement of the geometry"); + "number of refinement steps of the geometry"); } // Note that in the previous constructor we only passed the MPI @@ -578,6 +560,13 @@ namespace Step69 // mesh generated by GridGenerator::hyper_cube_with_cylindrical_hole(). // We refer to step-49, step-53, and step-54 for an overview how to // create advanced meshes. + // We first create 4 temporary (non distributed) coarse triangulations + // that we stitch together with the GridGenerator::merge_triangulation() + // function. We center the disc at $(0,0)$ with a diameter of + // disc_diameter. The lower left corner of the channel has + // coordinates (-disc_position, -height/2) and + // the upper right corner has (length-disc_position, + // height/2). template void Discretization::setup() { @@ -585,15 +574,6 @@ namespace Step69 triangulation.clear(); - // We first create 4 temporary (non distributed) coarse triangulations - // that we stitch together with the - // GridGenerator::merge_triangulation() function. We center the disc at - // $(0,0)$ with a diameter of disc_diameter. The lower - // left corner of the channel has coordinates - // (-disc_position, -height/2) and the upper - // right corner has (length-disc_position, - // height/2). - Triangulation tria1, tria2, tria3, tria4; GridGenerator::hyper_cube_with_cylindrical_hole( @@ -626,38 +606,38 @@ namespace Step69 // We have to fix up the left edge that is currently located at // $x=-$disc_diameter and has to be shifted to - // $x=-$disc_position: + // $x=-$disc_position. As a last step the boundary has to + // be colorized with do_nothing on the right, + // dirichlet on the left and slip on the + // upper and lower outer boundaries and the obstacle. - for (auto cell : triangulation.active_cell_iterators()) - for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; ++v) - { - auto &vertex = cell->vertex(v); - if (vertex[0] <= -disc_diameter + 1.e-6) - vertex[0] = -disc_position; - } + for (const auto &cell : triangulation.active_cell_iterators()) + { + for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; ++v) + { + auto &vertex = cell->vertex(v); + if (vertex[0] <= -disc_diameter + 1.e-6) + vertex[0] = -disc_position; + } + } - // As a last step the boundary has to be colorized with - // Boundary::do_nothing on the right, - // Boundary::dirichlet on the left and - // Boundary::slip on the upper and lower outer boundaries - // and the obstacle: - for (auto cell : triangulation.active_cell_iterators()) + for (const auto &cell : triangulation.active_cell_iterators()) { for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f) { const auto face = cell->face(f); - if (!face->at_boundary()) - continue; - - const auto center = face->center(); + if (face->at_boundary()) + { + const auto center = face->center(); - if (center[0] > length - disc_position - 1.e-6) - face->set_boundary_id(Boundary::do_nothing); - else if (center[0] < -disc_position + 1.e-6) - face->set_boundary_id(Boundary::dirichlet); - else - face->set_boundary_id(Boundary::slip); + if (center[0] > length - disc_position - 1.e-6) + face->set_boundary_id(do_nothing); + else if (center[0] < -disc_position + 1.e-6) + face->set_boundary_id(dirichlet); + else + face->set_boundary_id(slip); + } } } @@ -670,7 +650,7 @@ namespace Step69 // than initializing the corresponding class members in the // initialization list. template - OfflineData::OfflineData(const MPI_Comm & mpi_communicator, + OfflineData::OfflineData(const MPI_Comm mpi_communicator, TimerOutput & computing_timer, const Discretization &discretization, const std::string & subsection) @@ -704,9 +684,10 @@ namespace Step69 DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant); n_locally_relevant = locally_relevant.n_elements(); - partitioner.reset(new Utilities::MPI::Partitioner(locally_owned, - locally_relevant, - mpi_communicator)); + partitioner = + std::make_shared(locally_owned, + locally_relevant, + mpi_communicator); } const auto dofs_per_cell = discretization->finite_element.dofs_per_cell; @@ -773,7 +754,7 @@ namespace Step69 std::vector dof_indices(dofs_per_cell); - for (auto cell : dof_handler.active_cell_iterators()) + for (const auto &cell : dof_handler.active_cell_iterators()) { if (cell->is_artificial()) continue; @@ -830,7 +811,7 @@ namespace Step69 // syntactic // sugar for otherwise somewhat tedious code. - // The first function we introduce, get_entry, will be + // The first function we introduce, get_entry(), will be // used to read the value stored at the entry pointed by a // SparsityPattern iterator it of matrix. The // function works around a small deficiency in the SparseMatrix @@ -841,7 +822,7 @@ namespace Step69 // to the lack of an interface in the SparseMatrix for accessing the // element directly with a SparsityPattern iterator, we unfortunately // have to create a temporary SparseMatrix iterator. We simply hide - // this in the get_entry function. + // this in the get_entry() function. template DEAL_II_ALWAYS_INLINE inline typename Matrix::value_type @@ -853,8 +834,8 @@ namespace Step69 return matrix_iterator->value(); } - // The set_entry helper is the inverse operation of - // get_value: Given an iterator and a value, it sets the + // The set_entry() helper is the inverse operation of + // get_value(): Given an iterator and a value, it sets the // entry pointed to by the iterator in the matrix. template @@ -868,13 +849,13 @@ namespace Step69 matrix_iterator->value() = value; } - // gather_get_entry: we note that $\mathbf{c}_{ij} \in + // gather_get_entry(): we note that $\mathbf{c}_{ij} \in // \mathbb{R}^d$. If $d=2$ then $\mathbf{c}_{ij} = // [\mathbf{c}_{ij}^1,\mathbf{c}_{ij}^2]^\top$. Which basically implies // that we need one matrix per space dimension to store the // $\mathbf{c}_{ij}$ vectors. Similar observation follows for the // matrix $\mathbf{n}_{ij}$. The purpose of - // gather_get_entry is to retrieve those entries a store + // gather_get_entry() is to retrieve those entries a store // them into a Tensor<1, dim> for our convenience. template @@ -887,17 +868,17 @@ namespace Step69 return result; } - // gather (first interface): this first function + // gather() (first interface): this first function // signature, having three input arguments, will be used to retrieve // the individual components (i,l) of a matrix. The - // functionality of gather_get_entry and - // gather is very much the same, but their context is - // different: the function gather does not rely on an + // functionality of gather_get_entry() and + // gather() is very much the same, but their context is + // different: the function gather() does not rely on an // iterator (that actually knows the value pointed) but rather on the // 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 + // 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 // both $i$ and $j$ lie at the boundary. // @@ -918,7 +899,7 @@ namespace Step69 return result; } - // gather (second interface): this second function + // 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 // Tensor<1,problem_dimension> for our convenience. @@ -933,7 +914,7 @@ namespace Step69 return result; } - // scatter: this function has three input arguments, the + // scatter(): this function has three input arguments, the // 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 @@ -962,18 +943,18 @@ namespace Step69 // definition of // - Scratch data (i.e. input info required to carry out computations): in // this case it is scratch_data. - // - The worker: in the case it is local_assemble_system that + // - The worker: in the case it is local_assemble_system() that // actually computes the local (i.e. current cell) contributions from the // scratch data. // - A copy data: a struct that contains all the local assembly // contributions, in this case CopyData(). // - A copy data routine: in this case it is - // copy_local_to_global in charge of actually coping these + // copy_local_to_global() in charge of actually coping these // local contributions into the global objects (matrices and/or vectors) // // Most the following lines are spent in the definition of the worker - // local_assemble_system and the copy data routine - // copy_local_to_global. There is not much to say about the + // local_assemble_system() and the copy data routine + // copy_local_to_global(). There is not much to say about the // WorkStream framework since the vast majority of ideas are reasonably // well-documented in step-9, step-13 and step-32 among others. // @@ -982,7 +963,7 @@ namespace Step69 // // @f{align*} // \widehat{\boldsymbol{\nu}}_i := - // \frac{\boldsymbol{\nu}_i}{|\boldsymbol{\nu}_i|} \ \text{where} \ + // \frac{\boldsymbol{\nu}_i}{|\boldsymbol{\nu}_i|} \ \text{where} // \boldsymbol{\nu}_i := \sum_{T \subset \text{supp}(\phi_i)} // \sum_{F \subset \partial T \cap \partial \Omega} // \sum_{\mathbf{x}_{q,F}} \nu(\mathbf{x}_{q,F}) @@ -1010,7 +991,9 @@ namespace Step69 unsigned int dofs_per_cell = discretization->finite_element.dofs_per_cell; unsigned int n_q_points = discretization->quadrature.size(); - /* This is the implementation of the scratch data required by WorkStream */ + // What follows is the implementation of the scratch data required by + // WorkStream + MeshWorker::ScratchData scratch_data( discretization->mapping, discretization->finite_element, @@ -1025,7 +1008,6 @@ namespace Step69 computing_timer, "offline_data - assemble lumped mass matrix, and c_ij"); - /* This is the implementation of the "worker" required by WorkStream */ const auto local_assemble_system = [&](const auto &cell, auto & scratch, auto & copy) { @@ -1056,8 +1038,8 @@ namespace Step69 return partitioner->global_to_local(index); }); - /* We compute the local contributions for the lumped mass - matrix entries m_i and and vectors c_ij */ + // We compute the local contributions for the lumped mass matrix + // entries m_i and and vectors c_ij in the usual fashion: for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { const auto JxW = fe_values.JxW(q_point); @@ -1079,9 +1061,9 @@ namespace Step69 } /* j */ } /* q */ - /* Now we have to compute the boundary normals. Note that the - following loop does not do much unless the element - has faces on the boundary of the domain */ + // Now we have to compute the boundary normals. Note that the + // following loop does not do much unless the element has faces on + // the boundary of the domain. for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) { const auto face = cell->face(f); @@ -1100,14 +1082,15 @@ 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", store - it and pass it to the copy data routine. The proper - normalization requires an additional loop on nodes.*/ + // 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", + // store it and pass it to the copy data routine. The + // proper normalization requires an additional loop on + // nodes. Tensor<1, dim> normal; - if (id == Boundary::slip) + if (id == slip) { for (unsigned int q = 0; q < n_face_q_points; ++q) normal += fe_face_values.normal_vector(q) * @@ -1132,7 +1115,8 @@ namespace Step69 } }; - /* This is the copy data routine for WorkStream */ + // Last, we provide a copy_local_to_global function as required for + // the WorkStream const auto copy_local_to_global = [&](const auto ©) { const auto &is_artificial = copy.is_artificial; const auto &local_dof_indices = copy.local_dof_indices; @@ -1165,7 +1149,7 @@ namespace Step69 cij_matrix[k].add(local_dof_indices, cell_cij_matrix[k]); nij_matrix[k].add(local_dof_indices, cell_cij_matrix[k]); } - }; /* end of the copy data routine */ + }; WorkStream::run(dof_handler.begin_active(), dof_handler.end(), @@ -1329,7 +1313,7 @@ namespace Step69 // // 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 local_assemble_system. + // of the worker local_assemble_system(). { TimerOutput::Scope t(computing_timer, @@ -1370,7 +1354,7 @@ namespace Step69 if (!face->at_boundary()) continue; - if (id != Boundary::slip) + if (id != slip) continue; const auto &fe_face_values = scratch.reinit(cell, f); @@ -1853,7 +1837,7 @@ namespace Step69 // any surprising code: template - TimeStep::TimeStep(const MPI_Comm & mpi_communicator, + TimeStep::TimeStep(const MPI_Comm mpi_communicator, TimerOutput & computing_timer, const OfflineData & offline_data, const InitialValues &initial_values, @@ -1894,7 +1878,7 @@ namespace Step69 // state U in place and return the chosen time-step size. template - double TimeStep::step(vector_type &U, double t) + double TimeStep::make_one_step(vector_type &U, double t) { // Declare a number of read-only references to various different // variables and data structures. We do this is mainly to have shorter @@ -2209,7 +2193,7 @@ namespace Step69 // On slip boundaries we remove the normal component of the // momentum: - if (id == Boundary::slip) + if (id == slip) { auto m = ProblemDescription::momentum(U_i); m -= 1. * (m * normal) * normal; @@ -2219,7 +2203,7 @@ namespace Step69 // On Dirichlet boundaries we enforce initial conditions // strongly: - else if (id == Boundary::dirichlet) + else if (id == dirichlet) { U_i = initial_values->initial_state(position, t + tau_max); } @@ -2261,7 +2245,7 @@ namespace Step69 template SchlierenPostprocessor::SchlierenPostprocessor( - const MPI_Comm & mpi_communicator, + const MPI_Comm mpi_communicator, TimerOutput & computing_timer, const OfflineData &offline_data, const std::string & subsection /*= "SchlierenPostprocessor"*/) @@ -2321,7 +2305,7 @@ namespace Step69 // positive function such as // $\omega_i(\mathbf{x}) \equiv 1$ (that would allow us to recover the usual // notion of mean value). But as usual, the goal is to reuse the off-line - // data as much as it could be possible. In sense this, the most natural + // data as much as possible. In this sense, the most natural // choice of weight is $\omega_i = \phi_i$. Inserting this choice of // weight and the expansion $r_h(\mathbf{x}) = \sum_{j \in \mathcal{V}} // r_j \phi_j(\mathbf{x})$ into $\mathbf{(*)}$ we get : @@ -2420,7 +2404,7 @@ namespace Step69 const auto &normal = std::get<0>(bnm_it->second); const auto &id = std::get<1>(bnm_it->second); - if (id == Boundary::slip) + if (id == slip) r_i -= 1. * (r_i * normal) * normal; else r_i = 0.; @@ -2507,9 +2491,9 @@ namespace Step69 // restart from an interrupted computation, or not. template - TimeLoop::TimeLoop(const MPI_Comm &mpi_comm) + TimeLoop::TimeLoop(const MPI_Comm mpi_communicator) : ParameterAcceptor("A - TimeLoop") - , mpi_communicator(mpi_comm) + , mpi_communicator(mpi_communicator) , computing_timer(mpi_communicator, timer_output, TimerOutput::never, @@ -2695,9 +2679,10 @@ namespace Step69 // and then perform a single forward Euler step. Note that the // state vector U is updated in place and that - // time_step.step() return the chosen step size. + // time_step.make_one_step() returns the chosen step + // size. - t += time_step.step(U, t); + t += time_step.make_one_step(U, t); // Post processing, generating output and writing out the current // state is a CPU and IO intensive task that we cannot afford to do @@ -2726,7 +2711,7 @@ namespace Step69 template typename TimeLoop::vector_type - TimeLoop::interpolate_initial_values(double t) + TimeLoop::interpolate_initial_values(const double t) { pcout << "TimeLoop::interpolate_initial_values(t = " << t << ")" << std::endl; @@ -2782,9 +2767,9 @@ namespace Step69 template void TimeLoop::output(const typename TimeLoop::vector_type &U, const std::string & name, - double t, - unsigned int cycle, - bool checkpoint) + const double t, + const unsigned int cycle, + const bool checkpoint) { pcout << "TimeLoop::output(t = " << t << ", checkpoint = " << checkpoint << ")" << std::endl;