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