From 1777848615ffc84e97c5d80529224fefe2d89032 Mon Sep 17 00:00:00 2001 From: tcclevenger Date: Wed, 21 Aug 2019 15:46:02 -0400 Subject: [PATCH] add custom constraints --- doc/news/changes/minor/20190913Clevenger | 5 + .../deal.II/multigrid/mg_constrained_dofs.h | 66 ++++++++++++- source/lac/affine_constraints.inst.in | 4 + source/multigrid/mg_transfer_matrix_free.cc | 26 ++++- source/multigrid/mg_transfer_prebuilt.cc | 8 ++ source/multigrid/multigrid.cc | 16 +++ tests/multigrid/transfer_matrix_free_13.cc | 98 +++++++++++++++++++ .../multigrid/transfer_matrix_free_13.output | 7 ++ 8 files changed, 228 insertions(+), 2 deletions(-) create mode 100644 doc/news/changes/minor/20190913Clevenger create mode 100644 tests/multigrid/transfer_matrix_free_13.cc create mode 100644 tests/multigrid/transfer_matrix_free_13.output diff --git a/doc/news/changes/minor/20190913Clevenger b/doc/news/changes/minor/20190913Clevenger new file mode 100644 index 0000000000..74c2b4116e --- /dev/null +++ b/doc/news/changes/minor/20190913Clevenger @@ -0,0 +1,5 @@ +New: MGConstrainedDoFs now allows the user to define constraints +on each multigrid level to be applied during prolongation in a +matrix-free transfer. +
+(Conrad Clevenger, 2019/09/13) diff --git a/include/deal.II/multigrid/mg_constrained_dofs.h b/include/deal.II/multigrid/mg_constrained_dofs.h index eaa5ea9f03..614ab15be4 100644 --- a/include/deal.II/multigrid/mg_constrained_dofs.h +++ b/include/deal.II/multigrid/mg_constrained_dofs.h @@ -103,6 +103,24 @@ public: const std::set &boundary_ids, const ComponentMask & component_mask = ComponentMask()); + /** + * Add user defined constraints to be used on level @p level. + * + * The user can call this function multiple times and any new, + * conflicting constraints will overwrite the previous constraints + * for that DoF. + * + * Before the transfer, the user defined constraints will be distributed + * to the source vector, and then any DoF index set using + * make_zero_boundary_constraints() will be overwritten with + * value zero. + * + * @note This is currently only implemented for MGTransferMatrixFree. + */ + void + add_user_constraints(const unsigned int level, + const AffineConstraints &constraints_on_level); + /** * Fill the internal data structures with information * about no normal flux boundary dofs. @@ -196,6 +214,16 @@ public: const AffineConstraints & get_level_constraint_matrix(const unsigned int level) const; + /** + * Return the user defined constraint matrix for a given level. These + * constraints are set using the function add_user_constraints() and + * should not contain constraints for DoF indices set in + * make_zero_boundary_constraints() as they will be overwritten during + * the transfer. + */ + const AffineConstraints & + get_user_constraint_matrix(const unsigned int level) const; + private: /** * The indices of boundary dofs for each level. @@ -213,6 +241,11 @@ private: * periodic boundary conditions for each level . */ std::vector> level_constraints; + + /** + * Constraint matrices defined by user. + */ + std::vector> user_constraints; }; @@ -223,12 +256,14 @@ MGConstrainedDoFs::initialize(const DoFHandler &dof) boundary_indices.clear(); refinement_edge_indices.clear(); level_constraints.clear(); + user_constraints.clear(); const unsigned int nlevels = dof.get_triangulation().n_global_levels(); // At this point level_constraint and refinement_edge_indices are empty. - level_constraints.resize(nlevels); refinement_edge_indices.resize(nlevels); + level_constraints.resize(nlevels); + user_constraints.resize(nlevels); for (unsigned int l = 0; l < nlevels; ++l) { IndexSet relevant_dofs; @@ -399,11 +434,31 @@ MGConstrainedDoFs::make_no_normal_flux_constraints( } +inline void +MGConstrainedDoFs::add_user_constraints( + const unsigned int level, + const AffineConstraints &constraints_on_level) +{ + AssertIndexRange(level, user_constraints.size()); + + // Get the relevant DoFs from level_constraints if + // the user constraint matrix has not been initialized + if (user_constraints[level].get_local_lines().size() == 0) + user_constraints[level].reinit(level_constraints[level].get_local_lines()); + + user_constraints[level].merge( + constraints_on_level, + AffineConstraints::MergeConflictBehavior::right_object_wins); + user_constraints[level].close(); +} + + inline void MGConstrainedDoFs::clear() { boundary_indices.clear(); refinement_edge_indices.clear(); + user_constraints.clear(); } @@ -487,6 +542,15 @@ MGConstrainedDoFs::get_level_constraint_matrix(const unsigned int level) const +inline const AffineConstraints & +MGConstrainedDoFs::get_user_constraint_matrix(const unsigned int level) const +{ + AssertIndexRange(level, user_constraints.size()); + return user_constraints[level]; +} + + + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/source/lac/affine_constraints.inst.in b/source/lac/affine_constraints.inst.in index 4523cad087..e3877ee611 100644 --- a/source/lac/affine_constraints.inst.in +++ b/source/lac/affine_constraints.inst.in @@ -344,4 +344,8 @@ for (T : DEAL_II_VEC_TEMPLATES) T &) const; template void dealii::AffineConstraints::distribute>( T &) const; + + template void dealii::AffineConstraints::distribute< + LinearAlgebra::distributed::T>( + LinearAlgebra::distributed::T &) const; } diff --git a/source/multigrid/mg_transfer_matrix_free.cc b/source/multigrid/mg_transfer_matrix_free.cc index 49557a66dd..769c532e8c 100644 --- a/source/multigrid/mg_transfer_matrix_free.cc +++ b/source/multigrid/mg_transfer_matrix_free.cc @@ -355,6 +355,30 @@ MGTransferMatrixFree::do_prolongate_add( Utilities::fixed_power(n_child_dofs_1d); constexpr unsigned int three_to_dim = Utilities::pow(3, dim); + // If we have user defined MG constraints, we must create + // a non-const, ghosted version of the source vector to distribute + // constraints. + const LinearAlgebra::distributed::Vector *to_use = &src; + LinearAlgebra::distributed::Vector new_src; + if (this->mg_constrained_dofs != nullptr && + this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1) + .get_local_lines() + .size() > 0) + { + LinearAlgebra::distributed::Vector copy_src(src); + + // Distribute any user defined constraints + this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1) + .distribute(copy_src); + + // Re-initialize new ghosted vector with correct constraints + new_src.reinit(copy_src); + new_src = copy_src; + new_src.update_ghost_values(); + + to_use = &new_src; + } + for (unsigned int cell = 0; cell < n_owned_level_cells[to_level - 1]; cell += vec_size) { @@ -382,7 +406,7 @@ MGTransferMatrixFree::do_prolongate_add( for (unsigned int k = 0; k < (dim > 2 ? degree_size : 1); ++k) for (unsigned int j = 0; j < (dim > 1 ? degree_size : 1); ++j) for (unsigned int i = 0; i < degree_size; ++i, ++m) - evaluation_data[m][v] = src.local_element( + evaluation_data[m][v] = to_use->local_element( indices[c * n_scalar_cell_dofs + k * n_child_dofs_1d * n_child_dofs_1d + j * n_child_dofs_1d + i]); diff --git a/source/multigrid/mg_transfer_prebuilt.cc b/source/multigrid/mg_transfer_prebuilt.cc index a65380efff..0eceb6b356 100644 --- a/source/multigrid/mg_transfer_prebuilt.cc +++ b/source/multigrid/mg_transfer_prebuilt.cc @@ -106,6 +106,14 @@ MGTransferPrebuilt::prolongate(const unsigned int to_level, Assert((to_level >= 1) && (to_level <= prolongation_matrices.size()), ExcIndexRange(to_level, 1, prolongation_matrices.size() + 1)); +#ifdef DEBUG + if (this->mg_constrained_dofs != nullptr) + Assert(this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1) + .get_local_lines() + .size() == 0, + ExcNotImplemented()); +#endif + prolongation_matrices[to_level - 1]->vmult(dst, src); } diff --git a/source/multigrid/multigrid.cc b/source/multigrid/multigrid.cc index ff37ab9eb5..310fc0da13 100644 --- a/source/multigrid/multigrid.cc +++ b/source/multigrid/multigrid.cc @@ -94,6 +94,14 @@ MGTransferBlock::prolongate(const unsigned int to_level, Assert(dst.n_blocks() == this->n_mg_blocks, ExcDimensionMismatch(dst.n_blocks(), this->n_mg_blocks)); +#ifdef DEBUG + if (this->mg_constrained_dofs != nullptr) + Assert(this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1) + .get_local_lines() + .size() == 0, + ExcNotImplemented()); +#endif + // Multiplicate with prolongation // matrix, but only those blocks // selected. @@ -280,6 +288,14 @@ MGTransferBlockSelect::prolongate(const unsigned int to_level, Assert((to_level >= 1) && (to_level <= prolongation_matrices.size()), ExcIndexRange(to_level, 1, prolongation_matrices.size() + 1)); +#ifdef DEBUG + if (this->mg_constrained_dofs != nullptr) + Assert(this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1) + .get_local_lines() + .size() == 0, + ExcNotImplemented()); +#endif + prolongation_matrices[to_level - 1] ->block(selected_block, selected_block) .vmult(dst, src); diff --git a/tests/multigrid/transfer_matrix_free_13.cc b/tests/multigrid/transfer_matrix_free_13.cc new file mode 100644 index 0000000000..dd1e39bc11 --- /dev/null +++ b/tests/multigrid/transfer_matrix_free_13.cc @@ -0,0 +1,98 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Check MGTransferMatrixFree with custom user constraints + +#include + +#include + +#include + +#include + +#include +#include + +#include "../tests.h" + +template +void +check() +{ + Triangulation tr(Triangulation::limit_level_difference_at_vertices); + GridGenerator::hyper_cube(tr, 0, 1, false); + + typename Triangulation::active_cell_iterator cell = tr.begin_active(); + cell->face(0)->set_all_boundary_ids(1); + + tr.refine_global(1); + + FE_Q fe = FE_Q(1); + + DoFHandler mgdof(tr); + mgdof.distribute_dofs(fe); + mgdof.distribute_mg_dofs(); + + MGConstrainedDoFs mg_constrained_dofs; + mg_constrained_dofs.initialize(mgdof); + + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_level_dofs(mgdof, 0, relevant_dofs); + AffineConstraints user_constraints; + user_constraints.reinit(relevant_dofs); + + typename DoFHandler::level_face_iterator face0 = mgdof.begin(0)->face(0); + std::vector face_dofs(fe.dofs_per_face); + face0->get_mg_dof_indices(0, face_dofs); + for (unsigned int i = 0; i < face_dofs.size(); ++i) + { + if (user_constraints.can_store_line(face_dofs[i])) + { + user_constraints.add_line(face_dofs[i]); + user_constraints.set_inhomogeneity(face_dofs[i], 5.0); + } + } + user_constraints.close(); + mg_constrained_dofs.add_user_constraints(0, user_constraints); + + // build matrix-free transfer + MGTransferMatrixFree transfer(mg_constrained_dofs); + transfer.build(mgdof); + + LinearAlgebra::distributed::Vector src_level_0(mgdof.n_dofs(0)); + for (unsigned int i = 0; i < mgdof.n_dofs(0); ++i) + deallog << src_level_0(i) << " "; + deallog << std::endl; + + LinearAlgebra::distributed::Vector dst_level_1(mgdof.n_dofs(1)); + transfer.prolongate(1, dst_level_1, src_level_0); + for (unsigned int i = 0; i < mgdof.n_dofs(1); ++i) + deallog << dst_level_1(i) << " "; + deallog << std::endl; +} + + +int +main() +{ + initlog(); + + check<2>(); + deallog << std::endl; + check<3>(); + deallog << std::endl; +} diff --git a/tests/multigrid/transfer_matrix_free_13.output b/tests/multigrid/transfer_matrix_free_13.output new file mode 100644 index 0000000000..f51ed58732 --- /dev/null +++ b/tests/multigrid/transfer_matrix_free_13.output @@ -0,0 +1,7 @@ + +DEAL::0.00000 0.00000 0.00000 0.00000 +DEAL::5.00000 2.50000 5.00000 2.50000 0.00000 0.00000 5.00000 2.50000 0.00000 +DEAL:: +DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +DEAL::5.00000 2.50000 5.00000 2.50000 5.00000 2.50000 5.00000 2.50000 0.00000 0.00000 0.00000 0.00000 5.00000 2.50000 5.00000 2.50000 0.00000 0.00000 5.00000 2.50000 5.00000 2.50000 0.00000 0.00000 5.00000 2.50000 0.00000 +DEAL:: -- 2.39.5