From: Ignacio Tomas (-EXP) Date: Sun, 16 Feb 2020 08:35:01 +0000 (-0700) Subject: Almost done X-Git-Tag: v9.2.0-rc1~443^2~12 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=98c4404cf1e18f0bb2f53c01b3258e2f00d06561;p=dealii.git Almost done --- diff --git a/examples/step-69/doc/intro.dox b/examples/step-69/doc/intro.dox index e24b23a18e..6fc90ea6e0 100644 --- a/examples/step-69/doc/intro.dox +++ b/examples/step-69/doc/intro.dox @@ -31,7 +31,7 @@ Euler's equations that is based on three ingredients: a of finite elements; a graph-viscosity stabilization based on a guaranteed upper bound of the local wave speed; and explicit 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 +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-MPI parallelization; efficient local numbering @@ -45,7 +45,7 @@ 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. +the second-order scheme discussed in @cite GuermondEtAl2018. @@ -77,7 +77,7 @@ $\mathbb{f}(\mathbf{u})$ is defined as @f} 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. +$p$ that, in general, is defined by a closed-form equation of state. In this tutorial we limit the discussion to the class of polytropic ideal gases for which the pressure is given by @f{align*} @@ -115,7 +115,7 @@ of the parabolic regularization @f} Such solutions, which are understood as the solution recovered in the zero-viscosity limit, are often refered to as viscosity solutions. -Global existence and uniqueness of such solutions is a widely open issue. +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 all $\mathbf{x} \in \Omega$ and $t \geq 0$ where @@ -165,11 +165,11 @@ formulation: Here, $\mathbb{V}_h$ is an appropriate finite element space, and $s_h(\cdot,\cdot)_{L^2(\Omega)}$ is some linear stabilization method (possibly complemented with some ad-hoc shock-capturing technique, see for -instance @cite GuermondErn2004 Chapter 5 and references therein). Most +instance Chapter 5 of @cite GuermondErn2004 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 -to provide some notion of global (integral) stabiliy, meaning an +to provide some notion of global (integral) stability, meaning an estimate of the form @f{align*} |\!|\!| \mathbf{u}_{h}(t) |\!|\!| \leq |\!|\!| \mathbf{u}_{h}(0) |\!|\!| @@ -181,7 +181,7 @@ conservation laws have been very popular since the mid eighties, in particular combined with SUPG-type stabilization and/or upwinding techniques (see the early work of @cite Brooks1982 and @cite Johnson1986). They have proven to be some of the best approaches for simulations in the subsonic -shockless regime and similarly benign regimes. +shockless regime and similarly benign situations. However, in the transonic and supersonic regime, and shock-hydrodynamics applications the use of variational schemes might be questionable. In fact, @@ -204,10 +204,10 @@ collocation-type scheme) that preserves constraints pointwise, i.e., @f} 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 +works on 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 +and that a large number of classes, modules and namespaces from deal.II can be adapted for such purpose. @@ -301,15 +301,15 @@ $\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 +- 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 -of the domain ) then: $\{\textbf{U}_j^n\}_{j \in \mathcal{I}(i)}$ should contain +of the domain) then: $\{\textbf{U}_j^n\}_{j \in \mathcal{I}(i)}$ should contain nine state vector elements (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() in the case -of deal.ii) only collects values for the local cell (just a subset of the +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: @@ -330,55 +330,60 @@ application of this kind of schemes, also called edge-based or graph-based finite element schemes (see for instance @cite Rainald2008 for a historical overview). -

Conservation properties and boundary conditions

