libMesh
petsc_linear_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2019 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 #include "libmesh/libmesh_common.h"
19 
20 #ifdef LIBMESH_HAVE_PETSC
21 
22 
23 // C++ includes
24 #include <string.h>
25 
26 // Local Includes
27 #include "libmesh/dof_map.h"
28 #include "libmesh/libmesh_logging.h"
29 #include "libmesh/petsc_linear_solver.h"
30 #include "libmesh/shell_matrix.h"
31 #include "libmesh/petsc_matrix.h"
32 #include "libmesh/petsc_preconditioner.h"
33 #include "libmesh/petsc_vector.h"
34 #include "libmesh/enum_to_string.h"
35 #include "libmesh/system.h"
36 #include "libmesh/petsc_auto_fieldsplit.h"
37 #include "libmesh/solver_configuration.h"
38 #include "libmesh/enum_preconditioner_type.h"
39 #include "libmesh/enum_solver_type.h"
40 #include "libmesh/enum_convergence_flags.h"
41 #include "libmesh/auto_ptr.h" // libmesh_make_unique
42 
43 namespace libMesh
44 {
45 
46 extern "C"
47 {
48 #if PETSC_RELEASE_LESS_THAN(3,0,1)
49  PetscErrorCode libmesh_petsc_preconditioner_setup (void * ctx)
50  {
51  Preconditioner<Number> * preconditioner = static_cast<Preconditioner<Number> *>(ctx);
52 
53  if (!preconditioner->initialized())
54  libmesh_error_msg("Preconditioner not initialized! Make sure you call init() before solve!");
55 
56  preconditioner->setup();
57 
58  return 0;
59  }
60 
61 
62  PetscErrorCode libmesh_petsc_preconditioner_apply(void * ctx, Vec x, Vec y)
63  {
64  Preconditioner<Number> * preconditioner = static_cast<Preconditioner<Number> *>(ctx);
65 
66  PetscVector<Number> x_vec(x, preconditioner->comm());
67  PetscVector<Number> y_vec(y, preconditioner->comm());
68 
69  preconditioner->apply(x_vec,y_vec);
70 
71  return 0;
72  }
73 #else
74  PetscErrorCode libmesh_petsc_preconditioner_setup (PC pc)
75  {
76  void * ctx;
77  PetscErrorCode ierr = PCShellGetContext(pc,&ctx);CHKERRQ(ierr);
78  Preconditioner<Number> * preconditioner = static_cast<Preconditioner<Number> *>(ctx);
79 
80  if (!preconditioner->initialized())
81  libmesh_error_msg("Preconditioner not initialized! Make sure you call init() before solve!");
82 
83  preconditioner->setup();
84 
85  return 0;
86  }
87 
88  PetscErrorCode libmesh_petsc_preconditioner_apply(PC pc, Vec x, Vec y)
89  {
90  void * ctx;
91  PetscErrorCode ierr = PCShellGetContext(pc,&ctx);CHKERRQ(ierr);
92  Preconditioner<Number> * preconditioner = static_cast<Preconditioner<Number> *>(ctx);
93 
94  PetscVector<Number> x_vec(x, preconditioner->comm());
95  PetscVector<Number> y_vec(y, preconditioner->comm());
96 
97  preconditioner->apply(x_vec,y_vec);
98 
99  return 0;
100  }
101 #endif
102 
103 #ifdef LIBMESH_ENABLE_DEPRECATED
104 #if PETSC_RELEASE_LESS_THAN(3,0,1)
105  PetscErrorCode __libmesh_petsc_preconditioner_setup (void * ctx)
106  {
107  libmesh_deprecated();
108  return libmesh_petsc_preconditioner_setup(ctx);
109  }
110 
111  PetscErrorCode __libmesh_petsc_preconditioner_apply(void * ctx, Vec x, Vec y)
112  {
113  libmesh_deprecated();
114  return libmesh_petsc_preconditioner_apply(ctx, x, y);
115  }
116 
117 #else
118  PetscErrorCode __libmesh_petsc_preconditioner_setup (PC pc)
119  {
120  libmesh_deprecated();
121  return libmesh_petsc_preconditioner_setup(pc);
122  }
123 
124  PetscErrorCode __libmesh_petsc_preconditioner_apply(PC pc, Vec x, Vec y)
125  {
126  libmesh_deprecated();
127  return libmesh_petsc_preconditioner_apply(pc, x, y);
128  }
129 #endif
130 #endif
131 } // end extern "C"
132 
133 /*----------------------- functions ----------------------------------*/
134 template <typename T>
135 PetscLinearSolver<T>::PetscLinearSolver(const libMesh::Parallel::Communicator & comm_in) :
136  LinearSolver<T>(comm_in),
137  _restrict_solve_to_is(nullptr),
138  _restrict_solve_to_is_complement(nullptr),
139  _subset_solve_mode(SUBSET_ZERO)
140 {
141  if (this->n_processors() == 1)
142  this->_preconditioner_type = ILU_PRECOND;
143  else
144  this->_preconditioner_type = BLOCK_JACOBI_PRECOND;
145 }
146 
147 
148 
149 template <typename T>
150 void PetscLinearSolver<T>::clear ()
151 {
152  if (this->initialized())
153  {
154  /* If we were restricted to some subset, this restriction must
155  be removed and the subset index set destroyed. */
156  if (_restrict_solve_to_is != nullptr)
157  {
158  PetscErrorCode ierr = LibMeshISDestroy(&_restrict_solve_to_is);
159  LIBMESH_CHKERR(ierr);
160  _restrict_solve_to_is = nullptr;
161  }
162 
163  if (_restrict_solve_to_is_complement != nullptr)
164  {
165  PetscErrorCode ierr = LibMeshISDestroy(&_restrict_solve_to_is_complement);
166  LIBMESH_CHKERR(ierr);
167  _restrict_solve_to_is_complement = nullptr;
168  }
169 
170  this->_is_initialized = false;
171 
172  PetscErrorCode ierr=0;
173 
174  ierr = LibMeshKSPDestroy(&_ksp);
175  LIBMESH_CHKERR(ierr);
176 
177  // Mimic PETSc default solver and preconditioner
178  this->_solver_type = GMRES;
179 
180  if (!this->_preconditioner)
181  {
182  if (this->n_processors() == 1)
183  this->_preconditioner_type = ILU_PRECOND;
184  else
185  this->_preconditioner_type = BLOCK_JACOBI_PRECOND;
186  }
187  }
188 }
189 
190 
191 
192 template <typename T>
193 void PetscLinearSolver<T>::init (const char * name)
194 {
195  // Initialize the data structures if not done so already.
196  if (!this->initialized())
197  {
198  this->_is_initialized = true;
199 
200  PetscErrorCode ierr=0;
201 
202  // Create the linear solver context
203  ierr = KSPCreate (this->comm().get(), &_ksp);
204  LIBMESH_CHKERR(ierr);
205 
206  if (name)
207  {
208  ierr = KSPSetOptionsPrefix(_ksp, name);
209  LIBMESH_CHKERR(ierr);
210  }
211 
212  // Create the preconditioner context
213  ierr = KSPGetPC (_ksp, &_pc);
214  LIBMESH_CHKERR(ierr);
215 
216  // Set user-specified solver and preconditioner types
217  this->set_petsc_solver_type();
218 
219  // If the SolverConfiguration object is provided, use it to set
220  // options during solver initialization.
221  if (this->_solver_configuration)
222  {
223  this->_solver_configuration->set_options_during_init();
224  }
225 
226  // Set the options from user-input
227  // Set runtime options, e.g.,
228  // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
229  // These options will override those specified above as long as
230  // KSPSetFromOptions() is called _after_ any other customization
231  // routines.
232 
233  ierr = KSPSetFromOptions (_ksp);
234  LIBMESH_CHKERR(ierr);
235 
236  // Not sure if this is necessary, or if it is already handled by KSPSetFromOptions?
237  // NOT NECESSARY!!!!
238  //ierr = PCSetFromOptions (_pc);
239  //LIBMESH_CHKERR(ierr);
240 
241  // Have the Krylov subspace method use our good initial guess
242  // rather than 0, unless the user requested a KSPType of
243  // preonly, which complains if asked to use initial guesses.
244 #if PETSC_VERSION_LESS_THAN(3,0,0) || !PETSC_RELEASE_LESS_THAN(3,4,0)
245  // Pre-3.0 and petsc-dev (as of October 2012) use non-const versions
246  KSPType ksp_type;
247 #else
248  const KSPType ksp_type;
249 #endif
250 
251  ierr = KSPGetType (_ksp, &ksp_type);
252  LIBMESH_CHKERR(ierr);
253 
254  if (strcmp(ksp_type, "preonly"))
255  {
256  ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
257  LIBMESH_CHKERR(ierr);
258  }
259 
260  // Notify PETSc of location to store residual history.
261  // This needs to be called before any solves, since
262  // it sets the residual history length to zero. The default
263  // behavior is for PETSc to allocate (internally) an array
264  // of size 1000 to hold the residual norm history.
265  ierr = KSPSetResidualHistory(_ksp,
266  PETSC_NULL, // pointer to the array which holds the history
267  PETSC_DECIDE, // size of the array holding the history
268  PETSC_TRUE); // Whether or not to reset the history for each solve.
269  LIBMESH_CHKERR(ierr);
270 
271  PetscPreconditioner<T>::set_petsc_preconditioner_type(this->_preconditioner_type,_pc);
272 
273  //If there is a preconditioner object we need to set the internal setup and apply routines
274  if (this->_preconditioner)
275  {
276  this->_preconditioner->init();
277  PCShellSetContext(_pc,(void *)this->_preconditioner);
278  PCShellSetSetUp(_pc,libmesh_petsc_preconditioner_setup);
279  PCShellSetApply(_pc,libmesh_petsc_preconditioner_apply);
280  }
281  }
282 }
283 
284 
285 template <typename T>
286 void PetscLinearSolver<T>::init (PetscMatrix<T> * matrix,
287  const char * name)
288 {
289  // Initialize the data structures if not done so already.
290  if (!this->initialized())
291  {
292  this->_is_initialized = true;
293 
294  PetscErrorCode ierr=0;
295 
296  // Create the linear solver context
297  ierr = KSPCreate (this->comm().get(), &_ksp);
298  LIBMESH_CHKERR(ierr);
299 
300  if (name)
301  {
302  ierr = KSPSetOptionsPrefix(_ksp, name);
303  LIBMESH_CHKERR(ierr);
304  }
305 
306  //ierr = PCCreate (this->comm().get(), &_pc);
307  // LIBMESH_CHKERR(ierr);
308 
309  // Create the preconditioner context
310  ierr = KSPGetPC (_ksp, &_pc);
311  LIBMESH_CHKERR(ierr);
312 
313  // Set operators. The input matrix works as the preconditioning matrix
314 #if PETSC_RELEASE_LESS_THAN(3,5,0)
315  ierr = KSPSetOperators(_ksp, matrix->mat(), matrix->mat(),DIFFERENT_NONZERO_PATTERN);
316 #else
317  ierr = KSPSetOperators(_ksp, matrix->mat(), matrix->mat());
318 #endif
319  LIBMESH_CHKERR(ierr);
320 
321  // Set user-specified solver and preconditioner types
322  this->set_petsc_solver_type();
323 
324  // If the SolverConfiguration object is provided, use it to set
325  // options during solver initialization.
326  if (this->_solver_configuration)
327  {
328  this->_solver_configuration->set_options_during_init();
329  }
330 
331  // Set the options from user-input
332  // Set runtime options, e.g.,
333  // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
334  // These options will override those specified above as long as
335  // KSPSetFromOptions() is called _after_ any other customization
336  // routines.
337 
338  ierr = KSPSetFromOptions (_ksp);
339  LIBMESH_CHKERR(ierr);
340 
341  // Not sure if this is necessary, or if it is already handled by KSPSetFromOptions?
342  // NOT NECESSARY!!!!
343  //ierr = PCSetFromOptions (_pc);
344  //LIBMESH_CHKERR(ierr);
345 
346  // Have the Krylov subspace method use our good initial guess
347  // rather than 0, unless the user requested a KSPType of
348  // preonly, which complains if asked to use initial guesses.
349 #if PETSC_VERSION_LESS_THAN(3,0,0) || !PETSC_RELEASE_LESS_THAN(3,4,0)
350  KSPType ksp_type;
351 #else
352  const KSPType ksp_type;
353 #endif
354 
355  ierr = KSPGetType (_ksp, &ksp_type);
356  LIBMESH_CHKERR(ierr);
357 
358  if (strcmp(ksp_type, "preonly"))
359  {
360  ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
361  LIBMESH_CHKERR(ierr);
362  }
363 
364  // Notify PETSc of location to store residual history.
365  // This needs to be called before any solves, since
366  // it sets the residual history length to zero. The default
367  // behavior is for PETSc to allocate (internally) an array
368  // of size 1000 to hold the residual norm history.
369  ierr = KSPSetResidualHistory(_ksp,
370  PETSC_NULL, // pointer to the array which holds the history
371  PETSC_DECIDE, // size of the array holding the history
372  PETSC_TRUE); // Whether or not to reset the history for each solve.
373  LIBMESH_CHKERR(ierr);
374 
375  PetscPreconditioner<T>::set_petsc_preconditioner_type(this->_preconditioner_type,_pc);
376  if (this->_preconditioner)
377  {
378  this->_preconditioner->set_matrix(*matrix);
379  this->_preconditioner->init();
380  PCShellSetContext(_pc,(void *)this->_preconditioner);
381  PCShellSetSetUp(_pc,libmesh_petsc_preconditioner_setup);
382  PCShellSetApply(_pc,libmesh_petsc_preconditioner_apply);
383  }
384  }
385 }
386 
387 
388 
389 template <typename T>
390 void
391 PetscLinearSolver<T>::init_names (const System & sys)
392 {
393  petsc_auto_fieldsplit(this->pc(), sys);
394 }
395 
396 
397 
398 template <typename T>
399 void
400 PetscLinearSolver<T>::restrict_solve_to (const std::vector<unsigned int> * const dofs,
401  const SubsetSolveMode subset_solve_mode)
402 {
403  PetscErrorCode ierr=0;
404 
405  /* The preconditioner (in particular if a default preconditioner)
406  will have to be reset. We call this->clear() to do that. This
407  call will also remove and free any previous subset that may have
408  been set before. */
409  this->clear();
410 
411  _subset_solve_mode = subset_solve_mode;
412 
413  if (dofs != nullptr)
414  {
415  PetscInt * petsc_dofs = nullptr;
416  ierr = PetscMalloc(dofs->size()*sizeof(PetscInt), &petsc_dofs);
417  LIBMESH_CHKERR(ierr);
418 
419  for (auto i : index_range(*dofs))
420  petsc_dofs[i] = (*dofs)[i];
421 
422  ierr = ISCreateLibMesh(this->comm().get(),
423  cast_int<PetscInt>(dofs->size()),
424  petsc_dofs, PETSC_OWN_POINTER,
425  &_restrict_solve_to_is);
426  LIBMESH_CHKERR(ierr);
427  }
428 }
429 
430 
431 
432 template <typename T>
433 std::pair<unsigned int, Real>
434 PetscLinearSolver<T>::solve (SparseMatrix<T> & matrix_in,
435  SparseMatrix<T> & precond_in,
436  NumericVector<T> & solution_in,
437  NumericVector<T> & rhs_in,
438  const double tol,
439  const unsigned int m_its)
440 {
441  LOG_SCOPE("solve()", "PetscLinearSolver");
442 
443  // Make sure the data passed in are really of Petsc types
444  PetscMatrix<T> * matrix = cast_ptr<PetscMatrix<T> *>(&matrix_in);
445  PetscMatrix<T> * precond = cast_ptr<PetscMatrix<T> *>(&precond_in);
446  PetscVector<T> * solution = cast_ptr<PetscVector<T> *>(&solution_in);
447  PetscVector<T> * rhs = cast_ptr<PetscVector<T> *>(&rhs_in);
448 
449  this->init (matrix);
450 
451  PetscErrorCode ierr=0;
452  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
453  PetscReal final_resid=0.;
454 
455  // Close the matrices and vectors in case this wasn't already done.
456  matrix->close ();
457  precond->close ();
458  solution->close ();
459  rhs->close ();
460 
461  // // If matrix != precond, then this means we have specified a
462  // // special preconditioner, so reset preconditioner type to PCMAT.
463  // if (matrix != precond)
464  // {
465  // this->_preconditioner_type = USER_PRECOND;
466  // this->set_petsc_preconditioner_type ();
467  // }
468 
469  Mat submat = nullptr;
470  Mat subprecond = nullptr;
471  Vec subrhs = nullptr;
472  Vec subsolution = nullptr;
473  VecScatter scatter = nullptr;
474  std::unique_ptr<PetscMatrix<Number>> subprecond_matrix;
475 
476  // Set operators. Also restrict rhs and solution vector to
477  // subdomain if necessary.
478  if (_restrict_solve_to_is != nullptr)
479  {
480  PetscInt is_local_size = this->_restrict_solve_to_is_local_size();
481 
482  ierr = VecCreate(this->comm().get(),&subrhs);
483  LIBMESH_CHKERR(ierr);
484  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
485  LIBMESH_CHKERR(ierr);
486  ierr = VecSetFromOptions(subrhs);
487  LIBMESH_CHKERR(ierr);
488 
489  ierr = VecCreate(this->comm().get(),&subsolution);
490  LIBMESH_CHKERR(ierr);
491  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
492  LIBMESH_CHKERR(ierr);
493  ierr = VecSetFromOptions(subsolution);
494  LIBMESH_CHKERR(ierr);
495 
496  ierr = LibMeshVecScatterCreate(rhs->vec(), _restrict_solve_to_is, subrhs, nullptr, &scatter);
497  LIBMESH_CHKERR(ierr);
498 
499  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
500  LIBMESH_CHKERR(ierr);
501  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
502  LIBMESH_CHKERR(ierr);
503 
504  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
505  LIBMESH_CHKERR(ierr);
506  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
507  LIBMESH_CHKERR(ierr);
508 
509  ierr = LibMeshCreateSubMatrix(matrix->mat(),
510  _restrict_solve_to_is,
511  _restrict_solve_to_is,
512 #if PETSC_VERSION_LESS_THAN(3,1,0)
513  PETSC_DECIDE,
514 #endif
515  MAT_INITIAL_MATRIX,
516  &submat);
517  LIBMESH_CHKERR(ierr);
518 
519  ierr = LibMeshCreateSubMatrix(precond->mat(),
520  _restrict_solve_to_is,
521  _restrict_solve_to_is,
522 #if PETSC_VERSION_LESS_THAN(3,1,0)
523  PETSC_DECIDE,
524 #endif
525  MAT_INITIAL_MATRIX,
526  &subprecond);
527  LIBMESH_CHKERR(ierr);
528 
529  /* Since removing columns of the matrix changes the equation
530  system, we will now change the right hand side to compensate
531  for this. Note that this is not necessary if \p SUBSET_ZERO
532  has been selected. */
533  if (_subset_solve_mode!=SUBSET_ZERO)
534  {
535  _create_complement_is(rhs_in);
536  PetscInt is_complement_local_size =
537  cast_int<PetscInt>(rhs_in.local_size()-is_local_size);
538 
539  Vec subvec1 = nullptr;
540  Mat submat1 = nullptr;
541  VecScatter scatter1 = nullptr;
542 
543  ierr = VecCreate(this->comm().get(),&subvec1);
544  LIBMESH_CHKERR(ierr);
545  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
546  LIBMESH_CHKERR(ierr);
547  ierr = VecSetFromOptions(subvec1);
548  LIBMESH_CHKERR(ierr);
549 
550  ierr = LibMeshVecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1, nullptr, &scatter1);
551  LIBMESH_CHKERR(ierr);
552 
553  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
554  LIBMESH_CHKERR(ierr);
555  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
556  LIBMESH_CHKERR(ierr);
557 
558  ierr = VecScale(subvec1,-1.0);
559  LIBMESH_CHKERR(ierr);
560 
561  ierr = LibMeshCreateSubMatrix(matrix->mat(),
562  _restrict_solve_to_is,
563  _restrict_solve_to_is_complement,
564 #if PETSC_VERSION_LESS_THAN(3,1,0)
565  PETSC_DECIDE,
566 #endif
567  MAT_INITIAL_MATRIX,
568  &submat1);
569  LIBMESH_CHKERR(ierr);
570 
571  ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
572  LIBMESH_CHKERR(ierr);
573 
574  ierr = LibMeshVecScatterDestroy(&scatter1);
575  LIBMESH_CHKERR(ierr);
576  ierr = LibMeshVecDestroy(&subvec1);
577  LIBMESH_CHKERR(ierr);
578  ierr = LibMeshMatDestroy(&submat1);
579  LIBMESH_CHKERR(ierr);
580  }
581 #if PETSC_RELEASE_LESS_THAN(3,5,0)
582  ierr = KSPSetOperators(_ksp, submat, subprecond,
583  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
584 #else
585  ierr = KSPSetOperators(_ksp, submat, subprecond);
586 
587  PetscBool ksp_reuse_preconditioner = this->same_preconditioner ? PETSC_TRUE : PETSC_FALSE;
588  ierr = KSPSetReusePreconditioner(_ksp, ksp_reuse_preconditioner);
589 #endif
590  LIBMESH_CHKERR(ierr);
591 
592  if (this->_preconditioner)
593  {
594  subprecond_matrix = libmesh_make_unique<PetscMatrix<Number>>(subprecond, this->comm());
595  this->_preconditioner->set_matrix(*subprecond_matrix);
596  this->_preconditioner->init();
597  }
598  }
599  else
600  {
601 #if PETSC_RELEASE_LESS_THAN(3,5,0)
602  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
603  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
604 #else
605  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat());
606 
607  PetscBool ksp_reuse_preconditioner = this->same_preconditioner ? PETSC_TRUE : PETSC_FALSE;
608  ierr = KSPSetReusePreconditioner(_ksp, ksp_reuse_preconditioner);
609 #endif
610  LIBMESH_CHKERR(ierr);
611 
612  if (this->_preconditioner)
613  {
614  this->_preconditioner->set_matrix(matrix_in);
615  this->_preconditioner->init();
616  }
617  }
618 
619  // Set the tolerances for the iterative solver. Use the user-supplied
620  // tolerance for the relative residual & leave the others at default values.
621  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
622  PETSC_DEFAULT, max_its);
623  LIBMESH_CHKERR(ierr);
624 
625  // Allow command line options to override anything set programmatically.
626  ierr = KSPSetFromOptions(_ksp);
627  LIBMESH_CHKERR(ierr);
628 
629  // If the SolverConfiguration object is provided, use it to override
630  // solver options.
631  if (this->_solver_configuration)
632  {
633  this->_solver_configuration->configure_solver();
634  }
635 
636  // Solve the linear system
637  if (_restrict_solve_to_is != nullptr)
638  {
639  ierr = KSPSolve (_ksp, subrhs, subsolution);
640  LIBMESH_CHKERR(ierr);
641  }
642  else
643  {
644  ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
645  LIBMESH_CHKERR(ierr);
646  }
647 
648  // Get the number of iterations required for convergence
649  ierr = KSPGetIterationNumber (_ksp, &its);
650  LIBMESH_CHKERR(ierr);
651 
652  // Get the norm of the final residual to return to the user.
653  ierr = KSPGetResidualNorm (_ksp, &final_resid);
654  LIBMESH_CHKERR(ierr);
655 
656  if (_restrict_solve_to_is != nullptr)
657  {
658  switch(_subset_solve_mode)
659  {
660  case SUBSET_ZERO:
661  ierr = VecZeroEntries(solution->vec());
662  LIBMESH_CHKERR(ierr);
663  break;
664 
665  case SUBSET_COPY_RHS:
666  ierr = VecCopy(rhs->vec(),solution->vec());
667  LIBMESH_CHKERR(ierr);
668  break;
669 
670  case SUBSET_DONT_TOUCH:
671  /* Nothing to do here. */
672  break;
673 
674  default:
675  libmesh_error_msg("Invalid subset solve mode = " << _subset_solve_mode);
676  }
677  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
678  LIBMESH_CHKERR(ierr);
679  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
680  LIBMESH_CHKERR(ierr);
681 
682  ierr = LibMeshVecScatterDestroy(&scatter);
683  LIBMESH_CHKERR(ierr);
684 
685  if (this->_preconditioner)
686  {
687  // Before subprecond_matrix gets cleaned up, we should give
688  // the _preconditioner a different matrix.
689  this->_preconditioner->set_matrix(matrix_in);
690  this->_preconditioner->init();
691  }
692 
693  ierr = LibMeshVecDestroy(&subsolution);
694  LIBMESH_CHKERR(ierr);
695  ierr = LibMeshVecDestroy(&subrhs);
696  LIBMESH_CHKERR(ierr);
697  ierr = LibMeshMatDestroy(&submat);
698  LIBMESH_CHKERR(ierr);
699  ierr = LibMeshMatDestroy(&subprecond);
700  LIBMESH_CHKERR(ierr);
701  }
702 
703  // return the # of its. and the final residual norm.
704  return std::make_pair(its, final_resid);
705 }
706 
707 template <typename T>
708 std::pair<unsigned int, Real>
709 PetscLinearSolver<T>::adjoint_solve (SparseMatrix<T> & matrix_in,
710  NumericVector<T> & solution_in,
711  NumericVector<T> & rhs_in,
712  const double tol,
713  const unsigned int m_its)
714 {
715  LOG_SCOPE("solve()", "PetscLinearSolver");
716 
717  // Make sure the data passed in are really of Petsc types
718  PetscMatrix<T> * matrix = cast_ptr<PetscMatrix<T> *>(&matrix_in);
719  // Note that the matrix and precond matrix are the same
720  PetscMatrix<T> * precond = cast_ptr<PetscMatrix<T> *>(&matrix_in);
721  PetscVector<T> * solution = cast_ptr<PetscVector<T> *>(&solution_in);
722  PetscVector<T> * rhs = cast_ptr<PetscVector<T> *>(&rhs_in);
723 
724  this->init (matrix);
725 
726  PetscErrorCode ierr=0;
727  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
728  PetscReal final_resid=0.;
729 
730  // Close the matrices and vectors in case this wasn't already done.
731  matrix->close ();
732  precond->close ();
733  solution->close ();
734  rhs->close ();
735 
736  Mat submat = nullptr;
737  Mat subprecond = nullptr;
738  Vec subrhs = nullptr;
739  Vec subsolution = nullptr;
740  VecScatter scatter = nullptr;
741  std::unique_ptr<PetscMatrix<Number>> subprecond_matrix;
742 
743  // Set operators. Also restrict rhs and solution vector to
744  // subdomain if necessary.
745  if (_restrict_solve_to_is != nullptr)
746  {
747  PetscInt is_local_size = this->_restrict_solve_to_is_local_size();
748 
749  ierr = VecCreate(this->comm().get(),&subrhs);
750  LIBMESH_CHKERR(ierr);
751  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
752  LIBMESH_CHKERR(ierr);
753  ierr = VecSetFromOptions(subrhs);
754  LIBMESH_CHKERR(ierr);
755 
756  ierr = VecCreate(this->comm().get(),&subsolution);
757  LIBMESH_CHKERR(ierr);
758  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
759  LIBMESH_CHKERR(ierr);
760  ierr = VecSetFromOptions(subsolution);
761  LIBMESH_CHKERR(ierr);
762 
763  ierr = LibMeshVecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs, nullptr, &scatter);
764  LIBMESH_CHKERR(ierr);
765 
766  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
767  LIBMESH_CHKERR(ierr);
768  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
769  LIBMESH_CHKERR(ierr);
770 
771  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
772  LIBMESH_CHKERR(ierr);
773  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
774  LIBMESH_CHKERR(ierr);
775 
776  ierr = LibMeshCreateSubMatrix(matrix->mat(),
777  _restrict_solve_to_is,
778  _restrict_solve_to_is,
779 #if PETSC_VERSION_LESS_THAN(3,1,0)
780  PETSC_DECIDE,
781 #endif
782  MAT_INITIAL_MATRIX,
783  &submat);
784  LIBMESH_CHKERR(ierr);
785 
786  ierr = LibMeshCreateSubMatrix(precond->mat(),
787  _restrict_solve_to_is,
788  _restrict_solve_to_is,
789 #if PETSC_VERSION_LESS_THAN(3,1,0)
790  PETSC_DECIDE,
791 #endif
792  MAT_INITIAL_MATRIX,
793  &subprecond);
794  LIBMESH_CHKERR(ierr);
795 
796  /* Since removing columns of the matrix changes the equation
797  system, we will now change the right hand side to compensate
798  for this. Note that this is not necessary if \p SUBSET_ZERO
799  has been selected. */
800  if (_subset_solve_mode!=SUBSET_ZERO)
801  {
802  _create_complement_is(rhs_in);
803  PetscInt is_complement_local_size =
804  cast_int<PetscInt>(rhs_in.local_size()-is_local_size);
805 
806  Vec subvec1 = nullptr;
807  Mat submat1 = nullptr;
808  VecScatter scatter1 = nullptr;
809 
810  ierr = VecCreate(this->comm().get(),&subvec1);
811  LIBMESH_CHKERR(ierr);
812  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
813  LIBMESH_CHKERR(ierr);
814  ierr = VecSetFromOptions(subvec1);
815  LIBMESH_CHKERR(ierr);
816 
817  ierr = LibMeshVecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1, nullptr, &scatter1);
818  LIBMESH_CHKERR(ierr);
819 
820  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
821  LIBMESH_CHKERR(ierr);
822  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
823  LIBMESH_CHKERR(ierr);
824 
825  ierr = VecScale(subvec1,-1.0);
826  LIBMESH_CHKERR(ierr);
827 
828  ierr = LibMeshCreateSubMatrix(matrix->mat(),
829  _restrict_solve_to_is,
830  _restrict_solve_to_is_complement,
831 #if PETSC_VERSION_LESS_THAN(3,1,0)
832  PETSC_DECIDE,
833 #endif
834  MAT_INITIAL_MATRIX,
835  &submat1);
836  LIBMESH_CHKERR(ierr);
837 
838  ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
839  LIBMESH_CHKERR(ierr);
840 
841  ierr = LibMeshVecScatterDestroy(&scatter1);
842  LIBMESH_CHKERR(ierr);
843  ierr = LibMeshVecDestroy(&subvec1);
844  LIBMESH_CHKERR(ierr);
845  ierr = LibMeshMatDestroy(&submat1);
846  LIBMESH_CHKERR(ierr);
847  }
848 #if PETSC_RELEASE_LESS_THAN(3,5,0)
849  ierr = KSPSetOperators(_ksp, submat, subprecond,
850  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
851 #else
852  ierr = KSPSetOperators(_ksp, submat, subprecond);
853 
854  PetscBool ksp_reuse_preconditioner = this->same_preconditioner ? PETSC_TRUE : PETSC_FALSE;
855  ierr = KSPSetReusePreconditioner(_ksp, ksp_reuse_preconditioner);
856 #endif
857  LIBMESH_CHKERR(ierr);
858 
859  if (this->_preconditioner)
860  {
861  subprecond_matrix = libmesh_make_unique<PetscMatrix<Number>>(subprecond, this->comm());
862  this->_preconditioner->set_matrix(*subprecond_matrix);
863  this->_preconditioner->init();
864  }
865  }
866  else
867  {
868 #if PETSC_RELEASE_LESS_THAN(3,5,0)
869  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
870  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
871 #else
872  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat());
873 
874  PetscBool ksp_reuse_preconditioner = this->same_preconditioner ? PETSC_TRUE : PETSC_FALSE;
875  ierr = KSPSetReusePreconditioner(_ksp, ksp_reuse_preconditioner);
876 #endif
877  LIBMESH_CHKERR(ierr);
878 
879  if (this->_preconditioner)
880  {
881  this->_preconditioner->set_matrix(matrix_in);
882  this->_preconditioner->init();
883  }
884  }
885 
886  // Set the tolerances for the iterative solver. Use the user-supplied
887  // tolerance for the relative residual & leave the others at default values.
888  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
889  PETSC_DEFAULT, max_its);
890  LIBMESH_CHKERR(ierr);
891 
892  // Allow command line options to override anything set programmatically.
893  ierr = KSPSetFromOptions(_ksp);
894  LIBMESH_CHKERR(ierr);
895 
896  // Solve the linear system
897  if (_restrict_solve_to_is != nullptr)
898  {
899  ierr = KSPSolveTranspose (_ksp, subrhs, subsolution);
900  LIBMESH_CHKERR(ierr);
901  }
902  else
903  {
904  ierr = KSPSolveTranspose (_ksp, rhs->vec(), solution->vec());
905  LIBMESH_CHKERR(ierr);
906  }
907 
908  // Get the number of iterations required for convergence
909  ierr = KSPGetIterationNumber (_ksp, &its);
910  LIBMESH_CHKERR(ierr);
911 
912  // Get the norm of the final residual to return to the user.
913  ierr = KSPGetResidualNorm (_ksp, &final_resid);
914  LIBMESH_CHKERR(ierr);
915 
916  if (_restrict_solve_to_is != nullptr)
917  {
918  switch(_subset_solve_mode)
919  {
920  case SUBSET_ZERO:
921  ierr = VecZeroEntries(solution->vec());
922  LIBMESH_CHKERR(ierr);
923  break;
924 
925  case SUBSET_COPY_RHS:
926  ierr = VecCopy(rhs->vec(),solution->vec());
927  LIBMESH_CHKERR(ierr);
928  break;
929 
930  case SUBSET_DONT_TOUCH:
931  /* Nothing to do here. */
932  break;
933 
934  default:
935  libmesh_error_msg("Invalid subset solve mode = " << _subset_solve_mode);
936  }
937  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
938  LIBMESH_CHKERR(ierr);
939  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
940  LIBMESH_CHKERR(ierr);
941 
942  ierr = LibMeshVecScatterDestroy(&scatter);
943  LIBMESH_CHKERR(ierr);
944 
945  if (this->_preconditioner)
946  {
947  // Before subprecond_matrix gets cleaned up, we should give
948  // the _preconditioner a different matrix.
949  this->_preconditioner->set_matrix(matrix_in);
950  this->_preconditioner->init();
951  }
952 
953  ierr = LibMeshVecDestroy(&subsolution);
954  LIBMESH_CHKERR(ierr);
955  ierr = LibMeshVecDestroy(&subrhs);
956  LIBMESH_CHKERR(ierr);
957  ierr = LibMeshMatDestroy(&submat);
958  LIBMESH_CHKERR(ierr);
959  ierr = LibMeshMatDestroy(&subprecond);
960  LIBMESH_CHKERR(ierr);
961  }
962 
963  // return the # of its. and the final residual norm.
964  return std::make_pair(its, final_resid);
965 }
966 
967 
968 template <typename T>
969 std::pair<unsigned int, Real>
970 PetscLinearSolver<T>::solve (const ShellMatrix<T> & shell_matrix,
971  NumericVector<T> & solution_in,
972  NumericVector<T> & rhs_in,
973  const double tol,
974  const unsigned int m_its)
975 {
976 
977 #if PETSC_VERSION_LESS_THAN(3,1,0)
978  if (_restrict_solve_to_is != nullptr)
979  libmesh_error_msg("The current implementation of subset solves with " \
980  << "shell matrices requires PETSc version 3.1 or above. " \
981  << "Older PETSc version do not support automatic " \
982  << "submatrix generation of shell matrices.");
983 #endif
984 
985  LOG_SCOPE("solve()", "PetscLinearSolver");
986 
987  // Make sure the data passed in are really of Petsc types
988  PetscVector<T> * solution = cast_ptr<PetscVector<T> *>(&solution_in);
989  PetscVector<T> * rhs = cast_ptr<PetscVector<T> *>(&rhs_in);
990 
991  this->init ();
992 
993  PetscErrorCode ierr=0;
994  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
995  PetscReal final_resid=0.;
996 
997  Mat submat = nullptr;
998  Vec subrhs = nullptr;
999  Vec subsolution = nullptr;
1000  VecScatter scatter = nullptr;
1001 
1002  // Close the matrices and vectors in case this wasn't already done.
1003  solution->close ();
1004  rhs->close ();
1005 
1006  // Prepare the matrix.
1007  Mat mat;
1008  ierr = MatCreateShell(this->comm().get(),
1009  rhs_in.local_size(),
1010  solution_in.local_size(),
1011  rhs_in.size(),
1012  solution_in.size(),
1013  const_cast<void *>(static_cast<const void *>(&shell_matrix)),
1014  &mat);
1015  /* Note that the const_cast above is only necessary because PETSc
1016  does not accept a const void *. Inside the member function
1017  _petsc_shell_matrix() below, the pointer is casted back to a
1018  const ShellMatrix<T> *. */
1019 
1020  LIBMESH_CHKERR(ierr);
1021  ierr = MatShellSetOperation(mat,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
1022  LIBMESH_CHKERR(ierr);
1023  ierr = MatShellSetOperation(mat,MATOP_MULT_ADD,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult_add));
1024  LIBMESH_CHKERR(ierr);
1025  ierr = MatShellSetOperation(mat,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
1026  LIBMESH_CHKERR(ierr);
1027 
1028  // Restrict rhs and solution vectors and set operators. The input
1029  // matrix works as the preconditioning matrix.
1030  if (_restrict_solve_to_is != nullptr)
1031  {
1032  PetscInt is_local_size = this->_restrict_solve_to_is_local_size();
1033 
1034  ierr = VecCreate(this->comm().get(),&subrhs);
1035  LIBMESH_CHKERR(ierr);
1036  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
1037  LIBMESH_CHKERR(ierr);
1038  ierr = VecSetFromOptions(subrhs);
1039  LIBMESH_CHKERR(ierr);
1040 
1041  ierr = VecCreate(this->comm().get(),&subsolution);
1042  LIBMESH_CHKERR(ierr);
1043  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
1044  LIBMESH_CHKERR(ierr);
1045  ierr = VecSetFromOptions(subsolution);
1046  LIBMESH_CHKERR(ierr);
1047 
1048  ierr = LibMeshVecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs, nullptr, &scatter);
1049  LIBMESH_CHKERR(ierr);
1050 
1051  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1052  LIBMESH_CHKERR(ierr);
1053  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1054  LIBMESH_CHKERR(ierr);
1055 
1056  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1057  LIBMESH_CHKERR(ierr);
1058  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1059  LIBMESH_CHKERR(ierr);
1060 
1061 #if !PETSC_VERSION_LESS_THAN(3,1,0)
1062  ierr = LibMeshCreateSubMatrix(mat,
1063  _restrict_solve_to_is,
1064  _restrict_solve_to_is,
1065  MAT_INITIAL_MATRIX,
1066  &submat);
1067  LIBMESH_CHKERR(ierr);
1068 #endif
1069 
1070  /* Since removing columns of the matrix changes the equation
1071  system, we will now change the right hand side to compensate
1072  for this. Note that this is not necessary if \p SUBSET_ZERO
1073  has been selected. */
1074  if (_subset_solve_mode!=SUBSET_ZERO)
1075  {
1076  _create_complement_is(rhs_in);
1077  PetscInt is_complement_local_size =
1078  cast_int<PetscInt>(rhs_in.local_size()-is_local_size);
1079 
1080  Vec subvec1 = nullptr;
1081  Mat submat1 = nullptr;
1082  VecScatter scatter1 = nullptr;
1083 
1084  ierr = VecCreate(this->comm().get(),&subvec1);
1085  LIBMESH_CHKERR(ierr);
1086  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
1087  LIBMESH_CHKERR(ierr);
1088  ierr = VecSetFromOptions(subvec1);
1089  LIBMESH_CHKERR(ierr);
1090 
1091  ierr = LibMeshVecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1, nullptr, &scatter1);
1092  LIBMESH_CHKERR(ierr);
1093 
1094  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1095  LIBMESH_CHKERR(ierr);
1096  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1097  LIBMESH_CHKERR(ierr);
1098 
1099  ierr = VecScale(subvec1,-1.0);
1100  LIBMESH_CHKERR(ierr);
1101 
1102 #if !PETSC_VERSION_LESS_THAN(3,1,0)
1103  ierr = LibMeshCreateSubMatrix(mat,
1104  _restrict_solve_to_is,
1105  _restrict_solve_to_is_complement,
1106  MAT_INITIAL_MATRIX,
1107  &submat1);
1108  LIBMESH_CHKERR(ierr);
1109 #endif
1110 
1111  // The following lines would be correct, but don't work
1112  // correctly in PETSc up to 3.1.0-p5. See discussion in
1113  // petsc-users of Nov 9, 2010.
1114  //
1115  // ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
1116  // LIBMESH_CHKERR(ierr);
1117  //
1118  // We workaround by using a temporary vector. Note that the
1119  // fix in PETsc 3.1.0-p6 uses a temporary vector internally,
1120  // so this is no effective performance loss.
1121  Vec subvec2 = nullptr;
1122  ierr = VecCreate(this->comm().get(),&subvec2);
1123  LIBMESH_CHKERR(ierr);
1124  ierr = VecSetSizes(subvec2,is_local_size,PETSC_DECIDE);
1125  LIBMESH_CHKERR(ierr);
1126  ierr = VecSetFromOptions(subvec2);
1127  LIBMESH_CHKERR(ierr);
1128  ierr = MatMult(submat1,subvec1,subvec2);
1129  LIBMESH_CHKERR(ierr);
1130  ierr = VecAXPY(subrhs,1.0,subvec2);
1131 
1132  ierr = LibMeshVecScatterDestroy(&scatter1);
1133  LIBMESH_CHKERR(ierr);
1134  ierr = LibMeshVecDestroy(&subvec1);
1135  LIBMESH_CHKERR(ierr);
1136  ierr = LibMeshMatDestroy(&submat1);
1137  LIBMESH_CHKERR(ierr);
1138  }
1139 #if PETSC_RELEASE_LESS_THAN(3,5,0)
1140  ierr = KSPSetOperators(_ksp, submat, submat,
1141  DIFFERENT_NONZERO_PATTERN);
1142 #else
1143  ierr = KSPSetOperators(_ksp, submat, submat);
1144 #endif
1145  LIBMESH_CHKERR(ierr);
1146  }
1147  else
1148  {
1149 #if PETSC_RELEASE_LESS_THAN(3,5,0)
1150  ierr = KSPSetOperators(_ksp, mat, mat,
1151  DIFFERENT_NONZERO_PATTERN);
1152 #else
1153  ierr = KSPSetOperators(_ksp, mat, mat);
1154 #endif
1155  LIBMESH_CHKERR(ierr);
1156  }
1157 
1158  // Set the tolerances for the iterative solver. Use the user-supplied
1159  // tolerance for the relative residual & leave the others at default values.
1160  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
1161  PETSC_DEFAULT, max_its);
1162  LIBMESH_CHKERR(ierr);
1163 
1164  // Allow command line options to override anything set programmatically.
1165  ierr = KSPSetFromOptions(_ksp);
1166  LIBMESH_CHKERR(ierr);
1167 
1168  // Solve the linear system
1169  if (_restrict_solve_to_is != nullptr)
1170  {
1171  ierr = KSPSolve (_ksp, subrhs, subsolution);
1172  LIBMESH_CHKERR(ierr);
1173  }
1174  else
1175  {
1176  ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
1177  LIBMESH_CHKERR(ierr);
1178  }
1179 
1180  // Get the number of iterations required for convergence
1181  ierr = KSPGetIterationNumber (_ksp, &its);
1182  LIBMESH_CHKERR(ierr);
1183 
1184  // Get the norm of the final residual to return to the user.
1185  ierr = KSPGetResidualNorm (_ksp, &final_resid);
1186  LIBMESH_CHKERR(ierr);
1187 
1188  if (_restrict_solve_to_is != nullptr)
1189  {
1190  switch(_subset_solve_mode)
1191  {
1192  case SUBSET_ZERO:
1193  ierr = VecZeroEntries(solution->vec());
1194  LIBMESH_CHKERR(ierr);
1195  break;
1196 
1197  case SUBSET_COPY_RHS:
1198  ierr = VecCopy(rhs->vec(),solution->vec());
1199  LIBMESH_CHKERR(ierr);
1200  break;
1201 
1202  case SUBSET_DONT_TOUCH:
1203  /* Nothing to do here. */
1204  break;
1205 
1206  default:
1207  libmesh_error_msg("Invalid subset solve mode = " << _subset_solve_mode);
1208  }
1209  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1210  LIBMESH_CHKERR(ierr);
1211  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1212  LIBMESH_CHKERR(ierr);
1213 
1214  ierr = LibMeshVecScatterDestroy(&scatter);
1215  LIBMESH_CHKERR(ierr);
1216 
1217  ierr = LibMeshVecDestroy(&subsolution);
1218  LIBMESH_CHKERR(ierr);
1219  ierr = LibMeshVecDestroy(&subrhs);
1220  LIBMESH_CHKERR(ierr);
1221  ierr = LibMeshMatDestroy(&submat);
1222  LIBMESH_CHKERR(ierr);
1223  }
1224 
1225  // Destroy the matrix.
1226  ierr = LibMeshMatDestroy(&mat);
1227  LIBMESH_CHKERR(ierr);
1228 
1229  // return the # of its. and the final residual norm.
1230  return std::make_pair(its, final_resid);
1231 }
1232 
1233 
1234 
1235 template <typename T>
1236 std::pair<unsigned int, Real>
1237 PetscLinearSolver<T>::solve (const ShellMatrix<T> & shell_matrix,
1238  const SparseMatrix<T> & precond_matrix,
1239  NumericVector<T> & solution_in,
1240  NumericVector<T> & rhs_in,
1241  const double tol,
1242  const unsigned int m_its)
1243 {
1244 
1245 #if PETSC_VERSION_LESS_THAN(3,1,0)
1246  if (_restrict_solve_to_is != nullptr)
1247  libmesh_error_msg("The current implementation of subset solves with " \
1248  << "shell matrices requires PETSc version 3.1 or above. " \
1249  << "Older PETSc version do not support automatic " \
1250  << "submatrix generation of shell matrices.");
1251 #endif
1252 
1253  LOG_SCOPE("solve()", "PetscLinearSolver");
1254 
1255  // Make sure the data passed in are really of Petsc types
1256  const PetscMatrix<T> * precond = cast_ptr<const PetscMatrix<T> *>(&precond_matrix);
1257  PetscVector<T> * solution = cast_ptr<PetscVector<T> *>(&solution_in);
1258  PetscVector<T> * rhs = cast_ptr<PetscVector<T> *>(&rhs_in);
1259 
1260  this->init ();
1261 
1262  PetscErrorCode ierr=0;
1263  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
1264  PetscReal final_resid=0.;
1265 
1266  Mat submat = nullptr;
1267  Mat subprecond = nullptr;
1268  Vec subrhs = nullptr;
1269  Vec subsolution = nullptr;
1270  VecScatter scatter = nullptr;
1271  std::unique_ptr<PetscMatrix<Number>> subprecond_matrix;
1272 
1273  // Close the matrices and vectors in case this wasn't already done.
1274  solution->close ();
1275  rhs->close ();
1276 
1277  // Prepare the matrix.
1278  Mat mat;
1279  ierr = MatCreateShell(this->comm().get(),
1280  rhs_in.local_size(),
1281  solution_in.local_size(),
1282  rhs_in.size(),
1283  solution_in.size(),
1284  const_cast<void *>(static_cast<const void *>(&shell_matrix)),
1285  &mat);
1286  /* Note that the const_cast above is only necessary because PETSc
1287  does not accept a const void *. Inside the member function
1288  _petsc_shell_matrix() below, the pointer is casted back to a
1289  const ShellMatrix<T> *. */
1290 
1291  LIBMESH_CHKERR(ierr);
1292  ierr = MatShellSetOperation(mat,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
1293  LIBMESH_CHKERR(ierr);
1294  ierr = MatShellSetOperation(mat,MATOP_MULT_ADD,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult_add));
1295  LIBMESH_CHKERR(ierr);
1296  ierr = MatShellSetOperation(mat,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
1297  LIBMESH_CHKERR(ierr);
1298 
1299  // Restrict rhs and solution vectors and set operators. The input
1300  // matrix works as the preconditioning matrix.
1301  if (_restrict_solve_to_is != nullptr)
1302  {
1303  PetscInt is_local_size = this->_restrict_solve_to_is_local_size();
1304 
1305  ierr = VecCreate(this->comm().get(),&subrhs);
1306  LIBMESH_CHKERR(ierr);
1307  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
1308  LIBMESH_CHKERR(ierr);
1309  ierr = VecSetFromOptions(subrhs);
1310  LIBMESH_CHKERR(ierr);
1311 
1312  ierr = VecCreate(this->comm().get(),&subsolution);
1313  LIBMESH_CHKERR(ierr);
1314  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
1315  LIBMESH_CHKERR(ierr);
1316  ierr = VecSetFromOptions(subsolution);
1317  LIBMESH_CHKERR(ierr);
1318 
1319  ierr = LibMeshVecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs, nullptr, &scatter);
1320  LIBMESH_CHKERR(ierr);
1321 
1322  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1323  LIBMESH_CHKERR(ierr);
1324  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1325  LIBMESH_CHKERR(ierr);
1326 
1327  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1328  LIBMESH_CHKERR(ierr);
1329  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1330  LIBMESH_CHKERR(ierr);
1331 
1332 #if !PETSC_VERSION_LESS_THAN(3,1,0)
1333  ierr = LibMeshCreateSubMatrix(mat,
1334  _restrict_solve_to_is,
1335  _restrict_solve_to_is,
1336  MAT_INITIAL_MATRIX,
1337  &submat);
1338  LIBMESH_CHKERR(ierr);
1339 
1340  ierr = LibMeshCreateSubMatrix(const_cast<PetscMatrix<T> *>(precond)->mat(),
1341  _restrict_solve_to_is,
1342  _restrict_solve_to_is,
1343  MAT_INITIAL_MATRIX,
1344  &subprecond);
1345  LIBMESH_CHKERR(ierr);
1346 #endif
1347 
1348  /* Since removing columns of the matrix changes the equation
1349  system, we will now change the right hand side to compensate
1350  for this. Note that this is not necessary if \p SUBSET_ZERO
1351  has been selected. */
1352  if (_subset_solve_mode!=SUBSET_ZERO)
1353  {
1354  _create_complement_is(rhs_in);
1355  PetscInt is_complement_local_size = rhs_in.local_size()-is_local_size;
1356 
1357  Vec subvec1 = nullptr;
1358  Mat submat1 = nullptr;
1359  VecScatter scatter1 = nullptr;
1360 
1361  ierr = VecCreate(this->comm().get(),&subvec1);
1362  LIBMESH_CHKERR(ierr);
1363  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
1364  LIBMESH_CHKERR(ierr);
1365  ierr = VecSetFromOptions(subvec1);
1366  LIBMESH_CHKERR(ierr);
1367 
1368  ierr = LibMeshVecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1, nullptr, &scatter1);
1369  LIBMESH_CHKERR(ierr);
1370 
1371  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1372  LIBMESH_CHKERR(ierr);
1373  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1374  LIBMESH_CHKERR(ierr);
1375 
1376  ierr = VecScale(subvec1,-1.0);
1377  LIBMESH_CHKERR(ierr);
1378 
1379 #if !PETSC_VERSION_LESS_THAN(3,1,0)
1380  ierr = LibMeshCreateSubMatrix(mat,
1381  _restrict_solve_to_is,
1382  _restrict_solve_to_is_complement,
1383  MAT_INITIAL_MATRIX,
1384  &submat1);
1385  LIBMESH_CHKERR(ierr);
1386 #endif
1387 
1388  // The following lines would be correct, but don't work
1389  // correctly in PETSc up to 3.1.0-p5. See discussion in
1390  // petsc-users of Nov 9, 2010.
1391  //
1392  // ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
1393  // LIBMESH_CHKERR(ierr);
1394  //
1395  // We workaround by using a temporary vector. Note that the
1396  // fix in PETsc 3.1.0-p6 uses a temporary vector internally,
1397  // so this is no effective performance loss.
1398  Vec subvec2 = nullptr;
1399  ierr = VecCreate(this->comm().get(),&subvec2);
1400  LIBMESH_CHKERR(ierr);
1401  ierr = VecSetSizes(subvec2,is_local_size,PETSC_DECIDE);
1402  LIBMESH_CHKERR(ierr);
1403  ierr = VecSetFromOptions(subvec2);
1404  LIBMESH_CHKERR(ierr);
1405  ierr = MatMult(submat1,subvec1,subvec2);
1406  LIBMESH_CHKERR(ierr);
1407  ierr = VecAXPY(subrhs,1.0,subvec2);
1408  LIBMESH_CHKERR(ierr);
1409 
1410  ierr = LibMeshVecScatterDestroy(&scatter1);
1411  LIBMESH_CHKERR(ierr);
1412  ierr = LibMeshVecDestroy(&subvec1);
1413  LIBMESH_CHKERR(ierr);
1414  ierr = LibMeshMatDestroy(&submat1);
1415  LIBMESH_CHKERR(ierr);
1416  }
1417 
1418 #if PETSC_RELEASE_LESS_THAN(3,5,0)
1419  ierr = KSPSetOperators(_ksp, submat, subprecond,
1420  DIFFERENT_NONZERO_PATTERN);
1421 #else
1422  ierr = KSPSetOperators(_ksp, submat, subprecond);
1423 #endif
1424  LIBMESH_CHKERR(ierr);
1425 
1426  if (this->_preconditioner)
1427  {
1428  subprecond_matrix = libmesh_make_unique<PetscMatrix<Number>>(subprecond, this->comm());
1429  this->_preconditioner->set_matrix(*subprecond_matrix);
1430  this->_preconditioner->init();
1431  }
1432  }
1433  else
1434  {
1435 #if PETSC_RELEASE_LESS_THAN(3,5,0)
1436  ierr = KSPSetOperators(_ksp, mat, const_cast<PetscMatrix<T> *>(precond)->mat(),
1437  DIFFERENT_NONZERO_PATTERN);
1438 #else
1439  ierr = KSPSetOperators(_ksp, mat, const_cast<PetscMatrix<T> *>(precond)->mat());
1440 #endif
1441  LIBMESH_CHKERR(ierr);
1442 
1443  if (this->_preconditioner)
1444  {
1445  this->_preconditioner->set_matrix(const_cast<SparseMatrix<Number> &>(precond_matrix));
1446  this->_preconditioner->init();
1447  }
1448  }
1449 
1450  // Set the tolerances for the iterative solver. Use the user-supplied
1451  // tolerance for the relative residual & leave the others at default values.
1452  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
1453  PETSC_DEFAULT, max_its);
1454  LIBMESH_CHKERR(ierr);
1455 
1456  // Allow command line options to override anything set programmatically.
1457  ierr = KSPSetFromOptions(_ksp);
1458  LIBMESH_CHKERR(ierr);
1459 
1460  // Solve the linear system
1461  if (_restrict_solve_to_is != nullptr)
1462  {
1463  ierr = KSPSolve (_ksp, subrhs, subsolution);
1464  LIBMESH_CHKERR(ierr);
1465  }
1466  else
1467  {
1468  ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
1469  LIBMESH_CHKERR(ierr);
1470  }
1471 
1472  // Get the number of iterations required for convergence
1473  ierr = KSPGetIterationNumber (_ksp, &its);
1474  LIBMESH_CHKERR(ierr);
1475 
1476  // Get the norm of the final residual to return to the user.
1477  ierr = KSPGetResidualNorm (_ksp, &final_resid);
1478  LIBMESH_CHKERR(ierr);
1479 
1480  if (_restrict_solve_to_is != nullptr)
1481  {
1482  switch(_subset_solve_mode)
1483  {
1484  case SUBSET_ZERO:
1485  ierr = VecZeroEntries(solution->vec());
1486  LIBMESH_CHKERR(ierr);
1487  break;
1488 
1489  case SUBSET_COPY_RHS:
1490  ierr = VecCopy(rhs->vec(),solution->vec());
1491  LIBMESH_CHKERR(ierr);
1492  break;
1493 
1494  case SUBSET_DONT_TOUCH:
1495  /* Nothing to do here. */
1496  break;
1497 
1498  default:
1499  libmesh_error_msg("Invalid subset solve mode = " << _subset_solve_mode);
1500  }
1501  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1502  LIBMESH_CHKERR(ierr);
1503  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1504  LIBMESH_CHKERR(ierr);
1505 
1506  ierr = LibMeshVecScatterDestroy(&scatter);
1507  LIBMESH_CHKERR(ierr);
1508 
1509  if (this->_preconditioner)
1510  {
1511  // Before subprecond_matrix gets cleaned up, we should give
1512  // the _preconditioner a different matrix.
1513  this->_preconditioner->set_matrix(const_cast<SparseMatrix<Number> &>(precond_matrix));
1514  this->_preconditioner->init();
1515  }
1516 
1517  ierr = LibMeshVecDestroy(&subsolution);
1518  LIBMESH_CHKERR(ierr);
1519  ierr = LibMeshVecDestroy(&subrhs);
1520  LIBMESH_CHKERR(ierr);
1521  ierr = LibMeshMatDestroy(&submat);
1522  LIBMESH_CHKERR(ierr);
1523  ierr = LibMeshMatDestroy(&subprecond);
1524  LIBMESH_CHKERR(ierr);
1525  }
1526 
1527  // Destroy the matrix.
1528  ierr = LibMeshMatDestroy(&mat);
1529  LIBMESH_CHKERR(ierr);
1530 
1531  // return the # of its. and the final residual norm.
1532  return std::make_pair(its, final_resid);
1533 }
1534 
1535 
1536 
1537 template <typename T>
1538 void PetscLinearSolver<T>::get_residual_history(std::vector<double> & hist)
1539 {
1540  PetscErrorCode ierr = 0;
1541  PetscInt its = 0;
1542 
1543  // Fill the residual history vector with the residual norms
1544  // Note that GetResidualHistory() does not copy any values, it
1545  // simply sets the pointer p. Note that for some Krylov subspace
1546  // methods, the number of residuals returned in the history
1547  // vector may be different from what you are expecting. For
1548  // example, TFQMR returns two residual values per iteration step.
1549  PetscReal * p;
1550  ierr = KSPGetResidualHistory(_ksp, &p, &its);
1551  LIBMESH_CHKERR(ierr);
1552 
1553  // Check for early return
1554  if (its == 0) return;
1555 
1556  // Create space to store the result
1557  hist.resize(its);
1558 
1559  // Copy history into the vector provided by the user.
1560  for (PetscInt i=0; i<its; ++i)
1561  {
1562  hist[i] = *p;
1563  p++;
1564  }
1565 }
1566 
1567 
1568 
1569 
1570 template <typename T>
1571 Real PetscLinearSolver<T>::get_initial_residual()
1572 {
1573  PetscErrorCode ierr = 0;
1574  PetscInt its = 0;
1575 
1576  // Fill the residual history vector with the residual norms
1577  // Note that GetResidualHistory() does not copy any values, it
1578  // simply sets the pointer p. Note that for some Krylov subspace
1579  // methods, the number of residuals returned in the history
1580  // vector may be different from what you are expecting. For
1581  // example, TFQMR returns two residual values per iteration step.
1582  PetscReal * p;
1583  ierr = KSPGetResidualHistory(_ksp, &p, &its);
1584  LIBMESH_CHKERR(ierr);
1585 
1586  // Check no residual history
1587  if (its == 0)
1588  {
1589  libMesh::err << "No iterations have been performed, returning 0." << std::endl;
1590  return 0.;
1591  }
1592 
1593  // Otherwise, return the value pointed to by p.
1594  return *p;
1595 }
1596 
1597 
1598 
1599 
1600 template <typename T>
1601 void PetscLinearSolver<T>::set_petsc_solver_type()
1602 {
1603  PetscErrorCode ierr = 0;
1604 
1605  switch (this->_solver_type)
1606  {
1607 
1608  case CG:
1609  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPCG));
1610  LIBMESH_CHKERR(ierr);
1611  return;
1612 
1613  case CR:
1614  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPCR));
1615  LIBMESH_CHKERR(ierr);
1616  return;
1617 
1618  case CGS:
1619  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPCGS));
1620  LIBMESH_CHKERR(ierr);
1621  return;
1622 
1623  case BICG:
1624  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPBICG));
1625  LIBMESH_CHKERR(ierr);
1626  return;
1627 
1628  case TCQMR:
1629  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPTCQMR));
1630  LIBMESH_CHKERR(ierr);
1631  return;
1632 
1633  case TFQMR:
1634  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPTFQMR));
1635  LIBMESH_CHKERR(ierr);
1636  return;
1637 
1638  case LSQR:
1639  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPLSQR));
1640  LIBMESH_CHKERR(ierr);
1641  return;
1642 
1643  case BICGSTAB:
1644  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPBCGS));
1645  LIBMESH_CHKERR(ierr);
1646  return;
1647 
1648  case MINRES:
1649  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPMINRES));
1650  LIBMESH_CHKERR(ierr);
1651  return;
1652 
1653  case GMRES:
1654  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPGMRES));
1655  LIBMESH_CHKERR(ierr);
1656  return;
1657 
1658  case RICHARDSON:
1659  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPRICHARDSON));
1660  LIBMESH_CHKERR(ierr);
1661  return;
1662 
1663  case CHEBYSHEV:
1664 #if defined(LIBMESH_HAVE_PETSC) && PETSC_VERSION_LESS_THAN(3,3,0)
1665  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPCHEBYCHEV));
1666  LIBMESH_CHKERR(ierr);
1667  return;
1668 #else
1669  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPCHEBYSHEV));
1670  LIBMESH_CHKERR(ierr);
1671  return;
1672 #endif
1673 
1674 
1675  default:
1676  libMesh::err << "ERROR: Unsupported PETSC Solver: "
1677  << Utility::enum_to_string(this->_solver_type) << std::endl
1678  << "Continuing with PETSC defaults" << std::endl;
1679  }
1680 }
1681 
1682 
1683 
1684 template <typename T>
1685 LinearConvergenceReason PetscLinearSolver<T>::get_converged_reason() const
1686 {
1687  KSPConvergedReason reason;
1688  KSPGetConvergedReason(_ksp, &reason);
1689 
1690  switch(reason)
1691  {
1692 #if !PETSC_VERSION_LESS_THAN(3,2,0)
1693  case KSP_CONVERGED_RTOL_NORMAL : return CONVERGED_RTOL_NORMAL;
1694  case KSP_CONVERGED_ATOL_NORMAL : return CONVERGED_ATOL_NORMAL;
1695 #endif
1696  case KSP_CONVERGED_RTOL : return CONVERGED_RTOL;
1697  case KSP_CONVERGED_ATOL : return CONVERGED_ATOL;
1698  case KSP_CONVERGED_ITS : return CONVERGED_ITS;
1699  case KSP_CONVERGED_CG_NEG_CURVE : return CONVERGED_CG_NEG_CURVE;
1700  case KSP_CONVERGED_CG_CONSTRAINED : return CONVERGED_CG_CONSTRAINED;
1701  case KSP_CONVERGED_STEP_LENGTH : return CONVERGED_STEP_LENGTH;
1702  case KSP_CONVERGED_HAPPY_BREAKDOWN : return CONVERGED_HAPPY_BREAKDOWN;
1703  case KSP_DIVERGED_NULL : return DIVERGED_NULL;
1704  case KSP_DIVERGED_ITS : return DIVERGED_ITS;
1705  case KSP_DIVERGED_DTOL : return DIVERGED_DTOL;
1706  case KSP_DIVERGED_BREAKDOWN : return DIVERGED_BREAKDOWN;
1707  case KSP_DIVERGED_BREAKDOWN_BICG : return DIVERGED_BREAKDOWN_BICG;
1708  case KSP_DIVERGED_NONSYMMETRIC : return DIVERGED_NONSYMMETRIC;
1709  case KSP_DIVERGED_INDEFINITE_PC : return DIVERGED_INDEFINITE_PC;
1710 #if PETSC_VERSION_LESS_THAN(3,4,0)
1711  case KSP_DIVERGED_NAN : return DIVERGED_NAN;
1712 #else
1713  case KSP_DIVERGED_NANORINF : return DIVERGED_NAN;
1714 #endif
1715  case KSP_DIVERGED_INDEFINITE_MAT : return DIVERGED_INDEFINITE_MAT;
1716  case KSP_CONVERGED_ITERATING : return CONVERGED_ITERATING;
1717 #if !PETSC_VERSION_LESS_THAN(3,7,0)
1718 // PETSc-3.7.0 to 3.10.x
1719 #if PETSC_VERSION_LESS_THAN(3,11,0) && PETSC_VERSION_RELEASE
1720  case KSP_DIVERGED_PCSETUP_FAILED : return DIVERGED_PCSETUP_FAILED;
1721 // PETSc master and future PETSc
1722 #else
1723  case KSP_DIVERGED_PC_FAILED : return DIVERGED_PCSETUP_FAILED;
1724 #endif
1725 #endif
1726  default :
1727  libMesh::err << "Unknown convergence flag!" << std::endl;
1728  return UNKNOWN_FLAG;
1729  }
1730 }
1731 
1732 
1733 template <typename T>
1734 PetscErrorCode PetscLinearSolver<T>::_petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest)
1735 {
1736  /* Get the matrix context. */
1737  PetscErrorCode ierr=0;
1738  void * ctx;
1739  ierr = MatShellGetContext(mat,&ctx);
1740 
1741  /* Get user shell matrix object. */
1742  const ShellMatrix<T> & shell_matrix = *static_cast<const ShellMatrix<T> *>(ctx);
1743  CHKERRABORT(shell_matrix.comm().get(), ierr);
1744 
1745  /* Make \p NumericVector instances around the vectors. */
1746  PetscVector<T> arg_global(arg, shell_matrix.comm());
1747  PetscVector<T> dest_global(dest, shell_matrix.comm());
1748 
1749  /* Call the user function. */
1750  shell_matrix.vector_mult(dest_global,arg_global);
1751 
1752  return ierr;
1753 }
1754 
1755 
1756 
1757 template <typename T>
1758 PetscErrorCode PetscLinearSolver<T>::_petsc_shell_matrix_mult_add(Mat mat, Vec arg, Vec add, Vec dest)
1759 {
1760  /* Get the matrix context. */
1761  PetscErrorCode ierr=0;
1762  void * ctx;
1763  ierr = MatShellGetContext(mat,&ctx);
1764 
1765  /* Get user shell matrix object. */
1766  const ShellMatrix<T> & shell_matrix = *static_cast<const ShellMatrix<T> *>(ctx);
1767  CHKERRABORT(shell_matrix.comm().get(), ierr);
1768 
1769  /* Make \p NumericVector instances around the vectors. */
1770  PetscVector<T> arg_global(arg, shell_matrix.comm());
1771  PetscVector<T> dest_global(dest, shell_matrix.comm());
1772  PetscVector<T> add_global(add, shell_matrix.comm());
1773 
1774  if (add!=arg)
1775  {
1776  arg_global = add_global;
1777  }
1778 
1779  /* Call the user function. */
1780  shell_matrix.vector_mult_add(dest_global,arg_global);
1781 
1782  return ierr;
1783 }
1784 
1785 
1786 
1787 template <typename T>
1788 PetscErrorCode PetscLinearSolver<T>::_petsc_shell_matrix_get_diagonal(Mat mat, Vec dest)
1789 {
1790  /* Get the matrix context. */
1791  PetscErrorCode ierr=0;
1792  void * ctx;
1793  ierr = MatShellGetContext(mat,&ctx);
1794 
1795  /* Get user shell matrix object. */
1796  const ShellMatrix<T> & shell_matrix = *static_cast<const ShellMatrix<T> *>(ctx);
1797  CHKERRABORT(shell_matrix.comm().get(), ierr);
1798 
1799  /* Make \p NumericVector instances around the vector. */
1800  PetscVector<T> dest_global(dest, shell_matrix.comm());
1801 
1802  /* Call the user function. */
1803  shell_matrix.get_diagonal(dest_global);
1804 
1805  return ierr;
1806 }
1807 
1808 
1809 
1810 //------------------------------------------------------------------
1811 // Explicit instantiations
1812 template class PetscLinearSolver<Number>;
1813 
1814 } // namespace libMesh
1815 
1816 
1817 
1818 #endif // #ifdef LIBMESH_HAVE_PETSC
libMesh::ShellMatrix
Generic shell matrix, i.e.
Definition: eigen_preconditioner.h:36
libMesh::System
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:100
libMesh::DIVERGED_BREAKDOWN_BICG
Definition: enum_convergence_flags.h:49
libMesh::CONVERGED_STEP_LENGTH
Definition: enum_convergence_flags.h:42
libMesh::UNKNOWN_FLAG
Definition: enum_convergence_flags.h:58
libMesh::DIVERGED_INDEFINITE_MAT
Definition: enum_convergence_flags.h:53
libMesh::ShellMatrix::get_diagonal
virtual void get_diagonal(NumericVector< T > &dest) const =0
Copies the diagonal part of the matrix into dest.
libMesh::BICGSTAB
Definition: enum_solver_type.h:42
libMesh::DIVERGED_BREAKDOWN
Definition: enum_convergence_flags.h:48
libMesh::DIVERGED_DTOL
Definition: enum_convergence_flags.h:47
libMesh::PetscPreconditioner
This class provides an interface to the suite of preconditioners available from PETSc.
Definition: petsc_preconditioner.h:63
libMesh::DIVERGED_ITS
Definition: enum_convergence_flags.h:46
libMesh::Preconditioner::apply
virtual void apply(const NumericVector< T > &x, NumericVector< T > &y)=0
Computes the preconditioned vector y based on input vector x.
libMesh::BICG
Definition: enum_solver_type.h:41
libMesh::CR
Definition: enum_solver_type.h:37
libMesh::index_range
IntRange< std::size_t > index_range(const std::vector< T > &vec)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:106
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
Definition: factoryfunction.C:55
libMesh::PetscMatrix::mat
Mat mat()
Definition: petsc_matrix.h:281
libMesh::RICHARDSON
Definition: enum_solver_type.h:50
libMesh::MINRES
Definition: enum_solver_type.h:43
libMesh::Preconditioner::initialized
bool initialized() const
Definition: preconditioner.h:106
libMesh::DIVERGED_NULL
Definition: enum_convergence_flags.h:45
libMesh::ParallelObject::comm
const Parallel::Communicator & comm() const
Definition: parallel_object.h:94
libMesh::Preconditioner::setup
virtual void setup()
This is called every time the "operator might have changed".
Definition: preconditioner.h:132
libMesh::CONVERGED_CG_CONSTRAINED
Definition: enum_convergence_flags.h:41
libMesh::SUBSET_COPY_RHS
Definition: enum_subset_solve_mode.h:40
libMesh::ierr
ierr
Definition: petsc_dm_wrapper.C:72
libMesh::CONVERGED_ITERATING
Definition: enum_convergence_flags.h:56
libMesh::DIVERGED_NONSYMMETRIC
Definition: enum_convergence_flags.h:50
libMesh::SparseMatrix
Generic sparse matrix.
Definition: vector_fe_ex5.C:45
libMesh::SUBSET_ZERO
Definition: enum_subset_solve_mode.h:37
libMesh::NumericVector::size
virtual numeric_index_type size() const =0
PetscBool
PetscTruth PetscBool
Definition: petsc_macro.h:73
libMesh::NumericVector
Provides a uniform interface to vector storage schemes for different linear algebra libraries.
Definition: vector_fe_ex5.C:43
libMesh::CONVERGED_RTOL
Definition: enum_convergence_flags.h:37
libMesh::TriangleWrapper::init
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
libMesh::libmesh_petsc_preconditioner_apply
PetscErrorCode libmesh_petsc_preconditioner_apply(PC pc, Vec x, Vec y)
Definition: petsc_linear_solver.C:88
libMesh::CG
Definition: enum_solver_type.h:34
libMesh::LinearConvergenceReason
LinearConvergenceReason
Linear solver convergence flags (taken from the PETSc flags).
Definition: enum_convergence_flags.h:33
libMesh::TFQMR
Definition: enum_solver_type.h:40
libMesh::libmesh_petsc_preconditioner_setup
PetscErrorCode libmesh_petsc_preconditioner_setup(void *ctx)
This function is called by PETSc to initialize the preconditioner.
Definition: petsc_linear_solver.C:49
libMesh::CGS
Definition: enum_solver_type.h:36
libMesh::PetscVector::close
virtual void close() override
Calls the NumericVector's internal assembly routines, ensuring that the values are consistent across ...
Definition: petsc_vector.h:821
libMesh::ShellMatrix::vector_mult
virtual void vector_mult(NumericVector< T > &dest, const NumericVector< T > &arg) const =0
Multiplies the matrix with arg and stores the result in dest.
libMesh::__libmesh_petsc_preconditioner_setup
PetscErrorCode __libmesh_petsc_preconditioner_setup(void *ctx)
This function is called by PETSc to initialize the preconditioner.
Definition: petsc_linear_solver.C:105
libMesh::Preconditioner< Number >
libMesh::PetscLinearSolver
This class provides an interface to PETSc iterative solvers that is compatible with the libMesh Linea...
Definition: petsc_linear_solver.h:107
libMesh::TCQMR
Definition: enum_solver_type.h:39
libMesh::libmesh_petsc_preconditioner_setup
PetscErrorCode libmesh_petsc_preconditioner_setup(PC pc)
Definition: petsc_linear_solver.C:74
libMesh::BLOCK_JACOBI_PRECOND
Definition: enum_preconditioner_type.h:36
libMesh::PetscLinearSolver::PetscLinearSolver
PetscLinearSolver(const libMesh::Parallel::Communicator &comm_in)
Constructor.
Definition: petsc_linear_solver.C:135
libMesh::NumericVector::local_size
virtual numeric_index_type local_size() const =0
libMesh::PetscVector
This class provides a nice interface to PETSc's Vec object.
Definition: petsc_vector.h:72
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh::CONVERGED_HAPPY_BREAKDOWN
Definition: enum_convergence_flags.h:43
libMesh::GMRES
Definition: enum_solver_type.h:44
libMesh::CONVERGED_RTOL_NORMAL
Definition: enum_convergence_flags.h:35
libMesh::DIVERGED_INDEFINITE_PC
Definition: enum_convergence_flags.h:51
PETSC_OWN_POINTER
Definition: petsc_macro.h:103
libMesh::initialized
bool initialized()
Checks that library initialization has been done.
Definition: libmesh.C:265
libMesh::CONVERGED_CG_NEG_CURVE
Definition: enum_convergence_flags.h:40
libMesh::CONVERGED_ATOL_NORMAL
Definition: enum_convergence_flags.h:36
libMesh::ReferenceElem::get
const Elem & get(const ElemType type_in)
Definition: reference_elem.C:237
libMesh::SUBSET_DONT_TOUCH
Definition: enum_subset_solve_mode.h:45
libMesh::CONVERGED_ITS
Definition: enum_convergence_flags.h:39
libMesh::libMeshPrivateData::_is_initialized
bool _is_initialized
Flag that tells if init() has been called.
Definition: libmesh.C:246
libMesh::PetscMatrix
This class provides a nice interface to the PETSc C-based data structures for parallel,...
Definition: petsc_matrix.h:87
libMesh::LinearSolver
This base class can be inherited from to provide interfaces to linear solvers from different packages...
Definition: linear_solver.h:69
libMesh::CONVERGED_ATOL
Definition: enum_convergence_flags.h:38
libMesh::LSQR
Definition: enum_solver_type.h:45
libMesh::__libmesh_petsc_preconditioner_apply
PetscErrorCode __libmesh_petsc_preconditioner_apply(void *ctx, Vec x, Vec y)
This function is called by PETSc to actually apply the preconditioner.
Definition: petsc_linear_solver.C:111
libMesh::SubsetSolveMode
SubsetSolveMode
Definition: enum_subset_solve_mode.h:35
libMesh::DIVERGED_PCSETUP_FAILED
Definition: enum_convergence_flags.h:54
libMesh::DIVERGED_NAN
Definition: enum_convergence_flags.h:52
libMesh::ctx
void * ctx
Definition: petsc_dm_wrapper.C:71
libMesh::PetscMatrix::close
virtual void close() override
Calls the SparseMatrix's internal assembly routines, ensuring that the values are consistent across p...
Definition: petsc_matrix.C:989
libMesh::libmesh_petsc_preconditioner_apply
PetscErrorCode libmesh_petsc_preconditioner_apply(void *ctx, Vec x, Vec y)
This function is called by PETSc to actually apply the preconditioner.
Definition: petsc_linear_solver.C:62
libMesh::ShellMatrix::vector_mult_add
virtual void vector_mult_add(NumericVector< T > &dest, const NumericVector< T > &arg) const =0
Multiplies the matrix with arg and adds the result to dest.
libMesh::PetscVector::vec
Vec vec()
Definition: petsc_vector.h:335
libMesh::err
OStreamProxy err
libMesh::CHEBYSHEV
Definition: enum_solver_type.h:51
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:121
libMesh::CHKERRQ
CHKERRQ(ierr)
ierr
PetscErrorCode ierr
Definition: petscdmlibmeshimpl.C:1028
libMesh::petsc_auto_fieldsplit
void petsc_auto_fieldsplit(PC my_pc, const System &sys)
Definition: petsc_auto_fieldsplit.C:58
libMesh::Quality::name
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
libMesh::ILU_PRECOND
Definition: enum_preconditioner_type.h:43
Generated on Sat Jan 25 2020 12:07:04 for libMesh by   doxygen 1.8.16