left by the outward motion of material at this location.
-<h3>Using imhomogeneous constraints for implementing Dirichlet boundary conditions</h3>
+<h3>Implementation</h3>
+
+<h4>Using imhomogeneous constraints for implementing Dirichlet boundary conditions</h4>
In all the previous tutorial programs, we used the ConstraintMatrix merely
for handling hanging node constraints (with exception of step-11). However,
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.
<a name="constraint-matrix">
-Other ConstraintMatrix-related optimizations
+<h4>Using ConstraintMatrix for performance optimizations</h4>
</a>
Usually, the sparse matrix contains a substantial amount of elements that
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 <tt>false</tt> 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 <tt>false</tt> 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 <tt>condense()</tt>
-function on the sparsity pattern.
+This functions obviates, by the way, also the call to the
+<tt>condense()</tt> function on the sparsity pattern.
-<h3>Implementation</h3>
+<h4>Performance optimizations</h4>
-The program developed below has seen a lot of TLC. We have run it over
-
and over under profiling tools (mainly <a
-href="http://www.valgrind.org/">valgrind</a>'s cachegrind and
-
callgrind tools, as well as the KDE <a
+The program developed below has seen a lot of TLC. We have run it over and
+over under profiling tools (mainly <a
+href="http://www.valgrind.org/">valgrind</a>'s cachegrind and callgrind
+tools, as well as the KDE <a
href="http://kcachegrind.sourceforge.net/">KCachegrind</a> 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 <a href="#improved-solver">section on improved solvers in the
-results part</a> 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 <a
+href="#improved-solver">section on improved solvers in the results part</a>
+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,
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 <a href="#constraint-matrix">the section on new features of
-constraint matrix</a>.
+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 <a href="#constraint-matrix">the section on using advanced
+features of the ConstraintMatrix class</a>.
A profile of how many CPU instructions are spent at the various
different places in the program during refinement cycles
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<double>::vmult call in the middle, and another
-box including multiplications with SparseILU and SparseMatrix stemming from
-the solve for <i>U</i>). 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 <i>U</i>). 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
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
@image html step-22.3d.mesh-5.png
</td>
</tr>
-
- <tr>
- <td ALIGN="center">
- @image html step-22.3d.mesh-6.png
- </td>
-
- <td ALIGN="center">
- </td>
- </tr>
</table>
Again, they show essentially the location of singularities introduced
@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.
-
-<h5>Combining block preconditioner and multigrid</h5>
-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.
+
+<h5>Combining block preconditioner and multigrid</h5> 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.
<h5>No block matrices and vectors</h5>
Another possibility that can be taken into account is to not set up a block
projects / dealii-svn.git / commitdiff
