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 d919a3b8544b0ff24e57b8147d23900ea06c62cb..4f810cdfadcb63b1311712e258b22ff2cfa2be82 100644 (file)
--- 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 (file)
index 0000000..5ece096
--- /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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/lac/petsc_precondition.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/petsc_sparse_matrix.h>
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+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<double> 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" <<std::endl;
+
+    deallog << "prec A2, matrix A, A2 residual : ";
+    A2.residual(residual, solution, right_hand_side);
+    if (residual.l2_norm()<1.e-6)
+      deallog << residual.l2_norm() << std::endl;
+    else
+      deallog << "OK" << std::endl;
+  }
+
+}
+

diff --git a/tests/petsc/different_matrix_preconditioner.output b/tests/petsc/different_matrix_preconditioner.output
new file mode 100644 (file)
index 0000000..9beab9f
--- /dev/null
+++ b/tests/petsc/different_matrix_preconditioner.output
@@ -0,0 +1,4 @@
+
+DEAL::prec A , matrix A : OK
+DEAL::prec A2, matrix A, A residual : OK
+DEAL::prec A2, matrix A, A2 residual : OK