]> https://gitweb.dealii.org/ - dealii.git/commitdiff
step-9: Modify the refinement cycle and linear solver. 6900/head
authorDavid Wells <drwells@email.unc.edu>
Sun, 15 Jul 2018 16:37:23 +0000 (12:37 -0400)
committerDavid Wells <drwells@email.unc.edu>
Fri, 3 Aug 2018 02:39:05 +0000 (22:39 -0400)
1. The different refinement strategy seems to work better with the higher degree
   finite element.
2. Bicgstab has trouble converging (and does not converge in the max norm) for
   finer grids. GMRES is more consistent for this problem.
3. I added some extra console output to show the residual (which does not always
   converge with Bicgstab) and the number of DoFs.

examples/step-9/doc/intro.dox
examples/step-9/doc/results.dox
examples/step-9/step-9.cc

index 1bd0412c5c75065afe616d9cbf9d93cb0bb5c3d5..e8c4c6ca52f04e69d0926a3489dddc13c2d4e2e7 100644 (file)
@@ -180,8 +180,8 @@ as system matrix. We will assemble this matrix in the program.
 
 As the resulting matrix is no longer symmetric positive definite, we cannot
 use the usual Conjugate Gradient method (implemented in the
-SolverCG class) to solve the system. Instead, we use the BiCGStab (bi-conjugate
-gradient stabilized) method (implemented in SolverBicgstab) that is suitable
+SolverCG class) to solve the system. Instead, we use the GMRES (Generalized
+Minimum RESidual) method (implemented in SolverGMRES) that is suitable
 for problems of the kind we have here.
 
 
index a13a493a6f650023fb8b7d9ed240001df34a397a..487095b384922544b7504bbfc6b8d9b16532e221 100644 (file)
@@ -6,23 +6,55 @@ consist of the console output, some grid files, and the solution on
 each of these grids. First for the console output:
 @code
 Cycle 0:
-   Number of active cells:       256
-   Number of degrees of freedom: 289
+   Number of active cells:              64
+   Number of degrees of freedom:        1681
+   Iterations required for convergence: 298
+   Max norm of residual:                3.60316e-12
 Cycle 1:
-   Number of active cells:       643
-   Number of degrees of freedom: 793
+   Number of active cells:              124
+   Number of degrees of freedom:        3537
+   Iterations required for convergence: 415
+   Max norm of residual:                3.70682e-12
 Cycle 2:
-   Number of active cells:       1669
-   Number of degrees of freedom: 1950
+   Number of active cells:              247
+   Number of degrees of freedom:        6734
+   Iterations required for convergence: 543
+   Max norm of residual:                7.19716e-13
 Cycle 3:
-   Number of active cells:       4231
-   Number of degrees of freedom: 4923
+   Number of active cells:              502
+   Number of degrees of freedom:        14105
+   Iterations required for convergence: 666
+   Max norm of residual:                3.45628e-13
 Cycle 4:
-   Number of active cells:       10753
-   Number of degrees of freedom: 12175
+   Number of active cells:              1003
+   Number of degrees of freedom:        27462
+   Iterations required for convergence: 1064
+   Max norm of residual:                1.86495e-13
 Cycle 5:
-   Number of active cells:       27004
-   Number of degrees of freedom: 29810
+   Number of active cells:              1993
+   Number of degrees of freedom:        55044
+   Iterations required for convergence: 1251
+   Max norm of residual:                1.28765e-13
+Cycle 6:
+   Number of active cells:              3985
+   Number of degrees of freedom:        108492
+   Iterations required for convergence: 2035
+   Max norm of residual:                6.78085e-14
+Cycle 7:
+   Number of active cells:              7747
+   Number of degrees of freedom:        210612
+   Iterations required for convergence: 2187
+   Max norm of residual:                2.61457e-14
+Cycle 8:
+   Number of active cells:              15067
+   Number of degrees of freedom:        406907
+   Iterations required for convergence: 3079
+   Max norm of residual:                2.9932e-14
+Cycle 9:
+   Number of active cells:              29341
+   Number of degrees of freedom:        780591
+   Iterations required for convergence: 3913
+   Max norm of residual:                8.15689e-15
 @endcode
 
 Quite a number of cells are used on the finest level to resolve the features of
