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 :
19 :
20 : #ifndef LIBMESH_PETSC_LINEAR_SOLVER_H
21 : #define LIBMESH_PETSC_LINEAR_SOLVER_H
22 :
23 : #include "libmesh/libmesh_config.h"
24 :
25 : #ifdef LIBMESH_HAVE_PETSC
26 :
27 : #include "libmesh/petsc_macro.h"
28 : #include "libmesh/petsc_solver_exception.h"
29 : #include "libmesh/wrapped_petsc.h"
30 :
31 : // Petsc include files.
32 : #ifdef I
33 : # define LIBMESH_SAW_I
34 : #endif
35 : #include <petscksp.h>
36 : #ifndef LIBMESH_SAW_I
37 : # undef I // Avoid complex.h contamination
38 : #endif
39 :
40 : // Local includes
41 : #include "libmesh/linear_solver.h"
42 :
43 : // C++ includes
44 : #include <cstddef>
45 : #include <vector>
46 : #include <memory>
47 :
48 : //--------------------------------------------------------------------
49 : // Functions with C linkage to pass to PETSc. PETSc will call these
50 : // methods as needed for preconditioning
51 : //
52 : // Since they must have C linkage they have no knowledge of a namespace.
53 : // Give them an obscure name to avoid namespace pollution.
54 : extern "C"
55 : {
56 : /**
57 : * This function is called by PETSc to initialize the preconditioner.
58 : * ctx will hold the Preconditioner.
59 : */
60 : PetscErrorCode libmesh_petsc_preconditioner_setup (PC);
61 :
62 : /**
63 : * This function is called by PETSc to actually apply the preconditioner.
64 : * ctx will hold the Preconditioner.
65 : */
66 : PetscErrorCode libmesh_petsc_preconditioner_apply(PC, Vec x, Vec y);
67 :
68 : #if LIBMESH_ENABLE_DEPRECATED
69 : /**
70 : * This function is called by PETSc to initialize the preconditioner.
71 : * ctx will hold the Preconditioner.
72 : */
73 : PetscErrorCode __libmesh_petsc_preconditioner_setup (PC);
74 :
75 : /**
76 : * This function is called by PETSc to actually apply the preconditioner.
77 : * ctx will hold the Preconditioner.
78 : */
79 : PetscErrorCode __libmesh_petsc_preconditioner_apply(PC, Vec x, Vec y);
80 : #endif
81 : } // end extern "C"
82 :
83 :
84 : namespace libMesh
85 : {
86 :
87 : // forward declarations
88 : template <typename T> class PetscMatrixBase;
89 :
90 : /**
91 : * This class provides an interface to PETSc
92 : * iterative solvers that is compatible with the \p libMesh
93 : * \p LinearSolver<>
94 : *
95 : * \author Benjamin Kirk
96 : * \date 2002-2007
97 : */
98 : template <typename T>
99 : class PetscLinearSolver : public LinearSolver<T>
100 : {
101 : public:
102 : /**
103 : * Constructor. Initializes Petsc data structures
104 : */
105 : PetscLinearSolver (const libMesh::Parallel::Communicator & comm_in);
106 :
107 : /**
108 : * Destructor.
109 : */
110 45644 : virtual ~PetscLinearSolver () = default;
111 :
112 : /**
113 : * Release all memory and clear data structures.
114 : */
115 : virtual void clear () override;
116 :
117 : /**
118 : * Initialize data structures if not done so already.
119 : * Assigns a name, which is turned into an underscore-separated
120 : * prefix for the underlying KSP object.
121 : */
122 : virtual void init (const char * name = nullptr) override;
123 :
124 : /**
125 : * Initialize data structures if not done so already plus much more
126 : */
127 : void init (PetscMatrixBase<T> * matrix,
128 : const char * name = nullptr);
129 :
130 : /**
131 : * Apply names to the system to be solved. This sets an option
132 : * prefix from the system name and sets field names from the
133 : * system's variable names.
134 : *
135 : * Since field names are applied to DoF numberings, this method must
136 : * be called again after any System reinit.
137 : */
138 : virtual void init_names (const System &) override;
139 :
140 : /**
141 : * After calling this method, all successive solves will be
142 : * restricted to the given set of dofs, which must contain local
143 : * dofs on each processor only and not contain any duplicates. This
144 : * mode can be disabled by calling this method with \p dofs being a
145 : * \p nullptr.
146 : */
147 : virtual void restrict_solve_to (const std::vector<unsigned int> * const dofs,
148 : const SubsetSolveMode subset_solve_mode=SUBSET_ZERO) override;
149 :
150 : /**
151 : * Call the Petsc solver. It calls the method below, using the
152 : * same matrix for the system and preconditioner matrices.
153 : */
154 : virtual std::pair<unsigned int, Real>
155 403041 : solve (SparseMatrix<T> & matrix_in,
156 : NumericVector<T> & solution_in,
157 : NumericVector<T> & rhs_in,
158 : const std::optional<double> tol = std::nullopt,
159 : const std::optional<unsigned int> m_its = std::nullopt) override
160 : {
161 403041 : return this->solve(matrix_in, matrix_in, solution_in, rhs_in, tol, m_its);
162 : }
163 :
164 :
165 : /**
166 : * Call the Petsc solver. It calls the method below, using the
167 : * same matrix for the system and preconditioner matrices.
168 : */
169 : virtual std::pair<unsigned int, Real>
170 : adjoint_solve (SparseMatrix<T> & matrix_in,
171 : NumericVector<T> & solution_in,
172 : NumericVector<T> & rhs_in,
173 : const std::optional<double> tol = std::nullopt,
174 : const std::optional<unsigned int> m_its = std::nullopt) override;
175 :
176 : /**
177 : * This method allows you to call a linear solver while specifying
178 : * the matrix to use as the (left) preconditioning matrix.
179 : *
180 : * \note The linear solver will not compute a preconditioner in this
181 : * case, and will instead premultiply by the matrix you provide. In
182 : * PETSc, this is accomplished by calling
183 : * \code
184 : * PCSetType(_pc, PCMAT);
185 : * \endcode
186 : * before invoking KSPSolve().
187 : *
188 : * \note This functionality is not implemented in the LinearSolver
189 : * class since there is not a built-in analog to this method for
190 : * LASPACK. You could probably implement it by hand if you wanted.
191 : */
192 : virtual std::pair<unsigned int, Real>
193 : solve (SparseMatrix<T> & matrix,
194 : SparseMatrix<T> & preconditioner,
195 : NumericVector<T> & solution,
196 : NumericVector<T> & rhs,
197 : const std::optional<double> tol = std::nullopt,
198 : const std::optional<unsigned int> m_its = std::nullopt) override;
199 :
200 : /**
201 : * This function solves a system whose matrix is a shell matrix.
202 : */
203 : virtual std::pair<unsigned int, Real>
204 : solve (const ShellMatrix<T> & shell_matrix,
205 : NumericVector<T> & solution_in,
206 : NumericVector<T> & rhs_in,
207 : const std::optional<double> tol = std::nullopt,
208 : const std::optional<unsigned int> m_its = std::nullopt) override;
209 :
210 : /**
211 : * This function solves a system whose matrix is a shell matrix, but
212 : * a sparse matrix is used as preconditioning matrix, this allowing
213 : * other preconditioners than JACOBI.
214 : */
215 : virtual std::pair<unsigned int, Real>
216 : solve (const ShellMatrix<T> & shell_matrix,
217 : const SparseMatrix<T> & precond_matrix,
218 : NumericVector<T> & solution_in,
219 : NumericVector<T> & rhs_in,
220 : const std::optional<double> tol = std::nullopt,
221 : const std::optional<unsigned int> m_its = std::nullopt) override;
222 :
223 : /**
224 : * \returns The raw PETSc preconditioner context pointer.
225 : *
226 : * This allows you to specify the PCShellSetApply() and
227 : * PCShellSetSetUp() functions if you desire. Just don't do
228 : * anything crazy like calling libMeshPCDestroy() on the pointer!
229 : */
230 219514 : PC pc() { this->init(); return _pc; }
231 :
232 : /**
233 : * \returns The raw PETSc ksp context pointer.
234 : *
235 : * This is useful if you are for example setting a custom
236 : * convergence test with KSPSetConvergenceTest().
237 : */
238 : KSP ksp();
239 :
240 : /**
241 : * Fills the input vector with the sequence of residual norms
242 : * from the latest iterative solve.
243 : */
244 : void get_residual_history(std::vector<double> & hist);
245 :
246 : /**
247 : * \returns Just the initial residual for the solve just
248 : * completed with this interface.
249 : *
250 : * Use this method instead of the one above if you just want the
251 : * starting residual and not the entire history.
252 : */
253 : Real get_initial_residual();
254 :
255 : /**
256 : * \returns The solver's convergence flag
257 : */
258 : virtual LinearConvergenceReason get_converged_reason() const override;
259 :
260 : private:
261 :
262 : typedef PetscErrorCode (*ksp_solve_func_type)(KSP,Vec,Vec);
263 :
264 : /*
265 : * Helper function to run solve() and adjoint_solve()
266 : */
267 : virtual std::pair<unsigned int, Real>
268 : solve_common (SparseMatrix<T> & matrix_in,
269 : SparseMatrix<T> & precond_in,
270 : NumericVector<T> & solution_in,
271 : NumericVector<T> & rhs_in,
272 : const double rel_tol,
273 : const double abs_tol,
274 : const unsigned int m_its,
275 : ksp_solve_func_type solve_func);
276 :
277 : /*
278 : * Helper function to run ShellMatrix solve() with or without a
279 : * separate preconditioner matrix
280 : */
281 : virtual std::pair<unsigned int, Real>
282 : shell_solve_common (const ShellMatrix<T> & shell_matrix,
283 : PetscMatrixBase<T> * precond_matrix,
284 : NumericVector<T> & solution_in,
285 : NumericVector<T> & rhs_in,
286 : const double rel_tol,
287 : const double abs_tol,
288 : const unsigned int m_its);
289 :
290 : /*
291 : * Helper function for the helper functions
292 : */
293 : std::pair<unsigned int, Real>
294 : solve_base (SparseMatrix<T> * matrix,
295 : PetscMatrixBase<T> * precond,
296 : Mat mat,
297 : NumericVector<T> & solution_in,
298 : NumericVector<T> & rhs_in,
299 : const double rel_tol,
300 : const double abs_tol,
301 : const unsigned int m_its,
302 : ksp_solve_func_type solve_func);
303 :
304 :
305 : /**
306 : * Tells PETSC to use the user-specified solver stored in
307 : * \p _solver_type
308 : */
309 : void set_petsc_solver_type ();
310 :
311 : /**
312 : * Internal function if shell matrix mode is used.
313 : */
314 : static PetscErrorCode _petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest);
315 :
316 : /**
317 : * Internal function if shell matrix mode is used.
318 : */
319 : static PetscErrorCode _petsc_shell_matrix_mult_add(Mat mat, Vec arg, Vec add, Vec dest);
320 :
321 : /**
322 : * Internal function if shell matrix mode is used.
323 : */
324 : static PetscErrorCode _petsc_shell_matrix_get_diagonal(Mat mat, Vec dest);
325 :
326 : /**
327 : * Preconditioner context
328 : */
329 : PC _pc;
330 :
331 : /**
332 : * Krylov subspace context
333 : */
334 : WrappedPetsc<KSP> _ksp;
335 :
336 : /**
337 : * PETSc index set containing the dofs on which to solve (\p nullptr
338 : * means solve on all dofs).
339 : */
340 : WrappedPetsc<IS> _restrict_solve_to_is;
341 :
342 : /**
343 : * PETSc index set, complement to \p _restrict_solve_to_is. This
344 : * will be created on demand by the method \p
345 : * _create_complement_is().
346 : */
347 : WrappedPetsc<IS> _restrict_solve_to_is_complement;
348 :
349 : /**
350 : * \returns The local size of \p _restrict_solve_to_is.
351 : */
352 : PetscInt restrict_solve_to_is_local_size() const;
353 :
354 : /**
355 : * Creates \p _restrict_solve_to_is_complement to contain all
356 : * indices that are local in \p vec_in, except those that are
357 : * contained in \p _restrict_solve_to_is.
358 : */
359 : void create_complement_is (const NumericVector<T> & vec_in);
360 :
361 : /**
362 : * If restrict-solve-to-subset mode is active, this member decides
363 : * what happens with the dofs outside the subset.
364 : */
365 : SubsetSolveMode _subset_solve_mode;
366 : };
367 :
368 : } // namespace libMesh
369 :
370 :
371 : #endif // #ifdef LIBMESH_HAVE_PETSC
372 : #endif // LIBMESH_PETSC_LINEAR_SOLVER_H
|
