of finite elements; a graph-viscosity stabilization based on a
<i>guaranteed</i> 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
@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.
<a name="eulerequations"></a>
@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*}
@f}
Such solutions, which are understood as the solution recovered in the
zero-viscosity limit, are often refered to as <i>viscosity solutions</i>.
-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
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) |\!|\!|
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,
@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.
- $\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<dim, spacedim>.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:
<i>graph-based</i> finite element schemes (see for instance
@cite Rainald2008 for a historical overview).
-<h3>Conservation properties and boundary conditions</h3>
+<h3>Stable boundary conditions and conservation properties.</h3>
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
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.
+
+
// 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
// <code>distributed/tria.h</code> and
// <code>lac/la_parallel_vector.h</code>. Instead of a Trilinos, or PETSc
// @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 <code>Discretization</code> class, or a single
// method (such as the <code>step()</code> function of the
// <code>TimeStep</code> 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
// 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;
// The class <code>Discretization</code> 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
// (<code>length</code>, etc.) or the refinement level
// (<code>refinement</code>) from a parameter file. This requires us to
// constructor, and defer the creation of the mesh to the
// <code>setup()</code> method that can be called once all parameters are
// read-in via ParameterAcceptor::initialize().
-
template <int dim>
class Discretization : public ParameterAcceptor
{
//
// The class <code>OfflineData</code> contains pretty much all components
// of the discretization that do not evolve in time, in particular, the
- // DoFHandler, SparsityPattern, boundary maps, the lumped mass,
+ // DoFHandler, SparsityPattern, boundary maps, the lumped mass matrix,
// $\mathbf{c}_{ij}$ and $\mathbf{n}_{ij}$ matrices. Here, the term
// <i>offline</i> refers to the fact that all the class
// members of <code>OfflineData</code> have well-defined values
// either.
//
// We also compute and store a <code>boundary_normal_map</code> 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
// from ParameterAcceptor and follow the same idiom of providing a
// <code>setup()</code> (and <code>assemble()</code>) method as for the
// class Discretization.
-
template <int dim>
class OfflineData : public ParameterAcceptor
{
// $[\rho,\textbf{m},E]$.
//
// The purpose of the class members <code>component_names</code>,
- // <code>pressure</code>, and <code>speed_of_sound</code>, is evident
+ // <code>pressure</code>, and <code>speed_of_sound</code> is evident
// from their names. We also provide a function
// <code>compute_lambda_max</code>, that computes the wave speed estimate
// mentioned above, $\lambda_{max}(\mathbf{U},\mathbf{V},\mathbf{n})$,
// body is put in place for every invocation of the function, and no call
// (and code indirection) is generated. This is stronger than the
// <code>inline</code> 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. <code>DEAL_II_ALWAYS_INLINE</code> 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
// 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 <int dim>
class ProblemDescription
{
using rank1_type = Tensor<1, problem_dimension>;
using rank2_type = Tensor<1, problem_dimension, Tensor<1, dim>>;
- const static std::array<std::string, dim + 2> component_names;
+ const static std::array<std::string, problem_dimension> component_names;
static constexpr double gamma = 7. / 5.;
// 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 <code>setup()</code> method to be called (by-hand) after the
- // call to <code> ParameterAcceptor::initialize() </code> we provide an
+ // call to ParameterAcceptor::initialize() we provide an
// "implementation" for the class member
// <code>parse_parameters_call_back</code> which is automatically called when
- // invoking <code> ParameterAcceptor::initialize() </code> for every class
+ // invoking ParameterAcceptor::initialize() for every class
// that inherits from ParameterAceptor.
-
template <int dim>
class InitialValues : public ParameterAcceptor
{
// sets the proper partition and sparsity pattern for the temporary
// vector <code>temp</code> and the matrix <code>dij_matrix</code>
// respectively.
-
template <int dim>
class TimeStep : public ParameterAcceptor
{
// <code>schlieren</code>, 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.
-
+ // <a href="https://en.wikipedia.org/wiki/Schlieren">Schlieren</a>
+ // postprocessing is a standard method for enhancing the contrast of a
+ // visualization inspired by actual experimental X-ray and shadowgraphy
+ // techniques of visualization.
template <int dim>
class SchlierenPostprocessor : public ParameterAcceptor
{
// separate class <code>TimeLoop</code> that contains an object of every
// class and again reads in a number of parameters with the help of the
// ParameterAcceptor class.
-
template <int dim>
class TimeLoop : public ParameterAcceptor
{
// 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 <int dim>
Discretization<dim>::Discretization(const MPI_Comm & mpi_communicator,
TimerOutput & computing_timer,
// 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 <int dim>
void Discretization<dim>::setup()
{
// <code>Boundary::dirichlet</code> on the left and
// <code>Boundary::slip</code> 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)
// Not much is done in the constructor of <code>OfflineData</code> other
// than initializing the corresponding class members in the
// initialization list.
-
template <int dim>
OfflineData<dim>::OfflineData(const MPI_Comm & mpi_communicator,
TimerOutput & computing_timer,
// locally owned and locally relevant DOFs, and initialize a
// Utilities::MPI::Partitioner object that is needed for distributed
// vectors.
-
template <int dim>
void OfflineData<dim>::setup()
{
// <code>local_assemble_system</code> and the copy data routine
// <code>copy_local_to_global</code>. 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$,
on_subranges,
4096);
- // Finally, we normalize the vector stored in
+ // Finally, we normalize the vectors stored in
// <code>OfflineData<dim>::BoundaryNormalMap</code>. This operation has
// not been thread paralellized as it would neither illustrate any
// important concept nor lead to any noticeable speed gain.
}
}
- // 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
- // <code>local_assemble_system</code>.
+ // 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 <code>local_assemble_system</code>.
{
TimerOutput::Scope t(computing_timer,
// In this section we describe the implementation of the class members of
// the <code>ProblemDescription</code> 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
// <b>Step 4</b>: 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 <i>explicit
- // treatment of boundary conditions</i> 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 <code>on_subranges</code> 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
- // <code>OfflineData<dim>::assemble()</code>.
+ // which removes the normal component of $\mathbf{m}$.
{
TimerOutput::Scope time(computing_timer,
projects / dealii.git / commitdiff