From: Martin Kronbichler Date: Mon, 9 May 2022 15:01:49 +0000 (+0200) Subject: Add two new test cases for interleaved solver X-Git-Tag: v9.4.0-rc1~232^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=73ca892691daf6d2240037eaf48cc230265a0373;p=dealii.git Add two new test cases for interleaved solver --- diff --git a/tests/lac/solver_cg_interleave.cc b/tests/lac/solver_cg_interleave.cc new file mode 100644 index 0000000000..5047c738c2 --- /dev/null +++ b/tests/lac/solver_cg_interleave.cc @@ -0,0 +1,100 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 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. +// +// --------------------------------------------------------------------- + +// Check the path of SolverCG that can apply loops to sub-ranges of the +// matrix. That feature is used for matrix-free loops with increased data +// locality. + + +#include +#include +#include +#include + +#include "../tests.h" + + +struct MyDiagonalMatrix +{ + MyDiagonalMatrix(const LinearAlgebra::distributed::Vector &diagonal) + : diagonal(diagonal) + {} + + void + vmult(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const + { + dst = src; + dst.scale(diagonal); + } + + void + vmult( + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector & src, + const std::function &before, + const std::function &after) + const + { + before(0, dst.size()); + vmult(dst, src); + after(0, dst.size()); + } + + const LinearAlgebra::distributed::Vector &diagonal; +}; + + + +SolverControl::State +monitor_norm(const unsigned int iteration, + const double check_value, + const LinearAlgebra::distributed::Vector &) +{ + deallog << " CG estimated residual at iteration " << iteration << ": " + << check_value << std::endl; + return SolverControl::success; +} + + +int +main() +{ + initlog(); + + // Create diagonal matrix with entries between 1 and 30 + DiagonalMatrix> unit_matrix; + unit_matrix.get_vector().reinit(30); + unit_matrix.get_vector() = 1.0; + + LinearAlgebra::distributed::Vector matrix_entries(unit_matrix.m()); + for (unsigned int i = 0; i < unit_matrix.m(); ++i) + matrix_entries(i) = i + 1; + MyDiagonalMatrix matrix(matrix_entries); + + LinearAlgebra::distributed::Vector rhs(unit_matrix.m()), + sol(unit_matrix.m()); + rhs = 1.; + + deallog << "Solve with PreconditionIdentity: " << std::endl; + SolverControl control(30, 1e-4); + SolverCG> solver(control); + solver.connect(&monitor_norm); + solver.solve(matrix, sol, rhs, PreconditionIdentity()); + + deallog << "Solve with diagonal preconditioner: " << std::endl; + sol = 0; + solver.solve(matrix, sol, rhs, unit_matrix); +} diff --git a/tests/lac/solver_cg_interleave.output b/tests/lac/solver_cg_interleave.output new file mode 100644 index 0000000000..2485887418 --- /dev/null +++ b/tests/lac/solver_cg_interleave.output @@ -0,0 +1,53 @@ + +DEAL::Solve with PreconditionIdentity: +DEAL:cg::Starting value 5.47723 +DEAL:cg:: CG estimated residual at iteration 0: 5.47723 +DEAL:cg:: CG estimated residual at iteration 1: 3.05857 +DEAL:cg:: CG estimated residual at iteration 2: 2.21614 +DEAL:cg:: CG estimated residual at iteration 3: 1.69418 +DEAL:cg:: CG estimated residual at iteration 4: 1.30657 +DEAL:cg:: CG estimated residual at iteration 5: 0.998837 +DEAL:cg:: CG estimated residual at iteration 6: 0.750194 +DEAL:cg:: CG estimated residual at iteration 7: 0.550634 +DEAL:cg:: CG estimated residual at iteration 8: 0.393553 +DEAL:cg:: CG estimated residual at iteration 9: 0.273167 +DEAL:cg:: CG estimated residual at iteration 10: 0.183730 +DEAL:cg:: CG estimated residual at iteration 11: 0.119512 +DEAL:cg:: CG estimated residual at iteration 12: 0.0750441 +DEAL:cg:: CG estimated residual at iteration 13: 0.0454041 +DEAL:cg:: CG estimated residual at iteration 14: 0.0264187 +DEAL:cg:: CG estimated residual at iteration 15: 0.0147526 +DEAL:cg:: CG estimated residual at iteration 16: 0.00788820 +DEAL:cg:: CG estimated residual at iteration 17: 0.00402832 +DEAL:cg:: CG estimated residual at iteration 18: 0.00195897 +DEAL:cg:: CG estimated residual at iteration 19: 0.000904053 +DEAL:cg:: CG estimated residual at iteration 20: 0.000394320 +DEAL:cg:: CG estimated residual at iteration 21: 0.000161750 +DEAL:cg::Convergence step 22 value 6.20175e-05 +DEAL:cg:: CG estimated residual at iteration 22: 6.20175e-05 +DEAL::Solve with diagonal preconditioner: +DEAL:cg::Starting value 5.47723 +DEAL:cg:: CG estimated residual at iteration 0: 5.47723 +DEAL:cg:: CG estimated residual at iteration 1: 3.05857 +DEAL:cg:: CG estimated residual at iteration 2: 2.21614 +DEAL:cg:: CG estimated residual at iteration 3: 1.69418 +DEAL:cg:: CG estimated residual at iteration 4: 1.30657 +DEAL:cg:: CG estimated residual at iteration 5: 0.998837 +DEAL:cg:: CG estimated residual at iteration 6: 0.750194 +DEAL:cg:: CG estimated residual at iteration 7: 0.550634 +DEAL:cg:: CG estimated residual at iteration 8: 0.393553 +DEAL:cg:: CG estimated residual at iteration 9: 0.273167 +DEAL:cg:: CG estimated residual at iteration 10: 0.183730 +DEAL:cg:: CG estimated residual at iteration 11: 0.119512 +DEAL:cg:: CG estimated residual at iteration 12: 0.0750441 +DEAL:cg:: CG estimated residual at iteration 13: 0.0454041 +DEAL:cg:: CG estimated residual at iteration 14: 0.0264187 +DEAL:cg:: CG estimated residual at iteration 15: 0.0147526 +DEAL:cg:: CG estimated residual at iteration 16: 0.00788820 +DEAL:cg:: CG estimated residual at iteration 17: 0.00402832 +DEAL:cg:: CG estimated residual at iteration 18: 0.00195897 +DEAL:cg:: CG estimated residual at iteration 19: 0.000904053 +DEAL:cg:: CG estimated residual at iteration 20: 0.000394320 +DEAL:cg:: CG estimated residual at iteration 21: 0.000161750 +DEAL:cg::Convergence step 22 value 6.20175e-05 +DEAL:cg:: CG estimated residual at iteration 22: 6.20175e-05 diff --git a/tests/matrix_free/solver_cg_interleave.cc b/tests/matrix_free/solver_cg_interleave.cc new file mode 100644 index 0000000000..0063740579 --- /dev/null +++ b/tests/matrix_free/solver_cg_interleave.cc @@ -0,0 +1,234 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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. +// +// --------------------------------------------------------------------- + + + +// this tests the correctness of using the SolverCG class with interleaved +// vector operations, making use of the data locality options available by +// MatrixFree loops. + +#include +#include + +#include + +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include + +#include "../tests.h" + +#include "matrix_vector_mf.h" + + + +template +class HelmholtzOperator : public Subscriptor +{ +public: + using value_type = number; + using VectorType = LinearAlgebra::distributed::Vector; + + HelmholtzOperator(const MatrixFree &matrix_free) + : n_calls_vmult(0) + , data(matrix_free) + {} + + ~HelmholtzOperator() + { + if (n_calls_vmult > 0) + deallog << "Number of calls to special vmult: " << n_calls_vmult + << std::endl; + } + + void + initialize(const Mapping &mapping, DoFHandler &dof_handler) + { + n_calls_vmult = 0; + } + + void + vmult(VectorType &dst, const VectorType &src) const + { + const std::function< + void(const MatrixFree &, + VectorType &, + const VectorType &, + const std::pair &)> + function = helmholtz_operator; + data.cell_loop(function, dst, src, true); + for (unsigned int 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 + { + ++n_calls_vmult; + const std::function< + void(const MatrixFree &, + VectorType &, + const VectorType &, + const std::pair &)> + function = helmholtz_operator; + data.cell_loop( + function, dst, src, operation_before_loop, operation_after_loop); + for (unsigned int i : data.get_constrained_dofs()) + dst.local_element(i) = src.local_element(i); + } + +private: + mutable unsigned int n_calls_vmult; + MatrixFree data; +}; + + +template +struct MyDiagonalMatrix +{ + MyDiagonalMatrix(const LinearAlgebra::distributed::Vector &vec) + : vec(vec) + {} + + void + vmult(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const + { + dst = src; + dst.scale(vec); + } + + const LinearAlgebra::distributed::Vector &vec; +}; + + + +template +void +test(const unsigned int fe_degree) +{ + using VectorType = LinearAlgebra::distributed::Vector; + + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria); + tria.refine_global(6 - dim); + + MappingQ1 mapping; + FE_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + + IndexSet owned_set = dof.locally_owned_dofs(); + IndexSet relevant_set; + AffineConstraints constraints(relevant_set); + typename MatrixFree::AdditionalData addit_data; + addit_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + + { + DoFTools::extract_locally_relevant_dofs(dof, relevant_set); + constraints.close(); + + DoFRenumbering::matrix_free_data_locality(dof, constraints, addit_data); + } + + constraints.clear(); + DoFTools::extract_locally_relevant_dofs(dof, relevant_set); + constraints.close(); + + const QGauss<1> quadrature(dof.get_fe().degree + 1); + MatrixFree mf_data; + mf_data.reinit(mapping, dof, constraints, quadrature, addit_data); + + MatrixFreeTest mf(mf_data); + VectorType rhs, sol; + mf_data.initialize_dof_vector(rhs); + mf_data.initialize_dof_vector(sol); + rhs = 1.; + DiagonalMatrix preconditioner; + mf_data.initialize_dof_vector(preconditioner.get_vector()); + mf.vmult(preconditioner.get_vector(), rhs); + for (number &a : preconditioner.get_vector()) + if (a != 0.) + a = 1. / a; + else + a = 0.; + + // Step 1: solve with CG solver for a matrix that does not support the + // interleaving operation + { + SolverControl control(200, 1e-2 * rhs.l2_norm()); + SolverCG solver(control); + solver.solve(mf, sol, rhs, preconditioner); + deallog << "Norm of the solution: " << sol.l2_norm() << std::endl; + } + + // Step 2: solve with CG solver for a matrix that does support the + // interleaving operation + { + sol = 0; + HelmholtzOperator matrix(mf_data); + SolverControl control(200, 1e-2 * rhs.l2_norm()); + SolverCG solver(control); + solver.solve(matrix, sol, rhs, preconditioner); + deallog << "Norm of the solution: " << sol.l2_norm() << std::endl; + } + + // Step 3: solve with CG solver for a matrix that does support but a + // preconditioner that does not support the interleaving operation + { + sol = 0; + HelmholtzOperator matrix(mf_data); + MyDiagonalMatrix simple_diagonal(preconditioner.get_vector()); + SolverControl control(200, 1e-2 * rhs.l2_norm()); + SolverCG solver(control); + solver.solve(matrix, sol, rhs, simple_diagonal); + deallog << "Norm of the solution: " << sol.l2_norm() << std::endl; + } +} + + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + mpi_initlog(); + + test<2, double>(3); + test<3, double>(4); + test<3, double>(3); +} diff --git a/tests/matrix_free/solver_cg_interleave.with_p4est=true.mpirun=3.output b/tests/matrix_free/solver_cg_interleave.with_p4est=true.mpirun=3.output new file mode 100644 index 0000000000..7c60fc4942 --- /dev/null +++ b/tests/matrix_free/solver_cg_interleave.with_p4est=true.mpirun=3.output @@ -0,0 +1,31 @@ + +DEAL:cg::Starting value 49.0000 +DEAL:cg::Convergence step 51 value 0.429721 +DEAL::Norm of the solution: 11776.3 +DEAL:cg::Starting value 49.0000 +DEAL:cg::Convergence step 51 value 0.429721 +DEAL::Norm of the solution: 11776.3 +DEAL::Number of calls to special vmult: 51 +DEAL:cg::Starting value 49.0000 +DEAL:cg::Convergence step 51 value 0.429721 +DEAL::Norm of the solution: 11776.3 +DEAL:cg::Starting value 189.571 +DEAL:cg::Convergence step 61 value 1.71195 +DEAL::Norm of the solution: 684335. +DEAL:cg::Starting value 189.571 +DEAL:cg::Convergence step 61 value 1.71197 +DEAL::Norm of the solution: 684335. +DEAL::Number of calls to special vmult: 61 +DEAL:cg::Starting value 189.571 +DEAL:cg::Convergence step 61 value 1.71195 +DEAL::Norm of the solution: 684335. +DEAL:cg::Starting value 125.000 +DEAL:cg::Convergence step 39 value 1.10898 +DEAL::Norm of the solution: 196549. +DEAL:cg::Starting value 125.000 +DEAL:cg::Convergence step 39 value 1.10898 +DEAL::Norm of the solution: 196549. +DEAL::Number of calls to special vmult: 39 +DEAL:cg::Starting value 125.000 +DEAL:cg::Convergence step 39 value 1.10898 +DEAL::Norm of the solution: 196549.