]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Block Schur preconditioned solver on u and p using the BlockMatrixArray class.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 16 Apr 2008 06:43:10 +0000 (06:43 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 16 Apr 2008 06:43:10 +0000 (06:43 +0000)
git-svn-id: https://svn.dealii.org/trunk@15997 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-31/step-31.cc

index 5bf011105f61c3835172febdbeb24e0f622d73b7..033d8a4830bbf026b459d2a94ea81df7a16bad33 100644 (file)
 #include <lac/block_vector.h>
 #include <lac/full_matrix.h>
 #include <lac/block_sparse_matrix.h>
+#include <lac/solver_gmres.h>
 #include <lac/solver_cg.h>
 #include <lac/precondition.h>
 #include <lac/sparse_direct.h>
+#include <lac/sparse_ilu.h>
+#include <lac/block_matrix_array.h>
 
 #include <grid/tria.h>
 #include <grid/grid_generator.h>
 using namespace dealii;
 
 
-                                 
 template <int dim>
-class BoussinesqFlowProblem 
+ struct InnerPreconditioner;
+
+ template <>
+ struct InnerPreconditioner<2>
+ {
+     typedef SparseDirectUMFPACK type;
+ };
+
+ template <>
+ struct InnerPreconditioner<3>
+ {
+     typedef SparseILU<double> type;
+ };
+
+
+template <int dim>
+class BoussinesqFlowProblem
 {
   public:
     BoussinesqFlowProblem (const unsigned int degree);
     void run ();
-    
+
   private:
     void setup_dofs (const bool setup_matrices);
     void assemble_system ();
@@ -73,26 +91,27 @@ class BoussinesqFlowProblem
     void solve ();
     void output_results () const;
     void refine_mesh ();
-    
+
     const unsigned int   degree;
-    
+
     Triangulation<dim>   triangulation;
     FESystem<dim>        fe;
     DoFHandler<dim>      dof_handler;
 
     ConstraintMatrix     hanging_node_constraints;
-    
+
     BlockSparsityPattern      sparsity_pattern;
     BlockSparseMatrix<double> system_matrix;
 
     double time_step;
     unsigned int timestep_number;
+
     BlockVector<double> solution;
     BlockVector<double> old_solution;
     BlockVector<double> system_rhs;
 
-    boost::shared_ptr<SparseDirectUMFPACK> A_preconditioner;
+    boost::shared_ptr<typename InnerPreconditioner<dim>::type> A_preconditioner;
+    boost::shared_ptr<SparseILU<double> > Mp_preconditioner;
 
     bool rebuild_matrices;
     bool rebuild_preconditioner;
@@ -103,11 +122,11 @@ class BoussinesqFlowProblem
 
 
 template <int dim>
-class PressureBoundaryValues : public Function<dim> 
+class PressureBoundaryValues : public Function<dim>
 {
   public:
     PressureBoundaryValues () : Function<dim>(1) {}
-    
+
     virtual double value (const Point<dim>   &p,
                           const unsigned int  component = 0) const;
 };
@@ -116,7 +135,7 @@ class PressureBoundaryValues : public Function<dim>
 template <int dim>
 double
 PressureBoundaryValues<dim>::value (const Point<dim>  &/*p*/,
-                                    const unsigned int /*component*/) const 
+                                    const unsigned int /*component*/) const
 {
   return 0;
 }
@@ -124,11 +143,11 @@ PressureBoundaryValues<dim>::value (const Point<dim>  &/*p*/,
 
 
 template <int dim>
-class TemperatureBoundaryValues : public Function<dim> 
+class TemperatureBoundaryValues : public Function<dim>
 {
   public:
     TemperatureBoundaryValues () : Function<dim>(1) {}
-    
+
     virtual double value (const Point<dim>   &p,
                           const unsigned int  component = 0) const;
 };
@@ -138,7 +157,7 @@ class TemperatureBoundaryValues : public Function<dim>
 template <int dim>
 double
 TemperatureBoundaryValues<dim>::value (const Point<dim> &p,
-                                      const unsigned int /*component*/) const 
+                                      const unsigned int /*component*/) const
 {
 //TODO: leftover from olden times. replace by something sensible once we have
 //diffusion in the temperature field
@@ -152,15 +171,15 @@ TemperatureBoundaryValues<dim>::value (const Point<dim> &p,
 
 
 template <int dim>
-class InitialValues : public Function<dim> 
+class InitialValues : public Function<dim>
 {
   public:
     InitialValues () : Function<dim>(dim+2) {}
-    
+
     virtual double value (const Point<dim>   &p,
                           const unsigned int  component = 0) const;
 
-    virtual void vector_value (const Point<dim> &p, 
+    virtual void vector_value (const Point<dim> &p,
                                Vector<double>   &value) const;
 };
 
@@ -168,7 +187,7 @@ class InitialValues : public Function<dim>
 template <int dim>
 double
 InitialValues<dim>::value (const Point<dim>  &,
-                           const unsigned int) const 
+                           const unsigned int) const
 {
   return 0;
 }
@@ -177,7 +196,7 @@ InitialValues<dim>::value (const Point<dim>  &,
 template <int dim>
 void
 InitialValues<dim>::vector_value (const Point<dim> &p,
-                                  Vector<double>   &values) const 
+                                  Vector<double>   &values) const
 {
   for (unsigned int c=0; c<this->n_components; ++c)
     values(c) = InitialValues<dim>::value (p, c);
@@ -186,15 +205,15 @@ InitialValues<dim>::vector_value (const Point<dim> &p,
 
 
 template <int dim>
-class RightHandSide : public Function<dim> 
+class RightHandSide : public Function<dim>
 {
   public:
     RightHandSide () : Function<dim>(dim+2) {}
-    
+
     virtual double value (const Point<dim>   &p,
                           const unsigned int  component = 0) const;
 
-    virtual void vector_value (const Point<dim> &p, 
+    virtual void vector_value (const Point<dim> &p,
                                Vector<double>   &value) const;
 };
 
@@ -202,7 +221,7 @@ class RightHandSide : public Function<dim>
 template <int dim>
 double
 RightHandSide<dim>::value (const Point<dim>  &p,
-                           const unsigned int component) const 
+                           const unsigned int component) const
 {
   if (component == dim+1)
     return ((p.distance (Point<dim>(.3,.1)) < 1./32)
@@ -222,7 +241,7 @@ RightHandSide<dim>::value (const Point<dim>  &p,
 template <int dim>
 void
 RightHandSide<dim>::vector_value (const Point<dim> &p,
-                                  Vector<double>   &values) const 
+                                  Vector<double>   &values) const
 {
   for (unsigned int c=0; c<this->n_components; ++c)
     values(c) = RightHandSide<dim>::value (p, c);
@@ -256,7 +275,7 @@ InverseMatrix<Matrix,Preconditioner>::InverseMatrix (const Matrix &m,
 {}
 
 
-                                 
+
 template <class Matrix, class Preconditioner>
 void InverseMatrix<Matrix,Preconditioner>::vmult (Vector<double>       &dst,
                                                  const Vector<double> &src) const
@@ -277,45 +296,67 @@ void InverseMatrix<Matrix,Preconditioner>::vmult (Vector<double>       &dst,
 }
 
 
-
-template <class Preconditioner>
-class SchurComplement : public Subscriptor
+                               // This is the implementation
+                               // of the Schur complement
+                               // preconditioner as described
+                               // in the section on improved
+                               // solvers in step-22. See there
+                               // for more explanation of the
+                               // method.
+template <class PreconditionerA, class PreconditionerMp>
+class BlockSchurPreconditioner : public Subscriptor
 {
   public:
-    SchurComplement (const BlockSparseMatrix<double> &A,
-                     const InverseMatrix<SparseMatrix<double>,Preconditioner> &Minv);
+    BlockSchurPreconditioner (const BlockSparseMatrix<double>         &S,
+          const InverseMatrix<SparseMatrix<double>,PreconditionerMp>  &Mpinv,
+          const PreconditionerA &Apreconditioner);
 
-    void vmult (Vector<double>       &dst,
-                const Vector<double> &src) const;
+  void vmult (BlockVector<double>       &dst,
+              const BlockVector<double> &src) const;
 
   private:
     const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
-    const SmartPointer<const InverseMatrix<SparseMatrix<double>,Preconditioner> > m_inverse;
-    
-    mutable Vector<double> tmp1, tmp2;
-};
+    const SmartPointer<const InverseMatrix<SparseMatrix<double>,
+                       PreconditionerMp > > m_inverse;
+    const PreconditionerA &a_preconditioner;
 
+    mutable Vector<double> tmp;
 
+};
 
-template <class Preconditioner>
-SchurComplement<Preconditioner>::
-SchurComplement (const BlockSparseMatrix<double> &A,
-                 const InverseMatrix<SparseMatrix<double>,Preconditioner> &Minv)
+template <class PreconditionerA, class PreconditionerMp>
+BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::BlockSchurPreconditioner(
+          const BlockSparseMatrix<double>                            &S,
+          const InverseMatrix<SparseMatrix<double>,PreconditionerMp> &Mpinv,
+          const PreconditionerA &Apreconditioner
+          )
                 :
-                system_matrix (&A),
-                m_inverse (&Minv),
-                tmp1 (A.block(0,0).m()),
-                tmp2 (A.block(0,0).m())
-{}
-
+                system_matrix           (&S),
+                m_inverse               (&Mpinv),
+                a_preconditioner        (Apreconditioner),
+                tmp                     (S.block(1,1).m())
+{
+}
 
-template <class Preconditioner>
-void SchurComplement<Preconditioner>::vmult (Vector<double>       &dst,
-                                            const Vector<double> &src) const
+        // Now the interesting function, the multiplication of
+        // the preconditioner with a BlockVector.
+template <class PreconditionerA, class PreconditionerMp>
+void BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::vmult (
+                                     BlockVector<double>       &dst,
+                                     const BlockVector<double> &src) const
 {
-  system_matrix->block(0,1).vmult (tmp1, src);
-  m_inverse->vmult (tmp2, tmp1);
-  system_matrix->block(1,0).vmult (dst, tmp2);
+        // Form u_new = A^{-1} u
+  a_preconditioner.vmult (dst.block(0), src.block(0));
+        // Form tmp = - B u_new + p
+        // (<code>SparseMatrix::residual</code>
+        // does precisely this)
+  system_matrix->block(1,0).residual(tmp,
+                                     dst.block(0), src.block(1));
+        // Change sign in tmp.block(1)
+  tmp *= -1;
+        // Multiply by approximate Schur complement
+        // (i.e. a pressure mass matrix)
+  m_inverse->vmult (dst.block(1), tmp);
 }
 
 
@@ -338,8 +379,8 @@ BoussinesqFlowProblem<dim>::BoussinesqFlowProblem (const unsigned int degree)
 
 template <int dim>
 void BoussinesqFlowProblem<dim>::setup_dofs (const bool setup_matrices)
-{  
-  dof_handler.distribute_dofs (fe); 
+{
+  dof_handler.distribute_dofs (fe);
   DoFRenumbering::component_wise (dof_handler);
 
   hanging_node_constraints.clear ();
@@ -353,7 +394,7 @@ void BoussinesqFlowProblem<dim>::setup_dofs (const bool setup_matrices)
   hanging_node_constraints.close ();
 
   std::vector<unsigned int> dofs_per_component (dim+2);
-  DoFTools::count_dofs_per_component (dof_handler, dofs_per_component);  
+  DoFTools::count_dofs_per_component (dof_handler, dofs_per_component);
   const unsigned int n_u = dofs_per_component[0] * dim,
                      n_p = dofs_per_component[dim],
                      n_T = dofs_per_component[dim+1];
@@ -366,14 +407,14 @@ void BoussinesqFlowProblem<dim>::setup_dofs (const bool setup_matrices)
             << " (" << n_u << '+' << n_p << '+'<< n_T <<')'
             << std::endl
             << std::endl;
-  
+
   const unsigned int
     n_couplings = dof_handler.max_couplings_between_dofs();
 
   if (setup_matrices == true)
     {
       system_matrix.clear ();
-      
+
       sparsity_pattern.reinit (3,3);
       sparsity_pattern.block(0,0).reinit (n_u, n_u, n_couplings);
       sparsity_pattern.block(1,0).reinit (n_p, n_u, n_couplings);
@@ -384,29 +425,29 @@ void BoussinesqFlowProblem<dim>::setup_dofs (const bool setup_matrices)
       sparsity_pattern.block(0,2).reinit (n_u, n_T, n_couplings);
       sparsity_pattern.block(1,2).reinit (n_p, n_T, n_couplings);
       sparsity_pattern.block(2,2).reinit (n_T, n_T, n_couplings);
-  
-      sparsity_pattern.collect_sizes();    
 
-  
+      sparsity_pattern.collect_sizes();
+
+
       DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
       hanging_node_constraints.condense (sparsity_pattern);
       sparsity_pattern.compress();
-  
+
       system_matrix.reinit (sparsity_pattern);
     }
-                                   
+
   solution.reinit (3);
   solution.block(0).reinit (n_u);
   solution.block(1).reinit (n_p);
   solution.block(2).reinit (n_T);
   solution.collect_sizes ();
-  
+
   old_solution.reinit (3);
   old_solution.block(0).reinit (n_u);
   old_solution.block(1).reinit (n_p);
   old_solution.block(2).reinit (n_T);
   old_solution.collect_sizes ();
-  
+
   system_rhs.reinit (3);
   system_rhs.block(0).reinit (n_u);
   system_rhs.block(1).reinit (n_p);
@@ -417,14 +458,14 @@ void BoussinesqFlowProblem<dim>::setup_dofs (const bool setup_matrices)
 
 
 template <int dim>
-void BoussinesqFlowProblem<dim>::assemble_system () 
+void BoussinesqFlowProblem<dim>::assemble_system ()
 {
   if (rebuild_matrices == true)
     system_matrix=0;
 
   system_rhs=0;
-  
-  QGauss<dim>   quadrature_formula(degree+2); 
+
+  QGauss<dim>   quadrature_formula(degree+2);
   QGauss<dim-1> face_quadrature_formula(degree+2);
 
   FEValues<dim> fe_values (fe, quadrature_formula,
@@ -436,12 +477,12 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
                            update_gradients
                            :
                            UpdateFlags(0)));
-  FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula, 
+  FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
                                     update_values    | update_normal_vectors |
                                     update_quadrature_points  | update_JxW_values);
 
   const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
-  
+
   const unsigned int   n_q_points      = quadrature_formula.size();
   const unsigned int   n_face_q_points = face_quadrature_formula.size();
 
@@ -449,17 +490,17 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
   Vector<double>       local_rhs (dofs_per_cell);
 
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-  
+
   const PressureBoundaryValues<dim> pressure_boundary_values;
-  
+
   std::vector<double>               boundary_values (n_face_q_points);
-  
+
   std::vector<Vector<double> >      old_solution_values(n_q_points, Vector<double>(dim+2));
   std::vector<std::vector<Tensor<1,dim> > >  old_solution_grads(n_q_points,
                                                                 std::vector<Tensor<1,dim> > (dim+2));
 
   const double Raleigh_number = 10;
-  
+
   const FEValuesExtractors::Vector velocities (0);
   const FEValuesExtractors::Scalar pressure (dim);
   const FEValuesExtractors::Scalar temperature (dim+1);
@@ -468,7 +509,7 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
   for (; cell!=endc; ++cell)
-    { 
+    {
       fe_values.reinit (cell);
       local_matrix = 0;
       local_rhs = 0;
@@ -478,7 +519,7 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
       for (unsigned int q=0; q<n_q_points; ++q)
        {
          const double old_temperature = old_solution_values[q](dim+1);
-         
+
          for (unsigned int i=0; i<dofs_per_cell; ++i)
            {
 
@@ -491,9 +532,9 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
                  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);
-                 const double        phi_i_T       = fe_values[temperature].value (i, q); 
+                 const double        phi_i_T       = fe_values[temperature].value (i, q);
                  const Tensor<1,dim> grad_phi_i_T  = fe_values[temperature].gradient(i, q);
-            
+
                  for (unsigned int j=0; j<dofs_per_cell; ++j)
                    {
                      const SymmetricTensor<2,dim>
@@ -501,24 +542,24 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
                      const double        div_phi_j_u = fe_values[velocities].divergence (j, q);
                      const double        phi_j_p     = fe_values[pressure].value (j, q);
                      const double        phi_j_T     = fe_values[temperature].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
                                            + phi_i_T * phi_j_T)
-                                          * fe_values.JxW(q);     
+                                          * fe_values.JxW(q);
                    }
                }
-             
+
              const Point<dim> gravity (0,1);
-             
+
              local_rhs(i) += (Raleigh_number *
                               gravity * phi_i_u * old_temperature)*
                              fe_values.JxW(q);
           }
        }
-      
+
 
       for (unsigned int face_no=0;
            face_no<GeometryInfo<dim>::faces_per_cell;
@@ -526,12 +567,12 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
         if (cell->at_boundary(face_no))
           {
             fe_face_values.reinit (cell, face_no);
-            
+
             pressure_boundary_values
               .value_list (fe_face_values.get_quadrature_points(),
                            boundary_values);
 
-            for (unsigned int q=0; q<n_face_q_points; ++q) 
+            for (unsigned int q=0; q<n_face_q_points; ++q)
               for (unsigned int i=0; i<dofs_per_cell; ++i)
                 {
                   const Tensor<1,dim>
@@ -554,16 +595,16 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
                                 local_dof_indices[j],
                                 local_matrix(i,j));
        }
-      
+
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         system_rhs(local_dof_indices[i]) += local_rhs(i);
     }
 
   if (rebuild_matrices == true)
     hanging_node_constraints.condense (system_matrix);
-  
-  hanging_node_constraints.condense (system_rhs);  
-  
+
+  hanging_node_constraints.condense (system_rhs);
+
   if (rebuild_matrices == true)
     {
 //       std::map<unsigned int,double> boundary_values;
@@ -581,11 +622,11 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
 //         {
 //           std::cout << "Found cell and vertex: " << cell << ' '
 //                     << v << std::endl;
-             
+
 //           boundary_values[cell->vertex_dof_index(v,0)] = 0;
 //           break;
-//         } 
-      
+//         }
+
 //      std::vector<bool> component_mask (dim+2, true);
 //       component_mask[dim] = component_mask[dim+1] = false;
 //       VectorTools::interpolate_boundary_values (dof_handler,
@@ -597,7 +638,7 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
 //       MatrixTools::apply_boundary_values (boundary_values,
 //                                       system_matrix,
 //                                       solution,
-//                                       system_rhs);  
+//                                       system_rhs);
     }
 
   if (rebuild_preconditioner == true)
@@ -605,15 +646,29 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
       Assert (rebuild_matrices == true,
              ExcMessage ("There is no point in rebuilding the preconditioner "
                          "without a rebuilt matrix!"));
-      
+
       std::cout << "   Rebuilding preconditioner..." << std::flush;
-      
+
+                               // Rebuild the preconditioner
+                               // for the velocity-velocity
+                               // block (0,0)
       A_preconditioner
-       = boost::shared_ptr<SparseDirectUMFPACK>(new SparseDirectUMFPACK());
-      A_preconditioner->initialize (system_matrix.block(0,0));
+       = boost::shared_ptr<typename InnerPreconditioner<dim>::type>
+               (new typename InnerPreconditioner<dim>::type());
+      A_preconditioner->initialize (system_matrix.block(0,0),
+               typename InnerPreconditioner<dim>::type::AdditionalData());
+
+                               // Rebuild the preconditioner
+                               // for the pressure-pressure
+                               // block (1,1)
+      Mp_preconditioner
+       = boost::shared_ptr<SparseILU<double> >
+               (new SparseILU<double>);
+      Mp_preconditioner->initialize (system_matrix.block(0,0),
+                                    SparseILU<double>::AdditionalData());
 
       std::cout << std::endl;
-      
+
       rebuild_preconditioner = false;
     }
 
@@ -626,28 +681,28 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
 
 
 template <int dim>
-void BoussinesqFlowProblem<dim>::assemble_rhs_T () 
-{  
-  QGauss<dim>   quadrature_formula(degree+2); 
+void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
+{
+  QGauss<dim>   quadrature_formula(degree+2);
   QGauss<dim-1> face_quadrature_formula(degree+2);
-  FEValues<dim> fe_values (fe, quadrature_formula, 
+  FEValues<dim> fe_values (fe, quadrature_formula,
                            update_values    | update_gradients |
                            update_quadrature_points  | update_JxW_values);
-  FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula, 
+  FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
                                     update_values    | update_normal_vectors |
                                     update_quadrature_points  | update_JxW_values);
-  FESubfaceValues<dim> fe_subface_values (fe, face_quadrature_formula, 
+  FESubfaceValues<dim> fe_subface_values (fe, face_quadrature_formula,
                                          update_values    | update_normal_vectors |
                                          update_JxW_values);
-  FEFaceValues<dim> fe_face_values_neighbor (fe, face_quadrature_formula, 
+  FEFaceValues<dim> fe_face_values_neighbor (fe, face_quadrature_formula,
                                              update_values);
-  FESubfaceValues<dim> fe_subface_values_neighbor (fe, face_quadrature_formula, 
+  FESubfaceValues<dim> fe_subface_values_neighbor (fe, face_quadrature_formula,
                                                   update_values);
+
   const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
   const unsigned int   n_q_points      = quadrature_formula.size();
   const unsigned int   n_face_q_points = face_quadrature_formula.size();
-  
+
   Vector<double>       local_rhs (dofs_per_cell);
 
   std::vector<Vector<double> > old_solution_values(n_q_points, Vector<double>(dim+2));
@@ -660,14 +715,14 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
   std::vector<std::vector<Tensor<1,dim> > >
     present_solution_grads(n_q_points,
                           std::vector<Tensor<1,dim> >(dim+2));
-         
-  
+
+
   std::vector<double> neighbor_temperature (n_face_q_points);
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
   TemperatureBoundaryValues<dim> temperature_boundary_values;
   const FEValuesExtractors::Scalar temperature (dim+1);
-  
+
   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
@@ -681,7 +736,7 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
       fe_values.get_function_values (solution, present_solution_values);
       fe_values.get_function_gradients (solution, present_solution_grads);
 
-      for (unsigned int q=0; q<n_q_points; ++q) 
+      for (unsigned int q=0; q<n_q_points; ++q)
         for (unsigned int i=0; i<dofs_per_cell; ++i)
           {
             const double old_T = old_solution_values[q](dim+1);
@@ -691,13 +746,13 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
 
            double present_div_u = 0;
             for (unsigned int d=0; d<dim; ++d)
-              present_div_u += present_solution_grads[q][d][d];            
-           
+              present_div_u += present_solution_grads[q][d][d];
+
             const double        phi_i_T      = fe_values[temperature].value (i, q);
             const Tensor<1,dim> grad_phi_i_T = fe_values[temperature].gradient (i, q);
 
            const Point<dim> p = fe_values.quadrature_point(q);
-           
+
             local_rhs(i) += (time_step *
                              old_T *
                              (present_u *
@@ -715,8 +770,8 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
                             fe_values.JxW(q);
           }
 
-      
-//TODO: unify the code that actually does the assembly down below      
+
+//TODO: unify the code that actually does the assembly down below
       for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
            ++face_no)
         if (cell->at_boundary(face_no)
@@ -746,14 +801,14 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
              {
                const typename DoFHandler<dim>::active_cell_iterator
                  neighbor = cell->neighbor(face_no);
-               
+
                fe_face_values_neighbor.reinit (neighbor,
                                                cell->neighbor_of_neighbor(face_no));
-             
+
                fe_face_values_neighbor
                  .get_function_values (old_solution,
                                        old_solution_values_face_neighbor);
-             
+
                for (unsigned int q=0; q<n_face_q_points; ++q)
                  neighbor_temperature[q] = old_solution_values_face_neighbor[q](dim+1);
              }
@@ -768,7 +823,7 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
                                           fe_face_values.normal_vector(q);
 
                const bool is_outflow_q_point = (normal_flux >= 0);
-                                     
+
                for (unsigned int i=0; i<dofs_per_cell; ++i)
                  local_rhs(i) -= time_step *
                                  normal_flux *
@@ -803,11 +858,11 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
 
                  fe_face_values_neighbor.reinit (neighbor,
                                                  cell->neighbor_of_neighbor(face_no));
-             
+
                  fe_face_values_neighbor
                    .get_function_values (old_solution,
                                          old_solution_values_face_neighbor);
-                 
+
                  for (unsigned int q=0; q<n_face_q_points; ++q)
                    neighbor_temperature[q] = old_solution_values_face_neighbor[q](dim+1);
 
@@ -821,7 +876,7 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
                                                 fe_subface_values.normal_vector(q);
 
                      const bool is_outflow_q_point = (normal_flux >= 0);
-                                     
+
                      for (unsigned int i=0; i<dofs_per_cell; ++i)
                        local_rhs(i) -= time_step *
                                        normal_flux *
@@ -863,11 +918,11 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
              fe_subface_values_neighbor.reinit (neighbor,
                                                 neighbor_face_no,
                                                 neighbor_subface_no);
-             
+
              fe_subface_values_neighbor
                .get_function_values (old_solution,
                                      old_solution_values_face_neighbor);
-             
+
              for (unsigned int q=0; q<n_face_q_points; ++q)
                neighbor_temperature[q] = old_solution_values_face_neighbor[q](dim+1);
 
@@ -881,7 +936,7 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
                                             fe_face_values.normal_vector(q);
 
                  const bool is_outflow_q_point = (normal_flux >= 0);
-                                     
+
                  for (unsigned int i=0; i<dofs_per_cell; ++i)
                    local_rhs(i) -= time_step *
                                    normal_flux *
@@ -893,90 +948,87 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
                                    fe_face_values[temperature].value (i,q) *
                                    fe_face_values.JxW(q);
                }
-           }      
-      
+           }
+
       cell->get_dof_indices (local_dof_indices);
       for (unsigned int i=0; i<dofs_per_cell; ++i)
-        system_rhs(local_dof_indices[i]) += local_rhs(i);       
+        system_rhs(local_dof_indices[i]) += local_rhs(i);
     }
-} 
+}
 
 
 
 
 template <int dim>
-void BoussinesqFlowProblem<dim>::solve () 
+void BoussinesqFlowProblem<dim>::solve ()
 {
   solution = old_solution;
-  
-  const InverseMatrix<SparseMatrix<double>,SparseDirectUMFPACK>
-    A_inverse (system_matrix.block(0,0), *A_preconditioner);
-  Vector<double> tmp (solution.block(0).size());
-  Vector<double> schur_rhs (solution.block(1).size());
-  Vector<double> tmp2 (solution.block(2).size());
-  
+
+                               // Use the BlockMatrixArray structure
+                               // for extracting only the upper left
+                               // 2x2 blocks from the matrix that will
+                               // be used for the solution of the
+                               // blocked system.
   {
-    A_inverse.vmult (tmp, system_rhs.block(0));
-    system_matrix.block(1,0).vmult (schur_rhs, tmp);
-    schur_rhs -= system_rhs.block(1);
-
-    
-    SchurComplement<SparseDirectUMFPACK>
-      schur_complement (system_matrix, A_inverse);
-    
-    SolverControl solver_control (system_matrix.block(0,0).m(),
-                                  1e-6*schur_rhs.l2_norm());
-    SolverCG<>    cg (solver_control);
-    
-    PreconditionSSOR<> preconditioner;
-    preconditioner.initialize (system_matrix.block(1,1), 1.2);
-
-    InverseMatrix<SparseMatrix<double>,PreconditionSSOR<> >
-      m_inverse (system_matrix.block(1,1), preconditioner);
-    
-    try
-      {
-       cg.solve (schur_complement, solution.block(1), schur_rhs,
-                 m_inverse);
-      }
-    catch (...)
-      {
-       abort ();
-      }
+    GrowingVectorMemory<Vector<double> > simple_mem;
+    BlockMatrixArray<double> stokes_submatrix(2, 2, simple_mem);
 
-                                    // produce a consistent flow field
-    hanging_node_constraints.distribute (solution);
-  
-    std::cout << "   "
-              << solver_control.last_step()
-              << " CG Schur complement iterations for pressure."
-              << std::endl;    
-  }
+    stokes_submatrix.enter(system_matrix.block(0,0),0,0);
+    stokes_submatrix.enter(system_matrix.block(0,1),0,1);
+    stokes_submatrix.enter(system_matrix.block(1,0),1,0);
 
-  {
-    system_matrix.block(0,1).vmult (tmp, solution.block(1));
-    tmp *= -1;
-    tmp += system_rhs.block(0);
-    
-    A_inverse.vmult (solution.block(0), tmp);
+                               // Define some temporary vectors
+                               // for the solution process.
+    std::vector<unsigned int> block_sizes(2);
+    block_sizes[0] = solution.block(0).size();
+    block_sizes[1] = solution.block(1).size();
 
-                                    // produce a consistent pressure field
-    hanging_node_constraints.distribute (solution);
-  }
+    BlockVector<double> up_rhs(block_sizes);
+    BlockVector<double> up(block_sizes);
+
+    up_rhs.block(0) = system_rhs.block(0);
+    up_rhs.block(1) = system_rhs.block(1);
+
+                               // Set up inverse matrix for
+                               // pressure mass matrix
+    InverseMatrix<SparseMatrix<double>,SparseILU<double> >
+      mp_inverse (system_matrix.block(1,1), *Mp_preconditioner);
+
+                               // Set up block Schur preconditioner
+    BlockSchurPreconditioner<typename InnerPreconditioner<dim>::type,
+                           SparseILU<double> >
+      preconditioner (system_matrix, mp_inverse, *A_preconditioner);
+
+                               // Set up GMRES solver and
+                               // solve.
+    SolverControl solver_control (system_matrix.m(),
+                                    1e-6*system_rhs.l2_norm());
 
+    SolverGMRES<BlockVector<double> > gmres(solver_control,
+                       SolverGMRES<BlockVector<double> >::AdditionalData(100));
+
+    gmres.solve(stokes_submatrix, up, up_rhs,
+                preconditioner);
+
+                               // Produce a constistent solution field
+    hanging_node_constraints.distribute (up);
+
+    solution.block(0) = up.block(0);
+    solution.block(1) = up.block(1);
+  }
                                   // for DGQ1 needs to be /15
   time_step = GridTools::minimal_cell_diameter(triangulation) /
               std::max (get_maximal_velocity(), .05) / 2;
 
   assemble_rhs_T ();
   {
-    
+
     SolverControl solver_control (system_matrix.block(2,2).m(),
                                   1e-8*system_rhs.block(2).l2_norm());
     SolverCG<>   cg (solver_control);
     PreconditionJacobi<> preconditioner;
     preconditioner.initialize (system_matrix.block(2,2));
-    
+
     try
       {
        cg.solve (system_matrix.block(2,2), solution.block(2), system_rhs.block(2),
@@ -986,10 +1038,10 @@ void BoussinesqFlowProblem<dim>::solve ()
       {
        abort ();
       }
-       
+
                                     // produce a consistent temperature field
     hanging_node_constraints.distribute (solution);
-                
+
     std::cout << "   "
               << solver_control.last_step()
               << " CG iterations for temperature."
@@ -998,9 +1050,9 @@ void BoussinesqFlowProblem<dim>::solve ()
              << *std::max_element (solution.block(2).begin(),
                                    solution.block(2).end())
              << std::endl;
-  } 
+  }
 }
-                                 
+
 
 
 template <int dim>
@@ -1008,11 +1060,11 @@ void BoussinesqFlowProblem<dim>::output_results ()  const
 {
   if (timestep_number % 10 != 0)
     return;
-  
+
   std::vector<std::string> solution_names (dim, "velocity");
   solution_names.push_back ("p");
   solution_names.push_back ("T");
-  
+
   DataOut<dim> data_out;
 
   data_out.attach_dof_handler (dof_handler);
@@ -1023,13 +1075,13 @@ void BoussinesqFlowProblem<dim>::output_results ()  const
   for (unsigned int i=0; i<dim; ++i)
     data_component_interpretation[i]
       = DataComponentInterpretation::component_is_part_of_vector;
-  
+
   data_out.add_data_vector (solution, solution_names,
                            DataOut<dim>::type_dof_data,
                            data_component_interpretation);
-  
+
   data_out.build_patches ();
-  
+
   std::ostringstream filename;
   filename << "solution-" << Utilities::int_to_string(timestep_number, 4) << ".vtk";
 
@@ -1041,11 +1093,11 @@ void BoussinesqFlowProblem<dim>::output_results ()  const
 
 template <int dim>
 void
-BoussinesqFlowProblem<dim>::refine_mesh () 
+BoussinesqFlowProblem<dim>::refine_mesh ()
 {
   Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
 
-//TODO do this better      
+//TODO do this better
   DerivativeApproximation::approximate_gradient (dof_handler,
                                                 old_solution,
                                                 estimated_error_per_cell,
@@ -1069,7 +1121,7 @@ BoussinesqFlowProblem<dim>::refine_mesh ()
 
   Vector<double> x_old_solution (dof_handler.n_dofs());
   x_old_solution = old_solution;
-  
+
   soltrans.prepare_for_coarsening_and_refinement(x_old_solution);
 
   triangulation.execute_coarsening_and_refinement ();
@@ -1090,16 +1142,16 @@ template <int dim>
 double
 BoussinesqFlowProblem<dim>::get_maximal_velocity () const
 {
-  QGauss<dim>   quadrature_formula(degree+2); 
+  QGauss<dim>   quadrature_formula(degree+2);
   const unsigned int   n_q_points
     = quadrature_formula.size();
 
-  FEValues<dim> fe_values (fe, quadrature_formula, 
+  FEValues<dim> fe_values (fe, quadrature_formula,
                            update_values);
   std::vector<Vector<double> > solution_values(n_q_points,
                                                Vector<double>(dim+2));
   double max_velocity = 0;
-  
+
   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
@@ -1112,20 +1164,20 @@ BoussinesqFlowProblem<dim>::get_maximal_velocity () const
         {
           Tensor<1,dim> velocity;
           for (unsigned int i=0; i<dim; ++i)
-            velocity[i] = solution_values[q](i);          
-          
+            velocity[i] = solution_values[q](i);
+
           max_velocity = std::max (max_velocity,
                                    velocity.norm());
         }
     }
-  
+
   return max_velocity;
 }
 
 
 
 template <int dim>
-void BoussinesqFlowProblem<dim>::run () 
+void BoussinesqFlowProblem<dim>::run ()
 {
   switch (dim)
     {
@@ -1137,7 +1189,7 @@ void BoussinesqFlowProblem<dim>::run ()
 //     triangulation.set_boundary (0, boundary);
 
        GridGenerator::hyper_cube (triangulation);
-       
+
        triangulation.refine_global (6);
 
        break;
@@ -1147,10 +1199,10 @@ void BoussinesqFlowProblem<dim>::run ()
       {
        GridGenerator::hyper_shell (triangulation,
                                    Point<dim>(), 0.5, 1.0);
-       
+
        static HyperShellBoundary<dim> boundary;
        triangulation.set_boundary (0, boundary);
-       
+
        triangulation.refine_global (2);
 
        break;
@@ -1159,7 +1211,7 @@ void BoussinesqFlowProblem<dim>::run ()
       default:
            Assert (false, ExcNotImplemented());
     }
-  
+
 
   const bool do_adaptivity = false;
 
@@ -1172,7 +1224,7 @@ void BoussinesqFlowProblem<dim>::run ()
                            QGauss<dim>(degree+2),
                            InitialValues<dim>(),
                            old_solution);
-  
+
       for (unsigned int pre_refinement=0; pre_refinement<4-dim; ++pre_refinement)
        {
          refine_mesh ();
@@ -1194,29 +1246,29 @@ void BoussinesqFlowProblem<dim>::run ()
                            InitialValues<dim>(),
                            old_solution);
     }
-  
+
   timestep_number = 0;
   double time = 0;
-  
+
   do
-    { 
+    {
       std::cout << "Timestep " << timestep_number
                << ":  t=" << time
                << ", dt=" << time_step
-                << std::endl; 
+                << std::endl;
 
       std::cout << "   Assembling..." << std::endl;
-      assemble_system ();      
+      assemble_system ();
 
       std::cout << "   Solving..." << std::endl;
       solve ();
-      
+
       output_results ();
 
       time += time_step;
       ++timestep_number;
-   
-      old_solution = solution; 
+
+      old_solution = solution;
 
       std::cout << std::endl;
 
@@ -1227,9 +1279,9 @@ void BoussinesqFlowProblem<dim>::run ()
   while (time <= 5);
 }
 
-    
 
-int main () 
+
+int main ()
 {
   try
     {
@@ -1248,10 +1300,10 @@ int main ()
                 << "Aborting!" << std::endl
                 << "----------------------------------------------------"
                 << std::endl;
-      
+
       return 1;
     }
-  catch (...) 
+  catch (...)
     {
       std::cerr << std::endl << std::endl
                 << "----------------------------------------------------"

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.