1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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 #ifdef LIBMESH_FORWARD_DECLARE_ENUMS
31 namespace libMesh
32 {
33 enum SolverPackage : int;
34 }
35 #else
36 #include "libmesh/enum_solver_package.h"
37 #endif
38
39 // C++ includes
40 #include <cstddef>
41 #include <memory>
42
43 namespace libMesh
44 {
45
46 // forward declarations
47 template <typename T> class SparseMatrix;
48 template <typename T> class NumericVector;
49 template <typename T> class Preconditioner;
50 class SolverConfiguration;
51
52 /**
53 * This base class can be inherited from to provide interfaces to
54 * nonlinear solvers from different packages like PETSc and Trilinos.
55 *
56 * \author Benjamin Kirk
57 * \date 2005
58 */
59 template <typename T>
60 class NonlinearSolver : public ReferenceCountedObject<NonlinearSolver<T>>,
61 public ParallelObject
62 {
63 public:
64 /**
65 * The type of system
66 */
67 typedef NonlinearImplicitSystem sys_type;
68
69 /**
70 * Constructor. Initializes Solver data structures
71 */
72 explicit
73 NonlinearSolver (sys_type & s);
74
75 /**
76 * Destructor.
77 */
78 virtual ~NonlinearSolver ();
79
80 /**
81 * Builds a \p NonlinearSolver using the nonlinear solver package specified by
82 * \p solver_package
83 */
84 static std::unique_ptr<NonlinearSolver<T>> build(sys_type & s,
85 const SolverPackage solver_package = libMesh::default_solver_package());
86
87 /**
88 * \returns \p true if the data structures are
89 * initialized, false otherwise.
90 */
initialized()91 bool initialized () const { return _is_initialized; }
92
93 /**
94 * Release all memory and clear data structures.
95 */
clear()96 virtual void clear () {}
97
98 /**
99 * Initialize data structures if not done so already.
100 * May assign a name to the solver in some implementations
101 */
102 virtual void init (const char * name = nullptr) = 0;
103
104 /**
105 * Solves the nonlinear system.
106 */
107 virtual std::pair<unsigned int, Real> solve (SparseMatrix<T> &, // System Jacobian Matrix
108 NumericVector<T> &, // Solution vector
109 NumericVector<T> &, // Residual vector
110 const double, // Stopping tolerance
111 const unsigned int) = 0; // N. Iterations
112
113 /**
114 * Prints a useful message about why the latest nonlinear solve
115 * con(di)verged.
116 */
print_converged_reason()117 virtual void print_converged_reason() { libmesh_not_implemented(); }
118
119 /**
120 * Get the total number of linear iterations done in the last solve
121 */
122 virtual int get_total_linear_iterations() = 0;
123
124 /**
125 * \returns The current nonlinear iteration number if called
126 * *during* the solve(), for example by the user-specified residual
127 * or Jacobian function.
128 *
129 * Must be overridden in derived classes.
130 */
131 virtual unsigned get_current_nonlinear_iteration_number() const = 0;
132
133 /**
134 * Function that computes the residual \p R(X) of the nonlinear system
135 * at the input iterate \p X.
136 */
137 void (* residual) (const NumericVector<Number> & X,
138 NumericVector<Number> & R,
139 sys_type & S);
140
141 /**
142 * Object that computes the residual \p R(X) of the nonlinear system
143 * at the input iterate \p X.
144 */
145 NonlinearImplicitSystem::ComputeResidual * residual_object;
146
147 /**
148 * Object that computes the residual \p R(X) of the nonlinear system
149 * at the input iterate \p X for the purpose of forming a finite-differenced Jacobian.
150 */
151 NonlinearImplicitSystem::ComputeResidual * fd_residual_object;
152
153 /**
154 * Object that computes the residual \p R(X) of the nonlinear system
155 * at the input iterate \p X for the purpose of forming Jacobian-vector products
156 * via finite differencing.
157 */
158 NonlinearImplicitSystem::ComputeResidual * mffd_residual_object;
159
160 /**
161 * Function that computes the Jacobian \p J(X) of the nonlinear system
162 * at the input iterate \p X.
163 */
164 void (* jacobian) (const NumericVector<Number> & X,
165 SparseMatrix<Number> & J,
166 sys_type & S);
167
168 /**
169 * Object that computes the Jacobian \p J(X) of the nonlinear system
170 * at the input iterate \p X.
171 */
172 NonlinearImplicitSystem::ComputeJacobian * jacobian_object;
173
174 /**
175 * Function that computes either the residual \f$ R(X) \f$ or the
176 * Jacobian \f$ J(X) \f$ of the nonlinear system at the input
177 * iterate \f$ X \f$.
178 *
179 * \note Either \p R or \p J could be \p nullptr.
180 */
181 void (* matvec) (const NumericVector<Number> & X,
182 NumericVector<Number> * R,
183 SparseMatrix<Number> * J,
184 sys_type & S);
185
186 /**
187 * Object that computes either the residual \f$ R(X) \f$ or the
188 * Jacobian \f$ J(X) \f$ of the nonlinear system at the input
189 * iterate \f$ X \f$.
190 *
191 * \note Either \p R or \p J could be \p nullptr.
192 */
193 NonlinearImplicitSystem::ComputeResidualandJacobian * residual_and_jacobian_object;
194
195 /**
196 * Function that computes the lower and upper bounds \p XL and \p XU on the solution of the nonlinear system.
197 */
198 void (* bounds) (NumericVector<Number> & XL,
199 NumericVector<Number> & XU,
200 sys_type & S);
201 /**
202 * Object that computes the bounds vectors \f$ XL \f$ and \f$ XU \f$.
203 */
204 NonlinearImplicitSystem::ComputeBounds * bounds_object;
205
206 /**
207 * Function 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 void (* nullspace) (std::vector<NumericVector<Number> *> & sp, sys_type & S);
213
214 /**
215 * A callable object that computes a basis for the Jacobian's nullspace --
216 * the kernel or the "zero energy modes" -- that can be used in
217 * solving a degenerate problem iteratively, if the solver supports it
218 * (e.g., PETSc's KSP).
219 */
220 NonlinearImplicitSystem::ComputeVectorSubspace * nullspace_object;
221
222 /**
223 * Function that computes a basis for the transpose Jacobian's nullspace --
224 * when solving a degenerate problem iteratively, if the solver supports it
225 * (e.g., PETSc's KSP), it is used to remove contributions outside of R(jac)
226 */
227 void (* transpose_nullspace) (std::vector<NumericVector<Number> *> & sp, sys_type & S);
228
229 /**
230 * A callable object 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 NonlinearImplicitSystem::ComputeVectorSubspace * transpose_nullspace_object;
235
236 /**
237 * Function that computes a basis for the Jacobian's near nullspace --
238 * the set of "low energy modes" -- that can be used for AMG coarsening,
239 * if the solver supports it (e.g., ML, PETSc's GAMG).
240 */
241 void (* nearnullspace) (std::vector<NumericVector<Number> *> & sp, sys_type & S);
242
243 /**
244 * A callable object 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 NonlinearImplicitSystem::ComputeVectorSubspace * nearnullspace_object;
249
250 /**
251 * Customizable function pointer which users can attach to the
252 * solver. Gets called prior to every call to solve().
253 */
254 void (* user_presolve)(sys_type & S);
255
256 /**
257 * Function that performs a "check" on the Newton search direction
258 * and solution after each nonlinear step. See documentation for the
259 * NonlinearImplicitSystem::ComputePostCheck object for more
260 * information about the calling sequence.
261 */
262 void (* postcheck) (const NumericVector<Number> & old_soln,
263 NumericVector<Number> & search_direction,
264 NumericVector<Number> & new_soln,
265 bool & changed_search_direction,
266 bool & changed_new_soln,
267 sys_type & S);
268
269 /**
270 * A callable object that is executed after each nonlinear
271 * iteration. Allows the user to modify both the search direction
272 * and the solution vector in an application-specific way.
273 */
274 NonlinearImplicitSystem::ComputePostCheck * postcheck_object;
275
276 /**
277 * \returns A constant reference to the system we are solving.
278 */
system()279 const sys_type & system () const { return _system; }
280
281 /**
282 * \returns A writable reference to the system we are solving.
283 */
system()284 sys_type & system () { return _system; }
285
286 /**
287 * Attaches a Preconditioner object to be used during the linear solves.
288 */
289 void attach_preconditioner(Preconditioner<T> * preconditioner);
290
291 /**
292 * Maximum number of non-linear iterations.
293 */
294 unsigned int max_nonlinear_iterations;
295
296 /**
297 * Maximum number of function evaluations.
298 */
299 unsigned int max_function_evaluations;
300
301 /**
302 * The NonlinearSolver should exit after the residual is
303 * reduced to either less than absolute_residual_tolerance
304 * or less than relative_residual_tolerance times the
305 * initial residual.
306 *
307 * Users should increase any of these tolerances that they want to use for a
308 * stopping condition.
309 *
310 */
311 double absolute_residual_tolerance;
312 double relative_residual_tolerance;
313
314 /**
315 * The NonlinearSolver should exit if the residual becomes greater
316 * than the initial residual times the divergence_tolerance.
317 *
318 * Users should adjust this tolerances to prevent divergence of the
319 * NonlinearSolver.
320 */
321 double divergence_tolerance;
322
323 /**
324 * The NonlinearSolver should exit after the full nonlinear step norm is
325 * reduced to either less than absolute_step_tolerance
326 * or less than relative_step_tolerance times the largest
327 * nonlinear solution which has been seen so far.
328 *
329 * Users should increase any of these tolerances that they want to use for a
330 * stopping condition.
331 *
332 * \note Not all NonlinearSolvers support \p relative_step_tolerance!
333 */
334 double absolute_step_tolerance;
335 double relative_step_tolerance;
336
337 /**
338 * Each linear solver step should exit after \p max_linear_iterations
339 * is exceeded.
340 */
341 unsigned int max_linear_iterations;
342
343 /**
344 * Any required linear solves will at first be done with this tolerance;
345 * the NonlinearSolver may tighten the tolerance for later solves.
346 */
347 double initial_linear_tolerance;
348
349 /**
350 * The tolerance for linear solves is kept above this minimum
351 */
352 double minimum_linear_tolerance;
353
354 /**
355 * After a call to solve this will reflect whether or not the nonlinear
356 * solve was successful.
357 */
358 bool converged;
359
360 /**
361 * Set the solver configuration object.
362 */
363 void set_solver_configuration(SolverConfiguration & solver_configuration);
364
365 protected:
366 /**
367 * A reference to the system we are solving.
368 */
369 sys_type & _system;
370
371 /**
372 * Flag indicating if the data structures have been initialized.
373 */
374 bool _is_initialized;
375
376 /**
377 * Holds the Preconditioner object to be used for the linear solves.
378 */
379 Preconditioner<T> * _preconditioner;
380
381 /**
382 * Optionally store a SolverOptions object that can be used
383 * to set parameters like solver type, tolerances and iteration limits.
384 */
385 SolverConfiguration * _solver_configuration;
386 };
387
388
389
390
391 /*----------------------- inline functions ----------------------------------*/
392 template <typename T>
393 inline
NonlinearSolver(sys_type & s)394 NonlinearSolver<T>::NonlinearSolver (sys_type & s) :
395 ParallelObject (s),
396 residual (nullptr),
397 residual_object (nullptr),
398 fd_residual_object (nullptr),
399 mffd_residual_object (nullptr),
400 jacobian (nullptr),
401 jacobian_object (nullptr),
402 matvec (nullptr),
403 residual_and_jacobian_object (nullptr),
404 bounds (nullptr),
405 bounds_object (nullptr),
406 nullspace (nullptr),
407 nullspace_object (nullptr),
408 transpose_nullspace (nullptr),
409 transpose_nullspace_object (nullptr),
410 nearnullspace (nullptr),
411 nearnullspace_object (nullptr),
412 user_presolve (nullptr),
413 postcheck (nullptr),
414 postcheck_object (nullptr),
415 max_nonlinear_iterations(0),
416 max_function_evaluations(0),
417 absolute_residual_tolerance(0),
418 relative_residual_tolerance(0),
419 divergence_tolerance(0),
420 absolute_step_tolerance(0),
421 relative_step_tolerance(0),
422 max_linear_iterations(0),
423 initial_linear_tolerance(0),
424 minimum_linear_tolerance(0),
425 converged(false),
426 _system(s),
427 _is_initialized (false),
428 _preconditioner (nullptr),
429 _solver_configuration(nullptr)
430 {
431 }
432
433
434
435 template <typename T>
436 inline
~NonlinearSolver()437 NonlinearSolver<T>::~NonlinearSolver ()
438 {
439 this->clear ();
440 }
441
442
443 } // namespace libMesh
444
445
446 #endif // LIBMESH_NONLINEAR_SOLVER_H
447