From e74890e86de9570a19cc3abdd360c9669b5f1eee Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 17 Oct 2007 15:19:06 +0000 Subject: [PATCH] Use umfpack as a (very good) preconditioner. Improve screen output. git-svn-id: https://svn.dealii.org/trunk@15337 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-22/step-22.cc | 83 ++++++++++++++++++----------- 1 file changed, 53 insertions(+), 30 deletions(-) diff --git a/deal.II/examples/step-22/step-22.cc b/deal.II/examples/step-22/step-22.cc index 749dd6d007..13da6db65a 100644 --- a/deal.II/examples/step-22/step-22.cc +++ b/deal.II/examples/step-22/step-22.cc @@ -23,6 +23,7 @@ #include #include #include +#include #include #include @@ -47,8 +48,6 @@ #include #include -#include - using namespace dealii; @@ -85,6 +84,10 @@ class BoussinesqFlowProblem BlockVector solution; BlockVector old_solution; BlockVector system_rhs; + + boost::shared_ptr A_preconditioner; + + bool rebuild_preconditioner; }; @@ -264,40 +267,41 @@ extract_grad_T (const FEValuesBase &fe_values, -template +template class InverseMatrix : public Subscriptor { public: - InverseMatrix (const Matrix &m); + InverseMatrix (const Matrix &m, + const Preconditioner &preconditioner); void vmult (Vector &dst, const Vector &src) const; private: const SmartPointer matrix; + const Preconditioner &preconditioner; mutable GrowingVectorMemory<> vector_memory; }; -template -InverseMatrix::InverseMatrix (const Matrix &m) +template +InverseMatrix::InverseMatrix (const Matrix &m, + const Preconditioner &preconditioner) : - matrix (&m) + matrix (&m), + preconditioner (preconditioner) {} -template -void InverseMatrix::vmult (Vector &dst, - const Vector &src) const +template +void InverseMatrix::vmult (Vector &dst, + const Vector &src) const { SolverControl solver_control (src.size(), 1e-6*src.l2_norm()); SolverCG<> cg (solver_control, vector_memory); - PreconditionJacobi<> preconditioner; - preconditioner.initialize (*matrix); - dst = 0; try @@ -312,27 +316,29 @@ void InverseMatrix::vmult (Vector &dst, +template class SchurComplement : public Subscriptor { public: SchurComplement (const BlockSparseMatrix &A, - const InverseMatrix > &Minv); + const InverseMatrix,Preconditioner> &Minv); void vmult (Vector &dst, const Vector &src) const; private: const SmartPointer > system_matrix; - const SmartPointer > > m_inverse; + const SmartPointer,Preconditioner> > m_inverse; mutable Vector tmp1, tmp2; }; -SchurComplement:: +template +SchurComplement:: SchurComplement (const BlockSparseMatrix &A, - const InverseMatrix > &Minv) + const InverseMatrix,Preconditioner> &Minv) : system_matrix (&A), m_inverse (&Minv), @@ -341,8 +347,9 @@ SchurComplement (const BlockSparseMatrix &A, {} -void SchurComplement::vmult (Vector &dst, - const Vector &src) const +template +void SchurComplement::vmult (Vector &dst, + const Vector &src) const { system_matrix->block(0,1).vmult (tmp1, src); m_inverse->vmult (tmp2, tmp1); @@ -359,8 +366,9 @@ BoussinesqFlowProblem::BoussinesqFlowProblem (const unsigned int degree) FE_Q(degree), 1, FE_Q(degree), 1), dof_handler (triangulation), - n_refinement_steps (4), - time_step (0) + n_refinement_steps (6), + time_step (0), + rebuild_preconditioner (true) {} @@ -579,6 +587,19 @@ void BoussinesqFlowProblem::assemble_system () solution, system_rhs); } + + if (rebuild_preconditioner == true) + { + std::cout << " Rebuilding preconditioner..." << std::flush; + + A_preconditioner + = boost::shared_ptr(new SparseDirectUMFPACK()); + A_preconditioner->initialize (system_matrix.block(0,0)); + + std::cout << std::endl; + + rebuild_preconditioner = false; + } } @@ -715,8 +736,8 @@ void BoussinesqFlowProblem::assemble_rhs_T () template void BoussinesqFlowProblem::solve () { - const InverseMatrix > - A_inverse (system_matrix.block(0,0)); + const InverseMatrix,SparseDirectUMFPACK> + A_inverse (system_matrix.block(0,0), *A_preconditioner); Vector tmp (solution.block(0).size()); Vector schur_rhs (solution.block(1).size()); Vector tmp2 (solution.block(2).size()); @@ -728,7 +749,7 @@ void BoussinesqFlowProblem::solve () schur_rhs -= system_rhs.block(1); - SchurComplement + SchurComplement schur_complement (system_matrix, A_inverse); SolverControl solver_control (system_matrix.block(0,0).m(), @@ -889,22 +910,24 @@ void BoussinesqFlowProblem::run () do { std::cout << "Timestep " << timestep_number + << ": t=" << time + << ", dt=" << time_step << std::endl; - assemble_system (); + std::cout << " Assembling..." << std::endl; + assemble_system (); + std::cout << " Solving..." << std::endl; solve (); output_results (); time += time_step; ++timestep_number; - std::cout << " Now at t=" << time - << ", dt=" << time_step << '.' - << std::endl - << std::endl; + + std::cout << std::endl; } - while (time <= 50); + while (time <= 10); } -- 2.39.5