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