From 64807d9142d5e942403ab3082a70a6fb219188e6 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 26 Feb 2021 18:59:30 +0100 Subject: [PATCH] Add MGTransferGlobalCoarsening::interpolate_to_mg() --- .../multigrid/mg_transfer_global_coarsening.h | 76 +++++ .../mg_transfer_global_coarsening.templates.h | 304 +++++++++++++++--- .../multigrid/mg_transfer_matrix_free.h | 6 +- .../multigrid_a_01.cc | 30 +- .../multigrid_p_01.cc | 30 +- 5 files changed, 383 insertions(+), 63 deletions(-) diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.h index 6a5063ef10..ea86730e57 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.h @@ -122,6 +122,15 @@ public: */ void restrict_and_add(VectorType &dst, const VectorType &src) const; + + /** + * Perform interpolation of a solution vector from the fine level to the + * coarse level. This function is different from restriction, where a + * weighted residual is transferred to a coarser level (transposition of + * prolongation matrix). + */ + void + interpolate(VectorType &dst, const VectorType &src) const; }; @@ -184,6 +193,16 @@ public: restrict_and_add(LinearAlgebra::distributed::Vector & dst, const LinearAlgebra::distributed::Vector &src) const; + /** + * Perform interpolation of a solution vector from the fine level to the + * coarse level. This function is different from restriction, where a + * weighted residual is transferred to a coarser level (transposition of + * prolongation matrix). + */ + void + interpolate(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const; + private: /** * A multigrid transfer scheme. A multrigrid transfer class can have different @@ -241,6 +260,16 @@ private: */ AlignedVector> prolongation_matrix_1d; + /** + * Restriction matrix for non-tensor-product elements. + */ + AlignedVector> restriction_matrix; + + /** + * 1D restriction matrix for tensor-product elements. + */ + AlignedVector> restriction_matrix_1d; + /** * DoF indices of the coarse cells, expressed in indices local to the MPI * rank. @@ -366,6 +395,26 @@ public: OutVector & dst, const MGLevelObject &src) const; + /** + * Interpolate fine-mesh field @p src to each multigrid level in + * @p dof_handler and store the result in @p dst. This function is different + * from restriction, where a weighted residual is + * transferred to a coarser level (transposition of prolongation matrix). + * + * The argument @p dst has to be initialized with the correct size according + * to the number of levels of the triangulation. + * + * If an inner vector of @p dst is empty or has incorrect locally owned size, + * it will be resized to locally relevant degrees of freedom on each level. + * + * @note DoFHandler is not needed here, but is required by the interface. + */ + template + void + interpolate_to_mg(const DoFHandler &dof_handler, + MGLevelObject & dst, + const InVector & src) const; + private: /** * Collection of the two-level transfer operators. @@ -453,6 +502,33 @@ MGTransferGlobalCoarsening::copy_from_mg( dst.copy_locally_owned_data_from(src[src.max_level()]); } + + +template +template +void +MGTransferGlobalCoarsening::interpolate_to_mg( + const DoFHandler &dof_handler, + MGLevelObject & dst, + const InVector & src) const +{ + (void)dof_handler; + + const unsigned int min_level = transfer.min_level(); + const unsigned int max_level = transfer.max_level(); + + AssertDimension(min_level, dst.min_level()); + AssertDimension(max_level, dst.max_level()); + + for (unsigned int level = min_level; level <= max_level; ++level) + initialize_dof_vector(level, dst[level]); + + dst[transfer.max_level()].copy_locally_owned_data_from(src); + + for (unsigned int l = max_level; l > min_level; --l) + this->transfer[l].interpolate(dst[l - 1], dst[l]); +} + #endif DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h index da1c32f895..840cbb8a18 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h @@ -381,6 +381,34 @@ namespace internal return FETools::get_fe_by_name<1, 1>(fe_name); } + template + FullMatrix + get_restriction_matrix(const FiniteElement &fe, + const unsigned int child) + { + auto matrix = fe.get_restriction_matrix(child); + + for (unsigned int c_other = 0; c_other < child; ++c_other) + { + auto matrix_other = fe.get_restriction_matrix(c_other); + for (unsigned int i = 0; i < fe.dofs_per_cell; ++i) + { + if (fe.restriction_is_additive(i) == true) + continue; + + bool do_zero = false; + for (unsigned int j = 0; j < fe.dofs_per_cell; ++j) + if (matrix_other(i, j) != 0.) + do_zero = true; + + if (do_zero) + for (unsigned int j = 0; j < fe.dofs_per_cell; ++j) + matrix(i, j) = 0.0; + } + } + return matrix; + } + } // namespace class FineDoFHandlerViewCell @@ -1238,38 +1266,78 @@ namespace internal const unsigned int shift = fe->dofs_per_cell; const unsigned int n_child_dofs_1d = fe->dofs_per_cell * 2; - transfer.schemes[1].prolongation_matrix_1d.resize( - fe->dofs_per_cell * n_child_dofs_1d); - - for (unsigned int c = 0; c < GeometryInfo<1>::max_children_per_cell; - ++c) - for (unsigned int i = 0; i < fe->dofs_per_cell; ++i) - for (unsigned int j = 0; j < fe->dofs_per_cell; ++j) - transfer.schemes[1] - .prolongation_matrix_1d[i * n_child_dofs_1d + j + - c * shift] = - fe->get_prolongation_matrix(c)(renumbering[j], - renumbering[i]); + { + transfer.schemes[1].prolongation_matrix_1d.resize( + fe->dofs_per_cell * n_child_dofs_1d); + + for (unsigned int c = 0; + c < GeometryInfo<1>::max_children_per_cell; + ++c) + for (unsigned int i = 0; i < fe->dofs_per_cell; ++i) + for (unsigned int j = 0; j < fe->dofs_per_cell; ++j) + transfer.schemes[1] + .prolongation_matrix_1d[i * n_child_dofs_1d + j + + c * shift] = + fe->get_prolongation_matrix(c)(renumbering[j], + renumbering[i]); + } + { + transfer.schemes[1].restriction_matrix_1d.resize( + fe->dofs_per_cell * n_child_dofs_1d); + + for (unsigned int c = 0; + c < GeometryInfo<1>::max_children_per_cell; + ++c) + { + const auto matrix = get_restriction_matrix(*fe, c); + for (unsigned int i = 0; i < fe->dofs_per_cell; ++i) + for (unsigned int j = 0; j < fe->dofs_per_cell; ++j) + transfer.schemes[1] + .restriction_matrix_1d[i * n_child_dofs_1d + j + + c * shift] = + matrix(renumbering[i], renumbering[j]); + } + } } else { const auto & fe = fe_fine.base_element(0); const unsigned int n_dofs_per_cell = fe.n_dofs_per_cell(); - transfer.schemes[1].prolongation_matrix.resize( - n_dofs_per_cell * n_dofs_per_cell * - GeometryInfo::max_children_per_cell); - - for (unsigned int c = 0; - c < GeometryInfo::max_children_per_cell; - ++c) - for (unsigned int i = 0; i < n_dofs_per_cell; ++i) - for (unsigned int j = 0; j < n_dofs_per_cell; ++j) - transfer.schemes[1].prolongation_matrix - [i * n_dofs_per_cell * - GeometryInfo::max_children_per_cell + - j + c * n_dofs_per_cell] = - fe.get_prolongation_matrix(c)(j, i); + { + transfer.schemes[1].prolongation_matrix.resize( + n_dofs_per_cell * n_dofs_per_cell * + GeometryInfo::max_children_per_cell); + + for (unsigned int c = 0; + c < GeometryInfo::max_children_per_cell; + ++c) + for (unsigned int i = 0; i < n_dofs_per_cell; ++i) + for (unsigned int j = 0; j < n_dofs_per_cell; ++j) + transfer.schemes[1].prolongation_matrix + [i * n_dofs_per_cell * + GeometryInfo::max_children_per_cell + + j + c * n_dofs_per_cell] = + fe.get_prolongation_matrix(c)(j, i); + } + { + transfer.schemes[1].restriction_matrix.resize( + n_dofs_per_cell * n_dofs_per_cell * + GeometryInfo::max_children_per_cell); + + for (unsigned int c = 0; + c < GeometryInfo::max_children_per_cell; + ++c) + { + const auto matrix = get_restriction_matrix(fe, c); + for (unsigned int i = 0; i < n_dofs_per_cell; ++i) + for (unsigned int j = 0; j < n_dofs_per_cell; ++j) + transfer.schemes[1].restriction_matrix + [i * n_dofs_per_cell * + GeometryInfo::max_children_per_cell + + j + c * n_dofs_per_cell] = matrix(i, j); + } + } } } @@ -1658,16 +1726,33 @@ namespace internal fe_coarse->dofs_per_vertex; } - FullMatrix matrix(fe_fine->dofs_per_cell, - fe_coarse->dofs_per_cell); - FETools::get_projection_matrix(*fe_coarse, *fe_fine, matrix); - transfer.schemes[fe_index_no].prolongation_matrix_1d.resize( - fe_fine->dofs_per_cell * fe_coarse->dofs_per_cell); + { + FullMatrix matrix(fe_fine->dofs_per_cell, + fe_coarse->dofs_per_cell); + FETools::get_projection_matrix(*fe_coarse, *fe_fine, matrix); + transfer.schemes[fe_index_no].prolongation_matrix_1d.resize( + fe_fine->dofs_per_cell * fe_coarse->dofs_per_cell); + + for (unsigned int i = 0, k = 0; i < fe_coarse->dofs_per_cell; + ++i) + for (unsigned int j = 0; j < fe_fine->dofs_per_cell; ++j, ++k) + transfer.schemes[fe_index_no].prolongation_matrix_1d[k] = + matrix(renumbering_fine[j], renumbering_coarse[i]); + } - for (unsigned int i = 0, k = 0; i < fe_coarse->dofs_per_cell; ++i) - for (unsigned int j = 0; j < fe_fine->dofs_per_cell; ++j, ++k) - transfer.schemes[fe_index_no].prolongation_matrix_1d[k] = - matrix(renumbering_fine[j], renumbering_coarse[i]); + { + FullMatrix matrix(fe_coarse->dofs_per_cell, + fe_fine->dofs_per_cell); + FETools::get_projection_matrix(*fe_fine, *fe_coarse, matrix); + transfer.schemes[fe_index_no].restriction_matrix_1d.resize( + fe_fine->dofs_per_cell * fe_coarse->dofs_per_cell); + + for (unsigned int i = 0, k = 0; i < fe_coarse->dofs_per_cell; + ++i) + for (unsigned int j = 0; j < fe_fine->dofs_per_cell; ++j, ++k) + transfer.schemes[fe_index_no].restriction_matrix_1d[k] = + matrix(renumbering_coarse[i], renumbering_fine[j]); + } } else { @@ -1677,16 +1762,33 @@ namespace internal const auto &fe_coarse = dof_handler_coarse.get_fe(fe_index_pair.first).base_element(0); - FullMatrix matrix(fe_fine.dofs_per_cell, - fe_coarse.dofs_per_cell); - FETools::get_projection_matrix(fe_coarse, fe_fine, matrix); - transfer.schemes[fe_index_no].prolongation_matrix.resize( - fe_fine.dofs_per_cell * fe_coarse.dofs_per_cell); + { + FullMatrix matrix(fe_fine.dofs_per_cell, + fe_coarse.dofs_per_cell); + FETools::get_projection_matrix(fe_coarse, fe_fine, matrix); + transfer.schemes[fe_index_no].prolongation_matrix.resize( + fe_fine.dofs_per_cell * fe_coarse.dofs_per_cell); + + for (unsigned int i = 0, k = 0; i < fe_coarse.dofs_per_cell; + ++i) + for (unsigned int j = 0; j < fe_fine.dofs_per_cell; ++j, ++k) + transfer.schemes[fe_index_no].prolongation_matrix[k] = + matrix(j, i); + } - for (unsigned int i = 0, k = 0; i < fe_coarse.dofs_per_cell; ++i) - for (unsigned int j = 0; j < fe_fine.dofs_per_cell; ++j, ++k) - transfer.schemes[fe_index_no].prolongation_matrix[k] = - matrix(j, i); + { + FullMatrix matrix(fe_coarse.dofs_per_cell, + fe_fine.dofs_per_cell); + FETools::get_projection_matrix(fe_fine, fe_coarse, matrix); + transfer.schemes[fe_index_no].restriction_matrix.resize( + fe_fine.dofs_per_cell * fe_coarse.dofs_per_cell); + + for (unsigned int i = 0, k = 0; i < fe_coarse.dofs_per_cell; + ++i) + for (unsigned int j = 0; j < fe_fine.dofs_per_cell; ++j, ++k) + transfer.schemes[fe_index_no].restriction_matrix[k] = + matrix(i, j); + } } } @@ -2155,7 +2257,6 @@ MGTwoLevelTransfer>:: } // ----------------------------- coarse ---------------------------- - // write into dst vector { const unsigned int *indices = @@ -2181,6 +2282,119 @@ MGTwoLevelTransfer>:: +template +void +MGTwoLevelTransfer>:: + interpolate(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const +{ + using VectorizedArrayType = VectorizedArray; + + const unsigned int n_lanes = VectorizedArrayType::size(); + + this->vec_fine.copy_locally_owned_data_from(src); + this->vec_fine.update_ghost_values(); + + this->vec_coarse = 0.0; + + AlignedVector evaluation_data_fine; + AlignedVector evaluation_data_coarse; + + for (const auto &scheme : schemes) + { + // identity -> take short cut and work directly on global vectors + if (scheme.degree_fine == scheme.degree_coarse) + { + const unsigned int *indices_fine = + scheme.level_dof_indices_fine.data(); + const unsigned int *indices_coarse = + scheme.level_dof_indices_coarse.data(); + + for (unsigned int cell = 0; cell < scheme.n_coarse_cells; ++cell) + { + for (unsigned int i = 0; i < scheme.dofs_per_cell_fine; ++i) + this->vec_coarse.local_element(indices_coarse[i]) = + this->vec_fine.local_element(indices_fine[i]); + + indices_fine += scheme.dofs_per_cell_fine; + indices_coarse += scheme.dofs_per_cell_coarse; + } + + continue; + } + + // general case -> local restriction is needed + evaluation_data_fine.resize(scheme.dofs_per_cell_fine); + evaluation_data_coarse.resize(scheme.dofs_per_cell_fine); + + CellTransferFactory cell_transfer(scheme.degree_fine, + scheme.degree_coarse); + + const unsigned int n_scalar_dofs_fine = + scheme.dofs_per_cell_fine / n_components; + const unsigned int n_scalar_dofs_coarse = + scheme.dofs_per_cell_coarse / n_components; + + for (unsigned int cell = 0; cell < scheme.n_coarse_cells; cell += n_lanes) + { + const unsigned int n_lanes_filled = + (cell + n_lanes > scheme.n_coarse_cells) ? + (scheme.n_coarse_cells - cell) : + n_lanes; + + // read from source vector and weight + { + const unsigned int *indices = + &scheme.level_dof_indices_fine[cell * scheme.dofs_per_cell_fine]; + + for (unsigned int v = 0; v < n_lanes_filled; ++v) + { + for (unsigned int i = 0; i < scheme.dofs_per_cell_fine; ++i) + evaluation_data_fine[i][v] = + this->vec_fine.local_element(indices[i]); + + indices += scheme.dofs_per_cell_fine; + } + } + + // ------------------------------ fine ----------------------------- + for (int c = n_components - 1; c >= 0; --c) + { + CellRestrictor cell_restrictor( + scheme.restriction_matrix, + scheme.restriction_matrix_1d, + evaluation_data_fine.begin() + c * n_scalar_dofs_fine, + evaluation_data_coarse.begin() + c * n_scalar_dofs_coarse); + + if (scheme.restriction_matrix_1d.size() > 0) + cell_transfer.run(cell_restrictor); + else + cell_restrictor.run_full(n_scalar_dofs_fine, + n_scalar_dofs_coarse); + } + // ----------------------------- coarse ---------------------------- + + // write into dst vector + { + const unsigned int *indices = + &scheme + .level_dof_indices_coarse[cell * scheme.dofs_per_cell_coarse]; + for (unsigned int v = 0; v < n_lanes_filled; ++v) + { + for (unsigned int i = 0; i < scheme.dofs_per_cell_coarse; ++i) + this->vec_coarse.local_element(indices[i]) = + evaluation_data_coarse[i][v]; + indices += scheme.dofs_per_cell_coarse; + } + } + } + } + + dst.copy_locally_owned_data_from(this->vec_coarse); +} + + + template void MGTwoLevelTransfer>:: diff --git a/include/deal.II/multigrid/mg_transfer_matrix_free.h b/include/deal.II/multigrid/mg_transfer_matrix_free.h index be7d526ea5..6b620c0ce7 100644 --- a/include/deal.II/multigrid/mg_transfer_matrix_free.h +++ b/include/deal.II/multigrid/mg_transfer_matrix_free.h @@ -156,8 +156,10 @@ public: const LinearAlgebra::distributed::Vector &src) const override; /** - * Restrict fine-mesh field @p src to each multigrid level in @p dof_handler and - * store the result in @p dst. + * Interpolate fine-mesh field @p src to each multigrid level in + * @p dof_handler and store the result in @p dst. This function is different + * from restriction, where a weighted residual is + * transferred to a coarser level (transposition of prolongation matrix). * * The argument @p dst has to be initialized with the correct size according * to the number of levels of the triangulation. diff --git a/tests/multigrid-global-coarsening/multigrid_a_01.cc b/tests/multigrid-global-coarsening/multigrid_a_01.cc index 42f5eec640..def66ba1c2 100644 --- a/tests/multigrid-global-coarsening/multigrid_a_01.cc +++ b/tests/multigrid-global-coarsening/multigrid_a_01.cc @@ -128,16 +128,30 @@ test(const unsigned int n_refinements, << (do_simplex_mesh ? "tri " : "quad") << " " << solver_control.last_step() << std::endl; - DataOut data_out; + static unsigned int counter = 0; - data_out.attach_dof_handler(dof_handlers[max_level]); - data_out.add_data_vector(dst, "solution"); - data_out.build_patches(*mapping_, 2); + MGLevelObject results(min_level, max_level); - static unsigned int counter = 0; - std::ofstream output("test." + std::to_string(dim) + "." + - std::to_string(counter++) + ".vtk"); - data_out.write_vtk(output); + transfer.interpolate_to_mg(dof_handlers[max_level], results, dst); + + for (unsigned int l = min_level; l <= max_level; ++l) + { + DataOut data_out; + + data_out.attach_dof_handler(dof_handlers[l]); + data_out.add_data_vector( + results[l], + "solution", + DataOut_DoFData, dim>::DataVectorType::type_dof_data); + data_out.build_patches(*mapping_, 2); + + std::ofstream output("test." + std::to_string(dim) + "." + + std::to_string(counter) + "." + std::to_string(l) + + ".vtk"); + data_out.write_vtk(output); + } + + counter++; } int diff --git a/tests/multigrid-global-coarsening/multigrid_p_01.cc b/tests/multigrid-global-coarsening/multigrid_p_01.cc index 945ded4808..acbfb4918c 100644 --- a/tests/multigrid-global-coarsening/multigrid_p_01.cc +++ b/tests/multigrid-global-coarsening/multigrid_p_01.cc @@ -158,16 +158,30 @@ test(const unsigned int n_refinements, << (do_simplex_mesh ? "tri " : "quad") << " " << solver_control.last_step() << std::endl; - DataOut data_out; + static unsigned int counter = 0; - data_out.attach_dof_handler(dof_handlers[max_level]); - data_out.add_data_vector(dst, "solution"); - data_out.build_patches(*mapping_, 2); + MGLevelObject results(min_level, max_level); - static unsigned int counter = 0; - std::ofstream output("test." + std::to_string(dim) + "." + - std::to_string(counter++) + ".vtk"); - data_out.write_vtk(output); + transfer.interpolate_to_mg(dof_handlers[max_level], results, dst); + + for (unsigned int l = min_level; l <= max_level; ++l) + { + DataOut data_out; + + data_out.attach_dof_handler(dof_handlers[l]); + data_out.add_data_vector( + results[l], + "solution", + DataOut_DoFData, dim>::DataVectorType::type_dof_data); + data_out.build_patches(*mapping_, 2); + + std::ofstream output("test." + std::to_string(dim) + "." + + std::to_string(counter) + "." + std::to_string(l) + + ".vtk"); + data_out.write_vtk(output); + } + + counter++; } int -- 2.39.5