From 34c3dc3c4258c9472c88011af92aaa50009a9ef5 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 5 Jan 2010 13:52:18 +0000 Subject: [PATCH] Make it at least compile with a Schur complement based smoother. Doesn't run 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 | 151 +++++++++++++++++++++++++++- 1 file changed, 146 insertions(+), 5 deletions(-) diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index de3a193ba8..45f674a223 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -132,6 +132,7 @@ class StokesProblem std::vector > mg_dofs_per_component; std_cxx1x::shared_ptr::type> A_preconditioner; + std::vector::type> > mg_A_preconditioner; }; @@ -365,6 +366,7 @@ template void StokesProblem::setup_dofs () { A_preconditioner.reset (); + mg_A_preconditioner.resize (0); system_matrix.clear (); dof_handler.distribute_dofs (fe); @@ -668,6 +670,16 @@ void StokesProblem::assemble_multigrid () local_dof_indices, mg_interface_matrices[level]); } + + mg_A_preconditioner.resize (triangulation.n_levels()); + for (unsigned int level=0; level::type>(new typename InnerPreconditioner::type()); + mg_A_preconditioner[level] + ->initialize (mg_matrices[level].block(0,0), + typename InnerPreconditioner::type::AdditionalData()); + } } @@ -709,7 +721,127 @@ void StokesProblem::solve_block () << " block GMRES iterations"; } - template + +template +class SchurComplementSmoother +{ + public: + struct AdditionalData + { + const InnerPreconditioner *A_preconditioner; + }; + + void initialize (const BlockSparseMatrix &system_matrix, + const AdditionalData &data); + + void vmult (BlockVector &dst, + const BlockVector &src) const; + + void Tvmult (BlockVector &dst, + const BlockVector &src) const; + + void clear (); + + private: + SmartPointer > system_matrix; + SmartPointer A_preconditioner; +}; + + +template +void +SchurComplementSmoother:: +initialize (const BlockSparseMatrix &system_matrix, + const AdditionalData &data) +{ + this->system_matrix = &system_matrix; + this->A_preconditioner = data.A_preconditioner; +} + + + + +template +void +SchurComplementSmoother:: +vmult (BlockVector &dst, + const BlockVector &src) const +{ + const InverseMatrix,InnerPreconditioner> + A_inverse (system_matrix->block(0,0), *A_preconditioner); + Vector tmp (dst.block(0).size()); + + + { + Vector 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 + 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 preconditioner; + preconditioner.initialize (system_matrix->block(1,1), + SparseILU::AdditionalData()); + + InverseMatrix,SparseILU > + 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 +void +SchurComplementSmoother::clear () +{} + + + +template +void +SchurComplementSmoother:: +Tvmult (BlockVector &, + const BlockVector &) const +{ + Assert (false, ExcNotImplemented()); +} + + + +template void StokesProblem::solve_mg () { assemble_multigrid (); @@ -729,12 +861,21 @@ void StokesProblem::solve_mg () MGMatrix, BlockVector > mg_matrix(&mg_matrices); + typedef + SchurComplementSmoother::type> + Smoother; + + MGSmootherPrecondition, + Smoother, + BlockVector > + mg_smoother(mg_vector_memory); + +MGLevelObject +smoother_data (0, triangulation.n_levels()); - typedef PreconditionSSOR > REL; - MGSmootherRelaxation, REL, BlockVector > - mg_smoother(mg_vector_memory); +for (unsigned int level=0; level