From: kronbichler Date: Tue, 16 Sep 2008 14:06:41 +0000 (+0000) Subject: Make step-32 compliant with new TrilinosWrappers structure. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=60bd09c33beaf8205a7de6832d5d6155df73bb8d;p=dealii-svn.git Make step-32 compliant with new TrilinosWrappers structure. git-svn-id: https://svn.dealii.org/trunk@16835 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index c0c13e2f4a..b72f043873 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -189,8 +189,8 @@ namespace LinearSolvers const Preconditioner &preconditioner); - void vmult (TrilinosWrappers::Vector &dst, - const TrilinosWrappers::Vector &src) const; + void vmult (TrilinosWrappers::MPI::Vector &dst, + const TrilinosWrappers::MPI::Vector &src) const; private: const SmartPointer matrix; @@ -210,11 +210,11 @@ namespace LinearSolvers template void InverseMatrix::vmult ( - TrilinosWrappers::Vector &dst, - const TrilinosWrappers::Vector &src) const + TrilinosWrappers::MPI::Vector &dst, + const TrilinosWrappers::MPI::Vector &src) const { SolverControl solver_control (src.size(), 1e-6*src.l2_norm()); - SolverCG cg (solver_control); + SolverCG cg (solver_control); dst = 0; @@ -228,6 +228,8 @@ namespace LinearSolvers } } + + template class BlockSchurPreconditioner : public Subscriptor { @@ -237,14 +239,16 @@ namespace LinearSolvers const InverseMatrix &Mpinv, const PreconditionerA &Apreconditioner); - void vmult (TrilinosWrappers::BlockVector &dst, - const TrilinosWrappers::BlockVector &src) const; + void vmult (TrilinosWrappers::MPI::BlockVector &dst, + const TrilinosWrappers::MPI::BlockVector &src) const; private: const SmartPointer stokes_matrix; const SmartPointer > m_inverse; const PreconditionerA &a_preconditioner; + mutable TrilinosWrappers::MPI::Vector tmp; + }; @@ -257,16 +261,18 @@ namespace LinearSolvers : stokes_matrix (&S), m_inverse (&Mpinv), - a_preconditioner (Apreconditioner) + a_preconditioner (Apreconditioner), + tmp (stokes_matrix->block(1,1).matrix->RowMap()) {} + + template void BlockSchurPreconditioner::vmult ( - TrilinosWrappers::BlockVector &dst, - const TrilinosWrappers::BlockVector &src) const + TrilinosWrappers::MPI::BlockVector &dst, + const TrilinosWrappers::MPI::BlockVector &src) const { a_preconditioner.vmult (dst.block(0), src.block(0)); - TrilinosWrappers::Vector tmp (src.block(1), true); stokes_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1)); tmp *= -1; m_inverse->vmult (dst.block(1), tmp); @@ -332,8 +338,8 @@ class BoussinesqFlowProblem TrilinosWrappers::BlockSparseMatrix stokes_matrix; TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix; - TrilinosWrappers::BlockVector stokes_solution; - TrilinosWrappers::BlockVector stokes_rhs; + TrilinosWrappers::MPI::BlockVector stokes_solution; + TrilinosWrappers::MPI::BlockVector stokes_rhs; const unsigned int temperature_degree; @@ -346,10 +352,10 @@ class BoussinesqFlowProblem TrilinosWrappers::SparseMatrix temperature_stiffness_matrix; TrilinosWrappers::SparseMatrix temperature_matrix; - TrilinosWrappers::Vector temperature_solution; + TrilinosWrappers::MPI::Vector temperature_solution; TrilinosWrappers::Vector old_temperature_solution; TrilinosWrappers::Vector old_old_temperature_solution; - TrilinosWrappers::Vector temperature_rhs; + TrilinosWrappers::MPI::Vector temperature_rhs; double time_step; @@ -450,16 +456,13 @@ BoussinesqFlowProblem::get_extrapolated_temperature_range () const update_values); std::vector old_temperature_values(n_q_points); std::vector old_old_temperature_values(n_q_points); - - Vector localized_old_temperature_solution (old_temperature_solution); - Vector old_localized_old_temperature_solution (old_old_temperature_solution); double min_temperature = (1. + time_step/old_time_step) * old_temperature_solution.linfty_norm() + time_step/old_time_step * old_old_temperature_solution.linfty_norm(), - max_temperature = -min_temperature; + max_temperature = -min_temperature; typename DoFHandler::active_cell_iterator cell = temperature_dof_handler.begin_active(), @@ -468,8 +471,8 @@ BoussinesqFlowProblem::get_extrapolated_temperature_range () const if (cell->subdomain_id() == (unsigned int)trilinos_communicator.MyPID()) { fe_values.reinit (cell); - fe_values.get_function_values (localized_old_temperature_solution, old_temperature_values); - fe_values.get_function_values (old_localized_old_temperature_solution, old_old_temperature_values); + fe_values.get_function_values (old_temperature_solution, old_temperature_values); + fe_values.get_function_values (old_old_temperature_solution, old_old_temperature_values); for (unsigned int q=0; q::setup_dofs () BlockSparsityPattern stokes_sparsity_pattern (2,2); - //if (trilinos_communicator.MyPID() == 0) - { - BlockCompressedSetSparsityPattern csp (2,2); + BlockCompressedSetSparsityPattern csp (2,2); - csp.block(0,0).reinit (n_u, n_u); - 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.block(0,0).reinit (n_u, n_u); + 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 (); + csp.collect_sizes (); - Table<2,DoFTools::Coupling> coupling (dim+1, dim+1); + Table<2,DoFTools::Coupling> coupling (dim+1, dim+1); - for (unsigned int c=0; c::setup_dofs () BlockSparsityPattern stokes_preconditioner_sparsity_pattern (2,2); - //if (trilinos_communicator.MyPID() == 0) - { - BlockCompressedSetSparsityPattern csp (2,2); - - csp.block(0,0).reinit (n_u, n_u); - 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); - for (unsigned int c=0; c coupling (dim+1, dim+1); + for (unsigned int c=0; c::setup_dofs () SparsityPattern temperature_sparsity_pattern; - //if (trilinos_communicator.MyPID() == 0) - { - CompressedSetSparsityPattern csp (n_T, n_T); - DoFTools::make_sparsity_pattern (temperature_dof_handler, csp); - temperature_constraints.condense (csp); + CompressedSetSparsityPattern csp (n_T, n_T); + DoFTools::make_sparsity_pattern (temperature_dof_handler, csp); + temperature_constraints.condense (csp); - temperature_sparsity_pattern.copy_from (csp); - } + temperature_sparsity_pattern.copy_from (csp); temperature_matrix.reinit (temperature_partitioner, - temperature_sparsity_pattern); + temperature_sparsity_pattern); temperature_mass_matrix.reinit (temperature_partitioner, temperature_sparsity_pattern); temperature_stiffness_matrix.reinit (temperature_partitioner, - temperature_sparsity_pattern); + temperature_sparsity_pattern); } stokes_solution.reinit (stokes_partitioner); @@ -717,7 +711,6 @@ void BoussinesqFlowProblem::setup_dofs () temperature_solution.reinit (temperature_partitioner); old_temperature_solution.reinit (temperature_partitioner); old_old_temperature_solution.reinit (temperature_partitioner); - temperature_rhs.reinit (temperature_partitioner); } @@ -875,8 +868,6 @@ void BoussinesqFlowProblem::assemble_stokes_system () const FEValuesExtractors::Vector velocities (0); const FEValuesExtractors::Scalar pressure (dim); - Vector localized_old_temperature_solution (old_temperature_solution); - typename DoFHandler::active_cell_iterator cell = stokes_dof_handler.begin_active(), endc = stokes_dof_handler.end(); @@ -892,7 +883,7 @@ void BoussinesqFlowProblem::assemble_stokes_system () local_matrix = 0; local_rhs = 0; - temperature_fe_values.get_function_values (localized_old_temperature_solution, + temperature_fe_values.get_function_values (old_temperature_solution, old_temperature_values); for (unsigned int q=0; q::assemble_temperature_system () global_T_range = get_extrapolated_temperature_range(); const double global_Omega_diameter = GridTools::diameter (triangulation); - const Vector localized_old_temperature_solution (old_temperature_solution); - const Vector localized_old_old_temperature_solution (old_old_temperature_solution); - const BlockVector localized_stokes_solution (stokes_solution); + const TrilinosWrappers::BlockVector localized_stokes_solution (stokes_solution); typename DoFHandler::active_cell_iterator cell = temperature_dof_handler.begin_active(), @@ -1134,19 +1123,19 @@ void BoussinesqFlowProblem::assemble_temperature_system () temperature_fe_values.reinit (cell); stokes_fe_values.reinit (stokes_cell); - temperature_fe_values.get_function_values (localized_old_temperature_solution, + temperature_fe_values.get_function_values (old_temperature_solution, old_temperature_values); - temperature_fe_values.get_function_values (localized_old_old_temperature_solution, + temperature_fe_values.get_function_values (old_old_temperature_solution, old_old_temperature_values); - temperature_fe_values.get_function_gradients (localized_old_temperature_solution, + temperature_fe_values.get_function_gradients (old_temperature_solution, old_temperature_grads); - temperature_fe_values.get_function_gradients (localized_old_old_temperature_solution, + temperature_fe_values.get_function_gradients (old_old_temperature_solution, old_old_temperature_grads); - temperature_fe_values.get_function_hessians (localized_old_temperature_solution, + temperature_fe_values.get_function_hessians (old_temperature_solution, old_temperature_hessians); - temperature_fe_values.get_function_hessians (localized_old_old_temperature_solution, + temperature_fe_values.get_function_hessians (old_old_temperature_solution, old_old_temperature_hessians); temperature_right_hand_side.value_list (temperature_fe_values.get_quadrature_points(), @@ -1266,8 +1255,8 @@ void BoussinesqFlowProblem::solve () SolverControl solver_control (stokes_matrix.m(), 1e-6*stokes_rhs.l2_norm()); - SolverGMRES gmres(solver_control, - SolverGMRES::AdditionalData(100)); + SolverGMRES gmres(solver_control, + SolverGMRES::AdditionalData(100)); gmres.solve(stokes_matrix, stokes_solution, stokes_rhs, preconditioner); @@ -1276,7 +1265,7 @@ void BoussinesqFlowProblem::solve () << " GMRES iterations for Stokes subsystem." << std::endl; - BlockVector localized_stokes_solution (stokes_solution); + TrilinosWrappers::BlockVector localized_stokes_solution (stokes_solution); stokes_constraints.distribute (localized_stokes_solution); stokes_solution = localized_stokes_solution; } @@ -1295,7 +1284,7 @@ void BoussinesqFlowProblem::solve () SolverControl solver_control (temperature_matrix.m(), 1e-8*temperature_rhs.l2_norm()); - SolverCG cg (solver_control); + SolverCG cg (solver_control); TrilinosWrappers::PreconditionSSOR preconditioner (temperature_matrix, 1.2); @@ -1303,7 +1292,7 @@ void BoussinesqFlowProblem::solve () temperature_rhs, preconditioner); - Vector localized_temperature_solution (temperature_solution); + TrilinosWrappers::Vector localized_temperature_solution (temperature_solution); temperature_constraints.distribute (localized_temperature_solution); temperature_solution = localized_temperature_solution; @@ -1346,8 +1335,8 @@ void BoussinesqFlowProblem::output_results () const ExcInternalError()); Vector joint_solution (joint_dof_handler.n_dofs()); - BlockVector localized_stokes_solution (stokes_solution); - Vector localized_temperature_solution (temperature_solution); + TrilinosWrappers::BlockVector localized_stokes_solution (stokes_solution); + TrilinosWrappers::Vector localized_temperature_solution (temperature_solution); { std::vector local_joint_dof_indices (joint_fe.dofs_per_cell); @@ -1396,25 +1385,28 @@ void BoussinesqFlowProblem::output_results () const DataOut data_out; - data_out.attach_dof_handler (joint_dof_handler); + if (trilinos_communicator.MyPID() == 0) + { + data_out.attach_dof_handler (joint_dof_handler); - std::vector - data_component_interpretation - (dim+2, DataComponentInterpretation::component_is_scalar); - for (unsigned int i=0; i + data_component_interpretation + (dim+2, DataComponentInterpretation::component_is_scalar); + for (unsigned int i=0; i::type_dof_data, - data_component_interpretation); - data_out.build_patches (std::min(stokes_degree, temperature_degree)); + data_out.add_data_vector (joint_solution, joint_solution_names, + DataOut::type_dof_data, + data_component_interpretation); + data_out.build_patches (std::min(stokes_degree, temperature_degree)); - std::ostringstream filename; - filename << "solution-" << Utilities::int_to_string(timestep_number, 4) << ".vtk"; + std::ostringstream filename; + filename << "solution-" << Utilities::int_to_string(timestep_number, 4) << ".vtk"; - std::ofstream output (filename.str().c_str()); - data_out.write_vtk (output); + std::ofstream output (filename.str().c_str()); + data_out.write_vtk (output); + } } @@ -1425,7 +1417,7 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) { Vector estimated_error_per_cell (triangulation.n_active_cells()); - Vector localized_temperature_solution (temperature_solution); + TrilinosWrappers::Vector localized_temperature_solution (temperature_solution); KellyErrorEstimator::estimate (temperature_dof_handler, QGauss(temperature_degree+1), @@ -1442,11 +1434,11 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) cell != triangulation.end(); ++cell) cell->clear_refine_flag (); - std::vector > x_solution (2); + std::vector x_solution (2); x_solution[0] = temperature_solution; x_solution[1] = old_temperature_solution; - SolutionTransfer > soltrans(temperature_dof_handler); + SolutionTransfer soltrans(temperature_dof_handler); triangulation.prepare_coarsening_and_refinement(); soltrans.prepare_for_coarsening_and_refinement(x_solution); @@ -1454,7 +1446,7 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) triangulation.execute_coarsening_and_refinement (); setup_dofs (); - std::vector > tmp (2); + std::vector tmp (2); tmp[0] = temperature_solution; tmp[1] = temperature_solution; soltrans.interpolate(x_solution, tmp); 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 e50998fc2f..546378454d 100644 --- a/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h @@ -211,6 +211,21 @@ namespace TrilinosWrappers */ void compress (); + /** + * Returns the state of the + * matrix, i.e., whether + * compress() needs to be called + * after an operation requiring + * data exchange. Does only + * return non-true values when + * used in debug mode, + * since it is quite expensive to + * keep track of all operations + * that lead to the need for + * compress(). + */ + bool is_compressed () const; + /** * This function collects the * sizes of the sub-objects and @@ -221,10 +236,16 @@ namespace TrilinosWrappers * subobjects. You *must* call * this function each time after * you have changed the size of - * the sub-objects. Note that - * this is a collective - * operation, i.e., it needs to - * be called on all MPI processes. + * the sub-objects. Note that + * this is a collective + * operation, i.e., it needs to + * be called on all MPI + * processes. This command + * internally calls the method + * compress(), so you + * don't need to call that + * function in case you use + * collect_sizes(). */ void collect_sizes (); @@ -452,6 +473,24 @@ namespace TrilinosWrappers // ------------- inline and template functions ----------------- + inline + bool + BlockSparseMatrix::is_compressed () const + { + bool compressed = true; + for (unsigned int row=0; rowdebug mode, since + * it is quite expensive to + * keep track of all operations + * that lead to the need for + * compress(). + */ + bool is_compressed () const; + /** * Swap the contents of this * vector and the other vector @@ -319,6 +343,22 @@ namespace TrilinosWrappers + inline + bool + BlockVector::is_compressed () const + { + bool compressed = true; + for (unsigned int row=0; rowdebug mode, + * since it is quite expensive to + * keep track of all operations + * that lead to the need for + * compress(). + */ + bool is_compressed () const; + /** * This operator assigns a scalar * to a matrix. Since this does @@ -1384,6 +1399,13 @@ namespace TrilinosWrappers */ Epetra_CombineMode last_action; + /** + * A boolean variable to hold + * information on whether the + * vector is compressed or not. + */ + bool compressed; + public: /** * A sparse matrix object in @@ -1633,6 +1655,15 @@ namespace TrilinosWrappers (index < static_cast(end))); } + + + inline + bool + SparseMatrix::is_compressed () const + { + return compressed; + } + #endif // DOXYGEN } diff --git a/deal.II/lac/include/lac/trilinos_vector_base.h b/deal.II/lac/include/lac/trilinos_vector_base.h index fa8a0fdb0c..a24525734d 100644 --- a/deal.II/lac/include/lac/trilinos_vector_base.h +++ b/deal.II/lac/include/lac/trilinos_vector_base.h @@ -1,5 +1,5 @@ //--------------------------------------------------------------------------- -// $Id: trilinos_vector.h 16818 2008-09-12 14:30:24Z kronbichler $ +// $Id: trilinos_vector_base.h 16818 2008-09-12 14:30:24Z kronbichler $ // Version: $Name$ // // Copyright (C) 2008 by the deal.II authors @@ -313,9 +313,12 @@ namespace TrilinosWrappers * compress() has already been * called after an operation * requiring data exchange. Does - * only return non-false values + * only return non-true values * when used in debug - * mode. + * mode, since it is quite + * expensive to keep track of all + * operations that lead to a + * need for compress(). */ bool is_compressed () const; diff --git a/deal.II/lac/source/trilinos_block_vector.cc b/deal.II/lac/source/trilinos_block_vector.cc index 864f0b40de..a401e49e3d 100644 --- a/deal.II/lac/source/trilinos_block_vector.cc +++ b/deal.II/lac/source/trilinos_block_vector.cc @@ -62,6 +62,20 @@ namespace TrilinosWrappers + BlockVector & + BlockVector::operator = (const ::dealii::TrilinosWrappers::BlockVector &v) + { + Assert (n_blocks() == v.n_blocks(), + ExcDimensionMismatch(n_blocks(),v.n_blocks())); + + for (unsigned int i=0; in_blocks(); ++i) + this->components[i] = v.block(i); + + return *this; + } + + + BlockVector::~BlockVector () {} diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 5881615a48..38988ac32b 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -11,7 +11,6 @@ // //--------------------------------------------------------------------------- -#include #include #ifdef DEAL_II_USE_TRILINOS @@ -86,6 +85,7 @@ namespace TrilinosWrappers #endif col_map (row_map), last_action (Zero), + compressed (true), matrix (std::auto_ptr (new Epetra_FECrsMatrix(Copy, row_map, 0))) {} @@ -96,6 +96,7 @@ namespace TrilinosWrappers row_map (InputMap), col_map (row_map), last_action (Zero), + compressed (true), matrix (std::auto_ptr (new Epetra_FECrsMatrix(Copy, row_map, int(n_max_entries_per_row), false))) @@ -107,6 +108,7 @@ namespace TrilinosWrappers row_map (InputMap), col_map (row_map), last_action (Zero), + compressed (true), matrix (std::auto_ptr (new Epetra_FECrsMatrix(Copy, row_map, (int*)const_cast(&(n_entries_per_row[0])), @@ -120,6 +122,7 @@ namespace TrilinosWrappers row_map (InputRowMap), col_map (InputColMap), last_action (Zero), + compressed (true), matrix (std::auto_ptr (new Epetra_FECrsMatrix(Copy, row_map, col_map, int(n_max_entries_per_row), false))) @@ -132,6 +135,7 @@ namespace TrilinosWrappers row_map (InputRowMap), col_map (InputColMap), last_action (Zero), + compressed (true), matrix (std::auto_ptr (new Epetra_FECrsMatrix(Copy, row_map, col_map, (int*)const_cast(&(n_entries_per_row[0])), @@ -270,10 +274,12 @@ namespace TrilinosWrappers SparseMatrix::reinit (const SparseMatrix &sparse_matrix) { matrix.reset(); + row_map = sparse_matrix.row_map; col_map = sparse_matrix.col_map; matrix = std::auto_ptr(new Epetra_FECrsMatrix( *sparse_matrix.matrix)); + compressed = true; } @@ -354,6 +360,7 @@ namespace TrilinosWrappers &values[0], &row_indices[0]); } + compress(); } @@ -376,6 +383,7 @@ namespace TrilinosWrappers matrix = std::auto_ptr (new Epetra_FECrsMatrix(Copy, row_map, 0)); + compressed = true; } @@ -394,6 +402,8 @@ namespace TrilinosWrappers ierr = matrix->OptimizeStorage (); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + compressed = true; } @@ -437,6 +447,11 @@ namespace TrilinosWrappers if (last_action != Insert) last_action = Insert; +#ifdef DEBUG + if (in_local_range (i) == false) + compressed = false; +#endif + int trilinos_i = i; int trilinos_j = j; @@ -484,6 +499,11 @@ namespace TrilinosWrappers if (value == 0) return; +#ifdef DEBUG + if (in_local_range (i) == false) + compressed = false; +#endif + int trilinos_i = i; int trilinos_j = j; diff --git a/deal.II/lac/source/trilinos_vector_base.cc b/deal.II/lac/source/trilinos_vector_base.cc index 4d841242dd..dda97c4bca 100644 --- a/deal.II/lac/source/trilinos_vector_base.cc +++ b/deal.II/lac/source/trilinos_vector_base.cc @@ -1,5 +1,5 @@ //--------------------------------------------------------------------------- -// $Id: trilinos_vector.cc 16818 2008-09-12 14:30:24Z kronbichler $ +// $Id: trilinos_vector_base.cc 16818 2008-09-12 14:30:24Z kronbichler $ // Version: $Name$ // // Copyright (C) 2008 by the deal.II authors