step-57.h

Go to the documentation of this file.

 1);
 *         flow.run(4);
 *       }
 *     catch (std::exception &exc)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Exception on processing: " << std::endl
 *                   << exc.what() << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *     catch (...)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Unknown exception!" << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *     return 0;
 *   }
 * @endcode
<a name="step_57-Results"></a><h1>Results</h1>
 
 
Now we use the method we discussed above to solve Navier Stokes equations with
viscosity @f$1/400@f$ and @f$1/7500@f$.
 
<a name="step_57-Testcase1LowReynoldsNumber"></a><h3> Test case 1: Low Reynolds Number </h3>
 
 
In the first test case the viscosity is set to be @f$1/400@f$. As we discussed in the
introduction, the initial guess is the solution to the corresponding Stokes
problem. In the following table, the residuals at each Newton's iteration on
every mesh is shown. The data in the table shows that Newton's iteration
converges quadratically.
 
<table align="center" class="doxtable">
<tr>
    <th>@f$\mathrm{Re}=400@f$</th>
    <th colspan="2">Mesh0</th>
    <th colspan="2">Mesh1</th>
    <th colspan="2">Mesh2</th>
    <th colspan="2">Mesh3</th>
    <th colspan="2">Mesh4</th>
</tr>
<tr>
    <th>Newton iter   </th>
    <th>Residual      </th>
    <th>FGMRES        </th>
    <th>Residual      </th>
    <th>FGMRES        </th>
    <th>Residual      </th>
    <th>FGMRES        </th>
    <th>Residual      </th>
    <th>FGMRES        </th>
    <th>Residual      </th>
    <th>FGMRES        </th>
</tr>
<tr>
  <td>1</td>
  <td>3.7112e-03</td>
  <td>5</td>
  <td>6.4189e-03</td>
  <td>3</td>
  <td>2.4338e-03</td>
  <td>3</td>
  <td>1.0570e-03</td>
  <td>3</td>
  <td>4.9499e-04</td>
  <td>3</td>
</tr>
<tr>
  <td>2</td>
  <td>7.0849e-04</td>
  <td>5</td>
  <td>9.9458e-04</td>
  <td>5</td>
  <td>1.1409e-04</td>
  <td>6</td>
  <td>1.3544e-05</td>
  <td>6</td>
  <td>1.4171e-06</td>
  <td>6</td>
</tr>
<tr>
  <td>3</td>
  <td>1.9980e-05</td>
  <td>5</td>
  <td>4.5007e-05</td>
  <td>5</td>
  <td>2.9020e-08</td>
  <td>5</td>
  <td>4.4021e-10</td>
  <td>6</td>
  <td>6.3435e-11</td>
  <td>6</td>
</tr>
<tr>
  <td>4</td>
  <td>2.3165e-09</td>
  <td>6</td>
  <td>1.6891e-07</td>
  <td>5</td>
  <td>1.2338e-14</td>
  <td>7</td>
  <td>1.8506e-14</td>
  <td>8</td>
  <td>8.8563e-15</td>
  <td>8</td>
</tr>
<tr>
  <td>5</td>
  <td>1.2585e-13</td>
  <td>7</td>
  <td>1.4520e-11</td>
  <td>6</td>
  <td>1.9044e-13</td>
  <td>8</td>
  <td></td>
  <td></td>
  <td></td>
  <td></td>
</tr>
<tr>
  <td>6</td>
  <td></td>
  <td></td>
  <td>1.3998e-15</td>
  <td>8</td>
  <td></td>
  <td></td>
  <td></td>
  <td></td>
  <td></td>
  <td></td>
</tr>
</table>
 
 
 
 
 
 
The following figures show the sequence of generated grids. For the case
of @f$\mathrm{Re}=400@f$, the initial guess is obtained by solving Stokes on an
@f$8 \times 8@f$ mesh, and the mesh is refined adaptively. Between meshes, the
solution from the coarse mesh is interpolated to the fine mesh to be used as an
initial guess.
 
<table align="center">
  <tr>
    <td align="center">
      <img src="https://www.dealii.org/images/steps/developer/step-57.Re400_Mesh0.png" width="232px" alt="">
    </td>
    <td align="center">
      <img src="https://www.dealii.org/images/steps/developer/step-57.Re400_Mesh1.png" width="232px" alt="">
    </td>
    <td align="center">
      <img src="https://www.dealii.org/images/steps/developer/step-57.Re400_Mesh2.png" width="232px" alt="">
    </td>
  </tr>
  <tr>
    <td align="center">
      <img src="https://www.dealii.org/images/steps/developer/step-57.Re400_Mesh3.png" width="232px" alt="">
    </td>
    <td align="center">
      <img src="https://www.dealii.org/images/steps/developer/step-57.Re400_Mesh4.png" width="232px" alt="">
    </td>
  </tr>
