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_NONLINEAR_SOLVER_H
21 : #define LIBMESH_NONLINEAR_SOLVER_H
22 :
23 : // Local includes
24 : #include "libmesh/libmesh_common.h"
25 : #include "libmesh/reference_counted_object.h"
26 : #include "libmesh/nonlinear_implicit_system.h"
27 : #include "libmesh/libmesh.h"
28 : #include "libmesh/parallel_object.h"
29 :
30 : // C++ includes
31 : #include <cstddef>
32 : #include <memory>
33 :
34 : namespace libMesh
35 : {
36 :
37 : // forward declarations
38 : template <typename T> class SparseMatrix;
39 : template <typename T> class NumericVector;
40 : template <typename T> class Preconditioner;
41 : class SolverConfiguration;
42 : enum SolverPackage : int;
43 :
44 : /**
45 : * This base class can be inherited from to provide interfaces to
46 : * nonlinear solvers from different packages like PETSc and Trilinos.
47 : *
48 : * \author Benjamin Kirk
49 : * \date 2005
50 : */
51 : template <typename T>
52 : class NonlinearSolver : public ReferenceCountedObject<NonlinearSolver<T>>,
53 : public ParallelObject
54 : {
55 : public:
56 : /**
57 : * The type of system
58 : */
59 : typedef NonlinearImplicitSystem sys_type;
60 :
61 : /**
62 : * Constructor. Initializes Solver data structures
63 : */
64 : explicit
65 : NonlinearSolver (sys_type & s);
66 :
67 : /**
68 : * Destructor.
69 : */
70 : virtual ~NonlinearSolver ();
71 :
72 : /**
73 : * Builds a \p NonlinearSolver using the nonlinear solver package specified by
74 : * \p solver_package
75 : */
76 : static std::unique_ptr<NonlinearSolver<T>> build(sys_type & s,
77 : const SolverPackage solver_package = libMesh::default_solver_package());
78 :
79 : /**
80 : * \returns \p true if the data structures are
81 : * initialized, false otherwise.
82 : */
83 206214 : bool initialized () const { return _is_initialized; }
84 :
85 : /**
86 : * Release all memory and clear data structures.
87 : */
88 42 : virtual void clear () {}
89 :
90 : /**
91 : * Initialize data structures if not done so already.
92 : * May assign a name to the solver in some implementations
93 : */
94 : virtual void init (const char * name = nullptr) = 0;
95 :
96 : /**
97 : * Solves the nonlinear system.
98 : */
99 : virtual std::pair<unsigned int, Real> solve (SparseMatrix<T> &, // System Jacobian Matrix
100 : NumericVector<T> &, // Solution vector
101 : NumericVector<T> &, // Residual vector
102 : const double, // Stopping tolerance
103 : const unsigned int) = 0; // N. Iterations
104 :
105 : /**
106 : * Solves the nonlinear system using \p jac_in as the actual Jacobian operator and \p pre_in as
107 : * the preconditioning matrix -- which may be the same object (the common case) or genuinely
108 : * distinct (e.g. a matrix-free operator paired with an assembled preconditioning matrix). The
109 : * default implementation ignores \p jac_in and simply forwards to the single-matrix solve(),
110 : * for solver backends that do not distinguish between the two.
111 : */
112 0 : virtual std::pair<unsigned int, Real> solve (SparseMatrix<T> & /* jac_in */,
113 : SparseMatrix<T> & pre_in,
114 : NumericVector<T> & x_in,
115 : NumericVector<T> & r_in,
116 : const double tol,
117 : const unsigned int m_its)
118 : {
119 0 : return this->solve(pre_in, x_in, r_in, tol, m_its);
120 : }
121 :
122 : /**
123 : * Prints a useful message about why the latest nonlinear solve
124 : * con(di)verged.
125 : */
126 0 : virtual void print_converged_reason() { libmesh_not_implemented(); }
127 :
128 : /**
129 : * Get the total number of linear iterations done in the last solve
130 : */
131 : virtual int get_total_linear_iterations() = 0;
132 :
133 : /**
134 : * \returns The current nonlinear iteration number if called
135 : * *during* the solve(), for example by the user-specified residual
136 : * or Jacobian function.
137 : *
138 : * Must be overridden in derived classes.
139 : */
140 : virtual unsigned get_current_nonlinear_iteration_number() const = 0;
141 :
142 : /**
143 : * Function that computes the residual \p R(X) of the nonlinear system
144 : * at the input iterate \p X.
145 : */
146 : void (* residual) (const NumericVector<Number> & X,
147 : NumericVector<Number> & R,
148 : sys_type & S);
149 :
150 : /**
151 : * Object that computes the residual \p R(X) of the nonlinear system
152 : * at the input iterate \p X.
153 : */
154 : NonlinearImplicitSystem::ComputeResidual * residual_object;
155 :
156 : /**
157 : * Object that computes the residual \p R(X) of the nonlinear system
158 : * at the input iterate \p X for the purpose of forming a finite-differenced Jacobian.
159 : */
160 : NonlinearImplicitSystem::ComputeResidual * fd_residual_object;
161 :
162 : /**
163 : * Object that computes the residual \p R(X) of the nonlinear system
164 : * at the input iterate \p X for the purpose of forming Jacobian-vector products
165 : * via finite differencing.
166 : */
167 : NonlinearImplicitSystem::ComputeResidual * mffd_residual_object;
168 :
169 : /**
170 : * Function that computes the Jacobian \p J(X) of the nonlinear system
171 : * at the input iterate \p X.
172 : */
173 : void (* jacobian) (const NumericVector<Number> & X,
174 : SparseMatrix<Number> & J,
175 : sys_type & S);
176 :
177 : /**
178 : * Object that computes the Jacobian \p J(X) of the nonlinear system
179 : * at the input iterate \p X.
180 : */
181 : NonlinearImplicitSystem::ComputeJacobian * jacobian_object;
182 :
183 : /**
184 : * Function that computes either the residual \f$ R(X) \f$ or the
185 : * Jacobian \f$ J(X) \f$ of the nonlinear system at the input
186 : * iterate \f$ X \f$.
187 : *
188 : * \note Either \p R or \p J could be \p nullptr.
189 : */
190 : void (* matvec) (const NumericVector<Number> & X,
191 : NumericVector<Number> * R,
192 : SparseMatrix<Number> * J,
193 : sys_type & S);
194 :
195 : /**
196 : * Object that computes either the residual \f$ R(X) \f$ or the
197 : * Jacobian \f$ J(X) \f$ of the nonlinear system at the input
198 : * iterate \f$ X \f$.
199 : *
200 : * \note Either \p R or \p J could be \p nullptr.
201 : */
202 : NonlinearImplicitSystem::ComputeResidualandJacobian * residual_and_jacobian_object;
203 :
204 : /**
205 : * Function that computes the lower and upper bounds \p XL and \p XU on the solution of the nonlinear system.
206 : */
207 : void (* bounds) (NumericVector<Number> & XL,
208 : NumericVector<Number> & XU,
209 : sys_type & S);
210 : /**
211 : * Object that computes the bounds vectors \f$ XL \f$ and \f$ XU \f$.
212 : */
213 : NonlinearImplicitSystem::ComputeBounds * bounds_object;
214 :
215 : /**
216 : * Function that computes a basis for the Jacobian's nullspace --
217 : * the kernel or the "zero energy modes" -- that can be used in
218 : * solving a degenerate problem iteratively, if the solver supports it
219 : * (e.g., PETSc's KSP).
220 : */
221 : void (* nullspace) (std::vector<NumericVector<Number> *> & sp, sys_type & S);
222 :
223 : /**
224 : * A callable object that computes a basis for the Jacobian's nullspace --
225 : * the kernel or the "zero energy modes" -- that can be used in
226 : * solving a degenerate problem iteratively, if the solver supports it
227 : * (e.g., PETSc's KSP).
228 : */
229 : NonlinearImplicitSystem::ComputeVectorSubspace * nullspace_object;
230 :
231 : /**
232 : * Function that computes a basis for the transpose Jacobian's nullspace --
233 : * when solving a degenerate problem iteratively, if the solver supports it
234 : * (e.g., PETSc's KSP), it is used to remove contributions outside of R(jac)
235 : */
236 : void (* transpose_nullspace) (std::vector<NumericVector<Number> *> & sp, sys_type & S);
237 :
238 : /**
239 : * A callable object that computes a basis for the transpose Jacobian's nullspace --
240 : * when solving a degenerate problem iteratively, if the solver supports it
241 : * (e.g., PETSc's KSP), it is used to remove contributions outside of R(jac)
242 : */
243 : NonlinearImplicitSystem::ComputeVectorSubspace * transpose_nullspace_object;
244 :
245 : /**
246 : * Function that computes a basis for the Jacobian's near nullspace --
247 : * the set of "low energy modes" -- that can be used for AMG coarsening,
248 : * if the solver supports it (e.g., ML, PETSc's GAMG).
249 : */
250 : void (* nearnullspace) (std::vector<NumericVector<Number> *> & sp, sys_type & S);
251 :
252 : /**
253 : * A callable object that computes a basis for the Jacobian's near nullspace --
254 : * the set of "low energy modes" -- that can be used for AMG coarsening,
255 : * if the solver supports it (e.g., ML, PETSc's GAMG).
256 : */
257 : NonlinearImplicitSystem::ComputeVectorSubspace * nearnullspace_object;
258 :
259 : /**
260 : * Customizable function pointer which users can attach to the
261 : * solver. Gets called prior to every call to solve().
262 : */
263 : void (* user_presolve)(sys_type & S);
264 :
265 : /**
266 : * Function that performs a "check" on the Newton search direction
267 : * and solution after each nonlinear step. See documentation for the
268 : * NonlinearImplicitSystem::ComputePostCheck object for more
269 : * information about the calling sequence.
270 : */
271 : void (* postcheck) (const NumericVector<Number> & old_soln,
272 : NumericVector<Number> & search_direction,
273 : NumericVector<Number> & new_soln,
274 : bool & changed_search_direction,
275 : bool & changed_new_soln,
276 : sys_type & S);
277 :
278 : /**
279 : * A callable object that is executed after each nonlinear
280 : * iteration. Allows the user to modify both the search direction
281 : * and the solution vector in an application-specific way.
282 : */
283 : NonlinearImplicitSystem::ComputePostCheck * postcheck_object;
284 :
285 : NonlinearImplicitSystem::ComputePreCheck * precheck_object;
286 :
287 : /**
288 : * \returns A constant reference to the system we are solving.
289 : */
290 0 : const sys_type & system () const { return _system; }
291 :
292 : /**
293 : * \returns A writable reference to the system we are solving.
294 : */
295 407241 : sys_type & system () { return _system; }
296 :
297 : /**
298 : * Attaches a Preconditioner object to be used during the linear solves.
299 : */
300 : void attach_preconditioner(Preconditioner<T> * preconditioner);
301 :
302 : /**
303 : * Maximum number of non-linear iterations.
304 : */
305 : unsigned int max_nonlinear_iterations;
306 :
307 : /**
308 : * Maximum number of function evaluations.
309 : */
310 : unsigned int max_function_evaluations;
311 :
312 : /**
313 : * The NonlinearSolver should exit after the residual is
314 : * reduced to either less than absolute_residual_tolerance
315 : * or less than relative_residual_tolerance times the
316 : * initial residual.
317 : *
318 : * Users should increase any of these tolerances that they want to use for a
319 : * stopping condition.
320 : *
321 : */
322 : double absolute_residual_tolerance;
323 : double relative_residual_tolerance;
324 :
325 : /**
326 : * The NonlinearSolver should exit if the residual becomes greater
327 : * than the initial residual times the divergence_tolerance.
328 : *
329 : * Users should adjust this tolerances to prevent divergence of the
330 : * NonlinearSolver.
331 : */
332 : double divergence_tolerance;
333 :
334 : /**
335 : * The NonlinearSolver should exit after the full nonlinear step norm is
336 : * reduced to either less than absolute_step_tolerance
337 : * or less than relative_step_tolerance times the largest
338 : * nonlinear solution which has been seen so far.
339 : *
340 : * Users should increase any of these tolerances that they want to use for a
341 : * stopping condition.
342 : *
343 : * \note Not all NonlinearSolvers support \p relative_step_tolerance!
344 : */
345 : double absolute_step_tolerance;
346 : double relative_step_tolerance;
347 :
348 : /**
349 : * Each linear solver step should exit after \p max_linear_iterations
350 : * is exceeded.
351 : */
352 : unsigned int max_linear_iterations;
353 :
354 : /**
355 : * Any required linear solves will at first be done with this tolerance;
356 : * the NonlinearSolver may tighten the tolerance for later solves.
357 : */
358 : double initial_linear_tolerance;
359 :
360 : /**
361 : * The tolerance for linear solves is kept above this minimum
362 : */
363 : double minimum_linear_tolerance;
364 :
365 : /**
366 : * After a call to solve this will reflect whether or not the nonlinear
367 : * solve was successful.
368 : */
369 : bool converged;
370 :
371 : /**
372 : * Set the solver configuration object.
373 : */
374 : void set_solver_configuration(SolverConfiguration & solver_configuration);
375 :
376 : /**
377 : * Get the reuse_preconditioner flag
378 : */
379 : virtual bool reuse_preconditioner() const;
380 :
381 : /**
382 : * Set the reuse preconditioner flag
383 : */
384 : virtual void set_reuse_preconditioner(bool reuse);
385 :
386 : /**
387 : * Get the reuse_preconditioner_max_linear_its parameter
388 : */
389 : virtual unsigned int reuse_preconditioner_max_linear_its() const;
390 :
391 : /**
392 : * Set the reuse_preconditioner_max_linear_its parameter
393 : */
394 : virtual void set_reuse_preconditioner_max_linear_its(unsigned int i);
395 :
396 : /**
397 : * Immediately force a new preconditioner
398 : */
399 0 : virtual void force_new_preconditioner() {};
400 :
401 : /**
402 : * Enable (or disable; it is \p true by default) exact enforcement
403 : * of constraints at the solver level, correcting any constrained
404 : * DoF coefficients in \p current_local_solution as well as applying
405 : * nonlinear residual and Jacobian terms based on constraint
406 : * equations.
407 : *
408 : * This is probably only safe to disable if user code is setting
409 : * nonlinear residual and Jacobian terms based on constraint
410 : * equations at an element-by-element level, by combining the
411 : * \p asymmetric_constraint_rows option with the
412 : * \p residual_constrain_element_vector processing option in
413 : * \p DofMap.
414 : */
415 0 : virtual void set_exact_constraint_enforcement(bool enable)
416 : {
417 0 : _exact_constraint_enforcement = enable;
418 0 : }
419 :
420 0 : bool exact_constraint_enforcement()
421 : {
422 0 : return _exact_constraint_enforcement;
423 : }
424 :
425 : protected:
426 : /**
427 : * Whether we should reuse the linear preconditioner
428 : */
429 : bool _reuse_preconditioner;
430 :
431 : /**
432 : * Whether we should enforce exact constraints globally during a
433 : * solve.
434 : */
435 : bool _exact_constraint_enforcement;
436 :
437 : /**
438 : * Number of linear iterations to retain the preconditioner
439 : */
440 : unsigned int _reuse_preconditioner_max_linear_its;
441 :
442 : /**
443 : * A reference to the system we are solving.
444 : */
445 : sys_type & _system;
446 :
447 : /**
448 : * Flag indicating if the data structures have been initialized.
449 : */
450 : bool _is_initialized;
451 :
452 : /**
453 : * Holds the Preconditioner object to be used for the linear solves.
454 : */
455 : Preconditioner<T> * _preconditioner;
456 :
457 : /**
458 : * Optionally store a SolverOptions object that can be used
459 : * to set parameters like solver type, tolerances and iteration limits.
460 : */
461 : SolverConfiguration * _solver_configuration;
462 : };
463 :
464 :
465 :
466 :
467 : /*----------------------- inline functions ----------------------------------*/
468 : template <typename T>
469 : inline
470 1470 : NonlinearSolver<T>::NonlinearSolver (sys_type & s) :
471 : ParallelObject (s),
472 1386 : residual (nullptr),
473 1386 : residual_object (nullptr),
474 1386 : fd_residual_object (nullptr),
475 1386 : mffd_residual_object (nullptr),
476 1386 : jacobian (nullptr),
477 1386 : jacobian_object (nullptr),
478 1386 : matvec (nullptr),
479 1386 : residual_and_jacobian_object (nullptr),
480 1386 : bounds (nullptr),
481 1386 : bounds_object (nullptr),
482 1386 : nullspace (nullptr),
483 1386 : nullspace_object (nullptr),
484 1386 : transpose_nullspace (nullptr),
485 1386 : transpose_nullspace_object (nullptr),
486 1386 : nearnullspace (nullptr),
487 1386 : nearnullspace_object (nullptr),
488 1386 : user_presolve (nullptr),
489 1386 : postcheck (nullptr),
490 1386 : postcheck_object (nullptr),
491 1386 : precheck_object (nullptr),
492 1386 : max_nonlinear_iterations(0),
493 1386 : max_function_evaluations(0),
494 1386 : absolute_residual_tolerance(0),
495 1386 : relative_residual_tolerance(0),
496 1386 : divergence_tolerance(0),
497 1386 : absolute_step_tolerance(0),
498 1386 : relative_step_tolerance(0),
499 1386 : max_linear_iterations(0),
500 1386 : initial_linear_tolerance(0),
501 1386 : minimum_linear_tolerance(0),
502 1386 : converged(false),
503 1386 : _reuse_preconditioner(false),
504 1386 : _exact_constraint_enforcement(true),
505 1386 : _reuse_preconditioner_max_linear_its(0),
506 1386 : _system(s),
507 1386 : _is_initialized (false),
508 1386 : _preconditioner (nullptr),
509 1470 : _solver_configuration(nullptr)
510 : {
511 1470 : }
512 :
513 :
514 :
515 : template <typename T>
516 : inline
517 42 : NonlinearSolver<T>::~NonlinearSolver ()
518 : {
519 42 : this->NonlinearSolver::clear ();
520 42 : }
521 :
522 :
523 : } // namespace libMesh
524 :
525 :
526 : #endif // LIBMESH_NONLINEAR_SOLVER_H
|
