From: Martin Kronbichler Date: Wed, 8 Mar 2017 14:58:24 +0000 (+0100) Subject: Allow selection of components in MatrixFreeOperators::Base X-Git-Tag: v8.5.0-rc1~47^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b8aa42b84990fc072ddc452e72ff66c3be1612c3;p=dealii.git Allow selection of components in MatrixFreeOperators::Base --- diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h index a52c638595..d5c5dd1724 100644 --- a/include/deal.II/matrix_free/operators.h +++ b/include/deal.II/matrix_free/operators.h @@ -116,7 +116,49 @@ namespace MatrixFreeOperators * Currently, the only supported vectors are LinearAlgebra::distributed::Vector * and LinearAlgebra::distributed::BlockVector. * - * @author Denis Davydov, Daniel Arndt, 2016, 2017 + *

Selective use of blocks in MatrixFree

+ * + * MatrixFree allows to use several DoFHandler/ConstraintMatrix combinations + * by passing an std::vector with pointers to the respective objects into + * the MatrixFree::reinit function. This class supports to select only some + * of the blocks in the underlying MatrixFree object by optional integer + * lists that specify the chosen blocks. + * + * One application of selection is in problems with a Newton-type iteration + * or problems with inhomogeneous boundary conditions. In such a case, one + * has to deal with two different sets of constraints: One set of + * constraints applies to the solution vector which might include hanging + * node constraints or periodicity constraints but no constraints on + * inhomogeneous Dirichlet boundaries. Before the nonlinear iteration, the + * boundary values are set to the expected value in the vector representing + * the initial guess. In each iteration of the Newton method, a linear + * system subject to zero Dirichlet boundary conditions is solved that is + * then added to the initial guess. This setup can be realized by using a + * vector of two pointers pointing to the same DoFHandler object and a + * vector of two pointers to the two ConstraintMatrix objects. If the first + * constraint matrix is the one including the zero Dirichlet constraints, + * one would give an std::vector(1, 0) to the initialize() + * function, i.e., a vector of length 1 that selects exactly the first + * constraint matrix with index 0. + * + * A second application of using selectively constructing a matrix-free + * operator would be the setting of the step-32 tutorial program: This + * problem has three blocks, one for the velocity, one for the + * pressure, and one for temperature. The time lag scheme used for temporal + * evolution splits the temperature equation away from the Stokes system in + * velocity and pressure. However, there are cross terms like the velocity + * that enters the temperature advection-diffusion equation or the + * temperature that enters the right hand side of the velocity. In order to + * be sure that MatrixFree uses the same integer indexing to the different + * blocks, one needs to put all three blocks into the same MatrixFree + * object. However, when solving a linear system the operators involved + * either address the first two in the Stokes solver, or the last one for + * the temperature solver. In the former case, a BlockVector of two + * components would be selected with a vector selecting the blocks {0, 1} in + * MatrixFree, whereas in the latter, a non-block vector selecting the block + * {2} would be used. + * + * @author Denis Davydov, Daniel Arndt, Martin Kronbichler, 2016, 2017 */ template > class Base : public Subscriptor @@ -150,23 +192,61 @@ namespace MatrixFreeOperators /** * Initialize operator on fine scale. + * + * The optional selection vector allows to choose only some components + * from the underlying MatrixFree object, e.g. just a single one. The + * entry @p selected_row_blocks[i] in the vector chooses the DoFHandler + * and ConstraintMatrix object that was given as the + * @p selected_row_blocks[i]-th argument to the MatrixFree::reinit() call. + * Different arguments for rows and columns also make it possible to + * select non-diagonal blocks or rectangular blocks. If the row vector is + * empty, all components are selected, otherwise its size must be smaller + * or equal to MatrixFree::n_components() and all indices need to be + * unique and within the range of 0 and MatrixFree::n_components(). If the + * column selection vector is empty, it is taken the same as the row + * selection, defining a diagonal block. */ - void initialize (std_cxx11::shared_ptr > data); + void initialize (std_cxx11::shared_ptr > data, + const std::vector &selected_row_blocks = std::vector(), + const std::vector &selected_column_blocks = std::vector()); /** * Initialize operator on a level @p level for a single FiniteElement. + * + * The optional selection vector allows to choose only some components + * from the underlying MatrixFree object, e.g. just a single one. The + * entry @p selected_row_blocks[i] in the vector chooses the DoFHandler + * and ConstraintMatrix object that was given as the + * @p selected_row_blocks[i]-th argument to the MatrixFree::reinit() call. + * Since a multigrid operator is always associated to inverting a matrix + * and thus represents a diagonal block, the same vector for rows and + * columns is used as opposed to the non-level initialization function. If + * empty, all components are selected. */ void initialize (std_cxx11::shared_ptr > data, const MGConstrainedDoFs &mg_constrained_dofs, - const unsigned int level); + const unsigned int level, + const std::vector &selected_row_blocks = std::vector()); /** - * Initialize operator on a level @p level for multiple FiniteElement objects. + * Initialize operator on a level @p level for multiple FiniteElement + * objects. + * + * The optional selection vector allows to choose only some components + * from the underlying MatrixFree object, e.g. just a single one. The + * entry @p selected_row_blocks[i] in the vector chooses the DoFHandler + * and ConstraintMatrix object that was given as the + * @p selected_row_blocks[i]-th argument to the MatrixFree::reinit() call. + * Since a multigrid operator is always associated to inverting a matrix + * and thus represents a diagonal block, the same vector for rows and + * columns is used as opposed to the non-level initialization function. If + * empty, all components are selected. */ void initialize(std_cxx11::shared_ptr > data_, const std::vector &mg_constrained_dofs, - const unsigned int level); + const unsigned int level, + const std::vector &selected_row_blocks = std::vector()); /** * Return the dimension of the codomain (or range) space. @@ -293,6 +373,18 @@ namespace MatrixFreeOperators */ std_cxx11::shared_ptr > inverse_diagonal_entries; + /** + * A vector which defines the selection of sub-components of MatrixFree + * for the rows of the matrix representation. + */ + std::vector selected_rows; + + /** + * A vector which defines the selection of sub-components of MatrixFree + * for the columns of the matrix representation. + */ + std::vector selected_columns; + private: /** @@ -327,7 +419,8 @@ namespace MatrixFreeOperators * order to ensure that the cell loops will be able to access the ghost * indices with the correct local indices. */ - void adjust_ghost_range_if_necessary(const VectorType &vec) const; + void adjust_ghost_range_if_necessary(const VectorType &vec, + const bool is_row) const; }; @@ -480,6 +573,10 @@ namespace MatrixFreeOperators /** * This class implements the operation of the action of a mass matrix. * + * Note that this class only supports the non-blocked vector variant of the + * Base operator because only a single FEEvaluation object is used in the + * apply function. + * * @author Daniel Arndt, 2016 */ template > @@ -531,6 +628,10 @@ namespace MatrixFreeOperators * namely $ L_{ij} = \int_\Omega c(\mathbf x) \mathbf \nabla N_i(\mathbf x) \cdot \mathbf \nabla N_j(\mathbf x)\,d \mathbf x$, * where $c(\mathbf x)$ is the scalar heterogeneity coefficient. * + * Note that this class only supports the non-blocked vector variant of the + * Base operator because only a single FEEvaluation object is used in the + * apply function. + * * @author Denis Davydov, 2016 */ template > @@ -795,8 +896,8 @@ namespace MatrixFreeOperators Assert(data.get() != NULL, ExcNotInitialized()); typename Base::size_type total_size = 0; - for (unsigned int i=0; in_components(); ++i) - total_size += data->get_vector_partitioner(i)->size(); + for (unsigned int i=0; iget_vector_partitioner(selected_rows[i])->size(); return total_size; } @@ -806,7 +907,12 @@ namespace MatrixFreeOperators typename Base::size_type Base::n () const { - return m(); + Assert(data.get() != NULL, + ExcNotInitialized()); + typename Base::size_type total_size = 0; + for (unsigned int i=0; iget_vector_partitioner(selected_columns[i])->size(); + return total_size; } @@ -841,13 +947,15 @@ namespace MatrixFreeOperators { Assert(data.get() != NULL, ExcNotInitialized()); + AssertDimension(n_blocks(vec), selected_rows.size()); for (unsigned int i = 0; i < n_blocks(vec); ++i) { - if (!subblock(vec,i).partitioners_are_compatible - (*data->get_dof_info(i).vector_partitioner)) - data->initialize_dof_vector(subblock(vec,i)); - Assert(subblock(vec,i).partitioners_are_globally_compatible - (*data->get_dof_info(i).vector_partitioner), + const unsigned int index = selected_rows[i]; + if (!subblock(vec,index).partitioners_are_compatible + (*data->get_dof_info(index).vector_partitioner)) + data->initialize_dof_vector(subblock(vec,index)); + Assert(subblock(vec,index).partitioners_are_globally_compatible + (*data->get_dof_info(index).vector_partitioner), ExcInternalError()); } collect_sizes(vec); @@ -858,13 +966,50 @@ namespace MatrixFreeOperators template void Base:: - initialize (std_cxx11::shared_ptr::value_type> > data_) + initialize (std_cxx11::shared_ptr::value_type> > data_, + const std::vector &given_row_selection, + const std::vector &given_column_selection) { data = data_; + + selected_rows.clear(); + selected_columns.clear(); + if (given_row_selection.empty()) + for (unsigned int i=0; in_components(); ++i) + selected_rows.push_back(i); + else + { + for (unsigned int i=0; in_components()); + for (unsigned int j=0; jn_components()); + for (unsigned int j=0; jn_components()); + edge_constrained_indices.resize(selected_rows.size()); edge_constrained_values.clear(); - edge_constrained_values.resize(data->n_components()); + edge_constrained_values.resize(selected_rows.size()); have_interface_matrices = false; } @@ -874,11 +1019,12 @@ namespace MatrixFreeOperators void Base:: initialize (std_cxx11::shared_ptr::value_type> > data_, - const MGConstrainedDoFs &mg_constrained_dofs, - const unsigned int level) + const MGConstrainedDoFs &mg_constrained_dofs, + const unsigned int level, + const std::vector &given_row_selection) { std::vector mg_constrained_dofs_vector(1, mg_constrained_dofs); - initialize(data_, mg_constrained_dofs_vector, level); + initialize(data_, mg_constrained_dofs_vector, level, given_row_selection); } @@ -888,19 +1034,41 @@ namespace MatrixFreeOperators Base:: initialize (std_cxx11::shared_ptr > data_, const std::vector &mg_constrained_dofs, - const unsigned int level) + const unsigned int level, + const std::vector &given_row_selection) { AssertThrow (level != numbers::invalid_unsigned_int, ExcMessage("level is not set")); - AssertDimension(mg_constrained_dofs.size(), data_->n_components()); + + selected_rows.clear(); + selected_columns.clear(); + if (given_row_selection.empty()) + for (unsigned int i=0; in_components(); ++i) + selected_rows.push_back(i); + else + { + for (unsigned int i=0; in_components()); + for (unsigned int j=0; jn_components()); + edge_constrained_indices.resize(selected_rows.size()); edge_constrained_values.clear(); - edge_constrained_values.resize(data_->n_components()); + edge_constrained_values.resize(selected_rows.size()); data = data_; - for (unsigned int j = 0; j < data->n_components(); ++j) + for (unsigned int j = 0; j < selected_rows.size(); ++j) { if (data_->n_macro_cells() > 0) { @@ -915,7 +1083,7 @@ namespace MatrixFreeOperators edge_constrained_indices[j].clear(); edge_constrained_indices[j].reserve(interface_indices.size()); edge_constrained_values[j].resize(interface_indices.size()); - const IndexSet &locally_owned = data->get_dof_handler(j).locally_owned_mg_dofs(level); + const IndexSet &locally_owned = data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(level); for (unsigned int i=0; i & - constrained_dofs = data->get_constrained_dofs(j); + constrained_dofs = data->get_constrained_dofs(selected_rows[j]); for (unsigned int i=0; i void - Base::adjust_ghost_range_if_necessary(const VectorType &src) const + Base::adjust_ghost_range_if_necessary(const VectorType &src, + const bool is_row) const { typedef typename Base::value_type Number; for (unsigned int i=0; i done if (subblock(src,i).get_partitioner().get() == - data->get_dof_info(i).vector_partitioner.get()) + data->get_dof_info(mf_component).vector_partitioner.get()) continue; // If not, assert that the local ranges are the same and reset to the // current partitioner Assert(subblock(src,i).get_partitioner()->local_size() == - data->get_dof_info(i).vector_partitioner->local_size(), + data->get_dof_info(mf_component).vector_partitioner->local_size(), ExcMessage("The vector passed to the vmult() function does not have " "the correct size for compatibility with MatrixFree.")); @@ -1002,7 +1172,7 @@ namespace MatrixFreeOperators subblock(src,i).begin()); const Vector copy_vec = view_src_in; subblock(const_cast(src),i). - reinit(data->get_dof_info(i).vector_partitioner); + reinit(data->get_dof_info(mf_component).vector_partitioner); VectorView view_src_out(subblock(src,i).local_size(), subblock(src,i).begin()); static_cast&>(view_src_out) = copy_vec; @@ -1019,9 +1189,10 @@ namespace MatrixFreeOperators { AssertDimension(dst.size(), src.size()); AssertDimension(n_blocks(dst), n_blocks(src)); + AssertDimension(n_blocks(dst), selected_rows.size()); typedef typename Base::value_type Number; - adjust_ghost_range_if_necessary(src); - adjust_ghost_range_if_necessary(dst); + adjust_ghost_range_if_necessary(src, false); + adjust_ghost_range_if_necessary(dst, true); // set zero Dirichlet values on the input vector (and remember the src and // dst values because we need to reset them at the end) @@ -1045,7 +1216,7 @@ namespace MatrixFreeOperators for (unsigned int j=0; j & - constrained_dofs = data->get_constrained_dofs(j); + constrained_dofs = data->get_constrained_dofs(selected_rows[j]); for (unsigned int i=0; i::value_type Number; AssertDimension(dst.size(), src.size()); - adjust_ghost_range_if_necessary(src); - adjust_ghost_range_if_necessary(dst); + adjust_ghost_range_if_necessary(src, false); + adjust_ghost_range_if_necessary(dst, true); dst = Number(0.); @@ -1129,8 +1300,8 @@ namespace MatrixFreeOperators { typedef typename Base::value_type Number; AssertDimension(dst.size(), src.size()); - adjust_ghost_range_if_necessary(src); - adjust_ghost_range_if_necessary(dst); + adjust_ghost_range_if_necessary(src, false); + adjust_ghost_range_if_necessary(dst, true); dst = Number(0.); @@ -1360,7 +1531,7 @@ namespace MatrixFreeOperators const std::pair &cell_range) const { typedef typename Base::value_type Number; - FEEvaluation phi(data); + FEEvaluation phi(data, this->selected_rows[0]); for (unsigned int cell=cell_range.first; cell &cell_range) const { typedef typename Base::value_type Number; - FEEvaluation phi (data); + FEEvaluation phi (data, this->selected_rows[0]); for (unsigned int cell=cell_range.first; cell &cell_range) const { typedef typename Base::value_type Number; - FEEvaluation phi (data); + FEEvaluation phi (data, this->selected_rows[0]); for (unsigned int cell=cell_range.first; cell