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