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