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_CONTEXT_H
21 : #define LIBMESH_DIFF_CONTEXT_H
22 :
23 : // Local Includes
24 : #include "libmesh/dense_matrix.h"
25 : #include "libmesh/dense_submatrix.h"
26 : #include "libmesh/dense_subvector.h"
27 : #include "libmesh/dense_vector.h"
28 : #include "libmesh/id_types.h"
29 :
30 : // C++ includes
31 : #include <cstddef>
32 : #include <map>
33 : #include <vector>
34 : #include <memory>
35 :
36 : namespace libMesh
37 : {
38 :
39 : // Forward declarations
40 : template <typename T> class NumericVector;
41 : class System;
42 :
43 : /**
44 : * This class provides all data required for a physics package
45 : * (e.g. a DifferentiableSystem subclass) to perform local element
46 : * residual and jacobian integrations.
47 : *
48 : * This class is part of the new DifferentiableSystem framework,
49 : * which is still experimental. Users of this framework should
50 : * beware of bugs and future API changes.
51 : *
52 : * \author Roy H. Stogner
53 : * \date 2009
54 : */
55 663032 : class DiffContext
56 : {
57 : public:
58 :
59 : /**
60 : * Constructor. Optionally initializes required
61 : * data structures.
62 : */
63 : explicit
64 : DiffContext (const System &,
65 : bool allocate_local_matrices = true);
66 :
67 : /**
68 : * Destructor.
69 : */
70 : virtual ~DiffContext ();
71 :
72 : /**
73 : * Gives derived classes the opportunity to reinitialize data (FE objects in
74 : * FEMSystem, for example) needed for an interior integration at a new point
75 : * within a timestep
76 : */
77 0 : virtual void elem_reinit(Real /* theta */) {}
78 :
79 : /**
80 : * Gives derived classes the opportunity to reinitialize data needed
81 : * for a side integration at a new point within a timestep
82 : */
83 0 : virtual void elem_side_reinit(Real /* theta */) {}
84 :
85 : /**
86 : * Gives derived classes the opportunity to reinitialize data needed
87 : * for an edge integration at a new point within a timestep
88 : */
89 0 : virtual void elem_edge_reinit(Real /* theta */) {}
90 :
91 : /**
92 : * Gives derived classes the opportunity to reinitialize data needed
93 : * for nonlocal calculations at a new point within a timestep
94 : */
95 0 : virtual void nonlocal_reinit(Real /* theta */) {}
96 :
97 : /**
98 : * Number of variables in solution.
99 : */
100 1861633 : unsigned int n_vars() const
101 3723267 : { return cast_int<unsigned int>(_dof_indices_var.size()); }
102 :
103 : /**
104 : * Accessor for associated system.
105 : */
106 979778 : const System & get_system() const
107 9989832 : { return _system; }
108 :
109 : /**
110 : * Accessor for element solution.
111 : */
112 : const DenseVector<Number> & get_elem_solution() const
113 : { return _elem_solution; }
114 :
115 : /**
116 : * Non-const accessor for element solution.
117 : */
118 14144554 : DenseVector<Number> & get_elem_solution()
119 14144554 : { return _elem_solution; }
120 :
121 : /**
122 : * Accessor for element solution of a particular variable corresponding
123 : * to the variable index argument.
124 : */
125 34066666 : const DenseSubVector<Number> & get_elem_solution( unsigned int var ) const
126 : {
127 34066666 : libmesh_assert_greater(_elem_subsolutions.size(), var);
128 68133360 : return _elem_subsolutions[var];
129 : }
130 :
131 : /**
132 : * Accessor for element solution of a particular variable corresponding
133 : * to the variable index argument.
134 : */
135 6190106 : DenseSubVector<Number> & get_elem_solution( unsigned int var )
136 : {
137 6190106 : libmesh_assert_greater(_elem_subsolutions.size(), var);
138 12672723 : return _elem_subsolutions[var];
139 : }
140 :
141 : /**
142 : * Accessor for element solution rate of change w.r.t. time.
143 : */
144 : const DenseVector<Number> & get_elem_solution_rate() const
145 : { return _elem_solution_rate; }
146 :
147 : /**
148 : * Non-const accessor for element solution rate of change w.r.t.
149 : * time.
150 : */
151 9367340 : DenseVector<Number> & get_elem_solution_rate()
152 56101834 : { return _elem_solution_rate; }
153 :
154 : /**
155 : * Accessor for element solution rate for a particular variable
156 : * corresponding to the variable index argument.
157 : */
158 13085708 : const DenseSubVector<Number> & get_elem_solution_rate( unsigned int var ) const
159 : {
160 13085708 : libmesh_assert_greater(_elem_subsolution_rates.size(), var);
161 26171416 : return _elem_subsolution_rates[var];
162 : }
163 :
164 : /**
165 : * Accessor for element solution rate for a particular variable
166 : * corresponding to the variable index argument.
167 : */
168 3845484 : DenseSubVector<Number> & get_elem_solution_rate( unsigned int var )
169 : {
170 3845484 : libmesh_assert_greater(_elem_subsolution_rates.size(), var);
171 7691006 : return _elem_subsolution_rates[var];
172 : }
173 :
174 : /**
175 : * Accessor for element solution accel of change w.r.t. time.
176 : */
177 : const DenseVector<Number> & get_elem_solution_accel() const
178 : { return _elem_solution_accel; }
179 :
180 : /**
181 : * Non-const accessor for element solution accel of change w.r.t.
182 : * time.
183 : */
184 387088 : DenseVector<Number> & get_elem_solution_accel()
185 1409008 : { return _elem_solution_accel; }
186 :
187 : /**
188 : * Accessor for element solution accel for a particular variable
189 : * corresponding to the variable index argument.
190 : */
191 3369772 : const DenseSubVector<Number> & get_elem_solution_accel( unsigned int var ) const
192 : {
193 3369772 : libmesh_assert_greater(_elem_subsolution_accels.size(), var);
194 6739544 : return _elem_subsolution_accels[var];
195 : }
196 :
197 : /**
198 : * Accessor for element solution accel for a particular variable
199 : * corresponding to the variable index argument.
200 : */
201 369584 : DenseSubVector<Number> & get_elem_solution_accel( unsigned int var )
202 : {
203 369584 : libmesh_assert_greater(_elem_subsolution_accels.size(), var);
204 739168 : return _elem_subsolution_accels[var];
205 : }
206 :
207 : /**
208 : * Accessor for element fixed solution.
209 : */
210 : const DenseVector<Number> & get_elem_fixed_solution() const
211 : { return _elem_fixed_solution; }
212 :
213 : /**
214 : * Non-const accessor for element fixed solution.
215 : */
216 0 : DenseVector<Number> & get_elem_fixed_solution()
217 0 : { return _elem_fixed_solution; }
218 :
219 : /**
220 : * Accessor for element fixed solution of a particular variable corresponding
221 : * to the variable index argument.
222 : */
223 0 : const DenseSubVector<Number> & get_elem_fixed_solution( unsigned int var ) const
224 : {
225 0 : libmesh_assert_greater(_elem_fixed_subsolutions.size(), var);
226 0 : return _elem_fixed_subsolutions[var];
227 : }
228 :
229 : /**
230 : * Accessor for element fixed solution of a particular variable corresponding
231 : * to the variable index argument.
232 : */
233 0 : DenseSubVector<Number> & get_elem_fixed_solution( unsigned int var )
234 : {
235 0 : libmesh_assert_greater(_elem_fixed_subsolutions.size(), var);
236 0 : return _elem_fixed_subsolutions[var];
237 : }
238 :
239 : /**
240 : * Const accessor for element residual.
241 : */
242 : const DenseVector<Number> & get_elem_residual() const
243 : { return _elem_residual; }
244 :
245 : /**
246 : * Non-const accessor for element residual.
247 : */
248 11507464 : DenseVector<Number> & get_elem_residual()
249 79829471 : { return _elem_residual; }
250 :
251 : /**
252 : * Const accessor for element residual of a particular variable corresponding
253 : * to the variable index argument.
254 : */
255 : const DenseSubVector<Number> & get_elem_residual( unsigned int var ) const
256 : {
257 : libmesh_assert_greater(_elem_subresiduals.size(), var);
258 : return _elem_subresiduals[var];
259 : }
260 :
261 : /**
262 : * Non-const accessor for element residual of a particular variable corresponding
263 : * to the variable index argument.
264 : */
265 1210336 : DenseSubVector<Number> & get_elem_residual( unsigned int var )
266 : {
267 1210336 : libmesh_assert_greater(_elem_subresiduals.size(), var);
268 10087141 : return _elem_subresiduals[var];
269 : }
270 :
271 : /**
272 : * Const accessor for element Jacobian.
273 : */
274 : const DenseMatrix<Number> & get_elem_jacobian() const
275 : {
276 : libmesh_assert(_have_local_matrices);
277 : return _elem_jacobian;
278 : }
279 :
280 : /**
281 : * Non-const accessor for element Jacobian.
282 : */
283 7167736 : DenseMatrix<Number> & get_elem_jacobian()
284 : {
285 7167736 : libmesh_assert(_have_local_matrices);
286 50417051 : return _elem_jacobian;
287 : }
288 :
289 : /**
290 : * Const accessor for element Jacobian of particular variables corresponding
291 : * to the variable index arguments.
292 : */
293 : const DenseSubMatrix<Number> & get_elem_jacobian( unsigned int var1, unsigned int var2 ) const
294 : {
295 : libmesh_assert(_have_local_matrices);
296 : libmesh_assert_greater(_elem_subjacobians.size(), var1);
297 : libmesh_assert_greater(_elem_subjacobians[var1].size(), var2);
298 : return _elem_subjacobians[var1][var2];
299 : }
300 :
301 : /**
302 : * Non-const accessor for element Jacobian of particular variables corresponding
303 : * to the variable index arguments. Only available if \p
304 : * _have_local_matrices
305 : */
306 855657 : DenseSubMatrix<Number> & get_elem_jacobian( unsigned int var1, unsigned int var2 )
307 : {
308 855657 : libmesh_assert(_have_local_matrices);
309 855657 : libmesh_assert_greater(_elem_subjacobians.size(), var1);
310 855657 : libmesh_assert_greater(_elem_subjacobians[var1].size(), var2);
311 17955209 : return _elem_subjacobians[var1][var2];
312 : }
313 :
314 : /**
315 : * Const accessor for QoI vector.
316 : */
317 : const std::vector<Number> & get_qois() const
318 : { return _elem_qoi; }
319 :
320 : /**
321 : * Non-const accessor for QoI vector.
322 : */
323 1152 : std::vector<Number> & get_qois()
324 41425 : { return _elem_qoi; }
325 :
326 : /**
327 : * Const accessor for QoI derivatives.
328 : */
329 : const std::vector<DenseVector<Number>> & get_qoi_derivatives() const
330 : { return _elem_qoi_derivative; }
331 :
332 : /**
333 : * Non-const accessor for QoI derivatives.
334 : */
335 7506640 : std::vector<DenseVector<Number>> & get_qoi_derivatives()
336 47012822 : { return _elem_qoi_derivative; }
337 :
338 : /**
339 : * Const accessor for QoI derivative of a particular qoi and variable corresponding
340 : * to the index arguments.
341 : */
342 : const DenseSubVector<Number> & get_qoi_derivatives( std::size_t qoi, unsigned int var ) const
343 : {
344 : libmesh_assert_greater(_elem_qoi_subderivatives.size(), qoi);
345 : libmesh_assert_greater(_elem_qoi_subderivatives[qoi].size(), var);
346 : return _elem_qoi_subderivatives[qoi][var];
347 : }
348 :
349 : /**
350 : * Non-const accessor for QoI derivative of a particular qoi and variable corresponding
351 : * to the index arguments.
352 : */
353 18144 : DenseSubVector<Number> & get_qoi_derivatives( std::size_t qoi, unsigned int var )
354 : {
355 18144 : libmesh_assert_greater(_elem_qoi_subderivatives.size(), qoi);
356 18144 : libmesh_assert_greater(_elem_qoi_subderivatives[qoi].size(), var);
357 6231600 : return _elem_qoi_subderivatives[qoi][var];
358 : }
359 :
360 : /**
361 : * Accessor for element dof indices
362 : */
363 : const std::vector<dof_id_type> & get_dof_indices() const
364 : { return _dof_indices; }
365 :
366 : /**
367 : * Non-const accessor for element dof indices
368 : */
369 45023343 : std::vector<dof_id_type> & get_dof_indices()
370 201515786 : { return _dof_indices; }
371 :
372 : /**
373 : * Accessor for element dof indices of a particular variable corresponding
374 : * to the index argument.
375 : */
376 50719714 : const std::vector<dof_id_type> & get_dof_indices( unsigned int var ) const
377 : {
378 50719714 : libmesh_assert_greater(_dof_indices_var.size(), var);
379 583483476 : return _dof_indices_var[var];
380 : }
381 :
382 : /**
383 : * Accessor for element dof indices of a particular variable corresponding
384 : * to the index argument.
385 : */
386 16391562 : std::vector<dof_id_type> & get_dof_indices( unsigned int var )
387 : {
388 16391562 : libmesh_assert_greater(_dof_indices_var.size(), var);
389 149029780 : return _dof_indices_var[var];
390 : }
391 :
392 : /**
393 : * Total number of dof indices on the element
394 : */
395 : unsigned int n_dof_indices() const
396 : { return cast_int<unsigned int>(_dof_indices.size()); }
397 :
398 : /**
399 : * Total number of dof indices of the particular variable
400 : * corresponding to the index argument
401 : */
402 0 : unsigned int n_dof_indices( unsigned int var ) const
403 : {
404 0 : libmesh_assert_greater(_dof_indices_var.size(), var);
405 0 : return cast_int<unsigned int>(_dof_indices_var[var].size());
406 : }
407 :
408 : /**
409 : * Accessor for the time variable stored in the system class.
410 : */
411 7775952 : Real get_system_time() const
412 43285856 : { return system_time; }
413 :
414 : /**
415 : * Accessor for the time for which the current nonlinear_solution is defined.
416 : */
417 : Real get_time() const
418 : { return time; }
419 :
420 : /**
421 : * Set the time for which the current nonlinear_solution is defined.
422 : */
423 3887976 : void set_time( Real time_in )
424 43285856 : { time = time_in; }
425 :
426 : /**
427 : * The derivative of the current elem_solution w.r.t. the unknown
428 : * solution. Corresponding Jacobian contributions should be
429 : * multiplied by this amount, or may be skipped if
430 : * get_elem_solution_derivative() is 0.
431 : */
432 30 : Real get_elem_solution_derivative() const
433 60581274 : { return elem_solution_derivative; }
434 :
435 : /**
436 : * The derivative of the current elem_solution_rate w.r.t. the
437 : * unknown solution. Corresponding Jacobian contributions should be
438 : * multiplied by this amount, or may be skipped if
439 : * get_elem_solution_rate_derivative() is 0.
440 : */
441 6632080 : Real get_elem_solution_rate_derivative() const
442 67213324 : { return elem_solution_rate_derivative; }
443 :
444 : /**
445 : * The derivative of the current elem_solution_accel w.r.t. the
446 : * unknown solution. Corresponding Jacobian contributions should be
447 : * multiplied by this amount, or may be skipped if
448 : * get_elem_solution_accel_derivative() is 0.
449 : */
450 : Real get_elem_solution_accel_derivative() const
451 : { return elem_solution_accel_derivative; }
452 :
453 : /**
454 : * The derivative of the current fixed_elem_solution w.r.t. the
455 : * unknown solution. Corresponding Jacobian contributions should be
456 : * multiplied by this amount, or may be skipped if
457 : * get_fixed_elem_solution_derivative() is 0.
458 : */
459 : Real get_fixed_solution_derivative() const
460 : { return fixed_solution_derivative; }
461 :
462 : /**
463 : * Accessor for querying whether we need to do a primal
464 : * or adjoint solve
465 : */
466 : bool is_adjoint() const
467 : { return _is_adjoint; }
468 :
469 : /**
470 : * Accessor for setting whether we need to do a primal
471 : * or adjoint solve
472 : */
473 5460 : bool & is_adjoint()
474 5460 : { return _is_adjoint; }
475 :
476 : /**
477 : * For time-dependent problems, this is the time t for which the current
478 : * nonlinear_solution is defined.
479 : * FIXME - this needs to be tweaked mid-timestep by all transient solvers!
480 : */
481 : Real time;
482 :
483 : /**
484 : * This is the time stored in the System class at the time this context
485 : * was created, i.e. the time at the beginning of the current timestep.
486 : * This value gets set in the constructor and unlike DiffContext::time,
487 : * is not tweaked mid-timestep by transient solvers: it remains equal
488 : * to the value it was assigned at construction.
489 : */
490 : const Real system_time;
491 :
492 : /**
493 : * The derivative of elem_solution with respect to the current
494 : * nonlinear solution.
495 : */
496 : Real elem_solution_derivative;
497 :
498 : /**
499 : * The derivative of elem_solution_rate with respect to the current
500 : * nonlinear solution, for use by systems with non default
501 : * mass_residual terms.
502 : */
503 : Real elem_solution_rate_derivative;
504 :
505 : /**
506 : * The derivative of elem_solution_accel with respect to the current
507 : * nonlinear solution, for use by systems with non default
508 : * mass_residual terms.
509 : */
510 : Real elem_solution_accel_derivative;
511 :
512 : /**
513 : * The derivative of elem_fixed_solution with respect to the nonlinear
514 : * solution, for use by systems constructing jacobians with
515 : * elem_fixed_solution based methods
516 : */
517 : Real fixed_solution_derivative;
518 :
519 : /**
520 : * Points the _deltat member of this class at a timestep value
521 : * stored in the creating System, for example DiffSystem::deltat
522 : */
523 : void set_deltat_pointer(Real * dt);
524 :
525 : /**
526 : * \returns The value currently pointed to by this class's \p _deltat
527 : * member
528 : */
529 : Real get_deltat_value();
530 :
531 : /**
532 : * Adds a vector to the map of localized vectors. We can later evaluate interior_values,
533 : * interior_gradients and side_values for these fields these vectors represent.
534 : */
535 : void add_localized_vector (NumericVector<Number> & localized_vector, const System & sys);
536 :
537 : /**
538 : * Typedef for the localized_vectors iterator
539 : */
540 : typedef std::map<const NumericVector<Number> *, std::pair<DenseVector<Number>, std::vector<DenseSubVector<Number>>>>::iterator localized_vectors_iterator;
541 :
542 : /**
543 : * Return a reference to DenseVector localization of localized_vector
544 : * contained in the _localized_vectors map
545 : */
546 : DenseVector<Number> & get_localized_vector (const NumericVector<Number> & localized_vector);
547 :
548 : /**
549 : * const accessible version of get_localized_vector function
550 : */
551 : const DenseVector<Number> & get_localized_vector (const NumericVector<Number> & localized_vector) const;
552 :
553 : /**
554 : * Return a reference to DenseSubVector localization of localized_vector at variable var
555 : * contained in the _localized_vectors map
556 : */
557 : DenseSubVector<Number> & get_localized_subvector (const NumericVector<Number> & localized_vector, unsigned int var);
558 :
559 : /**
560 : * const accessible version of get_localized_subvector function
561 : */
562 : const DenseSubVector<Number> & get_localized_subvector (const NumericVector<Number> & localized_vector, unsigned int var) const;
563 :
564 : protected:
565 :
566 : /**
567 : * Contains pointers to vectors the user has asked to be localized, keyed with
568 : * pairs of element localized versions of that vector and per variable views
569 : */
570 :
571 : std::map<const NumericVector<Number> *, std::pair<DenseVector<Number>, std::vector<DenseSubVector<Number>>>> _localized_vectors;
572 :
573 : /**
574 : * Whether we have local matrices allocated/initialized
575 : */
576 : const bool _have_local_matrices;
577 :
578 : /**
579 : * Element by element components of nonlinear_solution
580 : * as adjusted by a time_solver
581 : */
582 : DenseVector<Number> _elem_solution;
583 : std::vector<DenseSubVector<Number>> _elem_subsolutions;
584 :
585 : /**
586 : * Element by element components of du/dt
587 : * as adjusted by a time_solver
588 : */
589 : DenseVector<Number> _elem_solution_rate;
590 : std::vector<DenseSubVector<Number>> _elem_subsolution_rates;
591 :
592 : /**
593 : * Element by element components of du/dt
594 : * as adjusted by a time_solver
595 : */
596 : DenseVector<Number> _elem_solution_accel;
597 : std::vector<DenseSubVector<Number>> _elem_subsolution_accels;
598 :
599 : /**
600 : * Element by element components of nonlinear_solution
601 : * at a fixed point in a timestep, for optional use by e.g.
602 : * stabilized methods
603 : */
604 : DenseVector<Number> _elem_fixed_solution;
605 : std::vector<DenseSubVector<Number>> _elem_fixed_subsolutions;
606 :
607 : /**
608 : * Element residual vector
609 : */
610 : DenseVector<Number> _elem_residual;
611 :
612 : /**
613 : * Element jacobian: derivatives of elem_residual with respect to
614 : * elem_solution. Only initialized if \p _have_local_matrices
615 : */
616 : DenseMatrix<Number> _elem_jacobian;
617 :
618 : /**
619 : * Element quantity of interest contributions
620 : */
621 : std::vector<Number> _elem_qoi;
622 :
623 : /**
624 : * Element quantity of interest derivative contributions
625 : */
626 : std::vector<DenseVector<Number>> _elem_qoi_derivative;
627 : std::vector<std::vector<DenseSubVector<Number>>> _elem_qoi_subderivatives;
628 :
629 : /**
630 : * Element residual subvectors and (if \p _have_local_matrices)
631 : * Jacobian submatrices
632 : */
633 : std::vector<DenseSubVector<Number>> _elem_subresiduals;
634 : std::vector<std::vector<DenseSubMatrix<Number>>> _elem_subjacobians;
635 :
636 : /**
637 : * Global Degree of freedom index lists
638 : */
639 : std::vector<dof_id_type> _dof_indices;
640 : std::vector<std::vector<dof_id_type>> _dof_indices_var;
641 :
642 : private:
643 :
644 : /**
645 : * Defaults to nullptr, can optionally be used to point to a timestep value
646 : * in the System-derived class responsible for creating this DiffContext.
647 : *
648 : * In DiffSystem's build_context() function, is assigned to point to
649 : * the deltat member of that class.
650 : *
651 : * Accessible via public get_deltat()/set_deltat() methods of this class.
652 : *
653 : * Always test for nullptr before using!
654 : */
655 : Real * _deltat;
656 :
657 : /**
658 : * A reference to the system this context is constructed with
659 : */
660 : const System & _system;
661 :
662 : /**
663 : * Is this context to be used for a primal or adjoint solve?
664 : */
665 : bool _is_adjoint;
666 :
667 : };
668 :
669 : } // namespace libMesh
670 :
671 :
672 : #endif // LIBMESH_DIFF_CONTEXT_H
|
