From: Matthias Maier Date: Thu, 26 Dec 2019 15:31:40 +0000 (-0600) Subject: update documentation, part II X-Git-Tag: v9.2.0-rc1~443^2~33 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ff24cb6a8d64ea7f570f506eef5b6f5025556c23;p=dealii.git update documentation, part II --- diff --git a/examples/step-69/doc/intro.dox b/examples/step-69/doc/intro.dox index ec8b1890b6..1989403057 100644 --- a/examples/step-69/doc/intro.dox +++ b/examples/step-69/doc/intro.dox @@ -3,8 +3,6 @@ and Ignacio Tomas (Sandia National Laboratories, Albuquerque). -@dealiiTutorialDOI{10.5281/zenodo.TODO,https://zenodo.org/badge/DOI/10.5281/zenodo.TODO.svg} - @note This tutorial step implements a first-order accurate guaranteed maximum wavespeed method based on a first-order graph viscosity for solving Euler's equations of gas dynamics @cite GuermondPopov2016. As @@ -14,6 +12,8 @@ high-performance implementation of a second-order accurate scheme that uses convex limiting techniques, and strong stability-preserving (SSP) time integration, see @cite GuermondEtAl2018. +@todo Add zenodo link +

Introduction

@@ -38,7 +38,9 @@ as elementary building blocks in higher-order schemes @cite GuermondEtAl2018. However, we hope that the reader still finds the tutorial step to be a good starting point (in particular with respect to the programming techniques) before jumping into full research codes such as -the second-order scheme @cite GuermondEtAl2018 maintained TODO. +the second-order scheme @cite GuermondEtAl2018. + +@todo Add link to repository and project of the second order code. @@ -173,8 +175,8 @@ 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 estimate -of the form +order 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) |\!|\!| @@ -190,15 +192,17 @@ 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 variational -schemes in such extreme regimes is the lack of pointwise stability. This stems -from the fact that a priori bounds on integrated quantities (e.g. -integrals of moments) have in general no implications on pointwise properties -of the solution. While some of these problems might be alleviated by the -(perpetual) chase of the right shock capturing scheme, finite difference-like -and finite volume schemes still have an edge in many regards. + +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 +variational schemes in such extreme regimes is the lack of pointwise +stability. This stems from the fact that a priori bounds on +integrated quantities (e.g. integrals of moments) have in general no +implications on pointwise properties of the solution. While some of these +problems might be alleviated by the (perpetual) chase of the right shock +capturing scheme, finite difference-like and finite volume schemes still +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 @@ -209,13 +213,13 @@ collocation-type scheme) that preserves constraints pointwise, i.e., \;\text{at every node}\;\mathbf{x}_i\;\text{of the mesh}. @f} -However, 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 to always work, all the time, no exception. This -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. +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 +to always work, all the time, no exception. This 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

@@ -232,15 +236,14 @@ spaces $\pmb{\mathbb{V}}_h := \{\mathbb{V}_h\}^{d+2}$. Let $\mathbf{u}_h \in \mathcal{V}} \mathbf{U}_i \phi_i$ where $\mathbf{U}_i \in \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 -context $\{\mathbf{x}_i\}_{i \in \mathcal{V}}$ be the set of all "support -points" (see @ref GlossSupport "this glossary entry") where $\mathbf{x}_i \in -\mathbb{R}^d$. Then each integer index $i \in \mathcal{V}$ -uniquely identifies a support point $\mathbf{x}_i$ and/or scalar-valued shape -function $\phi_i$. +@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 +\mathcal{V}$ uniquely identifies a support point $\mathbf{x}_i$, as well as +a scalar-valued shape function $\phi_i$. -With this notation 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} @@ -266,44 +269,40 @@ Where - $\textbf{n}_{ij} = \frac{\mathbf{c}_{ij}}{ \|\mathbf{c}_{ij}\|_{\ell^2} }$ The definition of $\lambda_{\text{max}} (\mathbf{U},\mathbf{V}, -\textbf{n})$ is far from trivial and we will postpone their definition in -order to focus on the computational/coding issues of this tutorial Step. -For the time being let's note that - - $m_i$ and $\mathbf{c}_{ij}$ do not evolve in time. It makes sense to - compute and store them once, and later recall them at very time step. - They are part of what we are going to call off-line data. +\textbf{n})$ is far from trivial and we will postpone the precise +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. - 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}), \lambda_{\text{max}} (\mathbf{U}_j^{n}, \mathbf{U}_i^{n}, - \textbf{n}_{ji}) \} \|\mathbf{c}_{ij}\|_{\ell^2} $ + \textbf{n}_{ji}) \} \|\mathbf{c}_{ij}\|_{\ell^2} $, which will + constitute the bulk of the computational cost. -<<<<<<< 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 +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)} := -\texttt{gather_cij_vectors}(\textbf{c}, \mathcal{I}(i)) \\ -&\ \ \ \ \{\textbf{U}_j^n\}_{j \in \mathcal{I}(i)} := -\texttt{gather_state_vectors}(\textbf{U}^n, \mathcal{I}(i)) \\ -&\ \ \ \ \ \textbf{U}_i^{n+1} := \mathbf{U}_i^{n} \\ -&\ \ \ \ \textbf{For } j \in \mathcal{I}(i) \\ -&\ \ \ \ \ \ \ \ \textbf{U}_i^{n+1} := \textbf{U}_i^{n+1} - \frac{\tau}{m_i} - \mathbb{f}(\mathbf{U}_j^{n})\cdot - \mathbf{c}_{ij} + d_{ij} \mathbf{U}_j^{n} \\ -&\ \ \ \ \textbf{EndFor} \\ -&\textbf{EndFor} +&\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)) \\ +&\ \ \ \ \ \textbf{U}_i^{n+1} \leftarrow \mathbf{U}_i^{n} \\ +&\ \ \ \ \textbf{for } j \in \mathcal{I}(i) \\ +&\ \ \ \ \ \ \ \ \texttt{compute } d_{ij} \\ +&\ \ \ \ \ \ \ \ \texttt{compute } \mathbb{f}(\mathbf{U}_j^{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)) \\ +&\textbf{end} @f} +<<<<<<< HEAD We note here that: - This "assembly" 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 @@ -329,22 +328,43 @@ 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 +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. + +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 +@f{align*} + \tau_n = c_{\text{cfl}}\,\min_{ + i\in\mathcal{V}}\left(\frac{m_i}{-2\,d_{ii}^{n}}\right), +@f} +where $0edge-based or +graph-based finite element schemes (see for instance @cite +Rainald2008 for a historical overview). +@todo Explain what to do for slip, dirichlet and do-nothing boundary +conditions.

