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_DIFF_SYSTEM_H
21 : #define LIBMESH_DIFF_SYSTEM_H
22 :
23 : // Local Includes
24 : #include "libmesh/diff_context.h"
25 : #include "libmesh/diff_physics.h"
26 : #include "libmesh/diff_qoi.h"
27 : #include "libmesh/implicit_system.h"
28 : #include "libmesh/time_solver.h"
29 :
30 : // C++ includes
31 : #include <memory>
32 : #include <stack>
33 :
34 : namespace libMesh
35 : {
36 :
37 : // Forward Declarations
38 : class TimeSolver;
39 :
40 : template <typename T> class NumericVector;
41 :
42 : /**
43 : * This class provides a specific system class. It aims
44 : * to generalize any system, linear or nonlinear, which
45 : * provides both a residual and a Jacobian.
46 : *
47 : * This class is part of the new DifferentiableSystem framework,
48 : * which is still experimental. Users of this framework should
49 : * beware of bugs and future API changes.
50 : *
51 : * \author Roy H. Stogner
52 : * \date 2006
53 : */
54 450 : class DifferentiableSystem : public ImplicitSystem,
55 : public virtual DifferentiablePhysics,
56 : public virtual DifferentiableQoI
57 : {
58 : public:
59 :
60 : /**
61 : * Constructor.
62 : */
63 : DifferentiableSystem (EquationSystems & es,
64 : const std::string & name,
65 : const unsigned int number);
66 :
67 : /**
68 : * Special functions.
69 : * - This class has the same restrictions as its base class.
70 : * - The destructor is defaulted out-of-line.
71 : */
72 : DifferentiableSystem (const DifferentiableSystem &) = delete;
73 : DifferentiableSystem & operator= (const DifferentiableSystem &) = delete;
74 : DifferentiableSystem (DifferentiableSystem &&) = default;
75 : DifferentiableSystem & operator= (DifferentiableSystem &&) = delete;
76 : virtual ~DifferentiableSystem ();
77 :
78 : /**
79 : * The type of system.
80 : */
81 : typedef DifferentiableSystem sys_type;
82 :
83 : /**
84 : * The type of the parent.
85 : */
86 : typedef ImplicitSystem Parent;
87 :
88 : /**
89 : * Clear all the data structures associated with
90 : * the system.
91 : */
92 : virtual void clear () override;
93 :
94 : /**
95 : * Reinitializes the member data fields associated with
96 : * the system, so that, e.g., \p assemble() may be used.
97 : */
98 : virtual void reinit () override;
99 :
100 : /**
101 : * Prepares \p matrix and \p rhs for matrix assembly.
102 : * Users should not reimplement this.
103 : * Note that in some cases only
104 : * \link current_local_solution \endlink is used during assembly,
105 : * and, therefore, if \link solution \endlink has been altered
106 : * without \link update() \endlink being called, then the
107 : * user must call \link update() \endlink before calling
108 : * this function.
109 : */
110 : virtual void assemble () override;
111 :
112 : /**
113 : * \returns A pointer to a linear solver appropriate for use in
114 : * adjoint and/or sensitivity solves
115 : */
116 : virtual LinearSolver<Number> * get_linear_solver() const override;
117 :
118 : /**
119 : * \returns An integer corresponding to the upper iteration count
120 : * limit and a Real corresponding to the convergence tolerance to
121 : * be used in linear adjoint and/or sensitivity solves
122 : */
123 : virtual std::pair<unsigned int, Real>
124 : get_linear_solve_parameters() const override;
125 :
126 : /**
127 : * Assembles a residual in \p rhs and/or a jacobian in \p matrix,
128 : * as requested.
129 : * Note that in some cases only
130 : * \link current_local_solution \endlink is used during assembly,
131 : * and, therefore, if \link solution \endlink has been altered
132 : * without \link update() \endlink being called, then the
133 : * user must call \link update() \endlink before calling
134 : * this function.
135 : */
136 : virtual void assembly (bool get_residual,
137 : bool get_jacobian,
138 : bool apply_heterogeneous_constraints = false,
139 : bool apply_no_constraints = false) override = 0;
140 :
141 : /**
142 : * Invokes the solver associated with the system. For steady state
143 : * solvers, this will find a root x where F(x) = 0. For transient
144 : * solvers, this will integrate dx/dt = F(x).
145 : */
146 : virtual void solve () override;
147 :
148 : /**
149 : * This function sets the _is_adjoint boolean member of TimeSolver to
150 : * true and then calls the adjoint_solve in implicit system
151 : */
152 : virtual std::pair<unsigned int, Real>
153 : adjoint_solve (const QoISet & qoi_indices = QoISet()) override;
154 :
155 : /**
156 : * We don't allow systems to be attached to each other
157 : */
158 0 : virtual std::unique_ptr<DifferentiablePhysics> clone_physics() override
159 : {
160 0 : libmesh_not_implemented();
161 : // dummy to avoid compiler warnings, not a real implementation
162 : return std::unique_ptr<DifferentiablePhysics>(nullptr);
163 : }
164 :
165 : /**
166 : * We don't allow systems to be attached to each other
167 : */
168 0 : virtual std::unique_ptr<DifferentiableQoI> clone() override
169 : {
170 0 : libmesh_not_implemented();
171 : // dummy to avoid compiler warnings, not a real implementation
172 : return std::unique_ptr<DifferentiableQoI>(nullptr);
173 : }
174 :
175 : /**
176 : * \returns A const reference to a DifferentiablePhysics object.
177 : *
178 : * \note If no external Physics objects have been attached, the default is
179 : * \p this.
180 : */
181 7569832 : const DifferentiablePhysics * get_physics() const
182 84593108 : { if (this->_diff_physics.empty())
183 84449492 : return this;
184 130560 : return this->_diff_physics.top().get(); }
185 :
186 : /**
187 : * \returns A reference to a DifferentiablePhysics object.
188 : *
189 : * \note If no external Physics object is attached, the default is
190 : * \p this.
191 : */
192 9444366 : DifferentiablePhysics * get_physics()
193 104827488 : { if (this->_diff_physics.empty())
194 104326176 : return this;
195 461184 : return this->_diff_physics.top().get(); }
196 :
197 : /**
198 : * Attach external Physics object.
199 : */
200 : void attach_physics( DifferentiablePhysics * physics_in )
201 : { this->_diff_physics.push(physics_in->clone_physics());
202 : this->_diff_physics.top()->init_physics(*this);}
203 :
204 : #ifdef LIBMESH_ENABLE_DEPRECATED
205 : /**
206 : * Swap current physics object with external object. This is
207 : * \deprecated for being too unsafe for easy use.
208 : */
209 : void swap_physics ( DifferentiablePhysics * & swap_physics );
210 : #endif // LIBMESH_ENABLE_DEPRECATED
211 :
212 : /**
213 : * Push a clone of a new physics object onto our stack, overriding
214 : * the current physics until the new physics is popped off again (or
215 : * until something else is pushed on top of it).
216 : */
217 : void push_physics ( DifferentiablePhysics & new_physics );
218 :
219 : /**
220 : * Pop a physics object off of our stack.
221 : */
222 : void pop_physics ();
223 :
224 : /**
225 : * \returns A const reference to a DifferentiableQoI object.
226 : *
227 : * \note If no external QoI object is attached, the default is \p this.
228 : */
229 : const DifferentiableQoI * get_qoi() const
230 : { if (this->_diff_qoi.empty())
231 : return this;
232 : return this->_diff_qoi.top().get(); }
233 :
234 : /**
235 : * \returns A reference to a DifferentiableQoI object.
236 : *
237 : * \note If no external QoI object is attached, the default is this.
238 : */
239 3500 : DifferentiableQoI * get_qoi()
240 117781 : { if (this->_diff_qoi.empty())
241 111981 : return this;
242 5636 : return this->_diff_qoi.top().get(); }
243 :
244 : /**
245 : * Attach external QoI object.
246 : */
247 : void attach_qoi( DifferentiableQoI * qoi_in );
248 :
249 : /**
250 : * A pointer to the solver object we're going to use.
251 : * This must be instantiated by the user before solving!
252 : */
253 : std::unique_ptr<TimeSolver> time_solver;
254 :
255 : /**
256 : * Sets the time_solver
257 : * FIXME: This code is a little dangerous as it transfers ownership
258 : * from the TimeSolver creator to this class. The user must no longer
259 : * access his original TimeSolver object after calling this function.
260 : */
261 : void set_time_solver(std::unique_ptr<TimeSolver> _time_solver)
262 : {
263 : time_solver.reset(_time_solver.release());
264 : }
265 :
266 : /**
267 : * \returns A pointer to the time solver attached to the calling system
268 : */
269 : TimeSolver & get_time_solver();
270 :
271 : /**
272 : * Non-const version of the above
273 : */
274 : const TimeSolver & get_time_solver() const;
275 :
276 : /**
277 : * For time-dependent problems, this is the amount delta t to advance the
278 : * solution in time.
279 : */
280 : Real deltat;
281 :
282 : /**
283 : * Builds a DiffContext object with enough information to do
284 : * evaluations on each element.
285 : *
286 : * For most problems, the default "Let FEMSystem build an * FEMContext"
287 : * reimplementation is correct; users who subclass FEMContext will need to
288 : * also reimplement this method to build it.
289 : */
290 : virtual std::unique_ptr<DiffContext> build_context();
291 :
292 : /**
293 : * Executes a postprocessing loop over all elements, and if
294 : * \p postprocess_sides is true over all sides.
295 : */
296 0 : virtual void postprocess () {}
297 :
298 : /**
299 : * Does any work that needs to be done on \p elem in a postprocessing loop.
300 : */
301 20482 : virtual void element_postprocess (DiffContext &) {}
302 :
303 : /**
304 : * Does any work that needs to be done on \p side of \p elem in a
305 : * postprocessing loop.
306 : */
307 1232 : virtual void side_postprocess (DiffContext &) {}
308 :
309 : /**
310 : * For a given second order (in time) variable var, this method will return
311 : * the index to the corresponding "dot" variable. For FirstOrderUnsteadySolver
312 : * classes, the "dot" variable would automatically be added and the returned
313 : * index will correspond to that variable. For SecondOrderUnsteadySolver classes,
314 : * this method will return var as there this is no "dot" variable per se, but
315 : * having this function allows one to use the interface to treat both
316 : * FirstOrderUnsteadySolver and SecondOrderUnsteadySolver simultaneously.
317 : */
318 : unsigned int get_second_order_dot_var( unsigned int var ) const;
319 :
320 : /**
321 : * Check for any first order vars that are also belong to FEFamily::SCALAR
322 : */
323 : bool have_first_order_scalar_vars() const;
324 :
325 : /**
326 : * Check for any second order vars that are also belong to FEFamily::SCALAR
327 : */
328 : bool have_second_order_scalar_vars() const;
329 :
330 : /**
331 : * If \p postprocess_sides is true (it is false by default), the
332 : * postprocessing loop will loop over all sides as well as all elements.
333 : */
334 : bool postprocess_sides;
335 :
336 : /**
337 : * Set print_residual_norms to true to print |U| whenever it is
338 : * used in an assembly() call
339 : */
340 : bool print_solution_norms;
341 :
342 : /**
343 : * Set print_solutions to true to print U whenever it is used in an
344 : * assembly() call
345 : */
346 : bool print_solutions;
347 :
348 : /**
349 : * Set print_residual_norms to true to print |F| whenever it is assembled.
350 : */
351 : bool print_residual_norms;
352 :
353 : /**
354 : * Set print_residuals to true to print F whenever it is assembled.
355 : */
356 : bool print_residuals;
357 :
358 : /**
359 : * Set print_jacobian_norms to true to print |J| whenever it is assembled.
360 : */
361 : bool print_jacobian_norms;
362 :
363 : /**
364 : * Set print_jacobians to true to print J whenever it is assembled.
365 : */
366 : bool print_jacobians;
367 :
368 : /**
369 : * Set print_element_solutions to true to print each U_elem input.
370 : */
371 : bool print_element_solutions;
372 :
373 : /**
374 : * Set print_element_residuals to true to print each R_elem contribution.
375 : */
376 : bool print_element_residuals;
377 :
378 : /**
379 : * Set print_element_jacobians to true to print each J_elem contribution.
380 : */
381 : bool print_element_jacobians;
382 :
383 : /**
384 : * set_constrain_in_solver to false to apply constraints only via
385 : * residual terms in the systems to be solved.
386 : */
387 : virtual void set_constrain_in_solver(bool enable);
388 :
389 19155839 : virtual bool get_constrain_in_solver()
390 : {
391 19155839 : return _constrain_in_solver;
392 : }
393 :
394 : protected:
395 :
396 : /**
397 : * Initializes the member data fields associated with
398 : * the system, so that, e.g., \p assemble() may be used.
399 : *
400 : * If the TimeSolver is a FirstOrderUnsteadySolver,
401 : * we check for second order in time variables and then add them
402 : * to the System as dot_<varname>. Then, during assembly, the
403 : * TimeSolver will populate the elem_accel vectors with the
404 : * dot_<varname> values so the user's element assembly function
405 : * can still treat the variable as a second order in time variable.
406 : */
407 : virtual void init_data () override;
408 :
409 : /**
410 : * Helper function to add "velocity" variables that are cousins to
411 : * second order-in-time variables in the DifferentiableSystem. This
412 : * function is only called if the TimeSolver is a FirstOrderUnsteadySolver.
413 : */
414 : void add_second_order_dot_vars();
415 :
416 : #ifdef LIBMESH_ENABLE_DIRICHLET
417 : /**
418 : * Helper function to and Dirichlet boundary conditions to "dot" variable
419 : * cousins of second order variables in the system. The function takes the
420 : * second order variable index, it's corresponding "dot" variable index and
421 : * then searches for DirichletBoundary objects for var_idx and then adds a
422 : * DirichletBoundary object for dot_var_idx using the same boundary ids and
423 : * functors for the var_idx DirichletBoundary.
424 : */
425 : void add_dot_var_dirichlet_bcs( unsigned int var_idx, unsigned int dot_var_idx);
426 : #endif
427 :
428 : /**
429 : * _constrain_in_solver defaults to true; if false then we apply
430 : * constraints only via residual terms in the systems to be solved.
431 : */
432 : bool _constrain_in_solver;
433 :
434 : private:
435 : /**
436 : * Stack of pointers to objects to use for physics assembly
437 : * evaluations. Physics assembly defaults to \p this for backwards
438 : * compatibility if the stack is empty; for the most flexibility
439 : * users should create separate physics objects.
440 : */
441 : std::stack<std::unique_ptr<DifferentiablePhysics>,
442 : std::vector<std::unique_ptr<DifferentiablePhysics>>> _diff_physics;
443 :
444 : /**
445 : * Pointer to object to use for quantity of interest assembly
446 : * evaluations. Defaults to \p this for backwards compatibility; in
447 : * the future users should create separate physics objects.
448 : */
449 : std::stack<std::unique_ptr<DifferentiableQoI>,
450 : std::vector<std::unique_ptr<DifferentiableQoI>>> _diff_qoi;
451 : };
452 :
453 : // --------------------------------------------------------------
454 : // DifferentiableSystem inline methods
455 : inline
456 15334 : TimeSolver & DifferentiableSystem::get_time_solver()
457 : {
458 15334 : libmesh_assert(time_solver.get());
459 15334 : libmesh_assert_equal_to (&(time_solver->system()), this);
460 15334 : return *time_solver;
461 : }
462 :
463 : inline
464 14932578 : const TimeSolver & DifferentiableSystem::get_time_solver() const
465 : {
466 14932578 : libmesh_assert(time_solver.get());
467 14932578 : libmesh_assert_equal_to (&(time_solver->system()), this);
468 14932578 : return *time_solver;
469 : }
470 :
471 : } // namespace libMesh
472 :
473 :
474 : #endif // LIBMESH_DIFF_SYSTEM_H
|