</table>
 
This picture is the graphical streamline result of lid-driven cavity with
@f$\mathrm{Re}=400@f$.
<img src="https://www.dealii.org/images/steps/developer/step-57.Re400_Streamline.png" alt="">
 
Then the solution is compared with a reference solution
from @cite Ghia1982 and the reference solution data can be found in the file "ref_2d_ghia_u.txt".
 
<img src="https://www.dealii.org/images/steps/developer/step-57.compare-Re400.svg" style="width:50%" alt="">
 
 
 
<a name="step_57-Testcase2HighReynoldsNumber"></a><h3> Test case 2: High Reynolds Number </h3>
 
 
Newton's iteration requires a good initial guess. However, the nonlinear term
dominates when the Reynolds number is large, so that the solution to the Stokes
equations may be far away from the exact solution. If the Stokes solution acts
as the initial guess, the convergence will be lost. The following picture
shows that the nonlinear iteration gets stuck and the residual no longer decreases
in further iterations.
 
<img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_loss_convergence.svg" style="width:50%" alt="">
 
The initial guess, therefore, has to be obtained via a continuation method
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
mesh is of size @f$32 \times 32@f$. After obtaining an initial guess, the mesh is
refined as in the previous test case. The following picture shows that at each
refinement Newton's iteration has quadratic convergence. 52 steps of Newton's
iterations are executed for solving this test case.
 
<img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_get_convergence.svg" style="width:50%" alt="">
 
We also show the residual from each step of Newton's iteration on every
mesh. The quadratic convergence is clearly visible in the table.
 
<table align="center" class="doxtable">
  <tr>
    <th>@f$\mathrm{Re}=7500@f$</th>
    <th colspan="2">Mesh0</th>
    <th colspan="2">Mesh1</th>
    <th colspan="2">Mesh2</th>
    <th colspan="2">Mesh3</th>
    <th colspan="2">Mesh4</th>
  </tr>
  <tr>
    <th>Newton iter   </th>
    <th>Residual      </th>
    <th>FGMRES        </th>
    <th>Residual      </th>
    <th>FGMRES        </th>
    <th>Residual      </th>
    <th>FGMRES        </th>
    <th>Residual      </th>
    <th>FGMRES        </th>
    <th>Residual      </th>
    <th>FGMRES        </th>
  </tr>
<tr>
  <td>1</td>
  <td>1.8922e-06</td>
  <td>6</td>
  <td>4.2506e-03</td>
  <td>3</td>
  <td>1.4299e-03</td>
  <td>3</td>
  <td>4.8793e-04</td>
  <td>2</td>
  <td>1.8998e-04</td>
  <td>2</td>
</tr>
<tr>
  <td>2</td>
  <td>3.1644e-09</td>
  <td>8</td>
  <td>1.3732e-03</td>
  <td>7</td>
  <td>4.1506e-04</td>
  <td>7</td>
  <td>9.1119e-05</td>
  <td>8</td>
  <td>1.3555e-05</td>
  <td>8</td>
</tr>
<tr>
  <td>3</td>
  <td>1.7611e-14</td>
  <td>9</td>
  <td>2.1946e-04</td>
  <td>6</td>
  <td>1.7881e-05</td>
  <td>6</td>
  <td>5.2678e-07</td>
  <td>7</td>
  <td>9.3739e-09</td>
  <td>7</td>
</tr>
<tr>
  <td>4</td>
  <td></td>
  <td></td>
  <td>8.8269e-06</td>
  <td>6</td>
  <td>6.8210e-09</td>
  <td>7</td>
  <td>2.2770e-11</td>
  <td>8</td>
  <td>1.2588e-13</td>
  <td>9</td>
</tr>
<tr>
  <td>5</td>
  <td></td>
  <td></td>
  <td>1.2974e-07</td>
  <td>7</td>
  <td>1.2515e-13</td>
  <td>9</td>
  <td>1.7801e-14</td>
  <td>1</td>
  <td></td>
  <td></td>
</tr>
<tr>
  <td>6</td>
  <td></td>
  <td></td>
  <td>4.4352e-11</td>
  <td>7</td>
  <td></td>
  <td></td>
  <td></td>
  <td></td>
  <td></td>
  <td></td>
