From d6b096bc248d2573d03544e3356a7f46cd4187d6 Mon Sep 17 00:00:00 2001 From: Marco Feder Date: Wed, 9 Oct 2024 16:47:12 +0200 Subject: [PATCH 1/2] Improve definition of agglomerates Update examples --- include/poly_utils.h | 89 +++++++++++++++++++ ...lly_distributed_poisson_sanity_check_01.cc | 20 ++--- ...lly_distributed_poisson_sanity_check_02.cc | 20 ++--- 3 files changed, 99 insertions(+), 30 deletions(-) diff --git a/include/poly_utils.h b/include/poly_utils.h index f14c5869..0c19d73b 100644 --- a/include/poly_utils.h +++ b/include/poly_utils.h @@ -616,6 +616,95 @@ namespace dealii::PolyUtils + /** + * Partition with METIS the locally owned regions of the given + * triangulation and insert agglomerates in the polytopic grid. + * + * @note The given triangulation must be a parallel::fullydistributed::Triangulation. This is + * required as the partitions generated by p4est, the partitioner for + * parallell::distributed::Triangulation, can generate discontinuous + * partitions which are not supported by the METIS partitioner. + * + */ + template + void + partition_locally_owned_regions( + AgglomerationHandler &agglomeration_handler, + const unsigned int n_partitions, + Triangulation &triangulation, + const SparsityTools::Partitioner partitioner) + { + AssertDimension(dim, spacedim); + Assert( + agglomeration_handler.n_agglomerates() == 0, + ExcMessage( + "The agglomerated grid must be empty upon calling this function.")); + Assert(n_partitions > 0, + ExcMessage("Invalid number of partitions, you provided " + + std::to_string(n_partitions))); + + auto parallel_triangulation = + dynamic_cast *>( + &triangulation); + Assert( + (parallel_triangulation != nullptr), + ExcMessage( + "Only fully distributed triangulations are supported. If you are using" + "a parallel::distributed::triangulation, you must convert it to a fully" + "distributed as explained in the documentation.")); + + // check for an easy return + if (n_partitions == 1) + { + for (const auto &cell : parallel_triangulation->active_cell_iterators()) + if (cell->is_locally_owned()) + agglomeration_handler.define_agglomerate({cell}); + return; + } + + // collect all locally owned cells + std::vector locally_owned_cells; + for (const auto &cell : triangulation.active_cell_iterators()) + if (cell->is_locally_owned()) + locally_owned_cells.push_back(cell->active_cell_index()); + + DynamicSparsityPattern cell_connectivity; + internal::get_face_connectivity_of_cells(*parallel_triangulation, + cell_connectivity, + locally_owned_cells); + + SparsityPattern sp_cell_connectivity; + sp_cell_connectivity.copy_from(cell_connectivity); + + // partition each locally owned connection graph and get + // back a vector of indices, one per degree + // of freedom (which is associated with a + // cell) + std::vector partition_indices( + parallel_triangulation->n_locally_owned_active_cells()); + SparsityTools::partition(sp_cell_connectivity, + n_partitions, + partition_indices, + partitioner); + + std::vector::active_cell_iterator>> + cells_per_partion_id; + cells_per_partion_id.resize(n_partitions); // number of agglomerates + + // finally loop over all cells and store the ones with same partition index + for (const auto &cell : parallel_triangulation->active_cell_iterators()) + if (cell->is_locally_owned()) + cells_per_partion_id[partition_indices[internal::get_index( + locally_owned_cells, cell->active_cell_index())]] + .push_back(cell); + + // All the cells with the same partition index will be merged together. + for (unsigned int i = 0; i < n_partitions; ++i) + agglomeration_handler.define_agglomerate(cells_per_partion_id[i]); + } + + + template std:: tuple, std::vector, std::vector, double> diff --git a/test/polydeal/fully_distributed_poisson_sanity_check_01.cc b/test/polydeal/fully_distributed_poisson_sanity_check_01.cc index f1d38683..580f2993 100644 --- a/test/polydeal/fully_distributed_poisson_sanity_check_01.cc +++ b/test/polydeal/fully_distributed_poisson_sanity_check_01.cc @@ -166,25 +166,15 @@ main(int argc, char *argv[]) const unsigned int n_local_agglomerates = 10; // number of agglomerates in each local subdomain - // Call the METIS partitioner to agglomerate within each processor. - PolyUtils::partition_locally_owned_regions(n_local_agglomerates, - tria_pft, - SparsityTools::Partitioner::metis); - // Setup the agglomeration handler. GridTools::Cache cached_tria(tria_pft); AgglomerationHandler ah(cached_tria); - // Agglomerate cells together based on their material id - std::vector::active_cell_iterator>> - cells_per_subdomain(n_local_agglomerates); - for (const auto &cell : tria_pft.active_cell_iterators()) - if (cell->is_locally_owned()) - cells_per_subdomain[cell->material_id()].push_back(cell); - - // Agglomerate elements with same id - for (std::size_t i = 0; i < cells_per_subdomain.size(); ++i) - ah.define_agglomerate(cells_per_subdomain[i]); + // Call the METIS partitioner to agglomerate within each processor. + PolyUtils::partition_locally_owned_regions(ah, + n_local_agglomerates, + tria_pft, + SparsityTools::Partitioner::metis); FE_DGQ<2> fe_dg(1); diff --git a/test/polydeal/fully_distributed_poisson_sanity_check_02.cc b/test/polydeal/fully_distributed_poisson_sanity_check_02.cc index 3561c789..2a017a30 100644 --- a/test/polydeal/fully_distributed_poisson_sanity_check_02.cc +++ b/test/polydeal/fully_distributed_poisson_sanity_check_02.cc @@ -148,10 +148,6 @@ main(int argc, char *argv[]) const unsigned int n_local_agglomerates = 10; // number of agglomerates in each local subdomain - // Call the METIS partitioner to agglomerate within each processor. - PolyUtils::partition_locally_owned_regions(n_local_agglomerates, - tria_pft, - SparsityTools::Partitioner::metis); // The rest of the program is identical to @@ -162,17 +158,11 @@ main(int argc, char *argv[]) GridTools::Cache cached_tria(tria_pft); AgglomerationHandler ah(cached_tria); - // Agglomerate cells together based on their material id - std::vector::active_cell_iterator>> - cells_per_subdomain(n_local_agglomerates); - for (const auto &cell : tria_pft.active_cell_iterators()) - if (cell->is_locally_owned()) - cells_per_subdomain[cell->material_id()].push_back(cell); - - // Agglomerate elements with same id - for (std::size_t i = 0; i < cells_per_subdomain.size(); ++i) - ah.define_agglomerate(cells_per_subdomain[i]); - + // Call the METIS partitioner to agglomerate within each processor. + PolyUtils::partition_locally_owned_regions(ah, + n_local_agglomerates, + tria_pft, + SparsityTools::Partitioner::metis); FE_DGQ<2> fe_dg(1); From 59cd0dff52607d959eb2a55768beef37acba65fd Mon Sep 17 00:00:00 2001 From: Marco Feder Date: Wed, 9 Oct 2024 18:18:33 +0200 Subject: [PATCH 2/2] Update examples --- examples/3D_piston.cc | 52 +++++----------------- examples/diffusion_reaction.cc | 55 +++++------------------ examples/poisson.cc | 81 ++++++++++++++++------------------ include/poly_utils.h | 4 +- 4 files changed, 62 insertions(+), 130 deletions(-) diff --git a/examples/3D_piston.cc b/examples/3D_piston.cc index 19288eb6..98b96c5f 100644 --- a/examples/3D_piston.cc +++ b/examples/3D_piston.cc @@ -423,18 +423,13 @@ DiffusionReactionProblem::setup_agglomerated_problem() ah = std::make_unique>(cached_tria); - // Agglomerate cells together based on their material id. - // The following vector stores cells with same material id, i.e.: - // v[material_id] = {cells with same id} - std::vector::active_cell_iterator>> - cells_per_material_id; - + double start, stop; + start = MPI_Wtime(); if (agglomeration_data.partitioning_strategy == PartitionerType::metis) { - cells_per_material_id.resize(agglomeration_data.agglomeration_parameter); - // Call the METIS partitioner to agglomerate within each processor. PolyUtils::partition_locally_owned_regions( + *ah, agglomeration_data.agglomeration_parameter, tria_pft, SparsityTools::Partitioner::metis); @@ -466,18 +461,8 @@ DiffusionReactionProblem::setup_agglomerated_problem() tree, agglomeration_data.agglomeration_parameter); const auto &agglomerates = csr_and_agglomerates.second; // vec> - std::size_t agglo_index = 0; - for (std::size_t i = 0; i < agglomerates.size(); ++i) - { - const auto &agglo = agglomerates[i]; - for (const auto &el : agglo) - el->set_material_id(agglo_index); - - ++agglo_index; // one agglomerate has been processed, increment - // counter - } - - cells_per_material_id.resize(agglo_index); + for (const auto &agglo : agglomerates) + ah->define_agglomerate(agglo); #ifdef AGGLO_DEBUG std::cout << "N agglomerates in process " @@ -492,22 +477,10 @@ DiffusionReactionProblem::setup_agglomerated_problem() // Irrespective of the partitioner strategy, define agglomerates based on the // material_id(s) attached to each cell. - - for (const auto &cell : tria_pft.active_cell_iterators()) - if (cell->is_locally_owned()) - cells_per_material_id[cell->material_id()].push_back(cell); - - // Agglomerate elements with same id - for (std::size_t i = 0; i < cells_per_material_id.size(); ++i) - ah->define_agglomerate(cells_per_material_id[i]); - - const unsigned int n_agglomerates = - agglomeration_data.partitioning_strategy == PartitionerType::rtree ? - Utilities::MPI::sum(cells_per_material_id.size(), comm) : - Utilities::MPI::n_mpi_processes(comm) * - agglomeration_data.agglomeration_parameter; - - pcout << "Total number of agglomerates: " << n_agglomerates << std::endl; + stop = MPI_Wtime(); + pcout << "Total number of agglomerates: " << ah->n_agglomerates() + << std::endl; + pcout << "Agglomeration performed in " << stop - start << "[s]" << std::endl; } @@ -922,7 +895,7 @@ main(int argc, char *argv[]) // Setup agglomeration data for rtree and METIS AgglomerationData rtree_data; rtree_data.partitioning_strategy = PartitionerType::rtree; - rtree_data.agglomeration_parameter = 2; // extraction level + rtree_data.agglomeration_parameter = 3; // extraction level AgglomerationData metis_data; metis_data.partitioning_strategy = PartitionerType::metis; @@ -931,9 +904,8 @@ main(int argc, char *argv[]) comm)); // number of local agglomerates const unsigned int reaction_coefficient = 0.; - for (const AgglomerationData &agglomeration_strategy : - {rtree_data, metis_data}) - for (unsigned int degree : {1, 2, 3, 4, 5, 6}) + for (const AgglomerationData &agglomeration_strategy : {metis_data}) + for (unsigned int degree : {1, 1, 1, 1}) { DiffusionReactionProblem problem(agglomeration_strategy, degree, diff --git a/examples/diffusion_reaction.cc b/examples/diffusion_reaction.cc index 6b8c9523..8bdb5a6a 100644 --- a/examples/diffusion_reaction.cc +++ b/examples/diffusion_reaction.cc @@ -17,8 +17,6 @@ // where Omega = [0,1]^3 // -#define HEX - #include #include #include @@ -26,8 +24,6 @@ #include #include -#include - #include #include @@ -311,16 +307,9 @@ class DiffusionReactionProblem void output_results() const; - const MPI_Comm comm; -#ifdef HEX - MappingQ1 mapping; - FE_DGQ fe_dg; -#else - MappingFE mapping; - FE_SimplexDGP fe_dg; -#endif - const unsigned int n_ranks; - + const MPI_Comm comm; + const unsigned int n_ranks; + FE_DGQ fe_dg; const unsigned int n_local_agglomerates; const double reaction_coefficient; parallel::fullydistributed::Triangulation tria_pft; @@ -351,13 +340,8 @@ DiffusionReactionProblem::DiffusionReactionProblem( const double reaction_coefficient, const MPI_Comm communicator) : comm(communicator) -#ifdef HEX - , mapping() -#else - , mapping(FE_SimplexP{1}) -#endif - , fe_dg(degree) , n_ranks(Utilities::MPI::n_mpi_processes(comm)) + , fe_dg(degree) , n_local_agglomerates(n_local_agglomerates) , reaction_coefficient(reaction_coefficient) , tria_pft(comm) @@ -375,15 +359,8 @@ DiffusionReactionProblem::make_fine_grid( const unsigned int n_global_refinements) { Triangulation tria; -#ifdef HEX GridGenerator::hyper_cube(tria); tria.refine_global(n_global_refinements); -#else - Triangulation tria_hex; - GridGenerator::hyper_cube(tria_hex, 0., 1.); - tria_hex.refine_global(n_global_refinements - 1); - GridGenerator::convert_hypercube_to_simplex_mesh(tria_hex, tria); -#endif // Partition serial triangulation: GridTools::partition_triangulation(n_ranks, tria); @@ -437,21 +414,13 @@ DiffusionReactionProblem::assemble_system() { constraints.close(); - const unsigned int quadrature_degree = fe_dg.get_degree() + 1; - const unsigned int face_quadrature_degree = fe_dg.get_degree() + 1; -#ifdef HEX + const unsigned int quadrature_degree = 2 * fe_dg.get_degree() + 1; + const unsigned int face_quadrature_degree = 2 * fe_dg.get_degree() + 1; ah->initialize_fe_values(QGauss(quadrature_degree), - update_gradients | update_JxW_values | - update_quadrature_points | update_JxW_values | - update_values, - QGauss(face_quadrature_degree)); -#else - ah->initialize_fe_values(QGaussSimplex(quadrature_degree), - update_gradients | update_JxW_values | - update_quadrature_points | update_JxW_values | - update_values, - QGaussSimplex(face_quadrature_degree)); -#endif + update_values | update_gradients | + update_JxW_values | update_quadrature_points, + QGauss(face_quadrature_degree), + update_JxW_values); ah->distribute_agglomerated_dofs(fe_dg); @@ -865,11 +834,7 @@ main(int argc, char *argv[]) // number of agglomerates in each local subdomain const unsigned int n_local_agglomerates = 10; const double reaction_coefficient = .5; -#ifdef HEX for (unsigned int degree : {1, 2, 3, 4}) -#else - for (unsigned int degree : {1, 2, 3}) // degree > 4 not implemented -#endif { DiffusionReactionProblem problem(n_local_agglomerates, degree, diff --git a/examples/poisson.cc b/examples/poisson.cc index f39f3da1..fe60311e 100644 --- a/examples/poisson.cc +++ b/examples/poisson.cc @@ -486,6 +486,7 @@ Poisson::Poisson(const GridType &grid_type, analytical_solution = std::make_unique>(); rhs_function = std::make_unique>(solution_type); + constraints.close(); } template @@ -546,11 +547,22 @@ Poisson::make_grid() GridTools::partition_triangulation(n_subdomains, tria, SparsityTools::Partitioner::metis); + + std::vector< + std::vector::active_cell_iterator>> + cells_per_subdomain(n_subdomains); + for (const auto &cell : tria.active_cell_iterators()) + cells_per_subdomain[cell->subdomain_id()].push_back(cell); + + // For every subdomain, agglomerate elements together + for (std::size_t i = 0; i < n_subdomains; ++i) + ah->define_agglomerate(cells_per_subdomain[i]); + + std::chrono::duration wctduration = (std::chrono::system_clock::now() - start); std::cout << "METIS built in " << wctduration.count() << " seconds [Wall Clock]" << std::endl; - std::cout << "N subdomains: " << n_subdomains << std::endl; } else if (partitioner_type == PartitionerType::rtree) { @@ -580,28 +592,17 @@ Poisson::make_grid() CellsAgglomerator agglomerator{tree, extraction_level}; - const auto agglomerates = agglomerator.extract_agglomerates(); + const auto vec_agglomerates = agglomerator.extract_agglomerates(); // ah->connect_hierarchy(agglomerator); // Flag elements for agglomeration - unsigned int agglo_index = 0; - for (unsigned int i = 0; i < agglomerates.size(); ++i) - { - const auto &agglo = agglomerates[i]; // i-th agglomerate - for (const auto &el : agglo) - { - el->set_subdomain_id(agglo_index); - } - ++agglo_index; - } + for (const auto &agglo : vec_agglomerates) + ah->define_agglomerate(agglo); std::chrono::duration wctduration = (std::chrono::system_clock::now() - start); std::cout << "R-tree agglomerates built in " << wctduration.count() << " seconds [Wall Clock]" << std::endl; - n_subdomains = agglo_index; - - std::cout << "N subdomains = " << n_subdomains << std::endl; // Check number of agglomerates if constexpr (dim == 2) @@ -634,34 +635,19 @@ Poisson::make_grid() { Assert(false, ExcMessage("Wrong partitioning.")); } - - - constraints.close(); + n_subdomains = ah->n_agglomerates(); + std::cout << "N subdomains = " << n_subdomains << std::endl; } template void Poisson::setup_agglomeration() { - if (partitioner_type != PartitionerType::no_partition) - { - std::vector< - std::vector::active_cell_iterator>> - cells_per_subdomain(n_subdomains); - for (const auto &cell : tria.active_cell_iterators()) - cells_per_subdomain[cell->subdomain_id()].push_back(cell); - - // For every subdomain, agglomerate elements together - for (std::size_t i = 0; i < cells_per_subdomain.size(); ++i) - ah->define_agglomerate(cells_per_subdomain[i]); - } - else + if (partitioner_type == PartitionerType::no_partition) { // No partitioning means that each cell is a master cell for (const auto &cell : tria.active_cell_iterators()) - { - ah->define_agglomerate({cell}); - } + ah->define_agglomerate({cell}); } ah->distribute_agglomerated_dofs(dg_fe); @@ -1036,23 +1022,32 @@ Poisson::output_results() PolyUtils::interpolate_to_fine_grid(*ah, interpolated_solution, solution, - true /*no_the_fly*/); + true /*on_the_fly*/); data_out.attach_dof_handler(ah->output_dh); data_out.add_data_vector(interpolated_solution, "u", DataOut::type_dof_data); - Vector agglomerated(tria.n_active_cells()); Vector agglo_idx(tria.n_active_cells()); - for (const auto &cell : tria.active_cell_iterators()) + + // Mark fine cells belonging to the same agglomerate. + for (const auto &polytope : ah->polytope_iterators()) { - agglomerated[cell->active_cell_index()] = - ah->get_relationships()[cell->active_cell_index()]; - agglo_idx[cell->active_cell_index()] = cell->subdomain_id(); + const types::global_cell_index polytope_index = polytope->index(); + const auto &patch_of_cells = polytope->get_agglomerate(); // fine cells + // Flag them + for (const auto &cell : patch_of_cells) + agglo_idx[cell->active_cell_index()] = polytope_index; } - data_out.add_data_vector(agglomerated, - "agglo_relationships", - DataOut::type_cell_data); + + // Old way, here just for completeness + // for (const auto &cell : tria.active_cell_iterators()) + // { + // agglomerated[cell->active_cell_index()] = + // ah->get_relationships()[cell->active_cell_index()]; + // agglo_idx[cell->active_cell_index()] = cell->subdomain_id(); + // } + data_out.add_data_vector(agglo_idx, "agglo_idx", DataOut::type_cell_data); diff --git a/include/poly_utils.h b/include/poly_utils.h index 0c19d73b..e2956ca2 100644 --- a/include/poly_utils.h +++ b/include/poly_utils.h @@ -650,8 +650,8 @@ namespace dealii::PolyUtils (parallel_triangulation != nullptr), ExcMessage( "Only fully distributed triangulations are supported. If you are using" - "a parallel::distributed::triangulation, you must convert it to a fully" - "distributed as explained in the documentation.")); + "a parallel::distributed::triangulation, you must convert it to a" + "fully distributed as explained in the documentation.")); // check for an easy return if (n_partitions == 1)