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_IMPLICIT_SYSTEM_H
21 : #define LIBMESH_IMPLICIT_SYSTEM_H
22 :
23 : // Local Includes
24 : #include "libmesh/explicit_system.h"
25 :
26 : // C++ includes
27 : #include <cstddef>
28 : #include <memory>
29 :
30 : namespace libMesh
31 : {
32 :
33 : // Forward declarations
34 : template <typename T> class LinearSolver;
35 : class StaticCondensation;
36 :
37 : /**
38 : * \brief Manages consistently variables, degrees of freedom, coefficient
39 : * vectors, and matrices for implicit systems.
40 : *
41 : * Implicit systems are characterized by the need to solve the (non-)linear
42 : * system Ax=b. This class provides, in addition to the ExplicitSystem class,
43 : * storage for sparse matrices. Hence this System provides means to manage a
44 : * right hand side vector and a sparse matrix. In addition, further matrices
45 : * can be managed using the ImplicitSystem::add_matrix method.
46 : *
47 : * This class provieds *no* means to solve the implicit system. This
48 : * functionality is provided, e.g., by the LinearImplicitSystem or the
49 : * NonlinearImplicitSystem class.
50 : *
51 : * \note Additional vectors/matrices can be added via parent class
52 : * interfaces.
53 : *
54 : * \author Benjamin S. Kirk
55 : * \date 2004
56 : */
57 1378 : class ImplicitSystem : public ExplicitSystem
58 : {
59 : public:
60 :
61 : /**
62 : * Constructor.
63 : */
64 : ImplicitSystem (EquationSystems & es,
65 : const std::string & name,
66 : const unsigned int number);
67 :
68 : /**
69 : * Special functions.
70 : * - This class has the same restrictions/defaults as its base class.
71 : * - The destructor is defaulted out-of-line.
72 : */
73 : ImplicitSystem (const ImplicitSystem &) = delete;
74 : ImplicitSystem & operator= (const ImplicitSystem &) = delete;
75 : ImplicitSystem (ImplicitSystem &&) = default;
76 : ImplicitSystem & operator= (ImplicitSystem &&) = delete;
77 : virtual ~ImplicitSystem ();
78 :
79 : /**
80 : * The type of system.
81 : */
82 : typedef ImplicitSystem sys_type;
83 :
84 : /**
85 : * \returns A reference to *this.
86 : */
87 : sys_type & system () { return *this; }
88 :
89 : /**
90 : * The type of the parent.
91 : */
92 : typedef ExplicitSystem Parent;
93 :
94 : /**
95 : * Clear all the data structures associated with
96 : * the system.
97 : */
98 : virtual void clear () override;
99 :
100 : /**
101 : * Prepares \p matrix and \p rhs for system assembly, then calls
102 : * user assembly function.
103 : * Can be overridden in derived classes.
104 : */
105 : virtual void assemble () override;
106 :
107 : /**
108 : * Avoids use of any cached data that might affect any solve result. Should
109 : * be overridden in derived systems.
110 : */
111 : virtual void disable_cache () override;
112 :
113 : /**
114 : * \returns \p "Implicit". Helps in identifying
115 : * the system type in an equation system file.
116 : */
117 0 : virtual std::string system_type () const override { return "Implicit"; }
118 :
119 : /**
120 : * \returns A dumb pointer to a local std::unique_ptr<LinearSolver>
121 : * member which can be used in adjoint and/or sensitivity solves.
122 : *
123 : * To mimic the previous behavior of this function, if the
124 : * linear_solver member is already initialized when this function is
125 : * called, this function first clears it. That is, no attempt is
126 : * made to reuse an existing LinearSolver object. The user MUST NOT
127 : * attempt to clean up this pointer, otherwise the std::unique_ptr
128 : * held by this class will likely become corrupted.
129 : *
130 : * This function is virtual so it can be overridden in derived
131 : * classes if necessary, however most will probably just want to use
132 : * the base class behavior.
133 : */
134 : virtual LinearSolver<Number> * get_linear_solver() const;
135 :
136 : /**
137 : * \returns An integer corresponding to the upper iteration count
138 : * limit and a Real corresponding to the convergence tolerance to
139 : * be used in linear adjoint and/or sensitivity solves
140 : */
141 : virtual std::pair<unsigned int, Real>
142 : get_linear_solve_parameters() const;
143 :
144 : #ifdef LIBMESH_ENABLE_DEPRECATED
145 : /**
146 : * Currently a no-op.
147 : *
148 : * \deprecated This function no longer needs to be called, since
149 : * get_linear_solver() no longer returns a heap-allocated dumb
150 : * pointer.
151 : */
152 : virtual void release_linear_solver(LinearSolver<Number> *) const;
153 : #endif // LIBMESH_ENABLE_DEPRECATED
154 :
155 : /**
156 : * Assembles a residual in \p rhs and/or a jacobian in \p matrix,
157 : * as requested.
158 : *
159 : * This is undefined in ImplicitSystem; subclasses each have their
160 : * own way of handling assembly.
161 : */
162 0 : virtual void assembly (bool /* get_residual */,
163 : bool /* get_jacobian */,
164 : bool /* apply_heterogeneous_constraints */ = false,
165 : bool /* apply_no_constraints */ = false)
166 0 : { libmesh_not_implemented(); }
167 :
168 : /**
169 : * Residual parameter derivative function.
170 : *
171 : * Uses finite differences by default.
172 : *
173 : * This will assemble the sensitivity rhs vectors to hold
174 : * -(partial R / partial p_i), making them ready to solve
175 : * the forward sensitivity equation.
176 : *
177 : * Can be overridden in derived classes.
178 : */
179 : virtual void assemble_residual_derivatives (const ParameterVector & parameters) override;
180 :
181 : /**
182 : * For explicit systems, just assemble and solve the system A*x=b.
183 : * Should be overridden in derived systems to provide a solver for the
184 : * system.
185 : */
186 0 : virtual void solve () override
187 0 : { libmesh_not_implemented(); }
188 :
189 : /**
190 : * Assembles & solves the linear system(s) (dR/du)*u_p = -dR/dp, for
191 : * those parameters contained within \p parameters.
192 : *
193 : * \returns A pair with the total number of linear iterations
194 : * performed and the (sum of the) final residual norms
195 : */
196 : virtual std::pair<unsigned int, Real>
197 : sensitivity_solve (const ParameterVector & parameters) override;
198 :
199 : /**
200 : * Assembles & solves the linear system(s) (dR/du)*u_w = sum(w_p*-dR/dp), for
201 : * those parameters p contained within \p parameters weighted by the
202 : * values w_p found within \p weights.
203 : *
204 : * \returns A pair with the total number of linear iterations
205 : * performed and the (sum of the) final residual norms
206 : */
207 : virtual std::pair<unsigned int, Real>
208 : weighted_sensitivity_solve (const ParameterVector & parameters,
209 : const ParameterVector & weights) override;
210 :
211 : /**
212 : * Assembles & solves the linear system (dR/du)^T*z = dq/du, for
213 : * those quantities of interest q specified by \p qoi_indices.
214 : *
215 : * Leave \p qoi_indices empty to solve all adjoint problems.
216 : *
217 : * \returns A pair with the total number of linear iterations
218 : * performed and the (sum of the) final residual norms
219 : */
220 : virtual std::pair<unsigned int, Real>
221 : adjoint_solve (const QoISet & qoi_indices = QoISet()) override;
222 :
223 : /**
224 : * Assembles & solves the linear system(s)
225 : * (dR/du)^T*z_w = sum(w_p*(d^2q/dudp - d^2R/dudp*z)), for those
226 : * parameters p contained within \p parameters, weighted by the
227 : * values w_p found within \p weights.
228 : *
229 : * Assumes that adjoint_solve has already calculated z for each qoi
230 : * in \p qoi_indices.
231 : *
232 : * \returns A pair with the total number of linear iterations
233 : * performed and the (sum of the) final residual norms
234 : */
235 : virtual std::pair<unsigned int, Real>
236 : weighted_sensitivity_adjoint_solve (const ParameterVector & parameters,
237 : const ParameterVector & weights,
238 : const QoISet & qoi_indices = QoISet()) override;
239 :
240 : /**
241 : * Solves for the derivative of each of the system's quantities of
242 : * interest q in \p qoi[qoi_indices] with respect to each parameter in
243 : * \p parameters, placing the result for qoi \p i and parameter \p j
244 : * into \p sensitivities[i][j].
245 : *
246 : * Uses adjoint_solve() and the adjoint sensitivity method.
247 : *
248 : * Currently uses finite differenced derivatives (partial q /
249 : * partial p) and (partial R / partial p).
250 : */
251 : virtual void adjoint_qoi_parameter_sensitivity (const QoISet & qoi_indices,
252 : const ParameterVector & parameters,
253 : SensitivityData & sensitivities) override;
254 :
255 : /**
256 : * Solves for the derivative of each of the system's quantities of
257 : * interest q in \p qoi[qoi_indices] with respect to each parameter in
258 : * \p parameters, placing the result for qoi \p i and parameter \p j
259 : * into \p sensitivities[i][j].
260 : *
261 : * Uses the forward sensitivity method.
262 : *
263 : * Currently uses finite differenced derivatives (partial q /
264 : * partial p) and (partial R / partial p).
265 : */
266 : virtual void forward_qoi_parameter_sensitivity (const QoISet & qoi_indices,
267 : const ParameterVector & parameters,
268 : SensitivityData & sensitivities) override;
269 :
270 : /**
271 : * For each of the system's quantities of interest q in
272 : * \p qoi[qoi_indices], and for a vector of parameters p, the
273 : * parameter sensitivity Hessian H_ij is defined as
274 : * H_ij = (d^2 q)/(d p_i d p_j)
275 : * This Hessian is the output of this method, where for each q_i,
276 : * H_jk is stored in \p hessian.second_derivative(i,j,k).
277 : *
278 : * Note that in some cases only
279 : * \link current_local_solution \endlink is used during assembly,
280 : * and, therefore, if \link solution \endlink has been altered
281 : * without \link update() \endlink being called, then the
282 : * user must call \link update() \endlink before calling
283 : * this function.
284 : */
285 : virtual void qoi_parameter_hessian(const QoISet & qoi_indices,
286 : const ParameterVector & parameters,
287 : SensitivityData & hessian) override;
288 :
289 : /**
290 : * For each of the system's quantities of interest q in
291 : * \p qoi[qoi_indices], and for a vector of parameters p, the
292 : * parameter sensitivity Hessian H_ij is defined as
293 : * H_ij = (d^2 q)/(d p_i d p_j)
294 : * The Hessian-vector product, for a vector v_k in parameter space, is
295 : * S_j = H_jk v_k
296 : * This product is the output of this method, where for each q_i,
297 : * S_j is stored in \p sensitivities[i][j].
298 : */
299 : virtual void qoi_parameter_hessian_vector_product(const QoISet & qoi_indices,
300 : const ParameterVector & parameters,
301 : const ParameterVector & vector,
302 : SensitivityData & product) override;
303 :
304 : /**
305 : * \returns A const reference to the system's primary matrix.
306 : */
307 : const SparseMatrix<Number> & get_system_matrix() const;
308 :
309 : /**
310 : * \returns A reference to the system's primary matrix.
311 : */
312 : SparseMatrix<Number> & get_system_matrix();
313 :
314 : /**
315 : * The system matrix. Implicit systems are characterized by
316 : * the need to solve the linear system Ax=b. This is the
317 : * system matrix A.
318 : *
319 : * Public access to this member variable will be deprecated in
320 : * the future! Use get_system_matrix() instead.
321 : */
322 : SparseMatrix<Number> * matrix;
323 :
324 : /**
325 : * By default, the system will zero out the matrix and the right hand side.
326 : * If this flag is false, it is the responsibility of the client code
327 : * to take care of setting these to zero before assembly begins
328 : */
329 : bool zero_out_matrix_and_rhs;
330 :
331 : /**
332 : * This class handles all the details of interfacing with various
333 : * linear algebra packages like PETSc or LASPACK. This is a public
334 : * member for backwards compatibility reasons, but in general it's
335 : * better to use get_linear_solver() to access this member, since
336 : * that function will also handle initialization if it hasn't
337 : * already been taken care of.
338 : */
339 : mutable std::unique_ptr<LinearSolver<Number>> linear_solver;
340 :
341 : virtual void create_static_condensation() override;
342 :
343 : /**
344 : * \returns The static condensation system matrix
345 : */
346 : StaticCondensation & get_static_condensation();
347 :
348 : protected:
349 : /**
350 : * Adds the system matrix
351 : */
352 : virtual void add_matrices() override;
353 :
354 : /**
355 : * Sets up the static condensation preconditioner for the supplied \p solver
356 : */
357 : template <typename T>
358 : void setup_static_condensation_preconditioner(T & solver);
359 :
360 : private:
361 : /**
362 : * Create the static condensation system matrix
363 : */
364 : void create_static_condensation_system_matrix();
365 :
366 : /**
367 : * The system matrix for static condensation problems
368 : */
369 : StaticCondensation * _sc_system_matrix;
370 : };
371 :
372 : inline
373 : StaticCondensation & ImplicitSystem::get_static_condensation()
374 : {
375 : libmesh_assert(_sc_system_matrix);
376 : return *_sc_system_matrix;
377 : }
378 :
379 : } // namespace libMesh
380 :
381 : #endif // LIBMESH_IMPLICIT_SYSTEM_H
|
