From fe088b7819cc69f8c08bcce710657478d27e41b0 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 12 Oct 2021 15:09:53 +0200 Subject: [PATCH] New test case --- .../parallel_multigrid_interleave.cc | 473 ++++++++++++++++++ ...4est=true.with_lapack=true.mpirun=2.output | 65 +++ 2 files changed, 538 insertions(+) create mode 100644 tests/matrix_free/parallel_multigrid_interleave.cc create mode 100644 tests/matrix_free/parallel_multigrid_interleave.with_p4est=true.with_lapack=true.mpirun=2.output diff --git a/tests/matrix_free/parallel_multigrid_interleave.cc b/tests/matrix_free/parallel_multigrid_interleave.cc new file mode 100644 index 0000000000..20b59fc2bc --- /dev/null +++ b/tests/matrix_free/parallel_multigrid_interleave.cc @@ -0,0 +1,473 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2014 - 2021 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. +// +// --------------------------------------------------------------------- + + +// Similar test as parallel_multigrid_mf, but using the functionality of +// PreconditionChebyshev to embed vector updates into the matrix-free loops of +// a suitable operator class. + +#include + +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include "../tests.h" + + +template +class LaplaceOperator : public Subscriptor +{ +public: + using value_type = number; + + LaplaceOperator() + : n_calls_vmult(0) + {} + + ~LaplaceOperator() + { + if (n_calls_vmult > 0) + deallog << "Number of calls to special vmult for Operator of size " << m() + << ": " << n_calls_vmult << std::endl; + } + + void + initialize(const Mapping & mapping, + const DoFHandler & dof_handler, + const std::set &dirichlet_boundaries, + const unsigned int level = numbers::invalid_unsigned_int) + { + n_calls_vmult = 0; + const QGauss<1> quad(dof_handler.get_fe().degree + 1); + typename MatrixFree::AdditionalData addit_data; + addit_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + addit_data.mg_level = level; + + // extract the constraints due to Dirichlet boundary conditions + AffineConstraints constraints; + Functions::ZeroFunction zero; + std::map *> functions; + for (std::set::const_iterator it = + dirichlet_boundaries.begin(); + it != dirichlet_boundaries.end(); + ++it) + functions[*it] = &zero; + if (level == numbers::invalid_unsigned_int) + VectorTools::interpolate_boundary_values(dof_handler, + functions, + constraints); + else + { + std::vector local_dofs; + typename DoFHandler::cell_iterator cell = dof_handler.begin(level), + endc = dof_handler.end(level); + for (; cell != endc; ++cell) + { + if (dof_handler.get_triangulation().locally_owned_subdomain() != + numbers::invalid_subdomain_id && + cell->level_subdomain_id() == numbers::artificial_subdomain_id) + continue; + const FiniteElement &fe = cell->get_fe(); + local_dofs.resize(fe.dofs_per_face); + + for (const unsigned int face_no : GeometryInfo::face_indices()) + if (cell->at_boundary(face_no) == true) + { + const typename DoFHandler::face_iterator face = + cell->face(face_no); + const types::boundary_id bi = face->boundary_id(); + if (functions.find(bi) != functions.end()) + { + face->get_mg_dof_indices(level, local_dofs); + for (unsigned int i = 0; i < fe.dofs_per_face; ++i) + constraints.add_line(local_dofs[i]); + } + } + } + } + constraints.close(); + + data.reinit(mapping, dof_handler, constraints, quad, addit_data); + + compute_inverse_diagonal(); + } + + void + vmult(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const + { + data.cell_loop(&LaplaceOperator::local_apply, this, dst, src, true); + for (unsigned int i : data.get_constrained_dofs()) + dst.local_element(i) = src.local_element(i); + } + + void + vmult(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::function + &operation_before_loop, + const std::function + &operation_after_loop) const + { + ++n_calls_vmult; + data.cell_loop(&LaplaceOperator::local_apply, + this, + dst, + src, + operation_before_loop, + operation_after_loop); + for (unsigned int i : data.get_constrained_dofs()) + dst.local_element(i) = src.local_element(i); + } + + void + Tvmult(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const + { + vmult(dst, src); + } + + void + vmult_add(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const + { + data.cell_loop(&LaplaceOperator::local_apply, this, dst, src); + for (unsigned int i : data.get_constrained_dofs()) + dst.local_element(i) += src.local_element(i); + } + + void + Tvmult_add(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const + { + vmult_add(dst, src); + } + + types::global_dof_index + m() const + { + return data.get_vector_partitioner()->size(); + } + + types::global_dof_index + n() const + { + return data.get_vector_partitioner()->size(); + } + + number + el(const unsigned int row, const unsigned int col) const + { + AssertThrow(false, + ExcMessage("Matrix-free does not allow for entry access")); + return number(); + } + + void + initialize_dof_vector( + LinearAlgebra::distributed::Vector &vector) const + { + if (!vector.partitioners_are_compatible( + *data.get_dof_info(0).vector_partitioner)) + data.initialize_dof_vector(vector); + Assert(vector.partitioners_are_globally_compatible( + *data.get_dof_info(0).vector_partitioner), + ExcInternalError()); + } + + const LinearAlgebra::distributed::Vector & + get_matrix_diagonal_inverse() const + { + Assert(inverse_diagonal_entries.size() > 0, ExcNotInitialized()); + return inverse_diagonal_entries; + } + + +private: + void + local_apply(const MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair &cell_range) const + { + FEEvaluation phi(data); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + phi.reinit(cell); + phi.read_dof_values(src); + phi.evaluate(EvaluationFlags::gradients); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(phi.get_gradient(q), q); + phi.integrate(EvaluationFlags::gradients); + phi.distribute_local_to_global(dst); + } + } + + void + compute_inverse_diagonal() + { + data.initialize_dof_vector(inverse_diagonal_entries); + unsigned int dummy = 0; + data.cell_loop(&LaplaceOperator::local_diagonal_cell, + this, + inverse_diagonal_entries, + dummy); + + for (unsigned int i = 0; i < inverse_diagonal_entries.local_size(); ++i) + if (std::abs(inverse_diagonal_entries.local_element(i)) > 1e-10) + inverse_diagonal_entries.local_element(i) = + 1. / inverse_diagonal_entries.local_element(i); + else + inverse_diagonal_entries.local_element(i) = 1.; + } + + void + local_diagonal_cell( + const MatrixFree & data, + LinearAlgebra::distributed::Vector &dst, + const unsigned int &, + const std::pair &cell_range) const + { + FEEvaluation phi(data); + + AlignedVector> local_diagonal_vector( + phi.dofs_per_cell); + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + phi.reinit(cell); + + for (unsigned int i = 0; i < phi.dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < phi.dofs_per_cell; ++j) + phi.begin_dof_values()[j] = VectorizedArray(); + phi.begin_dof_values()[i] = 1.; + phi.evaluate(EvaluationFlags::gradients); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(phi.get_gradient(q), q); + phi.integrate(EvaluationFlags::gradients); + local_diagonal_vector[i] = phi.begin_dof_values()[i]; + } + for (unsigned int i = 0; i < phi.dofs_per_cell; ++i) + phi.begin_dof_values()[i] = local_diagonal_vector[i]; + phi.distribute_local_to_global(dst); + } + } + + MatrixFree data; + LinearAlgebra::distributed::Vector inverse_diagonal_entries; + mutable unsigned int n_calls_vmult; +}; + + + +template +class MGTransferPrebuiltMF + : public MGTransferMatrixFree +{ +public: + MGTransferPrebuiltMF(const MGLevelObject &laplace, + const MGConstrainedDoFs & mg_constrained_dofs) + : MGTransferMatrixFree( + mg_constrained_dofs) + , laplace_operator(laplace){}; + + /** + * Overload copy_to_mg from MGTransferPrebuilt to get the vectors compatible + * with MatrixFree and bypass the crude initialization in MGTransferPrebuilt + */ + template + void + copy_to_mg( + const DoFHandler &mg_dof, + MGLevelObject< + LinearAlgebra::distributed::Vector> &dst, + const InVector &src) const + { + for (unsigned int level = dst.min_level(); level <= dst.max_level(); + ++level) + laplace_operator[level].initialize_dof_vector(dst[level]); + MGLevelGlobalTransfer>::copy_to_mg(mg_dof, dst, src); + } + +private: + const MGLevelObject &laplace_operator; +}; + + + +template +void +do_test(const DoFHandler &dof) +{ + deallog << "Testing " << dof.get_fe().get_name(); + deallog << std::endl; + deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl; + + const unsigned int fe_degree = dof.get_fe().degree; + MappingQ mapping(fe_degree + 1); + LaplaceOperator fine_matrix; + std::set dirichlet_boundaries; + dirichlet_boundaries.insert(0); + fine_matrix.initialize(mapping, dof, dirichlet_boundaries); + + LinearAlgebra::distributed::Vector in, sol; + fine_matrix.initialize_dof_vector(in); + fine_matrix.initialize_dof_vector(sol); + + // set constant rhs vector, except boundary + in = 1.; + { + std::map boundary_values; + VectorTools::interpolate_boundary_values(dof, + 0, + Functions::ZeroFunction(), + boundary_values); + for (const auto it : boundary_values) + if (dof.locally_owned_dofs().is_element(it.first)) + in(it.first) = 0.; + } + + // set up multigrid in analogy to step-37 + using LevelMatrixType = LaplaceOperator; + + MGLevelObject mg_matrices; + mg_matrices.resize(0, dof.get_triangulation().n_global_levels() - 1); + for (unsigned int level = 0; + level < dof.get_triangulation().n_global_levels(); + ++level) + { + mg_matrices[level].initialize(mapping, dof, dirichlet_boundaries, level); + } + + MGConstrainedDoFs mg_constrained_dofs; + mg_constrained_dofs.initialize(dof); + mg_constrained_dofs.make_zero_boundary_constraints(dof, {0}); + + MGTransferPrebuiltMF mg_transfer(mg_matrices, + mg_constrained_dofs); + mg_transfer.build(dof); + + using SMOOTHER = + PreconditionChebyshev>; + mg::SmootherRelaxation> + mg_smoother; + + MGLevelObject smoother_data; + smoother_data.resize(0, dof.get_triangulation().n_global_levels() - 1); + for (unsigned int level = 0; + level < dof.get_triangulation().n_global_levels(); + ++level) + { + smoother_data[level].smoothing_range = 15.; + smoother_data[level].degree = 5; + smoother_data[level].eig_cg_n_iterations = 15; + auto preconditioner = std::make_shared< + DiagonalMatrix>>(); + preconditioner->reinit(mg_matrices[level].get_matrix_diagonal_inverse()); + smoother_data[level].preconditioner = std::move(preconditioner); + } + + mg_smoother.initialize(mg_matrices, smoother_data); + MGCoarseGridApplySmoother> + mg_coarse(mg_smoother); + + mg::Matrix> mg_matrix(mg_matrices); + + Multigrid> mg( + mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother); + PreconditionMG, + MGTransferPrebuiltMF> + preconditioner(dof, mg, mg_transfer); + + { + // avoid output from inner (coarse-level) solver + deallog.depth_file(2); + ReductionControl control(30, 1e-20, 1e-7); + SolverCG> solver(control); + solver.solve(fine_matrix, sol, in, preconditioner); + } +} + + + +template +void +test(const unsigned int fe_degree) +{ + for (unsigned int i = 5; i < 7; ++i) + { + parallel::distributed::Triangulation tria( + MPI_COMM_WORLD, + Triangulation::limit_level_difference_at_vertices, + parallel::distributed::Triangulation< + dim>::construct_multigrid_hierarchy); + GridGenerator::hyper_cube(tria); + tria.refine_global(i - dim); + + FE_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + dof.distribute_mg_dofs(); + + do_test(dof); + } +} + + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + + mpi_initlog(); + + { + test<2, double>(1); + test<2, float>(3); + + test<3, double>(1); + test<3, float>(2); + } +} diff --git a/tests/matrix_free/parallel_multigrid_interleave.with_p4est=true.with_lapack=true.mpirun=2.output b/tests/matrix_free/parallel_multigrid_interleave.with_p4est=true.with_lapack=true.mpirun=2.output new file mode 100644 index 0000000000..34c42cc305 --- /dev/null +++ b/tests/matrix_free/parallel_multigrid_interleave.with_p4est=true.with_lapack=true.mpirun=2.output @@ -0,0 +1,65 @@ + +DEAL::Testing FE_Q<2>(1) +DEAL::Number of degrees of freedom: 81 +DEAL:cg::Starting value 7.00000 +DEAL:cg::Convergence step 3 value 1.44113e-07 +DEAL::Number of calls to special vmult for Operator of size 4: 15 +DEAL::Number of calls to special vmult for Operator of size 9: 27 +DEAL::Number of calls to special vmult for Operator of size 25: 27 +DEAL::Number of calls to special vmult for Operator of size 81: 27 +DEAL::Testing FE_Q<2>(1) +DEAL::Number of degrees of freedom: 289 +DEAL:cg::Starting value 15.0000 +DEAL:cg::Convergence step 3 value 1.49408e-06 +DEAL::Number of calls to special vmult for Operator of size 4: 15 +DEAL::Number of calls to special vmult for Operator of size 9: 27 +DEAL::Number of calls to special vmult for Operator of size 25: 27 +DEAL::Number of calls to special vmult for Operator of size 81: 27 +DEAL::Number of calls to special vmult for Operator of size 289: 27 +DEAL::Testing FE_Q<2>(3) +DEAL::Number of degrees of freedom: 625 +DEAL:cg::Starting value 23.0000 +DEAL:cg::Convergence step 4 value 6.30917e-08 +DEAL::Number of calls to special vmult for Operator of size 16: 20 +DEAL::Number of calls to special vmult for Operator of size 49: 36 +DEAL::Number of calls to special vmult for Operator of size 169: 36 +DEAL::Number of calls to special vmult for Operator of size 625: 36 +DEAL::Testing FE_Q<2>(3) +DEAL::Number of degrees of freedom: 2401 +DEAL:cg::Starting value 47.0000 +DEAL:cg::Convergence step 4 value 1.99839e-07 +DEAL::Number of calls to special vmult for Operator of size 16: 20 +DEAL::Number of calls to special vmult for Operator of size 49: 36 +DEAL::Number of calls to special vmult for Operator of size 169: 36 +DEAL::Number of calls to special vmult for Operator of size 625: 36 +DEAL::Number of calls to special vmult for Operator of size 2401: 36 +DEAL::Testing FE_Q<3>(1) +DEAL::Number of degrees of freedom: 125 +DEAL:cg::Starting value 5.19615 +DEAL:cg::Convergence step 3 value 5.16655e-09 +DEAL::Number of calls to special vmult for Operator of size 8: 15 +DEAL::Number of calls to special vmult for Operator of size 27: 27 +DEAL::Number of calls to special vmult for Operator of size 125: 27 +DEAL::Testing FE_Q<3>(1) +DEAL::Number of degrees of freedom: 729 +DEAL:cg::Starting value 18.5203 +DEAL:cg::Convergence step 3 value 1.01773e-06 +DEAL::Number of calls to special vmult for Operator of size 8: 15 +DEAL::Number of calls to special vmult for Operator of size 27: 27 +DEAL::Number of calls to special vmult for Operator of size 125: 27 +DEAL::Number of calls to special vmult for Operator of size 729: 27 +DEAL::Testing FE_Q<3>(2) +DEAL::Number of degrees of freedom: 729 +DEAL:cg::Starting value 18.5203 +DEAL:cg::Convergence step 4 value 5.65213e-09 +DEAL::Number of calls to special vmult for Operator of size 27: 20 +DEAL::Number of calls to special vmult for Operator of size 125: 36 +DEAL::Number of calls to special vmult for Operator of size 729: 36 +DEAL::Testing FE_Q<3>(2) +DEAL::Number of degrees of freedom: 4913 +DEAL:cg::Starting value 58.0948 +DEAL:cg::Convergence step 4 value 3.37720e-08 +DEAL::Number of calls to special vmult for Operator of size 27: 20 +DEAL::Number of calls to special vmult for Operator of size 125: 36 +DEAL::Number of calls to special vmult for Operator of size 729: 36 +DEAL::Number of calls to special vmult for Operator of size 4913: 36 -- 2.39.5