research computations you might want to consider exploring a corresponding
high-performance implementation of a second-order accurate scheme that uses
<i>convex limiting</i> techniques, and strong stability-preserving (SSP)
-time integration @cite GuermondEtAl2018. The repository can be found at
-TODO.
-
-
-<!-- ####################################################################### -->
-<!-- ####################################################################### -->
-<!-- ####################################################################### -->
+time integration, see @cite GuermondEtAl2018.
<h1>Introduction</h1>
the second-order scheme @cite GuermondEtAl2018 maintained TODO.
-<!-- ####################################################################### -->
-<!-- ####################################################################### -->
-<!-- ####################################################################### -->
-
-
<a name="eulerequations"></a>
<h3>Euler's equations of gas dynamics</h3>
specific heats</a>, and $\|\,.\|$ denotes the Euclidian norm.
-<!-- ####################################################################### -->
-<!-- ####################################################################### -->
-
-
<h4>Solution theory</h4>
Hyperbolic conservation laws, such as
$\mathbf{u}(\mathbf{x},t)$ remains in $\mathcal{B}$.
-<!-- ####################################################################### -->
-<!-- ####################################################################### -->
-
-
<h4>Variational versus collocation-type discretizations</h3>
Following Step-9, Step-12, and Step-33, at this point it might look tempting
schemes in Hilbert spaces and that a large number of classes, modules and
namespaces from deal.ii can be adapted for such purpose.
-<!-- ####################################################################### -->
-<!-- ####################################################################### -->
-<!-- ####################################################################### -->
<h3>Description of the scheme </h3>
\lambda_{\text{max}} (\mathbf{U}_j^{n}, \mathbf{U}_i^{n},
\textbf{n}_{ji}) \} \|\mathbf{c}_{ij}\|_{\ell^2} $
+<<<<<<< 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
@f{align*}
&\textbf{For } i \in \mathcal{V} \\
"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
-<!-- ####################################################################### -->
-<!-- ####################################################################### -->
-<!-- ####################################################################### -->
<h3>Implementation of the scheme </h3>
-
-
-
-
-
-
-
-
*/
// @sect3{Include files}
-// The set of include files is quite standard. The most intriguing part
+// 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 <code>la_parallel_vector.h</code> and the regular dealii
// \frac{\boldsymbol{\nu}_i}{|\boldsymbol{\nu}_i|}$ where
// $\boldsymbol{\nu}_i = \sum_{F \subset \text{supp}(\phi_i)}
// \sum_{\mathbf{x}_{q,F}} \nu(\mathbf{x}_{q,F})
- // \phi_i(\mathbf{x}_{q,F})$, here: $F \subset \partial \Omega$ denotes
- // faces of elements at the boundary of the domain, and $\mathbf{x}_{q,F}$
+ // \phi_i(\mathbf{x}_{q,F})$, here: $F \subset \partial \Omega$ denotes
+ // faces of elements at the boundary of the domain, and $\mathbf{x}_{q,F}$
// are quadrature points on such face.
- // Other more sophisticated definitions for $\nu_i$ are
+ // Other more sophisticated definitions for $\nu_i$ are
// possible but none of them have much influence in theory or practice.
// We remind the reader that <code>CopyData</code> includes the class member
// <code>local_boundary_normal_map</code> in order to store these local
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" and pass
+ /* 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" and pass
it to the copy data routine. */
Tensor<1, dim> normal;
if (id == Boundary::slip)
for (const auto &it : local_boundary_normal_map)
{
- auto &[normal, id, position] = boundary_normal_map[it.first];
- auto &[new_normal, new_id, new_position] = it.second;
+ auto &normal = std::get<0>(boundary_normal_map[it.first]);
+ auto &id = std::get<1>(boundary_normal_map[it.first]);
+ auto &position = std::get<2>(boundary_normal_map[it.first]);
+
+ const auto &new_normal = std::get<0>(it.second);
+ const auto &new_id = std::get<1>(it.second);
+ const auto &new_position = std::get<2>(it.second);
normal += new_normal;
id = std::max(id, new_id);
// contains a just copy of the matrix <code>cij_matrix</code>.
// That's not what we really
// want: we have to normalize its entries. In addition, we have not even
- // touched the entries of the matrix <code>norm_matrix</code> yet, and the
+ // touched the entries of the matrix <code>norm_matrix</code> yet, and the
// vectors stored in the map
// <code>OfflineData<dim>::BoundaryNormalMap</code> are not normalized.
//
// In principle, this is just offline data, it doesn't make much sense
// to over-optimize their computation, since their cost will get amortized
- // over the many time steps that we are going to use. However,
+ // over the many time steps that we are going to use. However,
// computing/storing the entries of the matrix
// <code>norm_matrix</code> and the normalization of <code>nij_matrix</code>
- // are perfect to illustrate thread-parallel node-loops:
- // - We want to visit every node $i$ in the mesh/sparsity graph,
+ // are perfect to illustrate thread-parallel node-loops:
+ // - We want to visit every node $i$ in the mesh/sparsity graph,
// - and for every such node we want to visit to every $j$ such that
// $\mathbf{c}_{ij} \not \equiv 0$.
//
- // From an algebraic point of view, this is equivalent to: visiting
+ // From an algebraic point of view, this is equivalent to: visiting
// every row in the matrix (equivalently sparsity
// pattern) and for each one of these rows execute a loop on the columns.
// Node-loops is a core theme of this tutorial step (see the pseudo-code
- // in the introduction) that will repeat over and over again. That's why
+ // in the introduction) that will repeat over and over again. That's why
// this is the right time to introduce them.
//
// We have the thread paralellization capability
// parallel::apply_to_subranges that is somehow more general than the
- // WorkStream framework. In particular, it can be used for our
+ // WorkStream framework. In particular, it can be used for our
// node-loops.
// This functionality requires four input arguments:
// - A begin iterator: <code>indices.begin()</code>
// element schemes when they are properly implemented.
//
// Finally, we normalize the vector stored in
- // <code>OfflineData<dim>::BoundaryNormalMap</code>. This operation has
- // not been thread paralellized as it would not illustrate any important
+ // <code>OfflineData<dim>::BoundaryNormalMap</code>. This operation has
+ // not been thread paralellized as it would not illustrate any important
// concept.
{
/* This is not thread parallelized, too bad! */
for (auto &it : boundary_normal_map)
{
- auto &[normal, id, _] = it.second;
+ auto &normal = std::get<0>(it.second);
normal /= (normal.norm() + std::numeric_limits<double>::epsilon());
}
}
// 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
+ // \cdot \boldsymbol{\nu}_i =0$ ) 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_j \, \mathrm{d}s$
//
// 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
+ // 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>.
{
// Now we define the implementation of <code>momentum</code>,
// <code>internal_energy</code>, <code>pressure</code>,
// <code>speed_of_sound</code>, and <code>f</code> (the flux of the system).
- // The functionality of each one of these functions is self-explanatory from
+ // The functionality of each one of these functions is self-explanatory from
// their names.
template <int dim>
return result;
}
- // The following function, <code>riemann_data_from_state</code>, takes the
- // full state $\mathbf{u} = [\rho,\mathbf{m},E]^\top$ defines a new
+ // The following function, <code>riemann_data_from_state</code>, takes the
+ // full state $\mathbf{u} = [\rho,\mathbf{m},E]^\top$ defines a new
// "projected state" defined as
//
// $\widetilde{\mathbf{u}} = [\rho,
// \mathbf{m} - (\mathbf{m}\cdot \mathbf{n}_{ij})\mathbf{n}_{ij},
// E - \tfrac{(\mathbf{m}\cdot \mathbf{n}_{ij})^2}{2\rho} ]^\top$
//
- // Projected states appear naturally when attempting to compute a maximum
+ // Projected states appear naturally when attempting to compute a maximum
// wavespeed appearing in Riemann problems.
namespace
return std::max(std::abs(u_i), std::abs(u_j)) + 5. * std::max(a_i, a_j);
}
- } /* End of namespace dedicated to the computation of the maximum wavespeed */
+ } // namespace
// Placeholder here.
if (i >= n_locally_owned)
continue;
- const auto &[normal, id, position] = it->second;
+ const auto &normal = std::get<0>(it->second);
+ const auto &id = std::get<1>(it->second);
+ const auto &position = std::get<2>(it->second);
/* Skip constrained degrees of freedom */
if (++sparsity.begin(i) == sparsity.end(i))
/* On boundary 2 enforce initial conditions: */
- if (id == Boundary::dirichlet)
+ else if (id == Boundary::dirichlet)
{
U_i = initial_values->initial_state(position, t + tau_max);
}
const auto bnm_it = boundary_normal_map.find(i);
if (bnm_it != boundary_normal_map.end())
{
- const auto [normal, id, _] = bnm_it->second;
+ const auto &normal = std::get<0>(bnm_it->second);
+ const auto &id = std::get<1>(bnm_it->second);
+
if (id == Boundary::slip)
- {
- r_i -= 1. * (r_i * normal) * normal;
- }
+ r_i -= 1. * (r_i * normal) * normal;
else
- {
- r_i = 0.;
- }
+ r_i = 0.;
}
const double m_i = lumped_mass_matrix.diag_element(i);
projects / dealii.git / commitdiff