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_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 1083299 : bool initialized () const { return _is_initialized; }
86 :
87 : /**
88 : * Release all memory and clear data structures.
89 : */
90 842 : 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 and set auxiliary preconditioner
100 : * data. For most packages this is a no-op; for PETSc this a) sets an option
101 : * prefix from the system name and sets field names from the system's
102 : * variable names and b) sets auxiliary data needed by preconditioners such
103 : * as hypre ams/ads.
104 : *
105 : * Since field names are applied to DoF numberings, this method must
106 : * be called again after any System reinit.
107 : */
108 3173 : virtual void init_systems (const System &) {}
109 :
110 : /**
111 : * \returns The type of solver to use.
112 : */
113 0 : SolverType solver_type () const { return _solver_type; }
114 :
115 : /**
116 : * Sets the type of solver to use.
117 : */
118 0 : void set_solver_type (const SolverType st)
119 0 : { _solver_type = st; }
120 :
121 : /**
122 : * \returns The type of preconditioner to use.
123 : */
124 : PreconditionerType preconditioner_type () const;
125 :
126 : /**
127 : * Sets the type of preconditioner to use.
128 : */
129 : void set_preconditioner_type (const PreconditionerType pct);
130 :
131 : /**
132 : * Attaches a Preconditioner object to be used
133 : */
134 : void attach_preconditioner(Preconditioner<T> * preconditioner);
135 :
136 : /**
137 : * Set the same_preconditioner flag, which indicates if we reuse the
138 : * same preconditioner for subsequent solves.
139 : */
140 : virtual void reuse_preconditioner(bool );
141 :
142 : /**
143 : * \returns \p same_preconditioner, which indicates if we reuse the
144 : * same preconditioner for subsequent solves.
145 : */
146 : bool get_same_preconditioner();
147 :
148 : /**
149 : * After calling this method, all successive solves will be
150 : * restricted to the given set of dofs, which must contain local
151 : * dofs on each processor only and not contain any duplicates. This
152 : * mode can be disabled by calling this method with \p dofs being a
153 : * \p nullptr.
154 : */
155 : virtual void restrict_solve_to (const std::vector<unsigned int> * const dofs,
156 : const SubsetSolveMode subset_solve_mode=SUBSET_ZERO);
157 :
158 : /**
159 : * This function calls the solver \p _solver_type preconditioned
160 : * with the \p _preconditioner_type preconditioner.
161 : *
162 : * \note This method will compute the preconditioner from the system
163 : * matrix.
164 : */
165 : virtual std::pair<unsigned int, Real> solve (SparseMatrix<T> &, // System Matrix
166 : NumericVector<T> &, // Solution vector
167 : NumericVector<T> &, // RHS vector
168 : const std::optional<double> tol = std::nullopt,
169 : const std::optional<unsigned int> m_its = std::nullopt) = 0; // N. Iterations
170 :
171 : /**
172 : * Function to solve the adjoint system.
173 : *
174 : * \note This method will compute the preconditioner from the system
175 : * matrix. This is not a pure virtual function and is defined
176 : * linear_solver.C
177 : */
178 : virtual std::pair<unsigned int, Real> adjoint_solve (SparseMatrix<T> &, // System Matrix
179 : NumericVector<T> &, // Solution vector
180 : NumericVector<T> &, // RHS vector
181 : const std::optional<double> tol = std::nullopt,
182 : const std::optional<unsigned int> m_its = std::nullopt); // N. Iterations
183 :
184 : /**
185 : * This function calls the solver
186 : * "_solver_type" preconditioned with the
187 : * "_preconditioner_type" preconditioner.
188 : */
189 : virtual std::pair<unsigned int, Real> solve (SparseMatrix<T> &, // System Matrix
190 : SparseMatrix<T> &, // Preconditioning Matrix
191 : NumericVector<T> &, // Solution vector
192 : NumericVector<T> &, // RHS vector
193 : const std::optional<double> tol = std::nullopt,
194 : const std::optional<unsigned int> m_its = std::nullopt) = 0; // N. Iterations
195 :
196 : /**
197 : * This function calls the solver "_solver_type" preconditioned with
198 : * the "_preconditioner_type" preconditioner. The preconditioning
199 : * matrix is used if it is provided, or the system matrix is used if
200 : * \p precond_matrix is null
201 : */
202 : std::pair<unsigned int, Real> solve (SparseMatrix<T> & matrix,
203 : SparseMatrix<T> * precond_matrix,
204 : NumericVector<T> &, // Solution vector
205 : NumericVector<T> &, // RHS vector
206 : const std::optional<double> tol = std::nullopt,
207 : const std::optional<unsigned int> m_its = std::nullopt); // N. Iterations
208 :
209 :
210 :
211 : /**
212 : * This function solves a system whose matrix is a shell matrix.
213 : */
214 : virtual std::pair<unsigned int, Real> solve (const ShellMatrix<T> & shell_matrix,
215 : NumericVector<T> &, // Solution vector
216 : NumericVector<T> &, // RHS vector
217 : const std::optional<double> tol = std::nullopt,
218 : const std::optional<unsigned int> m_its = std::nullopt) = 0; // N. Iterations
219 :
220 :
221 : /**
222 : * This function solves a system whose matrix is a shell matrix, but
223 : * a sparse matrix is used as preconditioning matrix, this allowing
224 : * other preconditioners than JACOBI.
225 : */
226 : virtual std::pair<unsigned int, Real> solve (const ShellMatrix<T> & shell_matrix,
227 : const SparseMatrix<T> & precond_matrix,
228 : NumericVector<T> &, // Solution vector
229 : NumericVector<T> &, // RHS vector
230 : const std::optional<double> tol = std::nullopt,
231 : const std::optional<unsigned int> m_its = std::nullopt) = 0; // N. Iterations
232 :
233 :
234 : /**
235 : * This function solves a system whose matrix is a shell matrix, but
236 : * an optional sparse matrix may be used as preconditioning matrix.
237 : */
238 : std::pair<unsigned int, Real> solve (const ShellMatrix<T> & matrix,
239 : const SparseMatrix<T> * precond_matrix,
240 : NumericVector<T> &, // Solution vector
241 : NumericVector<T> &, // RHS vector
242 : const std::optional<double> tol = std::nullopt,
243 : const std::optional<unsigned int> m_its = std::nullopt); // N. Iterations
244 :
245 :
246 : /**
247 : * Prints a useful message about why the latest linear solve
248 : * con(di)verged.
249 : */
250 : virtual void print_converged_reason() const;
251 :
252 : /**
253 : * \returns The solver's convergence flag
254 : */
255 : virtual LinearConvergenceReason get_converged_reason() const = 0;
256 :
257 : /**
258 : * Set the solver configuration object.
259 : */
260 : void set_solver_configuration(SolverConfiguration & solver_configuration);
261 :
262 : protected:
263 : /**
264 : * Enum stating which type of iterative solver to use.
265 : */
266 : SolverType _solver_type;
267 :
268 : /**
269 : * Enum stating with type of preconditioner to use.
270 : */
271 : PreconditionerType _preconditioner_type;
272 :
273 : /**
274 : * Flag indicating if the data structures have been initialized.
275 : */
276 : bool _is_initialized;
277 :
278 : /**
279 : * Holds the Preconditioner object to be used for the linear solves.
280 : */
281 : Preconditioner<T> * _preconditioner;
282 :
283 : /**
284 : * Boolean flag to indicate whether we want to use an identical
285 : * preconditioner to the previous solve. This can save
286 : * substantial work in the cases where the system matrix is
287 : * the same for successive solves.
288 : */
289 : bool same_preconditioner;
290 :
291 : /**
292 : * Optionally store a SolverOptions object that can be used
293 : * to set parameters like solver type, tolerances and iteration limits.
294 : */
295 : SolverConfiguration * _solver_configuration;
296 :
297 : /**
298 : * Get solver settings based on optional parameters and the solver
299 : * configuration object
300 : */
301 : double get_real_solver_setting(const std::string & setting_name,
302 : const std::optional<double> & setting,
303 : const std::optional<double> default_value = std::nullopt);
304 :
305 : /**
306 : * Get solver settings based on optional parameters and the solver
307 : * configuration object
308 : */
309 : int get_int_solver_setting(const std::string & setting_name,
310 : const std::optional<int> & setting,
311 : const std::optional<int> default_value = std::nullopt);
312 : };
313 :
314 :
315 :
316 :
317 : /*----------------------- inline functions ----------------------------------*/
318 : template <typename T>
319 : inline
320 842 : LinearSolver<T>::~LinearSolver ()
321 : {
322 842 : this->LinearSolver::clear ();
323 842 : }
324 :
325 : template <typename T>
326 : inline
327 0 : bool LinearSolver<T>::get_same_preconditioner()
328 : {
329 0 : return same_preconditioner;
330 : }
331 :
332 : template <typename T>
333 : inline
334 : std::pair<unsigned int, Real>
335 257234 : LinearSolver<T>::solve (SparseMatrix<T> & mat,
336 : SparseMatrix<T> * pc_mat,
337 : NumericVector<T> & sol,
338 : NumericVector<T> & rhs,
339 : const std::optional<double> tol,
340 : const std::optional<unsigned int> n_iter)
341 : {
342 257234 : if (pc_mat)
343 0 : return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
344 : else
345 257234 : return this->solve(mat, sol, rhs, tol, n_iter);
346 : }
347 :
348 :
349 : template <typename T>
350 : inline
351 : std::pair<unsigned int, Real>
352 70 : LinearSolver<T>::solve (const ShellMatrix<T> & mat,
353 : const SparseMatrix<T> * pc_mat,
354 : NumericVector<T> & sol,
355 : NumericVector<T> & rhs,
356 : const std::optional<double> tol,
357 : const std::optional<unsigned int> n_iter)
358 : {
359 70 : if (pc_mat)
360 70 : return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
361 : else
362 0 : return this->solve(mat, sol, rhs, tol, n_iter);
363 : }
364 :
365 :
366 : } // namespace libMesh
367 :
368 :
369 : #endif // LIBMESH_LINEAR_SOLVER_H
|
