From: David Wells Date: Sun, 15 Jul 2018 16:37:23 +0000 (-0400) Subject: step-9: Modify the refinement cycle and linear solver. X-Git-Tag: v9.1.0-rc1~853^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=aaca740e971b616b50a89683de1b4b7438f81d98;p=dealii.git step-9: Modify the refinement cycle and linear solver. 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. --- diff --git a/examples/step-9/doc/intro.dox b/examples/step-9/doc/intro.dox index 1bd0412c5c..e8c4c6ca52 100644 --- a/examples/step-9/doc/intro.dox +++ b/examples/step-9/doc/intro.dox @@ -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. diff --git a/examples/step-9/doc/results.dox b/examples/step-9/doc/results.dox index a13a493a6f..487095b384 100644 --- a/examples/step-9/doc/results.dox +++ b/examples/step-9/doc/results.dox @@ -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 diff --git a/examples/step-9/step-9.cc b/examples/step-9/step-9.cc index b5c648a416..7a19ff1ecb 100644 --- a/examples/step-9/step-9.cc +++ b/examples/step-9/step-9.cc @@ -27,7 +27,7 @@ #include #include #include -#include +#include #include #include #include @@ -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 void AdvectionProblem::solve() { - SolverControl solver_control(1000, 1e-12); - SolverBicgstab<> bicgstab(solver_control); - + SolverControl solver_control(std::max(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 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 GradientEstimation. 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 GradientEstimation. template void AdvectionProblem::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 void AdvectionProblem::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();