Reference documentation for deal.II version GIT 4e92747666 2022-12-02 08:40:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 
29 namespace 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,
54  ExcMessage(
55  "Floating point/arithmetic numbers have no derivatives."));
56  Assert(
58  ExcMessage(
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,
118  ExcMessage(
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(),
131  ExcMessage(
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,
148  ExcInternalError());
149 
150  if (index > 0)
151  {
152  Assert(
153  registered_marked_independent_variables[index - 1] == true,
154  ExcMessage(
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,
163  ExcMessage(
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,
214  ExcMessage(
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,
220  ExcInternalError());
221 
222  out = independent_variable_values[index];
223  }
224 
225 
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
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  }
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)
394  {
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;
467  }
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  }
479 
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  {
551  ExcInternalError());
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 
575 
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,
613  ExcMessage(
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,
624  ExcInternalError());
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  {
660  ExcInternalError());
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  }
674 
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,
685  ExcMessage(
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,
693  ExcMessage(
694  "Must be recording when registering dependent variables."));
695  }
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 
708 
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(),
729  ExcMessage(
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;
763  }
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(),
779  ExcMessage(
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  {
816  if ((ADNumberTraits<ad_type>::is_taped == true &&
817  this->taped_driver.keep_independent_values() == false) ||
819  {
820  Assert(
821  this->n_registered_independent_variables() ==
822  this->n_independent_variables(),
823  ExcMessage(
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,
832  ExcMessage(
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,
841  ExcMessage(
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  {
854  ExcInternalError());
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(),
878  ExcMessage(
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,
887  ExcMessage(
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,
903  ExcMessage(
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  {
917  ExcInternalError());
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  {
937  ExcMessage(
938  "Cannot computed function Hessian: AD number type does "
939  "not support the calculation of second order derivatives."));
940 
941  if ((ADNumberTraits<ad_type>::is_taped == true &&
942  this->taped_driver.keep_independent_values() == false))
943  {
944  Assert(
945  this->n_registered_independent_variables() ==
946  this->n_independent_variables(),
947  ExcMessage(
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,
956  ExcMessage(
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,
974  ExcMessage(
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()));
980 
981  this->taped_driver.hessian(this->active_tape_index(),
982  this->independent_variable_values,
983  hessian);
984  }
985  else
986  {
988  ExcInternalError());
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  {
1077  ExcInternalError());
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  {
1133  ExcInternalError());
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  {
1173  HelperBase<ADNumberTypeCode, ScalarType>::reset(n_independent_variables,
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(),
1195  ExcInternalError());
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(),
1285  ExcInternalError());
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  {
1388  ExcInternalError());
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  {
1447  ExcInternalError());
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  {
1527  ExcInternalError());
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  {
1722  ExcInternalError());
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  {
1780  ExcInternalError());
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
Definition: ad_helpers.cc:744
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
Definition: ad_helpers.h:850
void set_dof_values(const std::vector< scalar_type > &dof_values)
Definition: ad_helpers.cc:769
void register_dof_values(const std::vector< scalar_type > &dof_values)
Definition: ad_helpers.cc:721
CellLevelBase(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:710
void register_energy_functional(const ad_type &energy)
Definition: ad_helpers.cc:802
EnergyFunctional(const unsigned int n_independent_variables)
Definition: ad_helpers.cc:793
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
Definition: ad_helpers.h:1232
void compute_residual(Vector< scalar_type > &residual) const override
Definition: ad_helpers.cc:868
virtual void compute_linearization(FullMatrix< scalar_type > &linearization) const override
Definition: ad_helpers.cc:933
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:1225
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:643
typename AD::NumberTraits< ScalarType, ADNumberTypeCode >::scalar_type scalar_type
Definition: ad_helpers.h:177
void print_tape_stats(const typename Types< ad_type >::tape_index tape_index, std::ostream &stream) const
Definition: ad_helpers.cc:373
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:578
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:522
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:559
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:497
Types< ad_type >::tape_index active_tape_index() const
Definition: ad_helpers.cc:284
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:390
bool recorded_tape_requires_retaping(const typename Types< ad_type >::tape_index tape_index) const
Definition: ad_helpers.cc:484
std::vector< ad_type > dependent_variables
Definition: ad_helpers.h:749
std::size_t n_dependent_variables() const
Definition: ad_helpers.cc:262
void print_values(std::ostream &stream) const
Definition: ad_helpers.cc:359
void register_dependent_variable(const unsigned int index, const ad_type &func)
Definition: ad_helpers.cc:679
static void configure_tapeless_mode(const unsigned int n_independent_variables, const bool ensure_persistent_setting=true)
Definition: ad_helpers.cc:448
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:184
void activate_recorded_tape(const typename Types< ad_type >::tape_index tape_index)
Definition: ad_helpers.cc:474
bool is_symmetric_independent_variable(const unsigned int index) const
Definition: ad_helpers.cc:1192
PointLevelFunctionsBase(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:1155
void set_independent_variables(const std::vector< scalar_type > &values)
Definition: ad_helpers.cc:1296
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:2655
unsigned int n_symmetric_independent_variables() const
Definition: ad_helpers.cc:1206
void register_independent_variables(const std::vector< scalar_type > &values)
Definition: ad_helpers.cc:1220
void set_sensitivity_value(const unsigned int index, const bool symmetric_component, const scalar_type &value)
Definition: ad_helpers.cc:1274
const std::vector< ad_type > & get_sensitive_variables() const
Definition: ad_helpers.cc:1246
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:1168
ResidualLinearization(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:1006
virtual void compute_residual(Vector< scalar_type > &residual) const override
Definition: ad_helpers.cc:1032
virtual void compute_linearization(FullMatrix< scalar_type > &linearization) const override
Definition: ad_helpers.cc:1086
void register_residual_vector(const std::vector< ad_type > &residual)
Definition: ad_helpers.cc:1018
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
Definition: ad_helpers.h:3133
void register_dependent_variable(const ad_type &func)
Definition: ad_helpers.cc:1335
ScalarFunction(const unsigned int n_independent_variables)
Definition: ad_helpers.cc:1321
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)
void compute_gradient(Vector< scalar_type > &gradient) const
Definition: ad_helpers.cc:1398
void compute_hessian(FullMatrix< scalar_type > &hessian) const
Definition: ad_helpers.cc:1472
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:3126
void compute_values(Vector< scalar_type > &values) const
Definition: ad_helpers.cc:1677
void register_dependent_variables(const std::vector< ad_type > &funcs)
Definition: ad_helpers.cc:1661
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)
void compute_jacobian(FullMatrix< scalar_type > &jacobian) const
Definition: ad_helpers.cc:1733
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:3517
VectorFunction(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:1646
size_type n() const
size_type m() const
Definition: tensor.h:503
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcMessage(std::string arg1)
void set_tensor_entry(TensorType &t, const unsigned int unrolled_index, const NumberType &value)
Definition: ad_helpers.h:2436
std::enable_if_t<!(ADNumberTraits< ADNumberType >::type_code==NumberTypes::adolc_tapeless)> configure_tapeless_mode(const unsigned int)
Definition: ad_drivers.cc:1573
void set_dof_values(const DoFCellAccessor< dim, spacedim, lda > &cell, const Vector< number > &local_values, OutputVector &values, const bool perform_check)
static const unsigned int invalid_unsigned_int
Definition: types.h:206
static void initialize_global_environment(const unsigned int n_independent_variables)
Definition: ad_drivers.cc:1435
unsigned int tape_buffer_sizes
Definition: ad_drivers.h:105