</tr>
<tr>
  <td>7</td>
  <td></td>
  <td></td>
  <td>6.2863e-15</td>
  <td>9</td>
  <td></td>
  <td></td>
  <td></td>
  <td></td>
  <td></td>
  <td></td>
</tr>
</table>
 
 
 
 
 
 
The sequence of generated grids looks like this:
<table align="center">
  <tr>
    <td align="center">
      <img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_Mesh0.png" width="232px" alt="">
    </td>
    <td align="center">
      <img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_Mesh1.png" width="232px" alt="">
    </td>
    <td align="center">
      <img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_Mesh2.png" width="232px" alt="">
    </td>
  </tr>
  <tr>
    <td align="center">
      <img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_Mesh3.png" width="232px" alt="">
    </td>
    <td align="center">
      <img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_Mesh4.png" width="232px" alt="">
    </td>
  </tr>
</table>
We compare our solution with the reference solution from @cite Erturk2005 .
<img src="https://www.dealii.org/images/steps/developer/step-57.compare-Re7500.svg" style="width:50%" alt="">
The following picture presents the graphical result.
<img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_Streamline.png" alt="">
 
Furthermore, the error consists of the nonlinear error,
which decreases as we perform Newton iterations, and the discretization error,
which depends on the mesh size. That is why we have to refine the
mesh and repeat Newton's iteration on the next finer mesh. From the table above, we can
see that the final residual (nonlinear error) is below @f$10^{-12}@f$ on each mesh, but the
following picture shows us the difference between solutions on subsequently finer
meshes:
 
<img src="https://www.dealii.org/images/steps/developer/step-57.converge-Re7500.svg" style="width:50%" alt="">
 
 
<a name="step-57-extensions"></a>
 
<a name="step_57-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
 
 
<a name="step_57-Comparetoothersolvers"></a><h4>Compare to other solvers</h4>
 
 
It is easy to compare the currently implemented linear solver to just using
UMFPACK for the whole linear system. You need to remove the nullspace
containing the constant pressures and it is done in @ref step_56 "step-56". More interesting
is the comparison to other state of the art preconditioners like PCD. It turns
out that the preconditioner here is very competitive, as can be seen in the
paper @cite HeisterRapin2013.
 
The following table shows the timing results between our iterative approach
(FGMRES) compared to a direct solver (UMFPACK) for the whole system
with viscosity set to 1/400. Even though we use the same direct solver for
the velocity block in the iterative solver, it is considerably faster and
consumes less memory. This will be even more pronounced in 3d.
 
<table align="center" class="doxtable">
<tr>
  <th>Refinement Cycle</th>
  <th>DoFs</th>
  <th>Iterative: Total/s (Setup/s)</th>
  <th>Direct: Total/s (Setup/s)</th>
</tr>
<tr>
  <td>5</td>
  <td>9539</td>
  <td>0.10 (0.06)</td>
  <td>0.13 (0.12)</td>
</tr>
<tr>
  <td>6</td>
  <td>37507</td>
  <td>0.58 (0.37)</td>
  <td>1.03 (0.97)</td>
</tr>
<tr>
  <td>7</td>
  <td>148739</td>
  <td>3.59 (2.73)</td>
  <td>7.78 (7.53)</td>
</tr>
<tr>
  <td>8</td>
  <td>592387</td>
  <td>29.17 (24.94)</td>
  <td>(>4GB RAM)</td>
</tr>
</table>
 
 
<a name="step_57-3dcomputations"></a><h4>3d computations</h4>
 
 
The code is set up to also run in 3d. Of course the reference values are
different, see @cite Yang1998 for example. High resolution computations are not doable
with this example as is, because a direct solver for the velocity block does
not work well in 3d. Rather, a parallel solver based on algebraic or geometric
multigrid is needed -- see below.
 
<a name="step_57-Parallelization"></a><h4>Parallelization</h4>
 
 
For larger computations, especially in 3d, it is necessary to implement MPI
parallel solvers and preconditioners. A good starting point would be @ref step_55 "step-55",
which uses algebraic multigrid for the velocity block for the Stokes
equations. Another option would be to take a look at the list of codes
in the <a href="https://www.dealii.org/code-gallery.html">deal.II code
gallery</a>, which already contains parallel Navier-Stokes solvers.
 *
 *
<a name="step_57-PlainProg"></a>
<h1> The plain program</h1>
@include "step-57.cc"
*/