Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 : // Local Includes
24 : #include "libmesh/dof_map.h"
25 : #include "libmesh/libmesh_logging.h"
26 : #include "libmesh/petsc_linear_solver.h"
27 : #include "libmesh/shell_matrix.h"
28 : #include "libmesh/petsc_matrix.h"
29 : #include "libmesh/petsc_preconditioner.h"
30 : #include "libmesh/petsc_vector.h"
31 : #include "libmesh/enum_to_string.h"
32 : #include "libmesh/system.h"
33 : #include "libmesh/petsc_auto_fieldsplit.h"
34 : #include "libmesh/solver_configuration.h"
35 : #include "libmesh/enum_preconditioner_type.h"
36 : #include "libmesh/enum_solver_type.h"
37 : #include "libmesh/enum_convergence_flags.h"
38 :
39 : #if defined(LIBMESH_HAVE_PETSC_HYPRE) && PETSC_VERSION_LESS_THAN(3, 23, 0) && \
40 : !PETSC_VERSION_LESS_THAN(3, 12, 0) && defined(PETSC_HAVE_HYPRE_DEVICE)
41 : #include <HYPRE_utilities.h>
42 : #endif
43 :
44 : // C++ includes
45 : #include <memory>
46 : #include <string.h>
47 :
48 : namespace libMesh
49 : {
50 :
51 : extern "C"
52 : {
53 490 : PetscErrorCode libmesh_petsc_preconditioner_setup (PC pc)
54 : {
55 : PetscFunctionBegin;
56 :
57 : void * ctx;
58 490 : LibmeshPetscCallQ(PCShellGetContext(pc,&ctx));
59 490 : Preconditioner<Number> * preconditioner = static_cast<Preconditioner<Number> *>(ctx);
60 :
61 490 : libmesh_error_msg_if(!preconditioner->initialized(),
62 : "Preconditioner not initialized! Make sure you call init() before solve!");
63 :
64 490 : preconditioner->setup();
65 :
66 490 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
67 : }
68 :
69 1330 : PetscErrorCode libmesh_petsc_preconditioner_apply(PC pc, Vec x, Vec y)
70 : {
71 : PetscFunctionBegin;
72 :
73 : void * ctx;
74 1330 : LibmeshPetscCallQ(PCShellGetContext(pc,&ctx));
75 1330 : Preconditioner<Number> * preconditioner = static_cast<Preconditioner<Number> *>(ctx);
76 :
77 1406 : PetscVector<Number> x_vec(x, preconditioner->comm());
78 1368 : PetscVector<Number> y_vec(y, preconditioner->comm());
79 :
80 1330 : preconditioner->apply(x_vec,y_vec);
81 :
82 1368 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
83 1254 : }
84 :
85 : #ifdef LIBMESH_ENABLE_DEPRECATED
86 0 : PetscErrorCode __libmesh_petsc_preconditioner_setup (PC pc)
87 : {
88 : PetscFunctionBegin;
89 : libmesh_deprecated();
90 0 : PetscFunctionReturn(libmesh_petsc_preconditioner_setup(pc));
91 : }
92 :
93 0 : PetscErrorCode __libmesh_petsc_preconditioner_apply(PC pc, Vec x, Vec y)
94 : {
95 : PetscFunctionBegin;
96 : libmesh_deprecated();
97 0 : PetscFunctionReturn(libmesh_petsc_preconditioner_apply(pc, x, y));
98 : }
99 : #endif
100 : } // end extern "C"
101 :
102 :
103 :
104 : template <typename T>
105 22822 : PetscLinearSolver<T>::PetscLinearSolver(const libMesh::Parallel::Communicator & comm_in) :
106 : LinearSolver<T>(comm_in),
107 : _restrict_solve_to_is(nullptr),
108 : _restrict_solve_to_is_complement(nullptr),
109 22822 : _subset_solve_mode(SUBSET_ZERO)
110 : {
111 23496 : if (this->n_processors() == 1)
112 329 : this->_preconditioner_type = ILU_PRECOND;
113 : else
114 22493 : this->_preconditioner_type = BLOCK_JACOBI_PRECOND;
115 22822 : }
116 :
117 :
118 :
119 : template <typename T>
120 36631 : void PetscLinearSolver<T>::clear ()
121 : {
122 36631 : if (this->initialized())
123 : {
124 24701 : this->_is_initialized = false;
125 :
126 : // Calls specialized destroy() functions
127 24701 : if (_restrict_solve_to_is)
128 210 : _restrict_solve_to_is.reset_to_zero();
129 24701 : if (_restrict_solve_to_is_complement)
130 0 : _restrict_solve_to_is_complement.reset_to_zero();
131 :
132 : // Previously we only called KSPDestroy(), we did not reset _ksp
133 : // to nullptr, so that behavior is maintained here.
134 24701 : _ksp.destroy();
135 :
136 : // Mimic PETSc default solver and preconditioner
137 24701 : this->_solver_type = GMRES;
138 :
139 24701 : if (!this->_preconditioner)
140 : {
141 25417 : if (this->n_processors() == 1)
142 347 : this->_preconditioner_type = ILU_PRECOND;
143 : else
144 24214 : this->_preconditioner_type = BLOCK_JACOBI_PRECOND;
145 : }
146 : }
147 36631 : }
148 :
149 :
150 :
151 : template <typename T>
152 595811 : void PetscLinearSolver<T>::init (const char * name)
153 : {
154 : // Initialize the data structures if not done so already.
155 595811 : if (!this->initialized())
156 : {
157 44861 : this->_is_initialized = true;
158 :
159 44861 : LibmeshPetscCall(KSPCreate(this->comm().get(), _ksp.get()));
160 :
161 44861 : if (name)
162 490 : LibmeshPetscCall(KSPSetOptionsPrefix(_ksp, name));
163 :
164 : // Create the preconditioner context
165 44861 : LibmeshPetscCall(KSPGetPC(_ksp, &_pc));
166 :
167 : // Set user-specified solver and preconditioner types
168 44861 : this->set_petsc_solver_type();
169 :
170 : // If the SolverConfiguration object is provided, use it to set
171 : // options during solver initialization.
172 44861 : if (this->_solver_configuration)
173 : {
174 70 : this->_solver_configuration->set_options_during_init();
175 : }
176 :
177 : // Set the options from user-input
178 : // Set runtime options, e.g.,
179 : // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
180 : // These options will override those specified above as long as
181 : // KSPSetFromOptions() is called _after_ any other customization
182 : // routines.
183 44861 : LibmeshPetscCall(KSPSetFromOptions(_ksp));
184 :
185 : // Have the Krylov subspace method use our good initial guess
186 : // rather than 0, unless the user requested a KSPType of
187 : // preonly, which complains if asked to use initial guesses.
188 : KSPType ksp_type;
189 :
190 44861 : LibmeshPetscCall(KSPGetType(_ksp, &ksp_type));
191 :
192 44861 : if (strcmp(ksp_type, "preonly"))
193 44861 : LibmeshPetscCall(KSPSetInitialGuessNonzero(_ksp, PETSC_TRUE));
194 :
195 : // Notify PETSc of location to store residual history.
196 : // This needs to be called before any solves, since
197 : // it sets the residual history length to zero. The default
198 : // behavior is for PETSc to allocate (internally) an array
199 : // of size 1000 to hold the residual norm history.
200 44861 : LibmeshPetscCall(KSPSetResidualHistory(_ksp,
201 : LIBMESH_PETSC_NULLPTR, // pointer to the array which holds the history
202 : PETSC_DECIDE, // size of the array holding the history
203 : PETSC_TRUE)); // Whether or not to reset the history for each solve.
204 :
205 44861 : PetscPreconditioner<T>::set_petsc_preconditioner_type(this->_preconditioner_type,_pc);
206 :
207 : //If there is a preconditioner object we need to set the internal setup and apply routines
208 44861 : if (this->_preconditioner)
209 : {
210 210 : this->_preconditioner->init();
211 210 : LibmeshPetscCall(PCShellSetContext(_pc,(void *)this->_preconditioner));
212 210 : LibmeshPetscCall(PCShellSetSetUp(_pc,libmesh_petsc_preconditioner_setup));
213 210 : LibmeshPetscCall(PCShellSetApply(_pc,libmesh_petsc_preconditioner_apply));
214 : }
215 : }
216 595811 : }
217 :
218 :
219 : template <typename T>
220 427711 : void PetscLinearSolver<T>::init (PetscMatrixBase<T> * matrix,
221 : const char * name)
222 : {
223 : // Initialize the data structures if not done so already.
224 427711 : if (!this->initialized())
225 : {
226 786 : if (this->_preconditioner)
227 0 : this->_preconditioner->set_matrix(*matrix);
228 :
229 786 : this->init(name);
230 :
231 : // Set operators. The input matrix works as the preconditioning matrix
232 786 : LibmeshPetscCall(KSPSetOperators(_ksp, matrix->mat(), matrix->mat()));
233 : }
234 427711 : }
235 :
236 :
237 :
238 : template <typename T>
239 : void
240 219516 : PetscLinearSolver<T>::init_names (const System & sys)
241 : {
242 219516 : petsc_auto_fieldsplit(this->pc(), sys);
243 219516 : }
244 :
245 :
246 :
247 : template <typename T>
248 : void
249 420 : PetscLinearSolver<T>::restrict_solve_to (const std::vector<unsigned int> * const dofs,
250 : const SubsetSolveMode subset_solve_mode)
251 : {
252 : // The preconditioner (in particular if a default preconditioner)
253 : // will have to be reset. We call this->clear() to do that. This
254 : // call will also remove and free any previous subset that may have
255 : // been set before.
256 420 : this->clear();
257 :
258 420 : _subset_solve_mode = subset_solve_mode;
259 :
260 420 : if (dofs != nullptr)
261 : {
262 210 : PetscInt * petsc_dofs = nullptr;
263 216 : LibmeshPetscCall(PetscMalloc(dofs->size()*sizeof(PetscInt), &petsc_dofs));
264 :
265 65330 : for (auto i : index_range(*dofs))
266 65120 : petsc_dofs[i] = (*dofs)[i];
267 :
268 : // Create the IS
269 : // PETSc now takes over ownership of the "petsc_dofs"
270 : // array, so we don't have to worry about it any longer.
271 216 : LibmeshPetscCall(ISCreateGeneral(this->comm().get(),
272 : cast_int<PetscInt>(dofs->size()),
273 : petsc_dofs, PETSC_OWN_POINTER,
274 : _restrict_solve_to_is.get()));
275 : }
276 420 : }
277 :
278 :
279 :
280 : template <typename T>
281 : std::pair<unsigned int, Real>
282 403041 : PetscLinearSolver<T>::solve (SparseMatrix<T> & matrix_in,
283 : SparseMatrix<T> & precond_in,
284 : NumericVector<T> & solution_in,
285 : NumericVector<T> & rhs_in,
286 : const std::optional<double> tol,
287 : const std::optional<unsigned int> m_its)
288 : {
289 23584 : LOG_SCOPE("solve()", "PetscLinearSolver");
290 :
291 806082 : const double rel_tol = this->get_real_solver_setting("rel_tol", tol);
292 794290 : const double abs_tol = this->get_real_solver_setting("abs_tol",
293 : std::nullopt,
294 : static_cast<double>(PETSC_DEFAULT));
295 403041 : const double max_its = this->get_int_solver_setting("max_its", m_its);
296 :
297 379457 : return this->solve_common(matrix_in, precond_in, solution_in,
298 426625 : rhs_in, rel_tol, abs_tol, max_its, KSPSolve);
299 : }
300 :
301 : template <typename T>
302 : std::pair<unsigned int, Real>
303 24670 : PetscLinearSolver<T>::adjoint_solve (SparseMatrix<T> & matrix_in,
304 : NumericVector<T> & solution_in,
305 : NumericVector<T> & rhs_in,
306 : const std::optional<double> tol,
307 : const std::optional<unsigned int> m_its)
308 : {
309 1596 : LOG_SCOPE("adjoint_solve()", "PetscLinearSolver");
310 :
311 49340 : const double rel_tol = this->get_real_solver_setting("rel_tol", tol);
312 48542 : const double abs_tol = this->get_real_solver_setting("abs_tol",
313 : std::nullopt,
314 : static_cast<double>(PETSC_DEFAULT));
315 24670 : const double max_its = this->get_int_solver_setting("max_its", m_its);
316 :
317 : // Note that the matrix and precond matrix are the same
318 23074 : return this->solve_common(matrix_in, matrix_in, solution_in,
319 26266 : rhs_in, rel_tol, abs_tol, max_its, KSPSolveTranspose);
320 : }
321 :
322 :
323 : template <typename T>
324 : std::pair<unsigned int, Real>
325 427711 : PetscLinearSolver<T>::solve_common (SparseMatrix<T> & matrix_in,
326 : SparseMatrix<T> & precond_in,
327 : NumericVector<T> & solution_in,
328 : NumericVector<T> & rhs_in,
329 : const double rel_tol,
330 : const double abs_tol,
331 : const unsigned int m_its,
332 : ksp_solve_func_type solve_func)
333 : {
334 : // Make sure the data passed in are really of Petsc types
335 12590 : PetscMatrixBase<T> * matrix = cast_ptr<PetscMatrixBase<T> *>(&matrix_in);
336 12590 : PetscMatrixBase<T> * precond = cast_ptr<PetscMatrixBase<T> *>(&precond_in);
337 :
338 427711 : this->init (matrix);
339 :
340 : // Close the matrices in case this wasn't already done.
341 427711 : matrix->close ();
342 427711 : if (precond != matrix)
343 0 : precond->close ();
344 :
345 25180 : auto mat = matrix->mat();
346 :
347 : return this->solve_base
348 427711 : (matrix, precond, mat, solution_in, rhs_in, rel_tol, abs_tol, m_its, solve_func);
349 : }
350 :
351 :
352 : template <typename T>
353 : std::pair<unsigned int, Real>
354 0 : PetscLinearSolver<T>::solve (const ShellMatrix<T> & shell_matrix,
355 : NumericVector<T> & solution_in,
356 : NumericVector<T> & rhs_in,
357 : const std::optional<double> tol,
358 : const std::optional<unsigned int> m_its)
359 : {
360 0 : LOG_SCOPE("solve()", "PetscLinearSolver");
361 :
362 0 : const double rel_tol = this->get_real_solver_setting("rel_tol", tol);
363 0 : const double abs_tol = this->get_real_solver_setting("abs_tol",
364 : std::nullopt,
365 : static_cast<double>(PETSC_DEFAULT));
366 0 : const double max_its = this->get_int_solver_setting("max_its", m_its);
367 :
368 0 : return this->shell_solve_common(shell_matrix, nullptr, solution_in,
369 0 : rhs_in, rel_tol, abs_tol, max_its);
370 : }
371 :
372 :
373 :
374 : template <typename T>
375 : std::pair<unsigned int, Real>
376 70 : PetscLinearSolver<T>::solve (const ShellMatrix<T> & shell_matrix,
377 : const SparseMatrix<T> & precond_matrix,
378 : NumericVector<T> & solution_in,
379 : NumericVector<T> & rhs_in,
380 : const std::optional<double> tol,
381 : const std::optional<unsigned int> m_its)
382 : {
383 140 : const double rel_tol = this->get_real_solver_setting("rel_tol", tol);
384 138 : const double abs_tol = this->get_real_solver_setting("abs_tol",
385 : std::nullopt,
386 : static_cast<double>(PETSC_DEFAULT));
387 70 : const double max_its = this->get_int_solver_setting("max_its", m_its);
388 :
389 : // Make sure the data passed in are really of Petsc types
390 2 : const PetscMatrixBase<T> * precond = cast_ptr<const PetscMatrixBase<T> *>(&precond_matrix);
391 :
392 : return this->shell_solve_common
393 66 : (shell_matrix, const_cast<PetscMatrixBase<T> *>(precond), solution_in,
394 70 : rhs_in, rel_tol, abs_tol, max_its);
395 : }
396 :
397 :
398 :
399 : template <typename T>
400 : std::pair<unsigned int, Real>
401 70 : PetscLinearSolver<T>::shell_solve_common (const ShellMatrix<T> & shell_matrix,
402 : PetscMatrixBase<T> * precond,
403 : NumericVector<T> & solution_in,
404 : NumericVector<T> & rhs_in,
405 : const double rel_tol,
406 : const double abs_tol,
407 : const unsigned int m_its)
408 : {
409 4 : LOG_SCOPE("solve()", "PetscLinearSolver");
410 :
411 : // We don't really have a matrix for our matrix here
412 2 : SparseMatrix<T> * matrix = nullptr;
413 :
414 70 : this->init ();
415 :
416 : // Prepare the matrix.
417 4 : WrappedPetsc<Mat> mat;
418 70 : LibmeshPetscCall(MatCreateShell(this->comm().get(),
419 : rhs_in.local_size(),
420 : solution_in.local_size(),
421 : rhs_in.size(),
422 : solution_in.size(),
423 : const_cast<void *>(static_cast<const void *>(&shell_matrix)),
424 : mat.get()));
425 : // Note that the const_cast above is only necessary because PETSc
426 : // does not accept a const void *. Inside the member function
427 : // _petsc_shell_matrix() below, the pointer is casted back to a
428 : // const ShellMatrix<T> *.
429 :
430 70 : LibmeshPetscCall(MatShellSetOperation(mat, MATOP_MULT, reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult)));
431 70 : LibmeshPetscCall(MatShellSetOperation(mat, MATOP_MULT_ADD, reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult_add)));
432 70 : LibmeshPetscCall(MatShellSetOperation(mat, MATOP_GET_DIAGONAL, reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal)));
433 :
434 : return this->solve_base
435 140 : (matrix, precond, mat, solution_in, rhs_in, rel_tol, abs_tol, m_its, KSPSolve);
436 : }
437 :
438 :
439 :
440 : template <typename T>
441 : std::pair<unsigned int, Real>
442 427781 : PetscLinearSolver<T>::solve_base (SparseMatrix<T> * matrix,
443 : PetscMatrixBase<T> * precond,
444 : Mat mat,
445 : NumericVector<T> & solution_in,
446 : NumericVector<T> & rhs_in,
447 : const double rel_tol,
448 : const double abs_tol,
449 : const unsigned int m_its,
450 : ksp_solve_func_type solve_func)
451 : {
452 : // Make sure the data passed in are really of Petsc types
453 12592 : PetscVector<T> * solution = cast_ptr<PetscVector<T> *>(&solution_in);
454 12592 : PetscVector<T> * rhs = cast_ptr<PetscVector<T> *>(&rhs_in);
455 :
456 : // Close the vectors in case this wasn't already done.
457 427781 : solution->close ();
458 427781 : rhs->close ();
459 :
460 427781 : PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
461 427781 : PetscReal final_resid=0.;
462 :
463 415189 : std::unique_ptr<PetscMatrixBase<Number>> subprecond_matrix;
464 25184 : WrappedPetsc<Mat> subprecond;
465 :
466 25184 : WrappedPetsc<Mat> submat;
467 25184 : WrappedPetsc<Vec> subrhs;
468 25184 : WrappedPetsc<Vec> subsolution;
469 12592 : WrappedPetsc<VecScatter> scatter;
470 :
471 : // Restrict rhs and solution vectors and set operators. The input
472 : // matrix works as the preconditioning matrix.
473 427781 : if (_restrict_solve_to_is)
474 : {
475 210 : PetscInt is_local_size = this->restrict_solve_to_is_local_size();
476 :
477 210 : LibmeshPetscCall(VecCreate(this->comm().get(), subrhs.get()));
478 210 : LibmeshPetscCall(VecSetSizes(subrhs, is_local_size, PETSC_DECIDE));
479 210 : LibmeshPetscCall(VecSetFromOptions(subrhs));
480 :
481 210 : LibmeshPetscCall(VecCreate(this->comm().get(), subsolution.get()));
482 210 : LibmeshPetscCall(VecSetSizes(subsolution, is_local_size, PETSC_DECIDE));
483 210 : LibmeshPetscCall(VecSetFromOptions(subsolution));
484 :
485 210 : LibmeshPetscCall(VecScatterCreate(rhs->vec(), _restrict_solve_to_is, subrhs, nullptr, scatter.get()));
486 :
487 210 : VecScatterBeginEnd(this->comm(), scatter, rhs->vec(), subrhs, INSERT_VALUES, SCATTER_FORWARD);
488 210 : VecScatterBeginEnd(this->comm(), scatter, solution->vec(), subsolution, INSERT_VALUES, SCATTER_FORWARD);
489 :
490 210 : LibmeshPetscCall(LibMeshCreateSubMatrix(mat,
491 : _restrict_solve_to_is,
492 : _restrict_solve_to_is,
493 : MAT_INITIAL_MATRIX,
494 : submat.get()));
495 :
496 210 : if (precond)
497 210 : LibmeshPetscCall(LibMeshCreateSubMatrix(const_cast<PetscMatrixBase<T> *>(precond)->mat(),
498 : _restrict_solve_to_is,
499 : _restrict_solve_to_is,
500 : MAT_INITIAL_MATRIX,
501 : subprecond.get()));
502 :
503 : // Since removing columns of the matrix changes the equation
504 : // system, we will now change the right hand side to compensate
505 : // for this. Note that this is not necessary if \p SUBSET_ZERO
506 : // has been selected.
507 210 : if (_subset_solve_mode!=SUBSET_ZERO)
508 : {
509 0 : this->create_complement_is(rhs_in);
510 : PetscInt is_complement_local_size =
511 0 : cast_int<PetscInt>(rhs_in.local_size()-is_local_size);
512 :
513 0 : WrappedPetsc<Vec> subvec1;
514 0 : LibmeshPetscCall(VecCreate(this->comm().get(), subvec1.get()));
515 0 : LibmeshPetscCall(VecSetSizes(subvec1, is_complement_local_size, PETSC_DECIDE));
516 0 : LibmeshPetscCall(VecSetFromOptions(subvec1));
517 :
518 0 : WrappedPetsc<VecScatter> scatter1;
519 0 : LibmeshPetscCall(VecScatterCreate(rhs->vec(), _restrict_solve_to_is_complement, subvec1, nullptr, scatter1.get()));
520 :
521 0 : VecScatterBeginEnd(this->comm(), scatter1, _subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(), subvec1, INSERT_VALUES, SCATTER_FORWARD);
522 :
523 0 : LibmeshPetscCall(VecScale(subvec1, -1.0));
524 :
525 0 : WrappedPetsc<Mat> submat1;
526 0 : LibmeshPetscCall(LibMeshCreateSubMatrix(mat,
527 : _restrict_solve_to_is,
528 : _restrict_solve_to_is_complement,
529 : MAT_INITIAL_MATRIX,
530 : submat1.get()));
531 :
532 0 : LibmeshPetscCall(MatMultAdd(submat1, subvec1, subrhs, subrhs));
533 : }
534 :
535 210 : if (precond)
536 210 : LibmeshPetscCall(KSPSetOperators(_ksp, submat, subprecond));
537 : else
538 0 : LibmeshPetscCall(KSPSetOperators(_ksp, submat, submat));
539 :
540 210 : PetscBool ksp_reuse_preconditioner = this->same_preconditioner ? PETSC_TRUE : PETSC_FALSE;
541 210 : LibmeshPetscCall(KSPSetReusePreconditioner(_ksp, ksp_reuse_preconditioner));
542 :
543 210 : if (precond && this->_preconditioner)
544 : {
545 0 : subprecond_matrix = std::make_unique<PetscMatrix<Number>>(subprecond, this->comm());
546 0 : this->_preconditioner->set_matrix(*subprecond_matrix);
547 0 : this->_preconditioner->init();
548 : }
549 : }
550 : else
551 : {
552 427571 : PetscBool ksp_reuse_preconditioner = this->same_preconditioner ? PETSC_TRUE : PETSC_FALSE;
553 427571 : LibmeshPetscCall(KSPSetReusePreconditioner(_ksp, ksp_reuse_preconditioner));
554 :
555 427571 : if (precond)
556 427571 : LibmeshPetscCall(KSPSetOperators(_ksp, mat, const_cast<PetscMatrixBase<T> *>(precond)->mat()));
557 : else
558 0 : LibmeshPetscCall(KSPSetOperators(_ksp, mat, mat));
559 :
560 427571 : if (this->_preconditioner)
561 : {
562 210 : if (matrix)
563 6 : this->_preconditioner->set_matrix(*matrix);
564 0 : else if (precond)
565 0 : this->_preconditioner->set_matrix(const_cast<PetscMatrixBase<Number> &>(*precond));
566 :
567 210 : this->_preconditioner->init();
568 : }
569 : }
570 :
571 : // Set the tolerances for the iterative solver. Use the user-supplied
572 : // tolerance for the relative residual & leave the others at default values.
573 427781 : LibmeshPetscCall(KSPSetTolerances(_ksp, rel_tol, abs_tol,
574 : PETSC_DEFAULT, max_its));
575 :
576 : // Allow command line options to override anything set programmatically.
577 427781 : LibmeshPetscCall(KSPSetFromOptions(_ksp));
578 :
579 : #if defined(LIBMESH_HAVE_PETSC_HYPRE) && PETSC_VERSION_LESS_THAN(3, 23, 0) && \
580 : !PETSC_VERSION_LESS_THAN(3, 12, 0) && defined(PETSC_HAVE_HYPRE_DEVICE)
581 : {
582 : // Make sure hypre has been initialized
583 : LibmeshPetscCallExternal(HYPRE_Initialize);
584 : PetscScalar * dummyarray;
585 : PetscMemType mtype;
586 : LibmeshPetscCall(VecGetArrayAndMemType(solution->vec(), &dummyarray, &mtype));
587 : LibmeshPetscCall(VecRestoreArrayAndMemType(solution->vec(), &dummyarray));
588 : if (PetscMemTypeHost(mtype))
589 : LibmeshPetscCallExternal(HYPRE_SetMemoryLocation, HYPRE_MEMORY_HOST);
590 : }
591 : #endif
592 :
593 : // If the SolverConfiguration object is provided, use it to override
594 : // solver options.
595 427781 : if (this->_solver_configuration)
596 : {
597 70 : this->_solver_configuration->configure_solver();
598 : }
599 :
600 : // Solve the linear system
601 427781 : if (_restrict_solve_to_is)
602 210 : LibmeshPetscCall(solve_func (_ksp, subrhs, subsolution));
603 : else
604 427571 : LibmeshPetscCall(solve_func (_ksp, rhs->vec(), solution->vec()));
605 :
606 : // Get the number of iterations required for convergence
607 427781 : LibmeshPetscCall(KSPGetIterationNumber (_ksp, &its));
608 :
609 : // Get the norm of the final residual to return to the user.
610 427781 : LibmeshPetscCall(KSPGetResidualNorm (_ksp, &final_resid));
611 :
612 427781 : if (_restrict_solve_to_is)
613 : {
614 210 : switch(_subset_solve_mode)
615 : {
616 12 : case SUBSET_ZERO:
617 210 : LibmeshPetscCall(VecZeroEntries(solution->vec()));
618 6 : break;
619 :
620 0 : case SUBSET_COPY_RHS:
621 0 : LibmeshPetscCall(VecCopy(rhs->vec(),solution->vec()));
622 0 : break;
623 :
624 0 : case SUBSET_DONT_TOUCH:
625 : // Nothing to do here.
626 0 : break;
627 :
628 0 : default:
629 0 : libmesh_error_msg("Invalid subset solve mode = " << _subset_solve_mode);
630 : }
631 :
632 210 : VecScatterBeginEnd(this->comm(), scatter, subsolution, solution->vec(), INSERT_VALUES, SCATTER_REVERSE);
633 :
634 210 : if (precond && this->_preconditioner)
635 : {
636 : // Before subprecond_matrix gets cleaned up, we should give
637 : // the _preconditioner a different matrix.
638 0 : if (matrix)
639 0 : this->_preconditioner->set_matrix(*matrix);
640 : else
641 0 : this->_preconditioner->set_matrix(*precond);
642 :
643 0 : this->_preconditioner->init();
644 : }
645 : }
646 :
647 : // return the # of its. and the final residual norm.
648 881586 : return std::make_pair(its, final_resid);
649 402597 : }
650 :
651 :
652 :
653 : template <typename T>
654 280 : KSP PetscLinearSolver<T>::ksp()
655 : {
656 280 : this->init();
657 280 : return _ksp;
658 : }
659 :
660 : template <typename T>
661 0 : void PetscLinearSolver<T>::get_residual_history(std::vector<double> & hist)
662 : {
663 0 : PetscInt its = 0;
664 :
665 : // Fill the residual history vector with the residual norms
666 : // Note that GetResidualHistory() does not copy any values, it
667 : // simply sets the pointer p. Note that for some Krylov subspace
668 : // methods, the number of residuals returned in the history
669 : // vector may be different from what you are expecting. For
670 : // example, TFQMR returns two residual values per iteration step.
671 :
672 : // Recent versions of PETSc require the residual
673 : // history vector pointer to be declared as const.
674 : #if PETSC_VERSION_LESS_THAN(3,15,0)
675 : PetscReal * p;
676 : #else
677 : const PetscReal * p;
678 : #endif
679 :
680 0 : LibmeshPetscCall(KSPGetResidualHistory(_ksp, &p, &its));
681 :
682 : // Check for early return
683 0 : if (its == 0) return;
684 :
685 : // Create space to store the result
686 0 : hist.resize(its);
687 :
688 : // Copy history into the vector provided by the user.
689 0 : for (PetscInt i=0; i<its; ++i)
690 : {
691 0 : hist[i] = *p;
692 0 : p++;
693 : }
694 : }
695 :
696 :
697 :
698 :
699 : template <typename T>
700 0 : Real PetscLinearSolver<T>::get_initial_residual()
701 : {
702 0 : PetscInt its = 0;
703 :
704 : // Fill the residual history vector with the residual norms
705 : // Note that GetResidualHistory() does not copy any values, it
706 : // simply sets the pointer p. Note that for some Krylov subspace
707 : // methods, the number of residuals returned in the history
708 : // vector may be different from what you are expecting. For
709 : // example, TFQMR returns two residual values per iteration step.
710 :
711 : // Recent versions of PETSc require the residual
712 : // history vector pointer to be declared as const.
713 : #if PETSC_VERSION_LESS_THAN(3,15,0)
714 : PetscReal * p;
715 : #else
716 : const PetscReal * p;
717 : #endif
718 :
719 :
720 0 : LibmeshPetscCall(KSPGetResidualHistory(_ksp, &p, &its));
721 :
722 : // Check no residual history
723 0 : if (its == 0)
724 : {
725 0 : libMesh::err << "No iterations have been performed, returning 0." << std::endl;
726 0 : return 0.;
727 : }
728 :
729 : // Otherwise, return the value pointed to by p.
730 0 : return *p;
731 : }
732 :
733 :
734 :
735 :
736 : template <typename T>
737 44861 : void PetscLinearSolver<T>::set_petsc_solver_type()
738 : {
739 : #define CaseKSPSetType(SolverType, KSPT) \
740 : case SolverType: \
741 : LibmeshPetscCall(KSPSetType (_ksp, const_cast<KSPType>(KSPT))); \
742 : return;
743 :
744 44861 : switch (this->_solver_type)
745 : {
746 140 : CaseKSPSetType(CG, KSPCG)
747 0 : CaseKSPSetType(CR, KSPCR)
748 0 : CaseKSPSetType(CGS, KSPCGS)
749 0 : CaseKSPSetType(BICG, KSPBICG)
750 0 : CaseKSPSetType(TCQMR, KSPTCQMR)
751 0 : CaseKSPSetType(TFQMR, KSPTFQMR)
752 0 : CaseKSPSetType(LSQR, KSPLSQR)
753 0 : CaseKSPSetType(BICGSTAB, KSPBCGS)
754 0 : CaseKSPSetType(MINRES, KSPMINRES)
755 44651 : CaseKSPSetType(GMRES, KSPGMRES)
756 0 : CaseKSPSetType(RICHARDSON, KSPRICHARDSON)
757 0 : CaseKSPSetType(CHEBYSHEV, KSPCHEBYSHEV)
758 :
759 2 : default:
760 2 : libMesh::err << "ERROR: Unsupported PETSC Solver: "
761 140 : << Utility::enum_to_string(this->_solver_type) << std::endl
762 2 : << "Continuing with PETSC defaults" << std::endl;
763 : }
764 : }
765 :
766 :
767 :
768 : template <typename T>
769 160259 : LinearConvergenceReason PetscLinearSolver<T>::get_converged_reason() const
770 : {
771 : KSPConvergedReason reason;
772 160259 : LibmeshPetscCall(KSPGetConvergedReason(_ksp, &reason));
773 :
774 160259 : switch(reason)
775 : {
776 0 : case KSP_CONVERGED_RTOL_NORMAL : return CONVERGED_RTOL_NORMAL;
777 0 : case KSP_CONVERGED_ATOL_NORMAL : return CONVERGED_ATOL_NORMAL;
778 159839 : case KSP_CONVERGED_RTOL : return CONVERGED_RTOL;
779 420 : case KSP_CONVERGED_ATOL : return CONVERGED_ATOL;
780 0 : case KSP_CONVERGED_ITS : return CONVERGED_ITS;
781 : #if PETSC_VERSION_LESS_THAN(3,19,0)
782 0 : case KSP_CONVERGED_CG_NEG_CURVE : return CONVERGED_CG_NEG_CURVE;
783 : // This was deprecated for STEP_LENGTH
784 0 : case KSP_CONVERGED_CG_CONSTRAINED : return CONVERGED_CG_CONSTRAINED;
785 : #else
786 0 : case KSP_CONVERGED_NEG_CURVE : return CONVERGED_CG_NEG_CURVE;
787 : #endif
788 0 : case KSP_CONVERGED_STEP_LENGTH : return CONVERGED_STEP_LENGTH;
789 0 : case KSP_CONVERGED_HAPPY_BREAKDOWN : return CONVERGED_HAPPY_BREAKDOWN;
790 0 : case KSP_DIVERGED_NULL : return DIVERGED_NULL;
791 0 : case KSP_DIVERGED_ITS : return DIVERGED_ITS;
792 0 : case KSP_DIVERGED_DTOL : return DIVERGED_DTOL;
793 0 : case KSP_DIVERGED_BREAKDOWN : return DIVERGED_BREAKDOWN;
794 0 : case KSP_DIVERGED_BREAKDOWN_BICG : return DIVERGED_BREAKDOWN_BICG;
795 0 : case KSP_DIVERGED_NONSYMMETRIC : return DIVERGED_NONSYMMETRIC;
796 0 : case KSP_DIVERGED_INDEFINITE_PC : return DIVERGED_INDEFINITE_PC;
797 0 : case KSP_DIVERGED_NANORINF : return DIVERGED_NAN;
798 0 : case KSP_DIVERGED_INDEFINITE_MAT : return DIVERGED_INDEFINITE_MAT;
799 0 : case KSP_CONVERGED_ITERATING : return CONVERGED_ITERATING;
800 : #if !PETSC_VERSION_LESS_THAN(3,7,0)
801 : // PETSc-3.7.0 to 3.10.x
802 : #if PETSC_VERSION_LESS_THAN(3,11,0) && PETSC_VERSION_RELEASE
803 : case KSP_DIVERGED_PCSETUP_FAILED : return DIVERGED_PCSETUP_FAILED;
804 : // PETSc master and future PETSc
805 : #else
806 0 : case KSP_DIVERGED_PC_FAILED : return DIVERGED_PCSETUP_FAILED;
807 : #endif
808 : #endif
809 0 : default :
810 0 : libMesh::err << "Unknown convergence flag!" << std::endl;
811 0 : return UNKNOWN_FLAG;
812 : }
813 : }
814 :
815 :
816 : template <typename T>
817 988 : PetscErrorCode PetscLinearSolver<T>::_petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest)
818 : {
819 : PetscFunctionBegin;
820 :
821 : // Get the matrix context.
822 : void * ctx;
823 988 : PetscErrorCode ierr = MatShellGetContext(mat,&ctx);
824 :
825 : // Get user shell matrix object.
826 988 : const ShellMatrix<T> & shell_matrix = *static_cast<const ShellMatrix<T> *>(ctx);
827 988 : CHKERRABORT(shell_matrix.comm().get(), ierr);
828 :
829 : // Make \p NumericVector instances around the vectors.
830 1152 : PetscVector<T> arg_global(arg, shell_matrix.comm());
831 1070 : PetscVector<T> dest_global(dest, shell_matrix.comm());
832 :
833 : // Call the user function.
834 988 : shell_matrix.vector_mult(dest_global,arg_global);
835 :
836 1070 : PetscFunctionReturn(ierr);
837 824 : }
838 :
839 :
840 :
841 : template <typename T>
842 0 : PetscErrorCode PetscLinearSolver<T>::_petsc_shell_matrix_mult_add(Mat mat, Vec arg, Vec add, Vec dest)
843 : {
844 : PetscFunctionBegin;
845 :
846 : // Get the matrix context.
847 : void * ctx;
848 0 : PetscErrorCode ierr = MatShellGetContext(mat,&ctx);
849 :
850 : // Get user shell matrix object.
851 0 : const ShellMatrix<T> & shell_matrix = *static_cast<const ShellMatrix<T> *>(ctx);
852 0 : CHKERRABORT(shell_matrix.comm().get(), ierr);
853 :
854 : // Make \p NumericVector instances around the vectors.
855 0 : PetscVector<T> arg_global(arg, shell_matrix.comm());
856 0 : PetscVector<T> dest_global(dest, shell_matrix.comm());
857 0 : PetscVector<T> add_global(add, shell_matrix.comm());
858 :
859 0 : if (add!=arg)
860 : {
861 0 : arg_global = add_global;
862 : }
863 :
864 : // Call the user function.
865 0 : shell_matrix.vector_mult_add(dest_global,arg_global);
866 :
867 0 : PetscFunctionReturn(ierr);
868 0 : }
869 :
870 :
871 :
872 : template <typename T>
873 0 : PetscErrorCode PetscLinearSolver<T>::_petsc_shell_matrix_get_diagonal(Mat mat, Vec dest)
874 : {
875 : PetscFunctionBegin;
876 :
877 : // Get the matrix context.
878 : void * ctx;
879 0 : PetscErrorCode ierr = MatShellGetContext(mat,&ctx);
880 :
881 : // Get user shell matrix object.
882 0 : const ShellMatrix<T> & shell_matrix = *static_cast<const ShellMatrix<T> *>(ctx);
883 0 : CHKERRABORT(shell_matrix.comm().get(), ierr);
884 :
885 : // Make \p NumericVector instances around the vector.
886 0 : PetscVector<T> dest_global(dest, shell_matrix.comm());
887 :
888 : // Call the user function.
889 0 : shell_matrix.get_diagonal(dest_global);
890 :
891 0 : PetscFunctionReturn(ierr);
892 0 : }
893 :
894 :
895 :
896 : template <typename T>
897 : void
898 0 : PetscLinearSolver<T>::create_complement_is (const NumericVector<T> & vec_in)
899 : {
900 0 : libmesh_assert(_restrict_solve_to_is);
901 0 : if (!_restrict_solve_to_is_complement)
902 0 : LibmeshPetscCall(ISComplement(_restrict_solve_to_is,
903 : vec_in.first_local_index(),
904 : vec_in.last_local_index(),
905 : _restrict_solve_to_is_complement.get()));
906 0 : }
907 :
908 :
909 : template <typename T>
910 : PetscInt
911 210 : PetscLinearSolver<T>::restrict_solve_to_is_local_size() const
912 : {
913 6 : libmesh_assert(_restrict_solve_to_is);
914 :
915 : PetscInt s;
916 210 : LibmeshPetscCall(ISGetLocalSize(_restrict_solve_to_is, &s));
917 :
918 210 : return s;
919 : }
920 :
921 :
922 : //------------------------------------------------------------------
923 : // Explicit instantiations
924 : template class LIBMESH_EXPORT PetscLinearSolver<Number>;
925 :
926 : } // namespace libMesh
927 :
928 :
929 :
930 : #endif // #ifdef LIBMESH_HAVE_PETSC
|
