Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
step-69.h
Go to the documentation of this file.
1 
1462  * cell->get_dof_indices(dof_indices);
1463  * std::transform(dof_indices.begin(),
1464  * dof_indices.end(),
1465  * dof_indices.begin(),
1466  * [&](types::global_dof_index index) {
1467  * return partitioner->global_to_local(index);
1468  * });
1469  *
1470  * /* And simply add, for each dof, a coupling to all other "local"
1471  * * dofs on the cell: */
1472  * for (const auto dof : dof_indices)
1473  * dsp.add_entries(dof, dof_indices.begin(), dof_indices.end());
1474  * }
1475  *
1476  * sparsity_pattern.copy_from(dsp);
1477  *
1478  * lumped_mass_matrix.reinit(sparsity_pattern);
1479  * norm_matrix.reinit(sparsity_pattern);
1480  * for (auto &matrix : cij_matrix)
1481  * matrix.reinit(sparsity_pattern);
1482  * for (auto &matrix : nij_matrix)
1483  * matrix.reinit(sparsity_pattern);
1484  * }
1485  * }
1486  *
1487  * @endcode
1488  *
1489  * This concludes the setup of the DoFHandler and SparseMatrix objects.
1490  * Next, we have to assemble various matrices. We define a number of
1491  * helper functions and data structures in an anonymous namespace.
1492  *
1493 
1494  *
1495  *
1496  * @code
1497  * namespace
1498  * {
1499  * @endcode
1500  *
1501  * <code>CopyData</code> class that will be used to assemble the
1502  * offline data matrices using WorkStream. It acts as a container: it
1503  * is just a struct where WorkStream stores the local cell
1504  * contributions. Note that it also contains a class member
1505  * <code>local_boundary_normal_map</code> used to store the local
1506  * contributions required to compute the normals at the boundary.
1507  *
1508 
1509  *
1510  *
1511  * @code
1512  * template <int dim>
1513  * struct CopyData
1514  * {
1515  * bool is_artificial;
1516  * std::vector<types::global_dof_index> local_dof_indices;
1517  * typename OfflineData<dim>::BoundaryNormalMap local_boundary_normal_map;
1518  * FullMatrix<double> cell_lumped_mass_matrix;
1519  * std::array<FullMatrix<double>, dim> cell_cij_matrix;
1520  * };
1521  *
1522  * @endcode
1523  *
1524  * Next we introduce a number of helper functions that are all
1525  * concerned about reading and writing matrix and vector entries. They
1526  * are mainly motivated by providing slightly more efficient code and
1527  * <a href="https://en.wikipedia.org/wiki/Syntactic_sugar"> syntactic
1528  * sugar</a> for otherwise somewhat tedious code.
1529  *
1530 
1531  *
1532  * The first function we introduce, <code>get_entry()</code>, will be
1533  * used to read the value stored at the entry pointed by a
1534  * SparsityPattern iterator <code>it</code> of <code>matrix</code>. The
1535  * function works around a small deficiency in the SparseMatrix
1536  * interface: The SparsityPattern is concerned with all index
1537  * operations of the sparse matrix stored in CRS format. As such the
1538  * iterator already knows the global index of the corresponding matrix
1539  * entry in the low-level vector stored in the SparseMatrix object. Due
1540  * to the lack of an interface in the SparseMatrix for accessing the
1541  * element directly with a SparsityPattern iterator, we unfortunately
1542  * have to create a temporary SparseMatrix iterator. We simply hide
1543  * this in the <code>get_entry()</code> function.
1544  *
1545 
1546  *
1547  *
1548  * @code
1549  * template <typename IteratorType>
1551  * get_entry(const SparseMatrix<double> &matrix, const IteratorType &it)
1552  * {
1553  * const SparseMatrix<double>::const_iterator matrix_iterator(
1554  * &matrix, it->global_index());
1555  * return matrix_iterator->value();
1556  * }
1557  *
1558  * @endcode
1559  *
1560  * The <code>set_entry()</code> helper is the inverse operation of
1561  * <code>get_value()</code>: Given an iterator and a value, it sets the
1562  * entry pointed to by the iterator in the matrix.
1563  *
1564 
1565  *
1566  *
1567  * @code
1568  * template <typename IteratorType>
1569  * DEAL_II_ALWAYS_INLINE inline void
1570  * set_entry(SparseMatrix<double> & matrix,
1571  * const IteratorType & it,
1573  * {
1574  * SparseMatrix<double>::iterator matrix_iterator(&matrix,
1575  * it->global_index());
1576  * matrix_iterator->value() = value;
1577  * }
1578  *
1579  * @endcode
1580  *
1581  * <code>gather_get_entry()</code>: we note that @f$\mathbf{c}_{ij} \in
1582  * \mathbb{R}^d@f$. If @f$d=2@f$ then @f$\mathbf{c}_{ij} =
1583  * [\mathbf{c}_{ij}^1,\mathbf{c}_{ij}^2]^\top@f$. Which basically implies
1584  * that we need one matrix per space dimension to store the
1585  * @f$\mathbf{c}_{ij}@f$ vectors. Similar observation follows for the
1586  * matrix @f$\mathbf{n}_{ij}@f$. The purpose of
1587  * <code>gather_get_entry()</code> is to retrieve those entries and store
1588  * them into a <code>Tensor<1, dim></code> for our convenience.
1589  *
1590 
1591  *
1592  *
1593  * @code
1594  * template <std::size_t k, typename IteratorType>
1596  * gather_get_entry(const std::array<SparseMatrix<double>, k> &c_ij,
1597  * const IteratorType it)
1598  * {
1599  * Tensor<1, k> result;
1600  * for (unsigned int j = 0; j < k; ++j)
1601  * result[j] = get_entry(c_ij[j], it);
1602  * return result;
1603  * }
1604  *
1605  * @endcode
1606  *
1607  * <code>gather()</code> (first interface): this first function
1608  * signature, having three input arguments, will be used to retrieve
1609  * the individual components <code>(i,l)</code> of a matrix. The
1610  * functionality of <code>gather_get_entry()</code> and
1611  * <code>gather()</code> is very much the same, but their context is
1612  * different: the function <code>gather()</code> does not rely on an
1613  * iterator (that actually knows the value pointed to) but rather on the
1614  * indices <code>(i,l)</code> of the entry in order to retrieve its
1615  * actual value. We should expect <code>gather()</code> to be slightly
1616  * more expensive than <code>gather_get_entry()</code>. The use of
1617  * <code>gather()</code> will be limited to the task of computing the
1618  * algebraic viscosity @f$d_{ij}@f$ in the particular case that when
1619  * both @f$i@f$ and @f$j@f$ lie at the boundary.
1620  *
1621 
1622  *
1623  * @note The reader should be aware that accessing an arbitrary
1624  * <code>(i,l)</code> entry of a matrix (say for instance Trilinos or PETSc
1625  * matrices) is in general unacceptably expensive. Here is where we might
1626  * want to keep an eye on complexity: we want this operation to have
1627  * constant complexity, which is the case of the current implementation
1628  * using deal.II matrices.
1629  *
1630 
1631  *
1632  *
1633  * @code
1634  * template <std::size_t k>
1636  * gather(const std::array<SparseMatrix<double>, k> &n_ij,
1637  * const unsigned int i,
1638  * const unsigned int j)
1639  * {
1640  * Tensor<1, k> result;
1641  * for (unsigned int l = 0; l < k; ++l)
1642  * result[l] = n_ij[l](i, j);
1643  * return result;
1644  * }
1645  *
1646  * @endcode
1647  *
1648  * <code>gather()</code> (second interface): this second function
1649  * signature having two input arguments will be used to gather the
1650  * state at a node <code>i</code> and return it as a
1651  * <code>Tensor<1,problem_dimension></code> for our convenience.
1652  *
1653 
1654  *
1655  *
1656  * @code
1657  * template <std::size_t k>
1659  * gather(const std::array<LinearAlgebra::distributed::Vector<double>, k> &U,
1660  * const unsigned int i)
1661  * {
1662  * Tensor<1, k> result;
1663  * for (unsigned int j = 0; j < k; ++j)
1664  * result[j] = U[j].local_element(i);
1665  * return result;
1666  * }
1667  *
1668  * @endcode
1669  *
1670  * <code>scatter()</code>: this function has three input arguments, the
1671  * first one is meant to be a "global object" (say a locally owned or
1672  * locally relevant vector), the second argument which could be a
1673  * <code>Tensor<1,problem_dimension></code>, and the last argument
1674  * which represents a index of the global object. This function will be
1675  * primarily used to write the updated nodal values, stored as
1676  * <code>Tensor<1,problem_dimension></code>, into the global objects.
1677  *
1678 
1679  *
1680  *
1681  * @code
1682  * template <std::size_t k, int k2>
1683  * DEAL_II_ALWAYS_INLINE inline void
1684  * scatter(std::array<LinearAlgebra::distributed::Vector<double>, k> &U,
1685  * const Tensor<1, k2> & tensor,
1686  * const unsigned int i)
1687  * {
1688  * static_assert(k == k2,
1689  * "The dimensions of the input arguments must agree");
1690  * for (unsigned int j = 0; j < k; ++j)
1691  * U[j].local_element(i) = tensor[j];
1692  * }
1693  * } // namespace
1694  *
1695  * @endcode
1696  *
1697  * We are now in a position to assemble all matrices stored in
1698  * <code>OfflineData</code>: the lumped mass entries @f$m_i@f$, the
1699  * vector-valued matrices @f$\mathbf{c}_{ij}@f$ and @f$\mathbf{n}_{ij} =
1700  * \frac{\mathbf{c}_{ij}}{|\mathbf{c}_{ij}|}@f$, and the boundary normals
1701  * @f$\boldsymbol{\nu}_i@f$.
1702  *
1703 
1704  *
1705  * In order to exploit thread parallelization we use the WorkStream approach
1706  * detailed in the @ref threads "Parallel computing with multiple processors"
1707  * accessing shared memory. As customary this requires
1708  * definition of
1709  * - Scratch data (i.e. input info required to carry out computations): in
1710  * this case it is <code>scratch_data</code>.
1711  * - The worker: in our case this is the <code>local_assemble_system()</code>
1712  * function that
1713  * actually computes the local (i.e. current cell) contributions from the
1714  * scratch data.
1715  * - A copy data: a struct that contains all the local assembly
1716  * contributions, in this case <code>CopyData<dim>()</code>.
1717  * - A copy data routine: in this case it is
1718  * <code>copy_local_to_global()</code> in charge of actually coping these
1719  * local contributions into the global objects (matrices and/or vectors)
1720  *
1721 
1722  *
1723  * Most of the following lines are spent in the definition of the worker
1724  * <code>local_assemble_system()</code> and the copy data routine
1725  * <code>copy_local_to_global()</code>. There is not much to say about the
1726  * WorkStream framework since the vast majority of ideas are reasonably
1727  * well-documented in @ref step_9 "step-9", @ref step_13 "step-13" and @ref step_32 "step-32" among others.
1728  *
1729 
1730  *
1731  * Finally, assuming that @f$\mathbf{x}_i@f$ is a support point at the boundary,
1732  * the (nodal) normals are defined as:
1733  *
1734 
1735  *
1736  * @f{align*}
1737  * \widehat{\boldsymbol{\nu}}_i \dealcoloneq
1738  * \frac{\boldsymbol{\nu}_i}{|\boldsymbol{\nu}_i|} \ \text{ where }
1739  * \boldsymbol{\nu}_i \dealcoloneq \sum_{T \subset \text{supp}(\phi_i)}
1740  * \sum_{F \subset \partial T \cap \partial \Omega}
1741  * \sum_{\mathbf{x}_{q,F}} \nu(\mathbf{x}_{q,F})
1742  * \phi_i(\mathbf{x}_{q,F}).
1743  * @f}
1744  *
1745 
1746  *
1747  * Here @f$T@f$ denotes elements,
1748  * @f$\text{supp}(\phi_i)@f$ the support of the shape function @f$\phi_i@f$,
1749  * @f$F@f$ are faces of the element @f$T@f$, and @f$\mathbf{x}_{q,F}@f$
1750  * are quadrature points on such face. Note that this formula for
1751  * @f$\widehat{\boldsymbol{\nu}}_i@f$ is nothing else than some form of
1752  * weighted averaging. Other more sophisticated definitions for @f$\nu_i@f$
1753  * are possible but none of them have much influence in theory or practice.
1754  *
1755 
1756  *
1757  *
1758  * @code
1759  * template <int dim>
1761  * {
1762  * lumped_mass_matrix = 0.;
1763  * norm_matrix = 0.;
1764  * for (auto &matrix : cij_matrix)
1765  * matrix = 0.;
1766  * for (auto &matrix : nij_matrix)
1767  * matrix = 0.;
1768  *
1769  * unsigned int dofs_per_cell = discretization->finite_element.dofs_per_cell;
1770  * unsigned int n_q_points = discretization->quadrature.size();
1771  *
1772  * @endcode
1773  *
1774  * What follows is the initialization of the scratch data required by
1775  * WorkStream
1776  *
1777 
1778  *
1779  *
1780  * @code
1781  * MeshWorker::ScratchData<dim> scratch_data(
1782  * discretization->mapping,
1783  * discretization->finite_element,
1784  * discretization->quadrature,
1787  * discretization->face_quadrature,
1789  *
1790  * {
1791  * TimerOutput::Scope scope(
1792  * computing_timer,
1793  * "offline_data - assemble lumped mass matrix, and c_ij");
1794  *
1795  * const auto local_assemble_system =
1796  * [&](const typename DoFHandler<dim>::cell_iterator &cell,
1797  * MeshWorker::ScratchData<dim> & scratch,
1798  * CopyData<dim> & copy) {
1799  * copy.is_artificial = cell->is_artificial();
1800  * if (copy.is_artificial)
1801  * return;
1802  *
1803  * copy.local_boundary_normal_map.clear();
1804  * copy.cell_lumped_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
1805  * for (auto &matrix : copy.cell_cij_matrix)
1806  * matrix.reinit(dofs_per_cell, dofs_per_cell);
1807  *
1808  * const auto &fe_values = scratch.reinit(cell);
1809  *
1810  * copy.local_dof_indices.resize(dofs_per_cell);
1811  * cell->get_dof_indices(copy.local_dof_indices);
1812  *
1813  * std::transform(copy.local_dof_indices.begin(),
1814  * copy.local_dof_indices.end(),
1815  * copy.local_dof_indices.begin(),
1816  * [&](types::global_dof_index index) {
1817  * return partitioner->global_to_local(index);
1818  * });
1819  *
1820  * @endcode
1821  *
1822  * We compute the local contributions for the lumped mass matrix
1823  * entries @f$m_i@f$ and and vectors @f$c_{ij}@f$ in the usual fashion:
1824  *
1825  * @code
1826  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1827  * {
1828  * const auto JxW = fe_values.JxW(q_point);
1829  *
1830  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1831  * {
1832  * const auto value_JxW =
1833  * fe_values.shape_value(j, q_point) * JxW;
1834  * const auto grad_JxW = fe_values.shape_grad(j, q_point) * JxW;
1835  *
1836  * copy.cell_lumped_mass_matrix(j, j) += value_JxW;
1837  *
1838  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1839  * {
1840  * const auto value = fe_values.shape_value(i, q_point);
1841  * for (unsigned int d = 0; d < dim; ++d)
1842  * copy.cell_cij_matrix[d](i, j) += value * grad_JxW[d];
1843  *
1844  * } /* i */
1845  * } /* j */
1846  * } /* q */
1847  *
1848  * @endcode
1849  *
1850  * Now we have to compute the boundary normals. Note that the
1851  * following loop does not do much unless the element has faces on
1852  * the boundary of the domain.
1853  *
1854  * @code
1855  * for (auto f : GeometryInfo<dim>::face_indices())
1856  * {
1857  * const auto face = cell->face(f);
1858  * const auto id = face->boundary_id();
1859  *
1860  * if (!face->at_boundary())
1861  * continue;
1862  *
1863  * const auto &fe_face_values = scratch.reinit(cell, f);
1864  *
1865  * const unsigned int n_face_q_points =
1866  * fe_face_values.get_quadrature().size();
1867  *
1868  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1869  * {
1870  * if (!discretization->finite_element.has_support_on_face(j, f))
1871  * continue;
1872  *
1873  * @endcode
1874  *
1875  * Note that "normal" will only represent the contributions
1876  * from one of the faces in the support of the shape
1877  * function phi_j. So we cannot normalize this local
1878  * contribution right here, we have to take it "as is",
1879  * store it and pass it to the copy data routine. The
1880  * proper normalization requires an additional loop on
1881  * nodes. This is done in the copy function below.
1882  *
1883  * @code
1884  * Tensor<1, dim> normal;
1885  * if (id == Boundaries::free_slip)
1886  * {
1887  * for (unsigned int q = 0; q < n_face_q_points; ++q)
1888  * normal += fe_face_values.normal_vector(q) *
1889  * fe_face_values.shape_value(j, q);
1890  * }
1891  *
1892  * const auto index = copy.local_dof_indices[j];
1893  *
1894  * Point<dim> position;
1895  * for (unsigned int v = 0;
1896  * v < GeometryInfo<dim>::vertices_per_cell;
1897  * ++v)
1898  * if (cell->vertex_dof_index(v, 0) ==
1899  * partitioner->local_to_global(index))
1900  * {
1901  * position = cell->vertex(v);
1902  * break;
1903  * }
1904  *
1905  * const auto old_id =
1906  * std::get<1>(copy.local_boundary_normal_map[index]);
1907  * copy.local_boundary_normal_map[index] =
1908  * std::make_tuple(normal, std::max(old_id, id), position);
1909  * }
1910  * }
1911  * };
1912  *
1913  * @endcode
1914  *
1915  * Last, we provide a copy_local_to_global function as required for
1916  * the WorkStream
1917  *
1918  * @code
1919  * const auto copy_local_to_global = [&](const CopyData<dim> &copy) {
1920  * if (copy.is_artificial)
1921  * return;
1922  *
1923  * for (const auto &it : copy.local_boundary_normal_map)
1924  * {
1925  * std::get<0>(boundary_normal_map[it.first]) +=
1926  * std::get<0>(it.second);
1927  * std::get<1>(boundary_normal_map[it.first]) =
1928  * std::max(std::get<1>(boundary_normal_map[it.first]),
1929  * std::get<1>(it.second));
1930  * std::get<2>(boundary_normal_map[it.first]) = std::get<2>(it.second);
1931  * }
1932  *
1933  * lumped_mass_matrix.add(copy.local_dof_indices,
1934  * copy.cell_lumped_mass_matrix);
1935  *
1936  * for (int k = 0; k < dim; ++k)
1937  * {
1938  * cij_matrix[k].add(copy.local_dof_indices, copy.cell_cij_matrix[k]);
1939  * nij_matrix[k].add(copy.local_dof_indices, copy.cell_cij_matrix[k]);
1940  * }
1941  * };
1942  *
1943  * WorkStream::run(dof_handler.begin_active(),
1944  * dof_handler.end(),
1945  * local_assemble_system,
1946  * copy_local_to_global,
1947  * scratch_data,
1948  * CopyData<dim>());
1949  * }
1950  *
1951  * @endcode
1952  *
1953  * At this point in time we are done with the computation of @f$m_i@f$ and
1954  * @f$\mathbf{c}_{ij}@f$, but so far the matrix <code>nij_matrix</code>
1955  * contains just a copy of the matrix <code>cij_matrix</code>.
1956  * That's not what we really want: we have to normalize its entries. In
1957  * addition, we have not filled the entries of the matrix
1958  * <code>norm_matrix</code> and the vectors stored in the map
1959  * <code>OfflineData<dim>::BoundaryNormalMap</code> are not normalized.
1960  *
1961 
1962  *
1963  * In principle, this is just offline data, it doesn't make much sense
1964  * to over-optimize their computation, since their cost will get amortized
1965  * over the many time steps that we are going to use. However,
1966  * computing/storing the entries of the matrix
1967  * <code>norm_matrix</code> and the normalization of <code>nij_matrix</code>
1968  * are perfect to illustrate thread-parallel node-loops:
1969  * - we want to visit every node @f$i@f$ in the mesh/sparsity graph,
1970  * - and for every such node we want to visit to every @f$j@f$ such that
1971  * @f$\mathbf{c}_{ij} \not \equiv 0@f$.
1972  *
1973 
1974  *
1975  * From an algebraic point of view, this is equivalent to: visiting
1976  * every row in the matrix and for each one of these rows execute a loop on
1977  * the columns. Node-loops is a core theme of this tutorial step (see
1978  * the pseudo-code in the introduction) that will repeat over and over
1979  * again. That's why this is the right time to introduce them.
1980  *
1981 
1982  *
1983  * We have the thread parallelization capability
1984  * parallel::apply_to_subranges() that is somehow more general than the
1985  * WorkStream framework. In particular, parallel::apply_to_subranges() can
1986  * be used for our node-loops. This functionality requires four input
1987  * arguments which we explain in detail (for the specific case of our
1988  * thread-parallel node loops):
1989  * - The iterator <code>indices.begin()</code> points to a row index.
1990  * - The iterator <code>indices.end()</code> points to a numerically higher
1991  * row index.
1992  * - The function <code>on_subranges(i1,i2)</code> (where <code>i1</code>
1993  * and <code>i2</code> define a sub-range within the range spanned by
1994  * the end and begin iterators defined in the two previous bullets)
1995  * applies an operation to every iterator in such subrange. We may as
1996  * well call <code>on_subranges</code> the "worker".
1997  * - Grainsize: minimum number of iterators (in this case representing
1998  * rows) processed by each thread. We decided for a minimum of 4096
1999  * rows.
2000  *
2001 
2002  *
2003  * A minor caveat here is that the iterators <code>indices.begin()</code>
2004  * and <code>indices.end()</code> supplied to
2005  * parallel::apply_to_subranges() have to be random access iterators:
2006  * internally, parallel::apply_to_subranges() will break the range
2007  * defined by the <code>indices.begin()</code> and
2008  * <code>indices.end()</code> iterators into subranges (we want to be
2009  * able to read any entry in those subranges with constant complexity).
2010  * In order to provide such iterators we resort to
2011  * std_cxx20::ranges::iota_view.
2012  *
2013 
2014  *
2015  * The bulk of the following piece of code is spent defining
2016  * the "worker" <code>on_subranges</code>: i.e. the operation applied at
2017  * each row of the sub-range. Given a fixed <code>row_index</code>
2018  * we want to visit every column/entry in such row. In order to execute
2019  * such columns-loops we use
2020  * <a href="http://www.cplusplus.com/reference/algorithm/for_each/">
2021  * std::for_each</a>
2022  * from the standard library, where:
2023  * - <code>sparsity_pattern.begin(row_index)</code>
2024  * gives us an iterator starting at the first column of the row,
2025  * - <code>sparsity_pattern.end(row_index)</code> is an iterator pointing
2026  * at the last column of the row,
2027  * - the last argument required by `std::for_each` is the operation
2028  * applied at each nonzero entry (a lambda expression in this case)
2029  * of such row.
2030  *
2031 
2032  *
2033  * We note that, parallel::apply_to_subranges() will operate on
2034  * disjoint sets of rows (the subranges) and our goal is to write into
2035  * these rows. Because of the simple nature of the operations we want
2036  * to carry out (computation and storage of normals, and normalization
2037  * of the @f$\mathbf{c}_{ij}@f$ of entries) threads cannot conflict
2038  * attempting to write the same entry (we do not need a scheduler).
2039  *
2040  * @code
2041  * {
2042  * TimerOutput::Scope scope(computing_timer,
2043  * "offline_data - compute |c_ij|, and n_ij");
2044  *
2045  * const std_cxx20::ranges::iota_view<unsigned int, unsigned int> indices(
2046  * 0, n_locally_relevant);
2047  *
2048  * const auto on_subranges =
2049  * [&](
2050  * std_cxx20::ranges::iota_view<unsigned int, unsigned int>::iterator i1,
2051  * const std_cxx20::ranges::iota_view<unsigned int,
2052  * unsigned int>::iterator i2) {
2053  * for (const auto row_index :
2054  * std_cxx20::ranges::iota_view<unsigned int, unsigned int>(*i1,
2055  * *i2))
2056  * {
2057  * @endcode
2058  *
2059  * First column-loop: we compute and store the entries of the
2060  * matrix norm_matrix and write normalized entries into the
2061  * matrix nij_matrix:
2062  *
2063  * @code
2064  * std::for_each(
2065  * sparsity_pattern.begin(row_index),
2066  * sparsity_pattern.end(row_index),
2067  * [&](const ::SparsityPatternIterators::Accessor &jt) {
2068  * const auto c_ij = gather_get_entry(cij_matrix, &jt);
2069  * const double norm = c_ij.norm();
2070  *
2071  * set_entry(norm_matrix, &jt, norm);
2072  * for (unsigned int j = 0; j < dim; ++j)
2073  * set_entry(nij_matrix[j], &jt, c_ij[j] / norm);
2074  * });
2075  * }
2076  * };
2077  *
2078  * parallel::apply_to_subranges(indices.begin(),
2079  * indices.end(),
2080  * on_subranges,
2081  * 4096);
2082  *
2083  * @endcode
2084  *
2085  * Finally, we normalize the vectors stored in
2086  * <code>OfflineData<dim>::BoundaryNormalMap</code>. This operation has
2087  * not been thread parallelized as it would neither illustrate any
2088  * important concept nor lead to any noticeable speed gain.
2089  *
2090  * @code
2091  * for (auto &it : boundary_normal_map)
2092  * {
2093  * auto &normal = std::get<0>(it.second);
2094  * normal /= (normal.norm() + std::numeric_limits<double>::epsilon());
2095  * }
2096  * }
2097  *
2098  * @endcode
2099  *
2100  * As commented in the introduction (see section titled "Stable boundary
2101  * conditions and conservation properties") we use three different types of
2102  * boundary conditions: essential-like boundary conditions (we prescribe a
2103  * state in the left portion of our domain), outflow boundary conditions
2104  * (also called "do-nothing" boundary conditions) in the right portion of
2105  * the domain, and "reflecting" boundary conditions (also called "free
2106  * slip" boundary conditions). With these boundary conditions we should
2107  * not expect any form of conservation to hold.
2108  *
2109 
2110  *
2111  * However, if we were to use reflecting boundary conditions
2112  * @f$\mathbf{m} \cdot \boldsymbol{\nu}_i =0@f$ on the entirety of the
2113  * boundary we should preserve the density and total (mechanical)
2114  * energy. This requires us to modify the @f$\mathbf{c}_{ij}@f$ vectors at
2115  * the boundary as follows @cite GuermondEtAl2018 :
2116  *
2117 
2118  *
2119  * @f{align*}
2120  * \mathbf{c}_{ij} \, +\!\!= \int_{\partial \Omega}
2121  * (\boldsymbol{\nu}_j - \boldsymbol{\nu}(\mathbf{s})) \phi_i \phi_j \,
2122  * \mathrm{d}\mathbf{s} \ \text{whenever} \ \mathbf{x}_i \text{ and }
2123  * \mathbf{x}_j \text{ lie in the boundary.}
2124  * @f}
2125  *
2126 
2127  *
2128  * The ideas repeat themselves: we use WorkStream in order to compute
2129  * this correction, most of the following code is about the definition
2130  * of the worker <code>local_assemble_system()</code>.
2131  *
2132  * @code
2133  * {
2134  * TimerOutput::Scope scope(computing_timer,
2135  * "offline_data - fix slip boundary c_ij");
2136  *
2137  * const auto local_assemble_system =
2138  * [&](const typename DoFHandler<dim>::cell_iterator &cell,
2139  * MeshWorker::ScratchData<dim> & scratch,
2140  * CopyData<dim> & copy) {
2141  * copy.is_artificial = cell->is_artificial();
2142  *
2143  * if (copy.is_artificial)
2144  * return;
2145  *
2146  * for (auto &matrix : copy.cell_cij_matrix)
2147  * matrix.reinit(dofs_per_cell, dofs_per_cell);
2148  *
2149  * copy.local_dof_indices.resize(dofs_per_cell);
2150  * cell->get_dof_indices(copy.local_dof_indices);
2151  * std::transform(copy.local_dof_indices.begin(),
2152  * copy.local_dof_indices.end(),
2153  * copy.local_dof_indices.begin(),
2154  * [&](types::global_dof_index index) {
2155  * return partitioner->global_to_local(index);
2156  * });
2157  *
2158  * for (auto &matrix : copy.cell_cij_matrix)
2159  * matrix = 0.;
2160  *
2161  * for (auto f : GeometryInfo<dim>::face_indices())
2162  * {
2163  * const auto face = cell->face(f);
2164  * const auto id = face->boundary_id();
2165  *
2166  * if (!face->at_boundary())
2167  * continue;
2168  *
2169  * if (id != Boundaries::free_slip)
2170  * continue;
2171  *
2172  * const auto &fe_face_values = scratch.reinit(cell, f);
2173  *
2174  * const unsigned int n_face_q_points =
2175  * fe_face_values.get_quadrature().size();
2176  *
2177  * for (unsigned int q = 0; q < n_face_q_points; ++q)
2178  * {
2179  * const auto JxW = fe_face_values.JxW(q);
2180  * const auto normal_q = fe_face_values.normal_vector(q);
2181  *
2182  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
2183  * {
2184  * if (!discretization->finite_element.has_support_on_face(
2185  * j, f))
2186  * continue;
2187  *
2188  * const auto &normal_j = std::get<0>(
2189  * boundary_normal_map[copy.local_dof_indices[j]]);
2190  *
2191  * const auto value_JxW =
2192  * fe_face_values.shape_value(j, q) * JxW;
2193  *
2194  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2195  * {
2196  * const auto value = fe_face_values.shape_value(i, q);
2197  *
2198  * /* This is the correction of the boundary c_ij */
2199  * for (unsigned int d = 0; d < dim; ++d)
2200  * copy.cell_cij_matrix[d](i, j) +=
2201  * (normal_j[d] - normal_q[d]) * (value * value_JxW);
2202  * } /* i */
2203  * } /* j */
2204  * } /* q */
2205  * } /* f */
2206  * };
2207  *
2208  * const auto copy_local_to_global = [&](const CopyData<dim> &copy) {
2209  * if (copy.is_artificial)
2210  * return;
2211  *
2212  * for (int k = 0; k < dim; ++k)
2213  * cij_matrix[k].add(copy.local_dof_indices, copy.cell_cij_matrix[k]);
2214  * };
2215  *
2216  * WorkStream::run(dof_handler.begin_active(),
2217  * dof_handler.end(),
2218  * local_assemble_system,
2219  * copy_local_to_global,
2220  * scratch_data,
2221  * CopyData<dim>());
2222  * }
2223  * }
2224  *
2225  * @endcode
2226  *
2227  * At this point we are very much done with anything related to offline data.
2228  *
2229 
2230  *
2231  *
2232  * <a name="EquationofstateandapproximateRiemannsolver"></a>
2233  * <h4>Equation of state and approximate Riemann solver</h4>
2234  *
2235 
2236  *
2237  * In this section we describe the implementation of the class members of
2238  * the <code>ProblemDescription</code> class. Most of the code here is
2239  * specific to the compressible Euler's equations with an ideal gas law.
2240  * If we wanted to re-purpose @ref step_69 "step-69" for a different conservation law
2241  * (say for: instance the shallow water equation) most of the
2242  * implementation of this class would have to change. But most of the
2243  * other classes (in particular those defining loop structures) would
2244  * remain unchanged.
2245  *
2246 
2247  *
2248  * We start by implementing a number of small member functions for
2249  * computing <code>momentum</code>, <code>internal_energy</code>,
2250  * <code>pressure</code>, <code>speed_of_sound</code>, and the flux
2251  * <code>f</code> of the system. The functionality of each one of these
2252  * functions is self-explanatory from their names.
2253  *
2254 
2255  *
2256  *
2257  * @code
2258  * template <int dim>
2260  * ProblemDescription<dim>::momentum(const state_type &U)
2261  * {
2262  * Tensor<1, dim> result;
2263  * std::copy(&U[1], &U[1 + dim], &result[0]);
2264  * return result;
2265  * }
2266  *
2267  * template <int dim>
2268  * DEAL_II_ALWAYS_INLINE inline double
2269  * ProblemDescription<dim>::internal_energy(const state_type &U)
2270  * {
2271  * const double &rho = U[0];
2272  * const auto m = momentum(U);
2273  * const double &E = U[dim + 1];
2274  * return E - 0.5 * m.norm_square() / rho;
2275  * }
2276  *
2277  * template <int dim>
2278  * DEAL_II_ALWAYS_INLINE inline double
2279  * ProblemDescription<dim>::pressure(const state_type &U)
2280  * {
2281  * return (gamma - 1.) * internal_energy(U);
2282  * }
2283  *
2284  * template <int dim>
2285  * DEAL_II_ALWAYS_INLINE inline double
2286  * ProblemDescription<dim>::speed_of_sound(const state_type &U)
2287  * {
2288  * const double &rho = U[0];
2289  * const double p = pressure(U);
2290  *
2291  * return std::sqrt(gamma * p / rho);
2292  * }
2293  *
2294  * template <int dim>
2295  * DEAL_II_ALWAYS_INLINE inline typename ProblemDescription<dim>::flux_type
2296  * ProblemDescription<dim>::flux(const state_type &U)
2297  * {
2298  * const double &rho = U[0];
2299  * const auto m = momentum(U);
2300  * const auto p = pressure(U);
2301  * const double &E = U[dim + 1];
2302  *
2303  * flux_type result;
2304  *
2305  * result[0] = m;
2306  * for (unsigned int i = 0; i < dim; ++i)
2307  * {
2308  * result[1 + i] = m * m[i] / rho;
2309  * result[1 + i][i] += p;
2310  * }
2311  * result[dim + 1] = m / rho * (E + p);
2312  *
2313  * return result;
2314  * }
2315  *
2316  * @endcode
2317  *
2318  * Now we discuss the computation of @f$\lambda_{\text{max}}
2319  * (\mathbf{U}_i^{n},\mathbf{U}_j^{n}, \textbf{n}_{ij})@f$. The analysis
2320  * and derivation of sharp upper-bounds of maximum wavespeeds of Riemann
2321  * problems is a very technical endeavor and we cannot include an
2322  * advanced discussion about it in this tutorial. In this portion of the
2323  * documentation we will limit ourselves to sketch the main functionality
2324  * of our implementation functions and point to specific academic
2325  * references in order to help the (interested) reader trace the
2326  * source (and proper mathematical justification) of these ideas.
2327  *
2328 
2329  *
2330  * In general, obtaining a sharp guaranteed upper-bound on the maximum
2331  * wavespeed requires solving a quite expensive scalar nonlinear problem.
2332  * This is typically done with an iterative solver. In order to simplify
2333  * the presentation in this example step we decided not to include such
2334  * an iterative scheme. Instead, we will just use an initial guess as a
2335  * guess for an upper bound on the maximum wavespeed. More precisely,
2336  * equations (2.11) (3.7), (3.8) and (4.3) of @cite GuermondPopov2016b
2337  * are enough to define a guaranteed upper bound on the maximum
2338  * wavespeed. This estimate is returned by a call to the function
2339  * <code>lambda_max_two_rarefaction()</code>. At its core the
2340  * construction of such an upper bound uses the so-called two-rarefaction
2341  * approximation for the intermediate pressure @f$p^*@f$, see for instance
2342  * Equation (4.46), page 128 in @cite Toro2009.
2343  *
2344 
2345  *
2346  * The estimate returned by <code>lambda_max_two_rarefaction()</code> is
2347  * guaranteed to be an upper bound, it is in general quite sharp, and
2348  * overall sufficient for our purposes. However, for some specific
2349  * situations (in particular when one of states is close to vacuum
2350  * conditions) such an estimate will be overly pessimistic. That's why we
2351  * used a second estimate to avoid this degeneracy that will be invoked
2352  * by a call to the function <code>lambda_max_expansion()</code>. The most
2353  * important function here is <code>compute_lambda_max()</code> which
2354  * takes the minimum between the estimates returned by
2355  * <code>lambda_max_two_rarefaction()</code> and
2356  * <code>lambda_max_expansion()</code>.
2357  *
2358 
2359  *
2360  * We start again by defining a couple of helper functions:
2361  *
2362 
2363  *
2364  * The first function takes a state <code>U</code> and a unit vector
2365  * <code>n_ij</code> and computes the <i>projected</i> 1D state in
2366  * direction of the unit vector.
2367  *
2368  * @code
2369  * namespace
2370  * {
2371  * template <int dim>
2372  * DEAL_II_ALWAYS_INLINE inline std::array<double, 4> riemann_data_from_state(
2373  * const typename ProblemDescription<dim>::state_type U,
2374  * const Tensor<1, dim> & n_ij)
2375  * {
2376  * Tensor<1, 3> projected_U;
2377  * projected_U[0] = U[0];
2378  *
2379  * @endcode
2380  *
2381  * For this, we have to change the momentum to @f$\textbf{m}\cdot
2382  * n_{ij}@f$ and have to subtract the kinetic energy of the
2383  * perpendicular part from the total energy:
2384  *
2385  * @code
2386  * const auto m = ProblemDescription<dim>::momentum(U);
2387  * projected_U[1] = n_ij * m;
2388  *
2389  * const auto perpendicular_m = m - projected_U[1] * n_ij;
2390  * projected_U[2] = U[1 + dim] - 0.5 * perpendicular_m.norm_square() / U[0];
2391  *
2392  * @endcode
2393  *
2394  * We return the 1D state in <i>primitive</i> variables instead of
2395  * conserved quantities. The return array consists of density @f$\rho@f$,
2396  * velocity @f$u@f$, pressure @f$p@f$ and local speed of sound @f$a@f$:
2397  *
2398 
2399  *
2400  *
2401  * @code
2402  * return {{projected_U[0],
2403  * projected_U[1] / projected_U[0],
2404  * ProblemDescription<1>::pressure(projected_U),
2405  * ProblemDescription<1>::speed_of_sound(projected_U)}};
2406  * }
2407  *
2408  * @endcode
2409  *
2410  * At this point we also define two small functions that return the
2411  * positive and negative part of a double.
2412  *
2413 
2414  *
2415  *
2416  * @code
2417  * DEAL_II_ALWAYS_INLINE inline double positive_part(const double number)
2418  * {
2419  * return std::max(number, 0.);
2420  * }
2421  *
2422  *
2423  * DEAL_II_ALWAYS_INLINE inline double negative_part(const double number)
2424  * {
2425  * return -std::min(number, 0.);
2426  * }
2427  *
2428  * @endcode
2429  *
2430  * Next, we need two local wavenumbers that are defined in terms of a
2431  * primitive state @f$[\rho, u, p, a]@f$ and a given pressure @f$p^\ast@f$
2432  * @cite GuermondPopov2016 Eqn. (3.7):
2433  * @f{align*}
2434  * \lambda^- = u - a\,\sqrt{1 + \frac{\gamma+1}{2\gamma}
2435  * \left(\frac{p^\ast-p}{p}\right)_+}
2436  * @f}
2437  * Here, the @f$(\cdot)_{+}@f$ denotes the positive part of the given
2438  * argument.
2439  *
2440 
2441  *
2442  *
2443  * @code
2444  * DEAL_II_ALWAYS_INLINE inline double
2445  * lambda1_minus(const std::array<double, 4> &riemann_data,
2446  * const double p_star)
2447  * {
2448  * /* Implements formula (3.7) in Guermond-Popov-2016 */
2449  *
2450  * constexpr double gamma = ProblemDescription<1>::gamma;
2451  * const auto u = riemann_data[1];
2452  * const auto p = riemann_data[2];
2453  * const auto a = riemann_data[3];
2454  *
2455  * const double factor = (gamma + 1.0) / 2.0 / gamma;
2456  * const double tmp = positive_part((p_star - p) / p);
2457  * return u - a * std::sqrt(1.0 + factor * tmp);
2458  * }
2459  *
2460  * @endcode
2461  *
2462  * Analougously @cite GuermondPopov2016 Eqn. (3.8):
2463  * @f{align*}
2464  * \lambda^+ = u + a\,\sqrt{1 + \frac{\gamma+1}{2\gamma}
2465  * \left(\frac{p^\ast-p}{p}\right)_+}
2466  * @f}
2467  *
2468 
2469  *
2470  *
2471  * @code
2472  * DEAL_II_ALWAYS_INLINE inline double
2473  * lambda3_plus(const std::array<double, 4> &riemann_data, const double p_star)
2474  * {
2475  * /* Implements formula (3.8) in Guermond-Popov-2016 */
2476  *
2477  * constexpr double gamma = ProblemDescription<1>::gamma;
2478  * const auto u = riemann_data[1];
2479  * const auto p = riemann_data[2];
2480  * const auto a = riemann_data[3];
2481  *
2482  * const double factor = (gamma + 1.0) / 2.0 / gamma;
2483  * const double tmp = positive_part((p_star - p) / p);
2484  * return u + a * std::sqrt(1.0 + factor * tmp);
2485  * }
2486  *
2487  * @endcode
2488  *
2489  * All that is left to do is to compute the maximum of @f$\lambda^-@f$ and
2490  * @f$\lambda^+@f$ computed from the left and right primitive state
2491  * (@cite GuermondPopov2016 Eqn. (2.11)), where @f$p^\ast@f$ is given by
2492  * @cite GuermondPopov2016 Eqn (4.3):
2493  *
2494 
2495  *
2496  *
2497  * @code
2498  * DEAL_II_ALWAYS_INLINE inline double
2499  * lambda_max_two_rarefaction(const std::array<double, 4> &riemann_data_i,
2500  * const std::array<double, 4> &riemann_data_j)
2501  * {
2502  * constexpr double gamma = ProblemDescription<1>::gamma;
2503  * const auto u_i = riemann_data_i[1];
2504  * const auto p_i = riemann_data_i[2];
2505  * const auto a_i = riemann_data_i[3];
2506  * const auto u_j = riemann_data_j[1];
2507  * const auto p_j = riemann_data_j[2];
2508  * const auto a_j = riemann_data_j[3];
2509  *
2510  * const double numerator = a_i + a_j - (gamma - 1.) / 2. * (u_j - u_i);
2511  *
2512  * const double denominator =
2513  * a_i * std::pow(p_i / p_j, -1. * (gamma - 1.) / 2. / gamma) + a_j * 1.;
2514  *
2515  * /* Formula (4.3) in Guermond-Popov-2016 */
2516  *
2517  * const double p_star =
2518  * p_j * std::pow(numerator / denominator, 2. * gamma / (gamma - 1));
2519  *
2520  * const double lambda1 = lambda1_minus(riemann_data_i, p_star);
2521  * const double lambda3 = lambda3_plus(riemann_data_j, p_star);
2522  *
2523  * /* Formula (2.11) in Guermond-Popov-2016 */
2524  *
2525  * return std::max(positive_part(lambda3), negative_part(lambda1));
2526  * }
2527  *
2528  * @endcode
2529  *
2530  * We compute the second upper bound of the maximal wavespeed that is, in
2531  * general, not as sharp as the two-rarefaction estimate. But it will
2532  * save the day in the context of near vacuum conditions when the
2533  * two-rarefaction approximation might attain extreme values:
2534  * @f{align*}
2535  * \lambda_{\text{exp}} = \max(u_i,u_j) + 5. \max(a_i, a_j).
2536  * @f}
2537  * @note The constant 5.0 multiplying the maximum of the sound speeds
2538  * is <i>neither</i> an ad-hoc constant, <i>nor</i> a tuning parameter.
2539  * It defines an upper bound for any @f$\gamma \in [0,5/3]@f$. Do not play
2540  * with it!
2541  *
2542 
2543  *
2544  *
2545  * @code
2546  * DEAL_II_ALWAYS_INLINE inline double
2547  * lambda_max_expansion(const std::array<double, 4> &riemann_data_i,
2548  * const std::array<double, 4> &riemann_data_j)
2549  * {
2550  * const auto u_i = riemann_data_i[1];
2551  * const auto a_i = riemann_data_i[3];
2552  * const auto u_j = riemann_data_j[1];
2553  * const auto a_j = riemann_data_j[3];
2554  *
2555  * return std::max(std::abs(u_i), std::abs(u_j)) + 5. * std::max(a_i, a_j);
2556  * }
2557  * } // namespace
2558  *
2559  * @endcode
2560  *
2561  * The following is the main function that we are going to call in order to
2562  * compute @f$\lambda_{\text{max}} (\mathbf{U}_i^{n},\mathbf{U}_j^{n},
2563  * \textbf{n}_{ij})@f$. We simply compute both maximal wavespeed estimates
2564  * and return the minimum.
2565  *
2566 
2567  *
2568  *
2569  * @code
2570  * template <int dim>
2571  * DEAL_II_ALWAYS_INLINE inline double
2572  * ProblemDescription<dim>::compute_lambda_max(const state_type & U_i,
2573  * const state_type & U_j,
2574  * const Tensor<1, dim> &n_ij)
2575  * {
2576  * const auto riemann_data_i = riemann_data_from_state(U_i, n_ij);
2577  * const auto riemann_data_j = riemann_data_from_state(U_j, n_ij);
2578  *
2579  * const double lambda_1 =
2580  * lambda_max_two_rarefaction(riemann_data_i, riemann_data_j);
2581  *
2582  * const double lambda_2 =
2583  * lambda_max_expansion(riemann_data_i, riemann_data_j);
2584  *
2585  * return std::min(lambda_1, lambda_2);
2586  * }
2587  *
2588  * @endcode
2589  *
2590  * We conclude this section by defining static arrays
2591  * <code>component_names</code> that contain strings describing the
2592  * component names of our state vector. We have template specializations
2593  * for dimensions one, two and three, that are used later in DataOut for
2594  * naming the corresponding components:
2595  *
2596 
2597  *
2598  *
2599  * @code
2600  * template <>
2601  * const std::array<std::string, 3> ProblemDescription<1>::component_names{
2602  * {"rho", "m", "E"}};
2603  *
2604  * template <>
2605  * const std::array<std::string, 4> ProblemDescription<2>::component_names{
2606  * {"rho", "m_1", "m_2", "E"}};
2607  *
2608  * template <>
2609  * const std::array<std::string, 5> ProblemDescription<3>::component_names{
2610  * {"rho", "m_1", "m_2", "m_3", "E"}};
2611  *
2612  * @endcode
2613  *
2614  *
2615  * <a name="Initialvalues"></a>
2616  * <h4>Initial values</h4>
2617  *
2618 
2619  *
2620  * The last preparatory step, before we discuss the implementation of the
2621  * forward Euler scheme, is to briefly implement the `InitialValues` class.
2622  *
2623 
2624  *
2625  * In the constructor we initialize all parameters with default values,
2626  * declare all parameters for the `ParameterAcceptor` class and connect the
2627  * <code>parse_parameters_call_back</code> slot to the respective signal.
2628  *
2629 
2630  *
2631  * The <code>parse_parameters_call_back</code> slot will be invoked from
2632  * ParameterAceptor after the call to ParameterAcceptor::initialize(). In
2633  * that regard, its use is appropriate for situations where the
2634  * parameters have to be postprocessed (in some sense) or some
2635  * consistency condition between the parameters has to be checked.
2636  *
2637 
2638  *
2639  *
2640  * @code
2641  * template <int dim>
2642  * InitialValues<dim>::InitialValues(const std::string &subsection)
2643  * : ParameterAcceptor(subsection)
2644  * {
2645  * /* We wire up the slot InitialValues<dim>::parse_parameters_callback to
2646  * the ParameterAcceptor::parse_parameters_call_back signal: */
2648  * std::bind(&InitialValues<dim>::parse_parameters_callback, this));
2649  *
2650  * initial_direction[0] = 1.;
2651  * add_parameter("initial direction",
2652  * initial_direction,
2653  * "Initial direction of the uniform flow field");
2654  *
2655  * initial_1d_state[0] = ProblemDescription<dim>::gamma;
2656  * initial_1d_state[1] = 3.;
2657  * initial_1d_state[2] = 1.;
2658  * add_parameter("initial 1d state",
2659  * initial_1d_state,
2660  * "Initial 1d state (rho, u, p) of the uniform flow field");
2661  * }
2662  *
2663  * @endcode
2664  *
2665  * So far the constructor of <code>InitialValues</code> has defined
2666  * default values for the two private members
2667  * <code>initial_direction</code> and <code>initial_1d_state</code> and
2668  * added them to the parameter list. But we have not defined an
2669  * implementation of the only public member that we really care about,
2670  * which is <code>initial_state()</code> (the function that we are going to
2671  * call to actually evaluate the initial solution at the mesh nodes). At
2672  * the top of the function, we have to ensure that the provided initial
2673  * direction is not the zero vector.
2674  *
2675 
2676  *
2677  * @note As commented, we could have avoided using the method
2678  * <code>parse_parameters_call_back </code> and defined a class member
2679  * <code>setup()</code> in order to define the implementation of
2680  * <code>initial_state()</code>. But for illustrative purposes we want to
2681  * document a different way here and use the call back signal from
2683  *
2684 
2685  *
2686  *
2687  * @code
2688  * template <int dim>
2689  * void InitialValues<dim>::parse_parameters_callback()
2690  * {
2691  * AssertThrow(initial_direction.norm() != 0.,
2692  * ExcMessage(
2693  * "Initial shock front direction is set to the zero vector."));
2694  * initial_direction /= initial_direction.norm();
2695  *
2696  * @endcode
2697  *
2698  * Next, we implement the <code>initial_state</code> function object
2699  * with a lambda function computing a uniform flow field. For this we
2700  * have to translate a given primitive 1d state (density @f$\rho@f$,
2701  * velocity @f$u@f$, and pressure @f$p@f$) into a conserved n-dimensional state
2702  * (density @f$\rho@f$, momentum @f$\mathbf{m}@f$, and total energy @f$E@f$).
2703  *
2704 
2705  *
2706  *
2707  * @code
2708  * initial_state = [this](const Point<dim> & /*point*/, double /*t*/) {
2709  * const double rho = initial_1d_state[0];
2710  * const double u = initial_1d_state[1];
2711  * const double p = initial_1d_state[2];
2712  * static constexpr double gamma = ProblemDescription<dim>::gamma;
2713  *
2714  * state_type state;
2715  *
2716  * state[0] = rho;
2717  * for (unsigned int i = 0; i < dim; ++i)
2718  * state[1 + i] = rho * u * initial_direction[i];
2719  *
2720  * state[dim + 1] = p / (gamma - 1.) + 0.5 * rho * u * u;
2721  *
2722  * return state;
2723  * };
2724  * }
2725  *
2726  * @endcode
2727  *
2728  *
2729  * <a name="TheForwardEulerstep"></a>
2730  * <h4>The Forward Euler step</h4>
2731  *
2732 
2733  *
2734  * The constructor of the <code>%TimeStepping</code> class does not contain
2735  * any surprising code:
2736  *
2737 
2738  *
2739  *
2740  * @code
2741  * template <int dim>
2743  * const MPI_Comm mpi_communicator,
2744  * TimerOutput & computing_timer,
2745  * const OfflineData<dim> & offline_data,
2746  * const InitialValues<dim> &initial_values,
2747  * const std::string & subsection /*= "TimeStepping"*/)
2748  * : ParameterAcceptor(subsection)
2749  * , mpi_communicator(mpi_communicator)
2750  * , computing_timer(computing_timer)
2751  * , offline_data(&offline_data)
2752  * , initial_values(&initial_values)
2753  * {
2754  * cfl_update = 0.80;
2755  * add_parameter("cfl update",
2756  * cfl_update,
2757  * "Relative CFL constant used for update");
2758  * }
2759  *
2760  * @endcode
2761  *
2762  * In the class member <code>prepare()</code> we initialize the temporary
2763  * vector <code>temp</code> and the matrix <code>dij_matrix</code>. The
2764  * vector <code>temp</code> will be used to store the solution update
2765  * temporarily before its contents is swapped with the old vector.
2766  *
2767 
2768  *
2769  *
2770  * @code
2771  * template <int dim>
2772  * void TimeStepping<dim>::prepare()
2773  * {
2774  * TimerOutput::Scope scope(computing_timer,
2775  * "time_stepping - prepare scratch space");
2776  *
2777  * for (auto &it : temporary_vector)
2778  * it.reinit(offline_data->partitioner);
2779  *
2780  * dij_matrix.reinit(offline_data->sparsity_pattern);
2781  * }
2782  *
2783  * @endcode
2784  *
2785  * It is now time to implement the forward Euler step. Given a (writable
2786  * reference) to the old state <code>U</code> at time @f$t@f$ we update the
2787  * state <code>U</code> in place and return the chosen time-step size. We
2788  * first declare a number of read-only references to various different
2789  * variables and data structures. We do this is mainly to have shorter
2790  * variable names (e.g., <code>sparsity</code> instead of
2791  * <code>offline_data->sparsity_pattern</code>).
2792  *
2793 
2794  *
2795  *
2796  * @code
2797  * template <int dim>
2798  * double TimeStepping<dim>::make_one_step(vector_type &U, double t)
2799  * {
2800  * const auto &n_locally_owned = offline_data->n_locally_owned;
2801  * const auto &n_locally_relevant = offline_data->n_locally_relevant;
2802  *
2804  * indices_owned(0, n_locally_owned);
2806  * indices_relevant(0, n_locally_relevant);
2807  *
2808  * const auto &sparsity = offline_data->sparsity_pattern;
2809  *
2810  * const auto &lumped_mass_matrix = offline_data->lumped_mass_matrix;
2811  * const auto &norm_matrix = offline_data->norm_matrix;
2812  * const auto &nij_matrix = offline_data->nij_matrix;
2813  * const auto &cij_matrix = offline_data->cij_matrix;
2814  *
2815  * const auto &boundary_normal_map = offline_data->boundary_normal_map;
2816  *
2817  * @endcode
2818  *
2819  * <b>Step 1</b>: Computing the @f$d_{ij}@f$ graph viscosity matrix.
2820  *
2821 
2822  *
2823  * It is important to highlight that the viscosity matrix has to be
2824  * symmetric, i.e., @f$d_{ij} = d_{ji}@f$. In this regard we note here that
2825  * @f$\int_{\Omega} \nabla \phi_j \phi_i \, \mathrm{d}\mathbf{x}= -
2826  * \int_{\Omega} \nabla \phi_i \phi_j \, \mathrm{d}\mathbf{x}@f$ (or
2827  * equivalently @f$\mathbf{c}_{ij} = - \mathbf{c}_{ji}@f$) provided either
2828  * @f$\mathbf{x}_i@f$ or @f$\mathbf{x}_j@f$ is a support point located away
2829  * from the boundary. In this case we can check that
2830  * @f$\lambda_{\text{max}} (\mathbf{U}_i^{n}, \mathbf{U}_j^{n},
2831  * \textbf{n}_{ij}) = \lambda_{\text{max}} (\mathbf{U}_j^{n},
2832  * \mathbf{U}_i^{n},\textbf{n}_{ji})@f$ by construction, which guarantees
2833  * the property @f$d_{ij} = d_{ji}@f$.
2834  *
2835 
2836  *
2837  * However, if both support points @f$\mathbf{x}_i@f$ or @f$\mathbf{x}_j@f$
2838  * happen to lie on the boundary, then, the equalities @f$\mathbf{c}_{ij} =
2839  * - \mathbf{c}_{ji}@f$ and @f$\lambda_{\text{max}} (\mathbf{U}_i^{n},
2840  * \mathbf{U}_j^{n}, \textbf{n}_{ij}) = \lambda_{\text{max}}
2841  * (\mathbf{U}_j^{n}, \mathbf{U}_i^{n}, \textbf{n}_{ji})@f$ do not
2842  * necessarily hold true. The only mathematically safe solution for this
2843  * dilemma is to compute both of them @f$d_{ij}@f$ and @f$d_{ji}@f$ and
2844  * take the maximum.
2845  *
2846 
2847  *
2848  * Overall, the computation of @f$d_{ij}@f$ is quite expensive. In
2849  * order to save some computing time we exploit the fact that the viscosity
2850  * matrix has to be symmetric (as mentioned above): we only compute
2851  * the upper-triangular entries of @f$d_{ij}@f$ and copy the
2852  * corresponding entries to the lower-triangular counterpart.
2853  *
2854 
2855  *
2856  * We use again parallel::apply_to_subranges() for thread-parallel for
2857  * loops. Pretty much all the ideas for parallel traversal that we
2858  * introduced when discussing the assembly of the matrix
2859  * <code>norm_matrix</code> and the normalization of
2860  * <code>nij_matrix</code> above are used here again.
2861  *
2862 
2863  *
2864  * We define again a "worker" function <code>on_subranges</code> that
2865  * computes the viscosity @f$d_{ij}@f$ for a subrange [i1, i2) of column
2866  * indices:
2867  *
2868  * @code
2869  * {
2870  * TimerOutput::Scope scope(computing_timer,
2871  * "time_stepping - 1 compute d_ij");
2872  *
2873  * const auto on_subranges =
2874  * [&](
2876  * const std_cxx20::ranges::iota_view<unsigned int,
2877  * unsigned int>::iterator i2) {
2878  * for (const auto i :
2880  * *i2))
2881  * {
2882  * const auto U_i = gather(U, i);
2883  *
2884  * @endcode
2885  *
2886  * For a given column index i we iterate over the columns of the
2887  * sparsity pattern from <code>sparsity.begin(i)</code> to
2888  * <code>sparsity.end(i)</code>:
2889  *
2890  * @code
2891  * for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
2892  * {
2893  * const auto j = jt->column();
2894  *
2895  * @endcode
2896  *
2897  * We only compute @f$d_{ij}@f$ if @f$j < i@f$ (upper triangular
2898  * entries) and later copy the values over to @f$d_{ji}@f$.
2899  *
2900  * @code
2901  * if (j >= i)
2902  * continue;
2903  *
2904  * const auto U_j = gather(U, j);
2905  *
2906  * const auto n_ij = gather_get_entry(nij_matrix, jt);
2907  * const double norm = get_entry(norm_matrix, jt);
2908  *
2909  * const auto lambda_max =
2910  * ProblemDescription<dim>::compute_lambda_max(U_i, U_j, n_ij);
2911  *
2912  * double d = norm * lambda_max;
2913  *
2914  * @endcode
2915  *
2916  * If both support points happen to be at the boundary we
2917  * have to compute @f$d_{ji}@f$ as well and then take
2918  * @f$\max(d_{ij},d_{ji})@f$. After this we can finally set the
2919  * upper triangular and lower triangular entries.
2920  *
2921  * @code
2922  * if (boundary_normal_map.count(i) != 0 &&
2923  * boundary_normal_map.count(j) != 0)
2924  * {
2925  * const auto n_ji = gather(nij_matrix, j, i);
2926  * const auto lambda_max_2 =
2927  * ProblemDescription<dim>::compute_lambda_max(U_j,
2928  * U_i,
2929  * n_ji);
2930  * const double norm_2 = norm_matrix(j, i);
2931  *
2932  * d = std::max(d, norm_2 * lambda_max_2);
2933  * }
2934  *
2935  * set_entry(dij_matrix, jt, d);
2936  * dij_matrix(j, i) = d;
2937  * }
2938  * }
2939  * };
2940  *
2941  * parallel::apply_to_subranges(indices_relevant.begin(),
2942  * indices_relevant.end(),
2943  * on_subranges,
2944  * 4096);
2945  * }
2946  *
2947  * @endcode
2948  *
2949  * <b>Step 2</b>: Compute diagonal entries @f$d_{ii}@f$ and
2950  * @f$\tau_{\text{max}}@f$.
2951  *
2952 
2953  *
2954  * So far we have computed all off-diagonal entries of the matrix
2955  * <code>dij_matrix</code>. We still have to fill its diagonal entries
2956  * defined as @f$d_{ii}^n = - \sum_{j \in \mathcal{I}(i)\backslash \{i\}}
2957  * d_{ij}^n@f$. We use again parallel::apply_to_subranges() for this
2958  * purpose. While computing the @f$d_{ii}@f$s we also determine the
2959  * largest admissible time-step, which is defined as
2960  * \f[
2961  * \tau_n \dealcoloneq c_{\text{cfl}}\,\min_{i\in\mathcal{V}}
2962  * \left(\frac{m_i}{-2\,d_{ii}^{n}}\right) \, .
2963  * \f]
2964  * Note that the operation @f$\min_{i \in \mathcal{V}}@f$ is intrinsically
2965  * global, it operates on all nodes: first we have to take the minimum
2966  * over all threads (of a given node) and then we have to take the
2967  * minimum over all MPI processes. In the current implementation:
2968  * - We store <code>tau_max</code> (per node) as
2969  * <a
2970  * href="http://www.cplusplus.com/reference/atomic/atomic/"><code>std::atomic<double></code></a>.
2971  * The internal implementation of <code>std::atomic</code> will take
2972  * care of guarding any possible race condition when more than one
2973  * thread attempts to read and/or write <code>tau_max</code> at the
2974  * same time.
2975  * - In order to take the minimum over all MPI process we use the utility
2976  * function <code>Utilities::MPI::min</code>.
2977  *
2978 
2979  *
2980  *
2981  * @code
2982  * std::atomic<double> tau_max{std::numeric_limits<double>::infinity()};
2983  *
2984  * {
2985  * TimerOutput::Scope scope(computing_timer,
2986  * "time_stepping - 2 compute d_ii, and tau_max");
2987  *
2988  * @endcode
2989  *
2990  * on_subranges() will be executed on every thread individually. The
2991  * variable <code>tau_max_on_subrange</code> is thus stored thread
2992  * locally.
2993  *
2994 
2995  *
2996  *
2997  * @code
2998  * const auto on_subranges =
2999  * [&](
3000  * std_cxx20::ranges::iota_view<unsigned int, unsigned int>::iterator i1,
3001  * const std_cxx20::ranges::iota_view<unsigned int,
3002  * unsigned int>::iterator i2) {
3003  * double tau_max_on_subrange = std::numeric_limits<double>::infinity();
3004  *
3005  * for (const auto i :
3007  * *i2))
3008  * {
3009  * double d_sum = 0.;
3010  *
3011  * for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
3012  * {
3013  * const auto j = jt->column();
3014  *
3015  * if (j == i)
3016  * continue;
3017  *
3018  * d_sum -= get_entry(dij_matrix, jt);
3019  * }
3020  *
3021  * @endcode
3022  *
3023  * We store the negative sum of the d_ij entries at the
3024  * diagonal position
3025  *
3026  * @code
3027  * dij_matrix.diag_element(i) = d_sum;
3028  * @endcode
3029  *
3030  * and compute the maximal local time-step size
3031  * <code>tau</code>:
3032  *
3033  * @code
3034  * const double mass = lumped_mass_matrix.diag_element(i);
3035  * const double tau = cfl_update * mass / (-2. * d_sum);
3036  * tau_max_on_subrange = std::min(tau_max_on_subrange, tau);
3037  * }
3038  *
3039  * @endcode
3040  *
3041  * <code>tau_max_on_subrange</code> contains the largest possible
3042  * time-step size computed for the (thread local) subrange. At this
3043  * point we have to synchronize the value over all threads. This is
3044  * were we use the <a
3045  * href="http://www.cplusplus.com/reference/atomic/atomic/"><code>std::atomic<double></code></a>
3046  * <i>compare exchange</i> update mechanism:
3047  *
3048  * @code
3049  * double current_tau_max = tau_max.load();
3050  * while (current_tau_max > tau_max_on_subrange &&
3051  * !tau_max.compare_exchange_weak(current_tau_max,
3052  * tau_max_on_subrange))
3053  * ;
3054  * };
3055  *
3056  * parallel::apply_to_subranges(indices_relevant.begin(),
3057  * indices_relevant.end(),
3058  * on_subranges,
3059  * 4096);
3060  *
3061  * @endcode
3062  *
3063  * After all threads have finished we can simply synchronize the
3064  * value over all MPI processes:
3065  *
3066 
3067  *
3068  *
3069  * @code
3070  * tau_max.store(Utilities::MPI::min(tau_max.load(), mpi_communicator));
3071  *
3072  * @endcode
3073  *
3074  * This is a good point to verify that the computed
3075  * <code>tau_max</code> is indeed a valid floating point number.
3076  *
3077  * @code
3078  * AssertThrow(
3079  * !std::isnan(tau_max.load()) && !std::isinf(tau_max.load()) &&
3080  * tau_max.load() > 0.,
3081  * ExcMessage(
3082  * "I'm sorry, Dave. I'm afraid I can't do that. - We crashed."));
3083  * }
3084  *
3085  * @endcode
3086  *
3087  * <b>Step 3</b>: Perform update.
3088  *
3089 
3090  *
3091  * At this point, we have computed all viscosity coefficients @f$d_{ij}@f$
3092  * and we know the maximal admissible time-step size
3093  * @f$\tau_{\text{max}}@f$. This means we can now compute the update:
3094  *
3095 
3096  *
3097  * \f[
3098  * \mathbf{U}_i^{n+1} = \mathbf{U}_i^{n} - \frac{\tau_{\text{max}} }{m_i}
3099  * \sum_{j \in \mathcal{I}(i)} (\mathbb{f}(\mathbf{U}_j^{n}) -
3100  * \mathbb{f}(\mathbf{U}_i^{n})) \cdot \mathbf{c}_{ij} - d_{ij}
3101  * (\mathbf{U}_j^{n} - \mathbf{U}_i^{n})
3102  * \f]
3103  *
3104 
3105  *
3106  * This update formula is slightly different from what was discussed in
3107  * the introduction (in the pseudo-code). However, it can be shown that
3108  * both equations are algebraically equivalent (they will produce the
3109  * same numerical values). We favor this second formula since it has
3110  * natural cancellation properties that might help avoid numerical
3111  * artifacts.
3112  *
3113 
3114  *
3115  *
3116  * @code
3117  * {
3118  * TimerOutput::Scope scope(computing_timer,
3119  * "time_stepping - 3 perform update");
3120  *
3121  * const auto on_subranges =
3122  * [&](
3124  * const std_cxx20::ranges::iota_view<unsigned int,
3125  * unsigned int>::iterator i2) {
3126  * for (const auto i : boost::make_iterator_range(i1, i2))
3127  * {
3128  * Assert(i < n_locally_owned, ExcInternalError());
3129  *
3130  * const auto U_i = gather(U, i);
3131  *
3132  * const auto f_i = ProblemDescription<dim>::flux(U_i);
3133  * const double m_i = lumped_mass_matrix.diag_element(i);
3134  *
3135  * auto U_i_new = U_i;
3136  *
3137  * for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
3138  * {
3139  * const auto j = jt->column();
3140  *
3141  * const auto U_j = gather(U, j);
3142  * const auto f_j = ProblemDescription<dim>::flux(U_j);
3143  *
3144  * const auto c_ij = gather_get_entry(cij_matrix, jt);
3145  * const auto d_ij = get_entry(dij_matrix, jt);
3146  *
3147  * for (unsigned int k = 0; k < problem_dimension; ++k)
3148  * {
3149  * U_i_new[k] +=
3150  * tau_max / m_i *
3151  * (-(f_j[k] - f_i[k]) * c_ij + d_ij * (U_j[k] - U_i[k]));
3152  * }
3153  * }
3154  *
3155  * scatter(temporary_vector, U_i_new, i);
3156  * }
3157  * };
3158  *
3159  * parallel::apply_to_subranges(indices_owned.begin(),
3160  * indices_owned.end(),
3161  * on_subranges,
3162  * 4096);
3163  * }
3164  *
3165  * @endcode
3166  *
3167  * <b>Step 4</b>: Fix up boundary states.
3168  *
3169 
3170  *
3171  * As a last step in the Forward Euler method, we have to fix up all
3172  * boundary states. As discussed in the intro we
3173  * - advance in time satisfying no boundary condition at all,
3174  * - at the end of the time step enforce boundary conditions strongly
3175  * in a post-processing step.
3176  *
3177 
3178  *
3179  * Here, we compute the correction
3180  * \f[
3181  * \mathbf{m}_i \dealcoloneq \mathbf{m}_i - (\boldsymbol{\nu}_i \cdot
3182  * \mathbf{m}_i) \boldsymbol{\nu}_i,
3183  * \f]
3184  * which removes the normal component of @f$\mathbf{m}@f$.
3185  *
3186 
3187  *
3188  *
3189  * @code
3190  * {
3191  * TimerOutput::Scope scope(computing_timer,
3192  * "time_stepping - 4 fix boundary states");
3193  *
3194  * for (auto it : boundary_normal_map)
3195  * {
3196  * const auto i = it.first;
3197  *
3198  * @endcode
3199  *
3200  * We only iterate over the locally owned subset:
3201  *
3202  * @code
3203  * if (i >= n_locally_owned)
3204  * continue;
3205  *
3206  * const auto &normal = std::get<0>(it.second);
3207  * const auto &id = std::get<1>(it.second);
3208  * const auto &position = std::get<2>(it.second);
3209  *
3210  * auto U_i = gather(temporary_vector, i);
3211  *
3212  * @endcode
3213  *
3214  * On free slip boundaries we remove the normal component of the
3215  * momentum:
3216  *
3217  * @code
3218  * if (id == Boundaries::free_slip)
3219  * {
3220  * auto m = ProblemDescription<dim>::momentum(U_i);
3221  * m -= (m * normal) * normal;
3222  * for (unsigned int k = 0; k < dim; ++k)
3223  * U_i[k + 1] = m[k];
3224  * }
3225  *
3226  * @endcode
3227  *
3228  * On Dirichlet boundaries we enforce initial conditions
3229  * strongly:
3230  *
3231  * @code
3232  * else if (id == Boundaries::dirichlet)
3233  * {
3234  * U_i = initial_values->initial_state(position, t + tau_max);
3235  * }
3236  *
3237  * scatter(temporary_vector, U_i, i);
3238  * }
3239  * }
3240  *
3241  * @endcode
3242  *
3243  * <b>Step 5</b>: We now update the ghost layer over all MPI ranks,
3244  * swap the temporary vector with the solution vector <code>U</code>
3245  * (that will get returned by reference) and return the chosen
3246  * time-step size @f$\tau_{\text{max}}@f$:
3247  *
3248 
3249  *
3250  *
3251  * @code
3252  * for (auto &it : temporary_vector)
3253  * it.update_ghost_values();
3254  *
3255  * U.swap(temporary_vector);
3256  *
3257  * return tau_max;
3258  * }
3259  *
3260  * @endcode
3261  *
3262  *
3263  * <a name="Schlierenpostprocessing"></a>
3264  * <h4>Schlieren postprocessing</h4>
3265  *
3266 
3267  *
3268  * At various intervals we will output the current state <code>U</code>
3269  * of the solution together with a so-called Schlieren plot.
3270  * The constructor of the <code>SchlierenPostprocessor</code> class again
3271  * contains no surprises. We simply supply default values to and register
3272  * two parameters:
3273  * - schlieren_beta:
3274  * is an ad-hoc positive amplification factor in order to enhance the
3275  * contrast in the visualization. Its actual value is a matter of
3276  * taste.
3277  * - schlieren_index: is an integer indicating which component of the
3278  * state @f$[\rho, \mathbf{m},E]@f$ are we going to use in order to generate
3279  * the visualization.
3280  *
3281 
3282  *
3283  *
3284  * @code
3285  * template <int dim>
3286  * SchlierenPostprocessor<dim>::SchlierenPostprocessor(
3287  * const MPI_Comm mpi_communicator,
3288  * TimerOutput & computing_timer,
3289  * const OfflineData<dim> &offline_data,
3290  * const std::string & subsection /*= "SchlierenPostprocessor"*/)
3291  * : ParameterAcceptor(subsection)
3292  * , mpi_communicator(mpi_communicator)
3293  * , computing_timer(computing_timer)
3294  * , offline_data(&offline_data)
3295  * {
3296  * schlieren_beta = 10.;
3297  * add_parameter("schlieren beta",
3298  * schlieren_beta,
3299  * "Beta factor used in Schlieren-type postprocessor");
3300  *
3301  * schlieren_index = 0;
3302  * add_parameter("schlieren index",
3303  * schlieren_index,
3304  * "Use the corresponding component of the state vector for the "
3305  * "schlieren plot");
3306  * }
3307  *
3308  * @endcode
3309  *
3310  * Again, the <code>prepare()</code> function initializes two temporary
3311  * the vectors (<code>r</code> and <code>schlieren</code>).
3312  *
3313 
3314  *
3315  *
3316  * @code
3317  * template <int dim>
3318  * void SchlierenPostprocessor<dim>::prepare()
3319  * {
3320  * TimerOutput::Scope scope(computing_timer,
3321  * "schlieren_postprocessor - prepare scratch space");
3322  *
3323  * r.reinit(offline_data->n_locally_relevant);
3324  * schlieren.reinit(offline_data->partitioner);
3325  * }
3326  *
3327  * @endcode
3328  *
3329  * We now discuss the implementation of the class member
3330  * <code>SchlierenPostprocessor<dim>::compute_schlieren()</code>, which
3331  * basically takes a component of the state vector <code>U</code> and
3332  * computes the Schlieren indicator for such component (the formula of the
3333  * Schlieren indicator can be found just before the declaration of the class
3334  * <code>SchlierenPostprocessor</code>). We start by noting
3335  * that this formula requires the "nodal gradients" @f$\nabla r_j@f$.
3336  * However, nodal values of gradients are not defined for @f$\mathcal{C}^0@f$
3337  * finite element functions. More generally, pointwise values of
3338  * gradients are not defined for @f$W^{1,p}(\Omega)@f$ functions. The
3339  * simplest technique we can use to recover gradients at nodes is
3340  * weighted-averaging i.e.
3341  *
3342 
3343  *
3344  * \f[ \nabla r_j \dealcoloneq \frac{1}{\int_{S_i} \omega_i(\mathbf{x}) \,
3345  * \mathrm{d}\mathbf{x}}
3346  * \int_{S_i} r_h(\mathbf{x}) \omega_i(\mathbf{x}) \, \mathrm{d}\mathbf{x}
3347  * \ \ \ \ \ \mathbf{(*)} \f]
3348  *
3349 
3350  *
3351  * where @f$S_i@f$ is the support of the shape function @f$\phi_i@f$, and
3352  * @f$\omega_i(\mathbf{x})@f$ is the weight. The weight could be any
3353  * positive function such as
3354  * @f$\omega_i(\mathbf{x}) \equiv 1@f$ (that would allow us to recover the usual
3355  * notion of mean value). But as usual, the goal is to reuse the off-line
3356  * data as much as possible. In this sense, the most natural
3357  * choice of weight is @f$\omega_i = \phi_i@f$. Inserting this choice of
3358  * weight and the expansion @f$r_h(\mathbf{x}) = \sum_{j \in \mathcal{V}}
3359  * r_j \phi_j(\mathbf{x})@f$ into @f$\mathbf{(*)}@f$ we get :
3360  *
3361 
3362  *
3363  * \f[ \nabla r_j \dealcoloneq \frac{1}{m_i} \sum_{j \in \mathcal{I}(i)} r_j
3364  * \mathbf{c}_{ij} \ \ \ \ \ \mathbf{(**)} \, . \f]
3365  *
3366 
3367  *
3368  * Using this last formula we can recover averaged nodal gradients without
3369  * resorting to any form of quadrature. This idea aligns quite well with
3370  * the whole spirit of edge-based schemes (or algebraic schemes) where
3371  * we want to operate on matrices and vectors as directly as
3372  * it could be possible avoiding by all means assembly of bilinear
3373  * forms, cell-loops, quadrature, or any other
3374  * intermediate construct/operation between the input arguments (the state
3375  * from the previous time-step) and the actual matrices and vectors
3376  * required to compute the update.
3377  *
3378 
3379  *
3380  * The second thing to note is that we have to compute global minimum and
3381  * maximum @f$\max_j |\nabla r_j|@f$ and @f$\min_j |\nabla r_j|@f$. Following the
3382  * same ideas used to compute the time step size in the class member
3383  * <code>%TimeStepping%<dim%>::%step()</code> we define @f$\max_j |\nabla r_j|@f$
3384  * and @f$\min_j |\nabla r_j|@f$ as atomic doubles in order to resolve any
3385  * conflicts between threads. As usual, we use
3386  * <code>Utilities::MPI::max()</code> and
3387  * <code>Utilities::MPI::min()</code> to find the global maximum/minimum
3388  * among all MPI processes.
3389  *
3390 
3391  *
3392  * Finally, it is not possible to compute the Schlieren indicator in a single
3393  * loop over all nodes. The entire operation requires two loops over nodes:
3394  *
3395 
3396  *
3397  * - The first loop computes @f$|\nabla r_i|@f$ for all @f$i \in \mathcal{V}@f$ in
3398  * the mesh, and the bounds @f$\max_j |\nabla r_j|@f$ and
3399  * @f$\min_j |\nabla r_j|@f$.
3400  * - The second loop finally computes the Schlieren indicator using the
3401  * formula
3402  *
3403 
3404  *
3405  * \f[ \text{schlieren}[i] = e^{\beta \frac{ |\nabla r_i|
3406  * - \min_j |\nabla r_j| }{\max_j |\nabla r_j| - \min_j |\nabla r_j| } }
3407  * \, . \f]
3408  *
3409 
3410  *
3411  * This means that we will have to define two workers
3412  * <code>on_subranges</code> for each one of these stages.
3413  *
3414 
3415  *
3416  *
3417  * @code
3418  * template <int dim>
3419  * void SchlierenPostprocessor<dim>::compute_schlieren(const vector_type &U)
3420  * {
3421  * TimerOutput::Scope scope(
3422  * computing_timer, "schlieren_postprocessor - compute schlieren plot");
3423  *
3424  * const auto &sparsity = offline_data->sparsity_pattern;
3425  * const auto &lumped_mass_matrix = offline_data->lumped_mass_matrix;
3426  * const auto &cij_matrix = offline_data->cij_matrix;
3427  * const auto &boundary_normal_map = offline_data->boundary_normal_map;
3428  * const auto &n_locally_owned = offline_data->n_locally_owned;
3429  *
3430  * const auto indices =
3432  * n_locally_owned);
3433  *
3434  * @endcode
3435  *
3436  * We define the r_i_max and r_i_min in the current MPI process as
3437  * atomic doubles in order to avoid race conditions between threads:
3438  *
3439  * @code
3440  * std::atomic<double> r_i_max{0.};
3441  * std::atomic<double> r_i_min{std::numeric_limits<double>::infinity()};
3442  *
3443  * @endcode
3444  *
3445  * First loop: compute the averaged gradient at each node and the
3446  * global maxima and minima of the gradients.
3447  *
3448  * @code
3449  * {
3450  * const auto on_subranges =
3451  * [&](
3453  * const std_cxx20::ranges::iota_view<unsigned int,
3454  * unsigned int>::iterator i2) {
3455  * double r_i_max_on_subrange = 0.;
3456  * double r_i_min_on_subrange = std::numeric_limits<double>::infinity();
3457  *
3458  * for (const auto i : boost::make_iterator_range(i1, i2))
3459  * {
3460  * Assert(i < n_locally_owned, ExcInternalError());
3461  *
3462  * Tensor<1, dim> r_i;
3463  *
3464  * for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
3465  * {
3466  * const auto j = jt->column();
3467  *
3468  * if (i == j)
3469  * continue;
3470  *
3471  * const auto U_js = U[schlieren_index].local_element(j);
3472  * const auto c_ij = gather_get_entry(cij_matrix, jt);
3473  * r_i += c_ij * U_js;
3474  * }
3475  *
3476  * @endcode
3477  *
3478  * We fix up the gradient r_i at free slip boundaries similarly to
3479  * how we fixed up boundary states in the forward Euler step.
3480  * This avoids sharp, artificial gradients in the Schlieren
3481  * plot at free slip boundaries and is a purely cosmetic choice.
3482  *
3483 
3484  *
3485  *
3486  * @code
3487  * const auto bnm_it = boundary_normal_map.find(i);
3488  * if (bnm_it != boundary_normal_map.end())
3489  * {
3490  * const auto &normal = std::get<0>(bnm_it->second);
3491  * const auto &id = std::get<1>(bnm_it->second);
3492  *
3493  * if (id == Boundaries::free_slip)
3494  * r_i -= 1. * (r_i * normal) * normal;
3495  * else
3496  * r_i = 0.;
3497  * }
3498  *
3499  * @endcode
3500  *
3501  * We remind the reader that we are not interested in the nodal
3502  * gradients per se. We only want their norms in order to
3503  * compute the Schlieren indicator (weighted with the lumped
3504  * mass matrix @f$m_i@f$):
3505  *
3506  * @code
3507  * const double m_i = lumped_mass_matrix.diag_element(i);
3508  * r[i] = r_i.norm() / m_i;
3509  * r_i_max_on_subrange = std::max(r_i_max_on_subrange, r[i]);
3510  * r_i_min_on_subrange = std::min(r_i_min_on_subrange, r[i]);
3511  * }
3512  *
3513  * @endcode
3514  *
3515  * We compare the current_r_i_max and current_r_i_min (in the
3516  * current subrange) with r_i_max and r_i_min (for the current MPI
3517  * process) and update them if necessary:
3518  *
3519 
3520  *
3521  *
3522  * @code
3523  * double current_r_i_max = r_i_max.load();
3524  * while (current_r_i_max < r_i_max_on_subrange &&
3525  * !r_i_max.compare_exchange_weak(current_r_i_max,
3526  * r_i_max_on_subrange))
3527  * ;
3528  *
3529  * double current_r_i_min = r_i_min.load();
3530  * while (current_r_i_min > r_i_min_on_subrange &&
3531  * !r_i_min.compare_exchange_weak(current_r_i_min,
3532  * r_i_min_on_subrange))
3533  * ;
3534  * };
3535  *
3536  * parallel::apply_to_subranges(indices.begin(),
3537  * indices.end(),
3538  * on_subranges,
3539  * 4096);
3540  * }
3541  *
3542  * @endcode
3543  *
3544  * And synchronize <code>r_i_max</code> and <code>r_i_min</code> over
3545  * all MPI processes.
3546  *
3547 
3548  *
3549  *
3550  * @code
3551  * r_i_max.store(Utilities::MPI::max(r_i_max.load(), mpi_communicator));
3552  * r_i_min.store(Utilities::MPI::min(r_i_min.load(), mpi_communicator));
3553  *
3554  * @endcode
3555  *
3556  * Second loop: we now have the vector <code>r</code> and the scalars
3557  * <code>r_i_max</code> and <code>r_i_min</code> at our disposal. We
3558  * are thus in a position to actually compute the Schlieren indicator.
3559  *
3560 
3561  *
3562  *
3563  * @code
3564  * {
3565  * const auto on_subranges =
3566  * [&](
3568  * const std_cxx20::ranges::iota_view<unsigned int,
3569  * unsigned int>::iterator i2) {
3570  * for (const auto i : boost::make_iterator_range(i1, i2))
3571  * {
3572  * Assert(i < n_locally_owned, ExcInternalError());
3573  *
3574  * schlieren.local_element(i) =
3575  * 1. - std::exp(-schlieren_beta * (r[i] - r_i_min) /
3576  * (r_i_max - r_i_min));
3577  * }
3578  * };
3579  *
3580  * parallel::apply_to_subranges(indices.begin(),
3581  * indices.end(),
3582  * on_subranges,
3583  * 4096);
3584  * }
3585  *
3586  * @endcode
3587  *
3588  * And finally, exchange ghost elements.
3589  *
3590  * @code
3591  * schlieren.update_ghost_values();
3592  * }
3593  *
3594  * @endcode
3595  *
3596  *
3597  * <a name="Themainloop"></a>
3598  * <h4>The main loop</h4>
3599  *
3600 
3601  *
3602  * With all classes implemented it is time to create an instance of
3603  * <code>Discretization<dim></code>, <code>OfflineData<dim></code>,
3604  * <code>InitialValues<dim></code>, <code>TimeStepping<dim></code>, and
3605  * <code>SchlierenPostprocessor<dim></code>, and run the forward Euler
3606  * step in a loop.
3607  *
3608 
3609  *
3610  * In the constructor of <code>MainLoop<dim></code> we now initialize an
3611  * instance of all classes, and declare a number of parameters
3612  * controlling output. Most notable, we declare a boolean parameter
3613  * <code>resume</code> that will control whether the program attempts to
3614  * restart from an interrupted computation, or not.
3615  *
3616 
3617  *
3618  *
3619  * @code
3620  * template <int dim>
3621  * MainLoop<dim>::MainLoop(const MPI_Comm mpi_communicator)
3622  * : ParameterAcceptor("A - MainLoop")
3623  * , mpi_communicator(mpi_communicator)
3624  * , computing_timer(mpi_communicator,
3625  * timer_output,
3628  * , pcout(std::cout, Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
3629  * , discretization(mpi_communicator, computing_timer, "B - Discretization")
3630  * , offline_data(mpi_communicator,
3631  * computing_timer,
3632  * discretization,
3633  * "C - OfflineData")
3634  * , initial_values("D - InitialValues")
3635  * , time_stepping(mpi_communicator,
3636  * computing_timer,
3637  * offline_data,
3638  * initial_values,
3639  * "E - TimeStepping")
3640  * , schlieren_postprocessor(mpi_communicator,
3641  * computing_timer,
3642  * offline_data,
3643  * "F - SchlierenPostprocessor")
3644  * {
3645  * base_name = "test";
3646  * add_parameter("basename", base_name, "Base name for all output files");
3647  *
3648  * t_final = 4.;
3649  * add_parameter("final time", t_final, "Final time");
3650  *
3651  * output_granularity = 0.02;
3652  * add_parameter("output granularity",
3653  * output_granularity,
3654  * "time interval for output");
3655  *
3656  * #ifdef DEAL_II_WITH_THREADS
3657  * asynchronous_writeback = true;
3658  * #else
3659  * asynchronous_writeback = false;
3660  * #endif
3661  * add_parameter("asynchronous writeback",
3662  * asynchronous_writeback,
3663  * "Write out solution in a background thread performing IO");
3664  *
3665  * resume = false;
3666  * add_parameter("resume", resume, "Resume an interrupted computation.");
3667  * }
3668  *
3669  * @endcode
3670  *
3671  * We start by implementing a helper function <code>print_head()</code>
3672  * in an anonymous namespace that is used to output messages in the
3673  * terminal with some nice formatting.
3674  *
3675 
3676  *
3677  *
3678  * @code
3679  * namespace
3680  * {
3681  * void print_head(ConditionalOStream &pcout,
3682  * const std::string & header,
3683  * const std::string & secondary = "")
3684  * {
3685  * const auto header_size = header.size();
3686  * const auto padded_header = std::string((34 - header_size) / 2, ' ') +
3687  * header +
3688  * std::string((35 - header_size) / 2, ' ');
3689  *
3690  * const auto secondary_size = secondary.size();
3691  * const auto padded_secondary =
3692  * std::string((34 - secondary_size) / 2, ' ') + secondary +
3693  * std::string((35 - secondary_size) / 2, ' ');
3694  *
3695  * /* clang-format off */
3696  * pcout << std::endl;
3697  * pcout << " ####################################################" << std::endl;
3698  * pcout << " ######### #########" << std::endl;
3699  * pcout << " #########" << padded_header << "#########" << std::endl;
3700  * pcout << " #########" << padded_secondary << "#########" << std::endl;
3701  * pcout << " ######### #########" << std::endl;
3702  * pcout << " ####################################################" << std::endl;
3703  * pcout << std::endl;
3704  * /* clang-format on */
3705  * }
3706  * } // namespace
3707  *
3708  * @endcode
3709  *
3710  * With <code>print_head</code> in place it is now time to implement the
3711  * <code>MainLoop<dim>::run()</code> that contains the main loop of our
3712  * program.
3713  *
3714 
3715  *
3716  *
3717  * @code
3718  * template <int dim>
3719  * void MainLoop<dim>::run()
3720  * {
3721  * @endcode
3722  *
3723  * We start by reading in parameters and initializing all objects. We
3724  * note here that the call to ParameterAcceptor::initialize reads in
3725  * all parameters from the parameter file (whose name is given as a
3726  * string argument). ParameterAcceptor handles a global
3727  * ParameterHandler that is initialized with subsections and parameter
3728  * declarations for all class instances that are derived from
3729  * ParameterAceptor. The call to initialize enters the subsection for
3730  * each each derived class, and sets all variables that were added
3732  *
3733 
3734  *
3735  *
3736  * @code
3737  * pcout << "Reading parameters and allocating objects... " << std::flush;
3738  *
3739  * ParameterAcceptor::initialize("step-69.prm");
3740  * pcout << "done" << std::endl;
3741  *
3742  * @endcode
3743  *
3744  * Next we create the triangulation, assemble all matrices, set up
3745  * scratch space, and initialize the DataOut<dim> object:
3746  *
3747 
3748  *
3749  *
3750  * @code
3751  * {
3752  * print_head(pcout, "create triangulation");
3753  * discretization.setup();
3754  *
3755  * pcout << "Number of active cells: "
3756  * << discretization.triangulation.n_global_active_cells()
3757  * << std::endl;
3758  *
3759  * print_head(pcout, "compute offline data");
3760  * offline_data.setup();
3761  * offline_data.assemble();
3762  *
3763  * pcout << "Number of degrees of freedom: "
3764  * << offline_data.dof_handler.n_dofs() << std::endl;
3765  *
3766  * print_head(pcout, "set up time step");
3767  * time_stepping.prepare();
3768  * schlieren_postprocessor.prepare();
3769  * }
3770  *
3771  * @endcode
3772  *
3773  * We will store the current time and state in the variable
3774  * <code>t</code> and vector <code>U</code>:
3775  *
3776 
3777  *
3778  *
3779  * @code
3780  * double t = 0.;
3781  * unsigned int output_cycle = 0;
3782  *
3783  * print_head(pcout, "interpolate initial values");
3784  * vector_type U = interpolate_initial_values();
3785  *
3786  * @endcode
3787  *
3788  *
3789  * <a name="Resume"></a>
3790  * <h5>Resume</h5>
3791  *
3792 
3793  *
3794  * By default the boolean <code>resume</code> is set to false, i.e. the
3795  * following code snippet is not run. However, if
3796  * <code>resume==true</code> we indicate that we have indeed an
3797  * interrupted computation and the program shall restart by reading in
3798  * an old state consisting of <code>t</code>,
3799  * <code>output_cycle</code>, and <code>U</code> from a checkpoint
3800  * file. These checkpoint files will be created in the
3801  * <code>output()</code> routine discussed below.
3802  *
3803 
3804  *
3805  *
3806  * @code
3807  * if (resume)
3808  * {
3809  * print_head(pcout, "restore interrupted computation");
3810  *
3811  * const unsigned int i =
3812  * discretization.triangulation.locally_owned_subdomain();
3813  *
3814  * const std::string name = base_name + "-checkpoint-" +
3815  * Utilities::int_to_string(i, 4) + ".archive";
3816  * std::ifstream file(name, std::ios::binary);
3817  *
3818  * @endcode
3819  *
3820  * We use a <code>boost::archive</code> to store and read in the
3821  * contents the checkpointed state.
3822  *
3823 
3824  *
3825  *
3826  * @code
3827  * boost::archive::binary_iarchive ia(file);
3828  * ia >> t >> output_cycle;
3829  *
3830  * for (auto &it1 : U)
3831  * {
3832  * @endcode
3833  *
3834  * <code>it1</code> iterates over all components of the state
3835  * vector <code>U</code>. We read in every entry of the
3836  * component in sequence and update the ghost layer afterwards:
3837  *
3838  * @code
3839  * for (auto &it2 : it1)
3840  * ia >> it2;
3841  * it1.update_ghost_values();
3842  * }
3843  * }
3844  *
3845  * @endcode
3846  *
3847  * With either the initial state set up, or an interrupted state
3848  * restored it is time to enter the main loop:
3849  *
3850 
3851  *
3852  *
3853  * @code
3854  * output(U, base_name, t, output_cycle++);
3855  *
3856  * print_head(pcout, "enter main loop");
3857  *
3858  * for (unsigned int cycle = 1; t < t_final; ++cycle)
3859  * {
3860  * @endcode
3861  *
3862  * We first print an informative status message
3863  *
3864 
3865  *
3866  *
3867  * @code
3868  * std::ostringstream head;
3869  * std::ostringstream secondary;
3870  *
3871  * head << "Cycle " << Utilities::int_to_string(cycle, 6) << " ("
3872  * << std::fixed << std::setprecision(1) << t / t_final * 100
3873  * << "%)";
3874  * secondary << "at time t = " << std::setprecision(8) << std::fixed << t;
3875  *
3876  * print_head(pcout, head.str(), secondary.str());
3877  *
3878  * @endcode
3879  *
3880  * and then perform a single forward Euler step. Note that the
3881  * state vector <code>U</code> is updated in place and that
3882  * <code>time_stepping.make_one_step()</code> returns the chosen step
3883  * size.
3884  *
3885 
3886  *
3887  *
3888  * @code
3889  * t += time_stepping.make_one_step(U, t);
3890  *
3891  * @endcode
3892  *
3893  * Post processing, generating output and writing out the current
3894  * state is a CPU and IO intensive task that we cannot afford to do
3895  * every time step - in particular with explicit time stepping. We
3896  * thus only schedule output by calling the <code>output()</code>
3897  * function if we are past a threshold set by
3898  * <code>output_granularity</code>.
3899  *
3900 
3901  *
3902  *
3903  * @code
3904  * if (t > output_cycle * output_granularity)
3905  * {
3906  * output(U, base_name, t, output_cycle, true);
3907  * ++output_cycle;
3908  * }
3909  * }
3910  *
3911  * @endcode
3912  *
3913  * We wait for any remaining background output thread to finish before
3914  * printing a summary and exiting.
3915  *
3916  * @code
3917  * if (background_thread_state.valid())
3918  * background_thread_state.wait();
3919  *
3920  * computing_timer.print_summary();
3921  * pcout << timer_output.str() << std::endl;
3922  * }
3923  *
3924  * @endcode
3925  *
3926  * The <code>interpolate_initial_values</code> takes an initial time "t"
3927  * as input argument and populates a state vector <code>U</code> with the
3928  * help of the <code>InitialValues<dim>::initial_state</code> object.
3929  *
3930 
3931  *
3932  *
3933  * @code
3934  * template <int dim>
3935  * typename MainLoop<dim>::vector_type
3936  * MainLoop<dim>::interpolate_initial_values(const double t)
3937  * {
3938  * pcout << "MainLoop<dim>::interpolate_initial_values(t = " << t << ")"
3939  * << std::endl;
3940  * TimerOutput::Scope scope(computing_timer,
3941  * "main_loop - setup scratch space");
3942  *
3943  * vector_type U;
3944  *
3945  * for (auto &it : U)
3946  * it.reinit(offline_data.partitioner);
3947  *
3948  * constexpr auto problem_dimension =
3949  * ProblemDescription<dim>::problem_dimension;
3950  *
3951  * @endcode
3952  *
3953  * The function signature of
3954  * <code>InitialValues<dim>::initial_state</code> is not quite right
3955  * for VectorTools::interpolate(). We work around this issue by, first,
3956  * creating a lambda function that for a given position <code>x</code>
3957  * returns just the value of the <code>i</code>th component. This
3958  * lambda in turn is converted to a ::Function with the help of
3959  * the ScalarFunctionFromFunctionObject wrapper.
3960  *
3961 
3962  *
3963  *
3964  * @code
3965  * for (unsigned int i = 0; i < problem_dimension; ++i)
3966  * VectorTools::interpolate(offline_data.dof_handler,
3968  * [&](const Point<dim> &x) {
3969  * return initial_values.initial_state(x, t)[i];
3970  * }),
3971  * U[i]);
3972  *
3973  * for (auto &it : U)
3974  * it.update_ghost_values();
3975  *
3976  * return U;
3977  * }
3978  *
3979  * @endcode
3980  *
3981  *
3982  * <a name="Outputandcheckpointing"></a>
3983  * <h5>Output and checkpointing</h5>
3984  *
3985 
3986  *
3987  * Writing out the final vtk files is quite an IO intensive task that can
3988  * stall the main loop for a while. In order to avoid this we use an <a
3989  * href="https://en.wikipedia.org/wiki/Asynchronous_I/O">asynchronous
3990  * IO</a> strategy by creating a background thread that will perform IO
3991  * while the main loop is allowed to continue. In order for this to work
3992  * we have to be mindful of two things:
3993  * - Before running the <code>output_worker</code> thread, we have to create
3994  * a copy of the state vector <code>U</code>. We store it in the
3995  * vector <code>output_vector</code>.
3996  * - We have to avoid any MPI communication in the background thread,
3997  * otherwise the program might deadlock. This implies that we have to
3998  * run the postprocessing outside of the worker thread.
3999  *
4000 
4001  *
4002  *
4003  * @code
4004  * template <int dim>
4005  * void MainLoop<dim>::output(const typename MainLoop<dim>::vector_type &U,
4006  * const std::string & name,
4007  * const double t,
4008  * const unsigned int cycle,
4009  * const bool checkpoint)
4010  * {
4011  * pcout << "MainLoop<dim>::output(t = " << t
4012  * << ", checkpoint = " << checkpoint << ")" << std::endl;
4013  *
4014  * @endcode
4015  *
4016  * If the asynchronous writeback option is set we launch a background
4017  * thread performing all the slow IO to disc. In that case we have to
4018  * make sure that the background thread actually finished running. If
4019  * not, we have to wait to for it to finish. We launch said background
4020  * thread with <a
4021  * href="https://en.cppreference.com/w/cpp/thread/async"><code>std::async()</code></a>
4022  * that returns a <a
4023  * href="https://en.cppreference.com/w/cpp/thread/future"><code>std::future</code></a>
4024  * object. This <code>std::future</code> object contains the return
4025  * value of the function, which is in our case simply
4026  * <code>void</code>.
4027  *
4028 
4029  *
4030  *
4031  * @code
4032  * if (background_thread_state.valid())
4033  * {
4034  * TimerOutput::Scope timer(computing_timer, "main_loop - stalled output");
4035  * background_thread_state.wait();
4036  * }
4037  *
4038  * constexpr auto problem_dimension =
4039  * ProblemDescription<dim>::problem_dimension;
4040  *
4041  * @endcode
4042  *
4043  * At this point we make a copy of the state vector, run the schlieren
4044  * postprocessor, and run DataOut<dim>::build_patches() The actual
4045  * output code is standard: We create a DataOut instance, attach all
4046  * data vectors we want to output and call
4047  * DataOut<dim>::build_patches(). There is one twist, however. In order
4048  * to perform asynchronous IO on a background thread we create the
4049  * DataOut<dim> object as a shared pointer that we pass on to the
4050  * worker thread to ensure that once we exit this function and the
4051  * worker thread finishes the DataOut<dim> object gets destroyed again.
4052  *
4053 
4054  *
4055  *
4056  * @code
4057  * for (unsigned int i = 0; i < problem_dimension; ++i)
4058  * {
4059  * output_vector[i] = U[i];
4060  * output_vector[i].update_ghost_values();
4061  * }
4062  *
4063  * schlieren_postprocessor.compute_schlieren(output_vector);
4064  *
4065  * auto data_out = std::make_shared<DataOut<dim>>();
4066  *
4067  * data_out->attach_dof_handler(offline_data.dof_handler);
4068  *
4069  * const auto &component_names = ProblemDescription<dim>::component_names;
4070  *
4071  * for (unsigned int i = 0; i < problem_dimension; ++i)
4072  * data_out->add_data_vector(output_vector[i], component_names[i]);
4073  *
4074  * data_out->add_data_vector(schlieren_postprocessor.schlieren,
4075  * "schlieren_plot");
4076  *
4077  * data_out->build_patches(discretization.mapping,
4078  * discretization.finite_element.degree - 1);
4079  *
4080  * @endcode
4081  *
4082  * Next we create a lambda function for the background thread. We <a
4083  * href="https://en.cppreference.com/w/cpp/language/lambda">capture</a>
4084  * the <code>this</code> pointer as well as most of the arguments of
4085  * the output function by value so that we have access to them inside
4086  * the lambda function.
4087  *
4088  * @code
4089  * const auto output_worker = [this, name, t, cycle, checkpoint, data_out]() {
4090  * if (checkpoint)
4091  * {
4092  * @endcode
4093  *
4094  * We checkpoint the current state by doing the precise inverse
4095  * operation to what we discussed for the <a href="Resume">resume
4096  * logic</a>:
4097  *
4098 
4099  *
4100  *
4101  * @code
4102  * const unsigned int i =
4103  * discretization.triangulation.locally_owned_subdomain();
4104  * std::string filename =
4105  * name + "-checkpoint-" + Utilities::int_to_string(i, 4) + ".archive";
4106  *
4107  * std::ofstream file(filename, std::ios::binary | std::ios::trunc);
4108  *
4109  * boost::archive::binary_oarchive oa(file);
4110  * oa << t << cycle;
4111  * for (const auto &it1 : output_vector)
4112  * for (const auto &it2 : it1)
4113  * oa << it2;
4114  * }
4115  *
4116  * DataOutBase::VtkFlags flags(t,
4117  * cycle,
4118  * true,
4120  * data_out->set_flags(flags);
4121  *
4122  * data_out->write_vtu_with_pvtu_record(
4123  * "", name + "-solution", cycle, mpi_communicator, 6);
4124  * };
4125  *
4126  * @endcode
4127  *
4128  * If the asynchronous writeback option is set we launch a new
4129  * background thread with the help of
4130  * <a
4131  * href="https://en.cppreference.com/w/cpp/thread/async"><code>std::async</code></a>
4132  * function. The function returns a <a
4133  * href="https://en.cppreference.com/w/cpp/thread/future"><code>std::future</code></a>
4134  * object that we can use to query the status of the background thread.
4135  * At this point we can return from the <code>output()</code> function
4136  * and resume with the time stepping in the main loop - the thread will
4137  * run in the background.
4138  *
4139  * @code
4140  * if (asynchronous_writeback)
4141  * {
4142  * #ifdef DEAL_II_WITH_THREADS
4143  * background_thread_state = std::async(std::launch::async, output_worker);
4144  * #else
4145  * AssertThrow(
4146  * false,
4147  * ExcMessage(
4148  * "\"asynchronous_writeback\" was set to true but deal.II was built "
4149  * "without thread support (\"DEAL_II_WITH_THREADS=false\")."));
4150  * #endif
4151  * }
4152  * else
4153  * {
4154  * output_worker();
4155  * }
4156  * }
4157  *
4158  * } // namespace Step69
4159  *
4160  * @endcode
4161  *
4162  * And finally, the main function.
4163  *
4164 
4165  *
4166  *
4167  * @code
4168  * int main(int argc, char *argv[])
4169  * {
4170  * try
4171  * {
4172  * constexpr int dim = 2;
4173  *
4174  * using namespace dealii;
4175  * using namespace Step69;
4176  *
4177  * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
4178  *
4179  * MPI_Comm mpi_communicator(MPI_COMM_WORLD);
4180  * MainLoop<dim> main_loop(mpi_communicator);
4181  *
4182  * main_loop.run();
4183  * }
4184  * catch (std::exception &exc)
4185  * {
4186  * std::cerr << std::endl
4187  * << std::endl
4188  * << "----------------------------------------------------"
4189  * << std::endl;
4190  * std::cerr << "Exception on processing: " << std::endl
4191  * << exc.what() << std::endl
4192  * << "Aborting!" << std::endl
4193  * << "----------------------------------------------------"
4194  * << std::endl;
4195  * return 1;
4196  * }
4197  * catch (...)
4198  * {
4199  * std::cerr << std::endl
4200  * << std::endl
4201  * << "----------------------------------------------------"
4202  * << std::endl;
4203  * std::cerr << "Unknown exception!" << std::endl
4204  * << "Aborting!" << std::endl
4205  * << "----------------------------------------------------"
4206  * << std::endl;
4207  * return 1;
4208  * };
4209  * }
4210  * @endcode
4211 <a name="Results"></a>
4212 <a name="Results"></a><h1>Results</h1>
4213 
4214 
4215 Running the program with default parameters in release mode takes about 1
4216 minute on a 4 core machine (with hyperthreading):
4217 @verbatim
4218 # mpirun -np 4 ./step-69 | tee output
4219 Reading parameters and allocating objects... done
4220 
4221  ####################################################
4222  ######### #########
4223  ######### create triangulation #########
4224  ######### #########
4225  ####################################################
4226 
4227 Number of active cells: 36864
4228 
4229  ####################################################
4230  ######### #########
4231  ######### compute offline data #########
4232  ######### #########
4233  ####################################################
4234 
4235 Number of degrees of freedom: 37376
4236 
4237  ####################################################
4238  ######### #########
4239  ######### set up time step #########
4240  ######### #########
4241  ####################################################
4242 
4243  ####################################################
4244  ######### #########
4245  ######### interpolate initial values #########
4246  ######### #########
4247  ######### #########
4248  ####################################################
4249 
4250 TimeLoop<dim>::interpolate_initial_values(t = 0)
4251 TimeLoop<dim>::output(t = 0, checkpoint = 0)
4252 
4253  ####################################################
4254  ######### #########
4255  ######### enter main loop #########
4256  ######### #########
4257  ######### #########
4258  ####################################################
4259 
4260  ####################################################
4261  ######### #########
4262  ######### Cycle 000001 (0.0%) #########
4263  ######### at time t = 0.00000000 #########
4264  ######### #########
4265  ####################################################
4266 
4267 [...]
4268 
4269  ####################################################
4270  ######### #########
4271  ######### Cycle 007553 (100.0%) #########
4272  ######### at time t = 3.99984036 #########
4273  ######### #########
4274  ####################################################
4275 
4276 TimeLoop<dim>::output(t = 4.00038, checkpoint = 1)
4277 
4278 +------------------------------------------------------------------------+------------+------------+
4279 | Total CPU time elapsed since start | 357s | |
4280 | | | |
4281 | Section | no. calls | CPU time | % of total |
4282 +------------------------------------------------------------+-----------+------------+------------+
4283 | discretization - setup | 1 | 0.113s | 0% |
4284 | offline_data - assemble lumped mass matrix, and c_ij | 1 | 0.167s | 0% |
4285 | offline_data - compute |c_ij|, and n_ij | 1 | 0.00255s | 0% |
4286 | offline_data - create sparsity pattern and set up matrices | 1 | 0.0224s | 0% |
4287 | offline_data - distribute dofs | 1 | 0.0617s | 0% |
4288 | offline_data - fix slip boundary c_ij | 1 | 0.0329s | 0% |
4289 | schlieren_postprocessor - compute schlieren plot | 201 | 0.811s | 0.23% |
4290 | schlieren_postprocessor - prepare scratch space | 1 | 7.6e-05s | 0% |
4291 | time_loop - setup scratch space | 1 | 0.127s | 0% |
4292 | time_loop - stalled output | 200 | 0.000685s | 0% |
4293 | time_step - 1 compute d_ij | 7553 | 240s | 67% |
4294 | time_step - 2 compute d_ii, and tau_max | 7553 | 11.5s | 3.2% |
4295 | time_step - 3 perform update | 7553 | 101s | 28% |
4296 | time_step - 4 fix boundary states | 7553 | 0.724s | 0.2% |
4297 | time_step - prepare scratch space | 1 | 0.00245s | 0% |
4298 +------------------------------------------------------------+-----------+------------+------------+
4299 @endverbatim
4300 
4301 One thing that becomes evident is the fact that the program spends two
4302 thirds of the execution time computing the graph viscosity d_ij and about a
4303 third of the execution time in performing the update, where computing the
4304 flux @f$f(U)@f$ is the expensive operation. The preset default resolution is
4305 about 37k gridpoints, which amounts to about 148k spatial degrees of
4306 freedom in 2D. An animated schlieren plot of the solution looks as follows:
4307 
4308 <img src="https://www.dealii.org/images/steps/developer/step-69.coarse.gif" alt="" height="300">
4309 
4310 It is evident that 37k gridpoints for the first-order method is nowhere
4311 near the resolution needed to resolve any flow features. For comparison,
4312 here is a "reference" computation with a second-order method and about 9.5M
4313 gridpoints:
4314 
4315 <img src="https://www.dealii.org/images/steps/developer/step-69.2nd-order.t400.jpg" alt="" height="300">
4316 
4317 So, we give the first-order method a second chance and run it with about
4318 2.4M gridpoints on a small compute server:
4319 
4320 @verbatim
4321 # mpirun -np 16 ./step-69 | tee output
4322 
4323 [...]
4324 
4325  ####################################################
4326  ######### #########
4327  ######### Cycle 070216 (100.0%) #########
4328  ######### at time t = 3.99999231 #########
4329  ######### #########
4330  ####################################################
4331 
4332 TimeLoop<dim>::output(t = 4.00006, checkpoint = 1)
4333 
4334 [...]
4335 
4336 +------------------------------------------------------------------------+------------+------------+
4337 | Total wallclock time elapsed since start | 6.75e+03s | |
4338 | | | |
4339 | Section | no. calls | wall time | % of total |
4340 +------------------------------------------------------------+-----------+------------+------------+
4341 | discretization - setup | 1 | 1.97s | 0% |
4342 | offline_data - assemble lumped mass matrix, and c_ij | 1 | 1.19s | 0% |
4343 | offline_data - compute |c_ij|, and n_ij | 1 | 0.0172s | 0% |
4344 | offline_data - create sparsity pattern and set up matrices | 1 | 0.413s | 0% |
4345 | offline_data - distribute dofs | 1 | 1.05s | 0% |
4346 | offline_data - fix slip boundary c_ij | 1 | 0.252s | 0% |
4347 | schlieren_postprocessor - compute schlieren plot | 201 | 1.82s | 0% |
4348 | schlieren_postprocessor - prepare scratch space | 1 | 0.000497s | 0% |
4349 | time_loop - setup scratch space | 1 | 1.45s | 0% |
4350 | time_loop - stalled output | 200 | 0.00342s | 0% |
4351 | time_step - 1 compute d_ij | 70216 | 4.38e+03s | 65% |
4352 | time_step - 2 compute d_ii, and tau_max | 70216 | 419s | 6.2% |
4353 | time_step - 3 perform update | 70216 | 1.87e+03s | 28% |
4354 | time_step - 4 fix boundary states | 70216 | 24s | 0.36% |
4355 | time_step - prepare scratch space | 1 | 0.0227s | 0% |
4356 +------------------------------------------------------------+-----------+------------+------------+
4357 @endverbatim
4358 
4359 And with the following result:
4360 
4361 <img src="https://www.dealii.org/images/steps/developer/step-69.fine.gif" alt="" height="300">
4362 
4363 That's substantially better, although of course at the price of having run
4364 the code for roughly 2 hours on 16 cores.
4365  *
4366  *
4367 <a name="PlainProg"></a>
4368 <h1> The plain program</h1>
4369 @include "step-69.cc"
4370 */
SymmetricTensor::trace
constexpr Number trace(const SymmetricTensor< 2, dim, Number > &d)
Definition: symmetric_tensor.h:2705
Physics::Elasticity::Kinematics::F
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
std::exp
inline ::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5368
LAPACKSupport::diagonal
@ diagonal
Matrix is diagonal.
Definition: lapack_support.h:121
internal::QGaussLobatto::gamma
long double gamma(const unsigned int n)
Definition: quadrature_lib.cc:96
ParameterAcceptor::parse_parameters_call_back
boost::signals2::signal< void()> parse_parameters_call_back
Definition: parameter_acceptor.h:464
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition: fe_update_flags.h:122
TimerOutput::Scope
Definition: timer.h:554
TimerOutput::cpu_and_wall_times
@ cpu_and_wall_times
Definition: timer.h:653
Utilities::int_to_string
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:474
ScalarFunctionFromFunctionObject
Definition: function.h:727
dealii
Definition: namespace_dealii.h:25
LinearAlgebra::distributed::Vector< double >
TimeStepping::TimeStepping
Definition: time_stepping.h:96
Differentiation::SD::OptimizerType::lambda
@ lambda
WorkStream
Definition: work_stream.h:157
VectorizedArray::max
VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &x, const ::VectorizedArray< Number, width > &y)
Definition: vectorization.h:5466
parallel::transform
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Predicate &predicate, const unsigned int grainsize)
Definition: parallel.h:211
DataOutBase::vtk
@ vtk
Definition: data_out_base.h:1605
LAPACKSupport::V
static const char V
Definition: lapack_support.h:175
SparseMatrix
Definition: sparse_matrix.h:497
LAPACKSupport::U
static const char U
Definition: lapack_support.h:167
TimeStepping
Definition: time_stepping.h:36
GeometryInfo
Definition: geometry_info.h:1224
Physics::Elasticity::Kinematics::C
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
MPI_Comm
IteratorState::valid
@ valid
Iterator points to a valid object.
Definition: tria_iterator_base.h:38
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
update_values
@ update_values
Shape function values.
Definition: fe_update_flags.h:78
second
Point< 2 > second
Definition: grid_out.cc:4353
DEAL_II_ALWAYS_INLINE
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:99
DataOut::build_patches
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:1071
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
internal::p4est::functions
int(&) functions(const void *v1, const void *v2)
Definition: p4est_wrappers.cc:339
DoFHandler
Definition: dof_handler.h:205
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
ParameterAcceptor::add_parameter
void add_parameter(const std::string &entry, ParameterType &parameter, const std::string &documentation="", ParameterHandler &prm_=prm, const Patterns::PatternBase &pattern= *Patterns::Tools::Convert< ParameterType >::to_pattern())
Definition: parameter_acceptor.h:639
Utilities::CUDA::free
void free(T *&pointer)
Definition: cuda.h:96
WorkStream::run
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)
Definition: work_stream.h:1185
level
unsigned int level
Definition: grid_out.cc:4355
TimerOutput
Definition: timer.h:546
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
MeshWorker::ScratchData
Definition: scratch_data.h:232
DataOutBase::none
@ none
Definition: data_out_base.h:1557
VectorizedArray::exp
VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5368
SynchronousIterators::advance
void advance(std::tuple< I1, I2 > &t, const unsigned int n)
Definition: synchronous_iterator.h:146
LAPACKSupport::symmetric
@ symmetric
Matrix is symmetric.
Definition: lapack_support.h:115
Algorithms::Events::initial
const Event initial
Definition: event.cc:65
Tensor::norm
numbers::NumberTraits< Number >::real_type norm() const
ParameterAcceptor::initialize
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)
Definition: parameter_acceptor.cc:56
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
double
Tensor< 1, dim >
LocalIntegrators::Divergence::norm
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
update_gradients
@ update_gradients
Shape function gradients.
Definition: fe_update_flags.h:84
SymmetricTensor::sum
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
Physics::Elasticity::Kinematics::b
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:60
numbers::E
static constexpr double E
Definition: numbers.h:212
VectorTools::mean
@ mean
Definition: vector_tools_common.h:79
Threads::internal::call
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
Definition: thread_management.h:607
SparsityPattern
Definition: sparsity_pattern.h:865
LAPACKSupport::general
@ general
No special properties.
Definition: lapack_support.h:113
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:65
BlockIndices::swap
void swap(BlockIndices &u, BlockIndices &v)
Definition: block_indices.h:475
TimerOutput::never
@ never
Definition: timer.h:613
DataOutBase::VtkFlags
Definition: data_out_base.h:1095
LAPACKSupport::A
static const char A
Definition: lapack_support.h:155
unsigned int
value
static const bool value
Definition: dof_tools_constraints.cc:433
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:129
update_normal_vectors
@ update_normal_vectors
Normal vectors.
Definition: fe_update_flags.h:136
std::sqrt
inline ::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5412
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
Utilities::MPI::this_mpi_process
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
std_cxx20::ranges::iota_view
boost::integer_range< IncrementableType > iota_view
Definition: iota_view.h:48
VectorizedArray::sqrt
VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5412
internal::assemble
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition: loop.h:71
std::pow
inline ::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &x, const Number p)
Definition: vectorization.h:5428
ConditionalOStream
Definition: conditional_ostream.h:81
LAPACKSupport::zero
static const types::blas_int zero
Definition: lapack_support.h:179
Utilities::MPI::min
T min(const T &t, const MPI_Comm &mpi_communicator)
DataOutBase::VtkFlags::best_speed
@ best_speed
Definition: data_out_base.h:1142
Point< dim >
ParameterHandler
Definition: parameter_handler.h:845
FullMatrix< double >
FETools::interpolate
void interpolate(const DoFHandlerType1< dim, spacedim > &dof1, const InVector &u1, const DoFHandlerType2< dim, spacedim > &dof2, OutVector &u2)
VectorTools::interpolate
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask=ComponentMask())
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
std_cxx20
Definition: c++.h:183
MeshWorker::loop
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition: loop.h:443
DEAL_II_WITH_THREADS
#define DEAL_II_WITH_THREADS
Definition: config.h:64
first
Point< 2 > first
Definition: grid_out.cc:4352
ParameterAcceptor
Definition: parameter_acceptor.h:350
DataOut
Definition: data_out.h:148
make_iterator_range
IteratorRange< BaseIterator > make_iterator_range(const BaseIterator &begin, const typename identity< BaseIterator >::type &end)
Definition: iterator_range.h:296
gather
void gather(VectorizedArray< Number, width > &out, const std::array< Number *, width > &ptrs, const unsigned int offset)
Definition: vectorization.h:770
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
parallel
Definition: distributed.h:416
parallel::apply_to_subranges
void apply_to_subranges(const RangeType &begin, const typename identity< RangeType >::type &end, const Function &f, const unsigned int grainsize)
Definition: parallel.h:445
Utilities::MPI::max
T max(const T &t, const MPI_Comm &mpi_communicator)
internal::VectorOperations::copy
void copy(const T *begin, const T *end, U *dest)
Definition: vector_operations_internal.h:67
int
Utilities::MPI::MPI_InitFinalize
Definition: mpi.h:828