Reference documentation for deal.II version 9.3.3
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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
27namespace 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
82# include <deal.II/base/numbers.h>
83
86
87# include <complex>
88
90
91
92namespace 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
120namespace 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 typename std::enable_if<std::is_same<
151 SacadoNumber,
152 Sacado::Fad::DFad<typename SacadoNumber::value_type>>::value>::type>
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 typename std::enable_if<std::is_same<
171 SacadoNumber,
172 Sacado::Rad::ADvar<typename SacadoNumber::value_type>>::value>::type>
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<
192 Number,
193 typename std::enable_if<
194 std::is_arithmetic<typename std::decay<Number>::type>::value>::type>
195 {
196 static const unsigned int n_supported_derivative_levels = 0;
197 };
198
199
204 template <typename ScalarType>
205 struct ADNumberInfoFromEnum<
206 ScalarType,
208 typename std::enable_if<
209 std::is_floating_point<ScalarType>::value>::type>
210 {
211 static const bool is_taped = false;
212 using real_type = Sacado::Fad::DFad<ScalarType>;
213 using derivative_type =
214 typename SacadoNumberInfo<real_type>::derivative_type;
215 static const unsigned int n_supported_derivative_levels =
216 SacadoNumberInfo<real_type>::n_supported_derivative_levels;
217 };
218
219
224 template <typename ScalarType>
225 struct ADNumberInfoFromEnum<
226 ScalarType,
228 typename std::enable_if<
229 std::is_floating_point<ScalarType>::value>::type>
230 {
231 static const bool is_taped = false;
232 using real_type = Sacado::Fad::DFad<Sacado::Fad::DFad<ScalarType>>;
233 using derivative_type =
234 typename SacadoNumberInfo<real_type>::derivative_type;
235 static const unsigned int n_supported_derivative_levels =
236 SacadoNumberInfo<real_type>::n_supported_derivative_levels;
237 };
238
239
244 template <typename ScalarType>
245 struct ADNumberInfoFromEnum<
246 ScalarType,
248 typename std::enable_if<
249 std::is_floating_point<ScalarType>::value>::type>
250 {
251 static const bool is_taped = false;
252 using real_type = Sacado::Rad::ADvar<ScalarType>;
253 using derivative_type =
254 typename SacadoNumberInfo<real_type>::derivative_type;
255 static const unsigned int n_supported_derivative_levels =
256 SacadoNumberInfo<real_type>::n_supported_derivative_levels;
257 };
258
259
264 template <typename ScalarType>
265 struct ADNumberInfoFromEnum<
266 ScalarType,
268 typename std::enable_if<
269 std::is_floating_point<ScalarType>::value>::type>
270 {
271 static const bool is_taped = false;
272 using real_type = Sacado::Rad::ADvar<Sacado::Fad::DFad<ScalarType>>;
273 using derivative_type =
274 typename SacadoNumberInfo<real_type>::derivative_type;
275 static const unsigned int n_supported_derivative_levels =
276 SacadoNumberInfo<real_type>::n_supported_derivative_levels;
277 };
278
279
284 template <typename NumberType>
285 struct Marking<Sacado::Fad::DFad<NumberType>>
286 {
287 using ad_type =
288 typename SacadoNumberInfo<Sacado::Fad::DFad<NumberType>>::ad_type;
289 using derivative_type = typename SacadoNumberInfo<
290 Sacado::Fad::DFad<NumberType>>::derivative_type;
291 using scalar_type =
292 typename SacadoNumberInfo<Sacado::Fad::DFad<NumberType>>::scalar_type;
293
294 /*
295 * Initialize the state of an independent variable.
296 */
297 static void
298 independent_variable(const scalar_type &in,
299 const unsigned int index,
300 const unsigned int n_independent_variables,
301 ad_type & out)
302 {
303 // It is required that we first initialise the outer number before
304 // any of the nested ones.
305 out = ad_type(n_independent_variables, index, in);
306
307 // Initialize potential nested directional derivatives
308 Marking<derivative_type>::independent_variable(
309 in, index, n_independent_variables, out.val());
310 }
311
312 /*
313 * Initialize the state of a dependent variable.
314 */
315 static void
316 dependent_variable(ad_type &out, const ad_type &func)
317 {
318 out = func;
319 }
320 };
321
322
327 template <typename NumberType>
328 struct Marking<Sacado::Rad::ADvar<NumberType>>
329 {
330 using ad_type =
331 typename SacadoNumberInfo<Sacado::Rad::ADvar<NumberType>>::ad_type;
332 using derivative_type = typename SacadoNumberInfo<
333 Sacado::Rad::ADvar<NumberType>>::derivative_type;
334 using scalar_type = typename SacadoNumberInfo<
335 Sacado::Rad::ADvar<NumberType>>::scalar_type;
336
337 /*
338 * Initialize the state of an independent variable.
339 */
340 static void
341 independent_variable(const scalar_type &in,
342 const unsigned int index,
343 const unsigned int n_independent_variables,
344 ad_type & out)
345 {
346 // For Sacado::Rad::ADvar numbers, we have to initialize the
347 // ADNumber with an already fully-configured value. This means
348 // that if this nests another ADNumber then the nested number
349 // must already be setup and ready for use.
350
351 // Initialize potential nested directional derivatives
352 derivative_type derivative_initializer;
353 Marking<derivative_type>::independent_variable(
354 in, index, n_independent_variables, derivative_initializer);
355
356 // Initialize the outer ad_type
357 out = derivative_initializer;
358 }
359
360 /*
361 * Initialize the state of a dependent variable.
362 */
363 static void
364 dependent_variable(ad_type &out, const ad_type &func)
365 {
366 out = func;
367 }
368 };
369
370
376 template <typename NumberType>
377 struct ExtractData<Sacado::Fad::DFad<NumberType>>
378 {
379 using derivative_type = typename SacadoNumberInfo<
380 Sacado::Fad::DFad<NumberType>>::derivative_type;
381 using scalar_type =
382 typename SacadoNumberInfo<Sacado::Fad::DFad<NumberType>>::scalar_type;
383 using value_type =
384 typename SacadoNumberInfo<Sacado::Fad::DFad<NumberType>>::value_type;
385
389 static scalar_type
390 value(const Sacado::Fad::DFad<NumberType> &x)
391 {
392 return ExtractData<value_type>::value(x.val());
393 }
394
395
399 static unsigned int
400 n_directional_derivatives(const Sacado::Fad::DFad<NumberType> &x)
401 {
402 return x.size();
403 }
404
405
409 static derivative_type
410 directional_derivative(const Sacado::Fad::DFad<NumberType> &x,
411 const unsigned int direction)
412 {
413 if (x.hasFastAccess())
414 return x.fastAccessDx(direction);
415 else
416 return x.dx(direction);
417 }
418 };
419
420
426 template <typename NumberType>
427 struct ExtractData<Sacado::Rad::ADvar<NumberType>>
428 {
429 using derivative_type = typename SacadoNumberInfo<
430 Sacado::Rad::ADvar<NumberType>>::derivative_type;
431 using scalar_type = typename SacadoNumberInfo<
432 Sacado::Rad::ADvar<NumberType>>::scalar_type;
433 using value_type =
434 typename SacadoNumberInfo<Sacado::Rad::ADvar<NumberType>>::value_type;
435
439 static scalar_type
440 value(const Sacado::Rad::ADvar<NumberType> &x)
441 {
442 return ExtractData<value_type>::value(x.val());
443 }
444
445
449 static unsigned int
450 n_directional_derivatives(const Sacado::Rad::ADvar<NumberType> &)
451 {
452 // There are as many directional derivatives as there are
453 // independent variables, but each independent variable can
454 // only return one directional derivative.
455 return 1;
456 }
457
458
466 static derivative_type
467 directional_derivative(const Sacado::Rad::ADvar<NumberType> &x,
468 const unsigned int)
469 {
470 return x.adj();
471 }
472 };
473
474 } // namespace internal
475
476
477 /* -------------- NumberTypes::sacado_dfad -------------- */
478
479
486 template <typename ADNumberType>
487 struct ADNumberTraits<
488 ADNumberType,
489 typename std::enable_if<std::is_same<
490 ADNumberType,
491 Sacado::Fad::DFad<typename ADNumberType::scalar_type>>::value>::type>
492 : NumberTraits<typename ADNumberType::scalar_type,
493 NumberTypes::sacado_dfad>
494 {};
495
496
503 template <typename ADNumberType>
504 struct ADNumberTraits<
505 ADNumberType,
506 typename std::enable_if<std::is_same<
507 ADNumberType,
508 std::complex<Sacado::Fad::DFad<
509 typename ADNumberType::value_type::scalar_type>>>::value>::type>
510 : NumberTraits<
511 std::complex<typename ADNumberType::value_type::scalar_type>,
512 NumberTypes::sacado_dfad>
513 {};
514
515
520 template <>
521 struct NumberTraits<Sacado::Fad::DFad<float>, NumberTypes::sacado_dfad>
522 : NumberTraits<
523 typename ADNumberTraits<Sacado::Fad::DFad<float>>::scalar_type,
524 NumberTypes::sacado_dfad>
525 {};
526
527
532 template <>
533 struct NumberTraits<std::complex<Sacado::Fad::DFad<float>>,
535 : NumberTraits<typename ADNumberTraits<
536 std::complex<Sacado::Fad::DFad<float>>>::scalar_type,
537 NumberTypes::sacado_dfad>
538 {};
539
540
545 template <>
546 struct NumberTraits<Sacado::Fad::DFad<double>, NumberTypes::sacado_dfad>
547 : NumberTraits<
548 typename ADNumberTraits<Sacado::Fad::DFad<double>>::scalar_type,
549 NumberTypes::sacado_dfad>
550 {};
551
552
557 template <>
558 struct NumberTraits<std::complex<Sacado::Fad::DFad<double>>,
560 : NumberTraits<typename ADNumberTraits<
561 std::complex<Sacado::Fad::DFad<double>>>::scalar_type,
562 NumberTypes::sacado_dfad>
563 {};
564
565
566 /* -------------- NumberTypes::sacado_rad -------------- */
567
568
575 template <typename ADNumberType>
576 struct ADNumberTraits<
577 ADNumberType,
578 typename std::enable_if<std::is_same<
579 ADNumberType,
580 Sacado::Rad::ADvar<typename ADNumberType::ADVari::scalar_type>>::
581 value>::type>
582 : NumberTraits<typename ADNumberType::ADVari::scalar_type,
583 NumberTypes::sacado_rad>
584 {};
585
586
591 template <>
592 struct NumberTraits<Sacado::Rad::ADvar<float>, NumberTypes::sacado_rad>
593 : NumberTraits<
594 typename ADNumberTraits<Sacado::Rad::ADvar<float>>::scalar_type,
595 NumberTypes::sacado_rad>
596 {};
597
598
603 template <>
604 struct NumberTraits<Sacado::Rad::ADvar<double>, NumberTypes::sacado_rad>
605 : NumberTraits<
606 typename ADNumberTraits<Sacado::Rad::ADvar<double>>::scalar_type,
607 NumberTypes::sacado_rad>
608 {};
609
610
611# ifdef DEAL_II_TRILINOS_CXX_SUPPORTS_SACADO_COMPLEX_RAD
612
613
620 template <typename ADNumberType>
621 struct ADNumberTraits<
622 ADNumberType,
623 typename std::enable_if<std::is_same<
624 ADNumberType,
625 std::complex<Sacado::Rad::ADvar<typename ADNumberType::value_type::
626 ADVari::scalar_type>>>::value>::type>
627 : NumberTraits<
628 std::complex<typename ADNumberType::value_type::ADVari::scalar_type>,
629 NumberTypes::sacado_rad>
630 {};
631
632
637 template <>
638 struct NumberTraits<std::complex<Sacado::Rad::ADvar<float>>,
640 : NumberTraits<typename ADNumberTraits<
641 std::complex<Sacado::Rad::ADvar<float>>>::scalar_type,
642 NumberTypes::sacado_rad>
643 {};
644
645
650 template <>
651 struct NumberTraits<std::complex<Sacado::Rad::ADvar<double>>,
653 : NumberTraits<typename ADNumberTraits<
654 std::complex<Sacado::Rad::ADvar<double>>>::scalar_type,
655 NumberTypes::sacado_rad>
656 {};
657
658
659# endif
660
661
662 /* -------------- NumberTypes::sacado_dfad_dfad -------------- */
663
671 template <typename ADNumberType>
672 struct ADNumberTraits<
673 ADNumberType,
674 typename std::enable_if<
675 std::is_same<ADNumberType,
676 Sacado::Fad::DFad<Sacado::Fad::DFad<
677 typename ADNumberType::scalar_type>>>::value>::type>
678 : NumberTraits<typename ADNumberType::scalar_type,
679 NumberTypes::sacado_dfad_dfad>
680 {};
681
682
690 template <typename ADNumberType>
691 struct ADNumberTraits<
692 ADNumberType,
693 typename std::enable_if<std::is_same<
694 ADNumberType,
695 std::complex<Sacado::Fad::DFad<Sacado::Fad::DFad<
696 typename ADNumberType::value_type::scalar_type>>>>::value>::type>
697 : NumberTraits<
698 std::complex<typename ADNumberType::value_type::scalar_type>,
699 NumberTypes::sacado_dfad_dfad>
700 {};
701
702
707 template <>
708 struct NumberTraits<Sacado::Fad::DFad<Sacado::Fad::DFad<float>>,
710 : NumberTraits<typename ADNumberTraits<Sacado::Fad::DFad<
711 Sacado::Fad::DFad<float>>>::scalar_type,
712 NumberTypes::sacado_dfad_dfad>
713 {};
714
715
720 template <>
721 struct NumberTraits<
722 std::complex<Sacado::Fad::DFad<Sacado::Fad::DFad<float>>>,
724 : NumberTraits<typename ADNumberTraits<std::complex<Sacado::Fad::DFad<
725 Sacado::Fad::DFad<float>>>>::scalar_type,
726 NumberTypes::sacado_dfad_dfad>
727 {};
728
729
734 template <>
735 struct NumberTraits<Sacado::Fad::DFad<Sacado::Fad::DFad<double>>,
737 : NumberTraits<typename ADNumberTraits<Sacado::Fad::DFad<
738 Sacado::Fad::DFad<double>>>::scalar_type,
739 NumberTypes::sacado_dfad_dfad>
740 {};
741
742
747 template <>
748 struct NumberTraits<
749 std::complex<Sacado::Fad::DFad<Sacado::Fad::DFad<double>>>,
751 : NumberTraits<typename ADNumberTraits<std::complex<Sacado::Fad::DFad<
752 Sacado::Fad::DFad<double>>>>::scalar_type,
753 NumberTypes::sacado_dfad_dfad>
754 {};
755
756
757 /* -------------- NumberTypes::sacado_rad_dfad -------------- */
758
766 template <typename ADNumberType>
767 struct ADNumberTraits<
768 ADNumberType,
769 typename std::enable_if<std::is_same<
770 ADNumberType,
771 Sacado::Rad::ADvar<Sacado::Fad::DFad<
772 typename ADNumberType::ADVari::scalar_type>>>::value>::type>
773 : NumberTraits<typename ADNumberType::ADVari::scalar_type,
774 NumberTypes::sacado_rad_dfad>
775 {};
776
777
782 template <>
783 struct NumberTraits<Sacado::Rad::ADvar<Sacado::Fad::DFad<float>>,
785 : NumberTraits<typename ADNumberTraits<Sacado::Rad::ADvar<
786 Sacado::Fad::DFad<float>>>::scalar_type,
787 NumberTypes::sacado_rad_dfad>
788 {};
789
790
795 template <>
796 struct NumberTraits<Sacado::Rad::ADvar<Sacado::Fad::DFad<double>>,
798 : NumberTraits<typename ADNumberTraits<Sacado::Rad::ADvar<
799 Sacado::Fad::DFad<double>>>::scalar_type,
800 NumberTypes::sacado_rad_dfad>
801 {};
802
803
804# ifdef DEAL_II_TRILINOS_CXX_SUPPORTS_SACADO_COMPLEX_RAD
805
806
814 template <typename ADNumberType>
815 struct ADNumberTraits<
816 ADNumberType,
817 typename std::enable_if<std::is_same<
818 ADNumberType,
819 std::complex<Sacado::Rad::ADvar<Sacado::Fad::DFad<
820 typename ADNumberType::value_type::ADVari::scalar_type>>>>::value>::
821 type>
822 : NumberTraits<
823 std::complex<typename ADNumberType::value_type::ADVari::scalar_type>,
824 NumberTypes::sacado_rad_dfad>
825 {};
826
827
832 template <>
833 struct NumberTraits<
834 std::complex<Sacado::Rad::ADvar<Sacado::Fad::DFad<float>>>,
836 : NumberTraits<typename ADNumberTraits<std::complex<Sacado::Rad::ADvar<
837 Sacado::Fad::DFad<float>>>>::scalar_type,
838 NumberTypes::sacado_rad_dfad>
839 {};
840
841
846 template <>
847 struct NumberTraits<
848 std::complex<Sacado::Rad::ADvar<Sacado::Fad::DFad<double>>>,
850 : NumberTraits<typename ADNumberTraits<std::complex<Sacado::Rad::ADvar<
851 Sacado::Fad::DFad<double>>>>::scalar_type,
852 NumberTypes::sacado_rad_dfad>
853 {};
854
855
856# endif
857
858
859 /* -------------- Additional type traits -------------- */
860
861
862 template <typename NumberType>
863 struct is_sacado_dfad_number<
864 NumberType,
865 typename std::enable_if<
866 ADNumberTraits<typename std::decay<NumberType>::type>::type_code ==
867 NumberTypes::sacado_dfad ||
868 ADNumberTraits<typename std::decay<NumberType>::type>::type_code ==
869 NumberTypes::sacado_dfad_dfad>::type> : std::true_type
870 {};
871
872
873 template <typename NumberType>
874 struct is_sacado_dfad_number<
875 NumberType,
876 typename std::enable_if<std::is_same<
877 NumberType,
878 Sacado::Fad::Expr<typename NumberType::value_type>>::value>::type>
879 : std::true_type
880 {};
881
882
883 template <typename NumberType>
884 struct is_sacado_rad_number<
885 NumberType,
886 typename std::enable_if<
887 ADNumberTraits<typename std::decay<NumberType>::type>::type_code ==
888 NumberTypes::sacado_rad ||
889 ADNumberTraits<typename std::decay<NumberType>::type>::type_code ==
890 NumberTypes::sacado_rad_dfad>::type> : std::true_type
891 {};
892
893
894 template <typename NumberType>
895 struct is_sacado_rad_number<
896 NumberType,
897 typename std::enable_if<std::is_same<
898 NumberType,
899 Sacado::Rad::ADvari<Sacado::Fad::DFad<
900 typename NumberType::ADVari::scalar_type>>>::value>::type>
901 : std::true_type
902 {};
903
904
905 template <typename NumberType>
906 struct is_sacado_number<
907 NumberType,
908 typename std::enable_if<is_sacado_dfad_number<NumberType>::value ||
909 is_sacado_rad_number<NumberType>::value>::type>
910 : std::true_type
911 {};
912
913 } // namespace AD
914} // namespace Differentiation
915
916
917namespace numbers
918{
919 template <typename NumberType>
920 struct is_cuda_compatible<
921 NumberType,
922 typename std::enable_if<::Differentiation::AD::is_sacado_rad_number<
923 NumberType>::value>::type> : std::false_type
924 {};
925} // namespace numbers
926
927# endif // DOXYGEN
928
929
930
932
933
934#endif // DEAL_II_TRILINOS_WITH_SACADO
935
936#endif
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:416
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:454
STL namespace.