From 776f42be4ea941361dadae0f439e9dea99d8f345 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Tue, 17 Jan 2017 11:04:55 -0500 Subject: [PATCH] address some left-over comments --- examples/step-57/doc/intro.dox | 10 +- examples/step-57/doc/results.dox | 362 +++++++++++++++---------------- examples/step-57/step-57.cc | 11 +- 3 files changed, 189 insertions(+), 194 deletions(-) diff --git a/examples/step-57/doc/intro.dox b/examples/step-57/doc/intro.dox index be97c6f08d..396f4eac5a 100644 --- a/examples/step-57/doc/intro.dox +++ b/examples/step-57/doc/intro.dox @@ -125,9 +125,11 @@ Now, Newton's iteration can be used to solve for the update terms:
  • Initialization: Initial guess $u_0$ and $p_0$, tolerance $\tau$;
  • Linear solve to compute update term $\delta\textbf{u}^{k}$ and $\delta p^k$;
  • Update the approximation: $\textbf{u}^{k+1} = \textbf{u}^{k} + \delta\textbf{u}^{k}$ and $p^{k+1} = p^{k} + \delta p^{k}$; -
  • Check residual norm: $E^{k+1} = \|F(\mathbf{u}^{k+1}, p^{k+1})\|$. - If $E^{k+1} \leq \tau$, STOP. - If $E^{k+1} > \tau$, back to step 2. +
  • Check residual norm: $E^{k+1} = \|F(\mathbf{u}^{k+1}, p^{k+1})\|$: +

    Finding an Initial Guess

    @@ -236,7 +238,7 @@ Instead of solving the above system, we can solve the equivalent system \right) @f} with a parameter $\gamma$ and an invertible matrix W. Here -$\gamma B^TW^{-1}B & B^{T}$ is the Augmented Lagrangian term and +$\gamma B^TW^{-1}B$ is the Augmented Lagrangian term and see [1] for details. Denoting the system matrix of the new system by $G$ and the right-hand diff --git a/examples/step-57/doc/results.dox b/examples/step-57/doc/results.dox index f1475d646d..86d9c8c7a4 100644 --- a/examples/step-57/doc/results.dox +++ b/examples/step-57/doc/results.dox @@ -20,7 +20,7 @@ every mesh is shown. The data in the table shows that Newton's iteration converg Mesh4 - Newton's iter + Newton iter Residual FGMRES Residual @@ -33,112 +33,112 @@ every mesh is shown. The data in the table shows that Newton's iteration converg FGMRES - 1 - 7.40396e-3 - 3 - 1.05562e-3 - 3 - 4.94796e-4 - 3 - 2.5624e-4 - 2 - 1.26733e-4 - 2 + 1 + 7.40396e-3 + 3 + 1.05562e-3 + 3 + 4.94796e-4 + 3 + 2.5624e-4 + 2 + 1.26733e-4 + 2 - 2 - 3.86766e-3 - 4 - 1.3549e-5 - 3 - 1.41981e-6 - 3 - 1.29108e-6 - 4 - 6.14794e-7 - 4 + 2 + 3.86766e-3 + 4 + 1.3549e-5 + 3 + 1.41981e-6 + 3 + 1.29108e-6 + 4 + 6.14794e-7 + 4 - 3 - 1.60421e-3 - 4 - 1.24836e-9 - 3 - 9.11557e-11 - 3 - 3.35933e-11 - 3 - 5.86734e-11 - 2 + 3 + 1.60421e-3 + 4 + 1.24836e-9 + 3 + 9.11557e-11 + 3 + 3.35933e-11 + 3 + 5.86734e-11 + 2 - 4 - 9.26748e-4 - 4 - 2.75537e-14 - 4 - 1.39986e-14 - 5 - 2.18864e-14 - 5 - 3.38787e-14 - 5 + 4 + 9.26748e-4 + 4 + 2.75537e-14 + 4 + 1.39986e-14 + 5 + 2.18864e-14 + 5 + 3.38787e-14 + 5 - 5 - 1.34601e-5 - 4 -   -   -   -   -   -   -   -   + 5 + 1.34601e-5 + 4 +   +   +   +   +   +   +   +   - 6 - 2.5235e-8 - 5 -   -   -   -   -   -   -   -   + 6 + 2.5235e-8 + 5 +   +   +   +   +   +   +   +   - 7 - 1.38899e-12 - 4 -   -   -   -   -   -   -   -   + 7 + 1.38899e-12 + 4 +   +   +   +   +   +   +   +   - 8 - 4.68224e-15 - 4 -   -   -   -   -   -   -   -   + 8 + 4.68224e-15 + 4 +   +   +   +   +   +   +   +   -The following figures show the sequence of the generated grids. For the case +The following figures show the sequence of generated grids. For the case of Re=400, the initial guess is obtained by solving Stokes on an $8 \times 8$ 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. @@ -193,11 +193,10 @@ iterations are executed for solving this test case. -Also we show the residual from each step of Newton's iteration on every -mesh. The quadratic convergence is shown clearly in the table. +We also show the residual from each step of Newton's iteration on every +mesh. The quadratic convergence is clearly visible in the table. - @@ -206,9 +205,8 @@ mesh. The quadratic convergence is shown clearly in the table. - - + @@ -220,103 +218,96 @@ mesh. The quadratic convergence is shown clearly in the table. - - - - - - - - - - - - + + + + + + + + + + + - - - - - - - - - - - - + + + + + + + + + + + - - - - - - - - - - - - + + + + + + + + + + + - - - - - - - - - - - - + + + + + + + + + + + - - - - - - - - - - - - + + + + + + + + + + + - - - - - - - - - - - - + + + + + + + + + + + - - - - - - - - - - - - + + + + + + + + + + +
      Mesh0Mesh3 Mesh4
    Newton's iter Newton iter Residual FGMRES Residual Residual FGMRES
    1 1.89223e-6 6 4.2506e-3 3 1.42993e-3 3 4.87932e-4 2 1.89981e-04 2 1 1.89223e-6 6 4.2506e-3 3 1.42993e-3 3 4.87932e-4 2 1.89981e-04 2
    2 3.16439e-98 1.3732e-3 7 4.15062e-4 7 9.11191e-5 8 1.35553e-58 2 3.16439e-98 1.3732e-3 7 4.15062e-4 7 9.11191e-5 8 1.35553e-58
    3 1.7628e-149 2.19455e-4 6 1.78805e-5 6 5.26782e-7 7 9.37391e-9 7 3 1.7628e-149 2.19455e-4 6 1.78805e-5 6 5.26782e-7 7 9.37391e-9 7
    4     8.82693e-6 6 6.82096e-9 7 2.27696e-11 8 1.25899e-139 4     8.82693e-6 6 6.82096e-9 7 2.27696e-11 8 1.25899e-139
    5     1.29739e-77 1.25167e-13 9 1.76128e-14 10     5     1.29739e-77 1.25167e-13 9 1.76128e-14 10    
    6     4.43518e-117             6     4.43518e-117            
    7     6.42323e-15 9            7     6.42323e-15 9           
    The sequence of generated grids looks like this: @@ -377,40 +368,35 @@ consumes less memory. This will be even more pronounced in 3d. - - - - + + + + - - - - -
    RefDoFsIterative: Total/s (Setup/s)Direct: Total/s (Setup/s)RefDoFsIterative: Total/s (Setup/s)Direct: Total/s (Setup/s)
    5 9539 0.10 (0.06) 0.13 (0.12)
    6 37507 0.58 (0.37) 1.03 (0.97)
    7 148739 3.59 (2.73) 7.78 (7.53)
    8 592387 29.17 (24.94) (>4GB RAM)
    diff --git a/examples/step-57/step-57.cc b/examples/step-57/step-57.cc index 468df2e0f2..0760f493cf 100644 --- a/examples/step-57/step-57.cc +++ b/examples/step-57/step-57.cc @@ -189,7 +189,7 @@ namespace Step57 // Schur complement preconditioner is defined in this part. As discussed in // the introduction, the preconditioner in Krylov iterative methods is // implemented as a matrix-vector product operator. In practice, the Schur - // complement preconditioner is decomposed as a product of three matrices(as + // complement preconditioner is decomposed as a product of three matrices (as // presented in the first section). The $\tilde{A}^{-1}$ in the first factor // involves a solve for the linear system $\tilde{A}x=b$. Here we solve // this system via a direct solver for simplicity. The computation involved @@ -516,6 +516,13 @@ namespace Step57 { pressure_mass_matrix.reinit(sparsity_pattern.block(1,1)); pressure_mass_matrix.copy_from(system_matrix.block(1,1)); + + // Note that settings this pressure block to zero is not identical to + // not assembling anything in this block, because this operation here + // will (incorrectly) delete diagonal entries that come in from + // hanging node constraints for pressure DoFs. This means that our + // whole system matrix will have rows that are completely + // zero. Luckily, FGMRES handles these rows without any problem. system_matrix.block(1,1) = 0; } } @@ -545,7 +552,7 @@ namespace Step57 { const ConstraintMatrix &constraints_used = initial_step ? nonzero_constraints : zero_constraints; - SolverControl solver_control (system_matrix.m(),1e-4*system_rhs.l2_norm(), true); + SolverControl solver_control (system_matrix.m(), 1e-4*system_rhs.l2_norm(), true); SolverFGMRES > gmres(solver_control); SparseILU pmass_preconditioner; -- 2.39.5