From 0078676917fdaffe082eb4a27718877eb2d8e581 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Wed, 3 Sep 2008 11:08:14 +0000 Subject: [PATCH] Use Trilinos block structures on Stokes system now as well instead of the deal.II structures. Found an error in the implementation of TrilinosSparseMatrix. git-svn-id: https://svn.dealii.org/trunk@16721 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-31/step-31.cc | 165 +++++++++--------- .../lac/trilinos_block_sparse_matrix.h | 6 +- .../lac/include/lac/trilinos_block_vector.h | 8 +- deal.II/lac/source/trilinos_precondition.cc | 6 +- deal.II/lac/source/trilinos_sparse_matrix.cc | 8 +- 5 files changed, 96 insertions(+), 97 deletions(-) diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index c8ea6d0033..4cf90605a6 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -21,9 +21,9 @@ #include #include -#include +#include #include -#include +#include #include #include #include @@ -72,9 +72,6 @@ using namespace dealii; - // @sect3{Defining the AMG preconditioner} - - // @sect3{Equation data} // Again, the next stage in the program @@ -231,10 +228,6 @@ namespace EquationData } -namespace LinearSolvers -{ - - // @sect3{Linear solvers and preconditioners} @@ -248,6 +241,8 @@ namespace LinearSolvers // based preconditioners and solvers // have been introduced, with the // actual interface taken from step-22. +namespace LinearSolvers +{ // @sect4{The InverseMatrix class template} @@ -271,8 +266,9 @@ namespace LinearSolvers InverseMatrix (const Matrix &m, const Preconditioner &preconditioner); - void vmult (Vector &dst, - const Vector &src) const; + + void vmult (TrilinosWrappers::Vector &dst, + const TrilinosWrappers::Vector &src) const; private: const SmartPointer matrix; @@ -291,11 +287,12 @@ namespace LinearSolvers template - void InverseMatrix::vmult (Vector &dst, - const Vector &src) const + void InverseMatrix::vmult ( + TrilinosWrappers::Vector &dst, + const TrilinosWrappers::Vector &src) const { SolverControl solver_control (src.size(), 1e-6*src.l2_norm()); - SolverCG<> cg (solver_control); + SolverCG cg (solver_control); dst = 0; @@ -410,35 +407,36 @@ namespace LinearSolvers class BlockSchurPreconditioner : public Subscriptor { public: - BlockSchurPreconditioner (const BlockSparseMatrix &S, - const InverseMatrix,PreconditionerMp> &Mpinv, - const PreconditionerA &Apreconditioner); + BlockSchurPreconditioner ( + const TrilinosWrappers::BlockSparseMatrix &S, + const InverseMatrix &Mpinv, + const PreconditionerA &Apreconditioner); - void vmult (BlockVector &dst, - const BlockVector &src) const; + void vmult (TrilinosWrappers::BlockVector &dst, + const TrilinosWrappers::BlockVector &src) const; private: - const SmartPointer > stokes_matrix; - const SmartPointer, + const SmartPointer stokes_matrix; + const SmartPointer > m_inverse; const PreconditionerA &a_preconditioner; - mutable Vector tmp; + mutable TrilinosWrappers::Vector tmp; - }; +}; template BlockSchurPreconditioner:: - BlockSchurPreconditioner(const BlockSparseMatrix &S, - const InverseMatrix,PreconditionerMp> &Mpinv, - const PreconditionerA &Apreconditioner) + BlockSchurPreconditioner(const TrilinosWrappers::BlockSparseMatrix &S, + const InverseMatrix &Mpinv, + const PreconditionerA &Apreconditioner) : stokes_matrix (&S), m_inverse (&Mpinv), a_preconditioner (Apreconditioner), - tmp (S.block(1,1).m()) + tmp (stokes_matrix->block(1,1).row_map) {} @@ -463,8 +461,8 @@ namespace LinearSolvers // pressure vector. template void BlockSchurPreconditioner::vmult ( - BlockVector &dst, - const BlockVector &src) const + TrilinosWrappers::BlockVector &dst, + const TrilinosWrappers::BlockVector &src) const { a_preconditioner.vmult (dst.block(0), src.block(0)); stokes_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1)); @@ -526,30 +524,31 @@ class BoussinesqFlowProblem const double cell_diameter, const double old_time_step); - + Epetra_SerialComm trilinos_communicator; - - Triangulation triangulation; - const unsigned int stokes_degree; - FESystem stokes_fe; - DoFHandler stokes_dof_handler; - ConstraintMatrix stokes_constraints; - - BlockSparsityPattern stokes_sparsity_pattern; - BlockSparseMatrix stokes_matrix; - BlockSparsityPattern stokes_preconditioner_sparsity_pattern; - BlockSparseMatrix stokes_preconditioner_matrix; + Triangulation triangulation; - BlockVector stokes_solution; - BlockVector stokes_rhs; + const unsigned int stokes_degree; + FESystem stokes_fe; + DoFHandler stokes_dof_handler; + ConstraintMatrix stokes_constraints; + + std::vector stokes_partitioner; + BlockSparsityPattern stokes_sparsity_pattern; + TrilinosWrappers::BlockSparseMatrix stokes_matrix; + BlockSparsityPattern stokes_preconditioner_sparsity_pattern; + TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix; + + TrilinosWrappers::BlockVector stokes_solution; + TrilinosWrappers::BlockVector stokes_rhs; + + + const unsigned int temperature_degree; + FE_Q temperature_fe; + DoFHandler temperature_dof_handler; + ConstraintMatrix temperature_constraints; - - const unsigned int temperature_degree; - FE_Q temperature_fe; - DoFHandler temperature_dof_handler; - ConstraintMatrix temperature_constraints; - Epetra_Map temperature_partitioner; SparsityPattern temperature_sparsity_pattern; TrilinosWrappers::SparseMatrix temperature_mass_matrix; @@ -567,7 +566,7 @@ class BoussinesqFlowProblem unsigned int timestep_number; boost::shared_ptr Amg_preconditioner; - boost::shared_ptr > Mp_preconditioner; + boost::shared_ptr Mp_preconditioner; bool rebuild_stokes_matrix; bool rebuild_temperature_matrices; @@ -606,7 +605,7 @@ BoussinesqFlowProblem::BoussinesqFlowProblem () temperature_dof_handler (triangulation), temperature_partitioner (0, 0, trilinos_communicator), - + time_step (0), old_time_step (0), timestep_number (0), @@ -915,6 +914,13 @@ void BoussinesqFlowProblem::setup_dofs () // need to reassign the // system matrix structure to // the sparsity pattern. + stokes_partitioner.clear(); + { + Epetra_Map map_u(n_u, 0, trilinos_communicator); + stokes_partitioner.push_back (map_u); + Epetra_Map map_p(n_p, 0, trilinos_communicator); + stokes_partitioner.push_back (map_p); + } { stokes_matrix.clear (); @@ -924,7 +930,7 @@ void BoussinesqFlowProblem::setup_dofs () csp.block(0,1).reinit (n_u, n_p); csp.block(1,0).reinit (n_p, n_u); csp.block(1,1).reinit (n_p, n_p); - + csp.collect_sizes (); Table<2,DoFTools::Coupling> coupling (dim+1, dim+1); @@ -941,12 +947,13 @@ void BoussinesqFlowProblem::setup_dofs () coupling[c][d] = DoFTools::always; else coupling[c][d] = DoFTools::none; - + DoFTools::make_sparsity_pattern (stokes_dof_handler, coupling, csp); stokes_constraints.condense (csp); stokes_sparsity_pattern.copy_from (csp); - stokes_matrix.reinit (stokes_sparsity_pattern); + stokes_matrix.reinit (stokes_partitioner, stokes_sparsity_pattern); + stokes_matrix.collect_sizes(); } { @@ -975,7 +982,9 @@ void BoussinesqFlowProblem::setup_dofs () stokes_constraints.condense (csp); stokes_preconditioner_sparsity_pattern.copy_from (csp); - stokes_preconditioner_matrix.reinit (stokes_preconditioner_sparsity_pattern); + stokes_preconditioner_matrix.reinit (stokes_partitioner, + stokes_preconditioner_sparsity_pattern); + stokes_preconditioner_matrix.collect_sizes(); } temperature_partitioner = Epetra_Map (n_T, 0, trilinos_communicator); @@ -996,9 +1005,7 @@ void BoussinesqFlowProblem::setup_dofs () temperature_stiffness_matrix.reinit (temperature_partitioner, temperature_sparsity_pattern); } - - - + // As last action in this function, // we need to set the vectors // for the solution, the old @@ -1008,16 +1015,9 @@ void BoussinesqFlowProblem::setup_dofs () // three-block structure given // by velocity, pressure and // temperature. - stokes_solution.reinit (2); - stokes_solution.block(0).reinit (n_u); - stokes_solution.block(1).reinit (n_p); - stokes_solution.collect_sizes (); - - stokes_rhs.reinit (2); - stokes_rhs.block(0).reinit (n_u); - stokes_rhs.block(1).reinit (n_p); - stokes_rhs.collect_sizes (); - + stokes_solution.reinit (stokes_partitioner); + stokes_rhs.reinit (stokes_partitioner); + temperature_solution.reinit (temperature_partitioner); old_temperature_solution.reinit (temperature_partitioner); old_old_temperature_solution.reinit (temperature_partitioner); @@ -1080,6 +1080,7 @@ BoussinesqFlowProblem::assemble_stokes_preconditioner () local_dof_indices, stokes_preconditioner_matrix); } + stokes_preconditioner_matrix.compress(); } @@ -1132,11 +1133,10 @@ BoussinesqFlowProblem::build_stokes_preconditioner () // copied over to Trilinos. we need to // keep the (1,1) block, though - Mp_preconditioner = boost::shared_ptr > - (new SparseILU); - Mp_preconditioner->initialize (stokes_preconditioner_matrix.block(1,1), - SparseILU::AdditionalData()); - + Mp_preconditioner = boost::shared_ptr + (new TrilinosWrappers::PreconditionSSOR( + stokes_preconditioner_matrix.block(1,1),1.2)); + std::cout << std::endl; rebuild_stokes_preconditioner = false; @@ -1445,6 +1445,8 @@ void BoussinesqFlowProblem::assemble_stokes_system () local_dof_indices, stokes_rhs); } + stokes_matrix.compress(); + stokes_rhs.compress(); rebuild_stokes_matrix = false; @@ -1782,24 +1784,21 @@ void BoussinesqFlowProblem::solve () { // Set up inverse matrix for // pressure mass matrix - LinearSolvers::InverseMatrix,SparseILU > + LinearSolvers::InverseMatrix mp_inverse (stokes_preconditioner_matrix.block(1,1), *Mp_preconditioner); - // Set up block Schur preconditioner - /*BlockSchurPreconditioner::type, - SparseILU > - preconditioner (stokes_matrix, mp_inverse, *A_preconditioner);*/ LinearSolvers::BlockSchurPreconditioner > + TrilinosWrappers::PreconditionSSOR> preconditioner (stokes_matrix, mp_inverse, *Amg_preconditioner); // Set up GMRES solver and // solve. - SolverControl solver_control (stokes_matrix.m(), + SolverControl solver_control (20*stokes_matrix.m(), 1e-6*stokes_rhs.l2_norm()); - SolverGMRES > gmres(solver_control, - SolverGMRES >::AdditionalData(100)); + SolverGMRES gmres(solver_control, + SolverGMRES::AdditionalData(100)); stokes_solution = 0; gmres.solve(stokes_matrix, stokes_solution, stokes_rhs, preconditioner); @@ -2011,7 +2010,7 @@ void BoussinesqFlowProblem::run () const unsigned int initial_refinement = (dim == 2 ? 4 : 2); const unsigned int n_pre_refinement_steps = (dim == 2 ? 4 : 3); - + GridGenerator::hyper_cube (triangulation); triangulation.refine_global (initial_refinement); diff --git a/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h index 586a478be5..bddfb764e5 100644 --- a/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h @@ -250,7 +250,7 @@ namespace TrilinosWrappers * applicable if the matrix has * only one block column. */ - void vmult (BlockVector &dst, + void vmult (BlockVector &dst, const Vector &src) const; /** @@ -260,7 +260,7 @@ namespace TrilinosWrappers * applicable if the matrix has * only one block row. */ - void vmult (Vector &dst, + void vmult (Vector &dst, const BlockVector &src) const; /** @@ -312,7 +312,7 @@ namespace TrilinosWrappers * only one block. */ void Tvmult (Vector &dst, - const Vector &src) const; + const Vector &src) const; /** * Make the clear() function in the diff --git a/deal.II/lac/include/lac/trilinos_block_vector.h b/deal.II/lac/include/lac/trilinos_block_vector.h index d7b3c2d0e2..b944b98001 100644 --- a/deal.II/lac/include/lac/trilinos_block_vector.h +++ b/deal.II/lac/include/lac/trilinos_block_vector.h @@ -309,19 +309,19 @@ namespace TrilinosWrappers { unsigned int no_blocks = input_maps.size(); std::vector block_sizes (no_blocks); - + for (unsigned int i=0; iblock_indices.reinit (block_sizes); if (this->components.size() != this->n_blocks()) this->components.resize(this->n_blocks()); - + for (unsigned int i=0; in_blocks(); ++i) this->components[i].reinit(input_maps[i]); - + collect_sizes(); } diff --git a/deal.II/lac/source/trilinos_precondition.cc b/deal.II/lac/source/trilinos_precondition.cc index c980901f24..a079024d52 100755 --- a/deal.II/lac/source/trilinos_precondition.cc +++ b/deal.II/lac/source/trilinos_precondition.cc @@ -16,12 +16,10 @@ #include #include -#include - - - #ifdef DEAL_II_USE_TRILINOS +#include + DEAL_II_NAMESPACE_OPEN diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index b48139bc9f..1bea1a6616 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -766,11 +766,12 @@ namespace TrilinosWrappers void SparseMatrix::vmult_add (Vector &dst, - const Vector &src) const + const Vector &src) const { Assert (&src != &dst, ExcSourceEqualsDestination()); - Vector tmp = dst; + Vector tmp; + tmp = dst; vmult (dst, src); dst += tmp; } @@ -783,7 +784,8 @@ namespace TrilinosWrappers { Assert (&src != &dst, ExcSourceEqualsDestination()); - Vector tmp = dst; + Vector tmp; + tmp = dst; vmult (dst, src); dst += tmp; } -- 2.39.5