Reference documentation for deal.II version GIT d7aca55de5 2022-08-10 12:50:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
sacado_number_types.h
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.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_differentiation_ad_sacado_number_types_h
17 #define dealii_differentiation_ad_sacado_number_types_h
18 
19 #include <deal.II/base/config.h>
20 
21 #include <type_traits>
22 
23 
25 
26 
27 namespace Differentiation
28 {
29  namespace AD
30  {
36  template <typename NumberType, typename = void>
37  struct is_sacado_number : std::false_type
38  {};
39 
40 
46  template <typename NumberType, typename = void>
47  struct is_sacado_dfad_number : std::false_type
48  {};
49 
50 
56  template <typename NumberType, typename = void>
57  struct is_sacado_rad_number : std::false_type
58  {};
59 
60  } // namespace AD
61 } // namespace Differentiation
62 
63 
65 
66 
67 
68 #ifdef DEAL_II_TRILINOS_WITH_SACADO
69 
71 # include <Sacado.hpp>
72 // It appears that some versions of Trilinos do not directly or indirectly
73 // include all the headers for all forward and reverse Sacado AD types.
74 // So we directly include these both here as a precaution.
75 // Standard forward AD classes (templated)
76 # include <Sacado_Fad_DFad.hpp>
77 // Reverse AD classes (templated)
78 # include <Sacado_trad.hpp>
80 
81 # include <deal.II/base/exceptions.h>
82 # include <deal.II/base/numbers.h>
83 
86 
87 # include <complex>
88 
90 
91 
92 namespace Differentiation
93 {
94  namespace AD
95  {
96  namespace internal
97  {
104  template <typename SacadoNumber, typename = void>
106 
107  } // namespace internal
108 
109 
110 
111  } // namespace AD
112 } // namespace Differentiation
113 
114 
115 /* ----------- inline and template functions and specializations ----------- */
116 
117 
118 # ifndef DOXYGEN
119 
120 namespace Differentiation
121 {
122  namespace AD
123  {
124  namespace internal
125  {
126  // The documentation on Sacado numbers is pretty sparse and/or hard to
127  // navigate. As a point of reference, see
128  // https://trilinos.org/docs/dev/packages/sacado/doc/html/classSacado_1_1Fad_1_1SimpleFad.html
129  // for semi-applicable documentation for the Sacado::Fad::Dfad class.
130  // and the examples in
131  // https://github.com/trilinos/Trilinos/tree/master/packages/sacado/example
132  //
133  // If one dares to venture there, the relevant files for the classes
134  // supported here are:
135  //
136  // Forward-mode auto-differentiable types:
137  // https://github.com/trilinos/Trilinos/blob/master/packages/sacado/src/sacado_dfad_DFad.hpp
138  // https://github.com/trilinos/Trilinos/blob/master/packages/sacado/src/sacado_dfad_GeneralFad.hpp
139  //
140  // Reverse-mode auto-differentiable types:
141  // https://github.com/trilinos/Trilinos/blob/master/packages/sacado/src/Sacado_trad.hpp
142 
143 
147  template <typename SacadoNumber>
148  struct SacadoNumberInfo<
149  SacadoNumber,
150  std::enable_if_t<std::is_same<
151  SacadoNumber,
152  Sacado::Fad::DFad<typename SacadoNumber::value_type>>::value>>
153  {
154  using ad_type = SacadoNumber;
155  using scalar_type = typename ad_type::scalar_type;
156  using value_type = typename ad_type::value_type;
157  using derivative_type = typename ad_type::value_type;
158 
159  static const unsigned int n_supported_derivative_levels =
160  1 + SacadoNumberInfo<derivative_type>::n_supported_derivative_levels;
161  };
162 
163 
167  template <typename SacadoNumber>
168  struct SacadoNumberInfo<
169  SacadoNumber,
170  std::enable_if_t<std::is_same<
171  SacadoNumber,
172  Sacado::Rad::ADvar<typename SacadoNumber::value_type>>::value>>
173  {
174  using ad_type = SacadoNumber;
175  using scalar_type = typename ad_type::ADVari::scalar_type;
176  using value_type = typename ad_type::ADVari::value_type;
177  using derivative_type = typename ad_type::ADVari::value_type;
178 
179  static const unsigned int n_supported_derivative_levels =
180  1 + SacadoNumberInfo<derivative_type>::n_supported_derivative_levels;
181  };
182 
183 
190  template <typename Number>
191  struct SacadoNumberInfo<Number,
192  std::enable_if_t<std::is_arithmetic<
193  typename std::decay<Number>::type>::value>>
194  {
195  static const unsigned int n_supported_derivative_levels = 0;
196  };
197 
198 
203  template <typename ScalarType>
204  struct ADNumberInfoFromEnum<
205  ScalarType,
207  std::enable_if_t<std::is_floating_point<ScalarType>::value>>
208  {
209  static const bool is_taped = false;
210  using real_type = Sacado::Fad::DFad<ScalarType>;
211  using derivative_type =
212  typename SacadoNumberInfo<real_type>::derivative_type;
213  static const unsigned int n_supported_derivative_levels =
214  SacadoNumberInfo<real_type>::n_supported_derivative_levels;
215  };
216 
217 
222  template <typename ScalarType>
223  struct ADNumberInfoFromEnum<
224  ScalarType,
226  std::enable_if_t<std::is_floating_point<ScalarType>::value>>
227  {
228  static const bool is_taped = false;
229  using real_type = Sacado::Fad::DFad<Sacado::Fad::DFad<ScalarType>>;
230  using derivative_type =
231  typename SacadoNumberInfo<real_type>::derivative_type;
232  static const unsigned int n_supported_derivative_levels =
233  SacadoNumberInfo<real_type>::n_supported_derivative_levels;
234  };
235 
236 
241  template <typename ScalarType>
242  struct ADNumberInfoFromEnum<
243  ScalarType,
245  std::enable_if_t<std::is_floating_point<ScalarType>::value>>
246  {
247  static const bool is_taped = false;
248  using real_type = Sacado::Rad::ADvar<ScalarType>;
249  using derivative_type =
250  typename SacadoNumberInfo<real_type>::derivative_type;
251  static const unsigned int n_supported_derivative_levels =
252  SacadoNumberInfo<real_type>::n_supported_derivative_levels;
253  };
254 
255 
260  template <typename ScalarType>
261  struct ADNumberInfoFromEnum<
262  ScalarType,
264  std::enable_if_t<std::is_floating_point<ScalarType>::value>>
265  {
266  static const bool is_taped = false;
267  using real_type = Sacado::Rad::ADvar<Sacado::Fad::DFad<ScalarType>>;
268  using derivative_type =
269  typename SacadoNumberInfo<real_type>::derivative_type;
270  static const unsigned int n_supported_derivative_levels =
271  SacadoNumberInfo<real_type>::n_supported_derivative_levels;
272  };
273 
274 
279  template <typename NumberType>
280  struct Marking<Sacado::Fad::DFad<NumberType>>
281  {
282  using ad_type =
283  typename SacadoNumberInfo<Sacado::Fad::DFad<NumberType>>::ad_type;
284  using derivative_type = typename SacadoNumberInfo<
285  Sacado::Fad::DFad<NumberType>>::derivative_type;
286  using scalar_type =
287  typename SacadoNumberInfo<Sacado::Fad::DFad<NumberType>>::scalar_type;
288 
289  /*
290  * Initialize the state of an independent variable.
291  */
292  static void
293  independent_variable(const scalar_type &in,
294  const unsigned int index,
295  const unsigned int n_independent_variables,
296  ad_type & out)
297  {
298  // It is required that we first initialise the outer number before
299  // any of the nested ones.
300  out = ad_type(n_independent_variables, index, in);
301 
302  // Initialize potential nested directional derivatives
303  Marking<derivative_type>::independent_variable(
304  in, index, n_independent_variables, out.val());
305  }
306 
307  /*
308  * Initialize the state of a dependent variable.
309  */
310  static void
311  dependent_variable(ad_type &out, const ad_type &func)
312  {
313  out = func;
314  }
315  };
316 
317 
322  template <typename NumberType>
323  struct Marking<Sacado::Rad::ADvar<NumberType>>
324  {
325  using ad_type =
326  typename SacadoNumberInfo<Sacado::Rad::ADvar<NumberType>>::ad_type;
327  using derivative_type = typename SacadoNumberInfo<
328  Sacado::Rad::ADvar<NumberType>>::derivative_type;
329  using scalar_type = typename SacadoNumberInfo<
330  Sacado::Rad::ADvar<NumberType>>::scalar_type;
331 
332  /*
333  * Initialize the state of an independent variable.
334  */
335  static void
336  independent_variable(const scalar_type &in,
337  const unsigned int index,
338  const unsigned int n_independent_variables,
339  ad_type & out)
340  {
341  // For Sacado::Rad::ADvar numbers, we have to initialize the
342  // ADNumber with an already fully-configured value. This means
343  // that if this nests another ADNumber then the nested number
344  // must already be setup and ready for use.
345 
346  // Initialize potential nested directional derivatives
347  derivative_type derivative_initializer;
348  Marking<derivative_type>::independent_variable(
349  in, index, n_independent_variables, derivative_initializer);
350 
351  // Initialize the outer ad_type
352  out = derivative_initializer;
353  }
354 
355  /*
356  * Initialize the state of a dependent variable.
357  */
358  static void
359  dependent_variable(ad_type &out, const ad_type &func)
360  {
361  out = func;
362  }
363  };
364 
365 
371  template <typename NumberType>
372  struct ExtractData<Sacado::Fad::DFad<NumberType>>
373  {
374  using derivative_type = typename SacadoNumberInfo<
375  Sacado::Fad::DFad<NumberType>>::derivative_type;
376  using scalar_type =
377  typename SacadoNumberInfo<Sacado::Fad::DFad<NumberType>>::scalar_type;
378  using value_type =
379  typename SacadoNumberInfo<Sacado::Fad::DFad<NumberType>>::value_type;
380 
384  static scalar_type
385  value(const Sacado::Fad::DFad<NumberType> &x)
386  {
387  return ExtractData<value_type>::value(x.val());
388  }
389 
390 
394  static unsigned int
395  n_directional_derivatives(const Sacado::Fad::DFad<NumberType> &x)
396  {
397  return x.size();
398  }
399 
400 
404  static derivative_type
405  directional_derivative(const Sacado::Fad::DFad<NumberType> &x,
406  const unsigned int direction)
407  {
408  if (x.hasFastAccess())
409  return x.fastAccessDx(direction);
410  else
411  return x.dx(direction);
412  }
413  };
414 
415 
421  template <typename NumberType>
422  struct ExtractData<Sacado::Rad::ADvar<NumberType>>
423  {
424  using derivative_type = typename SacadoNumberInfo<
425  Sacado::Rad::ADvar<NumberType>>::derivative_type;
426  using scalar_type = typename SacadoNumberInfo<
427  Sacado::Rad::ADvar<NumberType>>::scalar_type;
428  using value_type =
429  typename SacadoNumberInfo<Sacado::Rad::ADvar<NumberType>>::value_type;
430 
434  static scalar_type
435  value(const Sacado::Rad::ADvar<NumberType> &x)
436  {
437  return ExtractData<value_type>::value(x.val());
438  }
439 
440 
444  static unsigned int
445  n_directional_derivatives(const Sacado::Rad::ADvar<NumberType> &)
446  {
447  // There are as many directional derivatives as there are
448  // independent variables, but each independent variable can
449  // only return one directional derivative.
450  return 1;
451  }
452 
453 
461  static derivative_type
462  directional_derivative(const Sacado::Rad::ADvar<NumberType> &x,
463  const unsigned int)
464  {
465  return x.adj();
466  }
467  };
468 
469  } // namespace internal
470 
471 
472  /* -------------- NumberTypes::sacado_dfad -------------- */
473 
474 
481  template <typename ADNumberType>
482  struct ADNumberTraits<
483  ADNumberType,
484  std::enable_if_t<std::is_same<
485  ADNumberType,
486  Sacado::Fad::DFad<typename ADNumberType::scalar_type>>::value>>
487  : NumberTraits<typename ADNumberType::scalar_type,
488  NumberTypes::sacado_dfad>
489  {};
490 
491 
498  template <typename ADNumberType>
499  struct ADNumberTraits<
500  ADNumberType,
501  std::enable_if_t<std::is_same<
502  ADNumberType,
503  std::complex<Sacado::Fad::DFad<
504  typename ADNumberType::value_type::scalar_type>>>::value>>
505  : NumberTraits<
506  std::complex<typename ADNumberType::value_type::scalar_type>,
507  NumberTypes::sacado_dfad>
508  {};
509 
510 
515  template <>
516  struct NumberTraits<Sacado::Fad::DFad<float>, NumberTypes::sacado_dfad>
517  : NumberTraits<
518  typename ADNumberTraits<Sacado::Fad::DFad<float>>::scalar_type,
519  NumberTypes::sacado_dfad>
520  {};
521 
522 
527  template <>
528  struct NumberTraits<std::complex<Sacado::Fad::DFad<float>>,
530  : NumberTraits<typename ADNumberTraits<
531  std::complex<Sacado::Fad::DFad<float>>>::scalar_type,
532  NumberTypes::sacado_dfad>
533  {};
534 
535 
540  template <>
541  struct NumberTraits<Sacado::Fad::DFad<double>, NumberTypes::sacado_dfad>
542  : NumberTraits<
543  typename ADNumberTraits<Sacado::Fad::DFad<double>>::scalar_type,
544  NumberTypes::sacado_dfad>
545  {};
546 
547 
552  template <>
553  struct NumberTraits<std::complex<Sacado::Fad::DFad<double>>,
555  : NumberTraits<typename ADNumberTraits<
556  std::complex<Sacado::Fad::DFad<double>>>::scalar_type,
557  NumberTypes::sacado_dfad>
558  {};
559 
560 
561  /* -------------- NumberTypes::sacado_rad -------------- */
562 
563 
570  template <typename ADNumberType>
571  struct ADNumberTraits<
572  ADNumberType,
573  std::enable_if_t<std::is_same<
574  ADNumberType,
575  Sacado::Rad::ADvar<typename ADNumberType::ADVari::scalar_type>>::value>>
576  : NumberTraits<typename ADNumberType::ADVari::scalar_type,
577  NumberTypes::sacado_rad>
578  {};
579 
580 
585  template <>
586  struct NumberTraits<Sacado::Rad::ADvar<float>, NumberTypes::sacado_rad>
587  : NumberTraits<
588  typename ADNumberTraits<Sacado::Rad::ADvar<float>>::scalar_type,
589  NumberTypes::sacado_rad>
590  {};
591 
592 
597  template <>
598  struct NumberTraits<Sacado::Rad::ADvar<double>, NumberTypes::sacado_rad>
599  : NumberTraits<
600  typename ADNumberTraits<Sacado::Rad::ADvar<double>>::scalar_type,
601  NumberTypes::sacado_rad>
602  {};
603 
604 
605 # ifdef DEAL_II_TRILINOS_CXX_SUPPORTS_SACADO_COMPLEX_RAD
606 
607 
614  template <typename ADNumberType>
615  struct ADNumberTraits<
616  ADNumberType,
617  std::enable_if_t<std::is_same<
618  ADNumberType,
619  std::complex<Sacado::Rad::ADvar<
620  typename ADNumberType::value_type::ADVari::scalar_type>>>::value>>
621  : NumberTraits<
622  std::complex<typename ADNumberType::value_type::ADVari::scalar_type>,
623  NumberTypes::sacado_rad>
624  {};
625 
626 
631  template <>
632  struct NumberTraits<std::complex<Sacado::Rad::ADvar<float>>,
634  : NumberTraits<typename ADNumberTraits<
635  std::complex<Sacado::Rad::ADvar<float>>>::scalar_type,
636  NumberTypes::sacado_rad>
637  {};
638 
639 
644  template <>
645  struct NumberTraits<std::complex<Sacado::Rad::ADvar<double>>,
647  : NumberTraits<typename ADNumberTraits<
648  std::complex<Sacado::Rad::ADvar<double>>>::scalar_type,
649  NumberTypes::sacado_rad>
650  {};
651 
652 
653 # endif
654 
655 
656  /* -------------- NumberTypes::sacado_dfad_dfad -------------- */
657 
665  template <typename ADNumberType>
666  struct ADNumberTraits<ADNumberType,
667  std::enable_if_t<std::is_same<
668  ADNumberType,
669  Sacado::Fad::DFad<Sacado::Fad::DFad<
670  typename ADNumberType::scalar_type>>>::value>>
671  : NumberTraits<typename ADNumberType::scalar_type,
672  NumberTypes::sacado_dfad_dfad>
673  {};
674 
675 
683  template <typename ADNumberType>
684  struct ADNumberTraits<
685  ADNumberType,
686  std::enable_if_t<std::is_same<
687  ADNumberType,
688  std::complex<Sacado::Fad::DFad<Sacado::Fad::DFad<
689  typename ADNumberType::value_type::scalar_type>>>>::value>>
690  : NumberTraits<
691  std::complex<typename ADNumberType::value_type::scalar_type>,
692  NumberTypes::sacado_dfad_dfad>
693  {};
694 
695 
700  template <>
701  struct NumberTraits<Sacado::Fad::DFad<Sacado::Fad::DFad<float>>,
703  : NumberTraits<typename ADNumberTraits<Sacado::Fad::DFad<
704  Sacado::Fad::DFad<float>>>::scalar_type,
705  NumberTypes::sacado_dfad_dfad>
706  {};
707 
708 
713  template <>
714  struct NumberTraits<
715  std::complex<Sacado::Fad::DFad<Sacado::Fad::DFad<float>>>,
717  : NumberTraits<typename ADNumberTraits<std::complex<Sacado::Fad::DFad<
718  Sacado::Fad::DFad<float>>>>::scalar_type,
719  NumberTypes::sacado_dfad_dfad>
720  {};
721 
722 
727  template <>
728  struct NumberTraits<Sacado::Fad::DFad<Sacado::Fad::DFad<double>>,
730  : NumberTraits<typename ADNumberTraits<Sacado::Fad::DFad<
731  Sacado::Fad::DFad<double>>>::scalar_type,
732  NumberTypes::sacado_dfad_dfad>
733  {};
734 
735 
740  template <>
741  struct NumberTraits<
742  std::complex<Sacado::Fad::DFad<Sacado::Fad::DFad<double>>>,
744  : NumberTraits<typename ADNumberTraits<std::complex<Sacado::Fad::DFad<
745  Sacado::Fad::DFad<double>>>>::scalar_type,
746  NumberTypes::sacado_dfad_dfad>
747  {};
748 
749 
750  /* -------------- NumberTypes::sacado_rad_dfad -------------- */
751 
759  template <typename ADNumberType>
760  struct ADNumberTraits<
761  ADNumberType,
762  std::enable_if_t<
763  std::is_same<ADNumberType,
764  Sacado::Rad::ADvar<Sacado::Fad::DFad<
765  typename ADNumberType::ADVari::scalar_type>>>::value>>
766  : NumberTraits<typename ADNumberType::ADVari::scalar_type,
767  NumberTypes::sacado_rad_dfad>
768  {};
769 
770 
775  template <>
776  struct NumberTraits<Sacado::Rad::ADvar<Sacado::Fad::DFad<float>>,
778  : NumberTraits<typename ADNumberTraits<Sacado::Rad::ADvar<
779  Sacado::Fad::DFad<float>>>::scalar_type,
780  NumberTypes::sacado_rad_dfad>
781  {};
782 
783 
788  template <>
789  struct NumberTraits<Sacado::Rad::ADvar<Sacado::Fad::DFad<double>>,
791  : NumberTraits<typename ADNumberTraits<Sacado::Rad::ADvar<
792  Sacado::Fad::DFad<double>>>::scalar_type,
793  NumberTypes::sacado_rad_dfad>
794  {};
795 
796 
797 # ifdef DEAL_II_TRILINOS_CXX_SUPPORTS_SACADO_COMPLEX_RAD
798 
799 
807  template <typename ADNumberType>
808  struct ADNumberTraits<
809  ADNumberType,
810  std::enable_if_t<std::is_same<
811  ADNumberType,
812  std::complex<Sacado::Rad::ADvar<Sacado::Fad::DFad<
813  typename ADNumberType::value_type::ADVari::scalar_type>>>>::value>>
814  : NumberTraits<
815  std::complex<typename ADNumberType::value_type::ADVari::scalar_type>,
816  NumberTypes::sacado_rad_dfad>
817  {};
818 
819 
824  template <>
825  struct NumberTraits<
826  std::complex<Sacado::Rad::ADvar<Sacado::Fad::DFad<float>>>,
828  : NumberTraits<typename ADNumberTraits<std::complex<Sacado::Rad::ADvar<
829  Sacado::Fad::DFad<float>>>>::scalar_type,
830  NumberTypes::sacado_rad_dfad>
831  {};
832 
833 
838  template <>
839  struct NumberTraits<
840  std::complex<Sacado::Rad::ADvar<Sacado::Fad::DFad<double>>>,
842  : NumberTraits<typename ADNumberTraits<std::complex<Sacado::Rad::ADvar<
843  Sacado::Fad::DFad<double>>>>::scalar_type,
844  NumberTypes::sacado_rad_dfad>
845  {};
846 
847 
848 # endif
849 
850 
851  /* -------------- Additional type traits -------------- */
852 
853 
854  template <typename NumberType>
855  struct is_sacado_dfad_number<
856  NumberType,
857  std::enable_if_t<
858  ADNumberTraits<typename std::decay<NumberType>::type>::type_code ==
859  NumberTypes::sacado_dfad ||
860  ADNumberTraits<typename std::decay<NumberType>::type>::type_code ==
861  NumberTypes::sacado_dfad_dfad>> : std::true_type
862  {};
863 
864 
865  template <typename NumberType>
866  struct is_sacado_dfad_number<
867  NumberType,
868  std::enable_if_t<std::is_same<
869  NumberType,
870  Sacado::Fad::Expr<typename NumberType::value_type>>::value>>
871  : std::true_type
872  {};
873 
874 
875  template <typename NumberType>
876  struct is_sacado_rad_number<
877  NumberType,
878  std::enable_if_t<
879  ADNumberTraits<typename std::decay<NumberType>::type>::type_code ==
880  NumberTypes::sacado_rad ||
881  ADNumberTraits<typename std::decay<NumberType>::type>::type_code ==
882  NumberTypes::sacado_rad_dfad>> : std::true_type
883  {};
884 
885 
886  template <typename NumberType>
887  struct is_sacado_rad_number<
888  NumberType,
889  std::enable_if_t<
890  std::is_same<NumberType,
891  Sacado::Rad::ADvari<Sacado::Fad::DFad<
892  typename NumberType::ADVari::scalar_type>>>::value>>
893  : std::true_type
894  {};
895 
896 
897  template <typename NumberType>
898  struct is_sacado_number<
899  NumberType,
900  std::enable_if_t<is_sacado_dfad_number<NumberType>::value ||
901  is_sacado_rad_number<NumberType>::value>>
902  : std::true_type
903  {};
904 
905  } // namespace AD
906 } // namespace Differentiation
907 
908 
909 namespace numbers
910 {
911  template <typename NumberType>
912  struct is_cuda_compatible<
913  NumberType,
914  std::enable_if_t<
915  ::Differentiation::AD::is_sacado_rad_number<NumberType>::value>>
916  : std::false_type
917  {};
918 } // namespace numbers
919 
920 # endif // DOXYGEN
921 
922 
923 
925 
926 
927 #endif // DEAL_II_TRILINOS_WITH_SACADO
928 
929 #endif
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:456
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:495