From: kronbichler Date: Mon, 20 Apr 2009 12:41:04 +0000 (+0000) Subject: Slight update in comments of step-22. Mention the changes on the news page. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6cb6449a3075b4c930f00321b746ab39acf07132;p=dealii-svn.git Slight update in comments of step-22. Mention the changes on the news page. git-svn-id: https://svn.dealii.org/trunk@18655 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/doxygen/tutorial/toc.html b/deal.II/doc/doxygen/tutorial/toc.html index e7c70e15b9..15d838d7df 100644 --- a/deal.II/doc/doxygen/tutorial/toc.html +++ b/deal.II/doc/doxygen/tutorial/toc.html @@ -210,13 +210,13 @@ searching for programs with a certain functionality. Step-22 Solving the Stokes equations of slow fluid flow on adaptive - meshes. More on Schur complement solvers + meshes. More on Schur complement solvers. Advanced use of the + ConstraintMatrix class. Step-23 - Finally a "real" time dependent problem, - the wave equation. + Finally a "real" time dependent problem, the wave equation. diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 0a6ac0c629..d9ac8f816e 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -158,6 +158,17 @@ inconvenience this causes.

General

    +
  1. +

    + Updated: The @ref step_22 "step-22" tutorial program now uses some + advanced features of the ConstraintMatrix such as imhomogeneities + for implementing Dirichlet boundary condition and condensation of + constraints on the fly. +
    + (Martin Kronbichler 2009/04/20) +

    +
    2. +

  2. New: There is now a program, contributed by Jean-Paul Pelteret, that can @@ -224,7 +235,7 @@ inconvenience this causes. (WB 2009/01/12)

    3. - +

  3. Updated: The @ref step_33 "step-33" tutorial program now uses the diff --git a/deal.II/examples/step-22/doc/intro.dox b/deal.II/examples/step-22/doc/intro.dox index e14b0f4990..2c3ea9d302 100644 --- a/deal.II/examples/step-22/doc/intro.dox +++ b/deal.II/examples/step-22/doc/intro.dox @@ -602,7 +602,9 @@ surface that sucks material all the way to the top surface to fill the gap left by the outward motion of material at this location. -

    Using imhomogeneous constraints for implementing Dirichlet boundary conditions

    +

    Implementation

    + +

    Using imhomogeneous constraints for implementing Dirichlet boundary conditions

    In all the previous tutorial programs, we used the ConstraintMatrix merely for handling hanging node constraints (with exception of step-11). However, @@ -650,13 +652,13 @@ side, it distributes the hanging node constraints and additionally implements (inhomogeneous) Dirichlet boundary conditions. That's nice, isn't it? -As results from this discussion, we can conclude that the ConstraintMatrix -provides an alternative to using MatrixTools::apply_boundary_values for -implementing Dirichlet boundary conditions. +We can conclude that the ConstraintMatrix provides an alternative to using +MatrixTools::apply_boundary_values for implementing Dirichlet boundary +conditions. -Other ConstraintMatrix-related optimizations +

    Using ConstraintMatrix for performance optimizations

     Usually, the sparse matrix contains a substantial amount of elements that @@ -670,54 +672,56 @@ considered in this tutorial program. Remember that matrix-vector products or preconditioners operate on all these elements that are zero, which is an inefficiency we have chosen to eliminate in the tutorial program. -Another advantage of directly resolving constrained degrees of freedom is -that we can avoid having all these entries that actually are zero in our -sparse matrix — we do not need them during matrix construction (as -opposed to the traditional algorithms, which first fill the matrix, and only -resolve constraints afterwards). The way we are going to do that is to pass -the information about constraints to the function that generates the -sparsity pattern, and then set a false argument specifying that we -do not intend to use constrained entries: +An advantage of directly resolving constrained degrees of freedom is that we +can avoid having all these entries that actually are zero in our sparse +matrix — we do not need them during matrix construction (as opposed to +the traditional algorithms, which first fill the matrix, and only resolve +constraints afterwards). This will save both memory and time when forming +matrix-vector products. The way we are going to do that is to pass the +information about constraints to the function that generates the sparsity +pattern, and then set a false argument specifying that we do not +intend to use constrained entries: @code DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern, constraints, false); @endcode -This functions saves, by the way, also the call to the condense() -function on the sparsity pattern. +This functions obviates, by the way, also the call to the +condense() function on the sparsity pattern. -

    Implementation

    +

    Performance optimizations

    -The program developed below has seen a lot of TLC. We have run it over -and over under profiling tools (mainly valgrind's cachegrind and -callgrind tools, as well as the KDE valgrind's cachegrind and callgrind +tools, as well as the KDE KCachegrind program for -visualization) to see where the bottlenecks are. This has paid off: -through this effort, the program has become almost twice as fast when +visualization) to see where the bottlenecks are. This has paid off: through +this effort, the program has become about four times as fast when considering the runtime of the refinement cycles zero through three, reducing the overall number of CPU instructions executed from -869,574,060,348 to 199,853,005,625. For higher refinement levels, the -gain is probably even larger since some algorithms that are not ${\cal -O}(N)$ have been eliminated. - -Essentially, there are currently two algorithms in the program that do -not scale linearly with the number of degrees of freedom: renumbering -of degrees of freedom, and the linear solver. As for the first, while -reordering degrees of freedom may not scale linearly, it is an -indispensible part of the overall algorithm as it greatly improves the -quality of the sparse ILU, easily making up for the time spent on -computing the renumbering; graphs and timings to demonstrate this are -shown in the documentation of the DoFRenumbering namespace, also -underlining the choice of the King reordering algorithm chosen below. - -As for the linear solver: as mentioned above, our implementation here -uses a Schur complement formulation. This is not necessarily the very -best choice but demonstrates various important techniques available in -deal.II. The question of which solver is best is again discussed in -the section on improved solvers in the -results part of this program, along with code showing alternative -solvers and a comparison of their results. +869,574,060,348 to 199,853,005,625. For higher refinement levels, the gain +is probably even larger since some algorithms that are not ${\cal O}(N)$ +have been eliminated. + +Essentially, there are currently two algorithms in the program that do not +scale linearly with the number of degrees of freedom: renumbering of degrees +of freedom (which is ${\val O}(N \log N)$, and the linear solver (which is +{\val O}(N^{4/3})$). As for the first, while reordering degrees of freedom +may not scale linearly, it is an indispensible part of the overall algorithm +as it greatly improves the quality of the sparse ILU, easily making up for +the time spent on computing the renumbering; graphs and timings to +demonstrate this are shown in the documentation of the DoFRenumbering +namespace, also underlining the choice of the King reordering algorithm +chosen below. + +As for the linear solver: as mentioned above, our implementation here uses a +Schur complement formulation. This is not necessarily the very best choice +but demonstrates various important techniques available in deal.II. The +question of which solver is best is again discussed in the section on improved solvers in the results part +of this program, along with code showing alternative solvers and a +comparison of their results. Apart from this, many other algorithms have been tested and improved during the creation of this program. For example, in building the sparsity pattern, @@ -727,14 +731,15 @@ large numbers of nonzero entries per row created by our discretization in 3d, leading to a quadratic behavior. Replacing the internal algorithms in deal.II to set many elements at a time, and using a BlockCompressedSimpleSparsityPattern as a better adapted data structure, -removed this bottleneck at the price of some more memory +removed this bottleneck at the price of a slightly higher memory consumption. Likewise, the implementation of the decomposition step in the SparseILU class was very inefficient and has been replaced by one that is -about 10 times faster. Small improvements were applied here and -there. Moreover, the ConstraintMatrix object has been used to eliminate a -lot of entries in the sparse matrix that are eventually going to be zero, -see the section on new features of -constraint matrix. +about 10 times faster. Even the vmult function of the SparseILU has been +improved to save about twenty precent of time. Small improvements were +applied here and there. Moreover, the ConstraintMatrix object has been used +to eliminate a lot of entries in the sparse matrix that are eventually going +to be zero, see the section on using advanced +features of the ConstraintMatrix class. A profile of how many CPU instructions are spent at the various different places in the program during refinement cycles @@ -744,16 +749,16 @@ zero through three in 3d is shown here: As can be seen, at this refinement level approximately three quarters of the instruction count is spent on the actual solver (the SparseILU::vmult calls -on the left, the SparseMatrix::vmult call in the middle, and another -box including multiplications with SparseILU and SparseMatrix stemming from -the solve for U). About one fifth of the instruction count is spent -on matrix assembly and sparse ILU computation (box in the lower right -corner) and the rest on other things. Since floating point operations such -as in the SparseILU::vmult calls typically take much longer than many of the -logical operations and table lookups in matrix assembly, the fraction of the -run time taken up by matrix assembly is actually significantly less than the -fraction of instructions, as will become apparent in the comparison we make -in the results section. +on the left, the SparseMatrix::vmult call in the middle for the Schur +complement solve, and another box representing the multiplications with +SparseILU and SparseMatrix in the solve for U). About one fifth of +the instruction count is spent on matrix assembly and sparse ILU computation +(box in the lower right corner) and the rest on other things. Since floating +point operations such as in the SparseILU::vmult calls typically take much +longer than many of the logical operations and table lookups in matrix +assembly, the fraction of the run time taken up by matrix assembly is +actually significantly less than the fraction of instructions, as will +become apparent in the comparison we make in the results section. For higher refinement levels, the boxes representing the solver as well as the blue box at the top right stemming from reordering algorithm are going diff --git a/deal.II/examples/step-22/doc/results.dox b/deal.II/examples/step-22/doc/results.dox index 3599cda871..00a5e21fc6 100644 --- a/deal.II/examples/step-22/doc/results.dox +++ b/deal.II/examples/step-22/doc/results.dox @@ -175,7 +175,7 @@ Refinement cycle 5 Again, we see that the number of outer iterations does not increase as we refine the mesh. Nevertheless, the compute time increases significantly: for each of the iterations above separately, it takes a -few seconds, a few seconds, 1min, 5min, 21min, and 1h35. This overall +few seconds, a few seconds, 30sec, 4min, 15min, and 1h18min. This overall superlinear (in the number of unknowns) increase in runtime is due to the fact that our inner solver is not ${\cal O}(N)$: a simple experiment shows that as we keep refining the mesh, the average number of @@ -218,15 +218,6 @@ look as follow: @image html step-22.3d.mesh-5.png - - - - @image html step-22.3d.mesh-6.png - - - - - Again, they show essentially the location of singularities introduced @@ -603,158 +594,147 @@ Let's first see the results in 2D: @code Refinement cycle 0 Number of active cells: 64 - Number of degrees of freedom: 679 (594+85) [0.013998 s] - Assembling... [0.004999 s] - Computing preconditioner... [0.003999 s] + Number of degrees of freedom: 679 (594+85) [0.005999 s] + Assembling... [0.002 s] + Computing preconditioner... [0.003 s] Solving... - Schur complement: 11 outer CG iterations for p [0.010999 s] - Block Schur preconditioner: 11 GMRES iterations [0.009998 s] - difference l_infty between solution vectors: 3.18714e-06 + Schur complement: 11 outer CG iterations for p [0.007999 s] + Block Schur preconditioner: 12 GMRES iterations [0.008998 s] + difference l_infty between solution vectors: 8.18909e-07 Refinement cycle 1 Number of active cells: 160 - Number of degrees of freedom: 1683 (1482+201) [0.040994 s] - Assembling... [0.013998 s] - Computing preconditioner... [0.016997 s] + Number of degrees of freedom: 1683 (1482+201) [0.013998 s] + Assembling... [0.005999 s] + Computing preconditioner... [0.012998 s] Solving... - Schur complement: 11 outer CG iterations for p [0.033995 s] - Block Schur preconditioner: 12 GMRES iterations [0.035995 s] - difference l_infty between solution vectors: 9.32671e-06 + Schur complement: 11 outer CG iterations for p [0.029995 s] + Block Schur preconditioner: 12 GMRES iterations [0.030995 s] + difference l_infty between solution vectors: 9.32504e-06 Refinement cycle 2 Number of active cells: 376 - Number of degrees of freedom: 3813 (3370+443) [0.099985 s] - Assembling... [0.033995 s] - Computing preconditioner... [0.052992 s] + Number of degrees of freedom: 3813 (3370+443) [0.031995 s] + Assembling... [0.014998 s] + Computing preconditioner... [0.044994 s] Solving... - Schur complement: 11 outer CG iterations for p [0.110983 s] - Block Schur preconditioner: 12 GMRES iterations [0.110983 s] - difference l_infty between solution vectors: 4.26989e-06 + Schur complement: 11 outer CG iterations for p [0.079987 s] + Block Schur preconditioner: 13 GMRES iterations [0.092986 s] + difference l_infty between solution vectors: 5.40689e-06 Refinement cycle 3 Number of active cells: 880 - Number of degrees of freedom: 8723 (7722+1001) [0.238963 s] - Assembling... [0.076989 s] - Computing preconditioner... [0.141978 s] + Number of degrees of freedom: 8723 (7722+1001) [0.074988 s] + Assembling... [0.035995 s] + Computing preconditioner... [0.110983 s] Solving... - Schur complement: 11 outer CG iterations for p [0.289956 s] - Block Schur preconditioner: 12 GMRES iterations [0.304953 s] - difference l_infty between solution vectors: 1.0266e-05 + Schur complement: 11 outer CG iterations for p [0.19697 s] + Block Schur preconditioner: 13 GMRES iterations [0.242963 s] + difference l_infty between solution vectors: 1.14676e-05 Refinement cycle 4 Number of active cells: 2008 - Number of degrees of freedom: 19383 (17186+2197) [0.561914 s] - Assembling... [0.170974 s] - Computing preconditioner... [0.409938 s] + Number of degrees of freedom: 19383 (17186+2197) [0.180973 s] + Assembling... [0.081987 s] + Computing preconditioner... [0.315952 s] Solving... - Schur complement: 11 outer CG iterations for p [0.740887 s] - Block Schur preconditioner: 13 GMRES iterations [0.837873 s] - difference l_infty between solution vectors: 3.13139e-05 + Schur complement: 11 outer CG iterations for p [0.673898 s] + Block Schur preconditioner: 13 GMRES iterations [0.778882 s] + difference l_infty between solution vectors: 3.13142e-05 Refinement cycle 5 Number of active cells: 4288 - Number of degrees of freedom: 40855 (36250+4605) [1.25581 s] - Assembling... [0.353946 s] - Computing preconditioner... [1.21681 s] + Number of degrees of freedom: 40855 (36250+4605) [0.386941 s] + Assembling... [0.171974 s] + Computing preconditioner... [0.766883 s] Solving... - Schur complement: 11 outer CG iterations for p [1.60576 s] - Block Schur preconditioner: 13 GMRES iterations [1.84572 s] - difference l_infty between solution vectors: 8.59663e-05 - -Refinement cycle 6 - Number of active cells: 8896 - Number of degrees of freedom: 83885 (74474+9411) [2.54961 s] - Assembling... [0.732889 s] - Computing preconditioner... [3.9704 s] - Solving... - Schur complement: 11 outer CG iterations for p [3.70344 s] - Block Schur preconditioner: 13 GMRES iterations [4.18436 s] - difference l_infty between solution vectors: 0.00022514 + Schur complement: 11 outer CG iterations for p [1.65275 s] + Block Schur preconditioner: 13 GMRES iterations [1.81372 s] + difference l_infty between solution vectors: 8.59668e-05 @endcode -We see that there is no huge difference in the solution time between -the block Schur complement preconditioner solver and the Schur -complement itself. The -reason is simple: we used a direct solve as preconditioner for $A$ - so -there is no substantial gain by avoiding the inner iterations. We see -that the number of iterations has slightly increased for GMRES, but all in all -the two choices are fairly similar. +We see that there is no huge difference in the solution time between the +block Schur complement preconditioner solver and the Schur complement +itself. The reason is simple: we used a direct solve as preconditioner for +$A$ - so we cannot expect any gain by avoiding the inner iterations. We see +that the number of iterations has slightly increased for GMRES, but all in +all the two choices are fairly similar. The picture of course changes in 3D: @code Refinement cycle 0 Number of active cells: 32 - Number of degrees of freedom: 1356 (1275+81) [0.336021 s] - Assembling... [0.204013 s] - Computing preconditioner... [0.128008 s] + Number of degrees of freedom: 1356 (1275+81) [0.025996 s] + Assembling... [0.056992 s] + Computing preconditioner... [0.027995 s] Solving... - Schur complement: 13 outer CG iterations for p [0.552034 s] - Block Schur preconditioner: 23 GMRES iterations [0.072005 s] - difference l_infty between solution vectors: 1.11101e-05 + Schur complement: 13 outer CG iterations for p [0.275958 s] + Block Schur preconditioner: 23 GMRES iterations [0.042994 s] + difference l_infty between solution vectors: 1.11307e-05 Refinement cycle 1 Number of active cells: 144 - Number of degrees of freedom: 5088 (4827+261) [0.79605 s] - Assembling... [0.428026 s] - Computing preconditioner... [0.30802 s] + Number of degrees of freedom: 5088 (4827+261) [0.102984 s] + Assembling... [0.254961 s] + Computing preconditioner... [0.161976 s] Solving... - Schur complement: 14 outer CG iterations for p [4.01625 s] - Block Schur preconditioner: 41 GMRES iterations [0.632039 s] - difference l_infty between solution vectors: 2.36202e-05 + Schur complement: 14 outer CG iterations for p [2.43963 s] + Block Schur preconditioner: 42 GMRES iterations [0.352946 s] + difference l_infty between solution vectors: 9.07409e-06 Refinement cycle 2 Number of active cells: 704 - Number of degrees of freedom: 22406 (21351+1055) [4.12826 s] - Assembling... [2.15614 s] - Computing preconditioner... [1.6721 s] + Number of degrees of freedom: 22406 (21351+1055) [0.52592 s] + Assembling... [1.24481 s] + Computing preconditioner... [0.948856 s] Solving... - Schur complement: 14 outer CG iterations for p [30.6179 s] - Block Schur preconditioner: 75 GMRES iterations [6.01638 s] - difference l_infty between solution vectors: 3.86168e-05 + Schur complement: 14 outer CG iterations for p [22.2056 s] + Block Schur preconditioner: 78 GMRES iterations [4.75928 s] + difference l_infty between solution vectors: 2.48042e-05 Refinement cycle 3 Number of active cells: 3168 - Number of degrees of freedom: 93176 (89043+4133) [18.8452 s] - Assembling... [9.47659 s] - Computing preconditioner... [7.47647 s] + Number of degrees of freedom: 93176 (89043+4133) [2.66759 s] + Assembling... [5.66014 s] + Computing preconditioner... [4.69529 s] Solving... - Schur complement: 15 outer CG iterations for p [228.366 s] - Block Schur preconditioner: 159 GMRES iterations [62.6799 s] - difference l_infty between solution vectors: 7.34297e-05 + Schur complement: 15 outer CG iterations for p [235.74 s] + Block Schur preconditioner: 162 GMRES iterations [63.7883 s] + difference l_infty between solution vectors: 5.62978e-05 Refinement cycle 4 Number of active cells: 11456 - Number of degrees of freedom: 327808 (313659+14149) [77.5688 s] - Assembling... [34.9422 s] - Computing preconditioner... [27.9577 s] + Number of degrees of freedom: 327808 (313659+14149) [12.0242 s] + Assembling... [20.2669 s] + Computing preconditioner... [17.3384 s] Solving... - Schur complement: 15 outer CG iterations for p [1040.28 s] - Block Schur preconditioner: 292 GMRES iterations [414.738 s] - difference l_infty between solution vectors: 0.00019614 + Schur complement: 15 outer CG iterations for p [817.287 s] + Block Schur preconditioner: 294 GMRES iterations [347.307 s] + difference l_infty between solution vectors: 0.000107536 Refinement cycle 5 Number of active cells: 45056 - Number of degrees of freedom: 1254464 (1201371+53093) [520.285 s] - Assembling... [136.909 s] - Computing preconditioner... [115.455 s] + Number of degrees of freedom: 1254464 (1201371+53093) [89.8533 s] + Assembling... [80.3588 s] + Computing preconditioner... [73.0849 s] Solving... - Schur complement: 14 outer CG iterations for p [6556.39 s] - Block Schur preconditioner: 590 GMRES iterations [3906.06 s] - difference l_infty between solution vectors: 0.000428885 + Schur complement: 14 outer CG iterations for p [4401.66 s] + Block Schur preconditioner: 587 GMRES iterations [3083.21 s] + difference l_infty between solution vectors: 0.00025531 @endcode Here, the block preconditioned solver is clearly superior to the Schur -complement, but the advantage gets less for more mesh points. This was expected -from the discussion above. It is also apparent that GMRES scales worse with the -problem size than CG (as explained above). -Nonetheless, the improvement by a factor of 3-5 for moderate problem sizes -is quite impressive. - -
    Combining block preconditioner and multigrid
    -An ultimate linear solver for this problem could be imagined as a combination of -an optimal preconditioner for $A$ (e.g. multigrid) and the block preconditioner -described above. +complement, but the advantage gets less for more mesh points. This is +because GMRES(k) scales worse with the problem size than CG, as we discussed +above. Nonetheless, the improvement by a factor of 3-5 for moderate problem +sizes is quite impressive. + +
    Combining block preconditioner and multigrid
    An ultimate linear +solver for this problem could be imagined as a combination of an optimal +preconditioner for $A$ (e.g. multigrid) and the block preconditioner +described above, which is the approach taken in the @ref step_31 "step-31" +tutorial program.
    No block matrices and vectors
    Another possibility that can be taken into account is to not set up a block diff --git a/deal.II/examples/step-22/doc/step-22.3d.mesh-6.png b/deal.II/examples/step-22/doc/step-22.3d.mesh-6.png deleted file mode 100644 index c1d5c55f19..0000000000 Binary files a/deal.II/examples/step-22/doc/step-22.3d.mesh-6.png and /dev/null differ diff --git a/deal.II/examples/step-22/step-22.cc b/deal.II/examples/step-22/step-22.cc index 5f7f81b35f..d1092b3ff5 100644 --- a/deal.II/examples/step-22/step-22.cc +++ b/deal.II/examples/step-22/step-22.cc @@ -934,11 +934,11 @@ void StokesProblem::assemble_system () // the indices in case we are // pointing into the empty part of // the local matrix. - cell->get_dof_indices (local_dof_indices); for (unsigned int i=0; iget_dof_indices (local_dof_indices); constraints.distribute_local_to_global (local_matrix, local_rhs, local_dof_indices, system_matrix, system_rhs);