From 8ea108901e364c16cdefad1c6a4c627d0d0ca885 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Wed, 1 Feb 2017 09:59:44 -0500 Subject: [PATCH] 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) --- examples/step-57/doc/results.dox | 391 +++++++++++++++---------------- examples/step-57/step-57.cc | 80 +++---- 2 files changed, 226 insertions(+), 245 deletions(-) diff --git a/examples/step-57/doc/results.dox b/examples/step-57/doc/results.dox index 86d9c8c7a4..42ac92872e 100644 --- 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. - - - +
 
+ + - - + + @@ -31,113 +31,92 @@ every mesh is shown. The data in the table shows that Newton's iteration converg - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Re=400 Mesh0 Mesh1 Mesh2 Mesh3 Mesh4
Newton iter Residual FGMRES FGMRES Residual FGMRES
1 7.40396e-33 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
3 1.60421e-34 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
5 1.34601e-54                
6 2.5235e-8 5                
7 1.38899e-12 4                
8 4.68224e-15 4                
13.7112e-0356.4189e-0332.4338e-0331.0570e-0334.9499e-043
27.0849e-045.0000e+009.9458e-0451.1409e-0461.3544e-0561.4171e-066
31.9980e-055.0000e+004.5007e-0552.9020e-0854.4021e-1066.3435e-116
42.3165e-096.0000e+001.6891e-0751.2338e-1471.8506e-1488.8563e-158
51.2585e-137.0000e+001.4520e-1161.9044e-138
61.3998e-158
+ + + + + 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. - +
- + @@ -218,98 +197,104 @@ mesh. The quadratic convergence is clearly visible in the table. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
 Re=7500 Mesh0 Mesh1 Mesh2Residual FGMRES
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
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
5     1.29739e-77 1.25167e-13 9 1.76128e-14 10    
6     4.43518e-117            
7     6.42323e-15 9           
11.8922e-0664.2506e-0331.4299e-0334.8793e-0421.8998e-042
23.1644e-0981.3732e-0374.1506e-0479.1119e-0581.3555e-058
31.7611e-1492.1946e-0461.7881e-0565.2678e-0779.3739e-097
48.8269e-0666.8210e-0972.2770e-1181.2588e-139
51.2974e-0771.2515e-1391.7801e-141
64.4352e-117
76.2863e-159
+ + + + + + The sequence of generated grids looks like this: diff --git a/examples/step-57/step-57.cc b/examples/step-57/step-57.cc index 0760f493cf..a54308a99f 100644 --- 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 - void StationaryNavierStokes::assemble_system(const bool initial_step, - const bool assemble_matrix, - const bool assemble_rhs) + void StationaryNavierStokes::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 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; qget_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 - void StationaryNavierStokes::assemble_matrix(const bool initial_step) + void StationaryNavierStokes::assemble_system(const bool initial_step) { - assemble_system(initial_step, true, false); + assemble(initial_step, true); } template void StationaryNavierStokes::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 -- 2.39.5