]> https://gitweb.dealii.org/ - dealii.git/commitdiff
A faster variant of assembly for step-22, added remark to this in step-20, and a...
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 19 Mar 2008 16:41:06 +0000 (16:41 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 19 Mar 2008 16:41:06 +0000 (16:41 +0000)
git-svn-id: https://svn.dealii.org/trunk@15910 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-20/doc/intro.dox
deal.II/examples/step-22/doc/results.dox
deal.II/examples/step-22/step-22.cc

index 595fa6a884ab703541586b59ac31087a036c340d..412d9e18fe48428f78428de16d62ac9445614ad2 100644 (file)
@@ -246,7 +246,10 @@ component. The whole mechanism is described in more detail in the
 
 In practice, it turns out that we can do a bit better if we evaluate the shape
 functions, their gradients and divergences only once per outermost loop, and
-store the result, as this saves us a few otherwise repeated computations.
+store the result, as this saves us a few otherwise repeated computations (it is
+possible to save even more repeated operations by calculating all relevant
+quantities in advance and then only inserting the results in the actual loop,
+see step-22 for a realization of that approach).
 The final result then looks like this, working in every space dimension:
 
 @code
index 94f130ea8ac0bcbfc19d3459dc988f87230c1107..b26b9e1f28be20540f822b16024b242e2f658f38 100644 (file)
@@ -352,6 +352,27 @@ demands on the linear algebra than standard elliptic problems as we have seen
 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 
@@ -592,73 +613,73 @@ 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.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
@@ -669,74 +690,58 @@ there won't be any 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.
 
-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
@@ -746,6 +751,11 @@ problem size than CG (as explained above).
 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
index 44b9334180222f68615286d15ee23a62a8d3f001..c5c70a1ec710eef3974a140a72f9e8b55ec3798e 100644 (file)
@@ -517,19 +517,15 @@ StokesProblem<dim>::StokesProblem (const unsigned int degree)
                                 // 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
@@ -564,7 +560,7 @@ void StokesProblem<dim>::setup_dofs ()
   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;
@@ -774,8 +770,53 @@ void StokesProblem<dim>::assemble_system ()
   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();
@@ -790,25 +831,20 @@ void StokesProblem<dim>::assemble_system ()
       
       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);     
 
                }
@@ -854,14 +890,30 @@ void StokesProblem<dim>::assemble_system ()
                                       // 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);
@@ -1017,18 +1069,18 @@ void StokesProblem<dim>::solve ()
                                     // 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());

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.