Loading [MathJax]/extensions/TeX/AMSsymbols.js
 Reference documentation for deal.II version 9.3.3
\(\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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
ad_helpers.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2016 - 2020 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
225
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();
245
246
247
248 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
249 unsigned int
251 const
252 {
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();
265 }
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: "
324 << "\n";
325 for (unsigned int i = 0; i < n_independent_variables(); i++)
326 stream << registered_independent_variable_values[i]
327 << (i < (n_independent_variables() - 1) ? "," : "");
328 stream << std::endl;
329
330 stream << "Independent variable values: "
331 << "\n";
332 print_values(stream);
333
334 stream << "Registered marked independent variables: "
335 << "\n";
336 for (unsigned int i = 0; i < n_independent_variables(); i++)
337 stream << registered_marked_independent_variables[i]
338 << (i < (n_independent_variables() - 1) ? "," : "")
339 << std::flush;
340 stream << std::endl;
342 stream << "Dependent variable values: "
343 << "\n";
344 for (unsigned int i = 0; i < n_dependent_variables(); i++)
345 stream << dependent_variables[i]
346 << (i < (n_dependent_variables() - 1) ? "," : "");
347 stream << std::endl;
348
349 stream << "Registered dependent variables: "
350 << "\n";
351 for (unsigned int i = 0; i < n_dependent_variables(); i++)
352 stream << registered_marked_dependent_variables[i]
353 << (i < (n_dependent_variables() - 1) ? "," : "");
354 stream << std::endl;
355
356 // Restore stream flags
357 stream.flags(stream_flags);
358 }
360
361
362 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
363 void
365 std::ostream &stream) const
367 for (unsigned int i = 0; i < n_independent_variables(); i++)
368 stream << independent_variable_values[i]
369 << (i < (n_independent_variables() - 1) ? "," : "")
370 << std::flush;
371 stream << std::endl;
372 }
373
374
375
376 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
377 void
379 const typename Types<ad_type>::tape_index tape_index,
380 std::ostream & stream) const
381 {
383 return;
384
385 Assert(is_registered_tape(tape_index),
386 ExcMessage("Tape number not registered"));
387
388 this->taped_driver.print_tape_stats(tape_index, stream);
389 }
390
391
392
393 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
394 void
396 const unsigned int n_independent_variables,
397 const unsigned int n_dependent_variables,
398 const bool clear_registered_tapes)
399 {
400 const unsigned int new_n_independent_variables =
401 (n_independent_variables != ::numbers::invalid_unsigned_int ?
402 n_independent_variables :
403 this->n_independent_variables());
404 const unsigned int new_n_dependent_variables =
405 (n_dependent_variables != ::numbers::invalid_unsigned_int ?
406 n_dependent_variables :
407 this->n_dependent_variables());
408
409 // Here we clear our vectors of AD data with a reallocation of memory
410 // (i.e. we *nuke* them entirely). Why we do this differs for each
411 // AD type:
412 // - ADOL-C taped mode must have their tapes fully cleared of marked data
413 // before the tapes can be overwritten.
414 // - ADOL-C tapeless mode must be configured for before any active live
415 // variables are created.
416 // - Reverse-mode Sacado numbers to must be destroyed to reset their
417 // accumulations.
418 // - Forward-mode Sacado numbers have no specific requirements, but it
419 // doesn't really hurt to perform this operation anyway.
420 {
421 std::vector<ad_type>().swap(independent_variables);
422 std::vector<ad_type>().swap(dependent_variables);
423 }
424
425 // Tapeless mode must be configured before any active live
426 // variables are created.
428 {
429 configure_tapeless_mode(new_n_independent_variables,
430 false /*ensure_persistent_setting*/);
431 }
433 taped_driver.reset(clear_registered_tapes);
434
435 independent_variable_values = std::vector<scalar_type>(
436 new_n_independent_variables,
438 registered_independent_variable_values =
439 std::vector<bool>(new_n_independent_variables, false);
440 registered_marked_independent_variables =
441 std::vector<bool>(new_n_independent_variables, false);
442 dependent_variables =
443 std::vector<ad_type>(new_n_dependent_variables,
445 registered_marked_dependent_variables =
446 std::vector<bool>(new_n_dependent_variables, false);
447 }
448
449
450
451 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
452 void
454 const unsigned int n_independent_variables,
455 const bool ensure_persistent_setting)
456 {
458 return;
459
460 // Try to safely initialize the global environment
462 n_independent_variables);
463
464 if (ensure_persistent_setting == true)
466 {
467 // In order to ensure that the settings remain for the entire
468 // duration of the simulation, we create a global live variable
469 // that doesn't go out of scope.
470 static ad_type num = 0.0;
471 (void)num;
472 }
473 }
474
475
476
477 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
478 void
480 const typename Types<ad_type>::tape_index tape_index)
481 {
482 activate_tape(tape_index, true /*read_mode*/);
483 }
484
485
486
487 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
488 bool
490 const typename Types<ad_type>::tape_index tape_index) const
491 {
493 return false;
494
495 return taped_driver.requires_retaping(tape_index);
496 }
497
498
499
500 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
501 bool
503 const
504 {
506 return false;
507
508 return taped_driver.last_action_requires_retaping();
509 }
510
511
512
513 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
514 void
516 {
518 return;
519
520 taped_driver.remove_tape(taped_driver.active_tape_index());
521 }
522
524
525 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
526 void
528 const typename Types<ad_type>::tape_index tape_index,
529 const bool read_mode)
530 {
532 {
534 ExcMessage("Invalid tape index"));
536 ExcMessage("Tape index exceeds maximum allowable value"));
537 taped_driver.activate_tape(tape_index);
538 reset_registered_independent_variables();
539
540 // A tape may have been defined by a different ADHelper, so in this
541 // case we ignore the fact that any dependent variables within the
542 // current data structure have not been marked as dependents
543 if (read_mode == true)
544 {
545 Assert(is_registered_tape(tape_index),
546 ExcMessage("Tape number not registered"));
547 reset_registered_dependent_variables(true);
548 Assert(n_registered_dependent_variables() ==
549 n_dependent_variables(),
550 ExcMessage("Not all dependent variables have been set!"));
551 }
552 }
553 else
554 {
557 }
558 }
559
560
561
562 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
563 void
565 const typename Types<ad_type>::tape_buffer_sizes obufsize,
566 const typename Types<ad_type>::tape_buffer_sizes lbufsize,
567 const typename Types<ad_type>::tape_buffer_sizes vbufsize,
568 const typename Types<ad_type>::tape_buffer_sizes tbufsize)
569 {
570 // When valid for the chosen AD number type, these values will be used the
571 // next time start_recording_operations() is called.
573 taped_driver.set_tape_buffer_sizes(obufsize,
574 lbufsize,
575 vbufsize,
576 tbufsize);
577 }
578
579
580
581 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
582 bool
584 const typename Types<ad_type>::tape_index tape_index,
585 const bool overwrite_tape,
586 const bool keep_independent_values)
587 {
588 // Define this here for clarity when this flag is used later.
589 const bool read_mode = false;
590
592 {
593 if (overwrite_tape != true)
594 {
595 Assert(is_recording() == false,
596 ExcMessage("Already recording..."));
597 }
598
599 // Check conditions to enable tracing
600 if (is_registered_tape(tape_index) == false || overwrite_tape == true)
601 {
602 // Setup the data structures for this class in the
603 // appropriate manner
604 activate_tape(tape_index, read_mode);
605
606 // Start taping
607 taped_driver.start_taping(active_tape_index(),
608 keep_independent_values);
609
610 // Clear the flags that state which independent and
611 // dependent variables have been registered
612 reset_registered_independent_variables();
613 reset_registered_dependent_variables();
614 }
615 else
616 {
617 Assert(is_recording() == false,
619 "Tape recording is unexpectedly still enabled."));
620
621 // Now we activate the pre-recorded tape so that its immediately
622 // available for use
623 activate_recorded_tape(tape_index);
624 }
625 }
626 else
627 {
628 Assert(ADNumberTraits<ad_type>::is_tapeless == true,
630
631 // Set the flag that states that we can safely mark dependent
632 // variables within this current phase of operations
633 tapeless_driver.allow_dependent_variable_marking();
634
635 // Dummy call to ensure that the intuitively correct
636 // value for the active tape (whether "valid" or not)
637 // is always returned to the user.
638 activate_tape(tape_index, read_mode);
639 }
640
641 return is_recording();
642 }
643
644
645
646 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
647 void
649 const bool write_tapes_to_file)
650 {
651 Assert(is_recording() == true, ExcMessage("Not currently recording..."));
652
653 // Double check that we've actually registered all DoFs
654 Assert(n_registered_independent_variables() == n_independent_variables(),
655 ExcMessage("Not all values of sensitivities have been recorded!"));
656
658 {
659 // Stop tracing
660 taped_driver.stop_taping(active_tape_index(), write_tapes_to_file);
661 }
662 else
663 {
666 // Double check that we've actually registered dependent variables
667 Assert(n_registered_dependent_variables() == n_dependent_variables(),
668 ExcMessage("Not all dependent variables have been set!"));
669
670 // By changing this flag, we ensure that the we can no longer
671 // legally alter the values of the dependent variables using
672 // set_dependent_variable(). This is important because the value of
673 // the tapeless independent variables are set and finalized when
674 // mark_dependent_variable() is called. So we cannot allow this to
675 // be done when not in the "recording" phase.
676 tapeless_driver.prevent_dependent_variable_marking();
677 }
678 }
679
680
681
682 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
683 void
685 const unsigned int index,
686 const ad_type & func)
687 {
688 Assert(index < n_dependent_variables(), ExcMessage("Index out of range"));
689 Assert(registered_marked_dependent_variables[index] == false,
691 "This dependent variable has already been registered."));
692
694 {
695 Assert(active_tape_index() != Numbers<ad_type>::invalid_tape_index,
696 ExcMessage("Invalid tape index"));
697 Assert(is_recording() == true,
699 "Must be recording when registering dependent variables."));
700 }
701
702 // Register the given dependent variable
703 internal::Marking<ad_type>::dependent_variable(dependent_variables[index],
704 func);
705 registered_marked_dependent_variables[index] = true;
706 }
707
708
710 /* -------------------- CellLevelBase -------------------- */
711
712
713
714 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
716 const unsigned int n_independent_variables,
717 const unsigned int n_dependent_variables)
718 : HelperBase<ADNumberTypeCode, ScalarType>(n_independent_variables,
719 n_dependent_variables)
720 {}
722
723
724 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
725 void
727 const std::vector<scalar_type> &dof_values)
728 {
729 // This is actually the same thing the set_independent_variable function,
730 // in the sense that we simply populate our array of independent values
731 // with a meaningful number. However, in this case we need to double check
732 // that we're not registering these variables twice
733 Assert(dof_values.size() == this->n_independent_variables(),
735 "Vector size does not match number of independent variables"));
736 for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
737 {
738 Assert(this->registered_independent_variable_values[i] == false,
739 ExcMessage("Independent variable value already registered."));
740 }
741 set_dof_values(dof_values);
742 }
743
744
745
746 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
747 const std::vector<
750 const
751 {
753 {
754 Assert(this->active_tape_index() !=
756 ExcMessage("Invalid tape index"));
757 }
758
759 // If necessary, initialize the internally stored vector of
760 // AD numbers that represents the independent variables
761 this->finalize_sensitive_independent_variables();
762 Assert(this->independent_variables.size() ==
763 this->n_independent_variables(),
764 ExcDimensionMismatch(this->independent_variables.size(),
765 this->n_independent_variables()));
766
767 return this->independent_variables;
768 }
769
771
772 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
773 void
775 const std::vector<scalar_type> &values)
776 {
778 {
779 Assert(this->active_tape_index() !=
781 ExcMessage("Invalid tape index"));
782 }
783 Assert(values.size() == this->n_independent_variables(),
785 "Vector size does not match number of independent variables"));
786 for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
788 i, values[i]);
789 }
790
791
792
793 /* ------------------ EnergyFunctional ------------------ */
794
795
796
797 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
799 const unsigned int n_independent_variables)
800 : CellLevelBase<ADNumberTypeCode, ScalarType>(n_independent_variables, 1)
801 {}
802
803
804
805 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
806 void
808 const ad_type &energy)
809 {
810 Assert(this->n_dependent_variables() == 1, ExcInternalError());
812 0, energy);
813 }
814
815
816
817 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
820 {
822 this->taped_driver.keep_independent_values() == false) ||
824 {
825 Assert(
826 this->n_registered_independent_variables() ==
827 this->n_independent_variables(),
829 "Not all values of sensitivities have been registered or subsequently set!"));
830 }
831 Assert(this->n_registered_dependent_variables() ==
832 this->n_dependent_variables(),
833 ExcMessage("Not all dependent variables have been registered."));
834
835 Assert(
836 this->n_dependent_variables() == 1,
838 "The EnergyFunctional class expects there to be only one dependent variable."));
839
841 {
842 Assert(this->active_tape_index() !=
844 ExcMessage("Invalid tape index"));
845 Assert(this->is_recording() == false,
847 "Cannot compute value while tape is being recorded."));
848 Assert(this->independent_variable_values.size() ==
849 this->n_independent_variables(),
850 ExcDimensionMismatch(this->independent_variable_values.size(),
851 this->n_independent_variables()));
852
853 return this->taped_driver.value(this->active_tape_index(),
854 this->independent_variable_values);
855 }
856 else
857 {
860 Assert(this->independent_variables.size() ==
861 this->n_independent_variables(),
862 ExcDimensionMismatch(this->independent_variables.size(),
863 this->n_independent_variables()));
864
865 return this->tapeless_driver.value(this->dependent_variables);
866 }
867 }
868
869
870
871 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
872 void
874 Vector<scalar_type> &gradient) const
875 {
877 this->taped_driver.keep_independent_values() == false) ||
879 {
880 Assert(
881 this->n_registered_independent_variables() ==
882 this->n_independent_variables(),
884 "Not all values of sensitivities have been registered or subsequently set!"));
885 }
886 Assert(this->n_registered_dependent_variables() ==
887 this->n_dependent_variables(),
888 ExcMessage("Not all dependent variables have been registered."));
889
890 Assert(
891 this->n_dependent_variables() == 1,
893 "The EnergyFunctional class expects there to be only one dependent variable."));
894
895 // We can neglect correctly initializing the entries as
896 // we'll be overwriting them immediately in the succeeding call to
897 // Drivers::gradient().
898 if (gradient.size() != this->n_independent_variables())
899 gradient.reinit(this->n_independent_variables(),
900 true /*omit_zeroing_entries*/);
901
903 {
904 Assert(this->active_tape_index() !=
906 ExcMessage("Invalid tape index"));
907 Assert(this->is_recording() == false,
909 "Cannot compute gradient while tape is being recorded."));
910 Assert(this->independent_variable_values.size() ==
911 this->n_independent_variables(),
912 ExcDimensionMismatch(this->independent_variable_values.size(),
913 this->n_independent_variables()));
914
915 this->taped_driver.gradient(this->active_tape_index(),
916 this->independent_variable_values,
917 gradient);
918 }
919 else
920 {
923 Assert(this->independent_variables.size() ==
924 this->n_independent_variables(),
925 ExcDimensionMismatch(this->independent_variables.size(),
926 this->n_independent_variables()));
927
928 this->tapeless_driver.gradient(this->independent_variables,
929 this->dependent_variables,
930 gradient);
931 }
932 }
933
934
935
936 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
937 void
939 FullMatrix<scalar_type> &hessian) const
940 {
943 "Cannot computed function Hessian: AD number type does "
944 "not support the calculation of second order derivatives."));
945
947 this->taped_driver.keep_independent_values() == false))
948 {
949 Assert(
950 this->n_registered_independent_variables() ==
951 this->n_independent_variables(),
953 "Not all values of sensitivities have been registered or subsequently set!"));
954 }
955 Assert(this->n_registered_dependent_variables() ==
956 this->n_dependent_variables(),
957 ExcMessage("Not all dependent variables have been registered."));
958
959 Assert(
960 this->n_dependent_variables() == 1,
962 "The EnergyFunctional class expects there to be only one dependent variable."));
963
964 // We can neglect correctly initializing the entries as
965 // we'll be overwriting them immediately in the succeeding call to
966 // Drivers::hessian().
967 if (hessian.m() != this->n_independent_variables() ||
968 hessian.n() != this->n_independent_variables())
969 hessian.reinit({this->n_independent_variables(),
970 this->n_independent_variables()},
971 true /*omit_default_initialization*/);
972
973 if (ADNumberTraits<ad_type>::is_taped == true)
974 {
975 Assert(this->active_tape_index() !=
977 ExcMessage("Invalid tape index"));
978 Assert(this->is_recording() == false,
980 "Cannot compute hessian while tape is being recorded."));
981 Assert(this->independent_variable_values.size() ==
982 this->n_independent_variables(),
983 ExcDimensionMismatch(this->independent_variable_values.size(),
984 this->n_independent_variables()));
985
986 this->taped_driver.hessian(this->active_tape_index(),
987 this->independent_variable_values,
988 hessian);
989 }
990 else
991 {
994 Assert(this->independent_variables.size() ==
995 this->n_independent_variables(),
996 ExcDimensionMismatch(this->independent_variables.size(),
997 this->n_independent_variables()));
998
999 this->tapeless_driver.hessian(this->independent_variables,
1000 this->dependent_variables,
1001 hessian);
1002 }
1003 }
1004
1005
1006 /* ------------------- ResidualLinearization ------------------- */
1007
1008
1009
1010 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
1012 const unsigned int n_independent_variables,
1013 const unsigned int n_dependent_variables)
1014 : CellLevelBase<ADNumberTypeCode, ScalarType>(n_independent_variables,
1015 n_dependent_variables)
1016 {}
1017
1018
1019
1020 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
1021 void
1023 register_residual_vector(const std::vector<ad_type> &residual)
1024 {
1025 Assert(residual.size() == this->n_dependent_variables(),
1026 ExcMessage(
1027 "Vector size does not match number of dependent variables"));
1028 for (unsigned int i = 0; i < this->n_dependent_variables(); ++i)
1030 i, residual[i]);
1031 }
1032
1033
1034
1035 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
1036 void
1039 {
1040 if ((ADNumberTraits<ad_type>::is_taped == true &&
1041 this->taped_driver.keep_independent_values() == false) ||
1043 {
1044 Assert(
1045 this->n_registered_independent_variables() ==
1046 this->n_independent_variables(),
1047 ExcMessage(
1048 "Not all values of sensitivities have been registered or subsequently set!"));
1049 }
1050 Assert(this->n_registered_dependent_variables() ==
1051 this->n_dependent_variables(),
1052 ExcMessage("Not all dependent variables have been registered."));
1053
1054 // We can neglect correctly initializing the entries as
1055 // we'll be overwriting them immediately in the succeeding call to
1056 // Drivers::values().
1057 if (values.size() != this->n_dependent_variables())
1058 values.reinit(this->n_dependent_variables(),
1059 true /*omit_zeroing_entries*/);
1060
1062 {
1063 Assert(this->active_tape_index() !=
1065 ExcMessage("Invalid tape index"));
1066 Assert(this->is_recording() == false,
1067 ExcMessage(
1068 "Cannot compute values while tape is being recorded."));
1069 Assert(this->independent_variable_values.size() ==
1070 this->n_independent_variables(),
1071 ExcDimensionMismatch(this->independent_variable_values.size(),
1072 this->n_independent_variables()));
1073
1074 this->taped_driver.values(this->active_tape_index(),
1075 this->n_dependent_variables(),
1076 this->independent_variable_values,
1077 values);
1078 }
1079 else
1080 {
1083 this->tapeless_driver.values(this->dependent_variables, values);
1084 }
1085 }
1086
1087
1088
1089 template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
1090 void
1092 FullMatrix<scalar_type> &jacobian) const
1093 {
1094 if ((ADNumberTraits<ad_type>::is_taped == true &&
1095 this->taped_driver.keep_independent_values() == false) ||
1097 {
1098 Assert(
1099 this->n_registered_independent_variables() ==
1100 this->n_independent_variables(),
1101 ExcMessage(
1102 "Not all values of sensitivities have been registered or subsequently set!"));
1103 }
1104 Assert(this->n_registered_dependent_variables() ==
1105 this->n_dependent_variables(),
1106 ExcMessage("Not all dependent variables have been registered."));
1107
1108 // We can neglect correctly initializing the entries as
1109 // we'll be overwriting them immediately in the succeeding call to
1110 // Drivers::jacobian().
1111 if (jacobian.m() != this->n_dependent_variables() ||
1112 jacobian.n() != this->n_independent_variables())
1113 jacobian.reinit({this->n_dependent_variables(),
1114 this->n_independent_variables()},
1115 true /*omit_default_initialization*/);
1116
1118 {
1119 Assert(this->active_tape_index() !=
1121 ExcMessage("Invalid tape index"));
1122 Assert(this->is_recording() == false,
1123 ExcMessage(
1124 "Cannot compute hessian while tape is being recorded."));
1125 Assert(this->independent_variable_values.size() ==
1126 this->n_independent_variables(),
1127 ExcDimensionMismatch(this->independent_variable_values.size(),
1128 this->n_independent_variables()));
1129
1130 this->taped_driver.jacobian(this->active_tape_index(),
1131 this->n_dependent_variables(),
1132 this->independent_variable_values,
1133 jacobian);
1134 }
1135 else
1136 {
1139 Assert(this->independent_variables.size() ==
1140 this->n_independent_variables(),
1141 ExcDimensionMismatch(this->independent_variables.size(),
1142 this->n_independent_variables()));
1143
1144 this->tapeless_driver.jacobian(this->independent_variables,
1145 this->dependent_variables,
1146 jacobian);
1147 }
1148 }
1149
1150
1151
1152 /* ----------------- PointLevelFunctionsBase ----------------- */
1153
1154
1155
1156 template <int dim,
1157 enum AD::NumberTypes ADNumberTypeCode,
1158 typename ScalarType>
1160 PointLevelFunctionsBase(const unsigned int n_independent_variables,
1161 const unsigned int n_dependent_variables)
1162 : HelperBase<ADNumberTypeCode, ScalarType>(n_independent_variables,
1163 n_dependent_variables)
1164 , symmetric_independent_variables(n_independent_variables, false)
1165 {}
1166
1167
1168
1169 template <int dim,
1170 enum AD::NumberTypes ADNumberTypeCode,
1171 typename ScalarType>
1172 void
1174 const unsigned int n_independent_variables,
1175 const unsigned int n_dependent_variables,
1176 const bool clear_registered_tapes)
1177 {
1179 n_dependent_variables,
1180 clear_registered_tapes);
1181
1182 const unsigned int new_n_independent_variables =
1183 (n_independent_variables != ::numbers::invalid_unsigned_int ?
1184 n_independent_variables :
1185 this->n_independent_variables());
1186 symmetric_independent_variables =
1187 std::vector<bool>(new_n_independent_variables, false);
1188 }
1189
1190
1191
1192 template <int dim,
1193 enum AD::NumberTypes ADNumberTypeCode,
1194 typename ScalarType>
1195 bool
1197 is_symmetric_independent_variable(const unsigned int index) const
1198 {
1199 Assert(index < symmetric_independent_variables.size(),
1201 return symmetric_independent_variables[index];
1202 }
1203
1204
1205
1206 template <int dim,
1207 enum AD::NumberTypes ADNumberTypeCode,
1208 typename ScalarType>
1209 unsigned int
1212 {
1213 return std::count(symmetric_independent_variables.begin(),
1214 symmetric_independent_variables.end(),
1215 true);
1216 }
1217
1218
1219
1220 template <int dim,
1221 enum AD::NumberTypes ADNumberTypeCode,
1222 typename ScalarType>
1223 void
1225 register_independent_variables(const std::vector<scalar_type> &values)
1226 {
1227 // This is actually the same thing the set_independent_variable function,
1228 // in the sense that we simply populate our array of independent values
1229 // with a meaningful number. However, in this case we need to double check
1230 // that we're not registering these variables twice
1231 Assert(values.size() == this->n_independent_variables(),
1232 ExcMessage(
1233 "Vector size does not match number of independent variables"));
1234 for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
1235 {
1236 Assert(this->registered_independent_variable_values[i] == false,
1237 ExcMessage("Independent variable value already registered."));
1238 }
1239 set_independent_variables(values);
1240 }
1241
1242
1243
1244 template <int dim,
1245 enum AD::NumberTypes ADNumberTypeCode,
1246 typename ScalarType>
1247 const std::vector<typename PointLevelFunctionsBase<dim,
1248 ADNumberTypeCode,
1249 ScalarType>::ad_type> &
1252 {
1254 {
1255 Assert(this->active_tape_index() !=
1257 ExcMessage("Invalid tape index"));
1258 }
1259
1260 // Just in case the user has not done so, we repeat the call to
1261 // initialize the internally stored vector of AD numbers that
1262 // represents the independent variables.
1263 this->finalize_sensitive_independent_variables();
1264 Assert(this->independent_variables.size() ==
1265 this->n_independent_variables(),
1266 ExcDimensionMismatch(this->independent_variables.size(),
1267 this->n_independent_variables()));
1268
1269 return this->independent_variables;
1270 }
1271
1272
1273
1274 template <int dim,
1275 enum AD::NumberTypes ADNumberTypeCode,
1276 typename ScalarType>
1277 void
1279 set_sensitivity_value(const unsigned int index,
1280 const bool symmetric_component,
1281 const scalar_type &value)
1282 {
1284 value);
1285 Assert(
1286 index < this->n_independent_variables(),
1287 ExcMessage(
1288 "Trying to set the symmetry flag of a non-existent independent variable."));
1289 Assert(index < symmetric_independent_variables.size(),
1291 symmetric_independent_variables[index] = symmetric_component;
1292 }
1293
1294
1295
1296 template <int dim,
1297 enum AD::NumberTypes ADNumberTypeCode,
1298 typename ScalarType>
1299 void
1301 set_independent_variables(const std::vector<scalar_type> &values)
1302 {
1304 {
1305 Assert(this->active_tape_index() !=
1307 ExcMessage("Invalid tape index"));
1308 }
1309 Assert(values.size() == this->n_independent_variables(),
1310 ExcMessage(
1311 "Vector size does not match number of independent variables"));
1312 for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
1314 i, values[i]);
1315 }
1316
1317
1318
1319 /* -------------------- ScalarFunction -------------------- */
1320
1321
1322
1323 template <int dim,
1324 enum AD::NumberTypes ADNumberTypeCode,
1325 typename ScalarType>
1327 const unsigned int n_independent_variables)
1328 : PointLevelFunctionsBase<dim, ADNumberTypeCode, ScalarType>(
1329 n_independent_variables,
1330 1)
1331 {}
1332
1333
1334
1335 template <int dim,
1336 enum AD::NumberTypes ADNumberTypeCode,
1337 typename ScalarType>
1338 void
1341 {
1342 Assert(this->n_dependent_variables() == 1, ExcInternalError());
1344 0, func);
1345 }
1346
1347
1348
1349 template <int dim,
1350 enum AD::NumberTypes ADNumberTypeCode,
1351 typename ScalarType>
1354 {
1355 if ((ADNumberTraits<ad_type>::is_taped == true &&
1356 this->taped_driver.keep_independent_values() == false) ||
1358 {
1359 Assert(
1360 this->n_registered_independent_variables() ==
1361 this->n_independent_variables(),
1362 ExcMessage(
1363 "Not all values of sensitivities have been registered or subsequently set!"));
1364 }
1365 Assert(this->n_registered_dependent_variables() ==
1366 this->n_dependent_variables(),
1367 ExcMessage("Not all dependent variables have been registered."));
1368
1369 Assert(
1370 this->n_dependent_variables() == 1,
1371 ExcMessage(
1372 "The ScalarFunction class expects there to be only one dependent variable."));
1373
1375 {
1376 Assert(this->active_tape_index() !=
1378 ExcMessage("Invalid tape index"));
1379 Assert(this->is_recording() == false,
1380 ExcMessage(
1381 "Cannot compute values while tape is being recorded."));
1382 Assert(this->independent_variable_values.size() ==
1383 this->n_independent_variables(),
1384 ExcDimensionMismatch(this->independent_variable_values.size(),
1385 this->n_independent_variables()));
1386
1387 return this->taped_driver.value(this->active_tape_index(),
1388 this->independent_variable_values);
1389 }
1390 else
1391 {
1394 return this->tapeless_driver.value(this->dependent_variables);
1395 }
1396 }
1397
1398
1399 template <int dim,
1400 enum AD::NumberTypes ADNumberTypeCode,
1401 typename ScalarType>
1402 void
1404 Vector<scalar_type> &gradient) const
1405 {
1406 if ((ADNumberTraits<ad_type>::is_taped == true &&
1407 this->taped_driver.keep_independent_values() == false) ||
1409 {
1410 Assert(
1411 this->n_registered_independent_variables() ==
1412 this->n_independent_variables(),
1413 ExcMessage(
1414 "Not all values of sensitivities have been registered or subsequently set!"));
1415 }
1416 Assert(this->n_registered_dependent_variables() ==
1417 this->n_dependent_variables(),
1418 ExcMessage("Not all dependent variables have been registered."));
1419
1420 Assert(
1421 this->n_dependent_variables() == 1,
1422 ExcMessage(
1423 "The ScalarFunction class expects there to be only one dependent variable."));
1424
1425 // We can neglect correctly initializing the entries as
1426 // we'll be overwriting them immediately in the succeeding call to
1427 // Drivers::gradient().
1428 if (gradient.size() != this->n_independent_variables())
1429 gradient.reinit(this->n_independent_variables(),
1430 true /*omit_zeroing_entries*/);
1431
1433 {
1434 Assert(this->active_tape_index() !=
1436 ExcMessage("Invalid tape index"));
1437 Assert(this->is_recording() == false,
1438 ExcMessage(
1439 "Cannot compute gradient while tape is being recorded."));
1440 Assert(this->independent_variable_values.size() ==
1441 this->n_independent_variables(),
1442 ExcDimensionMismatch(this->independent_variable_values.size(),
1443 this->n_independent_variables()));
1444
1445 this->taped_driver.gradient(this->active_tape_index(),
1446 this->independent_variable_values,
1447 gradient);
1448 }
1449 else
1450 {
1453 Assert(this->independent_variables.size() ==
1454 this->n_independent_variables(),
1455 ExcDimensionMismatch(this->independent_variables.size(),
1456 this->n_independent_variables()));
1457
1458 this->tapeless_driver.gradient(this->independent_variables,
1459 this->dependent_variables,
1460 gradient);
1461 }
1462
1463 // Account for symmetries of tensor components
1464 for (unsigned int i = 0; i < this->n_independent_variables(); i++)
1465 {
1466 if (this->is_symmetric_independent_variable(i) == true)
1467 gradient[i] *= 0.5;
1468 }
1469 }
1470
1471
1472
1473 template <int dim,
1474 enum AD::NumberTypes ADNumberTypeCode,
1475 typename ScalarType>
1476 void
1478 FullMatrix<scalar_type> &hessian) const
1479 {
1481 ExcMessage(
1482 "Cannot computed function Hessian: AD number type does "
1483 "not support the calculation of second order derivatives."));
1484
1485 if ((ADNumberTraits<ad_type>::is_taped == true &&
1486 this->taped_driver.keep_independent_values() == false))
1487 {
1488 Assert(
1489 this->n_registered_independent_variables() ==
1490 this->n_independent_variables(),
1491 ExcMessage(
1492 "Not all values of sensitivities have been registered or subsequently set!"));
1493 }
1494 Assert(this->n_registered_dependent_variables() ==
1495 this->n_dependent_variables(),
1496 ExcMessage("Not all dependent variables have been registered."));
1497
1498 Assert(
1499 this->n_dependent_variables() == 1,
1500 ExcMessage(
1501 "The ScalarFunction class expects there to be only one dependent variable."));
1502
1503 // We can neglect correctly initializing the entries as
1504 // we'll be overwriting them immediately in the succeeding call to
1505 // Drivers::hessian().
1506 if (hessian.m() != this->n_independent_variables() ||
1507 hessian.n() != this->n_independent_variables())
1508 hessian.reinit({this->n_independent_variables(),
1509 this->n_independent_variables()},
1510 true /*omit_default_initialization*/);
1511
1513 {
1514 Assert(this->active_tape_index() !=
1516 ExcMessage("Invalid tape index"));
1517 Assert(this->is_recording() == false,
1518 ExcMessage(
1519 "Cannot compute Hessian while tape is being recorded."));
1520 Assert(this->independent_variable_values.size() ==
1521 this->n_independent_variables(),
1522 ExcDimensionMismatch(this->independent_variable_values.size(),
1523 this->n_independent_variables()));
1524
1525 this->taped_driver.hessian(this->active_tape_index(),
1526 this->independent_variable_values,
1527 hessian);
1528 }
1529 else
1530 {
1533 Assert(this->independent_variables.size() ==
1534 this->n_independent_variables(),
1535 ExcDimensionMismatch(this->independent_variables.size(),
1536 this->n_independent_variables()));
1537
1538 this->tapeless_driver.hessian(this->independent_variables,
1539 this->dependent_variables,
1540 hessian);
1541 }
1542
1543 // Account for symmetries of tensor components
1544 for (unsigned int i = 0; i < this->n_independent_variables(); i++)
1545 for (unsigned int j = 0; j < i + 1; j++)
1546 {
1547 if (this->is_symmetric_independent_variable(i) == true &&
1548 this->is_symmetric_independent_variable(j) == true)
1549 {
1550 hessian[i][j] *= 0.25;
1551 if (i != j)
1552 hessian[j][i] *= 0.25;
1553 }
1554 else if ((this->is_symmetric_independent_variable(i) == true &&
1555 this->is_symmetric_independent_variable(j) == false) ||
1556 (this->is_symmetric_independent_variable(j) == true &&
1557 this->is_symmetric_independent_variable(i) == false))
1558 {
1559 hessian[i][j] *= 0.5;
1560 if (i != j)
1561 hessian[j][i] *= 0.5;
1562 }
1563 }
1564 }
1565
1566
1567
1568 template <int dim,
1569 enum AD::NumberTypes ADNumberTypeCode,
1570 typename ScalarType>
1571 Tensor<
1572 0,
1573 dim,
1577 const FEValuesExtractors::Scalar &extractor_row,
1578 const FEValuesExtractors::Scalar &extractor_col)
1579 {
1580 // NOTE: It is necessary to make special provision for the case when the
1581 // HessianType is scalar. Unfortunately Tensor<0,dim> does not provide
1582 // the function unrolled_to_component_indices!
1583 // NOTE: The order of components must be consistently defined throughout
1584 // this class.
1586
1587 // Get indexsets for the subblocks from which we wish to extract the
1588 // matrix values
1589 const std::vector<unsigned int> row_index_set(
1590 internal::extract_field_component_indices<dim>(extractor_row));
1591 const std::vector<unsigned int> col_index_set(
1592 internal::extract_field_component_indices<dim>(extractor_col));
1593 Assert(row_index_set.size() == 1, ExcInternalError());
1594 Assert(col_index_set.size() == 1, ExcInternalError());
1595
1597 0,
1598 hessian[row_index_set[0]][col_index_set[0]]);
1599
1600 return out;
1601 }
1602
1603
1604
1605 template <int dim,
1606 enum AD::NumberTypes ADNumberTypeCode,
1607 typename ScalarType>
1609 4,
1610 dim,
1614 const FullMatrix<scalar_type> & hessian,
1615 const FEValuesExtractors::SymmetricTensor<2> &extractor_row,
1616 const FEValuesExtractors::SymmetricTensor<2> &extractor_col)
1617 {
1618 // NOTE: The order of components must be consistently defined throughout
1619 // this class.
1620 // NOTE: We require a specialisation for rank-4 symmetric tensors because
1621 // they do not define their rank, and setting data using TableIndices is
1622 // somewhat specialised as well.
1624
1625 // Get indexsets for the subblocks from which we wish to extract the
1626 // matrix values
1627 const std::vector<unsigned int> row_index_set(
1628 internal::extract_field_component_indices<dim>(extractor_row));
1629 const std::vector<unsigned int> col_index_set(
1630 internal::extract_field_component_indices<dim>(extractor_col));
1631
1632 for (unsigned int r = 0; r < row_index_set.size(); ++r)
1633 for (unsigned int c = 0; c < col_index_set.size(); ++c)
1634 {
1636 out, r, c, hessian[row_index_set[r]][col_index_set[c]]);
1637 }
1638
1639 return out;
1640 }
1641
1642
1643
1644 /* -------------------- VectorFunction -------------------- */
1645
1646
1647
1648 template <int dim,
1649 enum AD::NumberTypes ADNumberTypeCode,
1650 typename ScalarType>
1652 const unsigned int n_independent_variables,
1653 const unsigned int n_dependent_variables)
1654 : PointLevelFunctionsBase<dim, ADNumberTypeCode, ScalarType>(
1655 n_independent_variables,
1656 n_dependent_variables)
1657 {}
1658
1659
1660
1661 template <int dim,
1662 enum AD::NumberTypes ADNumberTypeCode,
1663 typename ScalarType>
1664 void
1666 register_dependent_variables(const std::vector<ad_type> &funcs)
1667 {
1668 Assert(funcs.size() == this->n_dependent_variables(),
1669 ExcMessage(
1670 "Vector size does not match number of dependent variables"));
1671 for (unsigned int i = 0; i < this->n_dependent_variables(); ++i)
1673 i, funcs[i]);
1674 }
1675
1676
1677
1678 template <int dim,
1679 enum AD::NumberTypes ADNumberTypeCode,
1680 typename ScalarType>
1681 void
1684 {
1685 if ((ADNumberTraits<ad_type>::is_taped == true &&
1686 this->taped_driver.keep_independent_values() == false) ||
1688 {
1689 Assert(
1690 this->n_registered_independent_variables() ==
1691 this->n_independent_variables(),
1692 ExcMessage(
1693 "Not all values of sensitivities have been registered or subsequently set!"));
1694 }
1695 Assert(this->n_registered_dependent_variables() ==
1696 this->n_dependent_variables(),
1697 ExcMessage("Not all dependent variables have been registered."));
1698
1699 // We can neglect correctly initializing the entries as
1700 // we'll be overwriting them immediately in the succeeding call to
1701 // Drivers::values().
1702 if (values.size() != this->n_dependent_variables())
1703 values.reinit(this->n_dependent_variables(),
1704 true /*omit_zeroing_entries*/);
1705
1707 {
1708 Assert(this->active_tape_index() !=
1710 ExcMessage("Invalid tape index"));
1711 Assert(this->is_recording() == false,
1712 ExcMessage(
1713 "Cannot compute values while tape is being recorded."));
1714 Assert(this->independent_variable_values.size() ==
1715 this->n_independent_variables(),
1716 ExcDimensionMismatch(this->independent_variable_values.size(),
1717 this->n_independent_variables()));
1718
1719 this->taped_driver.values(this->active_tape_index(),
1720 this->n_dependent_variables(),
1721 this->independent_variable_values,
1722 values);
1723 }
1724 else
1725 {
1728 this->tapeless_driver.values(this->dependent_variables, values);
1729 }
1730 }
1731
1732
1733
1734 template <int dim,
1735 enum AD::NumberTypes ADNumberTypeCode,
1736 typename ScalarType>
1737 void
1739 FullMatrix<scalar_type> &jacobian) const
1740 {
1741 if ((ADNumberTraits<ad_type>::is_taped == true &&
1742 this->taped_driver.keep_independent_values() == false) ||
1744 {
1745 Assert(
1746 this->n_registered_independent_variables() ==
1747 this->n_independent_variables(),
1748 ExcMessage(
1749 "Not all values of sensitivities have been registered or subsequently set!"));
1750 }
1751 Assert(this->n_registered_dependent_variables() ==
1752 this->n_dependent_variables(),
1753 ExcMessage("Not all dependent variables have been registered."));
1754
1755 // We can neglect correctly initializing the entries as
1756 // we'll be overwriting them immediately in the succeeding call to
1757 // Drivers::jacobian().
1758 if (jacobian.m() != this->n_dependent_variables() ||
1759 jacobian.n() != this->n_independent_variables())
1760 jacobian.reinit({this->n_dependent_variables(),
1761 this->n_independent_variables()},
1762 true /*omit_default_initialization*/);
1763
1765 {
1766 Assert(this->active_tape_index() !=
1768 ExcMessage("Invalid tape index"));
1769 Assert(this->is_recording() == false,
1770 ExcMessage(
1771 "Cannot compute Jacobian while tape is being recorded."));
1772 Assert(this->independent_variable_values.size() ==
1773 this->n_independent_variables(),
1774 ExcDimensionMismatch(this->independent_variable_values.size(),
1775 this->n_independent_variables()));
1776
1777 this->taped_driver.jacobian(this->active_tape_index(),
1778 this->n_dependent_variables(),
1779 this->independent_variable_values,
1780 jacobian);
1781 }
1782 else
1783 {
1786 Assert(this->independent_variables.size() ==
1787 this->n_independent_variables(),
1788 ExcDimensionMismatch(this->independent_variables.size(),
1789 this->n_independent_variables()));
1790
1791 this->tapeless_driver.jacobian(this->independent_variables,
1792 this->dependent_variables,
1793 jacobian);
1794 }
1795
1796 for (unsigned int j = 0; j < this->n_independent_variables(); j++)
1797 {
1798 // Because we perform just a single differentiation
1799 // operation with respect to the "column" variables,
1800 // we only need to consider them for symmetry conditions.
1801 if (this->is_symmetric_independent_variable(j) == true)
1802 for (unsigned int i = 0; i < this->n_dependent_variables(); i++)
1803 jacobian[i][j] *= 0.5;
1804 }
1805 }
1806
1807
1808
1809 template <int dim,
1810 enum AD::NumberTypes ADNumberTypeCode,
1811 typename ScalarType>
1812 Tensor<
1813 0,
1814 dim,
1818 const FullMatrix<scalar_type> & jacobian,
1819 const FEValuesExtractors::Scalar &extractor_row,
1820 const FEValuesExtractors::Scalar &extractor_col)
1821 {
1822 // NOTE: It is necessary to make special provision for the case when the
1823 // HessianType is scalar. Unfortunately Tensor<0,dim> does not provide
1824 // the function unrolled_to_component_indices!
1825 // NOTE: The order of components must be consistently defined throughout
1826 // this class.
1828
1829 // Get indexsets for the subblocks from which we wish to extract the
1830 // matrix values
1831 const std::vector<unsigned int> row_index_set(
1832 internal::extract_field_component_indices<dim>(extractor_row));
1833 const std::vector<unsigned int> col_index_set(
1834 internal::extract_field_component_indices<dim>(extractor_col));
1835 Assert(row_index_set.size() == 1, ExcInternalError());
1836 Assert(col_index_set.size() == 1, ExcInternalError());
1837
1839 0,
1840 jacobian[row_index_set[0]][col_index_set[0]]);
1841
1842 return out;
1843 }
1844
1845
1846
1847 template <int dim,
1848 enum AD::NumberTypes ADNumberTypeCode,
1849 typename ScalarType>
1851 4,
1852 dim,
1856 const FullMatrix<scalar_type> & jacobian,
1857 const FEValuesExtractors::SymmetricTensor<2> &extractor_row,
1858 const FEValuesExtractors::SymmetricTensor<2> &extractor_col)
1859 {
1860 // NOTE: The order of components must be consistently defined throughout
1861 // this class.
1862 // NOTE: We require a specialisation for rank-4 symmetric tensors because
1863 // they do not define their rank, and setting data using TableIndices is
1864 // somewhat specialised as well.
1866
1867 // Get indexsets for the subblocks from which we wish to extract the
1868 // matrix values
1869 const std::vector<unsigned int> row_index_set(
1870 internal::extract_field_component_indices<dim>(extractor_row));
1871 const std::vector<unsigned int> col_index_set(
1872 internal::extract_field_component_indices<dim>(extractor_col));
1873
1874 for (unsigned int r = 0; r < row_index_set.size(); ++r)
1875 for (unsigned int c = 0; c < col_index_set.size(); ++c)
1876 {
1878 out, r, c, jacobian[row_index_set[r]][col_index_set[c]]);
1879 }
1880
1881 return out;
1882 }
1883
1884
1885 } // namespace AD
1886} // namespace Differentiation
1887
1888
1889/* --- Explicit instantiations --- */
1890# include "ad_helpers.inst"
1891
1892# ifdef DEAL_II_WITH_ADOLC
1893# include "ad_helpers.inst1"
1894# endif
1895# ifdef DEAL_II_TRILINOS_WITH_SACADO
1896# include "ad_helpers.inst2"
1897# endif
1898
1899
1901
1902#endif // defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO)
const std::vector< ad_type > & get_sensitive_dof_values() const
Definition: ad_helpers.cc:749
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
Definition: ad_helpers.h:851
void set_dof_values(const std::vector< scalar_type > &dof_values)
Definition: ad_helpers.cc:774
void register_dof_values(const std::vector< scalar_type > &dof_values)
Definition: ad_helpers.cc:726
CellLevelBase(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:715
void register_energy_functional(const ad_type &energy)
Definition: ad_helpers.cc:807
EnergyFunctional(const unsigned int n_independent_variables)
Definition: ad_helpers.cc:798
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
Definition: ad_helpers.h:1233
void compute_residual(Vector< scalar_type > &residual) const override
Definition: ad_helpers.cc:873
virtual void compute_linearization(FullMatrix< scalar_type > &linearization) const override
Definition: ad_helpers.cc:938
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:1226
unsigned int n_registered_dependent_variables() const
Definition: ad_helpers.cc:250
void stop_recording_operations(const bool write_tapes_to_file=false)
Definition: ad_helpers.cc:648
typename AD::NumberTraits< ScalarType, ADNumberTypeCode >::scalar_type scalar_type
Definition: ad_helpers.h:178
Types< ad_type >::tape_index active_tape_index() const
Definition: ad_helpers.cc:284
void print_tape_stats(const typename Types< ad_type >::tape_index tape_index, std::ostream &stream) const
Definition: ad_helpers.cc:378
std::size_t n_independent_variables() const
Definition: ad_helpers.cc:241
bool start_recording_operations(const typename Types< ad_type >::tape_index tape_index, const bool overwrite_tape=false, const bool keep_independent_values=true)
Definition: ad_helpers.cc:583
void mark_independent_variable(const unsigned int index, ad_type &out) const
Definition: ad_helpers.cc:142
bool is_registered_tape(const typename Types< ad_type >::tape_index tape_index) const
Definition: ad_helpers.cc:296
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)
Definition: ad_helpers.cc:527
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)
Definition: ad_helpers.cc:564
void print(std::ostream &stream) const
Definition: ad_helpers.cc:309
void set_sensitivity_value(const unsigned int index, const scalar_type &value)
Definition: ad_helpers.cc:105
bool active_tape_requires_retaping() const
Definition: ad_helpers.cc:502
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)
Definition: ad_helpers.cc:395
bool recorded_tape_requires_retaping(const typename Types< ad_type >::tape_index tape_index) const
Definition: ad_helpers.cc:489
std::vector< ad_type > dependent_variables
Definition: ad_helpers.h:750
std::size_t n_dependent_variables() const
Definition: ad_helpers.cc:262
void print_values(std::ostream &stream) const
Definition: ad_helpers.cc:364
void register_dependent_variable(const unsigned int index, const ad_type &func)
Definition: ad_helpers.cc:684
static void configure_tapeless_mode(const unsigned int n_independent_variables, const bool ensure_persistent_setting=true)
Definition: ad_helpers.cc:453
void initialize_non_sensitive_independent_variable(const unsigned int index, ad_type &out) const
Definition: ad_helpers.cc:205
typename AD::NumberTraits< ScalarType, ADNumberTypeCode >::ad_type ad_type
Definition: ad_helpers.h:185
void activate_recorded_tape(const typename Types< ad_type >::tape_index tape_index)
Definition: ad_helpers.cc:479
bool is_symmetric_independent_variable(const unsigned int index) const
Definition: ad_helpers.cc:1197
PointLevelFunctionsBase(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:1160
void set_independent_variables(const std::vector< scalar_type > &values)
Definition: ad_helpers.cc:1301
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:2649
unsigned int n_symmetric_independent_variables() const
Definition: ad_helpers.cc:1211
void register_independent_variables(const std::vector< scalar_type > &values)
Definition: ad_helpers.cc:1225
void set_sensitivity_value(const unsigned int index, const bool symmetric_component, const scalar_type &value)
Definition: ad_helpers.cc:1279
const std::vector< ad_type > & get_sensitive_variables() const
Definition: ad_helpers.cc:1251
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
Definition: ad_helpers.cc:1173
ResidualLinearization(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:1011
virtual void compute_residual(Vector< scalar_type > &residual) const override
Definition: ad_helpers.cc:1037
virtual void compute_linearization(FullMatrix< scalar_type > &linearization) const override
Definition: ad_helpers.cc:1091
void register_residual_vector(const std::vector< ad_type > &residual)
Definition: ad_helpers.cc:1023
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
Definition: ad_helpers.h:3095
void register_dependent_variable(const ad_type &func)
Definition: ad_helpers.cc:1340
ScalarFunction(const unsigned int n_independent_variables)
Definition: ad_helpers.cc:1326
void compute_gradient(Vector< scalar_type > &gradient) const
Definition: ad_helpers.cc:1403
void compute_hessian(FullMatrix< scalar_type > &hessian) const
Definition: ad_helpers.cc:1477
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:3088
void compute_values(Vector< scalar_type > &values) const
Definition: ad_helpers.cc:1682
void register_dependent_variables(const std::vector< ad_type > &funcs)
Definition: ad_helpers.cc:1666
void compute_jacobian(FullMatrix< scalar_type > &jacobian) const
Definition: ad_helpers.cc:1738
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:3479
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)
Definition: ad_helpers.cc:1651
size_type n() const
size_type m() const
Definition: tensor.h:472
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
std::enable_if<!(ADNumberTraits< ADNumberType >::type_code==NumberTypes::adolc_tapeless)>::type configure_tapeless_mode(const unsigned int)
Definition: ad_drivers.cc:1630
void set_tensor_entry(TensorType &t, const unsigned int unrolled_index, const NumberType &value)
Definition: ad_helpers.h:2437
static const unsigned int invalid_unsigned_int
Definition: types.h:196
static const Types< ADNumberType >::tape_index invalid_tape_index
Definition: ad_drivers.h:121
static void initialize_global_environment(const unsigned int n_independent_variables)
Definition: ad_drivers.cc:1491