Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 
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  {
65  configure_tapeless_mode(n_independent_variables,
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.
72  dependent_variables.resize(n_dependent_variables,
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: "
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;
341 
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  }
359 
360 
361 
362  template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
363  void
365  std::ostream &stream) const
366  {
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 
523 
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  {
556  ExcInternalError());
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,
618  ExcMessage(
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  {
629  ExcInternalError());
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  {
665  ExcInternalError());
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,
690  ExcMessage(
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,
698  ExcMessage(
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 
709 
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  {}
721 
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(),
734  ExcMessage(
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 
770 
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(),
784  ExcMessage(
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  {
821  if ((ADNumberTraits<ad_type>::is_taped == true &&
822  this->taped_driver.keep_independent_values() == false) ||
824  {
825  Assert(
826  this->n_registered_independent_variables() ==
827  this->n_independent_variables(),
828  ExcMessage(
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,
837  ExcMessage(
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,
846  ExcMessage(
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  {
859  ExcInternalError());
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  {
876  if ((ADNumberTraits<ad_type>::is_taped == true &&
877  this->taped_driver.keep_independent_values() == false) ||
879  {
880  Assert(
881  this->n_registered_independent_variables() ==
882  this->n_independent_variables(),
883  ExcMessage(
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,
892  ExcMessage(
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,
908  ExcMessage(
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  {
922  ExcInternalError());
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  {
942  ExcMessage(
943  "Cannot computed function Hessian: AD number type does "
944  "not support the calculation of second order derivatives."));
945 
946  if ((ADNumberTraits<ad_type>::is_taped == true &&
947  this->taped_driver.keep_independent_values() == false))
948  {
949  Assert(
950  this->n_registered_independent_variables() ==
951  this->n_independent_variables(),
952  ExcMessage(
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,
961  ExcMessage(
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 
974  {
975  Assert(this->active_tape_index() !=
977  ExcMessage("Invalid tape index"));
978  Assert(this->is_recording() == false,
979  ExcMessage(
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  {
993  ExcInternalError());
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
1038  Vector<scalar_type> &values) const
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  {
1082  ExcInternalError());
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  {
1138  ExcInternalError());
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  {
1178  HelperBase<ADNumberTypeCode, ScalarType>::reset(n_independent_variables,
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(),
1200  ExcInternalError());
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(),
1290  ExcInternalError());
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  {
1393  ExcInternalError());
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  {
1452  ExcInternalError());
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  {
1532  ExcInternalError());
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
1683  Vector<scalar_type> &values) const
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  {
1727  ExcInternalError());
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  {
1785  ExcInternalError());
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)
Differentiation::AD::ScalarFunction::scalar_type
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:3096
Differentiation::AD::is_tapeless_ad_number
Definition: ad_number_traits.h:272
Differentiation::AD::is_adolc_tapeless_number
Definition: adolc_number_types.h:62
Differentiation::AD::HelperBase::set_tape_buffer_sizes
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
SymmetricTensor
Definition: symmetric_tensor.h:611
Differentiation::AD::VectorFunction::compute_jacobian
void compute_jacobian(FullMatrix< scalar_type > &jacobian) const
Definition: ad_helpers.cc:1738
Differentiation::AD::VectorFunction::extract_jacobian_component
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)
Differentiation::AD::EnergyFunctional::compute_energy
scalar_type compute_energy() const
Definition: ad_helpers.cc:819
FEValuesExtractors::Scalar
Definition: fe_values_extractors.h:95
ad_helpers.h
FullMatrix::m
size_type m() const
Differentiation::AD::internal::set_tensor_entry
void set_tensor_entry(TensorType &t, const unsigned int unrolled_index, const NumberType &value)
Definition: ad_helpers.h:2441
TableBase< N, number >::reinit
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
FullMatrix::n
size_type n() const
Differentiation::AD::HelperBase::active_tape_index
Types< ad_type >::tape_index active_tape_index() const
Definition: ad_helpers.cc:284
Differentiation::AD::EnergyFunctional::compute_residual
void compute_residual(Vector< scalar_type > &residual) const override
Definition: ad_helpers.cc:873
Differentiation::AD::ScalarFunction::compute_hessian
void compute_hessian(FullMatrix< scalar_type > &hessian) const
Definition: ad_helpers.cc:1477
Differentiation::AD::ScalarFunction
Definition: ad_helpers.h:3087
Differentiation::AD::PointLevelFunctionsBase::PointLevelFunctionsBase
PointLevelFunctionsBase(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:1160
Differentiation::AD::HelperBase< ADNumberTypeCode, double >::ad_type
typename AD::NumberTraits< double, ADNumberTypeCode >::ad_type ad_type
Definition: ad_helpers.h:187
Differentiation::AD::HelperBase
Definition: ad_helpers.h:172
Differentiation::AD::internal::NumberType
Definition: numbers.h:653
FEValuesExtractors::SymmetricTensor
Definition: fe_values_extractors.h:199
Differentiation::AD::EnergyFunctional::register_energy_functional
void register_energy_functional(const ad_type &energy)
Definition: ad_helpers.cc:807
Differentiation::AD::VectorFunction::compute_values
void compute_values(Vector< scalar_type > &values) const
Definition: ad_helpers.cc:1682
Differentiation::AD::CellLevelBase::ad_type
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
Definition: ad_helpers.h:855
DataOutBase::none
@ none
Definition: data_out_base.h:1557
Differentiation::AD::is_taped_ad_number
Definition: ad_number_traits.h:260
Differentiation::AD::PointLevelFunctionsBase
Definition: ad_helpers.h:2640
Differentiation::AD::NumberTypes
NumberTypes
Definition: ad_number_types.h:36
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
Differentiation::AD::VectorFunction::scalar_type
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:3489
Tensor
Definition: tensor.h:450
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
ad_drivers.h
Differentiation
Definition: numbers.h:645
Differentiation::AD::HelperBase::register_dependent_variable
void register_dependent_variable(const unsigned int index, const ad_type &func)
Definition: ad_helpers.cc:684
Differentiation::AD::ResidualLinearization::compute_residual
virtual void compute_residual(Vector< scalar_type > &residual) const override
Definition: ad_helpers.cc:1037
Differentiation::AD::ResidualLinearization::compute_linearization
virtual void compute_linearization(FullMatrix< scalar_type > &linearization) const override
Definition: ad_helpers.cc:1091
Differentiation::AD::EnergyFunctional::scalar_type
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:1232
Differentiation::AD::ResidualLinearization::ResidualLinearization
ResidualLinearization(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:1011
Differentiation::AD::HelperBase< ADNumberTypeCode, double >::scalar_type
typename AD::NumberTraits< double, ADNumberTypeCode >::scalar_type scalar_type
Definition: ad_helpers.h:180
Differentiation::AD::EnergyFunctional::ad_type
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
Definition: ad_helpers.h:1239
unsigned int
Differentiation::AD::Types
Definition: ad_drivers.h:94
value
static const bool value
Definition: dof_tools_constraints.cc:433
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
Differentiation::AD::ScalarFunction::register_dependent_variable
void register_dependent_variable(const ad_type &func)
Definition: ad_helpers.cc:1340
Differentiation::AD::VectorFunction::VectorFunction
VectorFunction(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:1651
Differentiation::AD::ScalarFunction::ad_type
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
Definition: ad_helpers.h:3103
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
Differentiation::AD::CellLevelBase
Definition: ad_helpers.h:840
Differentiation::AD::ADNumberTraits
Definition: ad_number_traits.h:72
Differentiation::AD::ScalarFunction::extract_hessian_component
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)
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:191
Differentiation::AD::EnergyFunctional::compute_linearization
virtual void compute_linearization(FullMatrix< scalar_type > &linearization) const override
Definition: ad_helpers.cc:938
Differentiation::AD::VectorFunction::register_dependent_variables
void register_dependent_variables(const std::vector< ad_type > &funcs)
Definition: ad_helpers.cc:1666
Differentiation::AD::HelperBase::is_registered_tape
bool is_registered_tape(const typename Types< ad_type >::tape_index tape_index) const
Definition: ad_helpers.cc:296
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
ad_type
config.h
Differentiation::AD::TapelessDrivers
Definition: ad_drivers.h:522
FullMatrix
Definition: full_matrix.h:71
internal
Definition: aligned_vector.h:369
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
Differentiation::AD::ScalarFunction::compute_value
scalar_type compute_value() const
Definition: ad_helpers.cc:1353
Differentiation::AD::internal::configure_tapeless_mode
std::enable_if<!(ADNumberTraits< ADNumberType >::type_code==NumberTypes::adolc_tapeless)>::type configure_tapeless_mode(const unsigned int)
Definition: ad_drivers.cc:1630
Differentiation::AD::ScalarFunction::compute_gradient
void compute_gradient(Vector< scalar_type > &gradient) const
Definition: ad_helpers.cc:1403
Differentiation::AD::internal::Marking
Definition: ad_number_traits.h:151
Differentiation::AD::EnergyFunctional
Definition: ad_helpers.h:1224
Differentiation::AD::Numbers
Definition: ad_drivers.h:113
Differentiation::AD::ResidualLinearization::register_residual_vector
void register_residual_vector(const std::vector< ad_type > &residual)
Definition: ad_helpers.cc:1023
Differentiation::AD::HelperBase::is_recording
bool is_recording() const
Definition: ad_helpers.cc:271
Differentiation::AD::HelperBase::HelperBase
HelperBase(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:38