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