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_IMPLICIT_SYSTEM_H
21 : #define LIBMESH_NONLINEAR_IMPLICIT_SYSTEM_H
22 :
23 : // Local Includes
24 : #include "libmesh/implicit_system.h"
25 :
26 : // C++ includes
27 :
28 : namespace libMesh
29 : {
30 :
31 :
32 : // Forward declarations
33 : class DiffSolver;
34 : template<typename T> class NonlinearSolver;
35 :
36 :
37 : /**
38 : * \brief Manages consistently variables, degrees of freedom, coefficient
39 : * vectors, matrices and non-linear solvers for implicit systems.
40 : *
41 : * An implicit system is a system that requires the solution of a
42 : * system of equations. This class has the ability to create and use a
43 : * non-linear solver to solve the system of equations.
44 : *
45 : * The matrix NonlinearImplicitSystem::matrix and the vector
46 : * NonlinearImplicitSystem::rhs should be filled during assembly.
47 : *
48 : * \note Additional vectors/matrices can be added via parent class
49 : * interfaces.
50 : *
51 : * \author Benjamin S. Kirk
52 : * \date 2005
53 : */
54 194 : class NonlinearImplicitSystem : public ImplicitSystem
55 : {
56 : public:
57 :
58 : /**
59 : * Constructor.
60 : */
61 : NonlinearImplicitSystem (EquationSystems & es,
62 : const std::string & name,
63 : const unsigned int number);
64 :
65 : /**
66 : * Special functions.
67 : * - This class has the same restrictions/defaults as its base class.
68 : * - The destructor is defaulted out-of-line.
69 : */
70 : NonlinearImplicitSystem (const NonlinearImplicitSystem &) = delete;
71 : NonlinearImplicitSystem & operator= (const NonlinearImplicitSystem &) = delete;
72 : NonlinearImplicitSystem (NonlinearImplicitSystem &&) = default;
73 : NonlinearImplicitSystem & operator= (NonlinearImplicitSystem &&) = delete;
74 : virtual ~NonlinearImplicitSystem ();
75 :
76 : /**
77 : * The type of system.
78 : */
79 : typedef NonlinearImplicitSystem sys_type;
80 :
81 : /**
82 : * The type of the parent.
83 : */
84 : typedef ImplicitSystem Parent;
85 :
86 : /**
87 : * Abstract base class to be used to calculate the residual
88 : * of a nonlinear system.
89 : */
90 : class ComputeResidual
91 : {
92 : public:
93 : virtual ~ComputeResidual () = default;
94 : /**
95 : * Residual function. This function will be called to compute the
96 : * residual and must be implemented by the user in a derived class.
97 : */
98 : virtual void residual (const NumericVector<Number> & X,
99 : NumericVector<Number> & R,
100 : sys_type & S) = 0;
101 : };
102 :
103 :
104 : /**
105 : * Abstract base class to be used to calculate the Jacobian
106 : * of a nonlinear system.
107 : */
108 : class ComputeJacobian
109 : {
110 : public:
111 : virtual ~ComputeJacobian () = default;
112 :
113 : /**
114 : * Jacobian function. This function will be called to compute the
115 : * jacobian and must be implemented by the user in a derived class.
116 : */
117 : virtual void jacobian (const NumericVector<Number> & X,
118 : SparseMatrix<Number> & J,
119 : sys_type & S) = 0;
120 : };
121 :
122 :
123 : /**
124 : * Abstract base class to be used to calculate the bounds
125 : * on the degrees of freedom of a nonlinear system.
126 : */
127 : class ComputeBounds
128 : {
129 : public:
130 : virtual ~ComputeBounds () = default;
131 :
132 : /**
133 : * This function will be called to compute the bounds vector and
134 : * must be implemented by the user in a derived class.
135 : */
136 : virtual void bounds (NumericVector<Number> & XL,
137 : NumericVector<Number> & XU,
138 : sys_type & S) = 0;
139 : };
140 :
141 : /**
142 : * Callable abstract base class to be used as a callback to provide
143 : * the solver with a basis for the system's Jacobian's nullspace
144 : * (the kernel or the "zero energy modes") or near-nullspace
145 : * (the "low energy modes").
146 : *
147 : * A nullspace can be used to solve a degenerate problem iteratively
148 : * (e.g., with a Krylov subspace method). A near nullspace can be used
149 : * by an Algebraic Multigrid (AMG) preconditioner to construct smoothed-
150 : * aggregation-like coarse spaces.
151 : */
152 : class ComputeVectorSubspace
153 : {
154 : public:
155 : virtual ~ComputeVectorSubspace () = default;
156 :
157 : /**
158 : * This function will be called to compute the subspace basis
159 : * (e.g., nullspace or nearnullspace).
160 : * It must be implemented by the user in a derived class.
161 : */
162 : virtual void operator()(std::vector<NumericVector<Number> *> & sp,
163 : sys_type & s) = 0;
164 : };
165 :
166 : /**
167 : * Abstract base class to be used to calculate the residual and Jacobian
168 : * simultaneously of a nonlinear system.
169 : */
170 : class ComputeResidualandJacobian
171 : {
172 : public:
173 : virtual ~ComputeResidualandJacobian () = default;
174 :
175 : /**
176 : * Residual & Jacobian function, calculated simultaneously.
177 : * This function will be called to compute the residual and jacobian
178 : * simultaneously and must be implemented by the user in a derived class.
179 : */
180 : virtual void residual_and_jacobian (const NumericVector<Number> & X,
181 : NumericVector<Number> * R,
182 : SparseMatrix<Number> * J,
183 : sys_type & S) = 0;
184 : };
185 :
186 : /**
187 : * Abstract base class to be used for applying user modifications to
188 : * the solution vector and/or Newton update step after each
189 : * nonlinear step. A user should inherit from this class and define
190 : * a custom version of the postcheck() function.
191 : */
192 : class ComputePostCheck
193 : {
194 : public:
195 : virtual ~ComputePostCheck () = default;
196 :
197 : /**
198 : * This interface, which is inspired by PETSc's, passes the user:
199 : * .) A constant reference to the "old" solution vector (previous Newton iterate),
200 : * .) A pointer to the search direction (update) vector, and
201 : * .) The "proposed" new solution vector.
202 : * The user can then modify either the search direction or the new
203 : * solution vector according to their own application-specific
204 : * criteria. If they do, they must set the associated boolean
205 : * references appropriately.
206 : */
207 : virtual void postcheck (const NumericVector<Number> & old_soln,
208 : NumericVector<Number> & search_direction,
209 : NumericVector<Number> & new_soln,
210 : bool & changed_search_direction,
211 : bool & changed_new_soln,
212 : sys_type & S) = 0;
213 : };
214 :
215 : /**
216 : * Abstract base class to be used for applying user modifications to the Newton search direction
217 : * \emph before the solver's line search is called.
218 : */
219 : class ComputePreCheck
220 : {
221 : public:
222 : virtual ~ComputePreCheck () = default;
223 :
224 : /**
225 : * Abstract precheck method that users must override
226 : * @param precheck_soln This is the solution prior to the linesearch (which we have not
227 : * performed yet). This will correspond to the system's \p current_local_solution data member,
228 : * so it should hold all the information locally that a user might want to query from the vector
229 : * (e.g. possibly ghost entries)
230 : * @param search_direction The search direction that will be used in the upcoming line search.
231 : * \emph Note that this is the \emph negative of the solution update, at least when using PETSc
232 : * as the backend. The user may modify this vector
233 : * @param changed Whether the user has modified the \p search_direction
234 : * @param S The nonlinear implicit system
235 : */
236 : virtual void precheck(const NumericVector<Number> & precheck_soln,
237 : NumericVector<Number> & search_direction,
238 : bool & changed,
239 : sys_type & S) = 0;
240 : };
241 :
242 : /**
243 : * \returns A reference to *this.
244 : */
245 : sys_type & system () { return *this; }
246 :
247 : /**
248 : * Clear all the data structures associated with
249 : * the system.
250 : */
251 : virtual void clear () override;
252 :
253 : /**
254 : * Reinitializes the member data fields associated with
255 : * the system, so that, e.g., \p assemble() may be used.
256 : */
257 : virtual void reinit () override;
258 :
259 : /**
260 : * Assembles & solves the nonlinear system R(x) = 0.
261 : */
262 : virtual void solve () override;
263 :
264 : /**
265 : * \returns An integer corresponding to the upper iteration count
266 : * limit and a Real corresponding to the convergence tolerance to
267 : * be used in linear adjoint and/or sensitivity solves
268 : */
269 : virtual std::pair<unsigned int, Real>
270 : get_linear_solve_parameters() const override;
271 :
272 : /**
273 : * Assembles a residual in \p rhs and/or a jacobian in \p matrix,
274 : * as requested.
275 : */
276 : virtual void assembly(bool get_residual,
277 : bool get_jacobian,
278 : bool apply_heterogeneous_constraints = false,
279 : bool apply_no_constraints = false) override;
280 :
281 : /**
282 : * \returns \p "NonlinearImplicit". Helps in identifying
283 : * the system type in an equation system file.
284 : */
285 910 : virtual std::string system_type () const override { return "NonlinearImplicit"; }
286 :
287 : /**
288 : * The \p NonlinearSolver defines the default interface used to
289 : * solve the nonlinear_implicit system. This class handles all the
290 : * details of interfacing with various nonlinear algebra packages
291 : * like PETSc or LASPACK.
292 : */
293 : std::unique_ptr<NonlinearSolver<Number>> nonlinear_solver;
294 :
295 : /**
296 : * The \p DiffSolver defines an optional interface used to
297 : * solve the nonlinear_implicit system.
298 : */
299 : std::unique_ptr<DiffSolver> diff_solver;
300 :
301 : /**
302 : * \returns The number of iterations
303 : * taken for the most recent nonlinear solve.
304 : */
305 : unsigned int n_nonlinear_iterations() const { return _n_nonlinear_iterations; }
306 :
307 : /**
308 : * \returns The final residual for the nonlinear system solve.
309 : */
310 : Real final_nonlinear_residual() const { return _final_nonlinear_residual; }
311 :
312 : /**
313 : * If called *during* the solve(), for example by the user-specified
314 : * residual or Jacobian function, return the current nonlinear iteration
315 : * number.
316 : */
317 : unsigned get_current_nonlinear_iteration_number() const;
318 :
319 : virtual void create_static_condensation() override;
320 :
321 : protected:
322 :
323 : /**
324 : * Copies system parameters into nonlinear solver parameters
325 : */
326 : void set_solver_parameters();
327 :
328 : /**
329 : * The number of nonlinear iterations required to solve the nonlinear
330 : * system R(x)=0.
331 : */
332 : unsigned int _n_nonlinear_iterations;
333 :
334 : /**
335 : * The final residual for the nonlinear system R(x)
336 : */
337 : Real _final_nonlinear_residual;
338 : };
339 :
340 : } // namespace libMesh
341 :
342 : #endif // LIBMESH_NONLINEAR_IMPLICIT_SYSTEM_H
|
