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 558 : 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 8116134 : const DifferentiablePhysics * get_physics() const
182 91181278 : { if (this->_diff_physics.empty())
183 91037662 : 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 9910958 : DifferentiablePhysics * get_physics()
193 111164762 : { if (this->_diff_physics.empty())
194 110663450 : 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 : /**
205 : * Push a clone of a new physics object onto our stack, overriding
206 : * the current physics until the new physics is popped off again (or
207 : * until something else is pushed on top of it).
208 : */
209 : void push_physics ( DifferentiablePhysics & new_physics );
210 :
211 : /**
212 : * Pop a physics object off of our stack.
213 : */
214 : void pop_physics ();
215 :
216 : /**
217 : * \returns A const reference to a DifferentiableQoI object.
218 : *
219 : * \note If no external QoI object is attached, the default is \p this.
220 : */
221 : const DifferentiableQoI * get_qoi() const
222 : { if (this->_diff_qoi.empty())
223 : return this;
224 : return this->_diff_qoi.top().get(); }
225 :
226 : /**
227 : * \returns A reference to a DifferentiableQoI object.
228 : *
229 : * \note If no external QoI object is attached, the default is this.
230 : */
231 3500 : DifferentiableQoI * get_qoi()
232 117781 : { if (this->_diff_qoi.empty())
233 111981 : return this;
234 5636 : return this->_diff_qoi.top().get(); }
235 :
236 : /**
237 : * Attach external QoI object.
238 : */
239 : void attach_qoi( DifferentiableQoI * qoi_in );
240 :
241 : /**
242 : * A pointer to the solver object we're going to use.
243 : * This must be instantiated by the user before solving!
244 : */
245 : std::unique_ptr<TimeSolver> time_solver;
246 :
247 : /**
248 : * Sets the time_solver
249 : * FIXME: This code is a little dangerous as it transfers ownership
250 : * from the TimeSolver creator to this class. The user must no longer
251 : * access his original TimeSolver object after calling this function.
252 : */
253 : void set_time_solver(std::unique_ptr<TimeSolver> _time_solver)
254 : {
255 : time_solver.reset(_time_solver.release());
256 : }
257 :
258 : /**
259 : * \returns A pointer to the time solver attached to the calling system
260 : */
261 : TimeSolver & get_time_solver();
262 :
263 : /**
264 : * Non-const version of the above
265 : */
266 : const TimeSolver & get_time_solver() const;
267 :
268 : /**
269 : * For time-dependent problems, this is the amount delta t to advance the
270 : * solution in time.
271 : */
272 : Real deltat;
273 :
274 : /**
275 : * Builds a DiffContext object with enough information to do
276 : * evaluations on each element.
277 : *
278 : * For most problems, the default "Let FEMSystem build an * FEMContext"
279 : * reimplementation is correct; users who subclass FEMContext will need to
280 : * also reimplement this method to build it.
281 : */
282 : virtual std::unique_ptr<DiffContext> build_context();
283 :
284 : /**
285 : * Executes a postprocessing loop over all elements, and if
286 : * \p postprocess_sides is true over all sides.
287 : */
288 0 : virtual void postprocess () {}
289 :
290 : /**
291 : * Does any work that needs to be done on \p elem in a postprocessing loop.
292 : */
293 20482 : virtual void element_postprocess (DiffContext &) {}
294 :
295 : /**
296 : * Does any work that needs to be done on \p side of \p elem in a
297 : * postprocessing loop.
298 : */
299 1232 : virtual void side_postprocess (DiffContext &) {}
300 :
301 : /**
302 : * For a given second order (in time) variable var, this method will return
303 : * the index to the corresponding "dot" variable. For FirstOrderUnsteadySolver
304 : * classes, the "dot" variable would automatically be added and the returned
305 : * index will correspond to that variable. For SecondOrderUnsteadySolver classes,
306 : * this method will return var as there this is no "dot" variable per se, but
307 : * having this function allows one to use the interface to treat both
308 : * FirstOrderUnsteadySolver and SecondOrderUnsteadySolver simultaneously.
309 : */
310 : unsigned int get_second_order_dot_var( unsigned int var ) const;
311 :
312 : /**
313 : * Check for any first order vars that are also belong to FEFamily::SCALAR
314 : */
315 : bool have_first_order_scalar_vars() const;
316 :
317 : /**
318 : * Check for any second order vars that are also belong to FEFamily::SCALAR
319 : */
320 : bool have_second_order_scalar_vars() const;
321 :
322 : /**
323 : * If \p postprocess_sides is true (it is false by default), the
324 : * postprocessing loop will loop over all sides as well as all elements.
325 : */
326 : bool postprocess_sides;
327 :
328 : /**
329 : * Set print_residual_norms to true to print |U| whenever it is
330 : * used in an assembly() call
331 : */
332 : bool print_solution_norms;
333 :
334 : /**
335 : * Set print_solutions to true to print U whenever it is used in an
336 : * assembly() call
337 : */
338 : bool print_solutions;
339 :
340 : /**
341 : * Set print_residual_norms to true to print |F| whenever it is assembled.
342 : */
343 : bool print_residual_norms;
344 :
345 : /**
346 : * Set print_residuals to true to print F whenever it is assembled.
347 : */
348 : bool print_residuals;
349 :
350 : /**
351 : * Set print_jacobian_norms to true to print |J| whenever it is assembled.
352 : */
353 : bool print_jacobian_norms;
354 :
355 : /**
356 : * Set print_jacobians to true to print J whenever it is assembled.
357 : */
358 : bool print_jacobians;
359 :
360 : /**
361 : * Set print_element_solutions to true to print each U_elem input.
362 : */
363 : bool print_element_solutions;
364 :
365 : /**
366 : * Set print_element_residuals to true to print each R_elem contribution.
367 : */
368 : bool print_element_residuals;
369 :
370 : /**
371 : * Set print_element_jacobians to true to print each J_elem contribution.
372 : */
373 : bool print_element_jacobians;
374 :
375 : /**
376 : * set_constrain_in_solver to false to apply constraints only via
377 : * residual terms in the systems to be solved.
378 : */
379 : virtual void set_constrain_in_solver(bool enable);
380 :
381 20335877 : virtual bool get_constrain_in_solver()
382 : {
383 20335877 : return _constrain_in_solver;
384 : }
385 :
386 : protected:
387 :
388 : /**
389 : * Initializes the member data fields associated with
390 : * the system, so that, e.g., \p assemble() may be used.
391 : *
392 : * If the TimeSolver is a FirstOrderUnsteadySolver,
393 : * we check for second order in time variables and then add them
394 : * to the System as dot_<varname>. Then, during assembly, the
395 : * TimeSolver will populate the elem_accel vectors with the
396 : * dot_<varname> values so the user's element assembly function
397 : * can still treat the variable as a second order in time variable.
398 : */
399 : virtual void init_data () override;
400 :
401 : /**
402 : * Helper function to add "velocity" variables that are cousins to
403 : * second order-in-time variables in the DifferentiableSystem. This
404 : * function is only called if the TimeSolver is a FirstOrderUnsteadySolver.
405 : */
406 : void add_second_order_dot_vars();
407 :
408 : #ifdef LIBMESH_ENABLE_DIRICHLET
409 : /**
410 : * Helper function to and Dirichlet boundary conditions to "dot" variable
411 : * cousins of second order variables in the system. The function takes the
412 : * second order variable index, it's corresponding "dot" variable index and
413 : * then searches for DirichletBoundary objects for var_idx and then adds a
414 : * DirichletBoundary object for dot_var_idx using the same boundary ids and
415 : * functors for the var_idx DirichletBoundary.
416 : */
417 : void add_dot_var_dirichlet_bcs( unsigned int var_idx, unsigned int dot_var_idx);
418 : #endif
419 :
420 : /**
421 : * _constrain_in_solver defaults to true; if false then we apply
422 : * constraints only via residual terms in the systems to be solved.
423 : */
424 : bool _constrain_in_solver;
425 :
426 : private:
427 : /**
428 : * Stack of pointers to objects to use for physics assembly
429 : * evaluations. Physics assembly defaults to \p this for backwards
430 : * compatibility if the stack is empty; for the most flexibility
431 : * users should create separate physics objects.
432 : */
433 : std::stack<std::unique_ptr<DifferentiablePhysics>,
434 : std::vector<std::unique_ptr<DifferentiablePhysics>>> _diff_physics;
435 :
436 : /**
437 : * Pointer to object to use for quantity of interest assembly
438 : * evaluations. Defaults to \p this for backwards compatibility; in
439 : * the future users should create separate physics objects.
440 : */
441 : std::stack<std::unique_ptr<DifferentiableQoI>,
442 : std::vector<std::unique_ptr<DifferentiableQoI>>> _diff_qoi;
443 : };
444 :
445 : // --------------------------------------------------------------
446 : // DifferentiableSystem inline methods
447 : inline
448 27068 : TimeSolver & DifferentiableSystem::get_time_solver()
449 : {
450 27068 : libmesh_assert(time_solver.get());
451 27068 : libmesh_assert_equal_to (&(time_solver->system()), this);
452 27068 : return *time_solver;
453 : }
454 :
455 : inline
456 15878557 : const TimeSolver & DifferentiableSystem::get_time_solver() const
457 : {
458 15878557 : libmesh_assert(time_solver.get());
459 15878557 : libmesh_assert_equal_to (&(time_solver->system()), this);
460 15878557 : return *time_solver;
461 : }
462 :
463 : } // namespace libMesh
464 :
465 :
466 : #endif // LIBMESH_DIFF_SYSTEM_H
|
