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 1082171 : bool initialized () const { return _is_initialized; }
86 :
87 : /**
88 : * Release all memory and clear data structures.
89 : */
90 580 : 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 : * \returns The currently attached solver configuration object, if any.
259 : */
260 : SolverConfiguration * solver_configuration() const;
261 :
262 : /**
263 : * Set the solver configuration object. The solver does not take
264 : * ownership. A null pointer detaches the current configuration.
265 : */
266 : void set_solver_configuration(SolverConfiguration * solver_configuration);
267 :
268 : /**
269 : * Set the solver configuration object.
270 : */
271 0 : void set_solver_configuration(SolverConfiguration & solver_configuration)
272 : {
273 0 : this->set_solver_configuration(&solver_configuration);
274 0 : }
275 :
276 : protected:
277 : /**
278 : * Enum stating which type of iterative solver to use.
279 : */
280 : SolverType _solver_type;
281 :
282 : /**
283 : * Enum stating with type of preconditioner to use.
284 : */
285 : PreconditionerType _preconditioner_type;
286 :
287 : /**
288 : * Flag indicating if the data structures have been initialized.
289 : */
290 : bool _is_initialized;
291 :
292 : /**
293 : * Holds the Preconditioner object to be used for the linear solves.
294 : */
295 : Preconditioner<T> * _preconditioner;
296 :
297 : /**
298 : * Boolean flag to indicate whether we want to use an identical
299 : * preconditioner to the previous solve. This can save
300 : * substantial work in the cases where the system matrix is
301 : * the same for successive solves.
302 : */
303 : bool same_preconditioner;
304 :
305 : /**
306 : * Optionally store a SolverOptions object that can be used
307 : * to set parameters like solver type, tolerances and iteration limits.
308 : */
309 : SolverConfiguration * _solver_configuration;
310 :
311 : /**
312 : * Get solver settings based on optional parameters and the solver
313 : * configuration object
314 : */
315 : double get_real_solver_setting(const std::string & setting_name,
316 : const std::optional<double> & setting,
317 : const std::optional<double> default_value = std::nullopt);
318 :
319 : /**
320 : * Get solver settings based on optional parameters and the solver
321 : * configuration object
322 : */
323 : int get_int_solver_setting(const std::string & setting_name,
324 : const std::optional<int> & setting,
325 : const std::optional<int> default_value = std::nullopt);
326 : };
327 :
328 :
329 :
330 :
331 : /*----------------------- inline functions ----------------------------------*/
332 : template <typename T>
333 : inline
334 580 : LinearSolver<T>::~LinearSolver ()
335 : {
336 580 : this->LinearSolver::clear ();
337 580 : }
338 :
339 : template <typename T>
340 : inline
341 0 : bool LinearSolver<T>::get_same_preconditioner()
342 : {
343 0 : return same_preconditioner;
344 : }
345 :
346 : template <typename T>
347 : inline
348 : std::pair<unsigned int, Real>
349 257589 : LinearSolver<T>::solve (SparseMatrix<T> & mat,
350 : SparseMatrix<T> * pc_mat,
351 : NumericVector<T> & sol,
352 : NumericVector<T> & rhs,
353 : const std::optional<double> tol,
354 : const std::optional<unsigned int> n_iter)
355 : {
356 257589 : if (pc_mat)
357 0 : return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
358 : else
359 257589 : return this->solve(mat, sol, rhs, tol, n_iter);
360 : }
361 :
362 :
363 : template <typename T>
364 : inline
365 : std::pair<unsigned int, Real>
366 70 : LinearSolver<T>::solve (const ShellMatrix<T> & mat,
367 : const SparseMatrix<T> * pc_mat,
368 : NumericVector<T> & sol,
369 : NumericVector<T> & rhs,
370 : const std::optional<double> tol,
371 : const std::optional<unsigned int> n_iter)
372 : {
373 70 : if (pc_mat)
374 70 : return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
375 : else
376 0 : return this->solve(mat, sol, rhs, tol, n_iter);
377 : }
378 :
379 :
380 : } // namespace libMesh
381 :
382 :
383 : #endif // LIBMESH_LINEAR_SOLVER_H
|
