update step-57

- always assemble rhs when assembling matrix (so we get correct boundary
values from the ConstraintMatrix)
- this fixes having 0 iterations in the first Newton iteration
- restructure assemble* functions to make this more clear
- reorder some comments
- redo results section (neater tables)

diff --git a/examples/step-57/doc/results.dox b/examples/step-57/doc/results.dox
index 86d9c8c7a408a4a5539b44e16c6a0d5fdc5a014f..42ac92872e80083ce607685759d5d1890e788129 100644 (file)
--- a/examples/step-57/doc/results.dox
+++ b/examples/step-57/doc/results.dox
@@ -10,16 +10,16 @@ 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" border="1">
-  <tr>
-    <th>&nbsp;</th>
+<table align="center" class="doxtable">
+<tr>
+    <th>Re=400</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>
+</tr>
+<tr>
     <th>Newton iter   </th>
     <th>Residual      </th>
     <th>FGMRES        </th>
@@ -31,113 +31,92 @@ every mesh is shown. The data in the table shows that Newton's iteration converg
     <th>FGMRES        </th>
     <th>Residual      </th>
     <th>FGMRES        </th>
-  </tr>
-  <tr>
-    <td>1         </td>
-    <td>7.40396e-3</td>
-    <td>3         </td>
-    <td>1.05562e-3  </td>
-    <td>3         </td>
-    <td>4.94796e-4  </td>
-    <td>3         </td>
-    <td>2.5624e-4 </td>
-    <td>2         </td>
-    <td>1.26733e-4  </td>
-    <td>2         </td>
-  </tr>
-  <tr>
-    <td>2         </td>
-    <td>3.86766e-3  </td>
-    <td>4         </td>
-    <td>1.3549e-5  </td>
-    <td>3         </td>
-    <td>1.41981e-6  </td>
-    <td>3        </td>
-    <td>1.29108e-6  </td>
-    <td>4         </td>
-    <td>6.14794e-7  </td>
-    <td>4         </td>
-  </tr>
-  <tr>
-    <td>3         </td>
-    <td>1.60421e-3</td>
-    <td>4         </td>
-    <td>1.24836e-9  </td>
-    <td>3         </td>
-    <td>9.11557e-11  </td>
-    <td>3         </td>
-    <td>3.35933e-11 </td>
-    <td>3         </td>
-    <td>5.86734e-11 </td>
-    <td>2         </td>
-  </tr>
-  <tr>
-    <td>4         </td>
-    <td>9.26748e-4  </td>
-    <td>4         </td>
-    <td>2.75537e-14  </td>
-    <td>4         </td>
-    <td>1.39986e-14 </td>
-    <td>5         </td>
-    <td>2.18864e-14 </td>
-    <td>5         </td>
-    <td>3.38787e-14 </td>
-    <td>5         </td>
-  </tr>
-  <tr>
-    <td>5          </td>
-    <td>1.34601e-5</td>
-    <td>4          </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-  </tr>
-  <tr>
-    <td>6         </td>
-    <td>2.5235e-8    </td>
-    <td>5    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-  </tr>
-  <tr>
-    <td>7         </td>
-    <td>1.38899e-12    </td>
-    <td>4    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-  </tr>
-  <tr>
-    <td>8           </td>
-    <td>4.68224e-15 </td>
-    <td>4           </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;      </td>
-    <td>&nbsp;      </td>
-    <td>&nbsp;      </td>
-    <td>&nbsp;      </td>
-    <td>&nbsp;      </td>
-    <td>&nbsp;      </td>
-  </tr>
+</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.0000e+00</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.0000e+00</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.0000e+00</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.0000e+00</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 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
@@ -196,9 +175,9 @@ iterations are executed for solving this test case.
 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" border="1">
+<table align="center" class="doxtable">
   <tr>
-    <th>&nbsp;</th>
+    <th>Re=7500</th>
     <th colspan="2">Mesh0</th>
     <th colspan="2">Mesh1</th>
     <th colspan="2">Mesh2</th>
@@ -218,98 +197,104 @@ mesh. The quadratic convergence is clearly visible in the table.
     <th>Residual      </th>
     <th>FGMRES        </th>
   </tr>
