--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/timer.h>
+
+#include <deal.II/distributed/repartitioning_policy_tools.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/operators.h>
+#include <deal.II/matrix_free/tools.h>
+
+#include <deal.II/multigrid/mg_coarse.h>
+#include <deal.II/multigrid/mg_matrix.h>
+#include <deal.II/multigrid/mg_smoother.h>
+#include <deal.II/multigrid/mg_transfer_global_coarsening.h>
+#include <deal.II/multigrid/multigrid.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+#include <memory>
+
+#define ENABLE_MPI
+
+#include "performance_test_driver.h"
+
+using namespace dealii;
+
+dealii::ConditionalOStream debug_output(std::cout, false);
+
+
+
+template <int dim, typename number = double>
+class LaplaceOperator : public Subscriptor
+{
+public:
+ using value_type = number;
+ using VectorType = LinearAlgebra::distributed::Vector<number>;
+
+ LaplaceOperator(){};
+
+ void
+ initialize(const Mapping<dim> & mapping,
+ const DoFHandler<dim> & dof_handler,
+ const AffineConstraints<number> &constraints,
+ const DoFHandler<dim> & 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<dim, number>::AdditionalData mf_data;
+ mf_data.tasks_parallel_scheme =
+ MatrixFree<dim, number>::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<number> dg_constraints;
+
+ data.reinit(mapping,
+ std::vector<const DoFHandler<dim> *>{
+ {&dof_handler, &dg_dof_handler}},
+ std::vector<const AffineConstraints<number> *>{
+ {&constraints, &dg_constraints}},
+ std::vector<Quadrature<1>>{{quad, dg_quad}},
+ mf_data);
+ }
+
+ void
+ initialize(const Mapping<dim> & mapping,
+ const DoFHandler<dim> & dof_handler,
+ const AffineConstraints<number> &constraints)
+ {
+ const QGauss<1> quad(dof_handler.get_fe().degree + 1);
+ typename MatrixFree<dim, number>::AdditionalData mf_data;
+ mf_data.tasks_parallel_scheme =
+ MatrixFree<dim, number>::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<void(const unsigned int, const unsigned int)>
+ &operation_before_loop,
+ const std::function<void(const unsigned int, const unsigned int)>
+ &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<DiagonalMatrix<VectorType>>();
+ data.initialize_dof_vector(inverse_diagonal_entries->get_vector());
+ MatrixFreeTools::
+ compute_diagonal<dim, -1, 0, 1, number, VectorizedArray<number>>(
+ 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<DiagonalMatrix<VectorType>> &
+ get_matrix_diagonal_inverse() const
+ {
+ return inverse_diagonal_entries;
+ }
+
+ const MatrixFree<dim, number> &
+ get_matrix_free() const
+ {
+ return data;
+ }
+
+private:
+ void
+ local_apply(const MatrixFree<dim, number> & data,
+ VectorType & dst,
+ const VectorType & src,
+ const std::pair<unsigned int, unsigned int> &cell_range) const
+ {
+ FEEvaluation<dim, -1, 0, 1, number> 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<dim, number> data;
+ std::shared_ptr<DiagonalMatrix<VectorType>> inverse_diagonal_entries;
+};
+
+
+
+template <typename Number>
+void
+make_zero_mean(const std::vector<unsigned int> & constrained_dofs,
+ LinearAlgebra::distributed::Vector<Number> &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<Number>::epsilon() * vec.size(),
+ ExcInternalError());
+}
+
+
+
+// class to impose zero-mean constraint on coarse level
+template <class VectorType = LinearAlgebra::distributed::Vector<double>>
+class MGCoarseSolverSingular : public MGCoarseGridBase<VectorType>
+{
+public:
+ void
+ clear()
+ {
+ coarse_smooth = nullptr;
+ }
+
+ void
+ initialize(const MGSmootherBase<VectorType> &coarse_smooth,
+ const std::vector<unsigned int> & 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<const MGSmootherBase<VectorType>> coarse_smooth;
+ const std::vector<unsigned int> * 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 <int dim>
+class Solution : public Function<dim>
+{
+public:
+ Solution()
+ : Function<dim>(1)
+ {}
+
+ virtual double
+ value(const Point<dim> &p, const unsigned int = 0) const override
+ {
+ const Point<dim> x = Point<dim>(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 <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+ RightHandSide()
+ : Function<dim>(1)
+ {}
+
+ virtual double
+ value(const Point<dim> &p, const unsigned int = 0) const override
+ {
+ Solution<dim> sol;
+ return dim * 64. * numbers::PI * numbers::PI * sol.value(p);
+ }
+};
+
+
+
+template <int dim>
+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<dim> triangulation;
+ std::vector<std::shared_ptr<const Triangulation<dim>>> coarse_triangulations;
+ MappingQ<dim> mapping;
+ FE_DGQ<dim> dg_fe;
+ DoFHandler<dim> dg_dof_handler;
+ MGLevelObject<std::unique_ptr<FE_Q<dim>>> fes;
+ MGLevelObject<DoFHandler<dim>> dof_handlers;
+
+ LinearAlgebra::distributed::Vector<double> dg_rhs;
+ LinearAlgebra::distributed::Vector<double> rhs;
+ LinearAlgebra::distributed::Vector<double> solution;
+ LinearAlgebra::distributed::Vector<double> dg_solution;
+
+ LaplaceOperator<dim, double> system_matrix;
+ MGLevelObject<AffineConstraints<float>> level_constraints;
+ MGLevelObject<LaplaceOperator<dim, float>> level_matrices;
+ using VectorTypeMG = LinearAlgebra::distributed::Vector<float>;
+
+ using SmootherType =
+ PreconditionChebyshev<LaplaceOperator<dim, float>, VectorTypeMG>;
+ mg::SmootherRelaxation<SmootherType, VectorTypeMG> mg_smoother;
+
+ MGLevelObject<MGTwoLevelTransfer<dim, VectorTypeMG>> mg_transfers;
+ std::unique_ptr<MGTransferGlobalCoarsening<dim, VectorTypeMG>> mg_transfer;
+};
+
+
+
+template <int dim>
+LaplaceProblem<dim>::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 <int dim>
+void
+LaplaceProblem<dim>::setup_grid()
+{
+ GridGenerator::hyper_cube(triangulation, 0., 1.);
+ GridTools::transform([](const Point<dim> &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<dim>(0.3, 0.3, 0.3)) < 0.5)
+ cell->set_refine_flag();
+ triangulation.execute_coarsening_and_refinement();
+}
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::create_coarse_triangulations()
+{
+ coarse_triangulations =
+ MGTransferGlobalCoarseningTools::create_geometric_coarsening_sequence(
+ triangulation/*,
+ RepartitioningPolicyTools::MinimalGranularityPolicy<dim>(16)*/);
+}
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::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<unsigned int> 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<FE_Q<dim>>(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<dim> &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<float> &constraints = level_constraints[level];
+ constraints.reinit(relevant_dofs);
+ DoFTools::make_hanging_node_constraints(dof_h, constraints);
+ constraints.close();
+ typename MatrixFree<dim, float>::AdditionalData additional_data;
+ additional_data.tasks_parallel_scheme =
+ MatrixFree<dim, float>::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 <int dim>
+void
+LaplaceProblem<dim>::setup_matrix_free()
+{
+ AffineConstraints<double> 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 <int dim>
+void
+LaplaceProblem<dim>::setup_smoother()
+{
+ MGLevelObject<typename SmootherType::AdditionalData> 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<float>;
+ SolverCG<VectorType> solver(control);
+ internal::PreconditionChebyshevImplementation::EigenvalueTracker
+ eigenvalue_tracker;
+ solver.connect_eigenvalues_slot(
+ [&eigenvalue_tracker](const std::vector<double> &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 <int dim>
+void
+LaplaceProblem<dim>::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<MGTransferGlobalCoarsening<dim, VectorTypeMG>>(
+ mg_transfers, [&](const unsigned level, VectorTypeMG &vec) {
+ level_matrices[level].initialize_dof_vector(vec);
+ });
+}
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::compute_rhs()
+{
+ // interpolate to nodes
+ VectorTools::interpolate(mapping,
+ dg_dof_handler,
+ RightHandSide<dim>(),
+ 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<dim, -1, 0, 1, double> dg_eval(
+ system_matrix.get_matrix_free(), 1);
+ FEEvaluation<dim, -1, 0, 1, double> 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 <int dim>
+void
+LaplaceProblem<dim>::solve()
+{
+ MGCoarseSolverSingular<VectorTypeMG> mg_coarse;
+ mg_coarse.initialize(
+ mg_smoother, level_matrices[0].get_matrix_free().get_constrained_dofs());
+ mg::Matrix<VectorTypeMG> mg_matrix(level_matrices);
+
+ Multigrid<VectorTypeMG> mg(
+ mg_matrix, mg_coarse, *mg_transfer, mg_smoother, mg_smoother);
+ PreconditionMG<dim,
+ VectorTypeMG,
+ MGTransferGlobalCoarsening<dim, VectorTypeMG>>
+ preconditioner(dof_handlers.back(), mg, *mg_transfer);
+
+ SolverControl control(20, 1e-10 * rhs.l2_norm());
+ SolverCG<LinearAlgebra::distributed::Vector<double>> solver(control);
+
+ solver.solve(system_matrix, solution, rhs, preconditioner);
+ AssertThrow(control.last_step() < 10,
+ ExcMessage("Solve should converge in at most 10 iterations"));
+}
+
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::embed_solution_to_dg()
+{
+ make_zero_mean(system_matrix.get_matrix_free().get_constrained_dofs(),
+ solution);
+
+ FEEvaluation<dim, -1, 0, 1, double> dg_eval(system_matrix.get_matrix_free(),
+ 1);
+ MatrixFreeOperators::CellwiseInverseMassMatrix<dim, -1> inverse_mass(dg_eval);
+ FEEvaluation<dim, -1, 0, 1, double> 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<dim> 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<dim> 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 <int dim>
+Measurement
+LaplaceProblem<dim>::run()
+{
+ std::map<std::string, dealii::Timer> 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<float> 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<Metric, unsigned int, std::vector<std::string>>
+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();
+}
projects / dealii.git / commitdiff