deal.II version GIT relicensing-3507-gc50d9d6ce3 2025-06-12 13:20:01+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\}}\)
Loading...
Searching...
No Matches
sacado_number_types.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2017 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_differentiation_ad_sacado_number_types_h
16#define dealii_differentiation_ad_sacado_number_types_h
17
18#include <deal.II/base/config.h>
19
22
25
26#ifdef DEAL_II_TRILINOS_WITH_SACADO
27// It appears that some versions of Trilinos do not directly or indirectly
28// include all the headers for all forward and reverse Sacado AD types
29// in Sacado.hpp, so we also directly include these here as a precaution:
30# include <Sacado.hpp>
31# include <Sacado_Fad_DFad.hpp>
32# include <Sacado_trad.hpp>
33
34# include <complex>
35#endif
36
37#include <type_traits>
38
39
41
42
43namespace Differentiation
44{
45 namespace AD
46 {
52 template <typename NumberType, typename = void>
53 struct is_sacado_number : std::false_type
54 {};
55
56
62 template <typename NumberType, typename = void>
63 struct is_sacado_dfad_number : std::false_type
64 {};
65
66
72 template <typename NumberType, typename = void>
73 struct is_sacado_rad_number : std::false_type
74 {};
75
76 } // namespace AD
77} // namespace Differentiation
78
79
80
81#ifdef DEAL_II_TRILINOS_WITH_SACADO
82
83
84namespace Differentiation
85{
86 namespace AD
87 {
88 namespace internal
89 {
96 template <typename SacadoNumber, typename = void>
98
99 } // namespace internal
100
101
102
103 } // namespace AD
104} // namespace Differentiation
105
106
107/* ----------- inline and template functions and specializations ----------- */
108
109
110# ifndef DOXYGEN
111
112namespace Differentiation
113{
114 namespace AD
115 {
116 namespace internal
117 {
118 // The documentation on Sacado numbers is pretty sparse and/or hard to
119 // navigate. As a point of reference, see
120 // https://trilinos.org/docs/dev/packages/sacado/doc/html/classSacado_1_1Fad_1_1SimpleFad.html
121 // for semi-applicable documentation for the Sacado::Fad::Dfad class.
122 // and the examples in
123 // https://github.com/trilinos/Trilinos/tree/master/packages/sacado/example
124 //
125 // If one dares to venture there, the relevant files for the classes
126 // supported here are:
127 //
128 // Forward-mode auto-differentiable types:
129 // https://github.com/trilinos/Trilinos/blob/master/packages/sacado/src/sacado_dfad_DFad.hpp
130 // https://github.com/trilinos/Trilinos/blob/master/packages/sacado/src/sacado_dfad_GeneralFad.hpp
131 //
132 // Reverse-mode auto-differentiable types:
133 // https://github.com/trilinos/Trilinos/blob/master/packages/sacado/src/Sacado_trad.hpp
134
135
139 template <typename SacadoNumber>
140 struct SacadoNumberInfo<
141 SacadoNumber,
142 std::enable_if_t<
143 std::is_same_v<SacadoNumber,
144 Sacado::Fad::DFad<typename SacadoNumber::value_type>>>>
145 {
146 using ad_type = SacadoNumber;
147 using scalar_type = typename ad_type::scalar_type;
148 using value_type = typename ad_type::value_type;
149 using derivative_type = typename ad_type::value_type;
150
151 static const unsigned int n_supported_derivative_levels =
152 1 + SacadoNumberInfo<derivative_type>::n_supported_derivative_levels;
153 };
154
155
159 template <typename SacadoNumber>
160 struct SacadoNumberInfo<
161 SacadoNumber,
162 std::enable_if_t<std::is_same_v<
163 SacadoNumber,
164 Sacado::Rad::ADvar<typename SacadoNumber::value_type>>>>
165 {
166 using ad_type = SacadoNumber;
167 using scalar_type = typename ad_type::ADVari::scalar_type;
168 using value_type = typename ad_type::ADVari::value_type;
169 using derivative_type = typename ad_type::ADVari::value_type;
170
171 static const unsigned int n_supported_derivative_levels =
172 1 + SacadoNumberInfo<derivative_type>::n_supported_derivative_levels;
173 };
174
175
182 template <typename Number>
183 struct SacadoNumberInfo<
184 Number,
185 std::enable_if_t<std::is_arithmetic_v<std::decay_t<Number>>>>
186 {
187 static const unsigned int n_supported_derivative_levels = 0;
188 };
189
190
195 template <typename ScalarType>
196 struct ADNumberInfoFromEnum<
197 ScalarType,
199 std::enable_if_t<std::is_floating_point_v<ScalarType>>>
200 {
201 static const bool is_taped = false;
202 using real_type = Sacado::Fad::DFad<ScalarType>;
203 using derivative_type =
204 typename SacadoNumberInfo<real_type>::derivative_type;
205 static const unsigned int n_supported_derivative_levels =
206 SacadoNumberInfo<real_type>::n_supported_derivative_levels;
207 };
208
209
214 template <typename ScalarType>
215 struct ADNumberInfoFromEnum<
216 ScalarType,
218 std::enable_if_t<std::is_floating_point_v<ScalarType>>>
219 {
220 static const bool is_taped = false;
221 using real_type = Sacado::Fad::DFad<Sacado::Fad::DFad<ScalarType>>;
222 using derivative_type =
223 typename SacadoNumberInfo<real_type>::derivative_type;
224 static const unsigned int n_supported_derivative_levels =
225 SacadoNumberInfo<real_type>::n_supported_derivative_levels;
226 };
227
228
233 template <typename ScalarType>
234 struct ADNumberInfoFromEnum<
235 ScalarType,
237 std::enable_if_t<std::is_floating_point_v<ScalarType>>>
238 {
239 static const bool is_taped = false;
240 using real_type = Sacado::Rad::ADvar<ScalarType>;
241 using derivative_type =
242 typename SacadoNumberInfo<real_type>::derivative_type;
243 static const unsigned int n_supported_derivative_levels =
244 SacadoNumberInfo<real_type>::n_supported_derivative_levels;
245 };
246
247
252 template <typename ScalarType>
253 struct ADNumberInfoFromEnum<
254 ScalarType,
256 std::enable_if_t<std::is_floating_point_v<ScalarType>>>
257 {
258 static const bool is_taped = false;
259 using real_type = Sacado::Rad::ADvar<Sacado::Fad::DFad<ScalarType>>;
260 using derivative_type =
261 typename SacadoNumberInfo<real_type>::derivative_type;
262 static const unsigned int n_supported_derivative_levels =
263 SacadoNumberInfo<real_type>::n_supported_derivative_levels;
264 };
265
266
271 template <typename NumberType>
272 struct Marking<Sacado::Fad::DFad<NumberType>>
273 {
274 using ad_type =
275 typename SacadoNumberInfo<Sacado::Fad::DFad<NumberType>>::ad_type;
276 using derivative_type = typename SacadoNumberInfo<
277 Sacado::Fad::DFad<NumberType>>::derivative_type;
278 using scalar_type =
279 typename SacadoNumberInfo<Sacado::Fad::DFad<NumberType>>::scalar_type;
280
281 /*
282 * Initialize the state of an independent variable.
283 */
284 static void
285 independent_variable(const scalar_type &in,
286 const unsigned int index,
287 const unsigned int n_independent_variables,
288 ad_type &out)
289 {
290 // It is required that we first initialize the outer number before
291 // any of the nested ones.
292 out = ad_type(n_independent_variables, index, in);
293
294 // Initialize potential nested directional derivatives
295 Marking<derivative_type>::independent_variable(
296 in, index, n_independent_variables, out.val());
297 }
298
299 /*
300 * Initialize the state of a dependent variable.
301 */
302 static void
303 dependent_variable(ad_type &out, const ad_type &func)
304 {
305 out = func;
306 }
307 };
308
309
314 template <typename NumberType>
315 struct Marking<Sacado::Rad::ADvar<NumberType>>
316 {
317 using ad_type =
318 typename SacadoNumberInfo<Sacado::Rad::ADvar<NumberType>>::ad_type;
319 using derivative_type = typename SacadoNumberInfo<
320 Sacado::Rad::ADvar<NumberType>>::derivative_type;
321 using scalar_type = typename SacadoNumberInfo<
322 Sacado::Rad::ADvar<NumberType>>::scalar_type;
323
324 /*
325 * Initialize the state of an independent variable.
326 */
327 static void
328 independent_variable(const scalar_type &in,
329 const unsigned int index,
330 const unsigned int n_independent_variables,
331 ad_type &out)
332 {
333 // For Sacado::Rad::ADvar numbers, we have to initialize the
334 // ADNumber with an already fully-configured value. This means
335 // that if this nests another ADNumber then the nested number
336 // must already be set up and ready for use.
337
338 // Initialize potential nested directional derivatives
339 derivative_type derivative_initializer;
340 Marking<derivative_type>::independent_variable(
341 in, index, n_independent_variables, derivative_initializer);
342
343 // Initialize the outer ad_type
344 out = derivative_initializer;
345 }
346
347 /*
348 * Initialize the state of a dependent variable.
349 */
350 static void
351 dependent_variable(ad_type &out, const ad_type &func)
352 {
353 out = func;
354 }
355 };
356
357
363 template <typename NumberType>
364 struct ExtractData<Sacado::Fad::DFad<NumberType>>
365 {
366 using derivative_type = typename SacadoNumberInfo<
367 Sacado::Fad::DFad<NumberType>>::derivative_type;
368 using scalar_type =
369 typename SacadoNumberInfo<Sacado::Fad::DFad<NumberType>>::scalar_type;
370 using value_type =
371 typename SacadoNumberInfo<Sacado::Fad::DFad<NumberType>>::value_type;
372
376 static scalar_type
377 value(const Sacado::Fad::DFad<NumberType> &x)
378 {
379 return ExtractData<value_type>::value(x.val());
380 }
381
382
386 static unsigned int
387 n_directional_derivatives(const Sacado::Fad::DFad<NumberType> &x)
388 {
389 return x.size();
390 }
391
392
396 static derivative_type
397 directional_derivative(const Sacado::Fad::DFad<NumberType> &x,
398 const unsigned int direction)
399 {
400 if (x.hasFastAccess())
401 return x.fastAccessDx(direction);
402 else
403 return x.dx(direction);
404 }
405 };
406
407
413 template <typename NumberType>
414 struct ExtractData<Sacado::Rad::ADvar<NumberType>>
415 {
416 using derivative_type = typename SacadoNumberInfo<
417 Sacado::Rad::ADvar<NumberType>>::derivative_type;
418 using scalar_type = typename SacadoNumberInfo<
419 Sacado::Rad::ADvar<NumberType>>::scalar_type;
420 using value_type =
421 typename SacadoNumberInfo<Sacado::Rad::ADvar<NumberType>>::value_type;
422
426 static scalar_type
427 value(const Sacado::Rad::ADvar<NumberType> &x)
428 {
429 return ExtractData<value_type>::value(x.val());
430 }
431
432
436 static unsigned int
437 n_directional_derivatives(const Sacado::Rad::ADvar<NumberType> &)
438 {
439 // There are as many directional derivatives as there are
440 // independent variables, but each independent variable can
441 // only return one directional derivative.
442 return 1;
443 }
444
445
453 static derivative_type
454 directional_derivative(const Sacado::Rad::ADvar<NumberType> &x,
455 const unsigned int)
456 {
457 return x.adj();
458 }
459 };
460
461 } // namespace internal
462
463
464 /* -------------- NumberTypes::sacado_dfad -------------- */
465
466
473 template <typename ADNumberType>
474 struct ADNumberTraits<
475 ADNumberType,
476 std::enable_if_t<
477 std::is_same_v<ADNumberType,
478 Sacado::Fad::DFad<typename ADNumberType::scalar_type>>>>
479 : NumberTraits<typename ADNumberType::scalar_type,
480 NumberTypes::sacado_dfad>
481 {};
482
483
490 template <typename ADNumberType>
491 struct ADNumberTraits<
492 ADNumberType,
493 std::enable_if_t<
494 std::is_same_v<ADNumberType,
495 std::complex<Sacado::Fad::DFad<
496 typename ADNumberType::value_type::scalar_type>>>>>
497 : NumberTraits<
498 std::complex<typename ADNumberType::value_type::scalar_type>,
499 NumberTypes::sacado_dfad>
500 {};
501
502
507 template <>
508 struct NumberTraits<Sacado::Fad::DFad<float>, NumberTypes::sacado_dfad>
509 : NumberTraits<
510 typename ADNumberTraits<Sacado::Fad::DFad<float>>::scalar_type,
511 NumberTypes::sacado_dfad>
512 {};
513
514
519 template <>
520 struct NumberTraits<std::complex<Sacado::Fad::DFad<float>>,
522 : NumberTraits<typename ADNumberTraits<
523 std::complex<Sacado::Fad::DFad<float>>>::scalar_type,
524 NumberTypes::sacado_dfad>
525 {};
526
527
532 template <>
533 struct NumberTraits<Sacado::Fad::DFad<double>, NumberTypes::sacado_dfad>
534 : NumberTraits<
535 typename ADNumberTraits<Sacado::Fad::DFad<double>>::scalar_type,
536 NumberTypes::sacado_dfad>
537 {};
538
539
544 template <>
545 struct NumberTraits<std::complex<Sacado::Fad::DFad<double>>,
547 : NumberTraits<typename ADNumberTraits<
548 std::complex<Sacado::Fad::DFad<double>>>::scalar_type,
549 NumberTypes::sacado_dfad>
550 {};
551
552
553 /* -------------- NumberTypes::sacado_rad -------------- */
554
555
562 template <typename ADNumberType>
563 struct ADNumberTraits<
564 ADNumberType,
565 std::enable_if_t<std::is_same_v<
566 ADNumberType,
567 Sacado::Rad::ADvar<typename ADNumberType::ADVari::scalar_type>>>>
568 : NumberTraits<typename ADNumberType::ADVari::scalar_type,
569 NumberTypes::sacado_rad>
570 {};
571
572
577 template <>
578 struct NumberTraits<Sacado::Rad::ADvar<float>, NumberTypes::sacado_rad>
579 : NumberTraits<
580 typename ADNumberTraits<Sacado::Rad::ADvar<float>>::scalar_type,
581 NumberTypes::sacado_rad>
582 {};
583
584
589 template <>
590 struct NumberTraits<Sacado::Rad::ADvar<double>, NumberTypes::sacado_rad>
591 : NumberTraits<
592 typename ADNumberTraits<Sacado::Rad::ADvar<double>>::scalar_type,
593 NumberTypes::sacado_rad>
594 {};
595
596
597# ifdef DEAL_II_TRILINOS_CXX_SUPPORTS_SACADO_COMPLEX_RAD
598
599
606 template <typename ADNumberType>
607 struct ADNumberTraits<
608 ADNumberType,
609 std::enable_if_t<std::is_same_v<
610 ADNumberType,
611 std::complex<Sacado::Rad::ADvar<
612 typename ADNumberType::value_type::ADVari::scalar_type>>>>>
613 : NumberTraits<
614 std::complex<typename ADNumberType::value_type::ADVari::scalar_type>,
615 NumberTypes::sacado_rad>
616 {};
617
618
623 template <>
624 struct NumberTraits<std::complex<Sacado::Rad::ADvar<float>>,
626 : NumberTraits<typename ADNumberTraits<
627 std::complex<Sacado::Rad::ADvar<float>>>::scalar_type,
628 NumberTypes::sacado_rad>
629 {};
630
631
636 template <>
637 struct NumberTraits<std::complex<Sacado::Rad::ADvar<double>>,
639 : NumberTraits<typename ADNumberTraits<
640 std::complex<Sacado::Rad::ADvar<double>>>::scalar_type,
641 NumberTypes::sacado_rad>
642 {};
643
644
645# endif
646
647
648 /* -------------- NumberTypes::sacado_dfad_dfad -------------- */
649
657 template <typename ADNumberType>
658 struct ADNumberTraits<
659 ADNumberType,
660 std::enable_if_t<std::is_same_v<ADNumberType,
661 Sacado::Fad::DFad<Sacado::Fad::DFad<
662 typename ADNumberType::scalar_type>>>>>
663 : NumberTraits<typename ADNumberType::scalar_type,
664 NumberTypes::sacado_dfad_dfad>
665 {};
666
667
675 template <typename ADNumberType>
676 struct ADNumberTraits<
677 ADNumberType,
678 std::enable_if_t<
679 std::is_same_v<ADNumberType,
680 std::complex<Sacado::Fad::DFad<Sacado::Fad::DFad<
681 typename ADNumberType::value_type::scalar_type>>>>>>
682 : NumberTraits<
683 std::complex<typename ADNumberType::value_type::scalar_type>,
684 NumberTypes::sacado_dfad_dfad>
685 {};
686
687
692 template <>
693 struct NumberTraits<Sacado::Fad::DFad<Sacado::Fad::DFad<float>>,
695 : NumberTraits<typename ADNumberTraits<Sacado::Fad::DFad<
696 Sacado::Fad::DFad<float>>>::scalar_type,
697 NumberTypes::sacado_dfad_dfad>
698 {};
699
700
705 template <>
706 struct NumberTraits<
707 std::complex<Sacado::Fad::DFad<Sacado::Fad::DFad<float>>>,
709 : NumberTraits<typename ADNumberTraits<std::complex<Sacado::Fad::DFad<
710 Sacado::Fad::DFad<float>>>>::scalar_type,
711 NumberTypes::sacado_dfad_dfad>
712 {};
713
714
719 template <>
720 struct NumberTraits<Sacado::Fad::DFad<Sacado::Fad::DFad<double>>,
722 : NumberTraits<typename ADNumberTraits<Sacado::Fad::DFad<
723 Sacado::Fad::DFad<double>>>::scalar_type,
724 NumberTypes::sacado_dfad_dfad>
725 {};
726
727
732 template <>
733 struct NumberTraits<
734 std::complex<Sacado::Fad::DFad<Sacado::Fad::DFad<double>>>,
736 : NumberTraits<typename ADNumberTraits<std::complex<Sacado::Fad::DFad<
737 Sacado::Fad::DFad<double>>>>::scalar_type,
738 NumberTypes::sacado_dfad_dfad>
739 {};
740
741
742 /* -------------- NumberTypes::sacado_rad_dfad -------------- */
743
751 template <typename ADNumberType>
752 struct ADNumberTraits<ADNumberType,
753 std::enable_if_t<std::is_same_v<
754 ADNumberType,
755 Sacado::Rad::ADvar<Sacado::Fad::DFad<
756 typename ADNumberType::ADVari::scalar_type>>>>>
757 : NumberTraits<typename ADNumberType::ADVari::scalar_type,
758 NumberTypes::sacado_rad_dfad>
759 {};
760
761
766 template <>
767 struct NumberTraits<Sacado::Rad::ADvar<Sacado::Fad::DFad<float>>,
769 : NumberTraits<typename ADNumberTraits<Sacado::Rad::ADvar<
770 Sacado::Fad::DFad<float>>>::scalar_type,
771 NumberTypes::sacado_rad_dfad>
772 {};
773
774
779 template <>
780 struct NumberTraits<Sacado::Rad::ADvar<Sacado::Fad::DFad<double>>,
782 : NumberTraits<typename ADNumberTraits<Sacado::Rad::ADvar<
783 Sacado::Fad::DFad<double>>>::scalar_type,
784 NumberTypes::sacado_rad_dfad>
785 {};
786
787
788# ifdef DEAL_II_TRILINOS_CXX_SUPPORTS_SACADO_COMPLEX_RAD
789
790
798 template <typename ADNumberType>
799 struct ADNumberTraits<
800 ADNumberType,
801 std::enable_if_t<std::is_same_v<
802 ADNumberType,
803 std::complex<Sacado::Rad::ADvar<Sacado::Fad::DFad<
804 typename ADNumberType::value_type::ADVari::scalar_type>>>>>>
805 : NumberTraits<
806 std::complex<typename ADNumberType::value_type::ADVari::scalar_type>,
807 NumberTypes::sacado_rad_dfad>
808 {};
809
810
815 template <>
816 struct NumberTraits<
817 std::complex<Sacado::Rad::ADvar<Sacado::Fad::DFad<float>>>,
819 : NumberTraits<typename ADNumberTraits<std::complex<Sacado::Rad::ADvar<
820 Sacado::Fad::DFad<float>>>>::scalar_type,
821 NumberTypes::sacado_rad_dfad>
822 {};
823
824
829 template <>
830 struct NumberTraits<
831 std::complex<Sacado::Rad::ADvar<Sacado::Fad::DFad<double>>>,
833 : NumberTraits<typename ADNumberTraits<std::complex<Sacado::Rad::ADvar<
834 Sacado::Fad::DFad<double>>>>::scalar_type,
835 NumberTypes::sacado_rad_dfad>
836 {};
837
838
839# endif
840
841
842 /* -------------- Additional type traits -------------- */
843
844
845 template <typename NumberType>
846 struct is_sacado_dfad_number<
847 NumberType,
848 std::enable_if_t<ADNumberTraits<std::decay_t<NumberType>>::type_code ==
849 NumberTypes::sacado_dfad ||
850 ADNumberTraits<std::decay_t<NumberType>>::type_code ==
851 NumberTypes::sacado_dfad_dfad>> : std::true_type
852 {};
853
854
855 template <typename NumberType>
856 struct is_sacado_dfad_number<
857 NumberType,
858 std::enable_if_t<
859 std::is_same_v<NumberType,
860 Sacado::Fad::Expr<typename NumberType::value_type>>>>
861 : std::true_type
862 {};
863
864
865 template <typename NumberType>
866 struct is_sacado_rad_number<
867 NumberType,
868 std::enable_if_t<ADNumberTraits<std::decay_t<NumberType>>::type_code ==
869 NumberTypes::sacado_rad ||
870 ADNumberTraits<std::decay_t<NumberType>>::type_code ==
871 NumberTypes::sacado_rad_dfad>> : std::true_type
872 {};
873
874
875 template <typename NumberType>
876 struct is_sacado_rad_number<
877 NumberType,
878 std::enable_if_t<
879 std::is_same_v<NumberType,
880 Sacado::Rad::ADvari<Sacado::Fad::DFad<
881 typename NumberType::ADVari::scalar_type>>>>>
882 : std::true_type
883 {};
884
885
886 template <typename NumberType>
887 struct is_sacado_number<
888 NumberType,
889 std::enable_if_t<is_sacado_dfad_number<NumberType>::value ||
890 is_sacado_rad_number<NumberType>::value>>
891 : std::true_type
892 {};
893
894 } // namespace AD
895} // namespace Differentiation
896
897
898# endif // DOXYGEN
899
900
901#endif // DEAL_II_TRILINOS_WITH_SACADO
902
903
905
906#endif
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
STL namespace.