]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Make it at least compile with a Schur complement based smoother. Doesn't run
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 5 Jan 2010 13:52:18 +0000 (13:52 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 5 Jan 2010 13:52:18 +0000 (13:52 +0000)
right now, though.

git-svn-id: https://svn.dealii.org/trunk@20289 0785d39b-7218-0410-832d-ea1e28bc413d

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

index de3a193ba839c09f0cae820b158890d9f7728d73..45f674a223dd670505a187ff2ffd9adee181521c 100644 (file)
@@ -132,6 +132,7 @@ class StokesProblem
     std::vector<std::vector<unsigned int> >   mg_dofs_per_component;
                                  
     std_cxx1x::shared_ptr<typename InnerPreconditioner<dim>::type> A_preconditioner;
+    std::vector<std_cxx1x::shared_ptr<typename InnerPreconditioner<dim>::type> > mg_A_preconditioner;
 };
 
                        
@@ -365,6 +366,7 @@ template <int dim>
 void StokesProblem<dim>::setup_dofs ()
 {
   A_preconditioner.reset ();
+  mg_A_preconditioner.resize (0);
   system_matrix.clear ();
   
   dof_handler.distribute_dofs (fe);  
@@ -668,6 +670,16 @@ void StokesProblem<dim>::assemble_multigrid ()
                                   local_dof_indices,
                                   mg_interface_matrices[level]);
     }
+
+  mg_A_preconditioner.resize (triangulation.n_levels());
+  for (unsigned int level=0; level<triangulation.n_levels(); ++level)
+    {
+      mg_A_preconditioner[level]
+       = std_cxx1x::shared_ptr<typename InnerPreconditioner<dim>::type>(new typename InnerPreconditioner<dim>::type());
+      mg_A_preconditioner[level]
+       ->initialize (mg_matrices[level].block(0,0),
+                     typename InnerPreconditioner<dim>::type::AdditionalData());
+    }
 }
 
                                
@@ -709,7 +721,127 @@ void StokesProblem<dim>::solve_block ()
     << " block GMRES iterations";
 }
 
-  template <int dim>
+
+template <typename InnerPreconditioner>
+class SchurComplementSmoother
+{
+  public:
+    struct AdditionalData 
+    {
+       const InnerPreconditioner *A_preconditioner;
+    };
+
+    void initialize (const BlockSparseMatrix<double> &system_matrix,
+                    const AdditionalData            &data);
+    
+    void vmult (BlockVector<double> &dst,
+               const BlockVector<double> &src) const;
+
+    void Tvmult (BlockVector<double> &dst,
+                const BlockVector<double> &src) const;
+
+    void clear ();
+    
+  private:
+    SmartPointer<const BlockSparseMatrix<double> > system_matrix;
+    SmartPointer<const InnerPreconditioner>        A_preconditioner;
+};
+
+
+template <typename InnerPreconditioner>
+void
+SchurComplementSmoother<InnerPreconditioner>::
+initialize (const BlockSparseMatrix<double> &system_matrix,
+           const AdditionalData            &data)
+{
+  this->system_matrix    = &system_matrix;
+  this->A_preconditioner = data.A_preconditioner;
+}
+
+
+
+
+template <typename InnerPreconditioner>
+void
+SchurComplementSmoother<InnerPreconditioner>::
+vmult (BlockVector<double> &dst,
+       const BlockVector<double> &src) const
+{
+  const InverseMatrix<SparseMatrix<double>,InnerPreconditioner>
+    A_inverse (system_matrix->block(0,0), *A_preconditioner);
+  Vector<double> tmp (dst.block(0).size());
+  
+                                 
+  {
+    Vector<double> schur_rhs (dst.block(1).size());
+    A_inverse.vmult (tmp, src.block(0));
+    system_matrix->block(1,0).vmult (schur_rhs, tmp);
+    schur_rhs -= src.block(1);
+  
+    SchurComplement<InnerPreconditioner>
+      schur_complement (*system_matrix, A_inverse);
+    
+                                    // The usual control structures for
+                                    // the solver call are created...
+    SolverControl solver_control (dst.block(1).size(),
+                                 1e-6*schur_rhs.l2_norm());
+    SolverCG<>    cg (solver_control);
+    
+                                   
+                                  
+    SparseILU<double> preconditioner;
+    preconditioner.initialize (system_matrix->block(1,1), 
+                              SparseILU<double>::AdditionalData());
+  
+    InverseMatrix<SparseMatrix<double>,SparseILU<double> >
+      m_inverse (system_matrix->block(1,1), preconditioner);
+    
+                                   
+    cg.solve (schur_complement, dst.block(1), schur_rhs,
+             m_inverse);
+  
+// no constraints to be taken care of here
+    
+    std::cout << "  "
+             << solver_control.last_step()
+             << " CG Schur complement iterations in smoother"
+             << std::flush
+             << std::endl;    
+  }
+    
+                                
+  {
+    system_matrix->block(0,1).vmult (tmp, dst.block(1));
+    tmp *= -1;
+    tmp += src.block(0);
+  
+    A_inverse.vmult (dst.block(0), tmp);
+
+// no constraints here either
+  }
+}
+
+
+
+template <typename InnerPreconditioner>
+void
+SchurComplementSmoother<InnerPreconditioner>::clear ()
+{}
+
+
+
+template <typename InnerPreconditioner>
+void
+SchurComplementSmoother<InnerPreconditioner>::
+Tvmult (BlockVector<double> &,
+       const BlockVector<double> &) const
+{
+  Assert (false, ExcNotImplemented());
+}
+
+
+
+template <int dim>
 void StokesProblem<dim>::solve_mg () 
 {
   assemble_multigrid ();
@@ -729,12 +861,21 @@ void StokesProblem<dim>::solve_mg ()
   MGMatrix<BlockSparseMatrix<double>, BlockVector<double> >
     mg_matrix(&mg_matrices);
 
+  typedef
+    SchurComplementSmoother<typename InnerPreconditioner<dim>::type>
+    Smoother;
+    
+  MGSmootherPrecondition<BlockSparseMatrix<double>,
+    Smoother,
+    BlockVector<double> >
+  mg_smoother(mg_vector_memory);
+
+MGLevelObject<typename Smoother::AdditionalData>
+smoother_data (0, triangulation.n_levels());
 
-  typedef PreconditionSSOR<BlockSparseMatrix<double> > REL;
-  MGSmootherRelaxation<BlockSparseMatrix<double>, REL, BlockVector<double> >
-    mg_smoother(mg_vector_memory);
+for (unsigned int level=0; level<triangulation.n_levels(); ++level)
+  smoother_data[level].A_preconditioner = mg_A_preconditioner[level].get();
 
-  REL::AdditionalData smoother_data;
   mg_smoother.initialize(mg_matrices, smoother_data);
   mg_smoother.set_steps(2);
 

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.