From eb669968ad502d6e39dd47b83f2d5cc85e6412bb Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 28 Mar 2023 16:05:15 +0200 Subject: [PATCH] Add benchmark for multgrid with global coarsening --- tests/performance/timing_mg_glob_coarsen.cc | 863 ++++++++++++++++++ ...ds=1.mpirun=max.exclusive.release.run_only | 0 tests/performance/timing_step_37.cc | 3 - 3 files changed, 863 insertions(+), 3 deletions(-) create mode 100644 tests/performance/timing_mg_glob_coarsen.cc create mode 100644 tests/performance/timing_mg_glob_coarsen.threads=1.mpirun=max.exclusive.release.run_only diff --git a/tests/performance/timing_mg_glob_coarsen.cc b/tests/performance/timing_mg_glob_coarsen.cc new file mode 100644 index 0000000000..f93a487f61 --- /dev/null +++ b/tests/performance/timing_mg_glob_coarsen.cc @@ -0,0 +1,863 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 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. +// +// --------------------------------------------------------------------- + +// +// Description: +// +// A performance benchmark assessing a Poisson problem with the +// performance-oriented matrix-free framework. As opposed to the related +// timing_step_37 benchmark, this case uses the global-coarsening multigrid +// framework with p-multigrid and using a locally refined mesh with hanging +// nodes. It also uses a setup with multiple DoFHandler objects, imitating the +// projection from a related (higher-order) DG function space. +// +// Status: experimental +// + +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +#define ENABLE_MPI + +#include "performance_test_driver.h" + +using namespace dealii; + +dealii::ConditionalOStream debug_output(std::cout, false); + + + +template +class LaplaceOperator : public Subscriptor +{ +public: + using value_type = number; + using VectorType = LinearAlgebra::distributed::Vector; + + LaplaceOperator(){}; + + void + initialize(const Mapping & mapping, + const DoFHandler & dof_handler, + const AffineConstraints &constraints, + const DoFHandler & dg_dof_handler) + { + const QGauss<1> quad(dof_handler.get_fe().degree + 1); + const QGauss<1> dg_quad(dg_dof_handler.get_fe().degree + 1); + typename MatrixFree::AdditionalData mf_data; + mf_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + mf_data.mapping_update_flags |= update_quadrature_points; + mf_data.mapping_update_flags_inner_faces = + (update_gradients | update_JxW_values); + mf_data.mapping_update_flags_boundary_faces = + (update_gradients | update_JxW_values); + AffineConstraints dg_constraints; + + data.reinit(mapping, + std::vector *>{ + {&dof_handler, &dg_dof_handler}}, + std::vector *>{ + {&constraints, &dg_constraints}}, + std::vector>{{quad, dg_quad}}, + mf_data); + } + + void + initialize(const Mapping & mapping, + const DoFHandler & dof_handler, + const AffineConstraints &constraints) + { + const QGauss<1> quad(dof_handler.get_fe().degree + 1); + typename MatrixFree::AdditionalData mf_data; + mf_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + Assert(dof_handler.get_fe().dofs_per_vertex > 0, + ExcNotImplemented("Only continuous elements implemented")); + + data.reinit(mapping, dof_handler, constraints, quad, mf_data); + } + + void + vmult(VectorType &dst, const VectorType &src) const + { + data.cell_loop(&LaplaceOperator::local_apply, this, dst, src, true); + for (const auto i : data.get_constrained_dofs()) + dst.local_element(i) = src.local_element(i); + } + + void + vmult(VectorType & dst, + const VectorType &src, + const std::function + &operation_before_loop, + const std::function + &operation_after_loop) const + { + data.cell_loop(&LaplaceOperator::local_apply, + this, + dst, + src, + operation_before_loop, + operation_after_loop); + } + + void + Tvmult(VectorType &dst, const VectorType &src) const + { + vmult(dst, src); + } + + number + el(const types::global_dof_index, const types::global_dof_index) const + { + AssertThrow(false, ExcNotImplemented()); + return number(0.); + } + + types::global_dof_index + m() const + { + return data.get_vector_partitioner()->size(); + } + + types::global_dof_index + n() const + { + return data.get_vector_partitioner()->size(); + } + + void + initialize_dof_vector(VectorType & vector, + const unsigned int component = 0) const + { + data.initialize_dof_vector(vector, component); + } + + void + compute_inverse_diagonal() + { + inverse_diagonal_entries = std::make_shared>(); + data.initialize_dof_vector(inverse_diagonal_entries->get_vector()); + MatrixFreeTools:: + compute_diagonal>( + data, inverse_diagonal_entries->get_vector(), [](auto &eval) { + eval.evaluate(EvaluationFlags::gradients); + for (unsigned int q = 0; q < eval.n_q_points; ++q) + eval.submit_gradient(eval.get_gradient(q), q); + eval.integrate(EvaluationFlags::gradients); + }); + + for (number &entry : inverse_diagonal_entries->get_vector()) + if (std::abs(entry) > 1e-10) + entry = 1. / entry; + else + entry = 1.; + } + + const std::shared_ptr> & + get_matrix_diagonal_inverse() const + { + return inverse_diagonal_entries; + } + + const MatrixFree & + get_matrix_free() const + { + return data; + } + +private: + void + local_apply(const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &cell_range) const + { + FEEvaluation eval(data); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + eval.reinit(cell); + eval.gather_evaluate(src, EvaluationFlags::gradients); + for (unsigned int q = 0; q < eval.n_q_points; ++q) + eval.submit_gradient(eval.get_gradient(q), q); + eval.integrate_scatter(EvaluationFlags::gradients, dst); + } + } + + MatrixFree data; + std::shared_ptr> inverse_diagonal_entries; +}; + + + +template +void +make_zero_mean(const std::vector & constrained_dofs, + LinearAlgebra::distributed::Vector &vec) +{ + // set constrained entries to zero + for (const unsigned int index : constrained_dofs) + vec.local_element(index) = 0.; + + // rescale mean value computed among all vector entries to the vector size + // without constraints + const unsigned int n_unconstrained_dofs = + vec.locally_owned_size() - constrained_dofs.size(); + vec.add( + -vec.mean_value() * vec.size() / + Utilities::MPI::sum(n_unconstrained_dofs, vec.get_mpi_communicator())); + + // set constrained entries to zero again, this should now have zero mean + for (const unsigned int index : constrained_dofs) + vec.local_element(index) = 0.; + + Assert(std::abs(vec.mean_value()) < + std::numeric_limits::epsilon() * vec.size(), + ExcInternalError()); +} + + + +// class to impose zero-mean constraint on coarse level +template > +class MGCoarseSolverSingular : public MGCoarseGridBase +{ +public: + void + clear() + { + coarse_smooth = nullptr; + } + + void + initialize(const MGSmootherBase &coarse_smooth, + const std::vector & constrained_dofs) + { + this->coarse_smooth = &coarse_smooth; + this->constrained_dofs = &constrained_dofs; + } + + void + operator()(const unsigned int level, + VectorType & dst, + const VectorType & src) const override + { + src_copy.reinit(src, true); + src_copy.copy_locally_owned_data_from(src); + make_zero_mean(*constrained_dofs, src_copy); + coarse_smooth->apply(level, dst, src_copy); + make_zero_mean(*constrained_dofs, dst); + } + +private: + SmartPointer> coarse_smooth; + const std::vector * constrained_dofs; + + mutable VectorType src_copy; +}; + + +const Tensor<2, 3> deformation{ + {{1.05, 1e-3, 1e-2}, {1e-3, 1., -1e-3}, {1e-2, -1e-3, 0.95}}}; + +template +class Solution : public Function +{ +public: + Solution() + : Function(1) + {} + + virtual double + value(const Point &p, const unsigned int = 0) const override + { + const Point x = Point(invert(deformation) * p); + double value = 1.0; + for (unsigned int d = 0; d < dim; ++d) + value *= std::cos(8. * numbers::PI * x[d]); + return value; + } +}; + + +template +class RightHandSide : public Function +{ +public: + RightHandSide() + : Function(1) + {} + + virtual double + value(const Point &p, const unsigned int = 0) const override + { + Solution sol; + return dim * 64. * numbers::PI * numbers::PI * sol.value(p); + } +}; + + + +template +class LaplaceProblem +{ +public: + LaplaceProblem(const unsigned int degree); + + Measurement + run(); + +private: + void + setup_grid(); + void + create_coarse_triangulations(); + void + setup_dofs(); + void + setup_matrix_free(); + void + setup_smoother(); + void + setup_transfer(); + void + compute_rhs(); + void + solve(); + void + embed_solution_to_dg(); + + parallel::distributed::Triangulation triangulation; + std::vector>> coarse_triangulations; + MappingQ mapping; + FE_DGQ dg_fe; + DoFHandler dg_dof_handler; + MGLevelObject>> fes; + MGLevelObject> dof_handlers; + + LinearAlgebra::distributed::Vector dg_rhs; + LinearAlgebra::distributed::Vector rhs; + LinearAlgebra::distributed::Vector solution; + LinearAlgebra::distributed::Vector dg_solution; + + LaplaceOperator system_matrix; + MGLevelObject> level_constraints; + MGLevelObject> level_matrices; + using VectorTypeMG = LinearAlgebra::distributed::Vector; + + using SmootherType = + PreconditionChebyshev, VectorTypeMG>; + mg::SmootherRelaxation mg_smoother; + + MGLevelObject> mg_transfers; + std::unique_ptr> mg_transfer; +}; + + + +template +LaplaceProblem::LaplaceProblem(const unsigned int degree) +#ifdef DEAL_II_WITH_P4EST + : triangulation(MPI_COMM_WORLD) +#else + : triangulation() +#endif + , mapping(1) + , dg_fe(degree) + , dg_dof_handler(triangulation) +{} + + + +template +void +LaplaceProblem::setup_grid() +{ + GridGenerator::hyper_cube(triangulation, 0., 1.); + GridTools::transform([](const Point &p) { return deformation * p; }, + triangulation); + + switch (get_testing_environment()) + { + case TestingEnvironment::light: + triangulation.refine_global(3); + break; + case TestingEnvironment::medium: + triangulation.refine_global(4); + break; + case TestingEnvironment::heavy: + triangulation.refine_global(5); + break; + } + + for (const auto &cell : triangulation.active_cell_iterators()) + if (cell->is_locally_owned() && cell->center().norm() < 1.1) + cell->set_refine_flag(); + triangulation.execute_coarsening_and_refinement(); + for (const auto &cell : triangulation.active_cell_iterators()) + if (cell->is_locally_owned() && + cell->center().distance(Point(0.3, 0.3, 0.3)) < 0.5) + cell->set_refine_flag(); + triangulation.execute_coarsening_and_refinement(); +} + + +template +void +LaplaceProblem::create_coarse_triangulations() +{ + coarse_triangulations = + MGTransferGlobalCoarseningTools::create_geometric_coarsening_sequence( + triangulation/*, + RepartitioningPolicyTools::MinimalGranularityPolicy(16)*/); +} + + +template +void +LaplaceProblem::setup_dofs() +{ + dg_dof_handler.reinit(triangulation); + dg_dof_handler.distribute_dofs(dg_fe); + + // the solver uses ph-multigrid according to + // https://doi.org/10.1016/j.jcp.2020.109538 and + // https://doi.org/10.1145/3580314 + + // start by creating levels of continuous elements + std::vector p_levels({dg_fe.degree - 1}); + while (p_levels.back() > 2) + p_levels.push_back(std::max(p_levels.back() - 2, 2u)); + fes.resize(0, p_levels.size() - 1); + for (unsigned int level = 0; level < p_levels.size(); ++level) + fes[level] = + std::make_unique>(p_levels[p_levels.size() - 1 - level]); + + dof_handlers.resize(0, coarse_triangulations.size() - 1 + fes.max_level()); + level_constraints.resize(0, dof_handlers.max_level()); + for (unsigned int level = dof_handlers.min_level(); + level <= dof_handlers.max_level(); + ++level) + { + DoFHandler &dof_h = dof_handlers[level]; + dof_h.reinit( + *coarse_triangulations[std::min(level, + triangulation.n_global_levels() - 1)]); + if (level < coarse_triangulations.size()) + dof_h.distribute_dofs(*fes[0]); + else + dof_h.distribute_dofs(*fes[level + 1 - coarse_triangulations.size()]); + + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_dofs(dof_h, relevant_dofs); + AffineConstraints &constraints = level_constraints[level]; + constraints.reinit(relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_h, constraints); + constraints.close(); + typename MatrixFree::AdditionalData additional_data; + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + + DoFRenumbering::matrix_free_data_locality(dof_h, + constraints, + additional_data); + + // now create the final constraints object + DoFTools::extract_locally_relevant_dofs(dof_h, relevant_dofs); + constraints.clear(); + constraints.reinit(relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_h, constraints); + constraints.close(); + } +} + + + +template +void +LaplaceProblem::setup_matrix_free() +{ + AffineConstraints constraints_fine; + constraints_fine.reinit(level_constraints.back().get_local_lines()); + constraints_fine.copy_from(level_constraints.back()); + system_matrix.initialize(mapping, + dof_handlers.back(), + constraints_fine, + dg_dof_handler); + system_matrix.initialize_dof_vector(dg_rhs, 1); + system_matrix.initialize_dof_vector(dg_solution, 1); + system_matrix.initialize_dof_vector(rhs, 0); + system_matrix.initialize_dof_vector(solution, 0); + + level_matrices.resize(0, dof_handlers.max_level()); + for (unsigned int level = dof_handlers.min_level(); + level <= dof_handlers.max_level(); + ++level) + { + level_matrices[level].initialize(mapping, + dof_handlers[level], + level_constraints[level]); + } +} + + +template +void +LaplaceProblem::setup_smoother() +{ + MGLevelObject smoother_data( + 0, dof_handlers.max_level()); + for (unsigned int level = dof_handlers.min_level(); + level <= dof_handlers.max_level(); + ++level) + { + level_matrices[level].compute_inverse_diagonal(); + + // manually compute the eigenvalue estimate for Chebyshev because we + // need to be careful with the constrained indices + IterationNumberControl control(12, 1e-6, false, false); + + using VectorType = LinearAlgebra::distributed::Vector; + SolverCG solver(control); + internal::PreconditionChebyshevImplementation::EigenvalueTracker + eigenvalue_tracker; + solver.connect_eigenvalues_slot( + [&eigenvalue_tracker](const std::vector &eigenvalues) { + eigenvalue_tracker.slot(eigenvalues); + }); + + VectorType sol, rhs; + level_matrices[level].initialize_dof_vector(sol); + level_matrices[level].initialize_dof_vector(rhs); + + for (float &a : rhs) + a = (double)rand() / RAND_MAX; + make_zero_mean( + level_matrices[level].get_matrix_free().get_constrained_dofs(), rhs); + solver.solve(level_matrices[level], + sol, + rhs, + *level_matrices[level].get_matrix_diagonal_inverse()); + + if (level > 0) + { + smoother_data[level].smoothing_range = 15.; + smoother_data[level].degree = 4; + } + else + { + // Coarse level: Use MG smoother as solver (should use p-multigrid + // or AMG for complicated meshes) + smoother_data[level].smoothing_range = + eigenvalue_tracker.values.back() / + eigenvalue_tracker.values.front(); + smoother_data[0].degree = numbers::invalid_unsigned_int; + } + smoother_data[level].max_eigenvalue = eigenvalue_tracker.values.back(); + smoother_data[level].eig_cg_n_iterations = 0; + smoother_data[level].preconditioner = + level_matrices[level].get_matrix_diagonal_inverse(); + } + + mg_smoother.initialize(level_matrices, smoother_data); +} + + +template +void +LaplaceProblem::setup_transfer() +{ + mg_transfers.resize(0, dof_handlers.max_level()); + for (unsigned int level = 1; level <= dof_handlers.max_level(); ++level) + { + mg_transfers[level].reinit(dof_handlers[level], + dof_handlers[level - 1], + level_constraints[level], + level_constraints[level - 1]); + } + + mg_transfer = std::make_unique>( + mg_transfers, [&](const unsigned level, VectorTypeMG &vec) { + level_matrices[level].initialize_dof_vector(vec); + }); +} + + +template +void +LaplaceProblem::compute_rhs() +{ + // interpolate to nodes + VectorTools::interpolate(mapping, + dg_dof_handler, + RightHandSide(), + dg_rhs); + + // do the interpolation 10 times to get better significance in the numbers + for (unsigned int i = 0; i < 10; ++i) + { + rhs = 0.; + FEEvaluation dg_eval( + system_matrix.get_matrix_free(), 1); + FEEvaluation eval(system_matrix.get_matrix_free(), + 0); + for (unsigned int cell = 0; + cell < system_matrix.get_matrix_free().n_cell_batches(); + ++cell) + { + eval.reinit(cell); + dg_eval.reinit(cell); + dg_eval.gather_evaluate(dg_rhs, EvaluationFlags::values); + for (unsigned int q = 0; q < eval.n_q_points; ++q) + eval.submit_value(dg_eval.get_value(q), q); + eval.integrate_scatter(EvaluationFlags::values, rhs); + } + rhs.compress(VectorOperation::add); + + // since we use Neumann boundary conditions on the whole boundary, the + // right hand side must have zero mean value to ensure a solvable system + make_zero_mean(system_matrix.get_matrix_free().get_constrained_dofs(), + rhs); + } +} + + + +template +void +LaplaceProblem::solve() +{ + MGCoarseSolverSingular mg_coarse; + mg_coarse.initialize( + mg_smoother, level_matrices[0].get_matrix_free().get_constrained_dofs()); + mg::Matrix mg_matrix(level_matrices); + + Multigrid mg( + mg_matrix, mg_coarse, *mg_transfer, mg_smoother, mg_smoother); + PreconditionMG> + preconditioner(dof_handlers.back(), mg, *mg_transfer); + + SolverControl control(20, 1e-10 * rhs.l2_norm()); + SolverCG> solver(control); + + solver.solve(system_matrix, solution, rhs, preconditioner); + AssertThrow(control.last_step() < 10, + ExcMessage("Solve should converge in at most 10 iterations")); +} + + + +template +void +LaplaceProblem::embed_solution_to_dg() +{ + make_zero_mean(system_matrix.get_matrix_free().get_constrained_dofs(), + solution); + + FEEvaluation dg_eval(system_matrix.get_matrix_free(), + 1); + MatrixFreeOperators::CellwiseInverseMassMatrix inverse_mass(dg_eval); + FEEvaluation eval(system_matrix.get_matrix_free(), 0); + // to get better timings, run the evaluation 10 times + for (unsigned int i = 0; i < 10; ++i) + { + solution.update_ghost_values(); + for (unsigned int cell = 0; + cell < system_matrix.get_matrix_free().n_cell_batches(); + ++cell) + { + eval.reinit(cell); + dg_eval.reinit(cell); + eval.gather_evaluate(solution, EvaluationFlags::values); + inverse_mass.transform_from_q_points_to_basis( + 1, eval.begin_values(), dg_eval.begin_dof_values()); + dg_eval.set_dof_values(dg_solution); + } + solution.zero_out_ghost_values(); + } + + // compute error + double error = 0; + for (unsigned int cell = 0; + cell < system_matrix.get_matrix_free().n_cell_batches(); + ++cell) + { + dg_eval.reinit(cell); + dg_eval.gather_evaluate(dg_solution, EvaluationFlags::values); + Solution solution; + double local_error = 0; + for (unsigned int q : dg_eval.quadrature_point_indices()) + for (unsigned int v = 0; + v < + system_matrix.get_matrix_free().n_active_entries_per_cell_batch( + cell); + ++v) + { + Point quadrature_point; + for (unsigned int d = 0; d < dim; ++d) + quadrature_point[d] = dg_eval.quadrature_point(q)[d][v]; + local_error += + Utilities::fixed_power<2>(solution.value(quadrature_point) - + dg_eval.get_value(q)[v]) * + dg_eval.JxW(q)[v]; + } + error += local_error; + } + error = + std::sqrt(Utilities::MPI::sum(error, dg_solution.get_mpi_communicator())); + // do to the deformed mesh, the chosen right hand side and solution match + // only approximately - we request a tolerance of 1e-2 + AssertThrow(error < 1e-2, ExcMessage("Error should be less than 1e-2")); +} + + +template +Measurement +LaplaceProblem::run() +{ + std::map timer; + + timer["setup_grid"].start(); + setup_grid(); + timer["setup_grid"].stop(); + + timer["setup_coarse_grids"].start(); + create_coarse_triangulations(); + timer["setup_coarse_grids"].stop(); + + timer["setup_dofs"].start(); + setup_dofs(); + timer["setup_dofs"].stop(); + + timer["setup_matrix_free"].start(); + setup_matrix_free(); + timer["setup_matrix_free"].stop(); + + timer["setup_smoother"].start(); + setup_smoother(); + timer["setup_smoother"].stop(); + + timer["setup_transfer"].start(); + setup_transfer(); + timer["setup_transfer"].stop(); + + timer["compute_rhs"].start(); + compute_rhs(); + timer["compute_rhs"].stop(); + + timer["solve"].start(); + solve(); + timer["solve"].stop(); + + const unsigned int n_repeat = 50; + timer["matvec_double"].start(); + for (unsigned int t = 0; t < n_repeat; ++t) + system_matrix.vmult(rhs, solution); + timer["matvec_double"].stop(); + + LinearAlgebra::distributed::Vector vec1, vec2; + level_matrices[level_matrices.max_level()].initialize_dof_vector(vec1); + vec2.reinit(vec1); + timer["matvec_float"].start(); + for (unsigned int t = 0; t < n_repeat; ++t) + level_matrices[level_matrices.max_level()].vmult(vec2, vec1); + timer["matvec_float"].stop(); + + timer["embed_dg_and_error"].start(); + embed_solution_to_dg(); + timer["embed_dg_and_error"].stop(); + + debug_output << std::endl; + return {timer["setup_grid"].wall_time(), + timer["setup_coarse_grids"].wall_time(), + timer["setup_dofs"].wall_time(), + timer["setup_matrix_free"].wall_time(), + timer["setup_smoother"].wall_time(), + timer["setup_transfer"].wall_time(), + timer["compute_rhs"].wall_time(), + timer["solve"].wall_time(), + timer["matvec_double"].wall_time(), + timer["matvec_float"].wall_time(), + timer["embed_dg_and_error"].wall_time()}; +} + + +std::tuple> +describe_measurements() +{ + return {Metric::timing, + 5, + {"setup_grid", + "setup_coarse_grids", + "setup_dofs", + "setup_matrix_free", + "setup_smoother", + "setup_transfer", + "compute_rhs", + "solve", + "matvec_double", + "matvec_float", + "embed_dg_and_error"}}; +} + + +Measurement +perform_single_measurement() +{ + // run in 3d with degree 5, i.e., degree 4 for the FEM part making the + // actual solve + return LaplaceProblem<3>(5).run(); +} diff --git a/tests/performance/timing_mg_glob_coarsen.threads=1.mpirun=max.exclusive.release.run_only b/tests/performance/timing_mg_glob_coarsen.threads=1.mpirun=max.exclusive.release.run_only new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/performance/timing_step_37.cc b/tests/performance/timing_step_37.cc index ef35faec65..314b8f10e9 100644 --- a/tests/performance/timing_step_37.cc +++ b/tests/performance/timing_step_37.cc @@ -27,7 +27,6 @@ #include #include #include -#include #include #include @@ -38,7 +37,6 @@ #include #include -#include #include #include @@ -60,7 +58,6 @@ #include #include -#include #include #include -- 2.39.5