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