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_PHYSICS_H
21 : #define LIBMESH_DIFF_PHYSICS_H
22 :
23 : // Local Includes
24 : #include "libmesh/libmesh.h"
25 :
26 : // C++ includes
27 : #include <map>
28 : #include <memory>
29 : #include <set>
30 : #include <vector>
31 :
32 : namespace libMesh
33 : {
34 :
35 : // Forward Declarations
36 : class System;
37 : class DiffContext;
38 :
39 : /**
40 : * This class provides a specific system class. It aims
41 : * to generalize any system, linear or nonlinear, which
42 : * provides both a residual and a Jacobian.
43 : * For first order (in time) systems, the (nonlinear)
44 : * residual computed at each time step is
45 : * \f[ F(u) + G(u) - M(u,\dot{u})\dot{u} = 0 \f]
46 : * for unsteady TimeSolver and \f$ F(u) = 0\f$ for steady
47 : * TimeSolver. \f$F(u)\f$ is computed by element/side_time_derivative,
48 : * \f$ G(u) \f$ is computed using element/side_constraint, and
49 : * \f$ M(u,\dot{u})\dot{u} \f$ is computed using the mass_residual
50 : * methods.
51 : *
52 : * For second order (in time) systems, the (nonlinear)
53 : * residual computed at each time step is
54 : * \f[ M(u,\ddot{u})\ddot{u} + C(u,\dot{u})\dot{u} + F(u) + G(u) = 0 \f]
55 : * for unsteady TimeSolver and \f$ F(u) = 0\f$ for steady
56 : * TimeSolver. \f$F(u)\f$ is computed by element/side_time_derivative,
57 : * \f$G(u)\f$ is computed using element/side_constraint,
58 : * \f$C(u,\dot{u})\dot{u}\f$ is computed using the damping_residual methods
59 : * and \f$ -M(u,\ddot{u})\ddot{u}\f$ is computed using the mass_residual
60 : * methods. This is the sign convention used by the default implementation;
61 : * if the method is overridden, the user can choose any self-consistent sign
62 : * convention they wish.
63 : *
64 : * FEMContext provides methods for querying values of the solution \f${u}\f$,
65 : * its "rate" \f$\dot{u}\f$ and its "acceleration" \f$\ddot{u}\f$. Furthermore,
66 : * derivatives of each of these w.r.t the nonlinear iteration unknown (e.g. in
67 : * EulerSolver, the solution at the next time step \f$ u_{n+1} \f$) are provided
68 : * through DiffContext::get_elem_solution_derivative(),
69 : * DiffContext::get_elem_solution_rate_derivative(), and
70 : * DiffContext::get_elem_solution_accel_derivative(). The should be incorporated
71 : * into the Jacobian evaluations, if the Jacobian is being provided.
72 : *
73 : * \author Roy H. Stogner
74 : * \date 2006
75 : */
76 4364 : class DifferentiablePhysics
77 : {
78 : public:
79 :
80 : /**
81 : * Constructor. Optionally initializes required
82 : * data structures.
83 : */
84 828 : DifferentiablePhysics () :
85 804 : compute_internal_sides (false),
86 804 : _mesh_sys (nullptr),
87 804 : _mesh_x_var (libMesh::invalid_uint),
88 804 : _mesh_y_var (libMesh::invalid_uint),
89 828 : _mesh_z_var (libMesh::invalid_uint)
90 804 : {}
91 :
92 : /**
93 : * Destructor.
94 : */
95 : virtual ~DifferentiablePhysics ();
96 :
97 : /**
98 : * Copy of this object. User should override to copy any needed state.
99 : */
100 : virtual std::unique_ptr<DifferentiablePhysics> clone_physics() = 0;
101 :
102 : /**
103 : * Clear any data structures associated with the physics.
104 : */
105 : virtual void clear_physics ();
106 :
107 : /**
108 : * Initialize any data structures associated with the physics.
109 : */
110 : virtual void init_physics (const System & sys);
111 :
112 : /**
113 : * Adds the time derivative contribution on \p elem to elem_residual.
114 : * If this method receives request_jacobian = true, then it
115 : * should compute elem_jacobian and return true if possible. If
116 : * elem_jacobian has not been computed then the method should
117 : * return false.
118 : *
119 : * Users need to reimplement this for their particular PDE.
120 : *
121 : * To implement the physics model du/dt = F(u), the user should
122 : * examine u = elem_solution and add (F(u), phi_i) to elem_residual
123 : * in elem_time_derivative().
124 : */
125 0 : virtual bool element_time_derivative (bool request_jacobian,
126 : DiffContext &)
127 : {
128 0 : return request_jacobian;
129 : }
130 :
131 : /**
132 : * Adds the constraint contribution on \p elem to elem_residual.
133 : * If this method receives request_jacobian = true, then it
134 : * should compute elem_jacobian and return true if possible. If
135 : * elem_jacobian has not been computed then the method should
136 : * return false.
137 : *
138 : * Users may need to reimplement this for their particular PDE.
139 : *
140 : * To implement the constraint 0 = G(u), the user should
141 : * examine u = elem_solution and add (G(u), phi_i) to elem_residual
142 : * in elem_constraint().
143 : */
144 0 : virtual bool element_constraint (bool request_jacobian,
145 : DiffContext &)
146 : {
147 0 : return request_jacobian;
148 : }
149 :
150 : /**
151 : * \p compute_internal_sides is false by default, indicating that
152 : * side_* computations will only be done on boundary sides. If
153 : * compute_internal_sides is true, computations will be done
154 : * on sides between elements as well.
155 : */
156 : bool compute_internal_sides;
157 :
158 : /**
159 : * Adds the time derivative contribution on \p side of \p elem to
160 : * elem_residual.
161 : * If this method receives request_jacobian = true, then it
162 : * should compute elem_jacobian and return true if possible. If
163 : * elem_jacobian has not been computed then the method should
164 : * return false.
165 : *
166 : * Users may need to reimplement this for their particular PDE
167 : * depending on the boundary conditions.
168 : *
169 : * To implement a weak form of the source term du/dt = F(u) on
170 : * sides, such as might arise in a flux boundary condition, the user
171 : * should examine u = elem_solution and add (F(u), phi_i) boundary
172 : * integral contributions to elem_residual in side_constraint().
173 : */
174 0 : virtual bool side_time_derivative (bool request_jacobian,
175 : DiffContext &)
176 : {
177 0 : return request_jacobian;
178 : }
179 :
180 : /**
181 : * Adds the constraint contribution on \p side of \p elem to
182 : * elem_residual.
183 : * If this method receives request_jacobian = true, then it
184 : * should compute elem_jacobian and return true if possible. If
185 : * elem_jacobian has not been computed then the method should
186 : * return false.
187 : *
188 : * Users may need to reimplement this for their particular PDE
189 : * depending on the boundary conditions.
190 : *
191 : * To implement a weak form of the constraint 0 = G(u), the user
192 : * should examine u = elem_solution and add (G(u), phi_i) boundary
193 : * integral contributions to elem_residual in side_constraint().
194 : */
195 0 : virtual bool side_constraint (bool request_jacobian,
196 : DiffContext &)
197 : {
198 0 : return request_jacobian;
199 : }
200 :
201 : /**
202 : * Adds any nonlocal time derivative contributions (e.g. some
203 : * components of time derivatives in scalar variable equations) to
204 : * elem_residual
205 : *
206 : * If this method receives request_jacobian = true, then it
207 : * should also modify elem_jacobian and return true if possible. If
208 : * the Jacobian changes have not been computed then the method
209 : * should return false.
210 : *
211 : * Users may need to reimplement this for PDEs on systems to which
212 : * SCALAR variables have been added.
213 : */
214 0 : virtual bool nonlocal_time_derivative (bool request_jacobian,
215 : DiffContext &)
216 : {
217 0 : return request_jacobian;
218 : }
219 :
220 : /**
221 : * Adds any nonlocal constraint contributions (e.g. some
222 : * components of constraints in scalar variable equations) to
223 : * elem_residual
224 : *
225 : * If this method receives request_jacobian = true, then it
226 : * should also modify elem_jacobian and return true if possible. If
227 : * the Jacobian changes have not been computed then the method
228 : * should return false.
229 : *
230 : * Users may need to reimplement this for PDEs on systems to which
231 : * SCALAR variables with non-transient equations have been added.
232 : */
233 0 : virtual bool nonlocal_constraint (bool request_jacobian,
234 : DiffContext &)
235 : {
236 0 : return request_jacobian;
237 : }
238 :
239 : /**
240 : * Tells the DiffSystem that variable var is evolving with
241 : * respect to time. In general, the user's init() function
242 : * should call time_evolving() with order 1 for any variables which
243 : * behave like du/dt = F(u), with order 2 for any variables that
244 : * behave like d^2u/dt^2 = F(u), and should not call time_evolving()
245 : * for any variables which behave like 0 = G(u).
246 : *
247 : * Most derived systems will not have to reimplement this function; however
248 : * any system which reimplements mass_residual() may have to reimplement
249 : * time_evolving() to prepare data structures.
250 : */
251 : virtual void time_evolving (unsigned int var, unsigned int order);
252 :
253 : /**
254 : * \returns \p true iff variable \p var is evolving with
255 : * respect to time. In general, the user's init() function
256 : * should have set time_evolving() for any variables which
257 : * behave like du/dt = F(u), and should not call time_evolving()
258 : * for any variables which behave like 0 = G(u).
259 : */
260 1507262 : bool is_time_evolving (unsigned int var) const
261 : {
262 1507262 : libmesh_assert_less(var,_time_evolving.size());
263 1507262 : libmesh_assert( _time_evolving[var] == 0 ||
264 : _time_evolving[var] == 1 ||
265 : _time_evolving[var] == 2 );
266 18191156 : return _time_evolving[var];
267 : }
268 :
269 : /**
270 : * Adds a pseudo-convection contribution on \p elem to
271 : * elem_residual, if the nodes of \p elem are being translated by a
272 : * moving mesh.
273 : *
274 : * The library provides a basic implementation in
275 : * FEMPhysics::eulerian_residual()
276 : */
277 0 : virtual bool eulerian_residual (bool request_jacobian,
278 : DiffContext &)
279 : {
280 0 : return request_jacobian;
281 : }
282 :
283 : /**
284 : * Subtracts a mass vector contribution on \p elem from
285 : * elem_residual. For first-order-in-time problems, this
286 : * is the \f$ M(u,\dot{u})\dot{u} \f$ term. For
287 : * second-order-in-time problems, this is the
288 : * \f$ M(u,\ddot{u})\ddot{u} \f$ term. This method is only
289 : * called for UnsteadySolver-based TimeSolvers.
290 : *
291 : * If this method receives request_jacobian = true, then it
292 : * should compute elem_jacobian and return true if possible. If
293 : * elem_jacobian has not been computed then the method should
294 : * return false.
295 : *
296 : * Many first-order-in-time problems can use the reimplementation in
297 : * FEMPhysics::mass_residual which subtracts (du/dt,v) for each
298 : * transient variable u; users with more complicated transient
299 : * problems or second-order-in-time problems will need to reimplement
300 : * this themselves.
301 : */
302 0 : virtual bool mass_residual (bool request_jacobian,
303 : DiffContext &)
304 : {
305 0 : return request_jacobian;
306 : }
307 :
308 : /**
309 : * Subtracts a mass vector contribution on \p side of \p elem from
310 : * elem_residual.
311 : * If this method receives request_jacobian = true, then it
312 : * should compute elem_jacobian and return true if possible. If
313 : * elem_jacobian has not been computed then the method should
314 : * return false.
315 : *
316 : * For most problems, the default implementation of "do nothing"
317 : * is correct; users with boundary conditions including time
318 : * derivatives may need to reimplement this themselves.
319 : */
320 0 : virtual bool side_mass_residual (bool request_jacobian,
321 : DiffContext &)
322 : {
323 0 : return request_jacobian;
324 : }
325 :
326 : /**
327 : * Subtracts any nonlocal mass vector contributions (e.g. any time
328 : * derivative coefficients in scalar variable equations) from
329 : * elem_residual
330 : *
331 : * If this method receives request_jacobian = true, then it
332 : * should also modify elem_jacobian and return true if possible. If
333 : * the Jacobian changes have not been computed then the method
334 : * should return false.
335 : *
336 : * Many problems can use the reimplementation in
337 : * FEMPhysics::mass_residual which subtracts (du/dt,v) for each
338 : * transient scalar variable u; users with more complicated
339 : * transient scalar variable equations will need to reimplement this
340 : * themselves.
341 : */
342 : virtual bool nonlocal_mass_residual (bool request_jacobian,
343 : DiffContext & c);
344 :
345 : /**
346 : * Subtracts a damping vector contribution on \p elem from
347 : * elem_residual. This method is not used in first-order-in-time
348 : * problems. For second-order-in-time problems, this is the
349 : * \f$ C(u,\ddot{u})\ddot{u} \f$ term. This method is only
350 : * called for UnsteadySolver-based TimeSolvers.
351 : *
352 : * If this method receives request_jacobian = true, then it
353 : * should compute elem_jacobian and return true if possible. If
354 : * elem_jacobian has not been computed then the method should
355 : * return false.
356 : *
357 : * If the problem has no damping, the default "do-nothing" is correct.
358 : * Otherwise, this must be reimplemented.
359 : */
360 0 : virtual bool damping_residual (bool request_jacobian,
361 : DiffContext &)
362 : {
363 0 : return request_jacobian;
364 : }
365 :
366 : /**
367 : * Subtracts a damping vector contribution on \p side of \p elem from
368 : * elem_residual.
369 : * If this method receives request_jacobian = true, then it
370 : * should compute elem_jacobian and return true if possible. If
371 : * elem_jacobian has not been computed then the method should
372 : * return false.
373 : *
374 : * For most problems, the default implementation of "do nothing"
375 : * is correct; users with boundary conditions including first time
376 : * derivatives may need to reimplement this themselves.
377 : */
378 0 : virtual bool side_damping_residual (bool request_jacobian,
379 : DiffContext &)
380 : {
381 0 : return request_jacobian;
382 : }
383 :
384 : /**
385 : * Subtracts any nonlocal damping vector contributions (e.g. any
386 : * first time derivative coefficients in scalar variable equations) from
387 : * elem_residual
388 : *
389 : * If this method receives request_jacobian = true, then it
390 : * should also modify elem_jacobian and return true if possible. If
391 : * the Jacobian changes have not been computed then the method
392 : * should return false.
393 : */
394 0 : virtual bool nonlocal_damping_residual (bool request_jacobian,
395 : DiffContext &)
396 : {
397 0 : return request_jacobian;
398 : }
399 :
400 : /*
401 : * Prepares the result of a build_context() call for use.
402 : *
403 : * Most FEMSystem-based problems will need to reimplement this in order to
404 : * call FE::get_*() as their particular physics requires.
405 : */
406 0 : virtual void init_context(DiffContext &) {}
407 :
408 : /**
409 : * Tells the DifferentiablePhysics that system \p sys contains the
410 : * isoparametric Lagrangian variables which correspond to the
411 : * coordinates of mesh nodes, in problems where the mesh itself is
412 : * expected to move in time.
413 : *
414 : * The system with mesh coordinate data (which may be \p this system
415 : * itself, for fully coupled moving mesh problems) is currently
416 : * assumed to have new (end of time step) mesh coordinates stored in
417 : * solution, old (beginning of time step) mesh coordinates stored in
418 : * _old_nonlinear_solution, and constant velocity motion during each
419 : * time step.
420 : *
421 : * Activating this function ensures that local (but not neighbor!) element
422 : * geometry is correctly repositioned when evaluating element residuals.
423 : *
424 : * Currently \p sys must be \p *this for a tightly coupled moving
425 : * mesh problem or nullptr to stop mesh movement; loosely coupled
426 : * moving mesh problems are not implemented.
427 : *
428 : * This code is experimental. "Trust but verify, and not in that
429 : * order"
430 : */
431 : virtual void set_mesh_system(System * sys);
432 :
433 : /**
434 : * \returns A const reference to the system with variables corresponding to
435 : * mesh nodal coordinates, or nullptr if the mesh is fixed.
436 : * Useful for ALE calculations.
437 : */
438 : const System * get_mesh_system() const;
439 :
440 : /**
441 : * \returns A reference to the system with variables corresponding to
442 : * mesh nodal coordinates, or nullptr if the mesh is fixed.
443 : */
444 : System * get_mesh_system();
445 :
446 : /**
447 : * Tells the DifferentiablePhysics that variable \p var from the mesh system
448 : * should be used to update the x coordinate of mesh nodes, in problems where
449 : * the mesh itself is expected to move in time.
450 : *
451 : * The system with mesh coordinate data (which may be this system itself, for
452 : * fully coupled moving mesh problems) is currently assumed to have new (end
453 : * of time step) mesh coordinates stored in solution, old (beginning of time
454 : * step) mesh coordinates stored in _old_nonlinear_solution, and constant
455 : * velocity motion during each time step.
456 : *
457 : * Activating this function ensures that local (but not neighbor!) element
458 : * geometry is correctly repositioned when evaluating element residuals.
459 : */
460 : virtual void set_mesh_x_var(unsigned int var);
461 :
462 : /**
463 : * \returns The variable number corresponding to the
464 : * mesh x coordinate. Useful for ALE calculations.
465 : */
466 : unsigned int get_mesh_x_var() const;
467 :
468 : /**
469 : * Tells the DifferentiablePhysics that variable \p var from the mesh system
470 : * should be used to update the y coordinate of mesh nodes.
471 : */
472 : virtual void set_mesh_y_var(unsigned int var);
473 :
474 : /**
475 : * \returns The variable number corresponding to the
476 : * mesh y coordinate. Useful for ALE calculations.
477 : */
478 : unsigned int get_mesh_y_var() const;
479 :
480 : /**
481 : * Tells the DifferentiablePhysics that variable \p var from the mesh system
482 : * should be used to update the z coordinate of mesh nodes.
483 : */
484 : virtual void set_mesh_z_var(unsigned int var);
485 :
486 : /**
487 : * \returns The variable number corresponding to the
488 : * mesh z coordinate. Useful for ALE calculations.
489 : */
490 : unsigned int get_mesh_z_var() const;
491 :
492 : /**
493 : * This method simply combines element_time_derivative() and
494 : * eulerian_residual(), which makes its address useful as a
495 : * pointer-to-member-function when refactoring.
496 : */
497 : bool _eulerian_time_deriv (bool request_jacobian,
498 : DiffContext &);
499 :
500 0 : bool have_first_order_vars() const
501 0 : { return !_first_order_vars.empty(); }
502 :
503 : /**
504 : * \returns The set of first order in time variable indices. May be empty.
505 : */
506 0 : const std::set<unsigned int> & get_first_order_vars() const
507 0 : { return _first_order_vars; }
508 :
509 : bool is_first_order_var( unsigned int var ) const
510 : { return _first_order_vars.find(var) != _first_order_vars.end(); }
511 :
512 :
513 1703600 : bool have_second_order_vars() const
514 18934400 : { return !_second_order_vars.empty(); }
515 :
516 : /**
517 : * \returns The set of second order in time variable indices. May be empty.
518 : */
519 7018544 : const std::set<unsigned int> & get_second_order_vars() const
520 7018544 : { return _second_order_vars; }
521 :
522 471096 : bool is_second_order_var( unsigned int var ) const
523 471096 : { return _second_order_vars.find(var) != _second_order_vars.end(); }
524 :
525 :
526 : protected:
527 :
528 : /**
529 : * System from which to acquire moving mesh information
530 : */
531 : System * _mesh_sys;
532 :
533 : /**
534 : * Variables from which to acquire moving mesh information
535 : */
536 : unsigned int _mesh_x_var, _mesh_y_var, _mesh_z_var;
537 :
538 : /**
539 : * Stores unsigned int to tell us which variables are evolving
540 : * as first order in time (1), second order in time (2), or are
541 : * not time evolving (0).
542 : */
543 : std::vector<unsigned int> _time_evolving;
544 :
545 : /**
546 : * Variable indices for those variables that are first order in time.
547 : */
548 : std::set<unsigned int> _first_order_vars;
549 :
550 : /**
551 : * Variable indices for those variables that are second order in time.
552 : */
553 : std::set<unsigned int> _second_order_vars;
554 :
555 : /**
556 : * If the user adds any second order variables, then we need to also
557 : * cache the map to their corresponding dot variable that will
558 : * be added by this TimeSolver class.
559 : */
560 : std::map<unsigned int,unsigned int> _second_order_dot_vars;
561 :
562 : };
563 :
564 : // ------------------------------------------------------------
565 : // DifferentiablePhysics inline methods
566 :
567 :
568 : inline
569 0 : void DifferentiablePhysics::set_mesh_system(System * sys)
570 : {
571 : // For now we assume that we're doing fully coupled mesh motion
572 : // if (sys && sys != this)
573 : // libmesh_not_implemented();
574 :
575 : // For the foreseeable future we'll assume that we keep these
576 : // Systems in the same EquationSystems
577 : // libmesh_assert_equal_to (&this->get_equation_systems(),
578 : // &sys->get_equation_systems());
579 :
580 : // And for the immediate future this code may not even work
581 : libmesh_experimental();
582 :
583 0 : _mesh_sys = sys;
584 0 : }
585 :
586 :
587 :
588 : inline
589 0 : void DifferentiablePhysics::set_mesh_x_var (unsigned int var)
590 : {
591 0 : _mesh_x_var = var;
592 0 : }
593 :
594 :
595 :
596 : inline
597 0 : void DifferentiablePhysics::set_mesh_y_var (unsigned int var)
598 : {
599 0 : _mesh_y_var = var;
600 0 : }
601 :
602 :
603 :
604 : inline
605 0 : void DifferentiablePhysics::set_mesh_z_var (unsigned int var)
606 : {
607 0 : _mesh_z_var = var;
608 0 : }
609 :
610 :
611 :
612 : inline
613 : const System * DifferentiablePhysics::get_mesh_system() const
614 : {
615 : return _mesh_sys;
616 : }
617 :
618 : inline
619 14404 : System * DifferentiablePhysics::get_mesh_system()
620 : {
621 269920 : return _mesh_sys;
622 : }
623 :
624 : inline
625 14404 : unsigned int DifferentiablePhysics::get_mesh_x_var() const
626 : {
627 269920 : return _mesh_x_var;
628 : }
629 :
630 : inline
631 14404 : unsigned int DifferentiablePhysics::get_mesh_y_var() const
632 : {
633 269920 : return _mesh_y_var;
634 : }
635 :
636 : inline
637 14404 : unsigned int DifferentiablePhysics::get_mesh_z_var() const
638 : {
639 269920 : return _mesh_z_var;
640 : }
641 :
642 :
643 :
644 : } // namespace libMesh
645 :
646 :
647 : #endif // LIBMESH_DIFF_PHYSICS_H
|