Implementation of the scheme

+@todo Maybe comment on some key features of the implementation?! diff --git a/examples/step-69/step-69.cc b/examples/step-69/step-69.cc index 8823c0f2cb..948fb2dc44 100644 --- a/examples/step-69/step-69.cc +++ b/examples/step-69/step-69.cc @@ -1,6 +1,6 @@ /* --------------------------------------------------------------------- * - * Copyright (C) 2012 - 2019 by the deal.II authors + * Copyright (C) 2019 - 2020 by the deal.II authors * * This file is part of the deal.II library. * @@ -19,93 +19,115 @@ */ // @sect3{Include files} -// 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 -// sparse matrices sparse_matrix.h + +// The set of include files is quite standard. The most intriguing part is +// the fact that we will rely solely on deal.II data structures for MPI +// parallelization, in particular distributed::Triangulation and +// LinearAlgebra::distributed::Vector included through +// distributed/tria.h and +// 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. #include -#include #include #include #include #include #include #include + #include + #include #include #include + #include #include #include #include + #include #include + #include #include #include #include #include + #include + #include #include +// In addition to above deal.II specific includes, we also include four +// boost headers. The first two are for binary archives that we will use +// for implementing a checkpointing and restart mechanism. #include #include -#include + +// The last two boost header files are for creating custom iterator ranges +// over integer intervals. #include #include -#include - - -// @sect3{Declaration/s of the namespace Step69} +// @sect3{Class template declarations} +// Next we declare all data structures and class templates of the example +// step. namespace Step69 { - enum Boundary : dealii::types::boundary_id + using namespace dealii; + + // We start with an enum describing all possible boundary conditions + // encountered in this tutorial step. Such an enum allows us to refer to + // boundary types by a mnemonic (such as + // Boundary::do_nothing) rather than a numerical value. + enum Boundary : types::boundary_id { do_nothing = 0, slip = 2, dirichlet = 3, }; - - // @sect4{Declaration of Discretization class template} - - // The main goal of this class is to digest the input file and act as a - // "container" of members that may be changed/decided at run time (through the - // input file). It was natural to derive this class from - // dealii::ParameterAcceptor. This class is in charge of - // initializing mpi comunicator, geometry dimensions, triangulation, mapping, - // finite element, mapping, and quadratures. If we think of the class - // Discretization as a "container": it doesn't contain any - // memmory demanding class member such a dof_handlers, vectors or matrices. - // The most memmory thirsty class member is the - // dealii::parallel::distributed::Triangulation. - + // @sect4{class Discretization} + // + // The class Discretization contains all data structures + // concerning the mesh (triangulation) and discretization (mapping, + // finite element, quadrature) of the problem. We use the + // ParameterAcceptor class to automatically populate problem-specific + // parameters, such as the geometry information + // (immersed_disc_length, etc.) or the refinement level + // (refinement) from a parameter file. This requires us to + // split the initialization of data structures into two functions: We + // initialize everything that does not depend on parameters in the + // 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 dealii::ParameterAcceptor + class Discretization : public ParameterAcceptor { public: - Discretization(const MPI_Comm & mpi_communicator, - dealii::TimerOutput &computing_timer, - const std::string & subsection = "Discretization"); + Discretization(const MPI_Comm & mpi_communicator, + TimerOutput & computing_timer, + const std::string &subsection = "Discretization"); void setup(); const MPI_Comm &mpi_communicator; - dealii::parallel::distributed::Triangulation triangulation; + parallel::distributed::Triangulation triangulation; - const dealii::MappingQ mapping; - const dealii::FE_Q finite_element; - const dealii::QGauss quadrature; - const dealii::QGauss face_quadrature; + const MappingQ mapping; + const FE_Q finite_element; + const QGauss quadrature; + const QGauss face_quadrature; private: - dealii::TimerOutput &computing_timer; + TimerOutput &computing_timer; double immersed_disc_length; double immersed_disc_height; @@ -115,87 +137,110 @@ namespace Step69 unsigned int refinement; }; - // @sect4{Declaration of OfflineData class template} - - // The class OfflineData is initializes (and "owns") - // pretty much all the components of the discretization that - // do not evolve in time. In particular: dof_handler, sparsity - // patterns, boundary maps, lumped mass matrix, and other matrices - // and vectors that do not change in time are members of this class. - // The term "offline" here refers to the idea that all the class members - // of OfflineData are initialized and assigned values - // a "time step zero" and are not meant to be modified at any other later - // time step. For instance, the sparsity pattern should not - // change as we advance in time (we are not doing any form of adaptivity in - // space). Similarly, the entries of the vector - // lumped_mass_matrix should not be modified as we advance in - // time either. + // @sect4{class OfflineData} // - // Placeholder: Say something about BoundaryNormalMap. + // 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 + // 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 + // to be modified at any other later time step. For instance, the + // sparsity pattern should not change as we advance in time (we are not + // doing any form of adaptivity in space). Similarly, the entries of the + // lumped mass matrix should not be modified as we advance in time + // either. + // + // We also compute and store a boundary_normal_map that + // contains a map from a global index of type `types:global_dof_index` of + // a boundary degree of freedom to a tuple consisting of a normal vector, + // the boundary id, and the position associated with the degree of + // freedom. We actually have to compuate and store this geometric + // information in this class because we won't have access to geometric + // (or cell-based) information later on in the algebraic loops over the + // sparsity pattern. + // + // @note Even though this class currently does not have any parameters + // that could be read in from a parameter file we nevertheless dervie + // from ParameterAcceptor and follow the same idiom of providing a + // setup() (and assemble()) method as for the + // class Discretization. template - class OfflineData : public dealii::ParameterAcceptor + class OfflineData : public ParameterAcceptor { public: using BoundaryNormalMap = - std::map /* normal */, - dealii::types::boundary_id /* id */, - dealii::Point> /* position */>; + std::map, types::boundary_id, Point>>; OfflineData(const MPI_Comm & mpi_communicator, - dealii::TimerOutput & computing_timer, + TimerOutput & computing_timer, const Discretization &discretization, const std::string & subsection = "OfflineData"); void setup(); void assemble(); - dealii::DoFHandler dof_handler; + DoFHandler dof_handler; - std::shared_ptr partitioner; + std::shared_ptr partitioner; unsigned int n_locally_owned; unsigned int n_locally_relevant; - dealii::SparsityPattern sparsity_pattern; + SparsityPattern sparsity_pattern; BoundaryNormalMap boundary_normal_map; - dealii::SparseMatrix lumped_mass_matrix; - std::array, dim> cij_matrix; - std::array, dim> nij_matrix; - dealii::SparseMatrix norm_matrix; + SparseMatrix lumped_mass_matrix; + std::array, dim> cij_matrix; + std::array, dim> nij_matrix; + SparseMatrix norm_matrix; private: - const MPI_Comm & mpi_communicator; - dealii::TimerOutput &computing_timer; + const MPI_Comm &mpi_communicator; + TimerOutput & computing_timer; - dealii::SmartPointer> discretization; + SmartPointer> discretization; }; - // @sect4{Declaration of ProblemDescription class template} - - // Most of the implementations of the members of this class will be utility - // classes/functions specific for Euler's equations: - // - The type alias rank1_type will be used for the states - // $\mathbf{U}_i^n$ - // - The type alias rank2_type will be used for the fluxes - // $\mathbb{f}(\mathbf{U}_j^n)$. - // - The implementation of momentum will extract $\textbf{m}$ - // (out of the state vector $[\rho,\textbf{m},E]$) and store it in a - // Tensor<1, dim> for our convenience. - // - The implementation of internal_energy will compute - // $E - \frac{|\textbf{m}|^2}{2\rho}$ from the state vector + // @sect4{class ProblemDescription} + // + // The member functions of this class are utility functions 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 + // $\mathbb{f}(\mathbf{U}_j^n)$. + // - The momentum function extracts $\textbf{m}$ + // out of the state vector $[\rho,\textbf{m},E]$) and stores it in a + // Tensor<1, dim>. + // - The internal_energy function computes $E - + // \frac{|\textbf{m}|^2}{2\rho}$ from a given state vector // $[\rho,\textbf{m},E]$. // - // The purpose of the remaining class members component_names, - // pressure, and speed_of_sound, - // is evident from their names. Most notably, the last - // one compute_lambda_max is in charge of computing - // $\lambda_{max}(\mathbf{U},\mathbf{V},\mathbf{n})$ which is required - // to compute the first order viscosity $d_{ij}$ as detailed in the section - // Description of the scheme. + // The purpose of the class members component_names, + // pressure, and speed_of_sound, is evident + // from their names. We also provide a function + // compute_lambda_max, that computes the wave speed estimate + // mentioned above, $\lambda_{max}(\mathbf{U},\mathbf{V},\mathbf{n})$, + // which is used in the computation of the $d_{ij}$ matrix. + // + // @note The DEAL_II_ALWAYS_INLINE macro expands to a + // (compiler specific) pragma that ensures that the corresponding + // function defined in this class is always inlined, i.e., the function + // body is put in place for every invocation of the function, and no call + // (and code indirection) is generated. This is stronger than the + // inline keyword, which is more or less a (mild) suggestion + // to the compiler that the programmer things it would be beneficial to + // inline the function. DEAL_II_ALWAYS_INLINE should only be + // used rarely and with caution in situations such as this one, where we + // actually know (due to benchmarking) that inlining the function in + // question actually improves performance. template class ProblemDescription @@ -203,15 +248,14 @@ namespace Step69 public: static constexpr unsigned int problem_dimension = 2 + dim; - using rank1_type = dealii::Tensor<1, problem_dimension>; - using rank2_type = - dealii::Tensor<1, problem_dimension, dealii::Tensor<1, dim>>; + using rank1_type = Tensor<1, problem_dimension>; + using rank2_type = Tensor<1, problem_dimension, Tensor<1, dim>>; const static std::array component_names; static constexpr double gamma = 7. / 5.; - static DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, dim> + static DEAL_II_ALWAYS_INLINE inline Tensor<1, dim> momentum(const rank1_type U); static DEAL_II_ALWAYS_INLINE inline double @@ -225,34 +269,43 @@ namespace Step69 static DEAL_II_ALWAYS_INLINE inline rank2_type f(const rank1_type U); static DEAL_II_ALWAYS_INLINE inline double - compute_lambda_max(const rank1_type U_i, - const rank1_type U_j, - const dealii::Tensor<1, dim> &n_ij); + compute_lambda_max(const rank1_type U_i, + const rank1_type U_j, + const Tensor<1, dim> &n_ij); }; - // @sect4{Declaration of InitialValues class template} - - // Placeholder here + // @sect4{class InitialValues} + // + // The class InitialValues only public data type is a + // function initial_state that computes the initial state of + // a given point and time. For the purpose of this example step we simply + // implement a homogeneous uniform flow field for which the direction and + // a 1D primitive state (density, velocity, pressure) are read from the + // parameter file. + // + // Instead of implementing yet another setup() function we + // use a callback function parse_parameters_callback that + // can be hooked up to corresponding signal of the ParameterAcceptor, + // ParameterAcceptor::parse_parameters_call_back.connect(...). template - class InitialValues : public dealii::ParameterAcceptor + class InitialValues : public ParameterAcceptor { public: using rank1_type = typename ProblemDescription::rank1_type; InitialValues(const std::string &subsection = "InitialValues"); - std::function &point, double t)> - initial_state; + std::function &point, double t)> initial_state; private: void parse_parameters_callback(); - dealii::Tensor<1, dim> initial_direction; - dealii::Tensor<1, 3> initial_1d_state; + Tensor<1, dim> initial_direction; + Tensor<1, 3> initial_1d_state; }; - // @sect4{Declaration of TimeStep class template} + // @sect4{class TimeStep} // Placeholder here @@ -393,14 +446,7 @@ namespace Step69 vector_type output_vector; }; -} // namespace Step69 - - -// @sect3{Implementation of the classes in namespace Step69} - -namespace Step69 -{ - using namespace dealii; + // @sect3{Implementation of the classes in namespace Step69} // @sect4{Implementation of the members of the class Discretization}