Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
step-57.h
Go to the documentation of this file.
1  1);
1429  * flow.run(4);
1430  * }
1431  * catch (std::exception &exc)
1432  * {
1433  * std::cerr << std::endl
1434  * << std::endl
1435  * << "----------------------------------------------------"
1436  * << std::endl;
1437  * std::cerr << "Exception on processing: " << std::endl
1438  * << exc.what() << std::endl
1439  * << "Aborting!" << std::endl
1440  * << "----------------------------------------------------"
1441  * << std::endl;
1442  * return 1;
1443  * }
1444  * catch (...)
1445  * {
1446  * std::cerr << std::endl
1447  * << std::endl
1448  * << "----------------------------------------------------"
1449  * << std::endl;
1450  * std::cerr << "Unknown exception!" << std::endl
1451  * << "Aborting!" << std::endl
1452  * << "----------------------------------------------------"
1453  * << std::endl;
1454  * return 1;
1455  * }
1456  * return 0;
1457  * }
1458  * @endcode
1459 <a name="Results"></a><h1>Results</h1>
1460 
1461 
1462 Now we use the method we discussed above to solve Navier Stokes equations with
1463 viscosity @f$1/400@f$ and @f$1/7500@f$.
1464 
1465 <a name="Testcase1LowReynoldsNumber"></a><h3> Test case 1: Low Reynolds Number </h3>
1466 
1467 
1468 In the first test case the viscosity is set to be @f$1/400@f$. As we discussed in the
1469 introduction, the initial guess is the solution to the corresponding Stokes
1470 problem. In the following table, the residuals at each Newton's iteration on
1471 every mesh is shown. The data in the table shows that Newton's iteration
1472 converges quadratically.
1473 
1474 <table align="center" class="doxtable">
1475 <tr>
1476  <th>@f$\mathrm{Re}=400@f$</th>
1477  <th colspan="2">Mesh0</th>
1478  <th colspan="2">Mesh1</th>
1479  <th colspan="2">Mesh2</th>
1480  <th colspan="2">Mesh3</th>
1481  <th colspan="2">Mesh4</th>
1482 </tr>
1483 <tr>
1484  <th>Newton iter </th>
1485  <th>Residual </th>
1486  <th>FGMRES </th>
1487  <th>Residual </th>
1488  <th>FGMRES </th>
1489  <th>Residual </th>
1490  <th>FGMRES </th>
1491  <th>Residual </th>
1492  <th>FGMRES </th>
1493  <th>Residual </th>
1494  <th>FGMRES </th>
1495 </tr>
1496 <tr>
1497  <td>1</td>
1498  <td>3.7112e-03</td>
1499  <td>5</td>
1500  <td>6.4189e-03</td>
1501  <td>3</td>
1502  <td>2.4338e-03</td>
1503  <td>3</td>
1504  <td>1.0570e-03</td>
1505  <td>3</td>
1506  <td>4.9499e-04</td>
1507  <td>3</td>
1508 </tr>
1509 <tr>
1510  <td>2</td>
1511  <td>7.0849e-04</td>
1512  <td>5</td>
1513  <td>9.9458e-04</td>
1514  <td>5</td>
1515  <td>1.1409e-04</td>
1516  <td>6</td>
1517  <td>1.3544e-05</td>
1518  <td>6</td>
1519  <td>1.4171e-06</td>
1520  <td>6</td>
1521 </tr>
1522 <tr>
1523  <td>3</td>
1524  <td>1.9980e-05</td>
1525  <td>5</td>
1526  <td>4.5007e-05</td>
1527  <td>5</td>
1528  <td>2.9020e-08</td>
1529  <td>5</td>
1530  <td>4.4021e-10</td>
1531  <td>6</td>
1532  <td>6.3435e-11</td>
1533  <td>6</td>
1534 </tr>
1535 <tr>
1536  <td>4</td>
1537  <td>2.3165e-09</td>
1538  <td>6</td>
1539  <td>1.6891e-07</td>
1540  <td>5</td>
1541  <td>1.2338e-14</td>
1542  <td>7</td>
1543  <td>1.8506e-14</td>
1544  <td>8</td>
1545  <td>8.8563e-15</td>
1546  <td>8</td>
1547 </tr>
1548 <tr>
1549  <td>5</td>
1550  <td>1.2585e-13</td>
1551  <td>7</td>
1552  <td>1.4520e-11</td>
1553  <td>6</td>
1554  <td>1.9044e-13</td>
1555  <td>8</td>
1556  <td></td>
1557  <td></td>
1558  <td></td>
1559  <td></td>
1560 </tr>
1561 <tr>
1562  <td>6</td>
1563  <td></td>
1564  <td></td>
1565  <td>1.3998e-15</td>
1566  <td>8</td>
1567  <td></td>
1568  <td></td>
1569  <td></td>
1570  <td></td>
1571  <td></td>
1572  <td></td>
1573 </tr>
1574 </table>
1575 
1576 
1577 
1578 
1579 
1580 
1581 The following figures show the sequence of generated grids. For the case
1582 of @f$\mathrm{Re}=400@f$, the initial guess is obtained by solving Stokes on an
1583 @f$8 \times 8@f$ mesh, and the mesh is refined adaptively. Between meshes, the
1584 solution from the coarse mesh is interpolated to the fine mesh to be used as an
1585 initial guess.
1586 
1587 <table align="center">
1588  <tr>
1589  <td align="center">
1590  <img src="https://www.dealii.org/images/steps/developer/step-57.Re400_Mesh0.png" width="232px" alt="">
1591  </td>
1592  <td align="center">
1593  <img src="https://www.dealii.org/images/steps/developer/step-57.Re400_Mesh1.png" width="232px" alt="">
1594  </td>
1595  <td align="center">
1596  <img src="https://www.dealii.org/images/steps/developer/step-57.Re400_Mesh2.png" width="232px" alt="">
1597  </td>
1598  </tr>
1599  <tr>
1600  <td align="center">
1601  <img src="https://www.dealii.org/images/steps/developer/step-57.Re400_Mesh3.png" width="232px" alt="">
1602  </td>
1603  <td align="center">
1604  <img src="https://www.dealii.org/images/steps/developer/step-57.Re400_Mesh4.png" width="232px" alt="">
1605  </td>
1606  </tr>
1607 </table>
1608 
1609 This picture is the graphical streamline result of lid-driven cavity with
1610 @f$\mathrm{Re}=400@f$.
1611 <img src="https://www.dealii.org/images/steps/developer/step-57.Re400_Streamline.png" alt="">
1612 
1613 Then the solution is compared with a reference solution
1614 from [4] and the reference solution data can be found in the file "ref_2d_ghia_u.txt".
1615 
1616 <img src="https://www.dealii.org/images/steps/developer/step-57.compare-Re400.svg" style="width:50%" alt="">
1617 
1618 <a name="Testcase2HighReynoldsNumber"></a><h3> Test case 2: High Reynolds Number </h3>
1619 
1620 
1621 Newton's iteration requires a good initial guess. However, the nonlinear term
1622 dominates when the Reynolds number is large, so that the solution to the Stokes
1623 equations may be far away from the exact solution. If the Stokes solution acts
1624 as the initial guess, the convergence will be lost. The following picture
1625 shows that the nonlinear iteration gets stuck and the residual no longer decreases
1626 in further iterations.
1627 
1628 <img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_loss_convergence.svg" style="width:50%" alt="">
1629 
1630 The initial guess, therefore, has to be obtained via a continuation method
1631 which has been discussed in the introduction. Here the step size in the continuation method, that is @f$|\nu_{i}-\nu_{i+1}|@f$, is 2000 and the initial
1632 mesh is of size @f$32 \times 32@f$. After obtaining an initial guess, the mesh is
1633 refined as in the previous test case. The following picture shows that at each
1634 refinement Newton's iteration has quadratic convergence. 52 steps of Newton's
1635 iterations are executed for solving this test case.
1636 
1637 <img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_get_convergence.svg" style="width:50%" alt="">
1638 
1639 We also show the residual from each step of Newton's iteration on every
1640 mesh. The quadratic convergence is clearly visible in the table.
1641 
1642 <table align="center" class="doxtable">
1643  <tr>
1644  <th>@f$\mathrm{Re}=7500@f$</th>
1645  <th colspan="2">Mesh0</th>
1646  <th colspan="2">Mesh1</th>
1647  <th colspan="2">Mesh2</th>
1648  <th colspan="2">Mesh3</th>
1649  <th colspan="2">Mesh4</th>
1650  </tr>
1651  <tr>
1652  <th>Newton iter </th>
1653  <th>Residual </th>
1654  <th>FGMRES </th>
1655  <th>Residual </th>
1656  <th>FGMRES </th>
1657  <th>Residual </th>
1658  <th>FGMRES </th>
1659  <th>Residual </th>
1660  <th>FGMRES </th>
1661  <th>Residual </th>
1662  <th>FGMRES </th>
1663  </tr>
1664 <tr>
1665  <td>1</td>
1666  <td>1.8922e-06</td>
1667  <td>6</td>
1668  <td>4.2506e-03</td>
1669  <td>3</td>
1670  <td>1.4299e-03</td>
1671  <td>3</td>
1672  <td>4.8793e-04</td>
1673  <td>2</td>
1674  <td>1.8998e-04</td>
1675  <td>2</td>
1676 </tr>
1677 <tr>
1678  <td>2</td>
1679  <td>3.1644e-09</td>
1680  <td>8</td>
1681  <td>1.3732e-03</td>
1682  <td>7</td>
1683  <td>4.1506e-04</td>
1684  <td>7</td>
1685  <td>9.1119e-05</td>
1686  <td>8</td>
1687  <td>1.3555e-05</td>
1688  <td>8</td>
1689 </tr>
1690 <tr>
1691  <td>3</td>
1692  <td>1.7611e-14</td>
1693  <td>9</td>
1694  <td>2.1946e-04</td>
1695  <td>6</td>
1696  <td>1.7881e-05</td>
1697  <td>6</td>
1698  <td>5.2678e-07</td>
1699  <td>7</td>
1700  <td>9.3739e-09</td>
1701  <td>7</td>
1702 </tr>
1703 <tr>
1704  <td>4</td>
1705  <td></td>
1706  <td></td>
1707  <td>8.8269e-06</td>
1708  <td>6</td>
1709  <td>6.8210e-09</td>
1710  <td>7</td>
1711  <td>2.2770e-11</td>
1712  <td>8</td>
1713  <td>1.2588e-13</td>
1714  <td>9</td>
1715 </tr>
1716 <tr>
1717  <td>5</td>
1718  <td></td>
1719  <td></td>
1720  <td>1.2974e-07</td>
1721  <td>7</td>
1722  <td>1.2515e-13</td>
1723  <td>9</td>
1724  <td>1.7801e-14</td>
1725  <td>1</td>
1726  <td></td>
1727  <td></td>
1728 </tr>
1729 <tr>
1730  <td>6</td>
1731  <td></td>
1732  <td></td>
1733  <td>4.4352e-11</td>
1734  <td>7</td>
1735  <td></td>
1736  <td></td>
1737  <td></td>
1738  <td></td>
1739  <td></td>
1740  <td></td>
1741 </tr>
1742 <tr>
1743  <td>7</td>
1744  <td></td>
1745  <td></td>
1746  <td>6.2863e-15</td>
1747  <td>9</td>
1748  <td></td>
1749  <td></td>
1750  <td></td>
1751  <td></td>
1752  <td></td>
1753  <td></td>
1754 </tr>
1755 </table>
1756 
1757 
1758 
1759 
1760 
1761 
1762 The sequence of generated grids looks like this:
1763 <table align="center">
1764  <tr>
1765  <td align="center">
1766  <img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_Mesh0.png" width="232px" alt="">
1767  </td>
1768  <td align="center">
1769  <img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_Mesh1.png" width="232px" alt="">
1770  </td>
1771  <td align="center">
1772  <img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_Mesh2.png" width="232px" alt="">
1773  </td>
1774  </tr>
1775  <tr>
1776  <td align="center">
1777  <img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_Mesh3.png" width="232px" alt="">
1778  </td>
1779  <td align="center">
1780  <img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_Mesh4.png" width="232px" alt="">
1781  </td>
1782  </tr>
1783 </table>
1784 We compare our solution with reference solution from [5].
1785 <img src="https://www.dealii.org/images/steps/developer/step-57.compare-Re7500.svg" style="width:50%" alt="">
1786 The following picture presents the graphical result.
1787 <img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_Streamline.png" alt="">
1788 
1789 Furthermore, the error consists of the nonlinear error,
1790 which decreases as we perform Newton iterations, and the discretization error,
1791 which depends on the mesh size. That is why we have to refine the
1792 mesh and repeat Newton's iteration on the next finer mesh. From the table above, we can
1793 see that the residual (nonlinear error) is below 1e-12 on each mesh, but the
1794 following picture shows us the difference between solutions on subsequently finer
1795 meshes.
1796 
1797 <img src="https://www.dealii.org/images/steps/developer/step-57.converge-Re7500.svg" style="width:50%" alt="">
1798 
1799 <a name="extensions"></a>
1800 
1801 <a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1802 
1803 
1804 <a name="Comparetoothersolvers"></a><h4>Compare to other solvers</h4>
1805 
1806 
1807 It is easy to compare the currently implemented linear solver to just using
1808 UMFPACK for the whole linear system. You need to remove the nullspace
1809 containing the constant pressures and it is done in @ref step_56 "step-56". More interesting
1810 is the comparison to other state of the art preconditioners like PCD. It turns
1811 out that the preconditioner here is very competitive, as can be seen in the
1812 paper [2].
1813 
1814 The following table shows the timing results between our iterative approach
1815 (FGMRES) compared to a direct solver (UMFPACK) for the whole system
1816 with viscosity set to 1/400. Even though we use the same direct solver for
1817 the velocity block in the iterative solver, it is considerably faster and
1818 consumes less memory. This will be even more pronounced in 3d.
1819 
1820 <table align="center" class="doxtable">
1821 <tr>
1822  <th>Refinement Cycle</th>
1823  <th>DoFs</th>
1824  <th>Iterative: Total/s (Setup/s)</th>
1825  <th>Direct: Total/s (Setup/s)</th>
1826 </tr>
1827 <tr>
1828  <td>5</td>
1829  <td>9539</td>
1830  <td>0.10 (0.06)</td>
1831  <td>0.13 (0.12)</td>
1832 </tr>
1833 <tr>
1834  <td>6</td>
1835  <td>37507</td>
1836  <td>0.58 (0.37)</td>
1837  <td>1.03 (0.97)</td>
1838 </tr>
1839 <tr>
1840  <td>7</td>
1841  <td>148739</td>
1842  <td>3.59 (2.73)</td>
1843  <td>7.78 (7.53)</td>
1844 </tr>
1845 <tr>
1846  <td>8</td>
1847  <td>592387</td>
1848  <td>29.17 (24.94)</td>
1849  <td>(>4GB RAM)</td>
1850 </tr>
1851 </table>
1852 
1853 
1854 <a name="3dcomputations"></a><h4>3d computations</h4>
1855 
1856 
1857 The code is set up to also run in 3d. Of course the reference values are
1858 different, see [6] for example. High resolution computations are not doable
1859 with this example as is, because a direct solver for the velocity block does
1860 not work well in 3d. Rather, a parallel solver based on algebraic or geometric
1861 multigrid is needed. See below.
1862 
1863 <a name="Parallelization"></a><h4>Parallelization</h4>
1864 
1865 
1866 For larger computations, especially in 3d, it is necessary to implement MPI
1867 parallel solvers and preconditioners. A good starting point would be @ref step_55 "step-55",
1868 which uses algebraic multigrid for the velocity block for the Stokes
1869 equations. Another option would be to take a look at the list of codes
1870 in the <a href="https://www.dealii.org/code-gallery.html">deal.II code
1871 gallery</a>, which already contains parallel Navier-Stokes solvers.
1872  *
1873  *
1874 <a name="PlainProg"></a>
1875 <h1> The plain program</h1>
1876 @include "step-57.cc"
1877 */
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
Algorithms::Events::initial
const Event initial
Definition: event.cc:65
first
Point< 2 > first
Definition: grid_out.cc:4352
GridRefinement::refine
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
Definition: grid_refinement.cc:41