+

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) in the right portion of the domain, and -"reflecting boundary" conditions (also called "slip" boundary conditions) -@f{align*} - \mathbf{m}_j \cdot \boldsymbol{\nu}_j \equiv 0 \ \ - \mathbf{x}_j \in \partial\Omega \, . -@f} -in 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 code. In this portion of the documentation we will focus only -the "reflecting" boundary conditions which are somewhat more challenging. - -@note At the time of this writing (early 2020) it is accurate to say that -both analysis and implementation of boundary conditions for hyperbolic systems -of conservation is a widely open issue. Discussions about analysis and/or -implementation of boundary conditions in the academic literature is minimal to -non-existent. - -In this tutorial example we use the so-called "explicit treatment of boundary -conditions": -- Advance in time satisfying no boundary condition at all, -- At the end of the time step enforce boundary conditions strongly in a - post-processing step where we execute the projection +"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 +hyperbolic systems of conservation laws is an open issue. For the case of +variational formulations, stable boundary conditions are those leading to a +well-posed (coercive) bilinear form. But for general hyperbolic +systems of conservation laws (and for the algebraic formulation used in this +tutorial) coercivity has no applicability and/or meaning as a notion of +stability. In this tutorial step we will use preservation of the invariant set +as our main notion of stability which (at the very least) guarantees +well-posedness of the discrete problem. + +For the case of the reflecting boundary conditions we will proceed as follows: +- For every time step advance in time satisfying no boundary condition at all. +- Let $\partial\Omega^r$ be the portion of the boundary where we want to + enforce reflecting boundary conditions. At the end of the time step we enforce + 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) - \boldsymbol{\nu}_i \ \ \text{for all }\mathbf{x}_i \in \partial\Omega + \boldsymbol{\nu}_i \ \ \text{for all }\mathbf{x}_i \in \partial\Omega^r @f} - which removes the normal component of $\mathbf{m}$. Here the definition of - nodal normal $\boldsymbol{\nu}_i$ is very much arbitrary, but it is usually - computed with some form of averaging. - -At this point in time, the average finite element person might find this -approach questionable. Why would you want to do this? No doubt, when solving -parabolic, or elliptic equations, we typically enforce essential -(Dirichlet-like) boundary conditions by making them part of the approximation -space $\mathbb{V}$, and treat natural (e.g. Neumann and Robin) boundary + that removes the normal component of $\mathbf{m}$. Here the definition of + nodal normal $\boldsymbol{\nu}_i$ is very much arbitrary (there is no + unique definition) but it should be consistent upon refinement with the + underlying geometry. + +This is approach is usually called "explicit treatment of boundary conditions". +The well seasoned finite element person might find this approach questionable. +No doubt, when solving parabolic, or elliptic equations, we typically enforce +essential (Dirichlet-like) boundary conditions by making them part of the +approximation space $\mathbb{V}$, and treat natural (e.g. Neumann) boundary conditions as part of the variational formulation. We also know that explicit -treatment of boundary conditions (in the context of parabolic PDE) almost surely -leads to catastrophic consequences. However, in the context of nonlinear -hyperbolic equations: -- The most important reason: it is relatively easy to prove that (for the case -of reflecting boundary conditions) explicit treatment of boundary conditions is -not only conservative but also guarantees preservation of the property -$\mathbf{U}_i \in \mathcal{B}$ for all $i \in \mathcal{V}$. +treatment of boundary conditions (in the context of parabolic PDEs) almost +surely leads to catastrophic consequences. However, in the context of nonlinear +hyperbolic equations we have that: +- It is relatively easy to prove that (for the case of reflecting boundary +conditions) explicit treatment of boundary conditions is not only conservative +but also guarantees preservation of the property $\mathbf{U}_i \in \mathcal{B}$ +for all $i \in \mathcal{V}$ (well-posedness). This is perhaps the most +important reason to use explicit enforcement of boundary conditions. - To the best of our knowledge: we are not aware of any mathematical result proving that it is possible to guarantee the property $\mathbf{U}_i \in \mathcal{B}$ for all $i \in \mathcal{V}$ when using either direct enforcement of @@ -390,29 +395,35 @@ restrictive time step constraints. Dirichlet-like boundary conditions is stable under CFL conditions and does not introduce any loss in accuracy. -If we where to implement reflecting boundary conditions in the entirety of -the boundary, such implementation should be such that exact conservation of -density $\rho$ and mechanical energy $E$ +If $\mathbf{u}_t + \text{div} \, \mathbb{f}(\mathbf{u}) = \boldsymbol{0}$ +represents Euler's equation with reflecting boundary conditions in 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*} -\sum_{i \in \mathcal{V}} m_i \rho_i^{n+1}= \sum_{i \in \mathcal{V}} m_i -\rho_i^{n} \ \ \text{and} \ \ -\sum_{i \in \mathcal{V}} m_i E_i^{n+1}= \sum_{i \in \mathcal{V}} m_i -E_i^{n} +\int_{\Omega} \rho(\mathbf{x},t_2) \, \mathrm{d}\mathbf{x} = +\int_{\Omega} \rho(\mathbf{x},t_1) \, \mathrm{d}\mathbf{x} \ , \ \ +\int_{\Omega} E(\mathbf{x},t_2) \, \mathrm{d}\mathbf{x} = +\int_{\Omega} E(\mathbf{x},t_1) \, \mathrm{d}\mathbf{x} \ , \ \ +\int_{\Omega} \mathbf{m}(\mathbf{x},t_2) \, \mathrm{d}\mathbf{x} = +\int_{\Omega} \mathbf{m}(\mathbf{x},t_1) \, \mathrm{d}\mathbf{x} ++ \int_{t_1}^{t_2} \! \int_{\partial\Omega} p \boldsymbol{\nu} \, +\mathrm{d}\mathbf{s} \mathrm{d}t\, . @f} -is achieved (note that conservation of momentum is not a natural property -Euler's system with reflecting boundary conditions, that's why we did not -include it). Even though we will not use reflecting boundary conditions in the -entirety of the domain, we would like to know that our implementation of +Note that momentum is not a conserved quantity (interaction with walls leads to +momentum gain/loss). Even though we will not use reflecting boundary conditions +in the entirety of the domain, we would like to know that our implementation of boundary conditions is consistent with the conservation properties mentioned above. In order to guarantee such conservation property it is necessary to modify the values of the vectors $\mathbf{c}_{ij}$ as follows @f{align*} -\mathbf{c}_{ij} := \mathbf{c}_{ij} + \int_{\partial \Omega^r} -(\boldsymbol{\nu}_j - \boldsymbol{\nu}(s)) \phi_i \phi_j \, \mathrm{d}s -\ \ \text{whenever both }\mathbf{x}_i\text{ and }\mathbf{x}_j \text{ lie in the +\mathbf{c}_{ij} \dealcoloneq \mathbf{c}_{ij} + \int_{\partial \Omega^r} +(\boldsymbol{\nu}_j - \boldsymbol{\nu}(\mathbf{s})) \phi_i \phi_j \, +\mathrm{d}\mathbf{s} +\ \ \text{whenever both }\mathbf{x}_i\text{ and }\mathbf{x}_j \text{ lie on the boundary} @f} -where $\partial \Omega^r$ is the portion of the boundary where reflecting -boundary conditions are meant to be enforced. This 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. +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/step-69.cc b/examples/step-69/step-69.cc index 81f497c559..6cd9cb4bd8 100644 --- a/examples/step-69/step-69.cc +++ b/examples/step-69/step-69.cc @@ -31,7 +31,7 @@ // 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 +// parallelization, in particular parallel::distributed::Triangulation and // LinearAlgebra::distributed::Vector included through // distributed/tria.h and // lac/la_parallel_vector.h. Instead of a Trilinos, or PETSc @@ -84,15 +84,15 @@ // @sect3{Class template declarations} // -// We begin our actual implementation by declaring all classes with its +// We begin our actual implementation by declaring all classes with their // data structures and methods upfront. In contrast to previous example // steps we use a more fine-grained encapsulation of concepts, data // structures, and parameters into individual classes. A single class thus -// usually centers around either a single data structure (such as the +// 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 // TimeStep class). We typically declare parameter variables -// and scratch data object private and make methods and data structures +// 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 @@ -105,8 +105,7 @@ // 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. - +// mechanism can be found in step-59. namespace Step69 { using namespace dealii; @@ -127,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. As mentioned, we use - // ParameterAcceptor class to automatically populate problem-specific + // the ParameterAcceptor class to automatically populate problem-specific // parameters, such as the geometry information // (length, etc.) or the refinement level // (refinement) from a parameter file. This requires us to @@ -136,7 +135,6 @@ 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 { @@ -171,7 +169,7 @@ 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, + // DoFHandler, SparsityPattern, boundary maps, the lumped mass matrix, // $\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 @@ -184,7 +182,7 @@ namespace Step69 // 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 + // 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 compute and store this geometric @@ -197,7 +195,6 @@ namespace Step69 // from ParameterAcceptor and follow the same idiom of providing a // setup() (and assemble()) method as for the // class Discretization. - template class OfflineData : public ParameterAcceptor { @@ -253,7 +250,7 @@ namespace Step69 // $[\rho,\textbf{m},E]$. // // The purpose of the class members component_names, - // pressure, and speed_of_sound, is evident + // 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})$, @@ -265,13 +262,13 @@ namespace Step69 // 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 + // to the compiler that the programmer thinks 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. + // question improves performance. // - // Finally we note that: + // 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 @@ -283,7 +280,6 @@ namespace Step69 // 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. - template class ProblemDescription { @@ -295,7 +291,7 @@ namespace Step69 using rank1_type = Tensor<1, problem_dimension>; using rank2_type = Tensor<1, problem_dimension, Tensor<1, dim>>; - const static std::array component_names; + const static std::array component_names; static constexpr double gamma = 7. / 5.; @@ -330,16 +326,15 @@ namespace Step69 // It would be desirable to initialize the class in a single shot: // initialize/set the parameters and define the class members that // depend on these default parameters. However, since we do not know the - // actual final values for the parameters, this would be sort of - // meaningless an unsafe in general (we would like to have mechanisms to + // actual values for the parameters, this would be sort of + // meaningless and unsafe in general (we would like to have mechanisms to // check the consistency of the input parameters). Instead of defining // another setup() method to be called (by-hand) after the - // call to ParameterAcceptor::initialize() we provide an + // 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 + // invoking ParameterAcceptor::initialize() for every class // that inherits from ParameterAceptor. - template class InitialValues : public ParameterAcceptor { @@ -375,7 +370,6 @@ namespace Step69 // sets the proper partition and sparsity pattern for the temporary // vector temp and the matrix dij_matrix // respectively. - template class TimeStep : public ParameterAcceptor { @@ -422,12 +416,12 @@ namespace Step69 // 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 quantitiy. In practice // though, the density is a natural candidate, viz. $r := \rho$. - // Schlieren postprocessing is a standard method for enhancing the - // contrast of a visualization inspired by actual experimental X-ray and - // shadowgraphy techniques of visualization. - + // Schlieren + // postprocessing is a standard method for enhancing the contrast of a + // visualization inspired by actual experimental X-ray and shadowgraphy + // techniques of visualization. template class SchlierenPostprocessor : public ParameterAcceptor { @@ -472,7 +466,6 @@ namespace Step69 // separate class TimeLoop that contains an object of every // class and again reads in a number of parameters with the help of the // ParameterAcceptor class. - template class TimeLoop : public ParameterAcceptor { @@ -530,8 +523,7 @@ namespace Step69 // parameter file. Then, in the constructor body every parameter value is // initialized to a sensible default value and registered with the // ParameterAcceptor class with a call to - // ParameterAcceptor::add_parameter. - + // ParameterAcceptor::add_parameter(). template Discretization::Discretization(const MPI_Comm & mpi_communicator, TimerOutput & computing_timer, @@ -584,9 +576,8 @@ namespace Step69 // benchmark configuration consisting of a channel with a disc obstacle, see // @cite GuermondEtAl2018. We construct the geometry by modifying the // mesh generated by GridGenerator::hyper_cube_with_cylindrical_hole(). - // We refer to Step-49, Step-53, and Step-54 for an overview how to + // We refer to step-49, step-53, and step-54 for an overview how to // create advanced meshes. - template void Discretization::setup() { @@ -650,7 +641,6 @@ namespace Step69 // 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 (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f) @@ -679,7 +669,6 @@ namespace Step69 // Not much is done in the constructor of OfflineData other // than initializing the corresponding class members in the // initialization list. - template OfflineData::OfflineData(const MPI_Comm & mpi_communicator, TimerOutput & computing_timer, @@ -695,7 +684,6 @@ namespace Step69 // locally owned and locally relevant DOFs, and initialize a // Utilities::MPI::Partitioner object that is needed for distributed // vectors. - template void OfflineData::setup() { @@ -987,17 +975,19 @@ namespace Step69 // 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. + // well-documented in step-9, step-13 and step-32 among others. // // Finally, assuming that $\mathbf{x}_i$ is a support point at the boundary, - // the (nodal) normals are defined using averaging: + // the (nodal) normals are defined as: // - // $\widehat{\boldsymbol{\nu}}_i := - // \frac{\boldsymbol{\nu}_i}{|\boldsymbol{\nu}_i|}$ where - // $\boldsymbol{\nu}_i := \sum_{T \subset \text{supp}(\phi_i)} + // @f{align*} + // \widehat{\boldsymbol{\nu}}_i := + // \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}) - // \phi_i(\mathbf{x}_{q,F})$ + // \phi_i(\mathbf{x}_{q,F}) + // @f} // // here $T$ denotes elements, // $\text{supp}(\phi_i)$ the support of the shape function $\phi_i$, @@ -1303,7 +1293,7 @@ namespace Step69 on_subranges, 4096); - // Finally, we normalize the vector stored in + // Finally, we normalize the vectors stored in // OfflineData::BoundaryNormalMap. This operation has // not been thread paralellized as it would neither illustrate any // important concept nor lead to any noticeable speed gain. @@ -1315,18 +1305,31 @@ namespace Step69 } } - // In order to implement reflecting boundary conditions - // $\mathbf{m} \cdot \boldsymbol{\nu}_i =0$ (or equivalently $\mathbf{v} - // \cdot \boldsymbol{\nu}_i =0$ ) the vectors $\mathbf{c}_{ij}$ at the + // As commented in the introduction (see section titled "Stable boundary + // conditions and conservation properties") 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) in the right portion of + // the domain, and "reflecting" boundary conditions (also called "slip" + // boundary conditions). With these boundary conditions we should not expect + // any form of conservation to hold. + // + // However, if we were to use reflecting boundary conditions + // $\mathbf{m} \cdot \boldsymbol{\nu}_i =0$ in the entirety of the boundary + // we should preserve density and total mechanical energy. In order to be + // consistent which scenario the vectors $\mathbf{c}_{ij}$ at the // boundary have to be modified as: // - // $\mathbf{c}_{ij} \, +\!\!= \int_{\partial \Omega} - // (\boldsymbol{\nu}_j - \boldsymbol{\nu}(s)) \phi_i \phi_j \, \mathrm{d}s$ + // @f{align*} + // \mathbf{c}_{ij} \, +\!\!= \int_{\partial \Omega} + // (\boldsymbol{\nu}_j - \boldsymbol{\nu}(\mathbf{s})) \phi_i \phi_j \, + // \mathrm{d}\mathbf{s} \ \text{whenever} \ \mathbf{x}_i \text{ and } + // \mathbf{x}_j \text{ lie in the boundary.} + // @f} // - // Otherwise we will not be able to claim conservation. The ideas repeat - // themselves: we use Workstream in order to compute this correction, most - // of the following code is about the definition of the worker - // local_assemble_system. + // 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. { TimerOutput::Scope t(computing_timer, @@ -1434,7 +1437,7 @@ namespace Step69 // In this section we describe the implementation of the class members of // the ProblemDescription class. Most of the code here is // specific for compressible Euler's equations with an ideal gas law. - // If we wanted to re-purpose Step-69 for a different conservation law + // If we wanted to re-purpose step-69 for a different conservation law // (say for: instance the shallow water equation) most of the // implementation of this class would have to change. But most of the other // classes, however, (in particular those defining loop structures) would @@ -2172,42 +2175,18 @@ namespace Step69 // Step 4: Fix up boundary states. // As a last step in the Forward Euler method, we have to fix up all - // boundary states. This approach is an example of the explicit - // treatment of boundary conditions strategy: + // boundary states. As discussed in the intro we // - advance in time satisfying no boundary condition at all, // - at the end of the time step enforce boundary conditions strongly // in a post-processing step. // - // When solving parabolic, or elliptic equations, we typically enforce - // essential boundary conditions by making them part of the - // approximation space, and treat natural boundary conditions as part - // of the variational formulation. We also know that explicit treatment - // of boundary conditions (in the context of parabolic PDE) almost - // surely leads to catastrophic consequences. However, in the context - // of nonlinear hyperbolic equations there is enough numerical evidence - // suggesting that explicit treatment of essential boundary conditions - // is stable and does not introduce any loss in accuracy and - // convergence rates. In addition, it is relatively straightforward to - // prove that (for the case of reflecting boundary conditions) explicit - // treatment of boundary conditions is not only conservative but also - // guarantees preservation of the invariant set. We are not aware of - // any theoretical result showing that it is possible to provide such - // invariant-set guarantees when using either direct enforcement of - // boundary conditions into the approximation space, or weak - // enforcement using the Nitsche penalty method (which is for example - // widely used in discontinuous Galerkin schemes). - // // Here the worker on_subranges executes the correction // // \f[ // \mathbf{m}_i := \mathbf{m}_i - (\boldsymbol{\nu}_i \cdot \mathbf{m}_i) // \boldsymbol{\nu}_i, // \f] - // which removes the normal component of $\mathbf{m}$. We note that - // conservation is not just a consequence of this correction but also a - // consequence of modification of the $\mathbf{c}_{ij}$ coefficients at - // the boundary that we employed in - // OfflineData::assemble(). + // which removes the normal component of $\mathbf{m}$. { TimerOutput::Scope time(computing_timer,