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