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