From: Martin Kronbichler Date: Mon, 1 Feb 2010 12:41:43 +0000 (+0000) Subject: Use new features introduced in branch_distributed_grids even in mainline. X-Git-Tag: v8.0.0~6565 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e563e5e2523fb5c46733c191ece9f6cae87f05bd;p=dealii.git Use new features introduced in branch_distributed_grids even in mainline. git-svn-id: https://svn.dealii.org/trunk@20481 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/index_set.h b/deal.II/base/include/base/index_set.h index 3f4dc6b0bb..2ca1086285 100644 --- a/deal.II/base/include/base/index_set.h +++ b/deal.II/base/include/base/index_set.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2009 by the deal.II authors +// Copyright (C) 2009, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -150,6 +150,24 @@ class IndexSet */ unsigned int index_within_set (const unsigned int n) const; + /** + * Each index set can be + * represented as the union of a + * number of contiguous intervals + * of indices, where if necessary + * intervals may only consist of + * individual elements to + * represent isolated members of + * the index set. + * + * This function returns the + * minimal number of such + * intervals that are needed to + * represent the index set under + * consideration. + */ + unsigned int n_intervals () const; + /** * Compress the internal * representation by merging @@ -206,6 +224,14 @@ class IndexSet IndexSet get_view (const unsigned int begin, const unsigned int end) const; + /** + * Outputs a text representation of this + * IndexSet to the given stream. Used for + * testing. + */ + template + void print(STREAM &out) const; + #ifdef DEAL_II_USE_TRILINOS /** * Given an MPI communicator, @@ -655,6 +681,26 @@ IndexSet::operator != (const IndexSet &is) const return ranges != is.ranges; } +template +inline +void +IndexSet::print (STREAM &out) const +{ + compress(); + out << "{"; + std::vector::const_iterator p; + for (p = ranges.begin(); p != ranges.end(); ++p) + { + if (p->end-p->begin==1) + out << p->begin; + else + out << "[" << p->begin << "," << p->end-1 << "]"; + + if (p !=--ranges.end()) + out << ", "; + } + out << "}" << std::endl; +} DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/base/source/index_set.cc b/deal.II/base/source/index_set.cc index 385571375b..ac5e5cf71f 100644 --- a/deal.II/base/source/index_set.cc +++ b/deal.II/base/source/index_set.cc @@ -154,6 +154,15 @@ IndexSet::operator & (const IndexSet &is) const +unsigned int +IndexSet::n_intervals () const +{ + compress (); + return ranges.size(); +} + + + IndexSet IndexSet::get_view (const unsigned int begin, const unsigned int end) const diff --git a/deal.II/lac/include/lac/constraint_matrix.h b/deal.II/lac/include/lac/constraint_matrix.h index 789519f156..4b7f379544 100644 --- a/deal.II/lac/include/lac/constraint_matrix.h +++ b/deal.II/lac/include/lac/constraint_matrix.h @@ -1546,6 +1546,17 @@ class ConstraintMatrix : public Subscriptor << "to another DoF with number " << arg1 << ", which however is constrained by this object. This is not" << " allowed."); + /** + * Exception + * + * @ingroup Exceptions + */ + DeclException1 (ExcRowNotStoredHere, + int, + << "The index set given to this constraint matrix indicates " + << "constraints for degree of freedom " << arg1 + << " should not be stored by this object, but a constraint " + << "is being added."); private: @@ -2075,7 +2086,8 @@ ConstraintMatrix::calculate_line_index (const unsigned int line) const if (!local_lines.size()) return line; - Assert(local_lines.is_element(line), ExcInternalError()); + Assert(local_lines.is_element(line), + ExcRowNotStoredHere(line)); return local_lines.index_within_set(line); } diff --git a/deal.II/lac/include/lac/constraint_matrix.templates.h b/deal.II/lac/include/lac/constraint_matrix.templates.h index e956c626b4..37aa2128da 100644 --- a/deal.II/lac/include/lac/constraint_matrix.templates.h +++ b/deal.II/lac/include/lac/constraint_matrix.templates.h @@ -1077,7 +1077,7 @@ namespace internals // specialized sort and insert functions. struct GlobalRowsFromLocal { - GlobalRowsFromLocal (const unsigned int n_local_rows) + GlobalRowsFromLocal (const unsigned int n_local_rows) : total_row_indices (n_local_rows), n_active_rows (n_local_rows), n_inhomogeneous_rows (0){}; @@ -1110,9 +1110,9 @@ namespace internals [index_in_constraint]).second; }; bool have_indirect_rows () const { return data_cache.element_size; }; void insert_constraint (const unsigned int constrained_local_dof) - { --n_active_rows; + { --n_active_rows; total_row_indices[n_active_rows].local_row = constrained_local_dof; } - unsigned int last_constrained_local_row () + unsigned int last_constrained_local_row () { Assert (total_row_indices.back().global_row == numbers::invalid_unsigned_int, ExcInternalError()); return total_row_indices.back().local_row; }; @@ -1447,7 +1447,7 @@ namespace internals { if (value != 0.) { - Assert (col_ptr != 0, + Assert (col_ptr != 0, typename SparseMatrix::ExcInvalidIndex (row, column)); if (are_on_diagonal) { @@ -1491,7 +1491,7 @@ namespace internals const unsigned int loc_row = global_rows.local_row(i); #ifdef DEBUG - const unsigned int * col_ptr = sparsity.n_nonzero_elements() == 0 ? 0 : + const unsigned int * col_ptr = sparsity.n_nonzero_elements() == 0 ? 0 : &sparsity_struct[row_start[row]]; number * val_ptr = sparsity.n_nonzero_elements() == 0 ? 0 : &sparse_matrix->global_entry (row_start[row]); diff --git a/deal.II/lac/include/lac/petsc_solver.h b/deal.II/lac/include/lac/petsc_solver.h index 0aac5ad434..545ce0d0ba 100644 --- a/deal.II/lac/include/lac/petsc_solver.h +++ b/deal.II/lac/include/lac/petsc_solver.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2004, 2005, 2006, 2007, 2009 by the deal.II authors +// Copyright (C) 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -42,6 +42,12 @@ namespace PETScWrappers * classes simply set the right flags to select one solver or another, or to * set certain parameters for individual solvers. * + * Note: Repeated calls to solve() on a solver object with a Preconditioner + * must be used with care. The preconditioner is initialized in the first call + * to solve() and subsequent calls reuse the solver and preconditioner + * object. This is done for performance reasons. The solver and preconditioner + * can be reset by calling reset(). + * * One of the gotchas of PETSc is that -- in particular in MPI mode -- it * often does not produce very helpful error messages. In order to save * other users some time in searching a hard to track down error, here is @@ -101,7 +107,10 @@ namespace PETScWrappers * classes and the object passed as a * preconditioner, one of the linear * solvers and preconditioners of PETSc - * is chosen. + * is chosen. Repeated calls to + * solve() do not reconstruct the + * preconditioner for performance + * reasons. See class Documentation. */ void solve (const MatrixBase &A, @@ -110,6 +119,14 @@ namespace PETScWrappers const PreconditionerBase &preconditioner); + /** + * Resets the contained preconditioner + * and solver object. See class + * description for more details. + */ + virtual void reset(); + + /** * Access to object that controls * convergence. diff --git a/deal.II/lac/include/lac/trilinos_precondition.h b/deal.II/lac/include/lac/trilinos_precondition.h index edfd3415cd..8f25060460 100644 --- a/deal.II/lac/include/lac/trilinos_precondition.h +++ b/deal.II/lac/include/lac/trilinos_precondition.h @@ -1308,9 +1308,11 @@ namespace TrilinosWrappers */ AdditionalData (const bool elliptic = true, const bool higher_order_elements = false, + const unsigned int n_cycles = 1, + const bool w_cyle = false, const double aggregation_threshold = 1e-4, const std::vector > &constant_modes = std::vector > (1), - const unsigned int smoother_sweeps = 3, + const unsigned int smoother_sweeps = 2, const unsigned int smoother_overlap = 0, const bool output_details = false); @@ -1335,6 +1337,20 @@ namespace TrilinosWrappers */ bool higher_order_elements; + /** + * Defines how many multigrid cycles + * should be performed by the + * preconditioner. + */ + unsigned int n_cycles; + + /** + * Defines whether a w-cycle should be + * used instead of the standard setting + * of a v-cycle. + */ + bool w_cycle; + /** * This threshold tells the AMG setup * how the coarsening should be diff --git a/deal.II/lac/include/lac/trilinos_vector.h b/deal.II/lac/include/lac/trilinos_vector.h index ec0a4e82e7..ae24942466 100644 --- a/deal.II/lac/include/lac/trilinos_vector.h +++ b/deal.II/lac/include/lac/trilinos_vector.h @@ -51,6 +51,7 @@ namespace TrilinosWrappers */ namespace MPI { + class BlockVector; /** * This class implements a wrapper to use the Trilinos distributed @@ -232,6 +233,9 @@ namespace TrilinosWrappers const bool fast = false, const bool allow_different_maps = false); + void reinit (const BlockVector &v, + const bool import_data = false); + /** * Set all components of the * vector to the given number @p diff --git a/deal.II/lac/source/constraint_matrix.cc b/deal.II/lac/source/constraint_matrix.cc index 4969878ad6..a4da438d13 100644 --- a/deal.II/lac/source/constraint_matrix.cc +++ b/deal.II/lac/source/constraint_matrix.cc @@ -264,10 +264,19 @@ void ConstraintMatrix::close () // this line (including // ones that we have // appended in this go - // around) + // around) and see whether + // they are further + // constrained. ignore + // elements that we don't + // store on the current + // processor unsigned int entry = 0; while (entry < line->entries.size()) - if (is_constrained (line->entries[entry].first)) + if (((local_lines.size() == 0) + || + (local_lines.is_element(line->entries[entry].first))) + && + is_constrained (line->entries[entry].first)) { // ok, this entry is // further @@ -504,20 +513,27 @@ void ConstraintMatrix::close () #ifdef DEBUG // if in debug mode: check that no // dof is constraint to another dof - // that is also constrained + // that is also + // constrained. exclude dofs from + // this check whose constraint + // lines would not be store on the + // local processor for (std::vector::const_iterator line=lines.begin(); line!=lines.end(); ++line) for (std::vector >::const_iterator entry=line->entries.begin(); entry!=line->entries.end(); ++entry) - { - // make sure that - // entry->first is not the - // index of a line itself - const bool is_circle = is_constrained(entry->first); - Assert (is_circle == false, - ExcDoFConstrainedToConstrainedDoF(line->line, entry->first)); - }; + if ((local_lines.size() == 0) + || + (local_lines.is_element(entry->first))) + { + // make sure that + // entry->first is not the + // index of a line itself + const bool is_circle = is_constrained(entry->first); + Assert (is_circle == false, + ExcDoFConstrainedToConstrainedDoF(line->line, entry->first)); + } #endif sorted = true; @@ -1919,6 +1935,7 @@ void ConstraintMatrix::condense (BlockCompressedSimpleSparsityPattern &sparsity) // on the local range of the vector. But // we need access to values that the // local nodes are constrained to. + template<> void ConstraintMatrix::distribute (TrilinosWrappers::MPI::Vector &vec) const @@ -1992,6 +2009,98 @@ ConstraintMatrix::distribute (TrilinosWrappers::MPI::Vector &vec) const } } + + +template<> +void +ConstraintMatrix::distribute (TrilinosWrappers::MPI::BlockVector &vec) const +{ + IndexSet my_indices (vec.size()); + for (unsigned int block=0; block::const_iterator constraint_iterator; + ConstraintLine index_comparison; + index_comparison.line = vec.block(block).local_range().first + +vec.get_block_indices().block_start(block); + const constraint_iterator begin_my_constraints = + std::lower_bound (lines.begin(),lines.end(),index_comparison); + + index_comparison.line = vec.block(block).local_range().second + +vec.get_block_indices().block_start(block); + + const constraint_iterator end_my_constraints + = std::lower_bound(lines.begin(),lines.end(),index_comparison); + + // Here we search all the indices that we + // need to have read-access to - the local + // nodes and all the nodes that the + // constraints indicate. No caching done + // yet. would need some more clever data + // structures for doing that. + const std::pair + local_range = vec.block(block).local_range(); + + my_indices.add_range (local_range.first, local_range.second); + + std::set individual_indices; + for (constraint_iterator it = begin_my_constraints; + it != end_my_constraints; ++it) + for (unsigned int i=0; ientries.size(); ++i) + if ((it->entries[i].first < local_range.first) + || + (it->entries[i].first >= local_range.second)) + individual_indices.insert (it->entries[i].first); + + my_indices.add_indices (individual_indices.begin(), + individual_indices.end()); + } + +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI + const Epetra_MpiComm *mpi_comm + = dynamic_cast(&vec.block(0).trilinos_vector().Comm()); + + Assert (mpi_comm != 0, ExcInternalError()); + + TrilinosWrappers::MPI::Vector vec_distribute + (my_indices.make_trilinos_map (mpi_comm->Comm(), true)); +#else + TrilinosWrappers::MPI::Vector vec_distribute + (my_indices.make_trilinos_map (MPI_COMM_WORLD, true))); +#endif + + // here we import the data + vec_distribute.reinit(vec,true); + + for (unsigned int block=0; block::const_iterator constraint_iterator; + ConstraintLine index_comparison; + index_comparison.line = vec.block(block).local_range().first + +vec.get_block_indices().block_start(block); + const constraint_iterator begin_my_constraints = + std::lower_bound (lines.begin(),lines.end(),index_comparison); + + index_comparison.line = vec.block(block).local_range().second + +vec.get_block_indices().block_start(block); + + const constraint_iterator end_my_constraints + = std::lower_bound(lines.begin(),lines.end(),index_comparison); + + for (constraint_iterator it = begin_my_constraints; + it != end_my_constraints; ++it) + { + // fill entry in line + // next_constraint.line by adding the + // different contributions + double new_value = it->inhomogeneity; + for (unsigned int i=0; ientries.size(); ++i) + new_value += (vec_distribute(it->entries[i].first) * + it->entries[i].second); + vec(it->line) = new_value; + } + } +} + #endif diff --git a/deal.II/lac/source/petsc_matrix_base.cc b/deal.II/lac/source/petsc_matrix_base.cc index 868e7261b0..f19b3944e8 100644 --- a/deal.II/lac/source/petsc_matrix_base.cc +++ b/deal.II/lac/source/petsc_matrix_base.cc @@ -620,7 +620,7 @@ namespace PETScWrappers std::cout << "*** PETSC-Matrix: num-allocs = " << info.mallocs << " ***" << std::endl; - return sizeof(*this) + info.memory; + return sizeof(*this) + static_cast(info.memory); } } diff --git a/deal.II/lac/source/petsc_solver.cc b/deal.II/lac/source/petsc_solver.cc index cc9bf46911..94ec618415 100644 --- a/deal.II/lac/source/petsc_solver.cc +++ b/deal.II/lac/source/petsc_solver.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2004, 2006, 2008, 2009 by the deal.II authors +// Copyright (C) 2004, 2006, 2008, 2009, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -198,8 +198,8 @@ namespace PETScWrappers # endif #endif - // destroy solver object - solver_data.reset (); + // do not destroy solver object +// solver_data.reset (); // in case of failure: throw // exception @@ -211,6 +211,14 @@ namespace PETScWrappers + void + SolverBase::reset() + { + solver_data.reset (); + } + + + SolverControl & SolverBase::control() const { diff --git a/deal.II/lac/source/trilinos_precondition.cc b/deal.II/lac/source/trilinos_precondition.cc index f0f337bbff..2be68d6be1 100644 --- a/deal.II/lac/source/trilinos_precondition.cc +++ b/deal.II/lac/source/trilinos_precondition.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2008 by the deal.II authors +// Copyright (C) 2008, 2009, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -488,6 +488,8 @@ namespace TrilinosWrappers PreconditionAMG::AdditionalData:: AdditionalData (const bool elliptic, const bool higher_order_elements, + const unsigned int n_cycles, + const bool w_cycle, const double aggregation_threshold, const std::vector > &constant_modes, const unsigned int smoother_sweeps, @@ -496,6 +498,8 @@ namespace TrilinosWrappers : elliptic (elliptic), higher_order_elements (higher_order_elements), + n_cycles (n_cycles), + w_cycle (w_cycle), aggregation_threshold (aggregation_threshold), constant_modes (constant_modes), smoother_sweeps (smoother_sweeps), @@ -513,8 +517,6 @@ namespace TrilinosWrappers multilevel_operator.release(); const unsigned int n_rows = matrix.m(); - const unsigned int constant_modes_dimension = - additional_data.constant_modes.size(); // Build the AMG preconditioner. Teuchos::ParameterList parameter_list; @@ -559,6 +561,13 @@ namespace TrilinosWrappers parameter_list.set("smoother: sweeps", static_cast(additional_data.smoother_sweeps)); + parameter_list.set("cycle applications", + static_cast(additional_data.n_cycles)); + if (additional_data.w_cycle == true) + parameter_list.set("prec type", "MGW"); + else + parameter_list.set("prec type", "MGV"); + parameter_list.set("smoother: Chebyshev alpha",10.); parameter_list.set("smoother: ifpack overlap", static_cast(additional_data.smoother_overlap)); @@ -572,23 +581,26 @@ namespace TrilinosWrappers const Epetra_Map & domain_map = matrix.domain_partitioner(); + const unsigned int constant_modes_dimension = + additional_data.constant_modes.size(); Epetra_MultiVector distributed_constant_modes (domain_map, constant_modes_dimension); if (constant_modes_dimension > 1) { - Assert (n_rows == additional_data.constant_modes[0].size(), - ExcDimensionMismatch(n_rows, - additional_data.constant_modes[0].size())); + const bool constant_modes_are_global = + additional_data.constant_modes[0].size() == n_rows; + const unsigned int n_relevant_rows = + constant_modes_are_global ? n_rows : additional_data.constant_modes[0].size(); + const unsigned int my_size = domain_map.NumMyElements(); + if (constant_modes_are_global == false) + Assert (n_relevant_rows == my_size, + ExcDimensionMismatch(n_relevant_rows, my_size)); Assert (n_rows == static_cast(distributed_constant_modes.GlobalLength()), ExcDimensionMismatch(n_rows, distributed_constant_modes.GlobalLength())); - const unsigned int my_size = domain_map.NumMyElements(); - Assert (my_size == static_cast(domain_map.MaxLID()+1), - ExcDimensionMismatch (my_size, domain_map.MaxLID()+1)); - // Reshape null space as a // contiguous vector of // doubles so that Trilinos @@ -596,7 +608,7 @@ namespace TrilinosWrappers for (unsigned int d=0; d(additional_data.constant_modes[d][global_row_id])); } diff --git a/deal.II/lac/source/trilinos_vector.cc b/deal.II/lac/source/trilinos_vector.cc index 00795b8c23..ed08c205ea 100644 --- a/deal.II/lac/source/trilinos_vector.cc +++ b/deal.II/lac/source/trilinos_vector.cc @@ -17,8 +17,10 @@ #ifdef DEAL_II_USE_TRILINOS # include +# include # include # include +# include DEAL_II_NAMESPACE_OPEN @@ -209,6 +211,73 @@ namespace TrilinosWrappers + void + Vector::reinit (const BlockVector &v, + const bool import_data) + { + // In case we do not allow to + // have different maps, this + // call means that we have to + // reset the vector. So clear + // the vector, initialize our + // map with the map in v, and + // generate the vector. + if (v.n_blocks() == 0) + return; + + // create a vector that holds all the elements + // contained in the block vector. need to + // manually create an Epetra_Map. + unsigned int n_elements = 0, added_elements = 0, block_offset = 0; + for (unsigned int block=0; block global_ids (n_elements, -1); + for (unsigned int block=0; block (new Epetra_FEVector (new_map)); + + Teuchos::RCP actual_vec = (import_data == true) ? + Teuchos::rcp (new Epetra_FEVector (new_map), true) : + Teuchos::rcp (&*vector, false); + + TrilinosScalar* entries = (*actual_vec)[0]; + block_offset = 0; + for (unsigned int block=0; block(actual_vec->GlobalLength()) + == v.size(), + ExcDimensionMismatch (actual_vec->GlobalLength(), + v.size())); + + Epetra_Import data_exchange (vector->Map(), actual_vec->Map()); + + const int ierr = vector->Import(*actual_vec, data_exchange, Insert); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + last_action = Insert; + } + + } + + + Vector & Vector::operator = (const Vector &v) { @@ -226,8 +295,9 @@ namespace TrilinosWrappers { Assert (vector->Map().SameAs(v.vector->Map()) == true, ExcMessage ("The Epetra maps in the assignment operator =" - " do not match, even though the local_range " - " seems to be the same. Check vector setup!")); + " do not match, even though the local_range" + " seems to be the same. Check vector setup or" + " use reinit()!")); const int ierr = vector->Update(1.0, *v.vector, 0.0); Assert (ierr == 0, ExcTrilinosError(ierr)); @@ -240,7 +310,21 @@ namespace TrilinosWrappers { Epetra_Import data_exchange (vector->Map(), v.vector->Map()); const int ierr = vector->Import(*v.vector, data_exchange, Insert); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + Assert (ierr == 0, ExcTrilinosError(ierr)); + } + else if (size() == v.size() && local_range().first>=v.local_range().first && + local_range().second<=v.local_range().second && + vector->Map().UniqueGIDs()) + { + for (unsigned int i=0; iMap().GID(i)), + ExcMessage ("The right hand side vector does not contain " + "all local elements present in the left hand " + "side vector in the assignment operation =, " + "which is not allowed. Use reinit().")); + (*vector)[0][i] = v(vector->Map().GID(i)); + } } else { diff --git a/deal.II/lac/source/trilinos_vector_base.cc b/deal.II/lac/source/trilinos_vector_base.cc index 5fba85de93..2cfee040e6 100644 --- a/deal.II/lac/source/trilinos_vector_base.cc +++ b/deal.II/lac/source/trilinos_vector_base.cc @@ -272,17 +272,9 @@ namespace TrilinosWrappers { // get a representation of the vector and // loop over all the elements - TrilinosScalar *start_ptr; - int leading_dimension; - int ierr = vector->ExtractView (&start_ptr, &leading_dimension); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); - - // TODO: This - // won't work in parallel like - // this. Find out a better way to - // this in that case. + TrilinosScalar *start_ptr = (*vector)[0]; const TrilinosScalar *ptr = start_ptr, - *eptr = start_ptr + size(); + *eptr = start_ptr + local_size(); bool flag = true; while (ptr != eptr) {