1506 * cell->get_dof_indices(dof_indices);
1507 * std::transform(dof_indices.begin(),
1508 * dof_indices.end(),
1509 * dof_indices.begin(),
1511 * return partitioner->global_to_local(index);
1516 *
for (
const auto dof : dof_indices)
1517 * dsp.add_entries(dof, dof_indices.
begin(), dof_indices.
end());
1520 * sparsity_pattern.copy_from(dsp);
1522 * lumped_mass_matrix.reinit(sparsity_pattern);
1523 * norm_matrix.reinit(sparsity_pattern);
1524 *
for (
auto &matrix : cij_matrix)
1526 *
for (
auto &matrix : nij_matrix)
1534 * Next, we have to
assemble various matrices. We define a number of
1535 * helper
functions and
data structures in an anonymous
namespace.
1545 * <code>CopyData</code>
class that will be used to
assemble the
1546 * offline
data matrices
using WorkStream. It acts as a container: it
1547 * is just a
struct where
WorkStream stores the local cell
1548 * contributions. Note that it also contains a
class member
1549 * <code>local_boundary_normal_map</code> used to store the local
1550 * contributions required to compute the normals at the boundary.
1556 *
template <
int dim>
1559 *
bool is_artificial;
1560 * std::vector<types::global_dof_index> local_dof_indices;
1561 *
typename OfflineData<dim>::BoundaryNormalMap local_boundary_normal_map;
1563 * std::array<FullMatrix<double>, dim> cell_cij_matrix;
1568 * Next we introduce a number of helper
functions that are all
1569 * concerned about reading and writing
matrix and vector entries. They
1570 * are mainly motivated by providing slightly more efficient code and
1571 * <a href=
"https://en.wikipedia.org/wiki/Syntactic_sugar"> syntactic
1572 * sugar</a>
for otherwise somewhat tedious code.
1576 * The
first function we introduce, <code>get_entry()</code>, will be
1577 * used to read the
value stored at the entry pointed by a
1579 * function works around a small deficiency in the
SparseMatrix
1581 * operations of the sparse
matrix stored in CRS format. As such the
1582 * iterator already knows the global
index of the corresponding
matrix
1584 * to the lack of an
interface in the
SparseMatrix for accessing the
1586 * have to create a temporary
SparseMatrix iterator. We simply hide
1587 *
this in the <code>get_entry()</code> function.
1593 *
template <
typename IteratorType>
1598 * &matrix, it->global_index());
1599 *
return matrix_iterator->value();
1604 * The <code>set_entry()</code> helper is the inverse operation of
1605 * <code>get_value()</code>: Given an iterator and a
value, it sets the
1606 * entry pointed to by the iterator in the
matrix.
1612 *
template <
typename IteratorType>
1615 *
const IteratorType &it,
1619 * it->global_index());
1620 * matrix_iterator->value() =
value;
1625 * <code>gather_get_entry()</code>: we note that @f$\mathbf{c}_{ij} \in
1626 * \mathbb{R}^
d@f$. If @f$d=2@f$ then @f$\mathbf{c}_{ij} =
1627 * [\mathbf{c}_{ij}^1,\mathbf{c}_{ij}^2]^\top@f$. Which basically implies
1628 * that we need one
matrix per space dimension to store the
1629 * @f$\mathbf{c}_{ij}@f$ vectors. Similar observation follows
for the
1630 *
matrix @f$\mathbf{n}_{ij}@f$. The purpose of
1631 * <code>gather_get_entry()</code> is to retrieve those entries and store
1638 *
template <std::
size_t k,
typename IteratorType>
1641 *
const IteratorType it)
1644 *
for (
unsigned int j = 0; j < k; ++j)
1645 * result[j] = get_entry(c_ij[j], it);
1652 * signature, having three input arguments, will be used to retrieve
1653 * the individual components <code>(i,
l)</code> of a
matrix. The
1654 * functionality of <code>gather_get_entry()</code> and
1655 * <code>
gather()</code> is very much the same, but their context is
1656 * different: the function <code>
gather()</code> does not rely on an
1657 * iterator (that actually knows the
value pointed to) but rather on the
1658 * indices <code>(i,
l)</code> of the entry in order to retrieve its
1659 * actual
value. We should expect <code>
gather()</code> to be slightly
1660 * more expensive than <code>gather_get_entry()</code>. The use of
1661 * <code>
gather()</code> will be limited to the task of computing the
1662 * algebraic viscosity @f$d_{ij}@f$ in the particular case that when
1663 * both @f$i@f$ and @f$j@f$ lie at the boundary.
1667 * @note The reader should be aware that accessing an arbitrary
1668 * <code>(i,
l)</code> entry of a
matrix (say for instance Trilinos or
PETSc
1669 * matrices) is in
general unacceptably expensive. Here is where we might
1670 * want to keep an eye on complexity: we want this operation to have
1671 *
constant complexity, which is the case of the current implementation
1672 * using deal.II matrices.
1678 * template <
std::size_t k>
1681 * const unsigned
int i,
1682 * const unsigned
int j)
1685 *
for (
unsigned int l = 0;
l < k; ++
l)
1686 * result[l] = n_ij[l](i, j);
1693 * signature having two input arguments will be used to
gather the
1694 * state at a node <code>i</code> and return it as a
1695 * <code>
Tensor<1,n_solution_variables></code> for our convenience.
1701 * template <
std::size_t k>
1704 * const unsigned
int i)
1707 *
for (
unsigned int j = 0; j < k; ++j)
1708 * result[j] = U[j].local_element(i);
1714 * <code>
scatter()</code>:
this function has three input arguments, the
1715 *
first one is meant to be a
"global object" (say a locally owned or
1716 * locally relevant vector), the
second argument which could be a
1718 * which represents a index of the global
object. This function will be
1719 * primarily used to write the updated nodal values, stored as
1726 * template <std::size_t k, int k2>
1730 *
const unsigned int i)
1732 *
static_assert(k == k2,
1733 *
"The dimensions of the input arguments must agree");
1734 *
for (
unsigned int j = 0; j < k; ++j)
1735 * U[j].local_element(i) = tensor[j];
1741 * We are now in a position to
assemble all matrices stored in
1742 * <code>OfflineData</code>: the lumped mass entries @f$m_i@f$, the
1743 * vector-valued matrices @f$\mathbf{c}_{ij}@f$ and @f$\mathbf{n}_{ij} =
1744 * \frac{\mathbf{c}_{ij}}{|\mathbf{c}_{ij}|}@f$, and the boundary normals
1745 * @f$\boldsymbol{\nu}_i@f$.
1749 * In order to exploit thread parallelization we use the
WorkStream approach
1750 * detailed in the @ref threads
"Parallel computing with multiple processors"
1751 * accessing shared memory. As customary
this requires
1753 * - Scratch
data (i.e. input info required to carry out computations): in
1754 * this case it is <code>scratch_data</code>.
1755 * - The worker: in our case this is the <code>local_assemble_system()</code>
1757 * actually computes the local (i.
e. current cell) contributions from the
1759 * - A
copy data: a struct that contains all the local assembly
1760 * contributions, in this case <code>CopyData<dim>()</code>.
1761 * - A
copy data routine: in this case it is
1762 * <code>copy_local_to_global()</code> in charge of actually copying these
1763 * local contributions into the global objects (matrices and/or vectors)
1767 * Most of the following lines are spent in the definition of the worker
1768 * <code>local_assemble_system()</code> and the
copy data routine
1769 * <code>copy_local_to_global()</code>. There is not much to say about the
1770 *
WorkStream framework since the vast majority of ideas are reasonably
1771 * well-documented in @ref step_9
"step-9", @ref step_13
"step-13" and @ref step_32
"step-32" among others.
1775 * Finally, assuming that @f$\mathbf{x}_i@f$ is a support
point at the boundary,
1776 * the (nodal) normals are defined as:
1781 * \widehat{\boldsymbol{\nu}}_i \dealcoloneq
1782 * \frac{\int_{\partial\Omega} \phi_i \widehat{\boldsymbol{\nu}} \,
1783 * \, \mathrm{
d}\mathbf{s}}{\big|\int_{\partial\Omega} \phi_i
1784 * \widehat{\boldsymbol{\nu}} \, \mathrm{
d}\mathbf{s}\big|}
1789 * We will compute the numerator of
this expression
first and store it in
1790 * <code>OfflineData<dim>::BoundaryNormalMap</code>. We will normalize these
1791 * vectors in a posterior
loop.
1797 *
template <
int dim>
1798 *
void OfflineData<dim>::assemble()
1800 * lumped_mass_matrix = 0.;
1802 *
for (
auto &matrix : cij_matrix)
1804 *
for (
auto &matrix : nij_matrix)
1807 *
unsigned int dofs_per_cell =
1808 * discretization->finite_element.n_dofs_per_cell();
1809 *
unsigned int n_q_points = discretization->quadrature.size();
1813 * What follows is the initialization of the scratch
data required by
1821 * discretization->mapping,
1822 * discretization->finite_element,
1823 * discretization->quadrature,
1826 * discretization->face_quadrature,
1832 *
"offline_data - assemble lumped mass matrix, and c_ij");
1834 *
const auto local_assemble_system =
1837 * CopyData<dim> &
copy) {
1838 *
copy.is_artificial = cell->is_artificial();
1839 *
if (
copy.is_artificial)
1842 *
copy.local_boundary_normal_map.clear();
1843 *
copy.cell_lumped_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
1844 *
for (
auto &matrix :
copy.cell_cij_matrix)
1847 *
const auto &fe_values = scratch.reinit(cell);
1849 *
copy.local_dof_indices.resize(dofs_per_cell);
1850 * cell->get_dof_indices(
copy.local_dof_indices);
1852 * std::transform(
copy.local_dof_indices.begin(),
1853 *
copy.local_dof_indices.end(),
1854 *
copy.local_dof_indices.begin(),
1856 * return partitioner->global_to_local(index);
1861 * We compute the local contributions
for the lumped mass
matrix
1862 * entries @f$m_i@f$ and and vectors @f$c_{ij}@f$ in the usual fashion:
1865 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1867 *
const auto JxW = fe_values.JxW(q_point);
1869 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1871 *
const auto value_JxW =
1872 * fe_values.shape_value(j, q_point) * JxW;
1873 *
const auto grad_JxW = fe_values.shape_grad(j, q_point) * JxW;
1875 *
copy.cell_lumped_mass_matrix(j, j) += value_JxW;
1877 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1879 *
const auto value = fe_values.shape_value(i, q_point);
1880 *
for (
unsigned int d = 0;
d < dim; ++
d)
1881 *
copy.cell_cij_matrix[d](i, j) +=
value * grad_JxW[
d];
1889 * Now we have to compute the boundary normals. Note that the
1890 * following
loop does not
do much unless the element has faces on
1891 * the boundary of the domain.
1894 *
for (
const auto f : cell->face_indices())
1896 * const auto face = cell->face(f);
1897 *
const auto id = face->boundary_id();
1899 *
if (!face->at_boundary())
1902 *
const auto &fe_face_values = scratch.reinit(cell, f);
1904 *
const unsigned int n_face_q_points =
1905 * fe_face_values.get_quadrature().size();
1907 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1909 *
if (!discretization->finite_element.has_support_on_face(j, f))
1914 * Note that
"normal" will only represent the contributions
1915 * from one of the faces in the support of the shape
1916 * function phi_j. So we cannot normalize
this local
1917 * contribution right here, we have to take it
"as is",
1918 * store it and pass it to the
copy data routine. The
1919 * proper normalization
requires an additional
loop on
1920 * nodes. This is done in the
copy function below.
1924 *
if (
id == Boundaries::free_slip)
1926 *
for (
unsigned int q = 0; q < n_face_q_points; ++q)
1927 * normal += fe_face_values.normal_vector(q) *
1928 * fe_face_values.shape_value(j, q);
1931 *
const auto index =
copy.local_dof_indices[j];
1935 * if (cell->vertex_dof_index(v, 0) ==
1936 * partitioner->local_to_global(
index))
1938 * position = cell->vertex(v);
1942 *
const auto old_id =
1943 * std::get<1>(
copy.local_boundary_normal_map[index]);
1944 *
copy.local_boundary_normal_map[
index] =
1945 * std::make_tuple(normal,
std::max(old_id,
id), position);
1952 * Last, we provide a copy_local_to_global function as required
for
1956 *
const auto copy_local_to_global = [&](
const CopyData<dim> &
copy) {
1957 *
if (
copy.is_artificial)
1960 *
for (
const auto &it :
copy.local_boundary_normal_map)
1962 *
std::get<0>(boundary_normal_map[it.
first]) +=
1964 * std::get<1>(boundary_normal_map[it.first]) =
1965 *
std::max(std::get<1>(boundary_normal_map[it.first]),
1966 * std::get<1>(it.second));
1967 * std::get<2>(boundary_normal_map[it.first]) = std::get<2>(it.second);
1970 * lumped_mass_matrix.add(
copy.local_dof_indices,
1971 *
copy.cell_lumped_mass_matrix);
1973 *
for (
int k = 0; k < dim; ++k)
1975 * cij_matrix[k].add(
copy.local_dof_indices,
copy.cell_cij_matrix[k]);
1976 * nij_matrix[k].add(
copy.local_dof_indices,
copy.cell_cij_matrix[k]);
1981 * dof_handler.end(),
1982 * local_assemble_system,
1983 * copy_local_to_global,
1990 * At
this point in time we are done with the computation of @f$m_i@f$ and
1991 * @f$\mathbf{c}_{ij}@f$, but so far the matrix <code>nij_matrix</code>
1992 * contains just a
copy of the matrix <code>cij_matrix</code>.
1993 * That
's not what we really want: we have to normalize its entries. In
1994 * addition, we have not filled the entries of the matrix
1995 * <code>norm_matrix</code> and the vectors stored in the map
1996 * <code>OfflineData<dim>::BoundaryNormalMap</code> are not normalized.
2000 * In principle, this is just offline data, it doesn't make much sense
2001 * to over-optimize their computation, since their cost will get amortized
2002 * over the many time steps that we are going to use. However,
2003 * computing/storing the entries of the
matrix
2004 * <code>norm_matrix</code> and the normalization of <code>nij_matrix</code>
2005 * are perfect to illustrate thread-
parallel node-loops:
2006 * - we want to visit every node @f$i@f$ in the mesh/sparsity graph,
2007 * - and
for every such node we want to visit to every @f$j@f$ such that
2008 * @f$\mathbf{c}_{ij} \not \equiv 0@f$.
2012 * From an algebraic
point of view,
this is equivalent to: visiting
2013 * every row in the
matrix and
for each one of these rows execute a
loop on
2014 * the columns. Node-loops is a core theme of
this tutorial step (see
2015 * the pseudo-code in the introduction) that will repeat over and over
2016 * again. That
's why this is the right time to introduce them.
2020 * We have the thread parallelization capability
2021 * parallel::apply_to_subranges() that is somehow more general than the
2022 * WorkStream framework. In particular, parallel::apply_to_subranges() can
2023 * be used for our node-loops. This functionality requires four input
2024 * arguments which we explain in detail (for the specific case of our
2025 * thread-parallel node loops):
2026 * - The iterator <code>indices.begin()</code> points to a row index.
2027 * - The iterator <code>indices.end()</code> points to a numerically higher
2029 * - The function <code>on_subranges(index_begin, index_end)</code>
2030 * (where <code>index_begin</code> and <code>index_end</code>
2031 * define a sub-range within the range spanned by
2032 * the begin and end iterators defined in the two previous bullets)
2033 * applies an operation to every iterator in such subrange. We may as
2034 * well call <code>on_subranges</code> the "worker".
2035 * - Grainsize: minimum number of iterators (in this case representing
2036 * rows) processed by each thread. We decided for a minimum of 4096
2041 * A minor caveat here is that the iterators <code>indices.begin()</code>
2042 * and <code>indices.end()</code> supplied to
2043 * parallel::apply_to_subranges() have to be random access iterators:
2044 * internally, parallel::apply_to_subranges() will break the range
2045 * defined by the <code>indices.begin()</code> and
2046 * <code>indices.end()</code> iterators into subranges (we want to be
2047 * able to read any entry in those subranges with constant complexity).
2048 * In order to provide such iterators we resort to
2049 * std_cxx20::ranges::iota_view.
2053 * The bulk of the following piece of code is spent defining
2054 * the "worker" <code>on_subranges</code>: i.e. the operation applied at
2055 * each row of the sub-range. Given a fixed <code>row_index</code>
2056 * we want to visit every column/entry in such row. In order to execute
2057 * such columns-loops we use
2058 * <a href="http://www.cplusplus.com/reference/algorithm/for_each/">
2060 * from the standard library, where:
2061 * - <code>sparsity_pattern.begin(row_index)</code>
2062 * gives us an iterator starting at the first column of the row,
2063 * - <code>sparsity_pattern.end(row_index)</code> is an iterator pointing
2064 * at the last column of the row,
2065 * - the last argument required by `std::for_each` is the operation
2066 * applied at each nonzero entry (a lambda expression in this case)
2071 * We note that, parallel::apply_to_subranges() will operate on
2072 * disjoint sets of rows (the subranges) and our goal is to write into
2073 * these rows. Because of the simple nature of the operations we want
2074 * to carry out (computation and storage of normals, and normalization
2075 * of the @f$\mathbf{c}_{ij}@f$ of entries) threads cannot conflict
2076 * attempting to write the same entry (we do not need a scheduler).
2080 * TimerOutput::Scope scope(computing_timer,
2081 * "offline_data - compute |c_ij|, and n_ij");
2083 * const std_cxx20::ranges::iota_view<unsigned int, unsigned int> indices(
2084 * 0, n_locally_relevant);
2086 * const auto on_subranges =
2087 * [&](const auto index_begin, const auto index_end) {
2088 * for (const auto row_index :
2089 * std_cxx20::ranges::iota_view<unsigned int, unsigned int>(
2090 * *index_begin, *index_end))
2094 * First column-loop: we compute and store the entries of the
2095 * matrix norm_matrix and write normalized entries into the
2096 * matrix nij_matrix:
2099 * std::for_each(sparsity_pattern.begin(row_index),
2100 * sparsity_pattern.end(row_index),
2101 * [&](const SparsityPatternIterators::Accessor &jt) {
2103 * gather_get_entry(cij_matrix, &jt);
2104 * const double norm = c_ij.norm();
2106 * set_entry(norm_matrix, &jt, norm);
2107 * for (unsigned int j = 0; j < dim; ++j)
2108 * set_entry(nij_matrix[j], &jt, c_ij[j] / norm);
2113 * parallel::apply_to_subranges(indices.begin(),
2120 * Finally, we normalize the vectors stored in
2121 * <code>OfflineData<dim>::BoundaryNormalMap</code>. This operation has
2122 * not been thread parallelized as it would neither illustrate any
2123 * important concept nor lead to any noticeable speed gain.
2126 * for (auto &it : boundary_normal_map)
2128 * auto &normal = std::get<0>(it.second);
2129 * normal /= (normal.norm() + std::numeric_limits<double>::epsilon());
2136 * At this point we are very much done with anything related to offline data.
2141 * <a name="step_69-EquationofstateandapproximateRiemannsolver"></a>
2142 * <h4>Equation of state and approximate Riemann solver</h4>
2146 * In this section we describe the implementation of the class members of
2147 * the <code>ProblemDescription</code> class. Most of the code here is
2148 * specific to the compressible Euler's equations with an ideal gas law.
2149 * If we wanted to re-purpose @ref step_69
"step-69" for a different conservation law
2150 * (say
for: instance the shallow water equation) most of the
2151 * implementation of
this class would have to change. But most of the
2152 * other classes (in particular those defining loop structures) would
2157 * We start by implementing a number of small member
functions for
2158 * computing <code>momentum</code>, <code>internal_energy</code>,
2159 * <code>pressure</code>, <code>speed_of_sound</code>, and the flux
2160 * <code>f</code> of the system. The functionality of each one of these
2161 *
functions is self-explanatory from their names.
2167 *
template <
int dim>
2169 * ProblemDescription<dim>::momentum(
const state_type &U)
2172 * std::copy_n(&U[1], dim, &result[0]);
2176 *
template <
int dim>
2178 * ProblemDescription<dim>::internal_energy(
const state_type &U)
2180 *
const double &rho = U[0];
2181 *
const auto m = momentum(U);
2182 *
const double &
E = U[dim + 1];
2183 *
return E - 0.5 * m.norm_square() / rho;
2186 *
template <
int dim>
2188 * ProblemDescription<dim>::pressure(
const state_type &U)
2190 *
return (gamma - 1.) * internal_energy(U);
2193 *
template <
int dim>
2195 * ProblemDescription<dim>::speed_of_sound(
const state_type &U)
2197 *
const double &rho = U[0];
2198 *
const double p = pressure(U);
2203 *
template <
int dim>
2205 * ProblemDescription<dim>::flux(
const state_type &U)
2207 *
const double &rho = U[0];
2208 *
const auto m = momentum(U);
2209 *
const auto p = pressure(U);
2210 *
const double &
E = U[dim + 1];
2215 *
for (
unsigned int i = 0; i < dim; ++i)
2217 * result[1 + i] = m * m[i] / rho;
2218 * result[1 + i][i] += p;
2220 * result[dim + 1] = m / rho * (
E + p);
2227 * Now we discuss the computation of @f$\lambda_{\text{
max}}
2228 * (\mathbf{U}_i^{n},\mathbf{U}_j^{n}, \textbf{n}_{ij})@f$. The analysis
2229 * and derivation of sharp upper-bounds of maximum wavespeeds of Riemann
2230 * problems is a very technical endeavor and we cannot include an
2231 * advanced discussion about it in
this tutorial. In
this portion of the
2232 * documentation we will limit ourselves to sketch the main functionality
2233 * of our implementation functions and point to specific academic
2234 * references in order to help the (interested) reader
trace the
2235 * source (and proper mathematical justification) of these ideas.
2239 * In
general, obtaining a sharp guaranteed upper-bound on the maximum
2240 * wavespeed
requires solving a quite expensive
scalar nonlinear problem.
2241 * This is typically done with an iterative solver. In order to simplify
2242 * the presentation in
this example step we decided not to include such
2243 * an iterative scheme. Instead, we will just use an
initial guess as a
2244 * guess
for an upper bound on the maximum wavespeed. More precisely,
2245 * equations (2.11) (3.7), (3.8) and (4.3) of @cite GuermondPopov2016b
2246 * are enough to define a guaranteed upper bound on the maximum
2247 * wavespeed. This estimate is returned by a call to the function
2248 * <code>lambda_max_two_rarefaction()</code>. At its core the
2249 * construction of such an upper bound uses the so-called two-rarefaction
2250 * approximation for the intermediate pressure @f$p^*@f$, see for instance
2251 * Equation (4.46), page 128 in @cite Toro2009.
2255 * The estimate returned by <code>lambda_max_two_rarefaction()</code> is
2256 * guaranteed to be an upper bound, it is in general quite sharp, and
2257 * overall sufficient for our purposes. However, for some specific
2258 * situations (in particular when one of states is close to vacuum
2259 * conditions) such an estimate will be overly pessimistic. That's why we
2260 * used a
second estimate to avoid this degeneracy that will be invoked
2261 * by a call to the function <code>lambda_max_expansion()</code>. The most
2262 * important function here is <code>compute_lambda_max()</code> which
2263 * takes the minimum between the estimates returned by
2264 * <code>lambda_max_two_rarefaction()</code> and
2265 * <code>lambda_max_expansion()</code>.
2269 * We start again by defining a couple of helper functions:
2273 * The
first function takes a state <code>U</code> and a unit vector
2274 * <code>n_ij</code> and computes the <i>projected</i> 1d state in
2275 * direction of the unit vector.
2280 *
template <
int dim>
2282 *
const typename ProblemDescription<dim>::state_type U,
2286 * projected_U[0] = U[0];
2290 * For
this, we have to change the momentum to @f$\textbf{m}\cdot
2291 * n_{ij}@f$ and have to subtract the kinetic energy of the
2292 * perpendicular part from the total energy:
2295 *
const auto m = ProblemDescription<dim>::momentum(U);
2296 * projected_U[1] = n_ij * m;
2298 *
const auto perpendicular_m = m - projected_U[1] * n_ij;
2299 * projected_U[2] = U[1 + dim] - 0.5 * perpendicular_m.norm_square() / U[0];
2303 * We
return the 1
d state in <i>primitive</i> variables instead of
2304 * conserved quantities. The
return array consists of density @f$\rho@f$,
2305 * velocity @f$u@f$, pressure @f$p@f$ and local speed of sound @f$a@f$:
2311 *
return {{projected_U[0],
2312 * projected_U[1] / projected_U[0],
2313 * ProblemDescription<1>::pressure(projected_U),
2314 * ProblemDescription<1>::speed_of_sound(projected_U)}};
2319 * At
this point we also define two small
functions that
return the
2339 * Next, we need two local wavenumbers that are defined in terms of a
2340 * primitive state @f$[\rho, u, p, a]@f$ and a given pressure @f$p^\ast@f$
2341 * @cite GuermondPopov2016 Eqn. (3.7):
2344 * \left(\frac{p^\ast-p}{p}\right)_+}
2346 * Here, the @f$(\cdot)_{+}@f$ denotes the
positive part of the given
2354 * lambda1_minus(
const std::array<double, 4> &riemann_data,
2355 *
const double p_star)
2359 *
constexpr double gamma = ProblemDescription<1>::gamma;
2360 *
const auto u = riemann_data[1];
2361 *
const auto p = riemann_data[2];
2362 *
const auto a = riemann_data[3];
2364 *
const double factor = (
gamma + 1.0) / 2.0 / gamma;
2365 *
const double tmp = positive_part((p_star - p) / p);
2366 *
return u - a *
std::sqrt(1.0 + factor * tmp);
2371 * Analougously @cite GuermondPopov2016 Eqn. (3.8):
2374 * \left(\frac{p^\ast-p}{p}\right)_+}
2382 * lambda3_plus(
const std::array<double, 4> &riemann_data,
const double p_star)
2386 *
constexpr double gamma = ProblemDescription<1>::gamma;
2387 *
const auto u = riemann_data[1];
2388 *
const auto p = riemann_data[2];
2389 *
const auto a = riemann_data[3];
2391 *
const double factor = (
gamma + 1.0) / 2.0 / gamma;
2392 *
const double tmp = positive_part((p_star - p) / p);
2393 *
return u + a *
std::sqrt(1.0 + factor * tmp);
2398 * All that is left to
do is to compute the maximum of @f$\lambda^-@f$ and
2399 * @f$\lambda^+@f$ computed from the left and right primitive state
2400 * (@cite GuermondPopov2016 Eqn. (2.11)), where @f$p^\ast@f$ is given by
2401 * @cite GuermondPopov2016 Eqn (4.3):
2408 * lambda_max_two_rarefaction(const
std::array<double, 4> &riemann_data_i,
2409 * const
std::array<double, 4> &riemann_data_j)
2411 * constexpr double
gamma = ProblemDescription<1>::
gamma;
2412 *
const auto u_i = riemann_data_i[1];
2413 *
const auto p_i = riemann_data_i[2];
2414 *
const auto a_i = riemann_data_i[3];
2415 *
const auto u_j = riemann_data_j[1];
2416 *
const auto p_j = riemann_data_j[2];
2417 *
const auto a_j = riemann_data_j[3];
2419 *
const double numerator = a_i + a_j - (
gamma - 1.) / 2. * (u_j - u_i);
2421 *
const double denominator =
2422 * a_i *
std::pow(p_i / p_j, -1. * (gamma - 1.) / 2. / gamma) + a_j * 1.;
2426 *
const double p_star =
2427 * p_j *
std::pow(numerator / denominator, 2. * gamma / (gamma - 1));
2429 *
const double lambda1 = lambda1_minus(riemann_data_i, p_star);
2430 *
const double lambda3 = lambda3_plus(riemann_data_j, p_star);
2434 *
return std::max(positive_part(lambda3), negative_part(lambda1));
2439 * We compute the
second upper bound of the maximal wavespeed that is, in
2440 *
general, not as sharp as the two-rarefaction estimate. But it will
2441 * save the day in the context of near vacuum conditions when the
2442 * two-rarefaction approximation might attain extreme
values:
2444 * \lambda_{\text{
exp}} = \max(u_i,u_j) + 5. \max(a_i, a_j).
2446 * @note The
constant 5.0 multiplying the maximum of the sound speeds
2447 * is <i>neither</i> an ad-hoc
constant, <i>nor</i> a tuning parameter.
2448 * It defines an upper bound
for any @f$\gamma \in [0,5/3]@f$. Do not play
2456 * lambda_max_expansion(
const std::array<double, 4> &riemann_data_i,
2457 *
const std::array<double, 4> &riemann_data_j)
2459 *
const auto u_i = riemann_data_i[1];
2460 *
const auto a_i = riemann_data_i[3];
2461 *
const auto u_j = riemann_data_j[1];
2462 *
const auto a_j = riemann_data_j[3];
2470 * The following is the main function that we are going to call in order to
2471 * compute @f$\lambda_{\text{
max}} (\mathbf{U}_i^{n},\mathbf{U}_j^{n},
2472 * \textbf{n}_{ij})@f$. We simply compute both maximal wavespeed estimates
2473 * and
return the minimum.
2479 * template <int dim>
2481 * ProblemDescription<dim>::compute_lambda_max(
const state_type &U_i,
2482 *
const state_type &U_j,
2485 *
const auto riemann_data_i = riemann_data_from_state(U_i, n_ij);
2486 *
const auto riemann_data_j = riemann_data_from_state(U_j, n_ij);
2488 *
const double lambda_1 =
2489 * lambda_max_two_rarefaction(riemann_data_i, riemann_data_j);
2491 *
const double lambda_2 =
2492 * lambda_max_expansion(riemann_data_i, riemann_data_j);
2494 *
return std::min(lambda_1, lambda_2);
2499 * We conclude
this section by defining
static arrays
2500 * <code>component_names</code> that contain strings describing the
2501 * component names of our state vector. We have
template specializations
2502 *
for dimensions one, two and three, that are used later in
DataOut for
2503 * naming the corresponding components:
2510 *
const std::array<std::string, 3> ProblemDescription<1>::component_names{
2511 * {
"rho",
"m",
"E"}};
2514 *
const std::array<std::string, 4> ProblemDescription<2>::component_names{
2515 * {
"rho",
"m_1",
"m_2",
"E"}};
2518 *
const std::array<std::string, 5> ProblemDescription<3>::component_names{
2519 * {
"rho",
"m_1",
"m_2",
"m_3",
"E"}};
2524 * <a name=
"step_69-Initialvalues"></a>
2525 * <h4>Initial
values</h4>
2529 * The last preparatory step, before we discuss the implementation of the
2530 * forward Euler scheme, is to briefly implement the `InitialValues`
class.
2534 * In the constructor we initialize all parameters with
default values,
2536 * <code>parse_parameters_call_back</code> slot to the respective signal.
2540 * The <code>parse_parameters_call_back</code> slot will be invoked from
2542 * that regard, its use is appropriate
for situations where the
2543 * parameters have to be postprocessed (in some sense) or some
2544 * consistency condition between the parameters has to be checked.
2550 *
template <
int dim>
2551 * InitialValues<dim>::InitialValues(
const std::string &subsection)
2557 * [&]() { this->parse_parameters_callback(); });
2559 * initial_direction[0] = 1.;
2560 * add_parameter(
"initial direction",
2561 * initial_direction,
2562 *
"Initial direction of the uniform flow field");
2564 * initial_1d_state[0] = ProblemDescription<dim>::gamma;
2565 * initial_1d_state[1] = 3.;
2566 * initial_1d_state[2] = 1.;
2567 * add_parameter(
"initial 1d state",
2569 *
"Initial 1d state (rho, u, p) of the uniform flow field");
2574 * So far the constructor of <code>InitialValues</code> has defined
2575 *
default values for the two
private members
2576 * <code>initial_direction</code> and <code>initial_1d_state</code> and
2577 * added them to the parameter list. But we have not defined an
2578 * implementation of the only
public member that we really care about,
2579 * which is <code>initial_state()</code> (the function that we are going to
2580 * call to actually evaluate the
initial solution at the mesh nodes). At
2581 * the top of the function, we have to ensure that the provided initial
2582 * direction is not the zero vector.
2586 * @note As commented, we could have avoided
using the method
2587 * <code>parse_parameters_call_back </code> and defined a
class member
2588 * <code>setup()</code> in order to define the implementation of
2589 * <code>initial_state()</code>. But
for illustrative purposes we want to
2590 * document a different way here and use the call back signal from
2597 *
template <
int dim>
2598 *
void InitialValues<dim>::parse_parameters_callback()
2602 *
"Initial shock front direction is set to the zero vector."));
2603 * initial_direction /= initial_direction.norm();
2607 * Next, we implement the <code>initial_state</code> function
object
2608 * with a
lambda function computing a uniform flow field. For
this we
2609 * have to translate a given primitive 1
d state (density @f$\rho@f$,
2610 * velocity @f$u@f$, and pressure @f$p@f$) into a conserved n-dimensional state
2611 * (density @f$\rho@f$, momentum @f$\mathbf{m}@f$, and total energy @f$E@f$).
2617 * initial_state = [
this](
const Point<dim> & ,
double ) {
2618 *
const double rho = initial_1d_state[0];
2619 *
const double u = initial_1d_state[1];
2620 *
const double p = initial_1d_state[2];
2621 *
static constexpr double gamma = ProblemDescription<dim>::gamma;
2626 *
for (
unsigned int i = 0; i < dim; ++i)
2627 * state[1 + i] = rho * u * initial_direction[i];
2629 * state[dim + 1] = p / (
gamma - 1.) + 0.5 * rho * u * u;
2638 * <a name=
"step_69-TheForwardEulerstep"></a>
2639 * <h4>The Forward Euler step</h4>
2643 * The constructor of the <code>%
TimeStepping</code>
class does not contain
2644 * any surprising code:
2650 *
template <
int dim>
2654 *
const OfflineData<dim> &offline_data,
2655 *
const InitialValues<dim> &initial_values,
2656 *
const std::string &subsection )
2658 * , mpi_communicator(mpi_communicator)
2659 * , computing_timer(computing_timer)
2660 * , offline_data(&offline_data)
2661 * , initial_values(&initial_values)
2663 * cfl_update = 0.80;
2664 * add_parameter(
"cfl update",
2666 *
"Relative CFL constant used for update");
2671 * In the
class member <code>prepare()</code> we initialize the temporary
2672 * vector <code>temp</code> and the matrix <code>dij_matrix</code>. The
2673 * vector <code>temp</code> will be used to store the solution update
2674 * temporarily before its contents is swapped with the old vector.
2680 *
template <
int dim>
2681 *
void TimeStepping<dim>::prepare()
2684 *
"time_stepping - prepare scratch space");
2686 *
for (
auto &it : temporary_vector)
2687 * it.
reinit(offline_data->partitioner);
2689 * dij_matrix.reinit(offline_data->sparsity_pattern);
2694 * It is now time to implement the forward Euler step. Given a (writable
2695 * reference) to the old state <code>U</code> at time @f$t@f$ we update the
2696 * state <code>U</code> in place and
return the chosen time-step
size. We
2697 *
first declare a number of read-only references to various different
2698 * variables and
data structures. We
do this is mainly to have shorter
2699 * variable names (
e.g., <code>sparsity</code> instead of
2700 * <code>offline_data->sparsity_pattern</code>).
2706 *
template <
int dim>
2707 *
double TimeStepping<dim>::make_one_step(vector_type &U,
const double t)
2709 *
const auto &n_locally_owned = offline_data->n_locally_owned;
2710 *
const auto &n_locally_relevant = offline_data->n_locally_relevant;
2713 * indices_owned(0, n_locally_owned);
2715 * indices_relevant(0, n_locally_relevant);
2717 *
const auto &sparsity = offline_data->sparsity_pattern;
2719 *
const auto &lumped_mass_matrix = offline_data->lumped_mass_matrix;
2720 *
const auto &norm_matrix = offline_data->norm_matrix;
2721 *
const auto &nij_matrix = offline_data->nij_matrix;
2722 *
const auto &cij_matrix = offline_data->cij_matrix;
2724 *
const auto &boundary_normal_map = offline_data->boundary_normal_map;
2728 * <
b>Step 1</
b>: Computing the @f$d_{ij}@f$ graph viscosity
matrix.
2732 * It is important to highlight that the viscosity
matrix has to be
2733 *
symmetric, i.e., @f$d_{ij} = d_{ji}@f$. In
this regard we note here that
2734 * @f$\int_{\Omega} \nabla \phi_j \phi_i \, \mathrm{
d}\mathbf{x}= -
2735 * \int_{\Omega} \nabla \phi_i \phi_j \, \mathrm{
d}\mathbf{x}@f$ (or
2736 * equivalently @f$\mathbf{c}_{ij} = - \mathbf{c}_{ji}@f$) provided either
2737 * @f$\mathbf{x}_i@f$ or @f$\mathbf{x}_j@f$ is a support
point located away
2738 * from the boundary. In
this case we can
check that
2739 * @f$\lambda_{\text{
max}} (\mathbf{U}_i^{n}, \mathbf{U}_j^{n},
2740 * \textbf{n}_{ij}) = \lambda_{\text{
max}} (\mathbf{U}_j^{n},
2741 * \mathbf{U}_i^{n},\textbf{n}_{ji})@f$ by construction, which guarantees
2742 * the property @f$d_{ij} = d_{ji}@f$.
2746 * However,
if both support points @f$\mathbf{x}_i@f$ or @f$\mathbf{x}_j@f$
2747 * happen to lie on the boundary, then, the equalities @f$\mathbf{c}_{ij} =
2748 * - \mathbf{c}_{ji}@f$ and @f$\lambda_{\text{
max}} (\mathbf{U}_i^{n},
2749 * \mathbf{U}_j^{n}, \textbf{n}_{ij}) = \lambda_{\text{
max}}
2750 * (\mathbf{U}_j^{n}, \mathbf{U}_i^{n}, \textbf{n}_{ji})@f$
do not
2751 * necessarily hold
true. The only mathematically safe solution
for this
2752 * dilemma is to compute both of them @f$d_{ij}@f$ and @f$d_{ji}@f$ and
2757 * Overall, the computation of @f$d_{ij}@f$ is quite expensive. In
2758 * order to save some computing time we exploit the fact that the viscosity
2760 * the upper-triangular entries of @f$d_{ij}@f$ and
copy the
2761 * corresponding entries to the lower-triangular counterpart.
2766 * loops. Pretty much all the ideas for
parallel traversal that we
2767 * introduced when discussing the assembly of the
matrix
2768 * <code>norm_matrix</code> and the normalization of
2769 * <code>nij_matrix</code> above are used here again.
2773 * We define again a
"worker" function <code>on_subranges</code> that
2774 * computes the viscosity @f$d_{ij}@f$ for a subrange `[*index_begin,
2775 * *index_end)` of column indices:
2780 *
"time_stepping - 1 compute d_ij");
2782 *
const auto on_subranges =
2783 * [&](
const auto index_begin,
const auto index_end) {
2784 *
for (
const auto i :
2786 * *index_begin, *index_end))
2788 * const auto U_i =
gather(U, i);
2792 * For a given column
index i we iterate over the columns of the
2793 * sparsity pattern from <code>sparsity.begin(i)</code> to
2794 * <code>sparsity.end(i)</code>:
2797 *
for (
auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
2799 *
const auto j = jt->column();
2803 * We only compute @f$d_{ij}@f$
if @f$j < i@f$ (upper triangular
2804 * entries) and later
copy the
values over to @f$d_{ji}@f$.
2810 *
const auto U_j =
gather(U, j);
2812 *
const auto n_ij = gather_get_entry(nij_matrix, jt);
2813 *
const double norm = get_entry(norm_matrix, jt);
2815 *
const auto lambda_max =
2816 * ProblemDescription<dim>::compute_lambda_max(U_i, U_j, n_ij);
2818 *
double d =
norm * lambda_max;
2822 * If both support points happen to be at the boundary we
2823 * have to compute @f$d_{ji}@f$ as well and then take
2824 * @f$\max(d_{ij},d_{ji})@f$. After
this we can
finally set the
2825 * upper triangular and lower triangular entries.
2828 * if (boundary_normal_map.count(i) != 0 &&
2829 * boundary_normal_map.count(j) != 0)
2831 *
const auto n_ji =
gather(nij_matrix, j, i);
2832 *
const auto lambda_max_2 =
2833 * ProblemDescription<dim>::compute_lambda_max(U_j,
2836 *
const double norm_2 = norm_matrix(j, i);
2838 *
d =
std::max(d, norm_2 * lambda_max_2);
2841 * set_entry(dij_matrix, jt, d);
2842 * dij_matrix(j, i) =
d;
2848 * indices_relevant.end(),
2855 * <
b>Step 2</
b>: Compute
diagonal entries @f$d_{ii}@f$ and
2856 * @f$\tau_{\text{
max}}@f$.
2861 * <code>dij_matrix</code>. We still have to fill its
diagonal entries
2862 * defined as @f$d_{ii}^n = - \sum_{j \in \mathcal{I}(i)\backslash \{i\}}
2864 * purpose. While computing the @f$d_{ii}@f$s we also determine the
2865 * largest admissible time-step, which is defined as
2867 * \tau_n \dealcoloneq c_{\text{cfl}}\,\min_{i\in\mathcal{V}}
2868 * \left(\frac{m_i}{-2\,d_{ii}^{n}}\right) \, .
2870 * Note that the operation @f$\min_{i \in \mathcal{V}}@f$ is intrinsically
2871 * global, it operates on all nodes:
first we have to take the minimum
2872 * over all threads (of a given node) and then we have to take the
2873 * minimum over all
MPI processes. In the current implementation:
2874 * - We store <code>tau_max</code> (per node) as
2876 * href=
"http://www.cplusplus.com/reference/atomic/atomic/"><code>std::atomic<double></code></a>.
2877 * The
internal implementation of <code>std::atomic</code> will take
2878 * care of guarding any possible race condition when more than one
2879 * thread attempts to read and/or write <code>tau_max</code> at the
2881 * - In order to take the minimum over all
MPI process we use the utility
2882 * function <code>Utilities::MPI::min</code>.
2888 * std::atomic<double> tau_max{std::numeric_limits<double>::infinity()};
2892 *
"time_stepping - 2 compute d_ii, and tau_max");
2896 * on_subranges() will be executed on every thread individually. The
2897 * variable <code>tau_max_on_subrange</code> is thus stored thread
2904 * const auto on_subranges =
2905 * [&](const auto index_begin, const auto index_end) {
2906 *
double tau_max_on_subrange = std::numeric_limits<double>::infinity();
2908 *
for (
const auto i :
2910 * *index_begin, *index_end))
2912 * double d_sum = 0.;
2914 *
for (
auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
2916 *
const auto j = jt->column();
2921 * d_sum -= get_entry(dij_matrix, jt);
2926 * We store the
negative sum of the d_ij entries at the
2930 * dij_matrix.diag_element(i) = d_sum;
2933 * and compute the maximal local time-step
size
2937 *
const double mass = lumped_mass_matrix.diag_element(i);
2938 *
const double tau = cfl_update * mass / (-2. * d_sum);
2939 * tau_max_on_subrange =
std::min(tau_max_on_subrange, tau);
2944 * <code>tau_max_on_subrange</code> contains the largest possible
2945 * time-step
size computed
for the (thread local) subrange. At
this
2946 *
point we have to synchronize the
value over all threads. This is
2947 * were we use the <a
2948 * href=
"http://www.cplusplus.com/reference/atomic/atomic/"><code>std::atomic<double></code></a>
2949 * <i>compare exchange</i> update mechanism:
2952 *
double current_tau_max = tau_max.load();
2953 *
while (current_tau_max > tau_max_on_subrange &&
2954 * !tau_max.compare_exchange_weak(current_tau_max,
2955 * tau_max_on_subrange))
2960 * indices_relevant.end(),
2966 * After all threads have finished we can simply synchronize the
2977 * This is a good
point to verify that the computed
2978 * <code>tau_max</code> is indeed a
valid floating
point number.
2982 * !std::isnan(tau_max.load()) && !std::isinf(tau_max.load()) &&
2983 * tau_max.load() > 0.,
2985 *
"I'm sorry, Dave. I'm afraid I can't do that. - We crashed."));
2990 * <
b>Step 3</
b>: Perform update.
2994 * At
this point, we have computed all viscosity coefficients @f$d_{ij}@f$
2995 * and we know the maximal admissible time-step
size
2996 * @f$\tau_{\text{
max}}@f$. This means we can now compute the update:
3001 * \mathbf{U}_i^{n+1} = \mathbf{U}_i^{n} - \frac{\tau_{\text{
max}} }{m_i}
3002 * \sum_{j \in \mathcal{I}(i)} (\mathbb{f}(\mathbf{U}_j^{n}) -
3003 * \mathbb{f}(\mathbf{U}_i^{n})) \cdot \mathbf{c}_{ij} - d_{ij}
3004 * (\mathbf{U}_j^{n} - \mathbf{U}_i^{n})
3009 * This update formula is slightly different from what was discussed in
3010 * the introduction (in the pseudo-code). However, it can be shown that
3011 * both equations are algebraically equivalent (they will produce the
3012 * same numerical values). We favor
this second formula since it has
3013 * natural cancellation properties that might help avoid numerical
3022 *
"time_stepping - 3 perform update");
3024 *
const auto on_subranges =
3025 * [&](
const auto index_begin,
const auto index_end) {
3026 *
for (
const auto i :
3029 *
Assert(i < n_locally_owned, ExcInternalError());
3031 *
const auto U_i =
gather(U, i);
3033 *
const auto f_i = ProblemDescription<dim>::flux(U_i);
3034 *
const double m_i = lumped_mass_matrix.diag_element(i);
3036 *
auto U_i_new = U_i;
3038 *
for (
auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
3040 *
const auto j = jt->column();
3042 *
const auto U_j =
gather(U, j);
3043 *
const auto f_j = ProblemDescription<dim>::flux(U_j);
3045 *
const auto c_ij = gather_get_entry(cij_matrix, jt);
3046 *
const auto d_ij = get_entry(dij_matrix, jt);
3048 *
for (
unsigned int k = 0; k < n_solution_variables; ++k)
3052 * (-(f_j[k] - f_i[k]) * c_ij + d_ij * (U_j[k] - U_i[k]));
3056 *
scatter(temporary_vector, U_i_new, i);
3061 * indices_owned.end(),
3068 * <
b>Step 4</
b>: Fix up boundary states.
3072 * As a last step in the Forward Euler method, we have to fix up all
3073 * boundary states. As discussed in the intro we
3074 * -
advance in time satisfying no boundary condition at all,
3075 * - at the
end of the time step enforce boundary conditions strongly
3076 * in a post-processing step.
3080 * Here, we compute the correction
3082 * \mathbf{m}_i \dealcoloneq \mathbf{m}_i - (\boldsymbol{\nu}_i \cdot
3083 * \mathbf{m}_i) \boldsymbol{\nu}_i,
3085 * which removes the normal component of @f$\mathbf{m}@f$.
3093 *
"time_stepping - 4 fix boundary states");
3095 *
for (
const auto &it : boundary_normal_map)
3097 * const auto i = it.
first;
3101 * We only iterate over the locally owned subset:
3104 *
if (i >= n_locally_owned)
3107 *
const auto &normal = std::get<0>(it.second);
3108 *
const auto &
id = std::get<1>(it.second);
3109 *
const auto &position = std::get<2>(it.second);
3111 *
auto U_i =
gather(temporary_vector, i);
3115 * On free slip boundaries we remove the normal component of the
3119 *
if (
id == Boundaries::free_slip)
3121 *
auto m = ProblemDescription<dim>::momentum(U_i);
3122 * m -= (m * normal) * normal;
3123 *
for (
unsigned int k = 0; k < dim; ++k)
3124 * U_i[k + 1] = m[k];
3129 * On Dirichlet boundaries we enforce
initial conditions
3133 *
else if (
id == Boundaries::dirichlet)
3135 * U_i = initial_values->initial_state(position, t + tau_max);
3138 *
scatter(temporary_vector, U_i, i);
3144 * <
b>Step 5</
b>: We now update the ghost layer over all
MPI ranks,
3145 *
swap the temporary vector with the solution vector <code>U</code>
3146 * (that will get returned by reference) and
return the chosen
3147 * time-step
size @f$\tau_{\text{
max}}@f$:
3153 *
for (
auto &it : temporary_vector)
3154 * it.update_ghost_values();
3156 * U.swap(temporary_vector);
3164 * <a name=
"step_69-Schlierenpostprocessing"></a>
3165 * <h4>Schlieren postprocessing</h4>
3169 * At various intervals we will output the current state <code>U</code>
3170 * of the solution together with a so-called Schlieren plot.
3171 * The constructor of the <code>SchlierenPostprocessor</code>
class again
3172 * contains no surprises. We simply supply
default values to and
register
3175 * is an ad-hoc
positive amplification factor in order to enhance the
3176 * contrast in the visualization. Its actual
value is a matter of
3178 * - schlieren_index: is an integer indicating which component of the
3179 * state @f$[\rho, \mathbf{m},
E]@f$ are we going to use in order to generate
3180 * the visualization.
3186 *
template <
int dim>
3187 * SchlierenPostprocessor<dim>::SchlierenPostprocessor(
3190 *
const OfflineData<dim> &offline_data,
3191 *
const std::string &subsection )
3193 * , mpi_communicator(mpi_communicator)
3194 * , computing_timer(computing_timer)
3195 * , offline_data(&offline_data)
3197 * schlieren_beta = 10.;
3198 * add_parameter(
"schlieren beta",
3200 *
"Beta factor used in Schlieren-type postprocessor");
3202 * schlieren_index = 0;
3203 * add_parameter(
"schlieren index",
3205 *
"Use the corresponding component of the state vector for the "
3206 *
"schlieren plot");
3211 * Again, the <code>prepare()</code> function initializes two temporary
3212 * the vectors (<code>r</code> and <code>schlieren</code>).
3218 *
template <
int dim>
3219 *
void SchlierenPostprocessor<dim>::prepare()
3222 *
"schlieren_postprocessor - prepare scratch space");
3224 * r.reinit(offline_data->n_locally_relevant);
3225 * schlieren.reinit(offline_data->partitioner);
3230 * We now discuss the implementation of the
class member
3231 * <code>SchlierenPostprocessor<dim>::compute_schlieren()</code>, which
3232 * basically takes a component of the state vector <code>U</code> and
3233 * computes the Schlieren indicator
for such component (the formula of the
3234 * Schlieren indicator can be found just before the declaration of the
class
3235 * <code>SchlierenPostprocessor</code>). We start by noting
3236 * that
this formula
requires the
"nodal gradients" @f$\nabla r_j@f$.
3237 * However, nodal
values of
gradients are not defined
for @f$\mathcal{
C}^0@f$
3240 * simplest technique we can use to recover gradients at nodes is
3241 * weighted-averaging i.e.
3245 * \f[ \nabla r_j \dealcoloneq \frac{1}{\int_{S_i} \omega_i(\mathbf{x}) \,
3246 * \mathrm{
d}\mathbf{x}}
3247 * \int_{S_i} r_h(\mathbf{x}) \omega_i(\mathbf{x}) \, \mathrm{
d}\mathbf{x}
3248 * \ \ \ \ \ \mathbf{(*)} \f]
3252 * where @f$S_i@f$ is the support of the shape function @f$\phi_i@f$, and
3253 * @f$\omega_i(\mathbf{x})@f$ is the weight. The weight could be any
3254 * positive function such as
3255 * @f$\omega_i(\mathbf{x}) \equiv 1@f$ (that would allow us to recover the usual
3256 * notion of mean value). But as usual, the goal is to reuse the off-line
3257 *
data as much as possible. In
this sense, the most natural
3258 * choice of weight is @f$\omega_i = \phi_i@f$. Inserting
this choice of
3259 * weight and the expansion @f$r_h(\mathbf{x}) = \sum_{j \in \mathcal{V}}
3260 * r_j \phi_j(\mathbf{x})@f$ into @f$\mathbf{(*)}@f$ we get :
3264 * \f[ \nabla r_j \dealcoloneq \frac{1}{m_i} \sum_{j \in \mathcal{I}(i)} r_j
3265 * \mathbf{c}_{ij} \ \ \ \ \ \mathbf{(**)} \, . \f]
3269 * Using
this last formula we can recover averaged nodal
gradients without
3270 * resorting to any form of quadrature. This idea aligns quite well with
3271 * the whole spirit of edge-based schemes (or algebraic schemes) where
3272 * we want to operate on matrices and vectors as directly as
3273 * it could be possible avoiding by all means assembly of bilinear
3274 * forms, cell-loops, quadrature, or any other
3275 * intermediate construct/operation between the input arguments (the state
3276 * from the previous time-step) and the actual matrices and vectors
3277 * required to compute the update.
3281 * The
second thing to note is that we have to compute global minimum and
3282 * maximum @f$\max_j |\nabla r_j|@f$ and @f$\min_j |\nabla r_j|@f$. Following the
3283 * same ideas used to compute the time step
size in the
class member
3284 * <code>%
TimeStepping\<dim>::%step()</code> we define @f$\max_j |\nabla r_j|@f$
3285 * and @f$\min_j |\nabla r_j|@f$ as atomic doubles in order to resolve any
3286 * conflicts between threads. As usual, we use
3289 * among all
MPI processes.
3293 * Finally, it is not possible to compute the Schlieren indicator in a single
3294 *
loop over all nodes. The entire operation
requires two loops over nodes:
3298 * - The
first loop computes @f$|\nabla r_i|@f$
for all @f$i \in \mathcal{V}@f$ in
3299 * the mesh, and the bounds @f$\max_j |\nabla r_j|@f$ and
3300 * @f$\min_j |\nabla r_j|@f$.
3301 * - The
second loop finally computes the Schlieren indicator
using the
3306 * \f[ \text{schlieren}[i] =
e^{\beta \frac{ |\nabla r_i|
3307 * - \min_j |\nabla r_j| }{\max_j |\nabla r_j| - \min_j |\nabla r_j| } }
3312 * This means that we will have to define two workers
3313 * <code>on_subranges</code>
for each one of these stages.
3319 *
template <
int dim>
3320 *
void SchlierenPostprocessor<dim>::compute_schlieren(
const vector_type &U)
3323 * computing_timer,
"schlieren_postprocessor - compute schlieren plot");
3325 *
const auto &sparsity = offline_data->sparsity_pattern;
3326 *
const auto &lumped_mass_matrix = offline_data->lumped_mass_matrix;
3327 *
const auto &cij_matrix = offline_data->cij_matrix;
3328 *
const auto &boundary_normal_map = offline_data->boundary_normal_map;
3329 *
const auto &n_locally_owned = offline_data->n_locally_owned;
3331 *
const auto indices =
3337 * We define the r_i_max and r_i_min in the current
MPI process as
3338 * atomic doubles in order to avoid race conditions between threads:
3341 * std::atomic<double> r_i_max{0.};
3342 * std::atomic<double> r_i_min{std::numeric_limits<double>::infinity()};
3346 * First
loop: compute the averaged
gradient at each node and the
3347 * global maxima and minima of the
gradients.
3351 *
const auto on_subranges =
3352 * [&](
const auto index_begin,
const auto index_end) {
3353 *
double r_i_max_on_subrange = 0.;
3354 *
double r_i_min_on_subrange = std::numeric_limits<double>::infinity();
3356 *
for (
const auto i :
3359 *
Assert(i < n_locally_owned, ExcInternalError());
3363 *
for (
auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
3365 *
const auto j = jt->column();
3370 *
const auto U_js = U[schlieren_index].local_element(j);
3371 *
const auto c_ij = gather_get_entry(cij_matrix, jt);
3372 * r_i += c_ij * U_js;
3377 * We fix up the
gradient r_i at free slip boundaries similarly to
3378 * how we fixed up boundary states in the forward Euler step.
3379 * This avoids sharp, artificial
gradients in the Schlieren
3380 * plot at free slip boundaries and is a purely cosmetic choice.
3386 *
const auto bnm_it = boundary_normal_map.find(i);
3387 *
if (bnm_it != boundary_normal_map.end())
3389 *
const auto &normal = std::get<0>(bnm_it->second);
3390 *
const auto &
id = std::get<1>(bnm_it->second);
3392 *
if (
id == Boundaries::free_slip)
3393 * r_i -= 1. * (r_i * normal) * normal;
3400 * We remind the reader that we are not interested in the nodal
3401 *
gradients per se. We only want their norms in order to
3402 * compute the Schlieren indicator (weighted with the lumped
3403 * mass matrix @f$m_i@f$):
3406 * const double m_i = lumped_mass_matrix.diag_element(i);
3407 * r[i] = r_i.norm() / m_i;
3408 * r_i_max_on_subrange =
std::max(r_i_max_on_subrange, r[i]);
3409 * r_i_min_on_subrange =
std::min(r_i_min_on_subrange, r[i]);
3414 * We compare the current_r_i_max and current_r_i_min (in the
3415 * current subrange) with r_i_max and r_i_min (
for the current
MPI
3416 * process) and update them
if necessary:
3422 *
double current_r_i_max = r_i_max.load();
3423 *
while (current_r_i_max < r_i_max_on_subrange &&
3424 * !r_i_max.compare_exchange_weak(current_r_i_max,
3425 * r_i_max_on_subrange))
3428 *
double current_r_i_min = r_i_min.load();
3429 *
while (current_r_i_min > r_i_min_on_subrange &&
3430 * !r_i_min.compare_exchange_weak(current_r_i_min,
3431 * r_i_min_on_subrange))
3443 * And synchronize <code>r_i_max</code> and <code>r_i_min</code> over
3444 * all
MPI processes.
3455 * Second
loop: we now have the vector <code>r</code> and the scalars
3456 * <code>r_i_max</code> and <code>r_i_min</code> at our disposal. We
3457 * are thus in a position to actually compute the Schlieren indicator.
3464 *
const auto on_subranges =
3465 * [&](
const auto index_begin,
const auto index_end) {
3466 *
for (
const auto i :
3469 *
Assert(i < n_locally_owned, ExcInternalError());
3471 * schlieren.local_element(i) =
3472 * 1. -
std::exp(-schlieren_beta * (r[i] - r_i_min) /
3473 * (r_i_max - r_i_min));
3485 * And
finally, exchange ghost elements.
3488 * schlieren.update_ghost_values();
3494 * <a name=
"step_69-Themainloop"></a>
3495 * <h4>The main
loop</h4>
3499 * With all classes implemented it is time to create an instance of
3500 * <code>Discretization<dim></code>, <code>OfflineData<dim></code>,
3501 * <code>InitialValues<dim></code>, <code>%
TimeStepping\<dim></code>, and
3502 * <code>SchlierenPostprocessor<dim></code>, and
run the forward Euler
3507 * In the constructor of <code>MainLoop<dim></code> we now initialize an
3508 * instance of all classes, and declare a number of parameters
3509 * controlling output. Most notable, we declare a
boolean parameter
3510 * <code>resume</code> that will control whether the program attempts to
3511 * restart from an interrupted computation, or not.
3517 *
template <
int dim>
3518 * MainLoop<dim>::MainLoop(
const MPI_Comm mpi_communicator)
3520 * , mpi_communicator(mpi_communicator)
3521 * , computing_timer(mpi_communicator,
3526 * , discretization(mpi_communicator, computing_timer,
"B - Discretization")
3527 * , offline_data(mpi_communicator,
3530 *
"C - OfflineData")
3531 * , initial_values(
"D - InitialValues")
3532 * , time_stepping(mpi_communicator,
3536 *
"E - TimeStepping")
3537 * , schlieren_postprocessor(mpi_communicator,
3540 *
"F - SchlierenPostprocessor")
3542 * base_name =
"test";
3543 * add_parameter(
"basename", base_name,
"Base name for all output files");
3546 * add_parameter(
"final time", t_final,
"Final time");
3548 * output_granularity = 0.02;
3549 * add_parameter(
"output granularity",
3550 * output_granularity,
3551 *
"time interval for output");
3553 * asynchronous_writeback =
true;
3554 * add_parameter(
"asynchronous writeback",
3555 * asynchronous_writeback,
3556 *
"Write out solution in a background thread performing IO");
3559 * add_parameter(
"resume", resume,
"Resume an interrupted computation.");
3564 * We start by implementing a helper function <code>print_head()</code>
3565 * in an anonymous
namespace that is used to output messages in the
3566 * terminal with some nice formatting.
3575 *
const std::string &header,
3576 *
const std::string &secondary =
"")
3578 *
const auto header_size = header.size();
3579 *
const auto padded_header = std::string((34 - header_size) / 2,
' ') +
3581 * std::string((35 - header_size) / 2,
' ');
3583 *
const auto secondary_size = secondary.size();
3584 *
const auto padded_secondary =
3585 * std::string((34 - secondary_size) / 2,
' ') + secondary +
3586 * std::string((35 - secondary_size) / 2,
' ');
3589 * pcout << std::endl;
3590 * pcout <<
" ####################################################" << std::endl;
3591 * pcout <<
" ######### #########" << std::endl;
3592 * pcout <<
" #########" << padded_header <<
"#########" << std::endl;
3593 * pcout <<
" #########" << padded_secondary <<
"#########" << std::endl;
3594 * pcout <<
" ######### #########" << std::endl;
3595 * pcout <<
" ####################################################" << std::endl;
3596 * pcout << std::endl;
3603 * With <code>print_head</code> in place it is now time to implement the
3604 * <code>MainLoop<dim>::run()</code> that contains the main
loop of our
3611 *
template <
int dim>
3612 *
void MainLoop<dim>::run()
3616 * We start by reading in parameters and initializing all objects. We
3618 * all parameters from the parameter file (whose name is given as a
3621 * declarations
for all
class instances that are derived from
3622 * ParameterAceptor. The call to initialize enters the subsection
for
3623 * each each derived
class, and sets all variables that were added
3630 * pcout <<
"Reading parameters and allocating objects... " << std::flush;
3633 * pcout <<
"done" << std::endl;
3638 * scratch space, and initialize the
DataOut<dim> object. All of these
3639 * operations are pretty standard and discussed in detail in the
3640 * Discretization and OfflineData classes.
3647 * print_head(pcout,
"create triangulation");
3649 * discretization.setup();
3652 * discretization.triangulation.load(base_name +
"-checkpoint.mesh");
3654 * discretization.triangulation.refine_global(discretization.refinement);
3656 * pcout <<
"Number of active cells: "
3657 * << discretization.triangulation.n_global_active_cells()
3660 * print_head(pcout,
"compute offline data");
3661 * offline_data.setup();
3662 * offline_data.assemble();
3664 * pcout <<
"Number of degrees of freedom: "
3665 * << offline_data.dof_handler.n_dofs() << std::endl;
3667 * print_head(pcout,
"set up time step");
3668 * time_stepping.prepare();
3669 * schlieren_postprocessor.prepare();
3674 * We will store the current time and state in the variable
3675 * <code>t</code> and vector <code>U</code>:
3682 *
unsigned int output_cycle = 0;
3685 *
for (
auto &it : U)
3686 * it.
reinit(offline_data.partitioner);
3691 * <a name=
"step_69-Resume"></a>
3696 * By
default the boolean <code>resume</code> is set to
false, i.e. the
3697 * following code snippet is not
run. However,
if
3698 * <code>resume==
true</code> we indicate that we have indeed an
3699 * interrupted computation and the program shall restart by reading in
3700 * an old state consisting of <code>t</code>,
3701 * <code>output_cycle</code>, and the state vector <code>U</code> from
3706 * A
this point we have already read in the stored refinement history
3707 * of our
parallel distributed mesh. What is missing are the actual
3708 * state vector <code>U</code>, the time and output cycle. We use the
3710 * distributed::Triangulation::load() /
3711 * distributed::Triangulation::save() mechanism to read in the state
3712 * vector. A separate <code>
boost::archive</code> is used to retrieve
3713 * <code>t</code> and <code>output_cycle</code>. The checkpoint files
3714 * will be created in the <code>output()</code> routine discussed
3723 * print_head(pcout,
"resume interrupted computation");
3727 * solution_transfer(offline_data.dof_handler);
3729 * std::vector<LinearAlgebra::distributed::Vector<double> *> vectors;
3730 * std::transform(U.begin(),
3732 * std::back_inserter(vectors),
3733 * [](
auto &it) { return ⁢ });
3734 * solution_transfer.deserialize(vectors);
3736 *
for (
auto &it : U)
3737 * it.update_ghost_values();
3739 * std::ifstream file(base_name +
"-checkpoint.metadata",
3740 * std::ios::binary);
3742 * boost::archive::binary_iarchive ia(file);
3743 * ia >> t >> output_cycle;
3747 * print_head(pcout,
"interpolate initial values");
3748 * U = interpolate_initial_values();
3753 * With either the
initial state set up, or an interrupted state
3754 * restored it is time to enter the main
loop:
3760 * output(U, base_name, t, output_cycle++);
3762 * print_head(pcout,
"enter main loop");
3764 *
unsigned int timestep_number = 1;
3765 *
while (t < t_final)
3769 * We
first print an informative status message
3775 * std::ostringstream head;
3776 * std::ostringstream secondary;
3780 * << std::fixed << std::setprecision(1) << t / t_final * 100
3782 * secondary <<
"at time t = " << std::setprecision(8) << std::fixed << t;
3784 * print_head(pcout, head.str(), secondary.str());
3788 * and then perform a single forward Euler step. Note that the
3789 * state vector <code>U</code> is updated in place and that
3790 * <code>time_stepping.make_one_step()</code> returns the chosen step
3797 * t += time_stepping.make_one_step(U, t);
3801 * Post processing, generating output and writing out the current
3802 * state is a CPU and IO intensive task that we cannot afford to
do
3803 * every time step - in particular with
explicit time stepping. We
3804 * thus only schedule output by calling the <code>output()</code>
3805 * function
if we are past a threshold set by
3806 * <code>output_granularity</code>.
3812 *
if (t > output_cycle * output_granularity)
3814 * checkpoint(U, base_name, t, output_cycle);
3815 * output(U, base_name, t, output_cycle);
3819 * ++timestep_number;
3824 * We wait
for any remaining background output thread to finish before
3825 * printing a summary and exiting.
3828 *
if (background_thread_state.valid())
3829 * background_thread_state.wait();
3831 * computing_timer.print_summary();
3832 * pcout << timer_output.str() << std::endl;
3837 * The <code>interpolate_initial_values</code> takes an
initial time
"t"
3838 * as input argument and populates a state vector <code>U</code> with the
3839 * help of the <code>InitialValues<dim>::initial_state</code>
object.
3845 *
template <
int dim>
3846 *
typename MainLoop<dim>::vector_type
3847 * MainLoop<dim>::interpolate_initial_values(
const double t)
3849 * pcout <<
"MainLoop<dim>::interpolate_initial_values(t = " << t <<
')'
3852 *
"main_loop - setup scratch space");
3856 *
for (
auto &it : U)
3857 * it.
reinit(offline_data.partitioner);
3859 *
constexpr auto n_solution_variables =
3860 * ProblemDescription<dim>::n_solution_variables;
3864 * The function signature of
3865 * <code>InitialValues<dim>::initial_state</code> is not quite right
3867 * creating a
lambda function that
for a given position <code>x</code>
3868 * returns just the
value of the <code>i</code>th component. This
3873 *
for (
unsigned int i = 0; i < n_solution_variables; ++i)
3877 * return initial_values.initial_state(x, t)[i];
3881 *
for (
auto &it : U)
3882 * it.update_ghost_values();
3890 * <a name=
"step_69-Outputandcheckpointing"></a>
3891 * <h5>Output and checkpointing</h5>
3895 * We checkpoint the current state by doing the precise inverse
3896 * operation to what we discussed
for the <a href=
"Resume">resume
3903 *
template <
int dim>
3904 *
void MainLoop<dim>::checkpoint(
const typename MainLoop<dim>::vector_type &U,
3905 *
const std::string &name,
3907 *
const unsigned int cycle)
3909 * print_head(pcout,
"checkpoint computation");
3913 * solution_transfer(offline_data.dof_handler);
3915 * std::vector<const LinearAlgebra::distributed::Vector<double> *> vectors;
3916 * std::transform(U.begin(),
3918 * std::back_inserter(vectors),
3919 * [](
auto &it) { return ⁢ });
3921 * solution_transfer.prepare_for_serialization(vectors);
3923 * discretization.triangulation.save(name +
"-checkpoint.mesh");
3927 * std::ofstream file(name +
"-checkpoint.metadata", std::ios::binary);
3928 * boost::archive::binary_oarchive oa(file);
3935 * Writing out the
final vtk files is quite an IO intensive task that can
3936 * stall the main
loop for a
while. In order to avoid
this we use an <a
3937 * href=
"https://en.wikipedia.org/wiki/Asynchronous_I/O">asynchronous
3938 * IO</a> strategy by creating a background thread that will perform IO
3939 *
while the main
loop is allowed to
continue. In order
for this to work
3940 * we have to be mindful of two things:
3941 * - Before running the <code>output_worker</code> thread, we have to create
3942 * a
copy of the state vector <code>U</code>. We store it in the
3943 * vector <code>output_vector</code>.
3944 * - We have to avoid any
MPI communication in the background thread,
3945 * otherwise the program might deadlock. This implies that we have to
3946 *
run the postprocessing
outside of the worker thread.
3952 *
template <
int dim>
3953 *
void MainLoop<dim>::output(
const typename MainLoop<dim>::vector_type &U,
3954 *
const std::string &name,
3956 *
const unsigned int cycle)
3958 * pcout <<
"MainLoop<dim>::output(t = " << t <<
')' << std::endl;
3962 * If the asynchronous writeback option is set we launch a background
3963 * thread performing all the slow IO to disc. In that
case we have to
3964 * make sure that the background thread actually finished running. If
3965 * not, we have to wait to
for it to finish. We launch said background
3967 * href=
"https://en.cppreference.com/w/cpp/thread/async"><code>std::async()</code></a>
3969 * href=
"https://en.cppreference.com/w/cpp/thread/future"><code>std::future</code></a>
3970 *
object. This <code>std::future</code>
object contains the
return
3971 *
value of the function, which is in our
case simply
3972 * <code>
void</code>.
3978 *
if (background_thread_state.valid())
3981 * background_thread_state.wait();
3984 *
constexpr auto n_solution_variables =
3985 * ProblemDescription<dim>::n_solution_variables;
3989 * At
this point we make a
copy of the state vector,
run the schlieren
3991 * output code is standard: We create a
DataOut instance, attach all
3992 *
data vectors we want to output and call
3993 *
DataOut<dim>::build_patches(). There is one twist, however. In order
3994 * to perform asynchronous IO on a background thread we create the
3995 *
DataOut<dim>
object as a shared pointer that we pass on to the
3996 * worker thread to ensure that once we exit this function and the
3997 * worker thread finishes the
DataOut<dim>
object gets destroyed again.
4003 * for (
unsigned int i = 0; i < n_solution_variables; ++i)
4005 * output_vector[i] = U[i];
4006 * output_vector[i].update_ghost_values();
4009 * schlieren_postprocessor.compute_schlieren(output_vector);
4011 * std::unique_ptr<DataOut<dim>> data_out = std::make_unique<DataOut<dim>>();
4014 *
for (
unsigned int i = 0; i < n_solution_variables; ++i)
4015 * data_out->add_data_vector(output_vector[i],
4016 * ProblemDescription<dim>::component_names[i]);
4018 * data_out->add_data_vector(schlieren_postprocessor.schlieren,
4019 *
"schlieren_plot");
4021 * data_out->build_patches(discretization.mapping,
4022 * discretization.finite_element.degree - 1);
4026 * Next we create a
lambda function
for the background thread. We <a
4027 * href=
"https://en.cppreference.com/w/cpp/language/lambda">capture</a>
4028 * the <code>
this</code> pointer as well as most of the arguments of
4029 * the output function by
value so that we have access to them
inside
4034 * The
first capture argument of the
lambda function, `data_out_copy`
4035 * in essence creates a local variable
inside the
lambda function into
4036 * which we
"move" the `data_out` variable from above. The way
this works
4037 * is that we create a `std::unique_ptr` above that points to the
DataOut
4038 *
object. But we have no use
for it any more after
this point: We really
4039 * just want to move ownership from the current function to the
lambda
4040 * function implemented in the following few lines. We could have done
4041 *
this by
using a `std::shared_ptr` above and giving the
lambda function
4042 * a
copy of that shared pointer; once the current function returns (but
4043 * maybe
while the lambda function is still running), our local shared
4044 * pointer would go out of scope and stop pointing at the actual
4045 * object, at which
point the
lambda function simply becomes the sole
4046 * owner. But
using the `std::unique_ptr` is conceptually cleaner as it
4047 * makes it clear that the current function
's `data_out` variable isn't
4048 * even pointing to the
object any more.
4051 *
auto output_worker =
4052 * [data_out_copy = std::move(data_out),
this, name, t, cycle]() {
4055 * data_out_copy->set_flags(flags);
4057 * data_out_copy->write_vtu_with_pvtu_record(
4058 *
"", name +
"-solution", cycle, mpi_communicator, 6);
4063 * If the asynchronous writeback option is set we launch a
new
4064 * background thread with the help of
4066 * href=
"https://en.cppreference.com/w/cpp/thread/async"><code>std::async</code></a>
4067 * function. The function returns a <a
4068 * href=
"https://en.cppreference.com/w/cpp/thread/future"><code>std::future</code></a>
4069 *
object that we can use to query the status of the background thread.
4070 * At
this point we can
return from the <code>output()</code> function
4071 * and resume with the time stepping in the main
loop - the thread will
4072 *
run in the background.
4075 *
if (asynchronous_writeback)
4077 * background_thread_state =
4078 * std::async(std::launch::async, std::move(output_worker));
4090 * And
finally, the main function.
4096 *
int main(
int argc,
char *argv[])
4100 *
constexpr int dim = 2;
4102 *
using namespace dealii;
4103 *
using namespace Step69;
4107 *
MPI_Comm mpi_communicator(MPI_COMM_WORLD);
4108 * MainLoop<dim> main_loop(mpi_communicator);
4112 *
catch (std::exception &exc)
4114 * std::cerr << std::endl
4116 * <<
"----------------------------------------------------"
4118 * std::cerr <<
"Exception on processing: " << std::endl
4119 * << exc.what() << std::endl
4120 * <<
"Aborting!" << std::endl
4121 * <<
"----------------------------------------------------"
4127 * std::cerr << std::endl
4129 * <<
"----------------------------------------------------"
4131 * std::cerr <<
"Unknown exception!" << std::endl
4132 * <<
"Aborting!" << std::endl
4133 * <<
"----------------------------------------------------"
4139<a name=
"step_69-Results"></a><h1>Results</h1>
4142Running the program with
default parameters in release mode takes about 1
4143minute on a 4 core machine (with hyperthreading):
4145# mpirun -np 4 ./step-69 | tee output
4146Reading parameters and allocating objects... done
4148 ####################################################
4152 ####################################################
4154Number of active cells: 36864
4156 ####################################################
4158 ######### compute offline
data #########
4160 ####################################################
4162Number of degrees of freedom: 37376
4164 ####################################################
4166 ######### set up time step #########
4168 ####################################################
4170 ####################################################
4175 ####################################################
4177TimeLoop<dim>::interpolate_initial_values(t = 0)
4178TimeLoop<dim>::output(t = 0, checkpoint = 0)
4180 ####################################################
4182 ######### enter main
loop #########
4185 ####################################################
4187 ####################################################
4189 ######### Cycle 000001 (0.0%) #########
4190 ######### at time t = 0.00000000 #########
4192 ####################################################
4196 ####################################################
4198 ######### Cycle 007553 (100.0%) #########
4199 ######### at time t = 3.99984036 #########
4201 ####################################################
4203TimeLoop<dim>::output(t = 4.00038, checkpoint = 1)
4205+------------------------------------------------------------------------+------------+------------+
4206| Total CPU time elapsed since start | 357s | |
4208| Section | no. calls | CPU time | % of total |
4209+------------------------------------------------------------+-----------+------------+------------+
4210| discretization - setup | 1 | 0.113s | 0% |
4211| offline_data -
assemble lumped mass
matrix, and c_ij | 1 | 0.167s | 0% |
4212| offline_data - compute |c_ij|, and n_ij | 1 | 0.00255s | 0% |
4213| offline_data - create sparsity pattern and set up matrices | 1 | 0.0224s | 0% |
4214| offline_data - distribute dofs | 1 | 0.0617s | 0% |
4215| offline_data - fix slip boundary c_ij | 1 | 0.0329s | 0% |
4216| schlieren_postprocessor - compute schlieren plot | 201 | 0.811s | 0.23% |
4217| schlieren_postprocessor - prepare scratch space | 1 | 7.6e-05s | 0% |
4218| time_loop - setup scratch space | 1 | 0.127s | 0% |
4219| time_loop - stalled output | 200 | 0.000685s | 0% |
4220| time_step - 1 compute d_ij | 7553 | 240s | 67% |
4221| time_step - 2 compute d_ii, and tau_max | 7553 | 11.5s | 3.2% |
4222| time_step - 3 perform update | 7553 | 101s | 28% |
4223| time_step - 4 fix boundary states | 7553 | 0.724s | 0.2% |
4224| time_step - prepare scratch space | 1 | 0.00245s | 0% |
4225+------------------------------------------------------------+-----------+------------+------------+
4228One thing that becomes evident is the fact that the program spends two
4229thirds of the execution time computing the graph viscosity d_ij and about a
4230third of the execution time in performing the update, where computing the
4231flux @f$f(U)@f$ is the expensive operation. The preset default resolution is
4232about 37k gridpoints, which amounts to about 148k spatial degrees of
4233freedom in 2D. An animated schlieren plot of the solution looks as follows:
4235<img src=
"https://www.dealii.org/images/steps/developer/step-69.coarse.gif" alt=
"" height=
"300">
4237It is evident that 37k gridpoints for the
first-order method is nowhere
4238near the resolution needed to resolve any flow features. For comparison,
4239here is a
"reference" computation with a
second-order method and about 9.5M
4240gridpoints (<a href=
"https://github.com/conservation-laws/ryujin">github
4243<img src=
"https://www.dealii.org/images/steps/developer/step-69.2nd-order.t400.jpg" alt=
"" height=
"300">
4245So, we give the
first-order method a
second chance and
run it with about
42462.4M gridpoints on a small compute server:
4249# mpirun -np 16 ./step-69 | tee output
4253 ####################################################
4255 ######### Cycle 070216 (100.0%) #########
4256 ######### at time t = 3.99999231 #########
4258 ####################################################
4260TimeLoop<dim>::output(t = 4.00006, checkpoint = 1)
4264+------------------------------------------------------------------------+------------+------------+
4265| Total wallclock time elapsed since start | 6.75e+03s | |
4267| Section | no. calls | wall time | % of total |
4268+------------------------------------------------------------+-----------+------------+------------+
4269| discretization - setup | 1 | 1.97s | 0% |
4270| offline_data -
assemble lumped mass
matrix, and c_ij | 1 | 1.19s | 0% |
4271| offline_data - compute |c_ij|, and n_ij | 1 | 0.0172s | 0% |
4272| offline_data - create sparsity pattern and set up matrices | 1 | 0.413s | 0% |
4273| offline_data - distribute dofs | 1 | 1.05s | 0% |
4274| offline_data - fix slip boundary c_ij | 1 | 0.252s | 0% |
4275| schlieren_postprocessor - compute schlieren plot | 201 | 1.82s | 0% |
4276| schlieren_postprocessor - prepare scratch space | 1 | 0.000497s | 0% |
4277| time_loop - setup scratch space | 1 | 1.45s | 0% |
4278| time_loop - stalled output | 200 | 0.00342s | 0% |
4279| time_step - 1 compute d_ij | 70216 | 4.38e+03s | 65% |
4280| time_step - 2 compute d_ii, and tau_max | 70216 | 419s | 6.2% |
4281| time_step - 3 perform update | 70216 | 1.87e+03s | 28% |
4282| time_step - 4 fix boundary states | 70216 | 24s | 0.36% |
4283| time_step - prepare scratch space | 1 | 0.0227s | 0% |
4284+------------------------------------------------------------+-----------+------------+------------+
4287And with the following result:
4289<img src=
"https://www.dealii.org/images/steps/developer/step-69.fine.gif" alt=
"" height=
"300">
4291That
's substantially better, although of course at the price of having run
4292the code for roughly 2 hours on 16 cores.
4296<a name="step-69-extensions"></a>
4297<a name="step_69-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
4300The program showcased here is really only first-order accurate, as
4301discussed above. The pictures above illustrate how much diffusion that
4302introduces and how far the solution is from one that actually resolves
4303the features we care about.
4305This can be fixed, but it would exceed what a *tutorial* is about.
4306Nevertheless, it is worth showing what one can achieve by adding a
4307second-order scheme. For example, here is a video computed with <a
4308href=https://conservation-laws.org/>the following research code</a>
4309that shows (with a different color scheme) a 2d simulation that corresponds
4310to the cases shown above:
4314 <iframe width="560" height="315" src="https://www.youtube.com/embed/xIwJZlsXpZ4"
4316 allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture"
4317 allowfullscreen></iframe>
4321This simulation was done with 38 million degrees of freedom
4322(continuous @f$Q_1@f$ finite elements) per component of the solution
4323vector. The exquisite detail of the solution is remarkable for these
4324kinds of simulations, including in the sub-sonic region behind the
4327One can also with relative ease further extend this to the 3d case:
4331 <iframe width="560" height="315" src="https://www.youtube.com/embed/vBCRAF_c8m8"
4333 allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture"
4334 allowfullscreen></iframe>
4338Solving this becomes expensive, however: The simulation was done with
43391,817 million degrees of freedom (continuous @f$Q_1@f$ finite elements)
4340per component (for a total of 9.09 billion spatial degrees of freedom)
4341and ran on 30,720 MPI ranks. The code achieved an average throughput of
4342969M grid points per second (0.04M gridpoints per second per CPU). The
4343front and back wall show a "Schlieren plot": the magnitude of the
4344gradient of the density on an exponential scale from white (low) to
4345black (high). All other cutplanes and the surface of the obstacle show
4346the magnitude of the vorticity on a white (low) - yellow (medium) -
4347red (high) scale. (The scales of the individual cutplanes have been
4348adjusted for a nicer visualization.)
4351<a name="step_69-PlainProg"></a>
4352<h1> The plain program</h1>
4353@include "step-69.cc"
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
virtual void build_patches(const unsigned int n_subdivisions=0)
static void initialize(const std::string &filename="", const std::string &output_filename="", const ParameterHandler::OutputStyle output_style_for_output_filename=ParameterHandler::Short, ParameterHandler &prm=ParameterAcceptor::prm, const ParameterHandler::OutputStyle output_style_for_filename=ParameterHandler::DefaultStyle)
boost::signals2::signal< void()> parse_parameters_call_back
void add_parameter(const std::string &entry, ParameterType ¶meter, const std::string &documentation="", ParameterHandler &prm_=prm, const Patterns::PatternBase &pattern= *Patterns::Tools::Convert< ParameterType >::to_pattern())
#define DEAL_II_ALWAYS_INLINE
#define Assert(cond, exc)
#define AssertThrow(cond, exc)
typename ActiveSelector::cell_iterator cell_iterator
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
IteratorRange< BaseIterator > make_iterator_range(const BaseIterator &begin, const std_cxx20::type_identity_t< BaseIterator > &end)
std::vector< index_type > data
@ valid
Iterator points to a valid object.
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
@ diagonal
Matrix is diagonal.
@ general
No special properties.
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > E(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
T max(const T &t, const MPI_Comm mpi_communicator)
T min(const T &t, const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
T scatter(const MPI_Comm comm, const std::vector< T > &objects_to_send, const unsigned int root_process=0)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
void run(const std::vector< std::vector< Iterator > > &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
long double gamma(const unsigned int n)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
void apply_to_subranges(const tbb::blocked_range< Iterator > &range, const Function &f)
void apply_to_subranges(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, const Function &f, const unsigned int grainsize)
boost::integer_range< IncrementableType > iota_view
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
void swap(ObserverPointer< T, P > &t1, ObserverPointer< T, Q > &t2)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
void advance(std::tuple< I1, I2 > &t, const unsigned int n)
void gather(VectorizedArray< Number, width > &out, const std::array< Number *, width > &ptrs, const unsigned int offset)
1.9.8