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 1466 : 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 : /**
145 : * Assembles a residual in \p rhs and/or a jacobian in \p matrix,
146 : * as requested.
147 : *
148 : * This is undefined in ImplicitSystem; subclasses each have their
149 : * own way of handling assembly.
150 : */
151 0 : virtual void assembly (bool /* get_residual */,
152 : bool /* get_jacobian */,
153 : bool /* apply_heterogeneous_constraints */ = false,
154 : bool /* apply_no_constraints */ = false)
155 0 : { libmesh_not_implemented(); }
156 :
157 : /**
158 : * Residual parameter derivative function.
159 : *
160 : * Uses finite differences by default.
161 : *
162 : * This will assemble the sensitivity rhs vectors to hold
163 : * -(partial R / partial p_i), making them ready to solve
164 : * the forward sensitivity equation.
165 : *
166 : * Can be overridden in derived classes.
167 : */
168 : virtual void assemble_residual_derivatives (const ParameterVector & parameters) override;
169 :
170 : /**
171 : * For explicit systems, just assemble and solve the system A*x=b.
172 : * Should be overridden in derived systems to provide a solver for the
173 : * system.
174 : */
175 0 : virtual void solve () override
176 0 : { libmesh_not_implemented(); }
177 :
178 : /**
179 : * Assembles & solves the linear system(s) (dR/du)*u_p = -dR/dp, for
180 : * those parameters contained within \p parameters.
181 : *
182 : * \returns A pair with the total number of linear iterations
183 : * performed and the (sum of the) final residual norms
184 : */
185 : virtual std::pair<unsigned int, Real>
186 : sensitivity_solve (const ParameterVector & parameters) override;
187 :
188 : /**
189 : * Assembles & solves the linear system(s) (dR/du)*u_w = sum(w_p*-dR/dp), for
190 : * those parameters p contained within \p parameters weighted by the
191 : * values w_p found within \p weights.
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 : weighted_sensitivity_solve (const ParameterVector & parameters,
198 : const ParameterVector & weights) override;
199 :
200 : /**
201 : * Assembles & solves the linear system (dR/du)^T*z = dq/du, for
202 : * those quantities of interest q specified by \p qoi_indices.
203 : *
204 : * Leave \p qoi_indices empty to solve all adjoint problems.
205 : *
206 : * \returns A pair with the total number of linear iterations
207 : * performed and the (sum of the) final residual norms
208 : */
209 : virtual std::pair<unsigned int, Real>
210 : adjoint_solve (const QoISet & qoi_indices = QoISet()) override;
211 :
212 : /**
213 : * Assembles & solves the linear system(s)
214 : * (dR/du)^T*z_w = sum(w_p*(d^2q/dudp - d^2R/dudp*z)), for those
215 : * parameters p contained within \p parameters, weighted by the
216 : * values w_p found within \p weights.
217 : *
218 : * Assumes that adjoint_solve has already calculated z for each qoi
219 : * in \p qoi_indices.
220 : *
221 : * \returns A pair with the total number of linear iterations
222 : * performed and the (sum of the) final residual norms
223 : */
224 : virtual std::pair<unsigned int, Real>
225 : weighted_sensitivity_adjoint_solve (const ParameterVector & parameters,
226 : const ParameterVector & weights,
227 : const QoISet & qoi_indices = QoISet()) override;
228 :
229 : /**
230 : * Solves for the derivative of each of the system's quantities of
231 : * interest q in \p qoi[qoi_indices] with respect to each parameter in
232 : * \p parameters, placing the result for qoi \p i and parameter \p j
233 : * into \p sensitivities[i][j].
234 : *
235 : * Uses adjoint_solve() and the adjoint sensitivity method.
236 : *
237 : * Currently uses finite differenced derivatives (partial q /
238 : * partial p) and (partial R / partial p).
239 : */
240 : virtual void adjoint_qoi_parameter_sensitivity (const QoISet & qoi_indices,
241 : const ParameterVector & parameters,
242 : SensitivityData & sensitivities) override;
243 :
244 : /**
245 : * Solves for the derivative of each of the system's quantities of
246 : * interest q in \p qoi[qoi_indices] with respect to each parameter in
247 : * \p parameters, placing the result for qoi \p i and parameter \p j
248 : * into \p sensitivities[i][j].
249 : *
250 : * Uses the forward sensitivity method.
251 : *
252 : * Currently uses finite differenced derivatives (partial q /
253 : * partial p) and (partial R / partial p).
254 : */
255 : virtual void forward_qoi_parameter_sensitivity (const QoISet & qoi_indices,
256 : const ParameterVector & parameters,
257 : SensitivityData & sensitivities) override;
258 :
259 : /**
260 : * For each of the system's quantities of interest q in
261 : * \p qoi[qoi_indices], and for a vector of parameters p, the
262 : * parameter sensitivity Hessian H_ij is defined as
263 : * H_ij = (d^2 q)/(d p_i d p_j)
264 : * This Hessian is the output of this method, where for each q_i,
265 : * H_jk is stored in \p hessian.second_derivative(i,j,k).
266 : *
267 : * Note that in some cases only
268 : * \link current_local_solution \endlink is used during assembly,
269 : * and, therefore, if \link solution \endlink has been altered
270 : * without \link update() \endlink being called, then the
271 : * user must call \link update() \endlink before calling
272 : * this function.
273 : */
274 : virtual void qoi_parameter_hessian(const QoISet & qoi_indices,
275 : const ParameterVector & parameters,
276 : SensitivityData & hessian) override;
277 :
278 : /**
279 : * For each of the system's quantities of interest q in
280 : * \p qoi[qoi_indices], and for a vector of parameters p, the
281 : * parameter sensitivity Hessian H_ij is defined as
282 : * H_ij = (d^2 q)/(d p_i d p_j)
283 : * The Hessian-vector product, for a vector v_k in parameter space, is
284 : * S_j = H_jk v_k
285 : * This product is the output of this method, where for each q_i,
286 : * S_j is stored in \p sensitivities[i][j].
287 : */
288 : virtual void qoi_parameter_hessian_vector_product(const QoISet & qoi_indices,
289 : const ParameterVector & parameters,
290 : const ParameterVector & vector,
291 : SensitivityData & product) override;
292 :
293 : /**
294 : * \returns A const reference to the system's primary matrix.
295 : */
296 : const SparseMatrix<Number> & get_system_matrix() const;
297 :
298 : /**
299 : * \returns A reference to the system's primary matrix.
300 : */
301 : SparseMatrix<Number> & get_system_matrix();
302 :
303 : /**
304 : * The system matrix. Implicit systems are characterized by
305 : * the need to solve the linear system Ax=b. This is the
306 : * system matrix A.
307 : *
308 : * Public access to this member variable will be deprecated in
309 : * the future! Use get_system_matrix() instead.
310 : */
311 : SparseMatrix<Number> * matrix;
312 :
313 : /**
314 : * By default, the system will zero out the matrix and the right hand side.
315 : * If this flag is false, it is the responsibility of the client code
316 : * to take care of setting these to zero before assembly begins
317 : */
318 : bool zero_out_matrix_and_rhs;
319 :
320 : /**
321 : * This class handles all the details of interfacing with various
322 : * linear algebra packages like PETSc or LASPACK. This is a public
323 : * member for backwards compatibility reasons, but in general it's
324 : * better to use get_linear_solver() to access this member, since
325 : * that function will also handle initialization if it hasn't
326 : * already been taken care of.
327 : */
328 : mutable std::unique_ptr<LinearSolver<Number>> linear_solver;
329 :
330 : virtual void create_static_condensation() override;
331 :
332 : /**
333 : * \returns The static condensation system matrix
334 : */
335 : StaticCondensation & get_static_condensation();
336 :
337 : protected:
338 : /**
339 : * Adds the system matrix
340 : */
341 : virtual void add_matrices() override;
342 :
343 : /**
344 : * Sets up the static condensation preconditioner for the supplied \p solver
345 : */
346 : template <typename T>
347 : void setup_static_condensation_preconditioner(T & solver);
348 :
349 : private:
350 : /**
351 : * Create the static condensation system matrix
352 : */
353 : void create_static_condensation_system_matrix();
354 :
355 : /**
356 : * The system matrix for static condensation problems
357 : */
358 : StaticCondensation * _sc_system_matrix;
359 : };
360 :
361 : inline
362 : StaticCondensation & ImplicitSystem::get_static_condensation()
363 : {
364 : libmesh_assert(_sc_system_matrix);
365 : return *_sc_system_matrix;
366 : }
367 :
368 : } // namespace libMesh
369 :
370 : #endif // LIBMESH_IMPLICIT_SYSTEM_H
|
