From: Peter Munch Date: Tue, 29 Dec 2020 07:50:32 +0000 (+0100) Subject: Add test for p-multigrid X-Git-Tag: v9.3.0-rc1~712^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F11433%2Fhead;p=dealii.git Add test for p-multigrid --- diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h index d0b7f9d435..99498540c7 100644 --- a/include/deal.II/matrix_free/tools.h +++ b/include/deal.II/matrix_free/tools.h @@ -91,7 +91,7 @@ namespace MatrixFreeTools n_components, Number, VectorizedArrayType> &) const, - CLASS * owning_class, + const CLASS * owning_class, const unsigned int dof_no = 0, const unsigned int quad_no = 0, const unsigned int first_selected_component = 0); @@ -609,7 +609,7 @@ namespace MatrixFreeTools n_components, Number, VectorizedArrayType> &) const, - CLASS * owning_class, + const CLASS * owning_class, const unsigned int dof_no, const unsigned int quad_no, const unsigned int first_selected_component) diff --git a/tests/multigrid-global-coarsening/multigrid_p_01.cc b/tests/multigrid-global-coarsening/multigrid_p_01.cc new file mode 100644 index 0000000000..623d35b62d --- /dev/null +++ b/tests/multigrid-global-coarsening/multigrid_p_01.cc @@ -0,0 +1,161 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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. +// +// --------------------------------------------------------------------- + + +/** + * Test p-multigrid for a uniformly refined mesh both for simplex and + * hypercube mesh. + */ + +#include "multigrid_util.h" + +template +void +test(const unsigned int n_refinements, + const unsigned int fe_degree_fine, + const bool do_simplex_mesh) +{ + using VectorType = LinearAlgebra::distributed::Vector; + + Triangulation tria; + if (do_simplex_mesh) + GridGenerator::subdivided_hyper_cube_with_simplices(tria, 2); + else + GridGenerator::subdivided_hyper_cube(tria, 2); + tria.refine_global(n_refinements); + + const auto level_degrees = MGTransferGlobalCoarseningTools::create_p_sequence( + fe_degree_fine, + MGTransferGlobalCoarseningTools::PolynomialSequenceType::bisect); + + const unsigned int min_level = 0; + const unsigned int max_level = level_degrees.size() - 1; + + MGLevelObject> dof_handlers(min_level, max_level, tria); + MGLevelObject> constraints(min_level, max_level); + MGLevelObject> transfers(min_level, + max_level); + MGLevelObject> operators(min_level, max_level); + + std::unique_ptr> mapping_; + + // set up levels + for (auto l = min_level; l <= max_level; ++l) + { + auto &dof_handler = dof_handlers[l]; + auto &constraint = constraints[l]; + auto &op = operators[l]; + + std::unique_ptr> fe; + std::unique_ptr> quad; + std::unique_ptr> mapping; + + if (do_simplex_mesh) + { + fe = std::make_unique>(level_degrees[l]); + quad = std::make_unique>(level_degrees[l] + 1); + mapping = std::make_unique>(Simplex::FE_P(1)); + } + else + { + fe = std::make_unique>(level_degrees[l]); + quad = std::make_unique>(level_degrees[l] + 1); + mapping = std::make_unique>(FE_Q(1)); + } + + if (l == max_level) + mapping_ = mapping->clone(); + + // set up dofhandler + dof_handler.distribute_dofs(*fe); + + // set up constraints + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs(dof_handler, + locally_relevant_dofs); + constraint.reinit(locally_relevant_dofs); + VectorTools::interpolate_boundary_values( + *mapping, dof_handler, 0, Functions::ZeroFunction(), constraint); + constraint.close(); + + // set up operator + op.reinit(*mapping, dof_handler, *quad, constraint); + } + + // set up transfer operator + for (unsigned int l = min_level; l < max_level; ++l) + transfers[l + 1].reinit_polynomial_transfer(dof_handlers[l + 1], + dof_handlers[l], + constraints[l + 1], + constraints[l]); + + MGTransferGlobalCoarsening, VectorType> transfer( + operators, transfers); + + GMGParameters mg_data; // TODO + + VectorType dst, src; + operators[max_level].initialize_dof_vector(dst); + operators[max_level].initialize_dof_vector(src); + + operators[max_level].rhs(src); + + ReductionControl solver_control( + mg_data.maxiter, mg_data.abstol, mg_data.reltol, false, false); + + mg_solve(solver_control, + dst, + src, + mg_data, + dof_handlers[max_level], + operators[max_level], + operators, + transfer); + + deallog << dim << " " << fe_degree_fine << " " << n_refinements << " " + << (do_simplex_mesh ? "tri " : "quad") << " " + << solver_control.last_step() << std::endl; + + DataOut data_out; + + data_out.attach_dof_handler(dof_handlers[max_level]); + data_out.add_data_vector(dst, "solution"); + data_out.build_patches(*mapping_, 2); + + static unsigned int counter = 0; + std::ofstream output("test." + std::to_string(dim) + "." + + std::to_string(counter++) + ".vtk"); + data_out.write_vtk(output); +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll all; + + deallog.precision(8); + + for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements) + for (unsigned int degree = 2; degree <= 4; ++degree) + test<2>(n_refinements, degree, false /*quadrilateral*/); + + return 0; // TODO: enable simplex test once MGTwoLevelTransfer works for + // simplex meshes + + for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements) + for (unsigned int degree = 2; degree <= 2; ++degree) + test<2>(n_refinements, degree, true /*triangle*/); +} diff --git a/tests/multigrid-global-coarsening/multigrid_p_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output b/tests/multigrid-global-coarsening/multigrid_p_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output new file mode 100644 index 0000000000..3031424278 --- /dev/null +++ b/tests/multigrid-global-coarsening/multigrid_p_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output @@ -0,0 +1,10 @@ + +DEAL:0::2 2 2 quad 3 +DEAL:0::2 3 2 quad 3 +DEAL:0::2 4 2 quad 3 +DEAL:0::2 2 3 quad 3 +DEAL:0::2 3 3 quad 3 +DEAL:0::2 4 3 quad 3 +DEAL:0::2 2 4 quad 3 +DEAL:0::2 3 4 quad 3 +DEAL:0::2 4 4 quad 3 diff --git a/tests/multigrid-global-coarsening/multigrid_util.h b/tests/multigrid-global-coarsening/multigrid_util.h new file mode 100644 index 0000000000..e3176b616c --- /dev/null +++ b/tests/multigrid-global-coarsening/multigrid_util.h @@ -0,0 +1,465 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii_tests_multigrid_util_h +#define dealii_tests_multigrid_util_h + +#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 +#include +#include + +#include "../tests.h" + +namespace MGTransferGlobalCoarseningTools +{ + enum class PolynomialSequenceType + { + bisect, + decrease_by_one, + go_to_one + }; + + unsigned int + generate_level_degree(const unsigned int previous_fe_degree, + const PolynomialSequenceType &p_sequence) + { + switch (p_sequence) + { + case PolynomialSequenceType::bisect: + return std::max(previous_fe_degree / 2, 1u); + case PolynomialSequenceType::decrease_by_one: + return std::max(previous_fe_degree - 1, 1u); + case PolynomialSequenceType::go_to_one: + return 1u; + default: + Assert(false, StandardExceptions::ExcNotImplemented()); + return 1u; + } + } + + std::vector + create_p_sequence(const unsigned int degree, + const PolynomialSequenceType &p_sequence) + { + std::vector degrees; + degrees.push_back(degree); + + unsigned int previous_fe_degree = degree; + while (previous_fe_degree > 1) + { + const unsigned int level_degree = + generate_level_degree(previous_fe_degree, p_sequence); + + degrees.push_back(level_degree); + previous_fe_degree = level_degree; + } + + std::reverse(degrees.begin(), degrees.end()); + + return degrees; + } + +} // namespace MGTransferGlobalCoarseningTools + +template +class Operator : public Subscriptor +{ +public: + using value_type = Number; + using number = Number; + using VectorType = LinearAlgebra::distributed::Vector; + + static const int dim = dim_; + + using FECellIntegrator = FEEvaluation; + + void + reinit(const Mapping & mapping, + const DoFHandler & dof_handler, + const Quadrature & quad, + const AffineConstraints &constraints) + { + // Clear internal data structures (if operator is reused). + this->system_matrix.clear(); + + // Copy the constrains, since they might be needed for computation of the + // system matrix later on. + this->constraints.copy_from(constraints); + + // Set up MatrixFree. At the quadrature points, we only need to evaluate + // the gradient of the solution and test with the gradient of the shape + // functions so that we only need to set the flag `update_gradients`. + typename MatrixFree::AdditionalData data; + data.mapping_update_flags = update_gradients; + + matrix_free.reinit(mapping, dof_handler, constraints, quad, data); + } + + virtual types::global_dof_index + m() const + { + return matrix_free.get_dof_handler().n_dofs(); + } + + Number + el(unsigned int, unsigned int) const + { + Assert(false, ExcNotImplemented()); + return 0; + } + + virtual void + initialize_dof_vector(VectorType &vec) const + { + matrix_free.initialize_dof_vector(vec); + } + + virtual void + vmult(VectorType &dst, const VectorType &src) const + { + this->matrix_free.cell_loop( + &Operator::do_cell_integral_range, this, dst, src, true); + } + + void + Tvmult(VectorType &dst, const VectorType &src) const + { + vmult(dst, src); + } + + void + compute_inverse_diagonal(VectorType &diagonal) const + { + // compute diagonal + MatrixFreeTools::compute_diagonal(matrix_free, + diagonal, + &Operator::do_cell_integral_local, + this); + + // and invert it + for (auto &i : diagonal) + i = (std::abs(i) > 1.0e-10) ? (1.0 / i) : 1.0; + } + + virtual const TrilinosWrappers::SparseMatrix & + get_system_matrix() const + { + // Check if matrix has already been set up. + if (system_matrix.m() == 0 && system_matrix.n() == 0) + { + // Set up sparsity pattern of system matrix. + const auto &dof_handler = this->matrix_free.get_dof_handler(); + + TrilinosWrappers::SparsityPattern dsp( + dof_handler.locally_owned_dofs(), + matrix_free.get_task_info().communicator); + + DoFTools::make_sparsity_pattern(dof_handler, dsp, this->constraints); + + dsp.compress(); + system_matrix.reinit(dsp); + + // Assemble system matrix. + MatrixFreeTools::compute_matrix(matrix_free, + constraints, + system_matrix, + &Operator::do_cell_integral_local, + this); + } + + return this->system_matrix; + } + + void + rhs(VectorType &b) const + { + const int dummy = 0; + + matrix_free.template cell_loop( + [](const auto &matrix_free, auto &dst, const auto &, const auto cells) { + FECellIntegrator phi(matrix_free, cells); + for (unsigned int cell = cells.first; cell < cells.second; ++cell) + { + phi.reinit(cell); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_value(1.0, q); + + phi.integrate_scatter(EvaluationFlags::values, dst); + } + }, + b, + dummy, + true); + } + +private: + void + do_cell_integral_local(FECellIntegrator &integrator) const + { + integrator.evaluate(EvaluationFlags::gradients); + + for (unsigned int q = 0; q < integrator.n_q_points; ++q) + integrator.submit_gradient(integrator.get_gradient(q), q); + + integrator.integrate(EvaluationFlags::gradients); + } + + void + do_cell_integral_global(FECellIntegrator &integrator, + VectorType & dst, + const VectorType &src) const + { + integrator.gather_evaluate(src, EvaluationFlags::gradients); + + for (unsigned int q = 0; q < integrator.n_q_points; ++q) + integrator.submit_gradient(integrator.get_gradient(q), q); + + integrator.integrate_scatter(EvaluationFlags::gradients, dst); + } + + void + do_cell_integral_range( + const MatrixFree & matrix_free, + VectorType & dst, + const VectorType & src, + const std::pair &range) const + { + FECellIntegrator integrator(matrix_free, range); + + for (unsigned cell = range.first; cell < range.second; ++cell) + { + integrator.reinit(cell); + + do_cell_integral_global(integrator, dst, src); + } + } + + MatrixFree matrix_free; + + AffineConstraints constraints; + + mutable TrilinosWrappers::SparseMatrix system_matrix; +}; + +struct GMGParameters +{ + struct CoarseSolverParameters + { + std::string type = "cg_with_amg"; // "cg"; + unsigned int maxiter = 10000; + double abstol = 1e-20; + double reltol = 1e-4; + unsigned int smoother_sweeps = 1; + unsigned int n_cycles = 1; + std::string smoother_type = "ILU"; + }; + + struct SmootherParameters + { + std::string type = "chebyshev"; + double smoothing_range = 20; + unsigned int degree = 5; + unsigned int eig_cg_n_iterations = 20; + }; + + SmootherParameters smoother; + CoarseSolverParameters coarse_solver; + + unsigned int maxiter = 10000; + double abstol = 1e-20; + double reltol = 1e-4; +}; + +template +static void +mg_solve(SolverControl & solver_control, + VectorType & dst, + const VectorType & src, + const GMGParameters & mg_data, + const DoFHandler & dof, + const SystemMatrixType & fine_matrix, + const MGLevelObject &mg_matrices, + const MGTransferType & mg_transfer) +{ + AssertThrow(mg_data.smoother.type == "chebyshev", ExcNotImplemented()); + + const unsigned int min_level = mg_matrices.min_level(); + const unsigned int max_level = mg_matrices.max_level(); + + using Number = typename VectorType::value_type; + using SmootherPreconditionerType = DiagonalMatrix; + using SmootherType = PreconditionChebyshev; + using PreconditionerType = PreconditionMG; + + // Initialize level operators. + mg::Matrix mg_matrix(mg_matrices); + + // Initialize smoothers. + MGLevelObject smoother_data(min_level, + max_level); + + for (unsigned int level = min_level; level <= max_level; level++) + { + smoother_data[level].preconditioner = + std::make_shared(); + mg_matrices[level].compute_inverse_diagonal( + smoother_data[level].preconditioner->get_vector()); + smoother_data[level].smoothing_range = mg_data.smoother.smoothing_range; + smoother_data[level].degree = mg_data.smoother.degree; + smoother_data[level].eig_cg_n_iterations = + mg_data.smoother.eig_cg_n_iterations; + } + + MGSmootherPrecondition mg_smoother; + mg_smoother.initialize(mg_matrices, smoother_data); + + // Initialize coarse-grid solver. + ReductionControl coarse_grid_solver_control(mg_data.coarse_solver.maxiter, + mg_data.coarse_solver.abstol, + mg_data.coarse_solver.reltol, + false, + false); + SolverCG coarse_grid_solver(coarse_grid_solver_control); + + PreconditionIdentity precondition_identity; + PreconditionChebyshev> + precondition_chebyshev; + +#ifdef DEAL_II_WITH_TRILINOS + TrilinosWrappers::PreconditionAMG precondition_amg; +#endif + + std::unique_ptr> mg_coarse; + + if (mg_data.coarse_solver.type == "cg") + { + // CG with identity matrix as preconditioner + + mg_coarse = + std::make_unique, + LevelMatrixType, + PreconditionIdentity>>( + coarse_grid_solver, mg_matrices[min_level], precondition_identity); + } + else if (mg_data.coarse_solver.type == "cg_with_chebyshev") + { + // CG with Chebyshev as preconditioner + + typename SmootherType::AdditionalData smoother_data; + + smoother_data.preconditioner = + std::make_shared>(); + mg_matrices[min_level].compute_inverse_diagonal( + smoother_data.preconditioner->get_vector()); + smoother_data.smoothing_range = mg_data.smoother.smoothing_range; + smoother_data.degree = mg_data.smoother.degree; + smoother_data.eig_cg_n_iterations = mg_data.smoother.eig_cg_n_iterations; + + precondition_chebyshev.initialize(mg_matrices[min_level], smoother_data); + + mg_coarse = std::make_unique< + MGCoarseGridIterativeSolver, + LevelMatrixType, + decltype(precondition_chebyshev)>>( + coarse_grid_solver, mg_matrices[min_level], precondition_chebyshev); + } + else if (mg_data.coarse_solver.type == "cg_with_amg") + { + // CG with AMG as preconditioner + +#ifdef DEAL_II_WITH_TRILINOS + TrilinosWrappers::PreconditionAMG::AdditionalData amg_data; + amg_data.smoother_sweeps = mg_data.coarse_solver.smoother_sweeps; + amg_data.n_cycles = mg_data.coarse_solver.n_cycles; + amg_data.smoother_type = mg_data.coarse_solver.smoother_type.c_str(); + + // CG with AMG as preconditioner + precondition_amg.initialize(mg_matrices[min_level].get_system_matrix(), + amg_data); + + mg_coarse = std::make_unique< + MGCoarseGridIterativeSolver, + LevelMatrixType, + decltype(precondition_amg)>>( + coarse_grid_solver, mg_matrices[min_level], precondition_amg); +#else + AssertThrow(false, ExcNotImplemented()); +#endif + } + else + { + AssertThrow(false, ExcNotImplemented()); + } + + // Create multigrid object. + Multigrid mg( + mg_matrix, *mg_coarse, mg_transfer, mg_smoother, mg_smoother); + + // Convert it to a preconditioner. + PreconditionerType preconditioner(dof, mg, mg_transfer); + + // Finally, solve. + SolverCG(solver_control) + .solve(fine_matrix, dst, src, preconditioner); +} + +#endif