Reference documentation for deal.II version 9.5.0
\(\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// Copyright (C) 2016 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE at
12// the top level of the deal.II distribution.
13//
14// ---------------------------------------------------------------------
15
16#include <deal.II/base/config.h>
17
18#if defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO)
19
22
23# include <type_traits>
24
25
27
28
29namespace Differentiation
30{
31 namespace AD
32 {
33 /* -------------------------- HelperBase -------------------------- */
34
35
36
37 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
39 const unsigned int n_independent_variables,
40 const unsigned int n_dependent_variables)
41 : independent_variable_values(
42 n_independent_variables,
43 ::internal::NumberType<scalar_type>::value(0.0))
44 , registered_independent_variable_values(n_independent_variables, false)
45 , registered_marked_independent_variables(n_independent_variables, false)
46 , registered_marked_dependent_variables(n_dependent_variables, false)
47 {
48 // We have enabled the compilation of this class for arithmetic
49 // types (i.e. ADNumberTypeCode == NumberTypes::none), but we
50 // can't actually do anything with them. Lets not advance any further
51 // and seemingly allow any operations that will not give any
52 // sensible results.
53 Assert(ADNumberTypeCode != NumberTypes::none,
55 "Floating point/arithmetic numbers have no derivatives."));
56 Assert(
59 "The AD number type does not support the calculation of any derivatives."));
60
61 // Tapeless mode must be configured before any active live
62 // variables are created.
64 {
66 false /*ensure_persistent_setting*/);
67 }
68
69 // For safety, we ensure that the entries in this vector *really* are
70 // initialized correctly by sending in the constructed zero-valued
71 // initializer.
74 0.0));
75 }
76
77
78
79 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
80 void
81 HelperBase<ADNumberTypeCode,
82 ScalarType>::reset_registered_independent_variables()
83 {
84 std::fill(registered_independent_variable_values.begin(),
85 registered_independent_variable_values.end(),
86 false);
87 }
88
89
90
91 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
92 void
95 {
96 std::fill(registered_marked_dependent_variables.begin(),
97 registered_marked_dependent_variables.end(),
98 flag);
99 }
100
101
102
103 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
104 void
106 const unsigned int index,
107 const scalar_type &value)
108 {
110 {
111 // A dummy call in case the user does not encapsulate a set of
112 // calls to tapeless ADHelpers with the initial and final calls to
113 // [start,stop]_recording_operations.
114 if (this->is_recording() == false)
115 start_recording_operations(1 /*tape index*/);
116
117 Assert(this->is_recording() == true,
119 "Cannot change the value of an independent variable "
120 "of the tapeless variety while this class is not set "
121 "in recording operations."));
122 }
124 {
125 Assert(this->active_tape_index() !=
127 ExcMessage("Invalid tape index"));
128 }
129 Assert(
130 index < n_independent_variables(),
132 "Trying to set the value of a non-existent independent variable."));
133
134 independent_variable_values[index] = value;
135 registered_independent_variable_values[index] = true;
136 }
137
138
139
140 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
141 void
143 const unsigned int index,
144 ad_type & out) const
145 {
146 Assert(index < n_independent_variables(), ExcInternalError());
147 Assert(registered_independent_variable_values[index] == true,
149
150 if (index > 0)
151 {
152 Assert(
153 registered_marked_independent_variables[index - 1] == true,
155 "Need to extract sensitivities in the order they're created."));
156 }
157
159 {
160 Assert(active_tape_index() != Numbers<ad_type>::invalid_tape_index,
161 ExcMessage("Invalid tape index"));
162 Assert(is_recording() == true,
164 "The marking of independent variables is only valid "
165 "during recording."));
166 }
167
169 independent_variable_values[index],
170 index,
171 this->n_independent_variables(),
172 out);
173 registered_marked_independent_variables[index] = true;
174 }
175
176
177
178 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
179 void
180 HelperBase<ADNumberTypeCode,
181 ScalarType>::finalize_sensitive_independent_variables() const
182 {
183 // Double check that we've actually registered all DoFs
184 Assert(n_registered_independent_variables() == n_independent_variables(),
185 ExcMessage("Not all values of sensitivities have been recorded!"));
186
187 // This should happen only once
188 if (this->independent_variables.size() == 0)
189 {
190 this->independent_variables.resize(
191 this->n_independent_variables(),
193
194 // Indicate the sensitivity that each entry represents
195 for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
196 this->mark_independent_variable(i, this->independent_variables[i]);
197 }
198 }
199
200
201
202 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
203 void
206 ad_type &out) const
207 {
209 {
210 Assert(active_tape_index() != Numbers<ad_type>::invalid_tape_index,
211 ExcMessage("Invalid tape index"));
212 }
213 Assert(is_recording() == false,
215 "The initialization of non-sensitive independent variables is "
216 "only valid outside of recording operations."));
217
218 Assert(index < n_independent_variables(), ExcInternalError());
219 Assert(registered_independent_variable_values[index] == true,
221
222 out = independent_variable_values[index];
223 }
224
226
227 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
228 unsigned int
229 HelperBase<ADNumberTypeCode,
230 ScalarType>::n_registered_independent_variables() const
231 {
232 return std::count(registered_independent_variable_values.begin(),
233 registered_independent_variable_values.end(),
234 true);
235 }
236
237
238
239 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
240 std::size_t
242 {
243 return independent_variable_values.size();
244 }
245
246
247
248 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
249 unsigned int
251 const
253 return std::count(registered_marked_dependent_variables.begin(),
254 registered_marked_dependent_variables.end(),
255 true);
256 }
257
258
259
260 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
261 std::size_t
263 {
264 return dependent_variables.size();
266
267
268
269 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
270 bool
272 {
274 return taped_driver.is_recording();
275 else
276 return tapeless_driver.is_dependent_variable_marking_allowed();
277 }
278
279
280
281 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
282 typename Types<
285 {
287 return taped_driver.active_tape_index();
288 else
289 return 1;
290 }
291
292
293
294 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
295 bool
297 const typename Types<ad_type>::tape_index tape_index) const
298 {
300 return taped_driver.is_registered_tape(tape_index);
301 else
302 return true;
303 }
304
305
306
307 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
308 void
310 {
311 // Store stream flags
312 const std::ios_base::fmtflags stream_flags(stream.flags());
313 // Set stream to print booleans as "true"/"false"
314 stream.setf(std::ios_base::boolalpha);
315
316 stream << "Active tape index: " << active_tape_index() << '\n';
317 stream << "Recording? " << is_recording() << '\n';
318 stream << std::flush;
319
321 taped_driver.print(stream);
322
323 stream << "Registered independent variables: " << '\n';
324 for (unsigned int i = 0; i < n_independent_variables(); ++i)
325 stream << registered_independent_variable_values[i]
326 << (i < (n_independent_variables() - 1) ? "," : "");
327 stream << std::endl;
328
329 stream << "Independent variable values: " << '\n';
330 print_values(stream);
331
332 stream << "Registered marked independent variables: " << '\n';
333 for (unsigned int i = 0; i < n_independent_variables(); ++i)
334 stream << registered_marked_independent_variables[i]
335 << (i < (n_independent_variables() - 1) ? "," : "")
336 << std::flush;
337 stream << std::endl;
338
339 stream << "Dependent variable values: " << '\n';
340 for (unsigned int i = 0; i < n_dependent_variables(); ++i)
341 stream << dependent_variables[i]
342 << (i < (n_dependent_variables() - 1) ? "," : "");
343 stream << std::endl;
344
345 stream << "Registered dependent variables: " << '\n';
346 for (unsigned int i = 0; i < n_dependent_variables(); ++i)
347 stream << registered_marked_dependent_variables[i]
348 << (i < (n_dependent_variables() - 1) ? "," : "");
349 stream << std::endl;
350
351 // Restore stream flags
352 stream.flags(stream_flags);
353 }
354
355
356
357 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
358 void
360 std::ostream &stream) const
361 {
362 for (unsigned int i = 0; i < n_independent_variables(); ++i)
363 stream << independent_variable_values[i]
364 << (i < (n_independent_variables() - 1) ? "," : "")
365 << std::flush;
366 stream << std::endl;
367 }
368
369
370
371 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
372 void
374 const typename Types<ad_type>::tape_index tape_index,
375 std::ostream & stream) const
376 {
378 return;
379
380 Assert(is_registered_tape(tape_index),
381 ExcMessage("Tape number not registered"));
382
383 this->taped_driver.print_tape_stats(tape_index, stream);
384 }
385
386
387
388 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
389 void
391 const unsigned int n_independent_variables,
392 const unsigned int n_dependent_variables,
393 const bool clear_registered_tapes)
395 const unsigned int new_n_independent_variables =
396 (n_independent_variables != ::numbers::invalid_unsigned_int ?
397 n_independent_variables :
398 this->n_independent_variables());
399 const unsigned int new_n_dependent_variables =
400 (n_dependent_variables != ::numbers::invalid_unsigned_int ?
401 n_dependent_variables :
402 this->n_dependent_variables());
403
404 // Here we clear our vectors of AD data with a reallocation of memory
405 // (i.e. we *nuke* them entirely). Why we do this differs for each
406 // AD type:
407 // - ADOL-C taped mode must have their tapes fully cleared of marked data
408 // before the tapes can be overwritten.
409 // - ADOL-C tapeless mode must be configured for before any active live
410 // variables are created.
411 // - Reverse-mode Sacado numbers to must be destroyed to reset their
412 // accumulations.
413 // - Forward-mode Sacado numbers have no specific requirements, but it
414 // doesn't really hurt to perform this operation anyway.
415 {
416 std::vector<ad_type>().swap(independent_variables);
417 std::vector<ad_type>().swap(dependent_variables);
418 }
419
420 // Tapeless mode must be configured before any active live
421 // variables are created.
423 {
424 configure_tapeless_mode(new_n_independent_variables,
425 false /*ensure_persistent_setting*/);
426 }
428 taped_driver.reset(clear_registered_tapes);
429
430 independent_variable_values = std::vector<scalar_type>(
431 new_n_independent_variables,
433 registered_independent_variable_values =
434 std::vector<bool>(new_n_independent_variables, false);
435 registered_marked_independent_variables =
436 std::vector<bool>(new_n_independent_variables, false);
437 dependent_variables =
438 std::vector<ad_type>(new_n_dependent_variables,
440 registered_marked_dependent_variables =
441 std::vector<bool>(new_n_dependent_variables, false);
442 }
443
444
445
446 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
447 void
449 const unsigned int n_independent_variables,
450 const bool ensure_persistent_setting)
451 {
453 return;
454
455 // Try to safely initialize the global environment
457 n_independent_variables);
458
459 if (ensure_persistent_setting == true)
461 {
462 // In order to ensure that the settings remain for the entire
463 // duration of the simulation, we create a global live variable
464 // that doesn't go out of scope.
465 static ad_type num = 0.0;
466 (void)num;
468 }
469
470
471
472 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
473 void
475 const typename Types<ad_type>::tape_index tape_index)
476 {
477 activate_tape(tape_index, true /*read_mode*/);
478 }
480
481
482 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
483 bool
485 const typename Types<ad_type>::tape_index tape_index) const
486 {
488 return false;
489
490 return taped_driver.requires_retaping(tape_index);
491 }
492
493
494
495 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
496 bool
498 const
499 {
501 return false;
502
503 return taped_driver.last_action_requires_retaping();
504 }
505
506
507
508 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
509 void
511 {
513 return;
514
515 taped_driver.remove_tape(taped_driver.active_tape_index());
516 }
517
518
519
520 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
521 void
523 const typename Types<ad_type>::tape_index tape_index,
524 const bool read_mode)
525 {
527 {
529 ExcMessage("Invalid tape index"));
531 ExcMessage("Tape index exceeds maximum allowable value"));
532 taped_driver.activate_tape(tape_index);
533 reset_registered_independent_variables();
534
535 // A tape may have been defined by a different ADHelper, so in this
536 // case we ignore the fact that any dependent variables within the
537 // current data structure have not been marked as dependents
538 if (read_mode == true)
539 {
540 Assert(is_registered_tape(tape_index),
541 ExcMessage("Tape number not registered"));
542 reset_registered_dependent_variables(true);
543 Assert(n_registered_dependent_variables() ==
544 n_dependent_variables(),
545 ExcMessage("Not all dependent variables have been set!"));
546 }
547 }
548 else
549 {
552 }
553 }
554
555
556
557 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
558 void
560 const typename Types<ad_type>::tape_buffer_sizes obufsize,
561 const typename Types<ad_type>::tape_buffer_sizes lbufsize,
562 const typename Types<ad_type>::tape_buffer_sizes vbufsize,
563 const typename Types<ad_type>::tape_buffer_sizes tbufsize)
564 {
565 // When valid for the chosen AD number type, these values will be used the
566 // next time start_recording_operations() is called.
568 taped_driver.set_tape_buffer_sizes(obufsize,
569 lbufsize,
570 vbufsize,
571 tbufsize);
572 }
573
574
576 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
577 bool
579 const typename Types<ad_type>::tape_index tape_index,
580 const bool overwrite_tape,
581 const bool keep_independent_values)
582 {
583 // Define this here for clarity when this flag is used later.
584 const bool read_mode = false;
585
587 {
588 if (overwrite_tape != true)
589 {
590 Assert(is_recording() == false,
591 ExcMessage("Already recording..."));
592 }
593
594 // Check conditions to enable tracing
595 if (is_registered_tape(tape_index) == false || overwrite_tape == true)
596 {
597 // Setup the data structures for this class in the
598 // appropriate manner
599 activate_tape(tape_index, read_mode);
600
601 // Start taping
602 taped_driver.start_taping(active_tape_index(),
603 keep_independent_values);
604
605 // Clear the flags that state which independent and
606 // dependent variables have been registered
607 reset_registered_independent_variables();
608 reset_registered_dependent_variables();
609 }
610 else
611 {
612 Assert(is_recording() == false,
614 "Tape recording is unexpectedly still enabled."));
615
616 // Now we activate the pre-recorded tape so that its immediately
617 // available for use
618 activate_recorded_tape(tape_index);
619 }
620 }
621 else
622 {
623 Assert(ADNumberTraits<ad_type>::is_tapeless == true,
625
626 // Set the flag that states that we can safely mark dependent
627 // variables within this current phase of operations
628 tapeless_driver.allow_dependent_variable_marking();
629
630 // Dummy call to ensure that the intuitively correct
631 // value for the active tape (whether "valid" or not)
632 // is always returned to the user.
633 activate_tape(tape_index, read_mode);
634 }
635
636 return is_recording();
637 }
638
639
640
641 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
642 void
644 const bool write_tapes_to_file)
645 {
646 Assert(is_recording() == true, ExcMessage("Not currently recording..."));
647
648 // Double check that we've actually registered all DoFs
649 Assert(n_registered_independent_variables() == n_independent_variables(),
650 ExcMessage("Not all values of sensitivities have been recorded!"));
651
653 {
654 // Stop tracing
655 taped_driver.stop_taping(active_tape_index(), write_tapes_to_file);
656 }
657 else
658 {
661 // Double check that we've actually registered dependent variables
662 Assert(n_registered_dependent_variables() == n_dependent_variables(),
663 ExcMessage("Not all dependent variables have been set!"));
664
665 // By changing this flag, we ensure that the we can no longer
666 // legally alter the values of the dependent variables using
667 // set_dependent_variable(). This is important because the value of
668 // the tapeless independent variables are set and finalized when
669 // mark_dependent_variable() is called. So we cannot allow this to
670 // be done when not in the "recording" phase.
671 tapeless_driver.prevent_dependent_variable_marking();
672 }
673 }
675
676
677 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
678 void
680 const unsigned int index,
681 const ad_type & func)
682 {
683 Assert(index < n_dependent_variables(), ExcMessage("Index out of range"));
684 Assert(registered_marked_dependent_variables[index] == false,
686 "This dependent variable has already been registered."));
687
689 {
690 Assert(active_tape_index() != Numbers<ad_type>::invalid_tape_index,
691 ExcMessage("Invalid tape index"));
692 Assert(is_recording() == true,
694 "Must be recording when registering dependent variables."));
696
697 // Register the given dependent variable
698 internal::Marking<ad_type>::dependent_variable(dependent_variables[index],
699 func);
700 registered_marked_dependent_variables[index] = true;
701 }
702
703
704
705 /* -------------------- CellLevelBase -------------------- */
706
707
709 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
711 const unsigned int n_independent_variables,
712 const unsigned int n_dependent_variables)
713 : HelperBase<ADNumberTypeCode, ScalarType>(n_independent_variables,
714 n_dependent_variables)
715 {}
716
717
718
719 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
720 void
722 const std::vector<scalar_type> &dof_values)
723 {
724 // This is actually the same thing the set_independent_variable function,
725 // in the sense that we simply populate our array of independent values
726 // with a meaningful number. However, in this case we need to double check
727 // that we're not registering these variables twice
728 Assert(dof_values.size() == this->n_independent_variables(),
730 "Vector size does not match number of independent variables"));
731 for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
732 {
733 Assert(this->registered_independent_variable_values[i] == false,
734 ExcMessage("Independent variable value already registered."));
735 }
736 set_dof_values(dof_values);
737 }
738
739
740
741 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
742 const std::vector<
745 const
746 {
748 {
749 Assert(this->active_tape_index() !=
751 ExcMessage("Invalid tape index"));
752 }
753
754 // If necessary, initialize the internally stored vector of
755 // AD numbers that represents the independent variables
756 this->finalize_sensitive_independent_variables();
757 Assert(this->independent_variables.size() ==
758 this->n_independent_variables(),
759 ExcDimensionMismatch(this->independent_variables.size(),
760 this->n_independent_variables()));
761
762 return this->independent_variables;
764
765
766
767 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
768 void
770 const std::vector<scalar_type> &values)
771 {
773 {
774 Assert(this->active_tape_index() !=
776 ExcMessage("Invalid tape index"));
777 }
778 Assert(values.size() == this->n_independent_variables(),
780 "Vector size does not match number of independent variables"));
781 for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
783 i, values[i]);
784 }
785
786
787
788 /* ------------------ EnergyFunctional ------------------ */
789
790
791
792 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
794 const unsigned int n_independent_variables)
795 : CellLevelBase<ADNumberTypeCode, ScalarType>(n_independent_variables, 1)
796 {}
797
798
799
800 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
801 void
803 const ad_type &energy)
804 {
805 Assert(this->n_dependent_variables() == 1, ExcInternalError());
807 0, energy);
808 }
809
810
811
812 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
815 {
817 this->taped_driver.keep_independent_values() == false) ||
819 {
820 Assert(
821 this->n_registered_independent_variables() ==
822 this->n_independent_variables(),
824 "Not all values of sensitivities have been registered or subsequently set!"));
825 }
826 Assert(this->n_registered_dependent_variables() ==
827 this->n_dependent_variables(),
828 ExcMessage("Not all dependent variables have been registered."));
829
830 Assert(
831 this->n_dependent_variables() == 1,
833 "The EnergyFunctional class expects there to be only one dependent variable."));
834
836 {
837 Assert(this->active_tape_index() !=
839 ExcMessage("Invalid tape index"));
840 Assert(this->is_recording() == false,
842 "Cannot compute value while tape is being recorded."));
843 Assert(this->independent_variable_values.size() ==
844 this->n_independent_variables(),
845 ExcDimensionMismatch(this->independent_variable_values.size(),
846 this->n_independent_variables()));
847
848 return this->taped_driver.value(this->active_tape_index(),
849 this->independent_variable_values);
850 }
851 else
852 {
855 Assert(this->independent_variables.size() ==
856 this->n_independent_variables(),
857 ExcDimensionMismatch(this->independent_variables.size(),
858 this->n_independent_variables()));
859
860 return this->tapeless_driver.value(this->dependent_variables);
861 }
862 }
863
864
865
866 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
867 void
869 Vector<scalar_type> &gradient) const
870 {
872 this->taped_driver.keep_independent_values() == false) ||
874 {
875 Assert(
876 this->n_registered_independent_variables() ==
877 this->n_independent_variables(),
879 "Not all values of sensitivities have been registered or subsequently set!"));
880 }
881 Assert(this->n_registered_dependent_variables() ==
882 this->n_dependent_variables(),
883 ExcMessage("Not all dependent variables have been registered."));
884
885 Assert(
886 this->n_dependent_variables() == 1,
888 "The EnergyFunctional class expects there to be only one dependent variable."));
889
890 // We can neglect correctly initializing the entries as
891 // we'll be overwriting them immediately in the succeeding call to
892 // Drivers::gradient().
893 if (gradient.size() != this->n_independent_variables())
894 gradient.reinit(this->n_independent_variables(),
895 true /*omit_zeroing_entries*/);
896
898 {
899 Assert(this->active_tape_index() !=
901 ExcMessage("Invalid tape index"));
902 Assert(this->is_recording() == false,
904 "Cannot compute gradient while tape is being recorded."));
905 Assert(this->independent_variable_values.size() ==
906 this->n_independent_variables(),
907 ExcDimensionMismatch(this->independent_variable_values.size(),
908 this->n_independent_variables()));
909
910 this->taped_driver.gradient(this->active_tape_index(),
911 this->independent_variable_values,
912 gradient);
913 }
914 else
915 {
918 Assert(this->independent_variables.size() ==
919 this->n_independent_variables(),
920 ExcDimensionMismatch(this->independent_variables.size(),
921 this->n_independent_variables()));
922
923 this->tapeless_driver.gradient(this->independent_variables,
924 this->dependent_variables,
925 gradient);
926 }
927 }
928
929
930
931 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
932 void
934 FullMatrix<scalar_type> &hessian) const
935 {
938 "Cannot computed function Hessian: AD number type does "
939 "not support the calculation of second order derivatives."));
940
942 this->taped_driver.keep_independent_values() == false))
943 {
944 Assert(
945 this->n_registered_independent_variables() ==
946 this->n_independent_variables(),
948 "Not all values of sensitivities have been registered or subsequently set!"));
949 }
950 Assert(this->n_registered_dependent_variables() ==
951 this->n_dependent_variables(),
952 ExcMessage("Not all dependent variables have been registered."));
953
954 Assert(
955 this->n_dependent_variables() == 1,
957 "The EnergyFunctional class expects there to be only one dependent variable."));
958
959 // We can neglect correctly initializing the entries as
960 // we'll be overwriting them immediately in the succeeding call to
961 // Drivers::hessian().
962 if (hessian.m() != this->n_independent_variables() ||
963 hessian.n() != this->n_independent_variables())
964 hessian.reinit({this->n_independent_variables(),
965 this->n_independent_variables()},
966 true /*omit_default_initialization*/);
967
969 {
970 Assert(this->active_tape_index() !=
972 ExcMessage("Invalid tape index"));
973 Assert(this->is_recording() == false,
975 "Cannot compute hessian while tape is being recorded."));
976 Assert(this->independent_variable_values.size() ==
977 this->n_independent_variables(),
978 ExcDimensionMismatch(this->independent_variable_values.size(),
979 this->n_independent_variables()));
981 this->taped_driver.hessian(this->active_tape_index(),
982 this->independent_variable_values,
983 hessian);
984 }
985 else
986 {
989 Assert(this->independent_variables.size() ==
990 this->n_independent_variables(),
991 ExcDimensionMismatch(this->independent_variables.size(),
992 this->n_independent_variables()));
993
994 this->tapeless_driver.hessian(this->independent_variables,
995 this->dependent_variables,
996 hessian);
997 }
998 }
999
1000
1001 /* ------------------- ResidualLinearization ------------------- */
1002
1003
1004
1005 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
1007 const unsigned int n_independent_variables,
1008 const unsigned int n_dependent_variables)
1009 : CellLevelBase<ADNumberTypeCode, ScalarType>(n_independent_variables,
1010 n_dependent_variables)
1011 {}
1012
1013
1014
1015 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
1016 void
1018 register_residual_vector(const std::vector<ad_type> &residual)
1019 {
1020 Assert(residual.size() == this->n_dependent_variables(),
1021 ExcMessage(
1022 "Vector size does not match number of dependent variables"));
1023 for (unsigned int i = 0; i < this->n_dependent_variables(); ++i)
1025 i, residual[i]);
1026 }
1027
1028
1029
1030 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
1031 void
1033 Vector<scalar_type> &values) const
1034 {
1035 if ((ADNumberTraits<ad_type>::is_taped == true &&
1036 this->taped_driver.keep_independent_values() == false) ||
1038 {
1039 Assert(
1040 this->n_registered_independent_variables() ==
1041 this->n_independent_variables(),
1042 ExcMessage(
1043 "Not all values of sensitivities have been registered or subsequently set!"));
1044 }
1045 Assert(this->n_registered_dependent_variables() ==
1046 this->n_dependent_variables(),
1047 ExcMessage("Not all dependent variables have been registered."));
1048
1049 // We can neglect correctly initializing the entries as
1050 // we'll be overwriting them immediately in the succeeding call to
1051 // Drivers::values().
1052 if (values.size() != this->n_dependent_variables())
1053 values.reinit(this->n_dependent_variables(),
1054 true /*omit_zeroing_entries*/);
1055
1057 {
1058 Assert(this->active_tape_index() !=
1060 ExcMessage("Invalid tape index"));
1061 Assert(this->is_recording() == false,
1062 ExcMessage(
1063 "Cannot compute values while tape is being recorded."));
1064 Assert(this->independent_variable_values.size() ==
1065 this->n_independent_variables(),
1066 ExcDimensionMismatch(this->independent_variable_values.size(),
1067 this->n_independent_variables()));
1068
1069 this->taped_driver.values(this->active_tape_index(),
1070 this->n_dependent_variables(),
1071 this->independent_variable_values,
1072 values);
1073 }
1074 else
1075 {
1078 this->tapeless_driver.values(this->dependent_variables, values);
1079 }
1080 }
1081
1082
1083
1084 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
1085 void
1087 FullMatrix<scalar_type> &jacobian) const
1088 {
1089 if ((ADNumberTraits<ad_type>::is_taped == true &&
1090 this->taped_driver.keep_independent_values() == false) ||
1092 {
1093 Assert(
1094 this->n_registered_independent_variables() ==
1095 this->n_independent_variables(),
1096 ExcMessage(
1097 "Not all values of sensitivities have been registered or subsequently set!"));
1098 }
1099 Assert(this->n_registered_dependent_variables() ==
1100 this->n_dependent_variables(),
1101 ExcMessage("Not all dependent variables have been registered."));
1102
1103 // We can neglect correctly initializing the entries as
1104 // we'll be overwriting them immediately in the succeeding call to
1105 // Drivers::jacobian().
1106 if (jacobian.m() != this->n_dependent_variables() ||
1107 jacobian.n() != this->n_independent_variables())
1108 jacobian.reinit({this->n_dependent_variables(),
1109 this->n_independent_variables()},
1110 true /*omit_default_initialization*/);
1111
1113 {
1114 Assert(this->active_tape_index() !=
1116 ExcMessage("Invalid tape index"));
1117 Assert(this->is_recording() == false,
1118 ExcMessage(
1119 "Cannot compute hessian while tape is being recorded."));
1120 Assert(this->independent_variable_values.size() ==
1121 this->n_independent_variables(),
1122 ExcDimensionMismatch(this->independent_variable_values.size(),
1123 this->n_independent_variables()));
1124
1125 this->taped_driver.jacobian(this->active_tape_index(),
1126 this->n_dependent_variables(),
1127 this->independent_variable_values,
1128 jacobian);
1129 }
1130 else
1131 {
1134 Assert(this->independent_variables.size() ==
1135 this->n_independent_variables(),
1136 ExcDimensionMismatch(this->independent_variables.size(),
1137 this->n_independent_variables()));
1138
1139 this->tapeless_driver.jacobian(this->independent_variables,
1140 this->dependent_variables,
1141 jacobian);
1142 }
1143 }
1144
1145
1146
1147 /* ----------------- PointLevelFunctionsBase ----------------- */
1148
1149
1150
1151 template <int dim,
1152 enum AD::NumberTypes ADNumberTypeCode,
1153 typename ScalarType>
1155 PointLevelFunctionsBase(const unsigned int n_independent_variables,
1156 const unsigned int n_dependent_variables)
1157 : HelperBase<ADNumberTypeCode, ScalarType>(n_independent_variables,
1158 n_dependent_variables)
1159 , symmetric_independent_variables(n_independent_variables, false)
1160 {}
1161
1162
1163
1164 template <int dim,
1165 enum AD::NumberTypes ADNumberTypeCode,
1166 typename ScalarType>
1167 void
1169 const unsigned int n_independent_variables,
1170 const unsigned int n_dependent_variables,
1171 const bool clear_registered_tapes)
1172 {
1174 n_dependent_variables,
1175 clear_registered_tapes);
1176
1177 const unsigned int new_n_independent_variables =
1178 (n_independent_variables != ::numbers::invalid_unsigned_int ?
1179 n_independent_variables :
1180 this->n_independent_variables());
1181 symmetric_independent_variables =
1182 std::vector<bool>(new_n_independent_variables, false);
1183 }
1184
1185
1186
1187 template <int dim,
1188 enum AD::NumberTypes ADNumberTypeCode,
1189 typename ScalarType>
1190 bool
1192 is_symmetric_independent_variable(const unsigned int index) const
1193 {
1194 Assert(index < symmetric_independent_variables.size(),
1196 return symmetric_independent_variables[index];
1197 }
1198
1199
1200
1201 template <int dim,
1202 enum AD::NumberTypes ADNumberTypeCode,
1203 typename ScalarType>
1204 unsigned int
1207 {
1208 return std::count(symmetric_independent_variables.begin(),
1209 symmetric_independent_variables.end(),
1210 true);
1211 }
1212
1213
1214
1215 template <int dim,
1216 enum AD::NumberTypes ADNumberTypeCode,
1217 typename ScalarType>
1218 void
1220 register_independent_variables(const std::vector<scalar_type> &values)
1221 {
1222 // This is actually the same thing the set_independent_variable function,
1223 // in the sense that we simply populate our array of independent values
1224 // with a meaningful number. However, in this case we need to double check
1225 // that we're not registering these variables twice
1226 Assert(values.size() == this->n_independent_variables(),
1227 ExcMessage(
1228 "Vector size does not match number of independent variables"));
1229 for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
1230 {
1231 Assert(this->registered_independent_variable_values[i] == false,
1232 ExcMessage("Independent variable value already registered."));
1233 }
1234 set_independent_variables(values);
1235 }
1236
1237
1238
1239 template <int dim,
1240 enum AD::NumberTypes ADNumberTypeCode,
1241 typename ScalarType>
1242 const std::vector<typename PointLevelFunctionsBase<dim,
1243 ADNumberTypeCode,
1244 ScalarType>::ad_type> &
1247 {
1249 {
1250 Assert(this->active_tape_index() !=
1252 ExcMessage("Invalid tape index"));
1253 }
1254
1255 // Just in case the user has not done so, we repeat the call to
1256 // initialize the internally stored vector of AD numbers that
1257 // represents the independent variables.
1258 this->finalize_sensitive_independent_variables();
1259 Assert(this->independent_variables.size() ==
1260 this->n_independent_variables(),
1261 ExcDimensionMismatch(this->independent_variables.size(),
1262 this->n_independent_variables()));
1263
1264 return this->independent_variables;
1265 }
1266
1267
1268
1269 template <int dim,
1270 enum AD::NumberTypes ADNumberTypeCode,
1271 typename ScalarType>
1272 void
1274 set_sensitivity_value(const unsigned int index,
1275 const bool symmetric_component,
1276 const scalar_type &value)
1277 {
1279 value);
1280 Assert(
1281 index < this->n_independent_variables(),
1282 ExcMessage(
1283 "Trying to set the symmetry flag of a non-existent independent variable."));
1284 Assert(index < symmetric_independent_variables.size(),
1286 symmetric_independent_variables[index] = symmetric_component;
1287 }
1288
1289
1290
1291 template <int dim,
1292 enum AD::NumberTypes ADNumberTypeCode,
1293 typename ScalarType>
1294 void
1296 set_independent_variables(const std::vector<scalar_type> &values)
1297 {
1299 {
1300 Assert(this->active_tape_index() !=
1302 ExcMessage("Invalid tape index"));
1303 }
1304 Assert(values.size() == this->n_independent_variables(),
1305 ExcMessage(
1306 "Vector size does not match number of independent variables"));
1307 for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
1309 i, values[i]);
1310 }
1311
1312
1313
1314 /* -------------------- ScalarFunction -------------------- */
1315
1316
1317
1318 template <int dim,
1319 enum AD::NumberTypes ADNumberTypeCode,
1320 typename ScalarType>
1322 const unsigned int n_independent_variables)
1323 : PointLevelFunctionsBase<dim, ADNumberTypeCode, ScalarType>(
1324 n_independent_variables,
1325 1)
1326 {}
1327
1328
1329
1330 template <int dim,
1331 enum AD::NumberTypes ADNumberTypeCode,
1332 typename ScalarType>
1333 void
1336 {
1337 Assert(this->n_dependent_variables() == 1, ExcInternalError());
1339 0, func);
1340 }
1341
1342
1343
1344 template <int dim,
1345 enum AD::NumberTypes ADNumberTypeCode,
1346 typename ScalarType>
1349 {
1350 if ((ADNumberTraits<ad_type>::is_taped == true &&
1351 this->taped_driver.keep_independent_values() == false) ||
1353 {
1354 Assert(
1355 this->n_registered_independent_variables() ==
1356 this->n_independent_variables(),
1357 ExcMessage(
1358 "Not all values of sensitivities have been registered or subsequently set!"));
1359 }
1360 Assert(this->n_registered_dependent_variables() ==
1361 this->n_dependent_variables(),
1362 ExcMessage("Not all dependent variables have been registered."));
1363
1364 Assert(
1365 this->n_dependent_variables() == 1,
1366 ExcMessage(
1367 "The ScalarFunction class expects there to be only one dependent variable."));
1368
1370 {
1371 Assert(this->active_tape_index() !=
1373 ExcMessage("Invalid tape index"));
1374 Assert(this->is_recording() == false,
1375 ExcMessage(
1376 "Cannot compute values while tape is being recorded."));
1377 Assert(this->independent_variable_values.size() ==
1378 this->n_independent_variables(),
1379 ExcDimensionMismatch(this->independent_variable_values.size(),
1380 this->n_independent_variables()));
1381
1382 return this->taped_driver.value(this->active_tape_index(),
1383 this->independent_variable_values);
1384 }
1385 else
1386 {
1389 return this->tapeless_driver.value(this->dependent_variables);
1390 }
1391 }
1392
1393
1394 template <int dim,
1395 enum AD::NumberTypes ADNumberTypeCode,
1396 typename ScalarType>
1397 void
1399 Vector<scalar_type> &gradient) const
1400 {
1401 if ((ADNumberTraits<ad_type>::is_taped == true &&
1402 this->taped_driver.keep_independent_values() == false) ||
1404 {
1405 Assert(
1406 this->n_registered_independent_variables() ==
1407 this->n_independent_variables(),
1408 ExcMessage(
1409 "Not all values of sensitivities have been registered or subsequently set!"));
1410 }
1411 Assert(this->n_registered_dependent_variables() ==
1412 this->n_dependent_variables(),
1413 ExcMessage("Not all dependent variables have been registered."));
1414
1415 Assert(
1416 this->n_dependent_variables() == 1,
1417 ExcMessage(
1418 "The ScalarFunction class expects there to be only one dependent variable."));
1419
1420 // We can neglect correctly initializing the entries as
1421 // we'll be overwriting them immediately in the succeeding call to
1422 // Drivers::gradient().
1423 if (gradient.size() != this->n_independent_variables())
1424 gradient.reinit(this->n_independent_variables(),
1425 true /*omit_zeroing_entries*/);
1426
1428 {
1429 Assert(this->active_tape_index() !=
1431 ExcMessage("Invalid tape index"));
1432 Assert(this->is_recording() == false,
1433 ExcMessage(
1434 "Cannot compute gradient while tape is being recorded."));
1435 Assert(this->independent_variable_values.size() ==
1436 this->n_independent_variables(),
1437 ExcDimensionMismatch(this->independent_variable_values.size(),
1438 this->n_independent_variables()));
1439
1440 this->taped_driver.gradient(this->active_tape_index(),
1441 this->independent_variable_values,
1442 gradient);
1443 }
1444 else
1445 {
1448 Assert(this->independent_variables.size() ==
1449 this->n_independent_variables(),
1450 ExcDimensionMismatch(this->independent_variables.size(),
1451 this->n_independent_variables()));
1452
1453 this->tapeless_driver.gradient(this->independent_variables,
1454 this->dependent_variables,
1455 gradient);
1456 }
1457
1458 // Account for symmetries of tensor components
1459 for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
1460 {
1461 if (this->is_symmetric_independent_variable(i) == true)
1462 gradient[i] *= 0.5;
1463 }
1464 }
1465
1466
1467
1468 template <int dim,
1469 enum AD::NumberTypes ADNumberTypeCode,
1470 typename ScalarType>
1471 void
1473 FullMatrix<scalar_type> &hessian) const
1474 {
1476 ExcMessage(
1477 "Cannot computed function Hessian: AD number type does "
1478 "not support the calculation of second order derivatives."));
1479
1480 if ((ADNumberTraits<ad_type>::is_taped == true &&
1481 this->taped_driver.keep_independent_values() == false))
1482 {
1483 Assert(
1484 this->n_registered_independent_variables() ==
1485 this->n_independent_variables(),
1486 ExcMessage(
1487 "Not all values of sensitivities have been registered or subsequently set!"));
1488 }
1489 Assert(this->n_registered_dependent_variables() ==
1490 this->n_dependent_variables(),
1491 ExcMessage("Not all dependent variables have been registered."));
1492
1493 Assert(
1494 this->n_dependent_variables() == 1,
1495 ExcMessage(
1496 "The ScalarFunction class expects there to be only one dependent variable."));
1497
1498 // We can neglect correctly initializing the entries as
1499 // we'll be overwriting them immediately in the succeeding call to
1500 // Drivers::hessian().
1501 if (hessian.m() != this->n_independent_variables() ||
1502 hessian.n() != this->n_independent_variables())
1503 hessian.reinit({this->n_independent_variables(),
1504 this->n_independent_variables()},
1505 true /*omit_default_initialization*/);
1506
1508 {
1509 Assert(this->active_tape_index() !=
1511 ExcMessage("Invalid tape index"));
1512 Assert(this->is_recording() == false,
1513 ExcMessage(
1514 "Cannot compute Hessian while tape is being recorded."));
1515 Assert(this->independent_variable_values.size() ==
1516 this->n_independent_variables(),
1517 ExcDimensionMismatch(this->independent_variable_values.size(),
1518 this->n_independent_variables()));
1519
1520 this->taped_driver.hessian(this->active_tape_index(),
1521 this->independent_variable_values,
1522 hessian);
1523 }
1524 else
1525 {
1528 Assert(this->independent_variables.size() ==
1529 this->n_independent_variables(),
1530 ExcDimensionMismatch(this->independent_variables.size(),
1531 this->n_independent_variables()));
1532
1533 this->tapeless_driver.hessian(this->independent_variables,
1534 this->dependent_variables,
1535 hessian);
1536 }
1537
1538 // Account for symmetries of tensor components
1539 for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
1540 for (unsigned int j = 0; j < i + 1; ++j)
1541 {
1542 if (this->is_symmetric_independent_variable(i) == true &&
1543 this->is_symmetric_independent_variable(j) == true)
1544 {
1545 hessian[i][j] *= 0.25;
1546 if (i != j)
1547 hessian[j][i] *= 0.25;
1548 }
1549 else if ((this->is_symmetric_independent_variable(i) == true &&
1550 this->is_symmetric_independent_variable(j) == false) ||
1551 (this->is_symmetric_independent_variable(j) == true &&
1552 this->is_symmetric_independent_variable(i) == false))
1553 {
1554 hessian[i][j] *= 0.5;
1555 if (i != j)
1556 hessian[j][i] *= 0.5;
1557 }
1558 }
1559 }
1560
1561
1562
1563 template <int dim,
1564 enum AD::NumberTypes ADNumberTypeCode,
1565 typename ScalarType>
1566 Tensor<
1567 0,
1568 dim,
1572 const FEValuesExtractors::Scalar &extractor_row,
1573 const FEValuesExtractors::Scalar &extractor_col)
1574 {
1575 // NOTE: It is necessary to make special provision for the case when the
1576 // HessianType is scalar. Unfortunately Tensor<0,dim> does not provide
1577 // the function unrolled_to_component_indices!
1578 // NOTE: The order of components must be consistently defined throughout
1579 // this class.
1581
1582 // Get indexsets for the subblocks from which we wish to extract the
1583 // matrix values
1584 const std::vector<unsigned int> row_index_set(
1585 internal::extract_field_component_indices<dim>(extractor_row));
1586 const std::vector<unsigned int> col_index_set(
1587 internal::extract_field_component_indices<dim>(extractor_col));
1588 Assert(row_index_set.size() == 1, ExcInternalError());
1589 Assert(col_index_set.size() == 1, ExcInternalError());
1590
1592 0,
1593 hessian[row_index_set[0]][col_index_set[0]]);
1594
1595 return out;
1596 }
1597
1598
1599
1600 template <int dim,
1601 enum AD::NumberTypes ADNumberTypeCode,
1602 typename ScalarType>
1604 4,
1605 dim,
1609 const FullMatrix<scalar_type> & hessian,
1610 const FEValuesExtractors::SymmetricTensor<2> &extractor_row,
1611 const FEValuesExtractors::SymmetricTensor<2> &extractor_col)
1612 {
1613 // NOTE: The order of components must be consistently defined throughout
1614 // this class.
1615 // NOTE: We require a specialisation for rank-4 symmetric tensors because
1616 // they do not define their rank, and setting data using TableIndices is
1617 // somewhat specialised as well.
1619
1620 // Get indexsets for the subblocks from which we wish to extract the
1621 // matrix values
1622 const std::vector<unsigned int> row_index_set(
1623 internal::extract_field_component_indices<dim>(extractor_row));
1624 const std::vector<unsigned int> col_index_set(
1625 internal::extract_field_component_indices<dim>(extractor_col));
1626
1627 for (unsigned int r = 0; r < row_index_set.size(); ++r)
1628 for (unsigned int c = 0; c < col_index_set.size(); ++c)
1629 {
1631 out, r, c, hessian[row_index_set[r]][col_index_set[c]]);
1632 }
1633
1634 return out;
1635 }
1636
1637
1638
1639 /* -------------------- VectorFunction -------------------- */
1640
1641
1642
1643 template <int dim,
1644 enum AD::NumberTypes ADNumberTypeCode,
1645 typename ScalarType>
1647 const unsigned int n_independent_variables,
1648 const unsigned int n_dependent_variables)
1649 : PointLevelFunctionsBase<dim, ADNumberTypeCode, ScalarType>(
1650 n_independent_variables,
1651 n_dependent_variables)
1652 {}
1653
1654
1655
1656 template <int dim,
1657 enum AD::NumberTypes ADNumberTypeCode,
1658 typename ScalarType>
1659 void
1661 register_dependent_variables(const std::vector<ad_type> &funcs)
1662 {
1663 Assert(funcs.size() == this->n_dependent_variables(),
1664 ExcMessage(
1665 "Vector size does not match number of dependent variables"));
1666 for (unsigned int i = 0; i < this->n_dependent_variables(); ++i)
1668 i, funcs[i]);
1669 }
1670
1671
1672
1673 template <int dim,
1674 enum AD::NumberTypes ADNumberTypeCode,
1675 typename ScalarType>
1676 void
1678 Vector<scalar_type> &values) const
1679 {
1680 if ((ADNumberTraits<ad_type>::is_taped == true &&
1681 this->taped_driver.keep_independent_values() == false) ||
1683 {
1684 Assert(
1685 this->n_registered_independent_variables() ==
1686 this->n_independent_variables(),
1687 ExcMessage(
1688 "Not all values of sensitivities have been registered or subsequently set!"));
1689 }
1690 Assert(this->n_registered_dependent_variables() ==
1691 this->n_dependent_variables(),
1692 ExcMessage("Not all dependent variables have been registered."));
1693
1694 // We can neglect correctly initializing the entries as
1695 // we'll be overwriting them immediately in the succeeding call to
1696 // Drivers::values().
1697 if (values.size() != this->n_dependent_variables())
1698 values.reinit(this->n_dependent_variables(),
1699 true /*omit_zeroing_entries*/);
1700
1702 {
1703 Assert(this->active_tape_index() !=
1705 ExcMessage("Invalid tape index"));
1706 Assert(this->is_recording() == false,
1707 ExcMessage(
1708 "Cannot compute values while tape is being recorded."));
1709 Assert(this->independent_variable_values.size() ==
1710 this->n_independent_variables(),
1711 ExcDimensionMismatch(this->independent_variable_values.size(),
1712 this->n_independent_variables()));
1713
1714 this->taped_driver.values(this->active_tape_index(),
1715 this->n_dependent_variables(),
1716 this->independent_variable_values,
1717 values);
1718 }
1719 else
1720 {
1723 this->tapeless_driver.values(this->dependent_variables, values);
1724 }
1725 }
1726
1727
1728
1729 template <int dim,
1730 enum AD::NumberTypes ADNumberTypeCode,
1731 typename ScalarType>
1732 void
1734 FullMatrix<scalar_type> &jacobian) const
1735 {
1736 if ((ADNumberTraits<ad_type>::is_taped == true &&
1737 this->taped_driver.keep_independent_values() == false) ||
1739 {
1740 Assert(
1741 this->n_registered_independent_variables() ==
1742 this->n_independent_variables(),
1743 ExcMessage(
1744 "Not all values of sensitivities have been registered or subsequently set!"));
1745 }
1746 Assert(this->n_registered_dependent_variables() ==
1747 this->n_dependent_variables(),
1748 ExcMessage("Not all dependent variables have been registered."));
1749
1750 // We can neglect correctly initializing the entries as
1751 // we'll be overwriting them immediately in the succeeding call to
1752 // Drivers::jacobian().
1753 if (jacobian.m() != this->n_dependent_variables() ||
1754 jacobian.n() != this->n_independent_variables())
1755 jacobian.reinit({this->n_dependent_variables(),
1756 this->n_independent_variables()},
1757 true /*omit_default_initialization*/);
1758
1760 {
1761 Assert(this->active_tape_index() !=
1763 ExcMessage("Invalid tape index"));
1764 Assert(this->is_recording() == false,
1765 ExcMessage(
1766 "Cannot compute Jacobian while tape is being recorded."));
1767 Assert(this->independent_variable_values.size() ==
1768 this->n_independent_variables(),
1769 ExcDimensionMismatch(this->independent_variable_values.size(),
1770 this->n_independent_variables()));
1771
1772 this->taped_driver.jacobian(this->active_tape_index(),
1773 this->n_dependent_variables(),
1774 this->independent_variable_values,
1775 jacobian);
1776 }
1777 else
1778 {
1781 Assert(this->independent_variables.size() ==
1782 this->n_independent_variables(),
1783 ExcDimensionMismatch(this->independent_variables.size(),
1784 this->n_independent_variables()));
1785
1786 this->tapeless_driver.jacobian(this->independent_variables,
1787 this->dependent_variables,
1788 jacobian);
1789 }
1790
1791 for (unsigned int j = 0; j < this->n_independent_variables(); ++j)
1792 {
1793 // Because we perform just a single differentiation
1794 // operation with respect to the "column" variables,
1795 // we only need to consider them for symmetry conditions.
1796 if (this->is_symmetric_independent_variable(j) == true)
1797 for (unsigned int i = 0; i < this->n_dependent_variables(); ++i)
1798 jacobian[i][j] *= 0.5;
1799 }
1800 }
1801
1802
1803
1804 template <int dim,
1805 enum AD::NumberTypes ADNumberTypeCode,
1806 typename ScalarType>
1807 Tensor<
1808 0,
1809 dim,
1813 const FullMatrix<scalar_type> & jacobian,
1814 const FEValuesExtractors::Scalar &extractor_row,
1815 const FEValuesExtractors::Scalar &extractor_col)
1816 {
1817 // NOTE: It is necessary to make special provision for the case when the
1818 // HessianType is scalar. Unfortunately Tensor<0,dim> does not provide
1819 // the function unrolled_to_component_indices!
1820 // NOTE: The order of components must be consistently defined throughout
1821 // this class.
1823
1824 // Get indexsets for the subblocks from which we wish to extract the
1825 // matrix values
1826 const std::vector<unsigned int> row_index_set(
1827 internal::extract_field_component_indices<dim>(extractor_row));
1828 const std::vector<unsigned int> col_index_set(
1829 internal::extract_field_component_indices<dim>(extractor_col));
1830 Assert(row_index_set.size() == 1, ExcInternalError());
1831 Assert(col_index_set.size() == 1, ExcInternalError());
1832
1834 0,
1835 jacobian[row_index_set[0]][col_index_set[0]]);
1836
1837 return out;
1838 }
1839
1840
1841
1842 template <int dim,
1843 enum AD::NumberTypes ADNumberTypeCode,
1844 typename ScalarType>
1846 4,
1847 dim,
1851 const FullMatrix<scalar_type> & jacobian,
1852 const FEValuesExtractors::SymmetricTensor<2> &extractor_row,
1853 const FEValuesExtractors::SymmetricTensor<2> &extractor_col)
1854 {
1855 // NOTE: The order of components must be consistently defined throughout
1856 // this class.
1857 // NOTE: We require a specialisation for rank-4 symmetric tensors because
1858 // they do not define their rank, and setting data using TableIndices is
1859 // somewhat specialised as well.
1861
1862 // Get indexsets for the subblocks from which we wish to extract the
1863 // matrix values
1864 const std::vector<unsigned int> row_index_set(
1865 internal::extract_field_component_indices<dim>(extractor_row));
1866 const std::vector<unsigned int> col_index_set(
1867 internal::extract_field_component_indices<dim>(extractor_col));
1868
1869 for (unsigned int r = 0; r < row_index_set.size(); ++r)
1870 for (unsigned int c = 0; c < col_index_set.size(); ++c)
1871 {
1873 out, r, c, jacobian[row_index_set[r]][col_index_set[c]]);
1874 }
1875
1876 return out;
1877 }
1878
1879
1880 } // namespace AD
1881} // namespace Differentiation
1882
1883
1884/* --- Explicit instantiations --- */
1885# include "ad_helpers.inst"
1886
1887# ifdef DEAL_II_WITH_ADOLC
1888# include "ad_helpers.inst1"
1889# endif
1890# ifdef DEAL_II_TRILINOS_WITH_SACADO
1891# include "ad_helpers.inst2"
1892# endif
1893
1894
1896
1897#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:850
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:177
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:94
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:38
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:749
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:184
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#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:213
static void initialize_global_environment(const unsigned int n_independent_variables)