deal.II version GIT relicensing-2233-g9bd70d8389 2024-12-08 16:10:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
ad_helpers.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2018 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#include <deal.II/base/config.h>
16
17#if defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO)
18
21
22# include <type_traits>
23
24
26
27
28namespace Differentiation
29{
30 namespace AD
31 {
32 /* -------------------------- HelperBase -------------------------- */
33
34
35
36 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
38 const unsigned int n_independent_variables,
39 const unsigned int n_dependent_variables)
40 : independent_variable_values(
41 n_independent_variables,
42 ::internal::NumberType<scalar_type>::value(0.0))
43 , registered_independent_variable_values(n_independent_variables, false)
44 , registered_marked_independent_variables(n_independent_variables, false)
45 , registered_marked_dependent_variables(n_dependent_variables, false)
46 {
47 // We have enabled the compilation of this class for arithmetic
48 // types (i.e. ADNumberTypeCode == NumberTypes::none), but we
49 // can't actually do anything with them. Lets not advance any further
50 // and seemingly allow any operations that will not give any
51 // sensible results.
52 Assert(ADNumberTypeCode != NumberTypes::none,
54 "Floating point/arithmetic numbers have no derivatives."));
55 Assert(
58 "The AD number type does not support the calculation of any derivatives."));
59
60 // Tapeless mode must be configured before any active live
61 // variables are created.
63 {
65 false /*ensure_persistent_setting*/);
66 }
67
68 // For safety, we ensure that the entries in this vector *really* are
69 // initialized correctly by sending in the constructed zero-valued
70 // initializer.
73 0.0));
74 }
75
76
77
78 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
79 void
80 HelperBase<ADNumberTypeCode,
81 ScalarType>::reset_registered_independent_variables()
82 {
83 std::fill(registered_independent_variable_values.begin(),
84 registered_independent_variable_values.end(),
85 false);
86 }
87
88
89
90 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
91 void
94 {
95 std::fill(registered_marked_dependent_variables.begin(),
96 registered_marked_dependent_variables.end(),
97 flag);
98 }
99
100
101
102 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
103 void
105 const unsigned int index,
106 const scalar_type &value)
107 {
109 {
110 // A dummy call in case the user does not encapsulate a set of
111 // calls to tapeless ADHelpers with the initial and final calls to
112 // [start,stop]_recording_operations.
113 if (this->is_recording() == false)
114 start_recording_operations(1 /*tape index*/);
115
116 Assert(this->is_recording() == true,
118 "Cannot change the value of an independent variable "
119 "of the tapeless variety while this class is not set "
120 "in recording operations."));
121 }
123 {
124 Assert(this->active_tape_index() !=
126 ExcMessage("Invalid tape index"));
127 }
128 Assert(
129 index < n_independent_variables(),
131 "Trying to set the value of a non-existent independent variable."));
132
133 independent_variable_values[index] = value;
134 registered_independent_variable_values[index] = true;
135 }
136
137
138
139 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
140 void
142 const unsigned int index,
143 ad_type &out) const
144 {
145 Assert(index < n_independent_variables(), ExcInternalError());
146 Assert(registered_independent_variable_values[index] == true,
148
149 if (index > 0)
150 {
151 Assert(
152 registered_marked_independent_variables[index - 1] == true,
154 "Need to extract sensitivities in the order they're created."));
155 }
156
158 {
159 Assert(active_tape_index() != Numbers<ad_type>::invalid_tape_index,
160 ExcMessage("Invalid tape index"));
161 Assert(is_recording() == true,
163 "The marking of independent variables is only valid "
164 "during recording."));
165 }
166
168 independent_variable_values[index],
169 index,
170 this->n_independent_variables(),
171 out);
172 registered_marked_independent_variables[index] = true;
173 }
174
175
176
177 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
178 void
179 HelperBase<ADNumberTypeCode,
180 ScalarType>::finalize_sensitive_independent_variables() const
181 {
182 // Double check that we've actually registered all DoFs
183 Assert(n_registered_independent_variables() == n_independent_variables(),
184 ExcMessage("Not all values of sensitivities have been recorded!"));
185
186 // This should happen only once
187 if (this->independent_variables.empty())
188 {
189 this->independent_variables.resize(
190 this->n_independent_variables(),
192
193 // Indicate the sensitivity that each entry represents
194 for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
195 this->mark_independent_variable(i, this->independent_variables[i]);
196 }
197 }
198
199
200
201 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
202 void
205 ad_type &out) const
206 {
208 {
209 Assert(active_tape_index() != Numbers<ad_type>::invalid_tape_index,
210 ExcMessage("Invalid tape index"));
211 }
212 Assert(is_recording() == false,
214 "The initialization of non-sensitive independent variables is "
215 "only valid outside of recording operations."));
216
217 Assert(index < n_independent_variables(), ExcInternalError());
218 Assert(registered_independent_variable_values[index] == true,
220
221 out = independent_variable_values[index];
222 }
223
225
226 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
227 unsigned int
228 HelperBase<ADNumberTypeCode,
229 ScalarType>::n_registered_independent_variables() const
230 {
231 return std::count(registered_independent_variable_values.begin(),
232 registered_independent_variable_values.end(),
233 true);
234 }
235
236
237
238 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
239 std::size_t
241 {
242 return independent_variable_values.size();
243 }
244
245
246
247 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
248 unsigned int
250 const
252 return std::count(registered_marked_dependent_variables.begin(),
253 registered_marked_dependent_variables.end(),
254 true);
255 }
256
257
258
259 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
260 std::size_t
262 {
263 return dependent_variables.size();
265
266
267
268 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
269 bool
271 {
273 return taped_driver.is_recording();
274 else
275 return tapeless_driver.is_dependent_variable_marking_allowed();
276 }
277
278
279
280 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
281 typename Types<
284 {
286 return taped_driver.active_tape_index();
287 else
288 return 1;
289 }
290
291
292
293 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
294 bool
296 const typename Types<ad_type>::tape_index tape_index) const
297 {
299 return taped_driver.is_registered_tape(tape_index);
300 else
301 return true;
302 }
303
304
305
306 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
307 void
309 {
310 // Store stream flags
311 const std::ios_base::fmtflags stream_flags(stream.flags());
312 // Set stream to print booleans as "true"/"false"
313 stream.setf(std::ios_base::boolalpha);
314
315 stream << "Active tape index: " << active_tape_index() << '\n';
316 stream << "Recording? " << is_recording() << '\n';
317 stream << std::flush;
318
320 taped_driver.print(stream);
321
322 stream << "Registered independent variables: " << '\n';
323 for (unsigned int i = 0; i < n_independent_variables(); ++i)
324 stream << registered_independent_variable_values[i]
325 << (i < (n_independent_variables() - 1) ? "," : "");
326 stream << std::endl;
327
328 stream << "Independent variable values: " << '\n';
329 print_values(stream);
330
331 stream << "Registered marked independent variables: " << '\n';
332 for (unsigned int i = 0; i < n_independent_variables(); ++i)
333 stream << registered_marked_independent_variables[i]
334 << (i < (n_independent_variables() - 1) ? "," : "")
335 << std::flush;
336 stream << std::endl;
337
338 stream << "Dependent variable values: " << '\n';
339 for (unsigned int i = 0; i < n_dependent_variables(); ++i)
340 stream << dependent_variables[i]
341 << (i < (n_dependent_variables() - 1) ? "," : "");
342 stream << std::endl;
343
344 stream << "Registered dependent variables: " << '\n';
345 for (unsigned int i = 0; i < n_dependent_variables(); ++i)
346 stream << registered_marked_dependent_variables[i]
347 << (i < (n_dependent_variables() - 1) ? "," : "");
348 stream << std::endl;
349
350 // Restore stream flags
351 stream.flags(stream_flags);
352 }
353
354
355
356 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
357 void
359 std::ostream &stream) const
360 {
361 for (unsigned int i = 0; i < n_independent_variables(); ++i)
362 stream << independent_variable_values[i]
363 << (i < (n_independent_variables() - 1) ? "," : "")
364 << std::flush;
365 stream << std::endl;
366 }
367
368
369
370 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
371 void
373 const typename Types<ad_type>::tape_index tape_index,
374 std::ostream &stream) const
375 {
377 return;
378
379 Assert(is_registered_tape(tape_index),
380 ExcMessage("Tape number not registered"));
381
382 this->taped_driver.print_tape_stats(tape_index, stream);
383 }
384
385
386
387 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
388 void
390 const unsigned int n_independent_variables,
391 const unsigned int n_dependent_variables,
392 const bool clear_registered_tapes)
394 const unsigned int new_n_independent_variables =
395 (n_independent_variables != ::numbers::invalid_unsigned_int ?
396 n_independent_variables :
397 this->n_independent_variables());
398 const unsigned int new_n_dependent_variables =
399 (n_dependent_variables != ::numbers::invalid_unsigned_int ?
400 n_dependent_variables :
401 this->n_dependent_variables());
402
403 // Here we clear our vectors of AD data with a reallocation of memory
404 // (i.e. we *nuke* them entirely). Why we do this differs for each
405 // AD type:
406 // - ADOL-C taped mode must have their tapes fully cleared of marked data
407 // before the tapes can be overwritten.
408 // - ADOL-C tapeless mode must be configured for before any active live
409 // variables are created.
410 // - Reverse-mode Sacado numbers to must be destroyed to reset their
411 // accumulations.
412 // - Forward-mode Sacado numbers have no specific requirements, but it
413 // doesn't really hurt to perform this operation anyway.
414 {
415 std::vector<ad_type>().swap(independent_variables);
416 std::vector<ad_type>().swap(dependent_variables);
417 }
418
419 // Tapeless mode must be configured before any active live
420 // variables are created.
422 {
423 configure_tapeless_mode(new_n_independent_variables,
424 false /*ensure_persistent_setting*/);
425 }
427 taped_driver.reset(clear_registered_tapes);
428
429 independent_variable_values = std::vector<scalar_type>(
430 new_n_independent_variables,
432 registered_independent_variable_values =
433 std::vector<bool>(new_n_independent_variables, false);
434 registered_marked_independent_variables =
435 std::vector<bool>(new_n_independent_variables, false);
436 dependent_variables =
437 std::vector<ad_type>(new_n_dependent_variables,
439 registered_marked_dependent_variables =
440 std::vector<bool>(new_n_dependent_variables, false);
441 }
442
443
444
445 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
446 void
448 const unsigned int n_independent_variables,
449 const bool ensure_persistent_setting)
450 {
452 return;
453
454 // Try to safely initialize the global environment
456 n_independent_variables);
457
458 if (ensure_persistent_setting == true)
460 {
461 // In order to ensure that the settings remain for the entire
462 // duration of the simulation, we create a global live variable
463 // that doesn't go out of scope.
464 static ad_type num = 0.0;
465 (void)num;
467 }
468
469
470
471 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
472 void
474 const typename Types<ad_type>::tape_index tape_index)
475 {
476 activate_tape(tape_index, true /*read_mode*/);
477 }
479
480
481 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
482 bool
484 const typename Types<ad_type>::tape_index tape_index) const
485 {
487 return false;
488
489 return taped_driver.requires_retaping(tape_index);
490 }
491
492
493
494 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
495 bool
497 const
498 {
500 return false;
501
502 return taped_driver.last_action_requires_retaping();
503 }
504
505
506
507 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
508 void
510 {
512 return;
513
514 taped_driver.remove_tape(taped_driver.active_tape_index());
515 }
516
517
518
519 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
520 void
522 const typename Types<ad_type>::tape_index tape_index,
523 const bool read_mode)
524 {
526 {
528 ExcMessage("Invalid tape index"));
530 ExcMessage("Tape index exceeds maximum allowable value"));
531 taped_driver.activate_tape(tape_index);
532 reset_registered_independent_variables();
533
534 // A tape may have been defined by a different ADHelper, so in this
535 // case we ignore the fact that any dependent variables within the
536 // current data structure have not been marked as dependents
537 if (read_mode == true)
538 {
539 Assert(is_registered_tape(tape_index),
540 ExcMessage("Tape number not registered"));
541 reset_registered_dependent_variables(true);
542 Assert(n_registered_dependent_variables() ==
543 n_dependent_variables(),
544 ExcMessage("Not all dependent variables have been set!"));
545 }
546 }
547 else
548 {
551 }
552 }
553
554
555
556 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
557 void
559 const typename Types<ad_type>::tape_buffer_sizes obufsize,
560 const typename Types<ad_type>::tape_buffer_sizes lbufsize,
561 const typename Types<ad_type>::tape_buffer_sizes vbufsize,
562 const typename Types<ad_type>::tape_buffer_sizes tbufsize)
563 {
564 // When valid for the chosen AD number type, these values will be used the
565 // next time start_recording_operations() is called.
567 taped_driver.set_tape_buffer_sizes(obufsize,
568 lbufsize,
569 vbufsize,
570 tbufsize);
571 }
572
573
575 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
576 bool
578 const typename Types<ad_type>::tape_index tape_index,
579 const bool overwrite_tape,
580 const bool keep_independent_values)
581 {
582 // Define this here for clarity when this flag is used later.
583 const bool read_mode = false;
584
586 {
587 if (overwrite_tape != true)
588 {
589 Assert(is_recording() == false,
590 ExcMessage("Already recording..."));
591 }
592
593 // Check conditions to enable tracing
594 if (is_registered_tape(tape_index) == false || overwrite_tape == true)
595 {
596 // Setup the data structures for this class in the
597 // appropriate manner
598 activate_tape(tape_index, read_mode);
599
600 // Start taping
601 taped_driver.start_taping(active_tape_index(),
602 keep_independent_values);
603
604 // Clear the flags that state which independent and
605 // dependent variables have been registered
606 reset_registered_independent_variables();
607 reset_registered_dependent_variables();
608 }
609 else
610 {
611 Assert(is_recording() == false,
613 "Tape recording is unexpectedly still enabled."));
614
615 // Now we activate the pre-recorded tape so that its immediately
616 // available for use
617 activate_recorded_tape(tape_index);
618 }
619 }
620 else
621 {
622 Assert(ADNumberTraits<ad_type>::is_tapeless == true,
624
625 // Set the flag that states that we can safely mark dependent
626 // variables within this current phase of operations
627 tapeless_driver.allow_dependent_variable_marking();
628
629 // Dummy call to ensure that the intuitively correct
630 // value for the active tape (whether "valid" or not)
631 // is always returned to the user.
632 activate_tape(tape_index, read_mode);
633 }
634
635 return is_recording();
636 }
637
638
639
640 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
641 void
643 const bool write_tapes_to_file)
644 {
645 Assert(is_recording() == true, ExcMessage("Not currently recording..."));
646
647 // Double check that we've actually registered all DoFs
648 Assert(n_registered_independent_variables() == n_independent_variables(),
649 ExcMessage("Not all values of sensitivities have been recorded!"));
650
652 {
653 // Stop tracing
654 taped_driver.stop_taping(active_tape_index(), write_tapes_to_file);
655 }
656 else
657 {
660 // Double check that we've actually registered dependent variables
661 Assert(n_registered_dependent_variables() == n_dependent_variables(),
662 ExcMessage("Not all dependent variables have been set!"));
663
664 // By changing this flag, we ensure that the we can no longer
665 // legally alter the values of the dependent variables using
666 // set_dependent_variable(). This is important because the value of
667 // the tapeless independent variables are set and finalized when
668 // mark_dependent_variable() is called. So we cannot allow this to
669 // be done when not in the "recording" phase.
670 tapeless_driver.prevent_dependent_variable_marking();
671 }
672 }
674
675
676 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
677 void
679 const unsigned int index,
680 const ad_type &func)
681 {
682 Assert(index < n_dependent_variables(), ExcMessage("Index out of range"));
683 Assert(registered_marked_dependent_variables[index] == false,
685 "This dependent variable has already been registered."));
686
688 {
689 Assert(active_tape_index() != Numbers<ad_type>::invalid_tape_index,
690 ExcMessage("Invalid tape index"));
691 Assert(is_recording() == true,
693 "Must be recording when registering dependent variables."));
695
696 // Register the given dependent variable
697 internal::Marking<ad_type>::dependent_variable(dependent_variables[index],
698 func);
699 registered_marked_dependent_variables[index] = true;
700 }
701
702
703
704 /* -------------------- CellLevelBase -------------------- */
705
706
708 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
710 const unsigned int n_independent_variables,
711 const unsigned int n_dependent_variables)
712 : HelperBase<ADNumberTypeCode, ScalarType>(n_independent_variables,
713 n_dependent_variables)
714 {}
715
716
717
718 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
719 void
721 const std::vector<scalar_type> &dof_values)
722 {
723 // This is actually the same thing the set_independent_variable function,
724 // in the sense that we simply populate our array of independent values
725 // with a meaningful number. However, in this case we need to double check
726 // that we're not registering these variables twice
727 Assert(dof_values.size() == this->n_independent_variables(),
729 "Vector size does not match number of independent variables"));
730 for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
731 {
732 Assert(this->registered_independent_variable_values[i] == false,
733 ExcMessage("Independent variable value already registered."));
734 }
735 set_dof_values(dof_values);
736 }
737
738
739
740 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
741 const std::vector<
744 const
745 {
747 {
748 Assert(this->active_tape_index() !=
750 ExcMessage("Invalid tape index"));
751 }
752
753 // If necessary, initialize the internally stored vector of
754 // AD numbers that represents the independent variables
755 this->finalize_sensitive_independent_variables();
756 Assert(this->independent_variables.size() ==
757 this->n_independent_variables(),
758 ExcDimensionMismatch(this->independent_variables.size(),
759 this->n_independent_variables()));
760
761 return this->independent_variables;
763
764
765
766 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
767 void
769 const std::vector<scalar_type> &values)
770 {
772 {
773 Assert(this->active_tape_index() !=
775 ExcMessage("Invalid tape index"));
776 }
777 Assert(values.size() == this->n_independent_variables(),
779 "Vector size does not match number of independent variables"));
780 for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
782 i, values[i]);
783 }
784
785
786
787 /* ------------------ EnergyFunctional ------------------ */
788
789
790
791 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
793 const unsigned int n_independent_variables)
794 : CellLevelBase<ADNumberTypeCode, ScalarType>(n_independent_variables, 1)
795 {}
796
797
798
799 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
800 void
808
809
810
811 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
814 {
816 this->taped_driver.keep_independent_values() == false) ||
818 {
819 Assert(
820 this->n_registered_independent_variables() ==
821 this->n_independent_variables(),
823 "Not all values of sensitivities have been registered or subsequently set!"));
824 }
825 Assert(this->n_registered_dependent_variables() ==
826 this->n_dependent_variables(),
827 ExcMessage("Not all dependent variables have been registered."));
828
829 Assert(
830 this->n_dependent_variables() == 1,
832 "The EnergyFunctional class expects there to be only one dependent variable."));
833
835 {
836 Assert(this->active_tape_index() !=
838 ExcMessage("Invalid tape index"));
839 Assert(this->is_recording() == false,
841 "Cannot compute value while tape is being recorded."));
842 Assert(this->independent_variable_values.size() ==
843 this->n_independent_variables(),
844 ExcDimensionMismatch(this->independent_variable_values.size(),
845 this->n_independent_variables()));
846
847 return this->taped_driver.value(this->active_tape_index(),
848 this->independent_variable_values);
849 }
850 else
851 {
854 Assert(this->independent_variables.size() ==
855 this->n_independent_variables(),
856 ExcDimensionMismatch(this->independent_variables.size(),
857 this->n_independent_variables()));
858
859 return this->tapeless_driver.value(this->dependent_variables);
860 }
861 }
862
863
864
865 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
866 void
868 Vector<scalar_type> &gradient) const
869 {
871 this->taped_driver.keep_independent_values() == false) ||
873 {
874 Assert(
875 this->n_registered_independent_variables() ==
876 this->n_independent_variables(),
878 "Not all values of sensitivities have been registered or subsequently set!"));
879 }
880 Assert(this->n_registered_dependent_variables() ==
881 this->n_dependent_variables(),
882 ExcMessage("Not all dependent variables have been registered."));
883
884 Assert(
885 this->n_dependent_variables() == 1,
887 "The EnergyFunctional class expects there to be only one dependent variable."));
888
889 // We can neglect correctly initializing the entries as
890 // we'll be overwriting them immediately in the succeeding call to
891 // Drivers::gradient().
892 if (gradient.size() != this->n_independent_variables())
893 gradient.reinit(this->n_independent_variables(),
894 true /*omit_zeroing_entries*/);
895
897 {
898 Assert(this->active_tape_index() !=
900 ExcMessage("Invalid tape index"));
901 Assert(this->is_recording() == false,
903 "Cannot compute gradient while tape is being recorded."));
904 Assert(this->independent_variable_values.size() ==
905 this->n_independent_variables(),
906 ExcDimensionMismatch(this->independent_variable_values.size(),
907 this->n_independent_variables()));
908
909 this->taped_driver.gradient(this->active_tape_index(),
910 this->independent_variable_values,
911 gradient);
912 }
913 else
914 {
917 Assert(this->independent_variables.size() ==
918 this->n_independent_variables(),
919 ExcDimensionMismatch(this->independent_variables.size(),
920 this->n_independent_variables()));
921
922 this->tapeless_driver.gradient(this->independent_variables,
923 this->dependent_variables,
924 gradient);
925 }
926 }
927
928
929
930 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
931 void
933 FullMatrix<scalar_type> &hessian) const
934 {
937 "Cannot computed function Hessian: AD number type does "
938 "not support the calculation of second order derivatives."));
939
941 this->taped_driver.keep_independent_values() == false))
942 {
943 Assert(
944 this->n_registered_independent_variables() ==
945 this->n_independent_variables(),
947 "Not all values of sensitivities have been registered or subsequently set!"));
948 }
949 Assert(this->n_registered_dependent_variables() ==
950 this->n_dependent_variables(),
951 ExcMessage("Not all dependent variables have been registered."));
952
953 Assert(
954 this->n_dependent_variables() == 1,
956 "The EnergyFunctional class expects there to be only one dependent variable."));
957
958 // We can neglect correctly initializing the entries as
959 // we'll be overwriting them immediately in the succeeding call to
960 // Drivers::hessian().
961 if (hessian.m() != this->n_independent_variables() ||
962 hessian.n() != this->n_independent_variables())
963 hessian.reinit({this->n_independent_variables(),
964 this->n_independent_variables()},
965 true /*omit_default_initialization*/);
966
968 {
969 Assert(this->active_tape_index() !=
971 ExcMessage("Invalid tape index"));
972 Assert(this->is_recording() == false,
974 "Cannot compute hessian while tape is being recorded."));
975 Assert(this->independent_variable_values.size() ==
976 this->n_independent_variables(),
977 ExcDimensionMismatch(this->independent_variable_values.size(),
978 this->n_independent_variables()));
980 this->taped_driver.hessian(this->active_tape_index(),
981 this->independent_variable_values,
982 hessian);
983 }
984 else
985 {
988 Assert(this->independent_variables.size() ==
989 this->n_independent_variables(),
990 ExcDimensionMismatch(this->independent_variables.size(),
991 this->n_independent_variables()));
992
993 this->tapeless_driver.hessian(this->independent_variables,
994 this->dependent_variables,
995 hessian);
996 }
997 }
998
999
1000 /* ------------------- ResidualLinearization ------------------- */
1001
1002
1003
1004 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
1006 const unsigned int n_independent_variables,
1007 const unsigned int n_dependent_variables)
1008 : CellLevelBase<ADNumberTypeCode, ScalarType>(n_independent_variables,
1009 n_dependent_variables)
1010 {}
1011
1012
1013
1014 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
1015 void
1017 register_residual_vector(const std::vector<ad_type> &residual)
1018 {
1019 Assert(residual.size() == this->n_dependent_variables(),
1020 ExcMessage(
1021 "Vector size does not match number of dependent variables"));
1022 for (unsigned int i = 0; i < this->n_dependent_variables(); ++i)
1024 i, residual[i]);
1025 }
1026
1027
1028
1029 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
1030 void
1032 Vector<scalar_type> &values) const
1033 {
1034 if ((ADNumberTraits<ad_type>::is_taped == true &&
1035 this->taped_driver.keep_independent_values() == false) ||
1037 {
1038 Assert(
1039 this->n_registered_independent_variables() ==
1040 this->n_independent_variables(),
1041 ExcMessage(
1042 "Not all values of sensitivities have been registered or subsequently set!"));
1043 }
1044 Assert(this->n_registered_dependent_variables() ==
1045 this->n_dependent_variables(),
1046 ExcMessage("Not all dependent variables have been registered."));
1047
1048 // We can neglect correctly initializing the entries as
1049 // we'll be overwriting them immediately in the succeeding call to
1050 // Drivers::values().
1051 if (values.size() != this->n_dependent_variables())
1052 values.reinit(this->n_dependent_variables(),
1053 true /*omit_zeroing_entries*/);
1054
1056 {
1057 Assert(this->active_tape_index() !=
1059 ExcMessage("Invalid tape index"));
1060 Assert(this->is_recording() == false,
1061 ExcMessage(
1062 "Cannot compute values while tape is being recorded."));
1063 Assert(this->independent_variable_values.size() ==
1064 this->n_independent_variables(),
1065 ExcDimensionMismatch(this->independent_variable_values.size(),
1066 this->n_independent_variables()));
1067
1068 this->taped_driver.values(this->active_tape_index(),
1069 this->n_dependent_variables(),
1070 this->independent_variable_values,
1071 values);
1072 }
1073 else
1074 {
1077 this->tapeless_driver.values(this->dependent_variables, values);
1078 }
1079 }
1080
1081
1082
1083 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
1084 void
1086 FullMatrix<scalar_type> &jacobian) const
1087 {
1088 if ((ADNumberTraits<ad_type>::is_taped == true &&
1089 this->taped_driver.keep_independent_values() == false) ||
1091 {
1092 Assert(
1093 this->n_registered_independent_variables() ==
1094 this->n_independent_variables(),
1095 ExcMessage(
1096 "Not all values of sensitivities have been registered or subsequently set!"));
1097 }
1098 Assert(this->n_registered_dependent_variables() ==
1099 this->n_dependent_variables(),
1100 ExcMessage("Not all dependent variables have been registered."));
1101
1102 // We can neglect correctly initializing the entries as
1103 // we'll be overwriting them immediately in the succeeding call to
1104 // Drivers::jacobian().
1105 if (jacobian.m() != this->n_dependent_variables() ||
1106 jacobian.n() != this->n_independent_variables())
1107 jacobian.reinit({this->n_dependent_variables(),
1108 this->n_independent_variables()},
1109 true /*omit_default_initialization*/);
1110
1112 {
1113 Assert(this->active_tape_index() !=
1115 ExcMessage("Invalid tape index"));
1116 Assert(this->is_recording() == false,
1117 ExcMessage(
1118 "Cannot compute hessian while tape is being recorded."));
1119 Assert(this->independent_variable_values.size() ==
1120 this->n_independent_variables(),
1121 ExcDimensionMismatch(this->independent_variable_values.size(),
1122 this->n_independent_variables()));
1123
1124 this->taped_driver.jacobian(this->active_tape_index(),
1125 this->n_dependent_variables(),
1126 this->independent_variable_values,
1127 jacobian);
1128 }
1129 else
1130 {
1133 Assert(this->independent_variables.size() ==
1134 this->n_independent_variables(),
1135 ExcDimensionMismatch(this->independent_variables.size(),
1136 this->n_independent_variables()));
1137
1138 this->tapeless_driver.jacobian(this->independent_variables,
1139 this->dependent_variables,
1140 jacobian);
1141 }
1142 }
1143
1144
1145
1146 /* ----------------- PointLevelFunctionsBase ----------------- */
1147
1148
1149
1150 template <int dim,
1151 enum AD::NumberTypes ADNumberTypeCode,
1152 typename ScalarType>
1154 PointLevelFunctionsBase(const unsigned int n_independent_variables,
1155 const unsigned int n_dependent_variables)
1156 : HelperBase<ADNumberTypeCode, ScalarType>(n_independent_variables,
1157 n_dependent_variables)
1158 , symmetric_independent_variables(n_independent_variables, false)
1159 {}
1160
1161
1162
1163 template <int dim,
1164 enum AD::NumberTypes ADNumberTypeCode,
1165 typename ScalarType>
1166 void
1168 const unsigned int n_independent_variables,
1169 const unsigned int n_dependent_variables,
1170 const bool clear_registered_tapes)
1171 {
1173 n_dependent_variables,
1174 clear_registered_tapes);
1175
1176 const unsigned int new_n_independent_variables =
1177 (n_independent_variables != ::numbers::invalid_unsigned_int ?
1178 n_independent_variables :
1179 this->n_independent_variables());
1180 symmetric_independent_variables =
1181 std::vector<bool>(new_n_independent_variables, false);
1182 }
1183
1184
1185
1186 template <int dim,
1187 enum AD::NumberTypes ADNumberTypeCode,
1188 typename ScalarType>
1189 bool
1191 is_symmetric_independent_variable(const unsigned int index) const
1192 {
1193 Assert(index < symmetric_independent_variables.size(),
1195 return symmetric_independent_variables[index];
1196 }
1197
1198
1199
1200 template <int dim,
1201 enum AD::NumberTypes ADNumberTypeCode,
1202 typename ScalarType>
1203 unsigned int
1206 {
1207 return std::count(symmetric_independent_variables.begin(),
1208 symmetric_independent_variables.end(),
1209 true);
1210 }
1211
1212
1213
1214 template <int dim,
1215 enum AD::NumberTypes ADNumberTypeCode,
1216 typename ScalarType>
1217 void
1219 register_independent_variables(const std::vector<scalar_type> &values)
1220 {
1221 // This is actually the same thing the set_independent_variable function,
1222 // in the sense that we simply populate our array of independent values
1223 // with a meaningful number. However, in this case we need to double check
1224 // that we're not registering these variables twice
1225 Assert(values.size() == this->n_independent_variables(),
1226 ExcMessage(
1227 "Vector size does not match number of independent variables"));
1228 for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
1229 {
1230 Assert(this->registered_independent_variable_values[i] == false,
1231 ExcMessage("Independent variable value already registered."));
1232 }
1233 set_independent_variables(values);
1234 }
1235
1236
1237
1238 template <int dim,
1239 enum AD::NumberTypes ADNumberTypeCode,
1240 typename ScalarType>
1241 const std::vector<typename PointLevelFunctionsBase<dim,
1242 ADNumberTypeCode,
1243 ScalarType>::ad_type> &
1246 {
1248 {
1249 Assert(this->active_tape_index() !=
1251 ExcMessage("Invalid tape index"));
1252 }
1253
1254 // Just in case the user has not done so, we repeat the call to
1255 // initialize the internally stored vector of AD numbers that
1256 // represents the independent variables.
1257 this->finalize_sensitive_independent_variables();
1258 Assert(this->independent_variables.size() ==
1259 this->n_independent_variables(),
1260 ExcDimensionMismatch(this->independent_variables.size(),
1261 this->n_independent_variables()));
1262
1263 return this->independent_variables;
1264 }
1265
1266
1267
1268 template <int dim,
1269 enum AD::NumberTypes ADNumberTypeCode,
1270 typename ScalarType>
1271 void
1273 set_sensitivity_value(const unsigned int index,
1274 const bool symmetric_component,
1275 const scalar_type &value)
1276 {
1278 value);
1279 Assert(
1280 index < this->n_independent_variables(),
1281 ExcMessage(
1282 "Trying to set the symmetry flag of a non-existent independent variable."));
1283 Assert(index < symmetric_independent_variables.size(),
1285 symmetric_independent_variables[index] = symmetric_component;
1286 }
1287
1288
1289
1290 template <int dim,
1291 enum AD::NumberTypes ADNumberTypeCode,
1292 typename ScalarType>
1293 void
1295 set_independent_variables(const std::vector<scalar_type> &values)
1296 {
1298 {
1299 Assert(this->active_tape_index() !=
1301 ExcMessage("Invalid tape index"));
1302 }
1303 Assert(values.size() == this->n_independent_variables(),
1304 ExcMessage(
1305 "Vector size does not match number of independent variables"));
1306 for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
1308 i, values[i]);
1309 }
1310
1311
1312
1313 /* -------------------- ScalarFunction -------------------- */
1314
1315
1316
1317 template <int dim,
1318 enum AD::NumberTypes ADNumberTypeCode,
1319 typename ScalarType>
1321 const unsigned int n_independent_variables)
1322 : PointLevelFunctionsBase<dim, ADNumberTypeCode, ScalarType>(
1323 n_independent_variables,
1324 1)
1325 {}
1326
1327
1328
1329 template <int dim,
1330 enum AD::NumberTypes ADNumberTypeCode,
1331 typename ScalarType>
1332 void
1340
1341
1342
1343 template <int dim,
1344 enum AD::NumberTypes ADNumberTypeCode,
1345 typename ScalarType>
1348 {
1349 if ((ADNumberTraits<ad_type>::is_taped == true &&
1350 this->taped_driver.keep_independent_values() == false) ||
1352 {
1353 Assert(
1354 this->n_registered_independent_variables() ==
1355 this->n_independent_variables(),
1356 ExcMessage(
1357 "Not all values of sensitivities have been registered or subsequently set!"));
1358 }
1359 Assert(this->n_registered_dependent_variables() ==
1360 this->n_dependent_variables(),
1361 ExcMessage("Not all dependent variables have been registered."));
1362
1363 Assert(
1364 this->n_dependent_variables() == 1,
1365 ExcMessage(
1366 "The ScalarFunction class expects there to be only one dependent variable."));
1367
1369 {
1370 Assert(this->active_tape_index() !=
1372 ExcMessage("Invalid tape index"));
1373 Assert(this->is_recording() == false,
1374 ExcMessage(
1375 "Cannot compute values while tape is being recorded."));
1376 Assert(this->independent_variable_values.size() ==
1377 this->n_independent_variables(),
1378 ExcDimensionMismatch(this->independent_variable_values.size(),
1379 this->n_independent_variables()));
1380
1381 return this->taped_driver.value(this->active_tape_index(),
1382 this->independent_variable_values);
1383 }
1384 else
1385 {
1388 return this->tapeless_driver.value(this->dependent_variables);
1389 }
1390 }
1391
1392
1393 template <int dim,
1394 enum AD::NumberTypes ADNumberTypeCode,
1395 typename ScalarType>
1396 void
1398 Vector<scalar_type> &gradient) const
1399 {
1400 if ((ADNumberTraits<ad_type>::is_taped == true &&
1401 this->taped_driver.keep_independent_values() == false) ||
1403 {
1404 Assert(
1405 this->n_registered_independent_variables() ==
1406 this->n_independent_variables(),
1407 ExcMessage(
1408 "Not all values of sensitivities have been registered or subsequently set!"));
1409 }
1410 Assert(this->n_registered_dependent_variables() ==
1411 this->n_dependent_variables(),
1412 ExcMessage("Not all dependent variables have been registered."));
1413
1414 Assert(
1415 this->n_dependent_variables() == 1,
1416 ExcMessage(
1417 "The ScalarFunction class expects there to be only one dependent variable."));
1418
1419 // We can neglect correctly initializing the entries as
1420 // we'll be overwriting them immediately in the succeeding call to
1421 // Drivers::gradient().
1422 if (gradient.size() != this->n_independent_variables())
1423 gradient.reinit(this->n_independent_variables(),
1424 true /*omit_zeroing_entries*/);
1425
1427 {
1428 Assert(this->active_tape_index() !=
1430 ExcMessage("Invalid tape index"));
1431 Assert(this->is_recording() == false,
1432 ExcMessage(
1433 "Cannot compute gradient while tape is being recorded."));
1434 Assert(this->independent_variable_values.size() ==
1435 this->n_independent_variables(),
1436 ExcDimensionMismatch(this->independent_variable_values.size(),
1437 this->n_independent_variables()));
1438
1439 this->taped_driver.gradient(this->active_tape_index(),
1440 this->independent_variable_values,
1441 gradient);
1442 }
1443 else
1444 {
1447 Assert(this->independent_variables.size() ==
1448 this->n_independent_variables(),
1449 ExcDimensionMismatch(this->independent_variables.size(),
1450 this->n_independent_variables()));
1451
1452 this->tapeless_driver.gradient(this->independent_variables,
1453 this->dependent_variables,
1454 gradient);
1455 }
1456
1457 // Account for symmetries of tensor components
1458 for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
1459 {
1460 if (this->is_symmetric_independent_variable(i) == true)
1461 gradient[i] *= 0.5;
1462 }
1463 }
1464
1465
1466
1467 template <int dim,
1468 enum AD::NumberTypes ADNumberTypeCode,
1469 typename ScalarType>
1470 void
1472 FullMatrix<scalar_type> &hessian) const
1473 {
1475 ExcMessage(
1476 "Cannot computed function Hessian: AD number type does "
1477 "not support the calculation of second order derivatives."));
1478
1479 if ((ADNumberTraits<ad_type>::is_taped == true &&
1480 this->taped_driver.keep_independent_values() == false))
1481 {
1482 Assert(
1483 this->n_registered_independent_variables() ==
1484 this->n_independent_variables(),
1485 ExcMessage(
1486 "Not all values of sensitivities have been registered or subsequently set!"));
1487 }
1488 Assert(this->n_registered_dependent_variables() ==
1489 this->n_dependent_variables(),
1490 ExcMessage("Not all dependent variables have been registered."));
1491
1492 Assert(
1493 this->n_dependent_variables() == 1,
1494 ExcMessage(
1495 "The ScalarFunction class expects there to be only one dependent variable."));
1496
1497 // We can neglect correctly initializing the entries as
1498 // we'll be overwriting them immediately in the succeeding call to
1499 // Drivers::hessian().
1500 if (hessian.m() != this->n_independent_variables() ||
1501 hessian.n() != this->n_independent_variables())
1502 hessian.reinit({this->n_independent_variables(),
1503 this->n_independent_variables()},
1504 true /*omit_default_initialization*/);
1505
1507 {
1508 Assert(this->active_tape_index() !=
1510 ExcMessage("Invalid tape index"));
1511 Assert(this->is_recording() == false,
1512 ExcMessage(
1513 "Cannot compute Hessian while tape is being recorded."));
1514 Assert(this->independent_variable_values.size() ==
1515 this->n_independent_variables(),
1516 ExcDimensionMismatch(this->independent_variable_values.size(),
1517 this->n_independent_variables()));
1518
1519 this->taped_driver.hessian(this->active_tape_index(),
1520 this->independent_variable_values,
1521 hessian);
1522 }
1523 else
1524 {
1527 Assert(this->independent_variables.size() ==
1528 this->n_independent_variables(),
1529 ExcDimensionMismatch(this->independent_variables.size(),
1530 this->n_independent_variables()));
1531
1532 this->tapeless_driver.hessian(this->independent_variables,
1533 this->dependent_variables,
1534 hessian);
1535 }
1536
1537 // Account for symmetries of tensor components
1538 for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
1539 for (unsigned int j = 0; j < i + 1; ++j)
1540 {
1541 if (this->is_symmetric_independent_variable(i) == true &&
1542 this->is_symmetric_independent_variable(j) == true)
1543 {
1544 hessian[i][j] *= 0.25;
1545 if (i != j)
1546 hessian[j][i] *= 0.25;
1547 }
1548 else if ((this->is_symmetric_independent_variable(i) == true &&
1549 this->is_symmetric_independent_variable(j) == false) ||
1550 (this->is_symmetric_independent_variable(j) == true &&
1551 this->is_symmetric_independent_variable(i) == false))
1552 {
1553 hessian[i][j] *= 0.5;
1554 if (i != j)
1555 hessian[j][i] *= 0.5;
1556 }
1557 }
1558 }
1559
1560
1561
1562 template <int dim,
1563 enum AD::NumberTypes ADNumberTypeCode,
1564 typename ScalarType>
1565 Tensor<
1566 0,
1567 dim,
1571 const FEValuesExtractors::Scalar &extractor_row,
1572 const FEValuesExtractors::Scalar &extractor_col)
1573 {
1574 // NOTE: It is necessary to make special provision for the case when the
1575 // HessianType is scalar. Unfortunately Tensor<0,dim> does not provide
1576 // the function unrolled_to_component_indices!
1577 // NOTE: The order of components must be consistently defined throughout
1578 // this class.
1580
1581 // Get indexsets for the subblocks from which we wish to extract the
1582 // matrix values
1583 const std::vector<unsigned int> row_index_set(
1584 internal::extract_field_component_indices<dim>(extractor_row));
1585 const std::vector<unsigned int> col_index_set(
1586 internal::extract_field_component_indices<dim>(extractor_col));
1587 Assert(row_index_set.size() == 1, ExcInternalError());
1588 Assert(col_index_set.size() == 1, ExcInternalError());
1589
1591 0,
1592 hessian[row_index_set[0]][col_index_set[0]]);
1593
1594 return out;
1595 }
1596
1597
1598
1599 template <int dim,
1600 enum AD::NumberTypes ADNumberTypeCode,
1601 typename ScalarType>
1603 4,
1604 dim,
1608 const FullMatrix<scalar_type> &hessian,
1609 const FEValuesExtractors::SymmetricTensor<2> &extractor_row,
1610 const FEValuesExtractors::SymmetricTensor<2> &extractor_col)
1611 {
1612 // NOTE: The order of components must be consistently defined throughout
1613 // this class.
1614 // NOTE: We require a specialisation for rank-4 symmetric tensors because
1615 // they do not define their rank, and setting data using TableIndices is
1616 // somewhat specialised as well.
1618
1619 // Get indexsets for the subblocks from which we wish to extract the
1620 // matrix values
1621 const std::vector<unsigned int> row_index_set(
1622 internal::extract_field_component_indices<dim>(extractor_row));
1623 const std::vector<unsigned int> col_index_set(
1624 internal::extract_field_component_indices<dim>(extractor_col));
1625
1626 for (unsigned int r = 0; r < row_index_set.size(); ++r)
1627 for (unsigned int c = 0; c < col_index_set.size(); ++c)
1628 {
1630 out, r, c, hessian[row_index_set[r]][col_index_set[c]]);
1631 }
1632
1633 return out;
1634 }
1635
1636
1637
1638 /* -------------------- VectorFunction -------------------- */
1639
1640
1641
1642 template <int dim,
1643 enum AD::NumberTypes ADNumberTypeCode,
1644 typename ScalarType>
1646 const unsigned int n_independent_variables,
1647 const unsigned int n_dependent_variables)
1648 : PointLevelFunctionsBase<dim, ADNumberTypeCode, ScalarType>(
1649 n_independent_variables,
1650 n_dependent_variables)
1651 {}
1652
1653
1654
1655 template <int dim,
1656 enum AD::NumberTypes ADNumberTypeCode,
1657 typename ScalarType>
1658 void
1660 register_dependent_variables(const std::vector<ad_type> &funcs)
1661 {
1662 Assert(funcs.size() == this->n_dependent_variables(),
1663 ExcMessage(
1664 "Vector size does not match number of dependent variables"));
1665 for (unsigned int i = 0; i < this->n_dependent_variables(); ++i)
1667 i, funcs[i]);
1668 }
1669
1670
1671
1672 template <int dim,
1673 enum AD::NumberTypes ADNumberTypeCode,
1674 typename ScalarType>
1675 void
1677 Vector<scalar_type> &values) const
1678 {
1679 if ((ADNumberTraits<ad_type>::is_taped == true &&
1680 this->taped_driver.keep_independent_values() == false) ||
1682 {
1683 Assert(
1684 this->n_registered_independent_variables() ==
1685 this->n_independent_variables(),
1686 ExcMessage(
1687 "Not all values of sensitivities have been registered or subsequently set!"));
1688 }
1689 Assert(this->n_registered_dependent_variables() ==
1690 this->n_dependent_variables(),
1691 ExcMessage("Not all dependent variables have been registered."));
1692
1693 // We can neglect correctly initializing the entries as
1694 // we'll be overwriting them immediately in the succeeding call to
1695 // Drivers::values().
1696 if (values.size() != this->n_dependent_variables())
1697 values.reinit(this->n_dependent_variables(),
1698 true /*omit_zeroing_entries*/);
1699
1701 {
1702 Assert(this->active_tape_index() !=
1704 ExcMessage("Invalid tape index"));
1705 Assert(this->is_recording() == false,
1706 ExcMessage(
1707 "Cannot compute values while tape is being recorded."));
1708 Assert(this->independent_variable_values.size() ==
1709 this->n_independent_variables(),
1710 ExcDimensionMismatch(this->independent_variable_values.size(),
1711 this->n_independent_variables()));
1712
1713 this->taped_driver.values(this->active_tape_index(),
1714 this->n_dependent_variables(),
1715 this->independent_variable_values,
1716 values);
1717 }
1718 else
1719 {
1722 this->tapeless_driver.values(this->dependent_variables, values);
1723 }
1724 }
1725
1726
1727
1728 template <int dim,
1729 enum AD::NumberTypes ADNumberTypeCode,
1730 typename ScalarType>
1731 void
1733 FullMatrix<scalar_type> &jacobian) const
1734 {
1735 if ((ADNumberTraits<ad_type>::is_taped == true &&
1736 this->taped_driver.keep_independent_values() == false) ||
1738 {
1739 Assert(
1740 this->n_registered_independent_variables() ==
1741 this->n_independent_variables(),
1742 ExcMessage(
1743 "Not all values of sensitivities have been registered or subsequently set!"));
1744 }
1745 Assert(this->n_registered_dependent_variables() ==
1746 this->n_dependent_variables(),
1747 ExcMessage("Not all dependent variables have been registered."));
1748
1749 // We can neglect correctly initializing the entries as
1750 // we'll be overwriting them immediately in the succeeding call to
1751 // Drivers::jacobian().
1752 if (jacobian.m() != this->n_dependent_variables() ||
1753 jacobian.n() != this->n_independent_variables())
1754 jacobian.reinit({this->n_dependent_variables(),
1755 this->n_independent_variables()},
1756 true /*omit_default_initialization*/);
1757
1759 {
1760 Assert(this->active_tape_index() !=
1762 ExcMessage("Invalid tape index"));
1763 Assert(this->is_recording() == false,
1764 ExcMessage(
1765 "Cannot compute Jacobian while tape is being recorded."));
1766 Assert(this->independent_variable_values.size() ==
1767 this->n_independent_variables(),
1768 ExcDimensionMismatch(this->independent_variable_values.size(),
1769 this->n_independent_variables()));
1770
1771 this->taped_driver.jacobian(this->active_tape_index(),
1772 this->n_dependent_variables(),
1773 this->independent_variable_values,
1774 jacobian);
1775 }
1776 else
1777 {
1780 Assert(this->independent_variables.size() ==
1781 this->n_independent_variables(),
1782 ExcDimensionMismatch(this->independent_variables.size(),
1783 this->n_independent_variables()));
1784
1785 this->tapeless_driver.jacobian(this->independent_variables,
1786 this->dependent_variables,
1787 jacobian);
1788 }
1789
1790 for (unsigned int j = 0; j < this->n_independent_variables(); ++j)
1791 {
1792 // Because we perform just a single differentiation
1793 // operation with respect to the "column" variables,
1794 // we only need to consider them for symmetry conditions.
1795 if (this->is_symmetric_independent_variable(j) == true)
1796 for (unsigned int i = 0; i < this->n_dependent_variables(); ++i)
1797 jacobian[i][j] *= 0.5;
1798 }
1799 }
1800
1801
1802
1803 template <int dim,
1804 enum AD::NumberTypes ADNumberTypeCode,
1805 typename ScalarType>
1806 Tensor<
1807 0,
1808 dim,
1812 const FullMatrix<scalar_type> &jacobian,
1813 const FEValuesExtractors::Scalar &extractor_row,
1814 const FEValuesExtractors::Scalar &extractor_col)
1815 {
1816 // NOTE: It is necessary to make special provision for the case when the
1817 // HessianType is scalar. Unfortunately Tensor<0,dim> does not provide
1818 // the function unrolled_to_component_indices!
1819 // NOTE: The order of components must be consistently defined throughout
1820 // this class.
1822
1823 // Get indexsets for the subblocks from which we wish to extract the
1824 // matrix values
1825 const std::vector<unsigned int> row_index_set(
1826 internal::extract_field_component_indices<dim>(extractor_row));
1827 const std::vector<unsigned int> col_index_set(
1828 internal::extract_field_component_indices<dim>(extractor_col));
1829 Assert(row_index_set.size() == 1, ExcInternalError());
1830 Assert(col_index_set.size() == 1, ExcInternalError());
1831
1833 0,
1834 jacobian[row_index_set[0]][col_index_set[0]]);
1835
1836 return out;
1837 }
1838
1839
1840
1841 template <int dim,
1842 enum AD::NumberTypes ADNumberTypeCode,
1843 typename ScalarType>
1845 4,
1846 dim,
1850 const FullMatrix<scalar_type> &jacobian,
1851 const FEValuesExtractors::SymmetricTensor<2> &extractor_row,
1852 const FEValuesExtractors::SymmetricTensor<2> &extractor_col)
1853 {
1854 // NOTE: The order of components must be consistently defined throughout
1855 // this class.
1856 // NOTE: We require a specialisation for rank-4 symmetric tensors because
1857 // they do not define their rank, and setting data using TableIndices is
1858 // somewhat specialised as well.
1860
1861 // Get indexsets for the subblocks from which we wish to extract the
1862 // matrix values
1863 const std::vector<unsigned int> row_index_set(
1864 internal::extract_field_component_indices<dim>(extractor_row));
1865 const std::vector<unsigned int> col_index_set(
1866 internal::extract_field_component_indices<dim>(extractor_col));
1867
1868 for (unsigned int r = 0; r < row_index_set.size(); ++r)
1869 for (unsigned int c = 0; c < col_index_set.size(); ++c)
1870 {
1872 out, r, c, jacobian[row_index_set[r]][col_index_set[c]]);
1873 }
1874
1875 return out;
1876 }
1877
1878
1879 } // namespace AD
1880} // namespace Differentiation
1881
1882
1883/* --- Explicit instantiations --- */
1884# include "ad_helpers.inst"
1885
1886# ifdef DEAL_II_WITH_ADOLC
1887# include "ad_helpers.inst1"
1888# endif
1889# ifdef DEAL_II_TRILINOS_WITH_SACADO
1890# include "ad_helpers.inst2"
1891# endif
1892
1893
1895
1896#endif // defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO)
const std::vector< ad_type > & get_sensitive_dof_values() const
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
Definition ad_helpers.h:849
void set_dof_values(const std::vector< scalar_type > &dof_values)
void register_dof_values(const std::vector< scalar_type > &dof_values)
CellLevelBase(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
void register_energy_functional(const ad_type &energy)
EnergyFunctional(const unsigned int n_independent_variables)
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
void compute_residual(Vector< scalar_type > &residual) const override
virtual void compute_linearization(FullMatrix< scalar_type > &linearization) const override
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
unsigned int n_registered_dependent_variables() const
void stop_recording_operations(const bool write_tapes_to_file=false)
typename AD::NumberTraits< ScalarType, ADNumberTypeCode >::scalar_type scalar_type
Definition ad_helpers.h:176
Types< ad_type >::tape_index active_tape_index() const
void print_tape_stats(const typename Types< ad_type >::tape_index tape_index, std::ostream &stream) const
std::size_t n_independent_variables() const
bool start_recording_operations(const typename Types< ad_type >::tape_index tape_index, const bool overwrite_tape=false, const bool keep_independent_values=true)
void mark_independent_variable(const unsigned int index, ad_type &out) const
bool is_registered_tape(const typename Types< ad_type >::tape_index tape_index) const
void reset_registered_dependent_variables(const bool flag=false)
Definition ad_helpers.cc:93
void activate_tape(const typename Types< ad_type >::tape_index tape_index, const bool read_mode)
void set_tape_buffer_sizes(const typename Types< ad_type >::tape_buffer_sizes obufsize=64 *1024 *1024, const typename Types< ad_type >::tape_buffer_sizes lbufsize=64 *1024 *1024, const typename Types< ad_type >::tape_buffer_sizes vbufsize=64 *1024 *1024, const typename Types< ad_type >::tape_buffer_sizes tbufsize=64 *1024 *1024)
void print(std::ostream &stream) const
void set_sensitivity_value(const unsigned int index, const scalar_type &value)
bool active_tape_requires_retaping() const
HelperBase(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition ad_helpers.cc:37
virtual void reset(const unsigned int n_independent_variables=::numbers::invalid_unsigned_int, const unsigned int n_dependent_variables=::numbers::invalid_unsigned_int, const bool clear_registered_tapes=true)
bool recorded_tape_requires_retaping(const typename Types< ad_type >::tape_index tape_index) const
std::vector< ad_type > dependent_variables
Definition ad_helpers.h:748
std::size_t n_dependent_variables() const
void print_values(std::ostream &stream) const
void register_dependent_variable(const unsigned int index, const ad_type &func)
static void configure_tapeless_mode(const unsigned int n_independent_variables, const bool ensure_persistent_setting=true)
void initialize_non_sensitive_independent_variable(const unsigned int index, ad_type &out) const
typename AD::NumberTraits< ScalarType, ADNumberTypeCode >::ad_type ad_type
Definition ad_helpers.h:183
void activate_recorded_tape(const typename Types< ad_type >::tape_index tape_index)
bool is_symmetric_independent_variable(const unsigned int index) const
PointLevelFunctionsBase(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
void set_independent_variables(const std::vector< scalar_type > &values)
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
unsigned int n_symmetric_independent_variables() const
void register_independent_variables(const std::vector< scalar_type > &values)
void set_sensitivity_value(const unsigned int index, const bool symmetric_component, const scalar_type &value)
const std::vector< ad_type > & get_sensitive_variables() const
virtual void reset(const unsigned int n_independent_variables=::numbers::invalid_unsigned_int, const unsigned int n_dependent_variables=::numbers::invalid_unsigned_int, const bool clear_registered_tapes=true) override
ResidualLinearization(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
virtual void compute_residual(Vector< scalar_type > &residual) const override
virtual void compute_linearization(FullMatrix< scalar_type > &linearization) const override
void register_residual_vector(const std::vector< ad_type > &residual)
static internal::ScalarFieldHessian< dim, scalar_type, ExtractorType_Row, ExtractorType_Col >::type extract_hessian_component(const FullMatrix< scalar_type > &hessian, const ExtractorType_Row &extractor_row, const ExtractorType_Col &extractor_col)
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
void register_dependent_variable(const ad_type &func)
ScalarFunction(const unsigned int n_independent_variables)
void compute_gradient(Vector< scalar_type > &gradient) const
void compute_hessian(FullMatrix< scalar_type > &hessian) const
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
void compute_values(Vector< scalar_type > &values) const
void register_dependent_variables(const std::vector< ad_type > &funcs)
void compute_jacobian(FullMatrix< scalar_type > &jacobian) const
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
static internal::VectorFieldJacobian< dim, scalar_type, ExtractorType_Row, ExtractorType_Col >::type extract_jacobian_component(const FullMatrix< scalar_type > &jacobian, const ExtractorType_Row &extractor_row, const ExtractorType_Col &extractor_col)
VectorFunction(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
size_type n() const
size_type m() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
void set_tensor_entry(TensorType &t, const unsigned int unrolled_index, const NumberType &value)
static const unsigned int invalid_unsigned_int
Definition types.h:220
static void initialize_global_environment(const unsigned int n_independent_variables)