-  <tr>
-    <td>1          </td>
-    <td>1.89223e-6  </td>
-    <td>6          </td>
-    <td>4.2506e-3  </td>
-    <td>3          </td>
-    <td>1.42993e-3  </td>
-    <td>3          </td>
-    <td>4.87932e-4  </td>
-    <td>2          </td>
-    <td>1.89981e-04  </td>
-    <td>2         </td>
-  </tr>
-  <tr>
-    <td>2         </td>
-    <td>3.16439e-9</td>
-    <td>8         </td>
-    <td>1.3732e-3 </td>
-    <td>7         </td>
-    <td>4.15062e-4 </td>
-    <td>7         </td>
-    <td>9.11191e-5 </td>
-    <td>8         </td>
-    <td>1.35553e-5</td>
-    <td>8         </td>
-  </tr>
-  <tr>
-    <td>3         </td>
-    <td>1.7628e-14</td>
-    <td>9         </td>
-    <td>2.19455e-4 </td>
-    <td>6         </td>
-    <td>1.78805e-5 </td>
-    <td>6         </td>
-    <td>5.26782e-7 </td>
-    <td>7         </td>
-    <td>9.37391e-9 </td>
-    <td>7         </td>
-  </tr>
-  <tr>
-    <td>4         </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>8.82693e-6 </td>
-    <td>6         </td>
-    <td>6.82096e-9 </td>
-    <td>7         </td>
-    <td>2.27696e-11 </td>
-    <td>8         </td>
-    <td>1.25899e-13</td>
-    <td>9        </td>
-  </tr>
-  <tr>
-    <td>5         </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>1.29739e-7</td>
-    <td>7         </td>
-    <td>1.25167e-13 </td>
-    <td>9        </td>
-    <td>1.76128e-14 </td>
-    <td>10         </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-  </tr>
-  <tr>
-    <td>6         </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>4.43518e-11</td>
-    <td>7         </td>
-    <td>&nbsp; </td>
-    <td>&nbsp;        </td>
-    <td>&nbsp; </td>
-    <td>&nbsp;        </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-  </tr>
-  <tr>
-    <td>7         </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-    <td>6.42323e-15 </td>
-    <td>9         </td>
-    <td>&nbsp;</td>
-    <td>&nbsp;         </td>
-    <td>&nbsp; </td>
-    <td>&nbsp;         </td>
-    <td>&nbsp;    </td>
-    <td>&nbsp;    </td>
-  </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>

diff --git a/examples/step-57/step-57.cc b/examples/step-57/step-57.cc
index 0760f493cf62fa8d1966f66e567e689553b33195..a54308a99f8391b574a8c3661e49bd4981d2818e 100644 (file)
--- a/examples/step-57/step-57.cc
+++ b/examples/step-57/step-57.cc
@@ -98,10 +98,9 @@ namespace Step57
   private:
     void setup_dofs();
     void initialize_system();
-    void assemble_system(const bool initial_step,
-                         const bool assemble_matrix,
-                         const bool assemble_rhs);
-    void assemble_matrix(const bool initial_step);
+    void assemble(const bool initial_step,
+                  const bool assemble_matrix);
+    void assemble_system(const bool initial_step);
     void assemble_rhs(const bool initial_step);
     void solve(const bool initial_step);
     void refine_mesh();
@@ -373,24 +372,23 @@ namespace Step57
     system_rhs.reinit (dofs_per_block);
   }
 
-  // @sect4{StationaryNavierStokes::assemble_system}
+  // @sect4{StationaryNavierStokes::assemble}
 
   // This function builds the system matrix and right hand side that we
-  // actually work on. "initial_step" is given for applying different
-  // constraints (nonzero for the initial step and zero for the others). The
-  // other two flags are to determine whether to assemble the system matrix
-  // or the right hand side vector, respectively.
+  // currently work on. The @p initial_step argument is used to determine
+  // which set of constraints we apply (nonzero for the initial step and zero
+  // for the others). The @p assemble_matrix flag determines whether to
+  // assemble the whole system or only the right hand side vector,
+  // respectively.
 
   template <int dim>
