Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improve definition of agglomerates #141

Merged
merged 2 commits into from
Oct 9, 2024
Merged

Improve definition of agglomerates #141

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d6b096b
Improve definition of agglomerates
fdrmrc Oct 9, 2024
59cd0df
Update examples
fdrmrc Oct 9, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
52 changes: 12 additions & 40 deletions examples/3D_piston.cc
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading