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_LINEAR_SOLVER_H
21 : #define LIBMESH_LINEAR_SOLVER_H
22 :
23 : // Local includes
24 : #include "libmesh/libmesh_common.h"
25 : #include "libmesh/enum_subset_solve_mode.h" // SUBSET_ZERO
26 : #include "libmesh/reference_counted_object.h"
27 : #include "libmesh/libmesh.h"
28 : #include "libmesh/parallel_object.h"
29 :
30 : // C++ includes
31 : #include <cstddef>
32 : #include <vector>
33 : #include <memory>
34 : #include <optional>
35 :
36 : namespace libMesh
37 : {
38 :
39 : // forward declarations
40 : template <typename T> class SparseMatrix;
41 : template <typename T> class NumericVector;
42 : template <typename T> class ShellMatrix;
43 : template <typename T> class Preconditioner;
44 : class System;
45 : class SolverConfiguration;
46 : enum SolverPackage : int;
47 : enum PreconditionerType : int;
48 : enum SolverType : int;
49 : enum LinearConvergenceReason : int;
50 :
51 : /**
52 : * This base class can be inherited from to provide interfaces to
53 : * linear solvers from different packages like PETSc and LASPACK.
54 : *
55 : * \author Benjamin Kirk
56 : * \date 2003
57 : */
58 : template <typename T>
59 : class LinearSolver : public ReferenceCountedObject<LinearSolver<T>>,
60 : public ParallelObject
61 : {
62 : public:
63 :
64 : /**
65 : * Constructor. Initializes Solver data structures
66 : */
67 : LinearSolver (const libMesh::Parallel::Communicator & comm_in);
68 :
69 : /**
70 : * Destructor.
71 : */
72 : virtual ~LinearSolver ();
73 :
74 : /**
75 : * Builds a \p LinearSolver using the linear solver package specified by
76 : * \p solver_package
77 : */
78 : static std::unique_ptr<LinearSolver<T>> build(const libMesh::Parallel::Communicator & comm_in,
79 : const SolverPackage solver_package = libMesh::default_solver_package());
80 :
81 : /**
82 : * \returns \p true if the data structures are
83 : * initialized, false otherwise.
84 : */
85 1036414 : bool initialized () const { return _is_initialized; }
86 :
87 : /**
88 : * Release all memory and clear data structures.
89 : */
90 674 : virtual void clear () {}
91 :
92 : /**
93 : * Initialize data structures if not done so already.
94 : * May assign a name to the solver in some implementations
95 : */
96 : virtual void init (const char * name = nullptr) = 0;
97 :
98 : /**
99 : * Apply names to the system to be solved. For most packages this
100 : * is a no-op; for PETSc this sets an option prefix from the system
101 : * name and sets field names from the system's variable names.
102 : *
103 : * Since field names are applied to DoF numberings, this method must
104 : * be called again after any System reinit.
105 : */
106 3143 : virtual void init_names (const System &) {}
107 :
108 : /**
109 : * \returns The type of solver to use.
110 : */
111 0 : SolverType solver_type () const { return _solver_type; }
112 :
113 : /**
114 : * Sets the type of solver to use.
115 : */
116 0 : void set_solver_type (const SolverType st)
117 0 : { _solver_type = st; }
118 :
119 : /**
120 : * \returns The type of preconditioner to use.
121 : */
122 : PreconditionerType preconditioner_type () const;
123 :
124 : /**
125 : * Sets the type of preconditioner to use.
126 : */
127 : void set_preconditioner_type (const PreconditionerType pct);
128 :
129 : /**
130 : * Attaches a Preconditioner object to be used
131 : */
132 : void attach_preconditioner(Preconditioner<T> * preconditioner);
133 :
134 : /**
135 : * Set the same_preconditioner flag, which indicates if we reuse the
136 : * same preconditioner for subsequent solves.
137 : */
138 : virtual void reuse_preconditioner(bool );
139 :
140 : /**
141 : * \returns \p same_preconditioner, which indicates if we reuse the
142 : * same preconditioner for subsequent solves.
143 : */
144 : bool get_same_preconditioner();
145 :
146 : /**
147 : * After calling this method, all successive solves will be
148 : * restricted to the given set of dofs, which must contain local
149 : * dofs on each processor only and not contain any duplicates. This
150 : * mode can be disabled by calling this method with \p dofs being a
151 : * \p nullptr.
152 : */
153 : virtual void restrict_solve_to (const std::vector<unsigned int> * const dofs,
154 : const SubsetSolveMode subset_solve_mode=SUBSET_ZERO);
155 :
156 : /**
157 : * This function calls the solver \p _solver_type preconditioned
158 : * with the \p _preconditioner_type preconditioner.
159 : *
160 : * \note This method will compute the preconditioner from the system
161 : * matrix.
162 : */
163 : virtual std::pair<unsigned int, Real> solve (SparseMatrix<T> &, // System Matrix
164 : NumericVector<T> &, // Solution vector
165 : NumericVector<T> &, // RHS vector
166 : const std::optional<double> tol = std::nullopt,
167 : const std::optional<unsigned int> m_its = std::nullopt) = 0; // N. Iterations
168 :
169 : /**
170 : * Function to solve the adjoint system.
171 : *
172 : * \note This method will compute the preconditioner from the system
173 : * matrix. This is not a pure virtual function and is defined
174 : * linear_solver.C
175 : */
176 : virtual std::pair<unsigned int, Real> adjoint_solve (SparseMatrix<T> &, // System Matrix
177 : NumericVector<T> &, // Solution vector
178 : NumericVector<T> &, // RHS vector
179 : const std::optional<double> tol = std::nullopt,
180 : const std::optional<unsigned int> m_its = std::nullopt); // N. Iterations
181 :
182 : /**
183 : * This function calls the solver
184 : * "_solver_type" preconditioned with the
185 : * "_preconditioner_type" preconditioner.
186 : */
187 : virtual std::pair<unsigned int, Real> solve (SparseMatrix<T> &, // System Matrix
188 : SparseMatrix<T> &, // Preconditioning Matrix
189 : NumericVector<T> &, // Solution vector
190 : NumericVector<T> &, // RHS vector
191 : const std::optional<double> tol = std::nullopt,
192 : const std::optional<unsigned int> m_its = std::nullopt) = 0; // N. Iterations
193 :
194 : /**
195 : * This function calls the solver "_solver_type" preconditioned with
196 : * the "_preconditioner_type" preconditioner. The preconditioning
197 : * matrix is used if it is provided, or the system matrix is used if
198 : * \p precond_matrix is null
199 : */
200 : std::pair<unsigned int, Real> solve (SparseMatrix<T> & matrix,
201 : SparseMatrix<T> * precond_matrix,
202 : NumericVector<T> &, // Solution vector
203 : NumericVector<T> &, // RHS vector
204 : const std::optional<double> tol = std::nullopt,
205 : const std::optional<unsigned int> m_its = std::nullopt); // N. Iterations
206 :
207 :
208 :
209 : /**
210 : * This function solves a system whose matrix is a shell matrix.
211 : */
212 : virtual std::pair<unsigned int, Real> solve (const ShellMatrix<T> & shell_matrix,
213 : NumericVector<T> &, // Solution vector
214 : NumericVector<T> &, // RHS vector
215 : const std::optional<double> tol = std::nullopt,
216 : const std::optional<unsigned int> m_its = std::nullopt) = 0; // N. Iterations
217 :
218 :
219 : /**
220 : * This function solves a system whose matrix is a shell matrix, but
221 : * a sparse matrix is used as preconditioning matrix, this allowing
222 : * other preconditioners than JACOBI.
223 : */
224 : virtual std::pair<unsigned int, Real> solve (const ShellMatrix<T> & shell_matrix,
225 : const SparseMatrix<T> & precond_matrix,
226 : NumericVector<T> &, // Solution vector
227 : NumericVector<T> &, // RHS vector
228 : const std::optional<double> tol = std::nullopt,
229 : const std::optional<unsigned int> m_its = std::nullopt) = 0; // N. Iterations
230 :
231 :
232 : /**
233 : * This function solves a system whose matrix is a shell matrix, but
234 : * an optional sparse matrix may be used as preconditioning matrix.
235 : */
236 : std::pair<unsigned int, Real> solve (const ShellMatrix<T> & matrix,
237 : const SparseMatrix<T> * precond_matrix,
238 : NumericVector<T> &, // Solution vector
239 : NumericVector<T> &, // RHS vector
240 : const std::optional<double> tol = std::nullopt,
241 : const std::optional<unsigned int> m_its = std::nullopt); // N. Iterations
242 :
243 :
244 : /**
245 : * Prints a useful message about why the latest linear solve
246 : * con(di)verged.
247 : */
248 : virtual void print_converged_reason() const;
249 :
250 : /**
251 : * \returns The solver's convergence flag
252 : */
253 : virtual LinearConvergenceReason get_converged_reason() const = 0;
254 :
255 : /**
256 : * Set the solver configuration object.
257 : */
258 : void set_solver_configuration(SolverConfiguration & solver_configuration);
259 :
260 : protected:
261 : /**
262 : * Enum stating which type of iterative solver to use.
263 : */
264 : SolverType _solver_type;
265 :
266 : /**
267 : * Enum stating with type of preconditioner to use.
268 : */
269 : PreconditionerType _preconditioner_type;
270 :
271 : /**
272 : * Flag indicating if the data structures have been initialized.
273 : */
274 : bool _is_initialized;
275 :
276 : /**
277 : * Holds the Preconditioner object to be used for the linear solves.
278 : */
279 : Preconditioner<T> * _preconditioner;
280 :
281 : /**
282 : * Boolean flag to indicate whether we want to use an identical
283 : * preconditioner to the previous solve. This can save
284 : * substantial work in the cases where the system matrix is
285 : * the same for successive solves.
286 : */
287 : bool same_preconditioner;
288 :
289 : /**
290 : * Optionally store a SolverOptions object that can be used
291 : * to set parameters like solver type, tolerances and iteration limits.
292 : */
293 : SolverConfiguration * _solver_configuration;
294 :
295 : /**
296 : * Get solver settings based on optional parameters and the solver
297 : * configuration object
298 : */
299 : double get_real_solver_setting(const std::string & setting_name,
300 : const std::optional<double> & setting,
301 : const std::optional<double> default_value = std::nullopt);
302 :
303 : /**
304 : * Get solver settings based on optional parameters and the solver
305 : * configuration object
306 : */
307 : int get_int_solver_setting(const std::string & setting_name,
308 : const std::optional<int> & setting,
309 : const std::optional<int> default_value = std::nullopt);
310 : };
311 :
312 :
313 :
314 :
315 : /*----------------------- inline functions ----------------------------------*/
316 : template <typename T>
317 : inline
318 674 : LinearSolver<T>::~LinearSolver ()
319 : {
320 674 : this->LinearSolver::clear ();
321 674 : }
322 :
323 : template <typename T>
324 : inline
325 0 : bool LinearSolver<T>::get_same_preconditioner()
326 : {
327 0 : return same_preconditioner;
328 : }
329 :
330 : template <typename T>
331 : inline
332 : std::pair<unsigned int, Real>
333 244177 : LinearSolver<T>::solve (SparseMatrix<T> & mat,
334 : SparseMatrix<T> * pc_mat,
335 : NumericVector<T> & sol,
336 : NumericVector<T> & rhs,
337 : const std::optional<double> tol,
338 : const std::optional<unsigned int> n_iter)
339 : {
340 244177 : if (pc_mat)
341 0 : return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
342 : else
343 244177 : return this->solve(mat, sol, rhs, tol, n_iter);
344 : }
345 :
346 :
347 : template <typename T>
348 : inline
349 : std::pair<unsigned int, Real>
350 70 : LinearSolver<T>::solve (const ShellMatrix<T> & mat,
351 : const SparseMatrix<T> * pc_mat,
352 : NumericVector<T> & sol,
353 : NumericVector<T> & rhs,
354 : const std::optional<double> tol,
355 : const std::optional<unsigned int> n_iter)
356 : {
357 70 : if (pc_mat)
358 70 : return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
359 : else
360 0 : return this->solve(mat, sol, rhs, tol, n_iter);
361 : }
362 :
363 :
364 : } // namespace libMesh
365 :
366 :
367 : #endif // LIBMESH_LINEAR_SOLVER_H
|
