From: Daniel Arndt Date: Wed, 16 Dec 2015 23:25:49 +0000 (+0100) Subject: Fix use of different matrices for linear system and preconditioner using PETSc X-Git-Tag: v8.4.0-rc2~150^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6ee4c8b3758e57a15ade223086ae9c12f5126a81;p=dealii.git Fix use of different matrices for linear system and preconditioner using PETSc --- diff --git a/source/lac/petsc_solver.cc b/source/lac/petsc_solver.cc index d919a3b854..4f810cdfad 100644 --- a/source/lac/petsc_solver.cc +++ b/source/lac/petsc_solver.cc @@ -83,18 +83,6 @@ namespace PETScWrappers ierr = KSPCreate (mpi_communicator, &solver_data->ksp); AssertThrow (ierr == 0, ExcPETScError(ierr)); - // set the matrices involved. the - // last argument is irrelevant here, - // since we use the solver only once - // anyway -#if DEAL_II_PETSC_VERSION_LT(3, 5, 0) - ierr = KSPSetOperators (solver_data->ksp, A, preconditioner, - SAME_PRECONDITIONER); -#else - ierr = KSPSetOperators (solver_data->ksp, A, preconditioner); -#endif - AssertThrow (ierr == 0, ExcPETScError(ierr)); - // let derived classes set the solver // type, and the preconditioning // object set the type of @@ -104,6 +92,18 @@ namespace PETScWrappers ierr = KSPSetPC (solver_data->ksp, preconditioner.get_pc()); AssertThrow (ierr == 0, ExcPETScError(ierr)); + // setting the preconditioner overwrites the used matrices. + // hence, we need to set the matrices after the preconditioner. +#if DEAL_II_PETSC_VERSION_LT(3, 5, 0) + // the last argument is irrelevant here, + // since we use the solver only once anyway + ierr = KSPSetOperators (solver_data->ksp, A, preconditioner, + SAME_PRECONDITIONER); +#else + ierr = KSPSetOperators (solver_data->ksp, A, preconditioner); +#endif + AssertThrow (ierr == 0, ExcPETScError(ierr)); + // then a convergence monitor // function. that function simply // checks with the solver_control diff --git a/tests/petsc/different_matrix_preconditioner.cc b/tests/petsc/different_matrix_preconditioner.cc new file mode 100644 index 0000000000..5ece096cc0 --- /dev/null +++ b/tests/petsc/different_matrix_preconditioner.cc @@ -0,0 +1,134 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2014 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// If different matrices for building the preconditioner and the linear system +// to be solved are used, check whether the correct equation is solved +// https://github.com/dealii/dealii/issues/1978 + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; + +int main(int argc, char **argv) +{ + dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-6); + + // create stiffness and mass matrix templates + + Triangulation<2> triangulation; + GridGenerator::hyper_cube(triangulation, 0, 1.); + triangulation.refine_global(5); + + FE_Q<2> fe_q(1); + DoFHandler<2> dof_handler; + dof_handler.initialize(triangulation, fe_q); + + QGauss<2> quadrature(2); + + SparsityPattern sparsity_pattern(dof_handler.n_dofs(), dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, sparsity_pattern); + sparsity_pattern.compress(); + + SparseMatrix A_tmp; + A_tmp.reinit(sparsity_pattern); + MatrixCreator::create_mass_matrix(dof_handler, quadrature, A_tmp); + + // convert to PETSc types + + PETScWrappers::SparseMatrix A; + PETScWrappers::SparseMatrix A2; + A.reinit(sparsity_pattern); + A2.reinit(sparsity_pattern); + for (auto it : A_tmp) + { + const auto row = it.row(); + const auto column = it.column(); + A.set(row, column, it.value()); + A2.set(row, column, 1000. * it.value()); // scale A2 by 1000. + } + A.compress(VectorOperation::insert); + A2.compress(VectorOperation::insert); + + const unsigned int n_dofs = dof_handler.n_dofs(); + + PETScWrappers::MPI::Vector right_hand_side (MPI_COMM_WORLD, n_dofs, n_dofs); + for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i) + right_hand_side[i] = 1.; + right_hand_side.compress(VectorOperation::insert); + + // solve + PETScWrappers::MPI::Vector solution(MPI_COMM_WORLD, n_dofs, n_dofs); + PETScWrappers::MPI::Vector residual(MPI_COMM_WORLD, n_dofs, n_dofs); + + { + SolverControl solver_control(n_dofs, 1e-6, false, false); + PETScWrappers::SolverGMRES solver(solver_control); + + PETScWrappers::PreconditionBoomerAMG preconditioner(A); + + solver.solve(A, solution, right_hand_side, preconditioner); + + deallog << "prec A , matrix A : "; + A.residual(residual, solution, right_hand_side); + if (residual.l2_norm()>1.e-6) + deallog << residual.l2_norm() << std::endl; + else + deallog << "OK" << std::endl; + } + solution = 0.; + { + SolverControl solver_control(n_dofs, 1e-6, false, false); + PETScWrappers::SolverGMRES solver(solver_control); + + PETScWrappers::PreconditionBoomerAMG preconditioner(A2); + + solver.solve(A, solution, right_hand_side, preconditioner); + + deallog << "prec A2, matrix A, A residual : "; + A.residual(residual, solution, right_hand_side); + if (residual.l2_norm()>1.e-6) + deallog << residual.l2_norm() << std::endl; + else + deallog << "OK" <