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