From a3322a2e337011f81f4adb2c7ace2b8f80ca7da0 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Thu, 28 Sep 2023 16:46:26 +0200 Subject: [PATCH] MGTransferMF: check compatibility of DoFHandlers --- doc/news/changes/minor/20231003Munch | 7 + .../multigrid/mg_transfer_global_coarsening.h | 120 ++++++-- .../mg_transfer_global_coarsening.templates.h | 282 ++++++++++++++++-- .../multigrid_p_04.cc | 250 ++++++++++++++++ ...id_p_04.mpirun=1.with_trilinos=true.output | 79 +++++ 5 files changed, 694 insertions(+), 44 deletions(-) create mode 100644 doc/news/changes/minor/20231003Munch create mode 100644 tests/multigrid-global-coarsening/multigrid_p_04.cc create mode 100644 tests/multigrid-global-coarsening/multigrid_p_04.mpirun=1.with_trilinos=true.output diff --git a/doc/news/changes/minor/20231003Munch b/doc/news/changes/minor/20231003Munch new file mode 100644 index 0000000000..8bee73ba61 --- /dev/null +++ b/doc/news/changes/minor/20231003Munch @@ -0,0 +1,7 @@ +Improved: MGTransferMF can now renumerate DoFs during +MGTransferMF::copy_to_mg(), MGTransferMF::copy_to_mg(), +and MGTransferMF::interpolate_to_mg() if the outer +solver and the multigrid preconditioner have been set up +with different DoFHandler instances. +
+(Peter Munch, Laura Prieto Saavedra, 2023/10/03) diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.h index 68521ff6b2..187a5c7e4e 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.h @@ -49,6 +49,9 @@ namespace RepartitioningPolicyTools template class Base; } + +template +class MGTransferMF; #endif @@ -706,7 +709,19 @@ private: */ unsigned int n_components; + /** + * Pointer to the DoFHandler object used during initialization. + */ + SmartPointer> dof_handler_fine; + + /** + * Muligird level used during initialization. + */ + unsigned int mg_level_fine; + friend class internal::MGTwoLevelTransferImplementation; + + friend class MGTransferMF; }; @@ -919,6 +934,16 @@ private: LinearAlgebra::distributed::Vector &dst, const LinearAlgebra::distributed::Vector &src) const; + /** + * Pointer to the DoFHandler object used during initialization. + */ + SmartPointer> dof_handler_fine; + + /** + * Multigrid level used during initialization. + */ + unsigned int mg_level_fine; + /** * Object to evaluate shape functions on one mesh on visited support points of * the other mesh. @@ -954,6 +979,8 @@ private: * point. */ std::vector level_dof_indices_fine_ptrs; + + friend class MGTransferMF; }; @@ -1001,7 +1028,7 @@ public: * * @note See also MGTransferMatrixFree. */ - MGTransferMF() = default; + MGTransferMF(); /** * @name Global coarsening. @@ -1073,10 +1100,20 @@ public: void initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs); + /** @} */ + + /** + * @name Global coarsening and local smoothing. + */ + /** @{ */ + /** * Actually build the information for the prolongation for each level. * - * @note See also MGTransferMatrixFree. + * @note In the case of global coarsening, you can pass into this function + * a @p dof_handler with different DoF numbering as the one used within the + * provided two-level transfer objects. In this case, vector entries are + * permuted during copy_to_mg(), copy_from_mg(), and interpolate_to_mg(). */ void build(const DoFHandler &dof_handler, @@ -1087,7 +1124,10 @@ public: * Same as above but taking a lambda for initializing vector instead of * partitioners. * - * @note See also MGTransferMatrixFree. + * @note In the case of global coarsening, you can pass into this function + * a @p dof_handler with different DoF numbering as the one used within the + * provided two-level transfer objects. In this case, vector entries are + * permuted during copy_to_mg(), copy_from_mg(), and interpolate_to_mg(). */ void build(const DoFHandler &dof_handler, @@ -1163,19 +1203,17 @@ public: */ template void - interpolate_to_mg(MGLevelObject &dst, const InVector &src) const; + interpolate_to_mg(const DoFHandler &dof_handler, + MGLevelObject &dst, + const InVector &src) const; /** - * Like the above function but with a user-provided DoFHandler as - * additional argument. However, this DoFHandler is not used internally, but - * is required to be able to use MGTransferMF and - * MGTransferMatrixFree as template argument. + * Interpolate fine-mesh field @p src to each multigrid level in + * @p dof_handler and store the result in @p dst. */ template void - interpolate_to_mg(const DoFHandler &dof_handler, - MGLevelObject &dst, - const InVector &src) const; + interpolate_to_mg(MGLevelObject &dst, const InVector &src) const; /** @} */ @@ -1225,6 +1263,21 @@ private: const DoFHandler &dof_handler, const SmartPointer &mg_constrained_dofs); + /** + * Retreave finest DoFHandler from two-level transfer objects. + */ + std::pair *, unsigned int> + get_dof_handler_fine() const; + + /** + * Initialize copy indices for MGTransferMF::copy_to_mg(), + * MGTransferMF::copy_to_mg(), and MGTransferMF::interpolate_to_mg() + * in the case of global coarsening. + */ + void + fill_and_communicate_copy_indices_global_coarsening( + const DoFHandler &dof_handler); + /** * Set references to two-level transfer operators to be used. */ @@ -1243,6 +1296,14 @@ private: const InVector &vector_reference, const bool omit_zeroing_entries) const; + /** + * Check that the internal DoFHandler is compatible with the external one + * provided by copy_to_mg(), copy_from_mg() and interpolate_to_mg() + * used, e.g., by PreconditionMG. + */ + void + assert_dof_handler(const DoFHandler &dof_handler_out) const; + /** * Internal transfer operator. * @@ -1374,6 +1435,9 @@ MGTransferMF::MGTransferMF( const std::function &initialize_dof_vector) { + this->transfer.clear(); + this->internal_transfer.clear(); + this->initialize_transfer_references(transfer); this->build(initialize_dof_vector); } @@ -1465,7 +1529,7 @@ MGTransferMF::copy_to_mg(const DoFHandler &dof_handler, MGLevelObject &dst, const InVector &src) const { - (void)dof_handler; + assert_dof_handler(dof_handler); for (unsigned int level = dst.min_level(); level <= dst.max_level(); ++level) { @@ -1524,7 +1588,7 @@ MGTransferMF::copy_from_mg( OutVector &dst, const MGLevelObject &src) const { - (void)dof_handler; + assert_dof_handler(dof_handler); if (this->perform_plain_copy) { @@ -1576,6 +1640,22 @@ void MGTransferMF::interpolate_to_mg(MGLevelObject &dst, const InVector &src) const { + DoFHandler dof_handler_dummy; + + this->interpolate_to_mg(dof_handler_dummy, dst, src); +} + + + +template +template +void +MGTransferMF::interpolate_to_mg(const DoFHandler &dof_handler, + MGLevelObject &dst, + const InVector &src) const +{ + assert_dof_handler(dof_handler); + const unsigned int min_level = transfer.min_level(); const unsigned int max_level = transfer.max_level(); @@ -1636,20 +1716,6 @@ MGTransferMF::interpolate_to_mg(MGLevelObject &dst, } } - - -template -template -void -MGTransferMF::interpolate_to_mg(const DoFHandler &dof_handler, - MGLevelObject &dst, - const InVector &src) const -{ - (void)dof_handler; - - this->interpolate_to_mg(dst, src); -} - #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 903f620d4b..16c1667062 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h @@ -1774,17 +1774,8 @@ namespace internal (mg_level_coarse + 1 == mg_level_fine), ExcNotImplemented()); - // if (mg_level_fine != numbers::invalid_unsigned_int) - // AssertIndexRange(mg_level_fine, - // MGTools::max_level_for_coarse_mesh( - // dof_handler_fine.get_triangulation()) + - // 1); - // - // if (mg_level_coarse != numbers::invalid_unsigned_int) - // AssertIndexRange(mg_level_coarse, - // MGTools::max_level_for_coarse_mesh( - // dof_handler_coarse.get_triangulation()) + - // 1); + transfer.dof_handler_fine = &dof_handler_fine; + transfer.mg_level_fine = mg_level_fine; std::unique_ptr> dof_handler_fine_view; @@ -2326,6 +2317,9 @@ namespace internal "(numbers::invalid_unsigned_int) or on refinement levels without " "hanging nodes.")); + transfer.dof_handler_fine = &dof_handler_fine; + transfer.mg_level_fine = mg_level_fine; + std::unique_ptr> dof_handler_fine_view; if (internal::p_transfer_involves_repartitioning(dof_handler_fine, @@ -4051,10 +4045,21 @@ MGTwoLevelTransferBase>:: +template +MGTransferMF::MGTransferMF() +{ + this->transfer.clear(); + this->internal_transfer.clear(); +} + + + template MGTransferMF::MGTransferMF( const MGConstrainedDoFs &mg_constrained_dofs) { + this->transfer.clear(); + this->internal_transfer.clear(); this->initialize_constraints(mg_constrained_dofs); } @@ -4102,6 +4107,135 @@ MGTransferMF::initialize_internal_transfer( +template +std::pair *, unsigned int> +MGTransferMF::get_dof_handler_fine() const +{ + if (this->transfer.n_levels() <= + 1) // single level: the information cannot be retreaved + return {nullptr, numbers::invalid_unsigned_int}; + + if (const auto t = dynamic_cast< + const MGTwoLevelTransfer> *>( + this->transfer[this->transfer.max_level()].get())) + { + return {t->dof_handler_fine, t->mg_level_fine}; + } + else if (const auto t = dynamic_cast> *>( + this->transfer[this->transfer.max_level()].get())) + { + return {t->dof_handler_fine, t->mg_level_fine}; + } + else + { + Assert(false, ExcNotImplemented()); + return {nullptr, numbers::invalid_unsigned_int}; + } +} + + + +template +void +MGTransferMF::fill_and_communicate_copy_indices_global_coarsening( + const DoFHandler &dof_handler_out) +{ + const auto [dof_handler_in, level_in] = get_dof_handler_fine(); + + if ((dof_handler_in == nullptr) || (dof_handler_in == &dof_handler_out)) + return; // nothing to do + + this->copy_indices.resize(1); + this->copy_indices[0].reinit(2, dof_handler_out.n_locally_owned_dofs()); + + std::vector dof_indices_in; + std::vector dof_indices_out; + + this->perform_plain_copy = true; + + const auto &is_out = (level_in == numbers::invalid_unsigned_int) ? + dof_handler_out.locally_owned_dofs() : + dof_handler_out.locally_owned_mg_dofs(level_in); + + const auto &is_in = (level_in == numbers::invalid_unsigned_int) ? + dof_handler_in->locally_owned_dofs() : + dof_handler_in->locally_owned_mg_dofs(level_in); + + internal::loop_over_active_or_level_cells( + dof_handler_in->get_triangulation(), level_in, [&](const auto &cell) { + const auto cell_id = cell->id(); + + Assert( + dof_handler_out.get_triangulation().contains_cell(cell_id), + ExcMessage( + "DoFHandler instances used for set up of MGTransferMF and copy_to_mg(), " + "copy_from_mg(), or interpolate_to_mg() are not compatible.")); + + if (level_in == numbers::invalid_unsigned_int) + { + const auto cell_in = cell->as_dof_handler_iterator(*dof_handler_in); + const auto cell_out = dof_handler_out.get_triangulation() + .create_cell_iterator(cell_id) + ->as_dof_handler_iterator(dof_handler_out); + + AssertDimension(cell_in->get_fe().n_dofs_per_cell(), + cell_out->get_fe().n_dofs_per_cell()); + + dof_indices_in.resize(cell_in->get_fe().n_dofs_per_cell()); + dof_indices_out.resize(cell_out->get_fe().n_dofs_per_cell()); + + cell_in->get_dof_indices(dof_indices_in); + cell_out->get_dof_indices(dof_indices_out); + } + else + { + const auto cell_in = + cell->as_dof_handler_level_iterator(*dof_handler_in); + const auto cell_out = + dof_handler_out.get_triangulation() + .create_cell_iterator(cell_id) + ->as_dof_handler_level_iterator(dof_handler_out); + + AssertDimension(cell_in->get_fe().n_dofs_per_cell(), + cell_out->get_fe().n_dofs_per_cell()); + + dof_indices_in.resize(cell_in->get_fe().n_dofs_per_cell()); + dof_indices_out.resize(cell_out->get_fe().n_dofs_per_cell()); + + cell_in->get_mg_dof_indices(dof_indices_in); + cell_out->get_mg_dof_indices(dof_indices_out); + } + + this->perform_plain_copy &= (dof_indices_in == dof_indices_out); + + for (unsigned int i = 0; i < dof_indices_in.size(); ++i) + if (is_out.is_element(dof_indices_out[i])) + this->copy_indices[0](1, + is_out.index_within_set(dof_indices_out[i])) = + is_in.index_within_set(dof_indices_in[i]); + }); + + + this->perform_plain_copy = + Utilities::MPI::max(this->perform_plain_copy ? 1 : 0, + dof_handler_out.get_communicator()) != 0; + + if (this->perform_plain_copy) + { + this->copy_indices.clear(); + } + else + { + this->perform_renumbered_plain_copy = true; + this->solution_copy_indices = this->copy_indices; + } +} + + + template void MGTransferMF::build( @@ -4183,10 +4317,22 @@ MGTransferMF::build( const std::vector> &external_partitioners) { - this->initialize_internal_transfer(dof_handler, this->mg_constrained_dofs); - this->initialize_transfer_references(internal_transfer); + const bool use_local_smoothing = + this->transfer.n_levels() == 0 || this->internal_transfer.n_levels() > 0; + + if (use_local_smoothing) + { + this->initialize_internal_transfer(dof_handler, + this->mg_constrained_dofs); + this->initialize_transfer_references(internal_transfer); + } + this->build(external_partitioners); - this->fill_and_communicate_copy_indices(dof_handler); + + if (use_local_smoothing) + this->fill_and_communicate_copy_indices(dof_handler); + else + this->fill_and_communicate_copy_indices_global_coarsening(dof_handler); } @@ -4198,10 +4344,22 @@ MGTransferMF::build( const std::function &initialize_dof_vector) { - this->initialize_internal_transfer(dof_handler, this->mg_constrained_dofs); - this->initialize_transfer_references(internal_transfer); + const bool use_local_smoothing = + this->transfer.n_levels() == 0 || this->internal_transfer.n_levels() > 0; + + if (use_local_smoothing) + { + this->initialize_internal_transfer(dof_handler, + this->mg_constrained_dofs); + this->initialize_transfer_references(internal_transfer); + } + this->build(initialize_dof_vector); - this->fill_and_communicate_copy_indices(dof_handler); + + if (use_local_smoothing) + this->fill_and_communicate_copy_indices(dof_handler); + else + this->fill_and_communicate_copy_indices_global_coarsening(dof_handler); } @@ -4240,6 +4398,93 @@ MGTransferMF::restrict_and_add(const unsigned int from_level, +template +void +MGTransferMF::assert_dof_handler( + const DoFHandler &dof_handler_out) const +{ +#ifndef DEBUG + (void)dof_handler_out; +#else + + const auto [dof_handler_in, level_in] = get_dof_handler_fine(); + + if ((dof_handler_out.n_dofs() == 0) || // dummy DoFHandler + (dof_handler_in == nullptr) || // single level + (dof_handler_in == &dof_handler_out // same DoFHandler + )) + return; // nothing to do + + if (this->perform_plain_copy) + { + // global-coarsening path: compare indices of cells + + std::vector dof_indices_in; + std::vector dof_indices_out; + + internal::loop_over_active_or_level_cells( + dof_handler_in->get_triangulation(), level_in, [&](const auto &cell) { + const auto cell_id = cell->id(); + + Assert( + dof_handler_out.get_triangulation().contains_cell(cell_id), + ExcMessage( + "DoFHandler instances used for set up of MGTransferMF and copy_to_mg(), " + "copy_from_mg(), or interpolate_to_mg() are not compatible.")); + + if (level_in == numbers::invalid_unsigned_int) + { + const auto cell_in = + cell->as_dof_handler_iterator(*dof_handler_in); + const auto cell_out = + dof_handler_out.get_triangulation() + .create_cell_iterator(cell_id) + ->as_dof_handler_iterator(dof_handler_out); + + AssertDimension(cell_in->get_fe().n_dofs_per_cell(), + cell_out->get_fe().n_dofs_per_cell()); + + dof_indices_in.resize(cell_in->get_fe().n_dofs_per_cell()); + dof_indices_out.resize(cell_out->get_fe().n_dofs_per_cell()); + + cell_in->get_dof_indices(dof_indices_in); + cell_out->get_dof_indices(dof_indices_out); + } + else + { + const auto cell_in = + cell->as_dof_handler_level_iterator(*dof_handler_in); + const auto cell_out = + dof_handler_out.get_triangulation() + .create_cell_iterator(cell_id) + ->as_dof_handler_level_iterator(dof_handler_out); + + AssertDimension(cell_in->get_fe().n_dofs_per_cell(), + cell_out->get_fe().n_dofs_per_cell()); + + dof_indices_in.resize(cell_in->get_fe().n_dofs_per_cell()); + dof_indices_out.resize(cell_out->get_fe().n_dofs_per_cell()); + + cell_in->get_mg_dof_indices(dof_indices_in); + cell_out->get_mg_dof_indices(dof_indices_out); + } + + Assert( + dof_indices_in == dof_indices_out, + ExcMessage( + "DoFHandler instances used for set up of MGTransferMF and copy_to_mg(), " + "copy_from_mg(), or interpolate_to_mg() are not compatible.")); + }); + } + else if (this->perform_renumbered_plain_copy) + { + // nothing to do + } +#endif +} + + + template std::size_t MGTransferMF::memory_consumption() const @@ -4798,6 +5043,9 @@ MGTwoLevelTransferNonNested>:: dof_handler_fine.get_fe().n_components() > 0, ExcNotImplemented()); + this->dof_handler_fine = &dof_handler_fine; + this->mg_level_fine = numbers::invalid_unsigned_int; + this->fine_element_is_continuous = dof_handler_fine.get_fe().n_dofs_per_vertex() > 0; diff --git a/tests/multigrid-global-coarsening/multigrid_p_04.cc b/tests/multigrid-global-coarsening/multigrid_p_04.cc new file mode 100644 index 0000000000..c8fa658c7f --- /dev/null +++ b/tests/multigrid-global-coarsening/multigrid_p_04.cc @@ -0,0 +1,250 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 - 2022 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. +// +// --------------------------------------------------------------------- + + +/** + * Similar to test multigrid_01 but the DoFs are renumbered in the DoFHandler + * used by the outer solver. + */ + +#include + +#include "multigrid_util.h" + +template +void +test(const unsigned int n_refinements, + const unsigned int fe_degree_fine, + const bool do_simplex_mesh, + const unsigned int mesh_type) +{ + using VectorType = LinearAlgebra::distributed::Vector; + + Triangulation tria; + if (do_simplex_mesh) + GridGenerator::subdivided_hyper_cube_with_simplices(tria, 2); + else + GridGenerator::subdivided_hyper_cube(tria, 2); + + if (mesh_type == 0) + { + tria.refine_global(n_refinements); + } + else if (mesh_type == 1) + { + for (unsigned int i = 1; i < n_refinements; ++i) + { + for (auto cell : tria.active_cell_iterators()) + if (cell->is_locally_owned()) + { + bool flag = true; + for (int d = 0; d < dim; ++d) + if (cell->center()[d] > 0.5) + flag = false; + if (flag) + cell->set_refine_flag(); + } + tria.execute_coarsening_and_refinement(); + } + } + else + AssertThrow(false, ExcNotImplemented()); + + //////////////////////////////////////////////////////////////////////// + + std::unique_ptr> fe; + std::unique_ptr> quad; + std::unique_ptr> mapping; + + if (do_simplex_mesh) + { + fe = std::make_unique>(fe_degree_fine); + quad = std::make_unique>(fe_degree_fine + 1); + mapping = std::make_unique>(FE_SimplexP(1)); + } + else + { + fe = std::make_unique>(fe_degree_fine); + quad = std::make_unique>(fe_degree_fine + 1); + mapping = std::make_unique>(FE_Q(1)); + } + + DoFHandler fine_dof_handler(tria); + fine_dof_handler.distribute_dofs(*fe); + DoFRenumbering::Cuthill_McKee(fine_dof_handler); + + AffineConstraints fine_constraints; + const IndexSet locally_relevant_dofs = + DoFTools::extract_locally_relevant_dofs(fine_dof_handler); + fine_constraints.reinit(locally_relevant_dofs); + VectorTools::interpolate_boundary_values(*mapping, + fine_dof_handler, + 0, + Functions::ZeroFunction(), + fine_constraints); + DoFTools::make_hanging_node_constraints(fine_dof_handler, fine_constraints); + fine_constraints.close(); + + // set up operator + Operator fine_operator; + fine_operator.reinit(*mapping, fine_dof_handler, *quad, fine_constraints); + + //////////////////////////////////////////////////////////////////////// + + const auto level_degrees = + MGTransferGlobalCoarseningTools::create_polynomial_coarsening_sequence( + fe_degree_fine, + MGTransferGlobalCoarseningTools::PolynomialCoarseningSequenceType:: + bisect); + + const unsigned int min_level = 0; + const unsigned int max_level = level_degrees.size() - 1; + + MGLevelObject> dof_handlers(min_level, max_level, tria); + MGLevelObject> constraints(min_level, max_level); + MGLevelObject> transfers(min_level, + max_level); + MGLevelObject> operators(min_level, max_level); + + std::unique_ptr> mapping_; + + // set up levels + for (auto l = min_level; l <= max_level; ++l) + { + auto &dof_handler = dof_handlers[l]; + auto &constraint = constraints[l]; + auto &op = operators[l]; + + std::unique_ptr> fe; + std::unique_ptr> quad; + std::unique_ptr> mapping; + + if (do_simplex_mesh) + { + fe = std::make_unique>(level_degrees[l]); + quad = std::make_unique>(level_degrees[l] + 1); + mapping = std::make_unique>(FE_SimplexP(1)); + } + else + { + fe = std::make_unique>(level_degrees[l]); + quad = std::make_unique>(level_degrees[l] + 1); + mapping = std::make_unique>(FE_Q(1)); + } + + if (l == max_level) + mapping_ = mapping->clone(); + + // set up dofhandler + dof_handler.distribute_dofs(*fe); + + // set up constraints + const IndexSet locally_relevant_dofs = + DoFTools::extract_locally_relevant_dofs(dof_handler); + constraint.reinit(locally_relevant_dofs); + VectorTools::interpolate_boundary_values( + *mapping, dof_handler, 0, Functions::ZeroFunction(), constraint); + DoFTools::make_hanging_node_constraints(dof_handler, constraint); + constraint.close(); + + // set up operator + op.reinit(*mapping, dof_handler, *quad, constraint); + } + + // set up transfer operator + for (unsigned int l = min_level; l < max_level; ++l) + transfers[l + 1].reinit(dof_handlers[l + 1], + dof_handlers[l], + constraints[l + 1], + constraints[l]); + + MGTransferGlobalCoarsening transfer; + transfer.initialize_two_level_transfers(transfers); + transfer.build(fine_dof_handler, [&](const auto l, auto &vec) { + operators[l].initialize_dof_vector(vec); + }); + + GMGParameters mg_data; // TODO + + VectorType dst, src; + fine_operator.initialize_dof_vector(dst); + fine_operator.initialize_dof_vector(src); + fine_operator.rhs(src); + + ReductionControl solver_control( + mg_data.maxiter, mg_data.abstol, mg_data.reltol, false, false); + + mg_solve(solver_control, + dst, + src, + mg_data, + fine_dof_handler, + fine_operator, + operators, + transfer); + + fine_constraints.distribute(dst); + + deallog << dim << ' ' << fe_degree_fine << ' ' << n_refinements << ' ' + << (do_simplex_mesh ? "tri " : "quad") << ' ' + << solver_control.last_step() << std::endl; + + static unsigned int counter = 0; + + MGLevelObject results(min_level, max_level); + + transfer.interpolate_to_mg(dof_handlers[max_level], results, dst); + + for (unsigned int l = min_level; l <= max_level; ++l) + { + deallog << "Norm interpolated solution on level " << l << ": " + << results[l].l2_norm() << std::endl; + DataOut data_out; + + data_out.attach_dof_handler(dof_handlers[l]); + data_out.add_data_vector( + results[l], + "solution", + DataOut_DoFData::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 +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll all; + + deallog.precision(8); + + for (unsigned int mesh_type = 0; mesh_type < 2; ++mesh_type) + { + for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements) + for (unsigned int degree = 2; degree <= 4; ++degree) + test<2>(n_refinements, degree, false /*quadrilateral*/, mesh_type); + + for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements) + for (unsigned int degree = 2; degree <= 2; ++degree) + test<2>(n_refinements, degree, true /*triangle*/, mesh_type); + } +} diff --git a/tests/multigrid-global-coarsening/multigrid_p_04.mpirun=1.with_trilinos=true.output b/tests/multigrid-global-coarsening/multigrid_p_04.mpirun=1.with_trilinos=true.output new file mode 100644 index 0000000000..f2980eab6f --- /dev/null +++ b/tests/multigrid-global-coarsening/multigrid_p_04.mpirun=1.with_trilinos=true.output @@ -0,0 +1,79 @@ + +DEAL:0::2 2 2 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.33877875 +DEAL:0::Norm interpolated solution on level 1: 0.66014529 +DEAL:0::2 3 2 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.33872299 +DEAL:0::Norm interpolated solution on level 1: 0.99016130 +DEAL:0::2 4 2 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.33877890 +DEAL:0::Norm interpolated solution on level 1: 0.66010393 +DEAL:0::Norm interpolated solution on level 2: 1.3201957 +DEAL:0::2 2 3 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.66459543 +DEAL:0::Norm interpolated solution on level 1: 1.3203627 +DEAL:0::2 3 3 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.66458798 +DEAL:0::Norm interpolated solution on level 1: 1.9805368 +DEAL:0::2 4 3 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.66459725 +DEAL:0::Norm interpolated solution on level 1: 1.3203574 +DEAL:0::Norm interpolated solution on level 2: 2.6407133 +DEAL:0::2 2 4 quad 3 +DEAL:0::Norm interpolated solution on level 0: 1.3225826 +DEAL:0::Norm interpolated solution on level 1: 2.6407347 +DEAL:0::2 3 4 quad 3 +DEAL:0::Norm interpolated solution on level 0: 1.3225816 +DEAL:0::Norm interpolated solution on level 1: 3.9611011 +DEAL:0::2 4 4 quad 3 +DEAL:0::Norm interpolated solution on level 0: 1.3225830 +DEAL:0::Norm interpolated solution on level 1: 2.6407340 +DEAL:0::Norm interpolated solution on level 2: 5.2814679 +DEAL:0::2 2 2 tri 3 +DEAL:0::Norm interpolated solution on level 0: 0.33806168 +DEAL:0::Norm interpolated solution on level 1: 0.66011744 +DEAL:0::2 2 3 tri 3 +DEAL:0::Norm interpolated solution on level 0: 0.66445481 +DEAL:0::Norm interpolated solution on level 1: 1.3203595 +DEAL:0::2 2 4 tri 3 +DEAL:0::Norm interpolated solution on level 0: 1.3225585 +DEAL:0::Norm interpolated solution on level 1: 2.6407344 +DEAL:0::2 2 2 quad 2 +DEAL:0::Norm interpolated solution on level 0: 0.13157424 +DEAL:0::Norm interpolated solution on level 1: 0.24449632 +DEAL:0::2 3 2 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.13144458 +DEAL:0::Norm interpolated solution on level 1: 0.35587970 +DEAL:0::2 4 2 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.13139421 +DEAL:0::Norm interpolated solution on level 1: 0.24404614 +DEAL:0::Norm interpolated solution on level 2: 0.46558561 +DEAL:0::2 2 3 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.23622379 +DEAL:0::Norm interpolated solution on level 1: 0.45479680 +DEAL:0::2 3 3 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.23657515 +DEAL:0::Norm interpolated solution on level 1: 0.66773399 +DEAL:0::2 4 3 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.23652211 +DEAL:0::Norm interpolated solution on level 1: 0.45456995 +DEAL:0::Norm interpolated solution on level 2: 0.87733284 +DEAL:0::2 2 4 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.43119420 +DEAL:0::Norm interpolated solution on level 1: 0.85199621 +DEAL:0::2 3 4 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.43119824 +DEAL:0::Norm interpolated solution on level 1: 1.2615628 +DEAL:0::2 4 4 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.43118187 +DEAL:0::Norm interpolated solution on level 1: 0.85199072 +DEAL:0::Norm interpolated solution on level 2: 1.6669458 +DEAL:0::2 2 2 tri 3 +DEAL:0::Norm interpolated solution on level 0: 0.13575160 +DEAL:0::Norm interpolated solution on level 1: 0.24280231 +DEAL:0::2 2 3 tri 3 +DEAL:0::Norm interpolated solution on level 0: 0.21915908 +DEAL:0::Norm interpolated solution on level 1: 0.41924192 +DEAL:0::2 2 4 tri 4 +DEAL:0::Norm interpolated solution on level 0: 0.40718520 +DEAL:0::Norm interpolated solution on level 1: 0.80533381 -- 2.39.5