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