index b5c648a41691a10b959f7392a0b00b10039664b1..7a19ff1ecb646f69e8327a43695fe908f7f4951c 100644 (file)
@@ -27,7 +27,7 @@
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/sparse_matrix.h>
 #include <deal.II/lac/dynamic_sparsity_pattern.h>
-#include <deal.II/lac/solver_bicgstab.h>
+#include <deal.II/lac/solver_gmres.h>
 #include <deal.II/lac/precondition.h>
 #include <deal.II/lac/affine_constraints.h>
 #include <deal.II/grid/tria.h>
@@ -759,27 +759,36 @@ namespace Step9
   // Here comes the linear solver routine. As the system is no longer
   // symmetric positive definite as in all the previous examples, we cannot
   // use the Conjugate Gradient method anymore. Rather, we use a solver that
-  // is tailored to nonsymmetric systems like the one at hand, the BiCGStab
-  // method. As preconditioner, we use the Jacobi method.
+  // is more general and does not rely on any special properties of the
+  // matrix: the GMRES method. GMRES, like the conjugate gradient method,
+  // requires a decent preconditioner: we use a Jacobi preconditioner here,
+  // which works well enough for this problem.
   template <int dim>
   void AdvectionProblem<dim>::solve()
   {
-    SolverControl    solver_control(1000, 1e-12);
-    SolverBicgstab<> bicgstab(solver_control);
-
+    SolverControl        solver_control(std::max<std::size_t>(1000,
+                                                       system_rhs.size() / 10),
+                                 1e-10 * system_rhs.l2_norm());
+    SolverGMRES<>        solver(solver_control);
     PreconditionJacobi<> preconditioner;
     preconditioner.initialize(system_matrix, 1.0);
+    solver.solve(system_matrix, solution, system_rhs, preconditioner);
+
+    Vector<double> residual(dof_handler.n_dofs());
 
-    bicgstab.solve(system_matrix, solution, system_rhs, preconditioner);
+    system_matrix.vmult(residual, solution);
+    residual -= system_rhs;
+    std::cout << "   Iterations required for convergence: "
+              << solver_control.last_step() << '\n'
+              << "   Max norm of residual:                "
+              << residual.linfty_norm() << '\n';
 
     hanging_node_constraints.distribute(solution);
   }
 
   // The following function refines the grid according to the quantity
   // described in the introduction. The respective computations are made in
-  // the class <code>GradientEstimation</code>. The only difference to
-  // previous examples is that we refine a little more aggressively (0.5
-  // instead of 0.3 of the number of cells).
+  // the class <code>GradientEstimation</code>.
   template <int dim>
   void AdvectionProblem<dim>::refine_grid()
   {
@@ -791,7 +800,7 @@ namespace Step9
 
     GridRefinement::refine_and_coarsen_fixed_number(triangulation,
                                                     estimated_error_per_cell,
-                                                    0.5,
+                                                    0.3,
                                                     0.03);
 
     triangulation.execute_coarsening_and_refinement();
@@ -850,14 +859,14 @@ namespace Step9
   template <int dim>
   void AdvectionProblem<dim>::run()
   {
-    for (unsigned int cycle = 0; cycle < 6; ++cycle)
+    for (unsigned int cycle = 0; cycle < 10; ++cycle)
       {
         std::cout << "Cycle " << cycle << ':' << std::endl;
 
         if (cycle == 0)
           {
             GridGenerator::hyper_cube(triangulation, -1, 1);
-            triangulation.refine_global(4);
+            triangulation.refine_global(3);
           }
         else
           {
@@ -865,13 +874,13 @@ namespace Step9
           }
 
 
-        std::cout << "   Number of active cells:       "
+        std::cout << "   Number of active cells:              "
                   << triangulation.n_active_cells() << std::endl;
 
         setup_system();
 
-        std::cout << "   Number of degrees of freedom: " << dof_handler.n_dofs()
-                  << std::endl;
+        std::cout << "   Number of degrees of freedom:        "
+                  << dof_handler.n_dofs() << std::endl;
 
         assemble_system();
         solve();

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.