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