Skip to content

Commit

Permalink
Improve definition of agglomerates (#141)
Browse files Browse the repository at this point in the history 
* Improve definition of agglomerates

Update examples

* Update examples
  • Loading branch information
@fdrmrc
fdrmrc authored Oct 9, 2024
1 parent da242ce commit c4e67e5
Show file tree
Hide file tree
Showing 6 changed files with 159 additions and 158 deletions.
52 changes: 12 additions & 40 deletions examples/3D_piston.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -423,18 +423,13 @@ DiffusionReactionProblem<dim>::setup_agglomerated_problem()

ah = std::make_unique<AgglomerationHandler<dim>>(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<std::vector<typename Triangulation<dim>::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);
Expand Down Expand Up @@ -466,18 +461,8 @@ DiffusionReactionProblem<dim>::setup_agglomerated_problem()
tree, agglomeration_data.agglomeration_parameter);
const auto &agglomerates = csr_and_agglomerates.second; // vec<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 "
Expand All @@ -492,22 +477,10 @@ DiffusionReactionProblem<dim>::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;
}


Expand Down Expand Up @@ -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;
Expand All @@ -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<dim> problem(agglomeration_strategy,
degree,
Expand Down
55 changes: 10 additions & 45 deletions examples/diffusion_reaction.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,17 +17,13 @@
// where Omega = [0,1]^3
//

#define HEX

#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/utilities.h>

#include <deal.II/distributed/fully_distributed_tria.h>
#include <deal.II/distributed/tria.h>

#include <deal.II/fe/mapping_fe.h>

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>

Expand Down Expand Up @@ -311,16 +307,9 @@ class DiffusionReactionProblem
void
output_results() const;

const MPI_Comm comm;
#ifdef HEX
MappingQ1<dim> mapping;
FE_DGQ<dim> fe_dg;
#else
MappingFE<dim> mapping;
FE_SimplexDGP<dim> fe_dg;
#endif
const unsigned int n_ranks;

const MPI_Comm comm;
const unsigned int n_ranks;
FE_DGQ<dim> fe_dg;
const unsigned int n_local_agglomerates;
const double reaction_coefficient;
parallel::fullydistributed::Triangulation<dim> tria_pft;
Expand Down Expand Up @@ -351,13 +340,8 @@ DiffusionReactionProblem<dim>::DiffusionReactionProblem(
const double reaction_coefficient,
const MPI_Comm communicator)
: comm(communicator)
#ifdef HEX
, mapping()
#else
, mapping(FE_SimplexP<dim>{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)
Expand All @@ -375,15 +359,8 @@ DiffusionReactionProblem<dim>::make_fine_grid(
const unsigned int n_global_refinements)
{
Triangulation<dim> tria;
#ifdef HEX
GridGenerator::hyper_cube(tria);
tria.refine_global(n_global_refinements);
#else
Triangulation<dim> 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);
Expand Down Expand Up @@ -437,21 +414,13 @@ DiffusionReactionProblem<dim>::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<dim>(quadrature_degree),
update_gradients | update_JxW_values |
update_quadrature_points | update_JxW_values |
update_values,
QGauss<dim - 1>(face_quadrature_degree));
#else
ah->initialize_fe_values(QGaussSimplex<dim>(quadrature_degree),
update_gradients | update_JxW_values |
update_quadrature_points | update_JxW_values |
update_values,
QGaussSimplex<dim - 1>(face_quadrature_degree));
#endif
update_values | update_gradients |
update_JxW_values | update_quadrature_points,
QGauss<dim - 1>(face_quadrature_degree),
update_JxW_values);

ah->distribute_agglomerated_dofs(fe_dg);

Expand Down Expand Up @@ -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<dim> problem(n_local_agglomerates,
degree,
Expand Down
81 changes: 38 additions & 43 deletions examples/poisson.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -486,6 +486,7 @@ Poisson<dim>::Poisson(const GridType &grid_type,
analytical_solution = std::make_unique<SolutionProductSine<dim>>();

rhs_function = std::make_unique<const RightHandSide<dim>>(solution_type);
constraints.close();
}

template <int dim>
Expand Down Expand Up @@ -546,11 +547,22 @@ Poisson<dim>::make_grid()
GridTools::partition_triangulation(n_subdomains,
tria,
SparsityTools::Partitioner::metis);

std::vector<
std::vector<typename Triangulation<dim>::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<double> 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)
{
Expand Down Expand Up @@ -580,28 +592,17 @@ Poisson<dim>::make_grid()

CellsAgglomerator<dim, decltype(tree)> 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<double> 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)
Expand Down Expand Up @@ -634,34 +635,19 @@ Poisson<dim>::make_grid()
{
Assert(false, ExcMessage("Wrong partitioning."));
}


constraints.close();
n_subdomains = ah->n_agglomerates();
std::cout << "N subdomains = " << n_subdomains << std::endl;
}

template <int dim>
void
Poisson<dim>::setup_agglomeration()
{
if (partitioner_type != PartitionerType::no_partition)
{
std::vector<
std::vector<typename Triangulation<dim>::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);
Expand Down Expand Up @@ -1036,23 +1022,32 @@ Poisson<dim>::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<dim>::type_dof_data);

Vector<float> agglomerated(tria.n_active_cells());
Vector<float> 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<dim>::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<dim>::type_cell_data);
Expand Down
Loading

0 comments on commit c4e67e5

Please sign in to comment.