in the early tutorial programs. The improvements are still rather
unsatisfactory, if one compares with an elliptic problem of similar size.
+<a name="improved-ilu">
+<h4>Better ILU decomposition by smart reordering</h4>
+</a>
+A first attempt to improve the speed of the linear solution process is to choose
+a dof reordering that makes the ILU being closer to a full LU decomposition, as
+already mentioned in the in-code comments. The DoFRenumbering namespace compares
+several choices for the renumbering of dofs for the Stokes equations. The best
+result regarding the computing time was found for the King ordering, which is
+accessed through the call DoFRenumbering::boost::king_ordering. With that
+program, the inner solver needs considerably less operations, e.g. about 62
+inner CG iterations for the inversion of $A$ at cycle 4 compared to about 75
+iterations with the standard Cuthill_McKee-algorithm. Also, the computing time
+at cycle 4 decreased from about 20 to 13 minutes for the <code>solve()</code>
+call. However, the King ordering (and the orderings provided by the
+DoFRenumbering::boost namespace in general) has a serious drawback - it uses
+much more memory than the in-build deal versions, since it acts on abstract
+graphs rather than the geometry provided by the triangulation. In the present
+case, the renumbering takes about 5 times as much memory, which yields an
+infeasible algorithm for the two last cycles in 3D with 1.2 and 4.6 million
+unknowns, respectively.
+
<h4>Better preconditioner for the inner CG solver</h4>
The first way to improve the situation would be to choose a preconditioner that
makes CG for the (0,0) matrix $A$ converge in a mesh-independent number of
@code
Refinement cycle 0
Number of active cells: 64
- Number of degrees of freedom: 679 (594+85) [0.022997 s]
- Assembling... [0.013998 s]
- Computing preconditioner... [0.004999 s]
+ Number of degrees of freedom: 679 (594+85) [0.013998 s]
+ Assembling... [0.004999 s]
+ Computing preconditioner... [0.003999 s]
Solving...
- Schur complement: 11 outer CG iterations for p [0.010998 s]
- Block Schur preconditioner: 11 GMRES iterations [0.009999 s]
+ 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
Refinement cycle 1
Number of active cells: 160
- Number of degrees of freedom: 1683 (1482+201) [0.072989 s]
- Assembling... [0.034994 s]
- Computing preconditioner... [0.017998 s]
+ Number of degrees of freedom: 1683 (1482+201) [0.040994 s]
+ Assembling... [0.013998 s]
+ Computing preconditioner... [0.016997 s]
Solving...
- Schur complement: 11 outer CG iterations for p [0.036994 s]
- Block Schur preconditioner: 12 GMRES iterations [0.036994 s]
- difference l_infty between solution vectors: 1.57149e-06
+ 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
Refinement cycle 2
Number of active cells: 376
- Number of degrees of freedom: 3813 (3370+443) [0.184972 s]
- Assembling... [0.083987 s]
- Computing preconditioner... [0.054992 s]
+ Number of degrees of freedom: 3813 (3370+443) [0.099985 s]
+ Assembling... [0.033995 s]
+ Computing preconditioner... [0.052992 s]
Solving...
- Schur complement: 11 outer CG iterations for p [0.115983 s]
- Block Schur preconditioner: 12 GMRES iterations [0.123981 s]
- difference l_infty between solution vectors: 4.25998e-06
+ 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
Refinement cycle 3
Number of active cells: 880
- Number of degrees of freedom: 8723 (7722+1001) [0.410937 s]
- Assembling... [0.19597 s]
- Computing preconditioner... [0.147978 s]
+ Number of degrees of freedom: 8723 (7722+1001) [0.238963 s]
+ Assembling... [0.076989 s]
+ Computing preconditioner... [0.141978 s]
Solving...
- Schur complement: 11 outer CG iterations for p [0.324951 s]
- Block Schur preconditioner: 12 GMRES iterations [0.32495 s]
- difference l_infty between solution vectors: 1.02911e-05
+ 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
Refinement cycle 4
Number of active cells: 2008
- Number of degrees of freedom: 19383 (17186+2197) [0.92086 s]
- Assembling... [0.439933 s]
- Computing preconditioner... [0.465929 s]
+ Number of degrees of freedom: 19383 (17186+2197) [0.561914 s]
+ Assembling... [0.170974 s]
+ Computing preconditioner... [0.409938 s]
Solving...
- Schur complement: 11 outer CG iterations for p [0.796879 s]
- Block Schur preconditioner: 13 GMRES iterations [0.896864 s]
- difference l_infty between solution vectors: 3.12965e-05
+ 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
Refinement cycle 5
Number of active cells: 4288
- Number of degrees of freedom: 40855 (36250+4605) [2.08068 s]
- Assembling... [0.936857 s]
- Computing preconditioner... [1.53277 s]
+ Number of degrees of freedom: 40855 (36250+4605) [1.25581 s]
+ Assembling... [0.353946 s]
+ Computing preconditioner... [1.21681 s]
Solving...
- Schur complement: 11 outer CG iterations for p [1.80873 s]
- Block Schur preconditioner: 13 GMRES iterations [2.01869 s]
- difference l_infty between solution vectors: 8.60583e-05
+ 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) [4.30634 s]
- Assembling... [1.93171 s]
- Computing preconditioner... [4.76627 s]
+ 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 [4.02939 s]
- Block Schur preconditioner: 13 GMRES iterations [4.5883 s]
- difference l_infty between solution vectors: 0.000225481
+ 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
@endcode
We see that there is no huge difference in the solution time between
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. For the computing times reported below,
-cycle 5 was run without the memory
-consuming options BlockCompressedSetSparsityPattern and
-DoFRenumbering::boost::king_ordering and Cuthill_McKee numbering and
-CompressedBlockSparsityPattern instead - calculations from cycle 0-4 with that
-setup are about 10-30% slower (in the setup_dofs and solver function) than with
-the one standard configuration as above.
+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.189971 s]
- Assembling... [0.46293 s]
- Computing preconditioner... [0.055991 s]
+ Number of degrees of freedom: 1356 (1275+81) [0.13498 s]
+ Assembling... [0.081987 s]
+ Computing preconditioner... [0.056992 s]
Solving...
- Schur complement: 13 outer CG iterations for p [0.307953 s]
- Block Schur preconditioner: 20 GMRES iterations [0.047993 s]
- difference l_infty between solution vectors: 1.03218e-05
+ Schur complement: 13 outer CG iterations for p [0.334949 s]
+ Block Schur preconditioner: 23 GMRES iterations [0.051992 s]
+ difference l_infty between solution vectors: 1.11101e-05
Refinement cycle 1
Number of active cells: 144
- Number of degrees of freedom: 5088 (4827+261) [1.05084 s]
- Assembling... [2.12268 s]
- Computing preconditioner... [0.307953 s]
+ Number of degrees of freedom: 5088 (4827+261) [0.739887 s]
+ Assembling... [0.39494 s]
+ Computing preconditioner... [0.310953 s]
Solving...
- Schur complement: 14 outer CG iterations for p [2.72559 s]
- Block Schur preconditioner: 37 GMRES iterations [0.422936 s]
- difference l_infty between solution vectors: 1.17651e-05
+ Schur complement: 14 outer CG iterations for p [2.74158 s]
+ Block Schur preconditioner: 41 GMRES iterations [0.451931 s]
+ difference l_infty between solution vectors: 2.36202e-05
Refinement cycle 2
Number of active cells: 704
- Number of degrees of freedom: 22406 (21351+1055) [5.60815 s]
- Assembling... [10.4114 s]
- Computing preconditioner... [1.59876 s]
+ Number of degrees of freedom: 22406 (21351+1055) [3.90041 s]
+ Assembling... [1.9877 s]
+ Computing preconditioner... [1.67175 s]
Solving...
- Schur complement: 14 outer CG iterations for p [18.7521 s]
- Block Schur preconditioner: 62 GMRES iterations [3.76743 s]
- difference l_infty between solution vectors: 2.46681e-05
+ Schur complement: 14 outer CG iterations for p [21.7047 s]
+ Block Schur preconditioner: 75 GMRES iterations [4.6303 s]
+ difference l_infty between solution vectors: 3.86168e-05
Refinement cycle 3
Number of active cells: 3168
- Number of degrees of freedom: 93176 (89043+4133) [23.6404 s]
- Assembling... [46.3899 s]
- Computing preconditioner... [7.27989 s]
+ Number of degrees of freedom: 93176 (89043+4133) [17.3564 s]
+ Assembling... [8.67068 s]
+ Computing preconditioner... [7.41487 s]
Solving...
- Schur complement: 15 outer CG iterations for p [205.028 s]
- Block Schur preconditioner: 130 GMRES iterations [47.8267 s]
- difference l_infty between solution vectors: 4.91733e-05
+ Schur complement: 15 outer CG iterations for p [253.15 s]
+ Block Schur preconditioner: 159 GMRES iterations [54.0888 s]
+ difference l_infty between solution vectors: 7.34297e-05
Refinement cycle 4
Number of active cells: 11456
- Number of degrees of freedom: 327808 (313659+14149) [93.4108 s]
- Assembling... [173.294 s]
- Computing preconditioner... [27.0479 s]
+ Number of degrees of freedom: 327808 (313659+14149) [73.0609 s]
+ Assembling... [31.9341 s]
+ Computing preconditioner... [27.6778 s]
Solving...
- Schur complement: 15 outer CG iterations for p [824.465 s]
- Block Schur preconditioner: 209 GMRES iterations [289.923 s]
- difference l_infty between solution vectors: 0.000185872
-
-Refinement cycle 5
- Number of active cells: 45056
- Number of degrees of freedom: 1254464 (1201371+53093) [2141.82 s]
- Assembling... [699.713 s]
- Computing preconditioner... [113.561 s]
- Solving...
- Schur complement: 14 outer CG iterations for p [6461.56 s]
- Block Schur preconditioner: 557 GMRES iterations [3640.34 s]
- difference l_infty between solution vectors: 0.000427821
+ Schur complement: 15 outer CG iterations for p [1190.48 s]
+ Block Schur preconditioner: 292 GMRES iterations [416.802 s]
+ difference l_infty between solution vectors: 0.00019614
@endcode
Here, the block preconditioned solver is clearly superior to the Schur
Nonetheless, the improvement by a factor of 3-5 for moderate problem sizes
is quite impressive.
+<h4>Combining block preconditioner and multigrid</h4>
+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.
+
<h4>No block matrices and vectors</h4>
Another possibility that can be taken into account is to not set up a block
system, but rather solve the system of velocity and pressure all at once. The
// the results we obtain with several
// of these algorithms based on the
// testcase discussed here in this
- // tutorial program. It turns out
- // that not the traditional
- // Cuthill-McKee algorithm already
- // used in some of the previous
- // tutorial programs is the best, but
- // the King ordering. deal.II
- // implements them by interfacing
- // with the Boost Graph Library,
- // which is part of the deal.II
- // distribution. The call to use it
- // below is
- // DoFRenumbering::boost::king_ordering.
- //
+ // tutorial program. Here, we will
+ // use the traditional Cuthill-McKee
+ // algorithm already used in some of
+ // the previous tutorial programs.
+ // In the
+ // <a href="#improved-ilu">section on improved ILU</a>
+ // we're going to discuss this issue
+ // in more detail.
+
// There is one more change compared
// to previous tutorial programs:
// There is no reason in sorting the
system_matrix.clear ();
dof_handler.distribute_dofs (fe);
- DoFRenumbering::boost::king_ordering (dof_handler);
+ DoFRenumbering::Cuthill_McKee (dof_handler);
std::vector<unsigned int> block_component (dim+1,0);
block_component[dim] = 1;
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
- // We can then start the loop over all
- // cells:
+ // As an extension over step-20
+ // and step-21, we include a few
+ // optimizations that make assembly
+ // faster for this particular problem.
+ // The improvements are based on the
+ // observation that we do a few
+ // calculations too many times when
+ // we do as in step-20: The
+ // symmetric gradient actually has
+ // <code>dofs_per_cell</code>
+ // different values per quadrature
+ // point, but we calculate it
+ // <code>dofs_per_cell*dofs_per_cell</code>
+ // times - for both the loop over
+ // <code>i</code> and the loop over
+ // <code>j</code>. So what we're
+ // going to do here is to avoid
+ // such double calculations by
+ // getting a vector of rank-2
+ // tensors (and similarly for
+ // the divergence and the basis
+ // function value on pressure)
+ // at the quadrature point prior
+ // to starting the loop over the
+ // dofs on the cell. First, we
+ // create the respective objects
+ // that will hold the respective
+ // values. Then, we start the
+ // loop over all cells and the loop
+ // over the quadrature points,
+ // where we first extract these
+ // values. There is one more
+ // optimization we implement here:
+ // the local matrix (as well as
+ // the global one) is going to
+ // be symmetric, since the all
+ // the operations involved are
+ // symmetric with respect to $i$
+ // and $j$. This is implemented by
+ // simply running the inner loop
+ // not to <code>dofs_per_cell</code>,
+ // but only up to <code>i</code>,
+ // the index of the outer loop.
+ std::vector<SymmetricTensor<2,dim> > phi_grads_u (dofs_per_cell);
+ std::vector<double> div_phi_u (dofs_per_cell);
+ std::vector<double> phi_p (dofs_per_cell);
+
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (unsigned int q=0; q<n_q_points; ++q)
{
+ for (unsigned int k=0; k<dofs_per_cell; ++k)
+ {
+ phi_grads_u[k] = fe_values[velocities].symmetric_gradient (k, q);
+ div_phi_u[k] = fe_values[velocities].divergence (k, q);
+ phi_p[k] = fe_values[pressure].value (k, q);
+ }
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
- const SymmetricTensor<2,dim>
- phi_i_grads_u = fe_values[velocities].symmetric_gradient (i, q);
- const Tensor<1,dim> phi_i_u = fe_values[velocities].value (i, q);
- const double div_phi_i_u = fe_values[velocities].divergence (i, q);
- const double phi_i_p = fe_values[pressure].value (i, q);
-
- for (unsigned int j=0; j<dofs_per_cell; ++j)
+ for (unsigned int j=0; j<=i; ++j)
{
- const SymmetricTensor<2,dim>
- phi_j_grads_u = fe_values[velocities].symmetric_gradient (j, q);
- const double div_phi_j_u = fe_values[velocities].divergence (j, q);
- const double phi_j_p = fe_values[pressure].value (j, q);
-
- local_matrix(i,j) += (phi_i_grads_u * phi_j_grads_u
- - div_phi_i_u * phi_j_p
- - phi_i_p * div_phi_j_u
- + phi_i_p * phi_j_p)
+ local_matrix(i,j) += (phi_grads_u[i] * phi_grads_u[j]
+ - div_phi_u[i] * phi_p[j]
+ - phi_p[i] * div_phi_u[j]
+ + phi_p[i] * phi_p[j])
* fe_values.JxW(q);
}
// terms constituting the pressure mass
// matrix are written into the correct
// position without any further
- // interaction:
+ // interaction. We have to be careful
+ // about one thing, though. We have
+ // only build up half of the local
+ // matrix because of symmetry, but
+ // we have to save the full system
+ // matrix in order to use the standard
+ // functions for the solution. Hence,
+ // we use the element below the
+ // diagonal if we happen to look
+ // above it.
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
- system_matrix.add (local_dof_indices[i],
- local_dof_indices[j],
- local_matrix(i,j));
+ {
+ if (j <= i)
+ system_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ local_matrix(i,j));
+ else
+ system_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ local_matrix(j,i));
+ }
for (unsigned int i=0; i<dofs_per_cell; ++i)
system_rhs(local_dof_indices[i]) += local_rhs(i);
// independently on the mesh size. This
// is precisely what we do here: We
// choose another ILU preconditioner
- // and take it along to the
+ // and take it along to the
// InverseMatrix object via the
// corresponding template parameter. A
// CG solver is then called within the
// vmult operation of the inverse matrix.
- //
- // An alternative that is cheaper to build,
- // but needs more iterations afterwards,
- // would be to choose a SSOR preconditioner
- // with factor 1.2. It needs about twice
- // the number of iterations, but the costs
- // for its generation are almost neglible.
+ //
+ // An alternative that is cheaper to build,
+ // but needs more iterations afterwards,
+ // would be to choose a SSOR preconditioner
+ // with factor 1.2. It needs about twice
+ // the number of iterations, but the costs
+ // for its generation are almost neglible.
SparseILU<double> preconditioner;
preconditioner.initialize (system_matrix.block(1,1),
SparseILU<double>::AdditionalData());
projects / dealii-svn.git / commitdiff