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