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
#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>
// 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()
{
GridRefinement::refine_and_coarsen_fixed_number(triangulation,
estimated_error_per_cell,
- 0.5,
+ 0.3,
0.03);
triangulation.execute_coarsening_and_refinement();
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
{
}
- 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();