-  void StationaryNavierStokes<dim>::assemble_system(const bool initial_step,
-                                                    const bool assemble_matrix,
-                                                    const bool assemble_rhs)
+  void StationaryNavierStokes<dim>::assemble(const bool initial_step,
+                                             const bool assemble_matrix)
   {
     if (assemble_matrix)
       system_matrix    = 0;
 
-    if (assemble_rhs)
-      system_rhs       = 0;
+    system_rhs       = 0;
 
     QGauss<dim>   quadrature_formula(degree+2);
 
@@ -445,12 +443,13 @@ namespace Step57
         fe_values[pressure].get_function_values(evaluation_point,
                                                 present_pressure_values);
 
-        // The assembly is similar to step-22. An additional term with gamma as a coefficient
-        // is the Augmented Lagrangian (AL), which is assembled via grad-div stabilization.
-        // As we discussed in the introduction, the bottom right block of the system matrix should be
-        // zero. Since the pressure mass matrix is used while creating the preconditioner,
-        // we assemble it here and then move it into a separate SparseMatrix at the end (same as in step-22).
-
+        // The assembly is similar to step-22. An additional term with gamma
+        // as a coefficient is the Augmented Lagrangian (AL), which is
+        // assembled via grad-div stabilization.  As we discussed in the
+        // introduction, the bottom right block of the system matrix should be
+        // zero. Since the pressure mass matrix is used while creating the
+        // preconditioner, we assemble it here and then move it into a
+        // separate SparseMatrix at the end (same as in step-22).
         for (unsigned int q=0; q<n_q_points; ++q)
           {
             for (unsigned int k=0; k<dofs_per_cell; ++k)
@@ -478,32 +477,29 @@ namespace Step57
                       }
                   }
 
-                if (assemble_rhs)
-                  {
-                    double present_velocity_divergence =  trace(present_velocity_gradients[q]);
-                    local_rhs(i) += ( - viscosity*scalar_product(present_velocity_gradients[q],grad_phi_u[i])
-                                      - present_velocity_gradients[q]*present_velocity_values[q]*phi_u[i]
-                                      + present_pressure_values[q]*div_phi_u[i]
-                                      + present_velocity_divergence*phi_p[i]
-                                      - gamma*present_velocity_divergence*div_phi_u[i])
-                                    * fe_values.JxW(q);
-                  }
+                double present_velocity_divergence =  trace(present_velocity_gradients[q]);
+                local_rhs(i) += ( - viscosity*scalar_product(present_velocity_gradients[q],grad_phi_u[i])
+                                  - present_velocity_gradients[q]*present_velocity_values[q]*phi_u[i]
+                                  + present_pressure_values[q]*div_phi_u[i]
+                                  + present_velocity_divergence*phi_p[i]
+                                  - gamma*present_velocity_divergence*div_phi_u[i])
+                                * fe_values.JxW(q);
               }
           }
 
         cell->get_dof_indices (local_dof_indices);
 
         const ConstraintMatrix &constraints_used = initial_step ? nonzero_constraints : zero_constraints;
-        // Finally we move pressure mass matrix into a separate matrix:
 
         if (assemble_matrix)
           {
             constraints_used.distribute_local_to_global(local_matrix,
+                                                        local_rhs,
                                                         local_dof_indices,
-                                                        system_matrix);
+                                                        system_matrix,
+                                                        system_rhs);
           }
-
-        if (assemble_rhs)
+        else
           {
             constraints_used.distribute_local_to_global(local_rhs,
                                                         local_dof_indices,
@@ -514,6 +510,8 @@ namespace Step57
 
     if (assemble_matrix)
       {
+        // Finally we move pressure mass matrix into a separate matrix:
+
         pressure_mass_matrix.reinit(sparsity_pattern.block(1,1));
         pressure_mass_matrix.copy_from(system_matrix.block(1,1));
 
@@ -528,16 +526,17 @@ namespace Step57
   }
 
   template <int dim>
-  void StationaryNavierStokes<dim>::assemble_matrix(const bool initial_step)
+  void StationaryNavierStokes<dim>::assemble_system(const bool initial_step)
   {
-    assemble_system(initial_step, true, false);
+    assemble(initial_step, true);
   }
 
   template <int dim>
   void StationaryNavierStokes<dim>::assemble_rhs(const bool initial_step)
   {
-    assemble_system(initial_step, false, true);
+    assemble(initial_step, false);
   }
+
   // @sect4{StationaryNavierStokes::solve}
   // In this function, we use FGMRES together with the block preconditioner,
   // which is defined at the beginning of the program, to solve the linear
@@ -660,8 +659,7 @@ namespace Step57
                 setup_dofs();
                 initialize_system();
                 evaluation_point = present_solution;
-                assemble_matrix(first_step);
-                assemble_rhs(first_step);
+                assemble_system(first_step);
                 solve(first_step);
                 present_solution = newton_update;
                 nonzero_constraints.distribute(present_solution);
@@ -678,9 +676,7 @@ namespace Step57
             else
               {
                 evaluation_point = present_solution;
-                assemble_matrix(first_step);
-                if (outer_iteration == 0)
-                  assemble_rhs(first_step);
+                assemble_system(first_step);
                 solve(first_step);
 
                 // To make sure our solution is getting close to the exact solution, we