deal.II version GIT relicensing-2167-g9622207b8f 2024-11-21 12:40:00+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
20#include <type_traits>
21
22
24
25
26namespace Differentiation
27{
28 namespace AD
29 {
35 template <typename NumberType, typename = void>
36 struct is_sacado_number : std::false_type
37 {};
38
39
45 template <typename NumberType, typename = void>
46 struct is_sacado_dfad_number : std::false_type
47 {};
48
49
55 template <typename NumberType, typename = void>
56 struct is_sacado_rad_number : std::false_type
57 {};
58
59 } // namespace AD
60} // namespace Differentiation
61
62
64
65
66
67#ifdef DEAL_II_TRILINOS_WITH_SACADO
68
70# include <deal.II/base/numbers.h>
71
74
75// It appears that some versions of Trilinos do not directly or indirectly
76// include all the headers for all forward and reverse Sacado AD types
77// in Sacado.hpp, so we also directly include these here as a precaution:
78# include <Sacado.hpp>
79# include <Sacado_Fad_DFad.hpp>
80# include <Sacado_trad.hpp>
81
82# include <complex>
83
85
86
87namespace Differentiation
88{
89 namespace AD
90 {
91 namespace internal
92 {
99 template <typename SacadoNumber, typename = void>
101
102 } // namespace internal
103
104
105
106 } // namespace AD
107} // namespace Differentiation
108
109
110/* ----------- inline and template functions and specializations ----------- */
111
112
113# ifndef DOXYGEN
114
115namespace Differentiation
116{
117 namespace AD
118 {
119 namespace internal
120 {
121 // The documentation on Sacado numbers is pretty sparse and/or hard to
122 // navigate. As a point of reference, see
123 // https://trilinos.org/docs/dev/packages/sacado/doc/html/classSacado_1_1Fad_1_1SimpleFad.html
124 // for semi-applicable documentation for the Sacado::Fad::Dfad class.
125 // and the examples in
126 // https://github.com/trilinos/Trilinos/tree/master/packages/sacado/example
127 //
128 // If one dares to venture there, the relevant files for the classes
129 // supported here are:
130 //
131 // Forward-mode auto-differentiable types:
132 // https://github.com/trilinos/Trilinos/blob/master/packages/sacado/src/sacado_dfad_DFad.hpp
133 // https://github.com/trilinos/Trilinos/blob/master/packages/sacado/src/sacado_dfad_GeneralFad.hpp
134 //
135 // Reverse-mode auto-differentiable types:
136 // https://github.com/trilinos/Trilinos/blob/master/packages/sacado/src/Sacado_trad.hpp
137
138
142 template <typename SacadoNumber>
143 struct SacadoNumberInfo<
144 SacadoNumber,
145 std::enable_if_t<
146 std::is_same_v<SacadoNumber,
147 Sacado::Fad::DFad<typename SacadoNumber::value_type>>>>
148 {
149 using ad_type = SacadoNumber;
150 using scalar_type = typename ad_type::scalar_type;
151 using value_type = typename ad_type::value_type;
152 using derivative_type = typename ad_type::value_type;
153
154 static const unsigned int n_supported_derivative_levels =
155 1 + SacadoNumberInfo<derivative_type>::n_supported_derivative_levels;
156 };
157
158
162 template <typename SacadoNumber>
163 struct SacadoNumberInfo<
164 SacadoNumber,
165 std::enable_if_t<std::is_same_v<
166 SacadoNumber,
167 Sacado::Rad::ADvar<typename SacadoNumber::value_type>>>>
168 {
169 using ad_type = SacadoNumber;
170 using scalar_type = typename ad_type::ADVari::scalar_type;
171 using value_type = typename ad_type::ADVari::value_type;
172 using derivative_type = typename ad_type::ADVari::value_type;
173
174 static const unsigned int n_supported_derivative_levels =
175 1 + SacadoNumberInfo<derivative_type>::n_supported_derivative_levels;
176 };
177
178
185 template <typename Number>
186 struct SacadoNumberInfo<
187 Number,
188 std::enable_if_t<std::is_arithmetic_v<std::decay_t<Number>>>>
189 {
190 static const unsigned int n_supported_derivative_levels = 0;
191 };
192
193
198 template <typename ScalarType>
199 struct ADNumberInfoFromEnum<
200 ScalarType,
202 std::enable_if_t<std::is_floating_point_v<ScalarType>>>
203 {
204 static const bool is_taped = false;
205 using real_type = Sacado::Fad::DFad<ScalarType>;
206 using derivative_type =
207 typename SacadoNumberInfo<real_type>::derivative_type;
208 static const unsigned int n_supported_derivative_levels =
209 SacadoNumberInfo<real_type>::n_supported_derivative_levels;
210 };
211
212
217 template <typename ScalarType>
218 struct ADNumberInfoFromEnum<
219 ScalarType,
221 std::enable_if_t<std::is_floating_point_v<ScalarType>>>
222 {
223 static const bool is_taped = false;
224 using real_type = Sacado::Fad::DFad<Sacado::Fad::DFad<ScalarType>>;
225 using derivative_type =
226 typename SacadoNumberInfo<real_type>::derivative_type;
227 static const unsigned int n_supported_derivative_levels =
228 SacadoNumberInfo<real_type>::n_supported_derivative_levels;
229 };
230
231
236 template <typename ScalarType>
237 struct ADNumberInfoFromEnum<
238 ScalarType,
240 std::enable_if_t<std::is_floating_point_v<ScalarType>>>
241 {
242 static const bool is_taped = false;
243 using real_type = Sacado::Rad::ADvar<ScalarType>;
244 using derivative_type =
245 typename SacadoNumberInfo<real_type>::derivative_type;
246 static const unsigned int n_supported_derivative_levels =
247 SacadoNumberInfo<real_type>::n_supported_derivative_levels;
248 };
249
250
255 template <typename ScalarType>
256 struct ADNumberInfoFromEnum<
257 ScalarType,
259 std::enable_if_t<std::is_floating_point_v<ScalarType>>>
260 {
261 static const bool is_taped = false;
262 using real_type = Sacado::Rad::ADvar<Sacado::Fad::DFad<ScalarType>>;
263 using derivative_type =
264 typename SacadoNumberInfo<real_type>::derivative_type;
265 static const unsigned int n_supported_derivative_levels =
266 SacadoNumberInfo<real_type>::n_supported_derivative_levels;
267 };
268
269
274 template <typename NumberType>
275 struct Marking<Sacado::Fad::DFad<NumberType>>
276 {
277 using ad_type =
278 typename SacadoNumberInfo<Sacado::Fad::DFad<NumberType>>::ad_type;
279 using derivative_type = typename SacadoNumberInfo<
280 Sacado::Fad::DFad<NumberType>>::derivative_type;
281 using scalar_type =
282 typename SacadoNumberInfo<Sacado::Fad::DFad<NumberType>>::scalar_type;
283
284 /*
285 * Initialize the state of an independent variable.
286 */
287 static void
288 independent_variable(const scalar_type &in,
289 const unsigned int index,
290 const unsigned int n_independent_variables,
291 ad_type &out)
292 {
293 // It is required that we first initialize the outer number before
294 // any of the nested ones.
295 out = ad_type(n_independent_variables, index, in);
296
297 // Initialize potential nested directional derivatives
298 Marking<derivative_type>::independent_variable(
299 in, index, n_independent_variables, out.val());
300 }
301
302 /*
303 * Initialize the state of a dependent variable.
304 */
305 static void
306 dependent_variable(ad_type &out, const ad_type &func)
307 {
308 out = func;
309 }
310 };
311
312
317 template <typename NumberType>
318 struct Marking<Sacado::Rad::ADvar<NumberType>>
319 {
320 using ad_type =
321 typename SacadoNumberInfo<Sacado::Rad::ADvar<NumberType>>::ad_type;
322 using derivative_type = typename SacadoNumberInfo<
323 Sacado::Rad::ADvar<NumberType>>::derivative_type;
324 using scalar_type = typename SacadoNumberInfo<
325 Sacado::Rad::ADvar<NumberType>>::scalar_type;
326
327 /*
328 * Initialize the state of an independent variable.
329 */
330 static void
331 independent_variable(const scalar_type &in,
332 const unsigned int index,
333 const unsigned int n_independent_variables,
334 ad_type &out)
335 {
336 // For Sacado::Rad::ADvar numbers, we have to initialize the
337 // ADNumber with an already fully-configured value. This means
338 // that if this nests another ADNumber then the nested number
339 // must already be set up and ready for use.
340
341 // Initialize potential nested directional derivatives
342 derivative_type derivative_initializer;
343 Marking<derivative_type>::independent_variable(
344 in, index, n_independent_variables, derivative_initializer);
345
346 // Initialize the outer ad_type
347 out = derivative_initializer;
348 }
349
350 /*
351 * Initialize the state of a dependent variable.
352 */
353 static void
354 dependent_variable(ad_type &out, const ad_type &func)
355 {
356 out = func;
357 }
358 };
359
360
366 template <typename NumberType>
367 struct ExtractData<Sacado::Fad::DFad<NumberType>>
368 {
369 using derivative_type = typename SacadoNumberInfo<
370 Sacado::Fad::DFad<NumberType>>::derivative_type;
371 using scalar_type =
372 typename SacadoNumberInfo<Sacado::Fad::DFad<NumberType>>::scalar_type;
373 using value_type =
374 typename SacadoNumberInfo<Sacado::Fad::DFad<NumberType>>::value_type;
375
379 static scalar_type
380 value(const Sacado::Fad::DFad<NumberType> &x)
381 {
382 return ExtractData<value_type>::value(x.val());
383 }
384
385
389 static unsigned int
390 n_directional_derivatives(const Sacado::Fad::DFad<NumberType> &x)
391 {
392 return x.size();
393 }
394
395
399 static derivative_type
400 directional_derivative(const Sacado::Fad::DFad<NumberType> &x,
401 const unsigned int direction)
402 {
403 if (x.hasFastAccess())
404 return x.fastAccessDx(direction);
405 else
406 return x.dx(direction);
407 }
408 };
409
410
416 template <typename NumberType>
417 struct ExtractData<Sacado::Rad::ADvar<NumberType>>
418 {
419 using derivative_type = typename SacadoNumberInfo<
420 Sacado::Rad::ADvar<NumberType>>::derivative_type;
421 using scalar_type = typename SacadoNumberInfo<
422 Sacado::Rad::ADvar<NumberType>>::scalar_type;
423 using value_type =
424 typename SacadoNumberInfo<Sacado::Rad::ADvar<NumberType>>::value_type;
425
429 static scalar_type
430 value(const Sacado::Rad::ADvar<NumberType> &x)
431 {
432 return ExtractData<value_type>::value(x.val());
433 }
434
435
439 static unsigned int
440 n_directional_derivatives(const Sacado::Rad::ADvar<NumberType> &)
441 {
442 // There are as many directional derivatives as there are
443 // independent variables, but each independent variable can
444 // only return one directional derivative.
445 return 1;
446 }
447
448
456 static derivative_type
457 directional_derivative(const Sacado::Rad::ADvar<NumberType> &x,
458 const unsigned int)
459 {
460 return x.adj();
461 }
462 };
463
464 } // namespace internal
465
466
467 /* -------------- NumberTypes::sacado_dfad -------------- */
468
469
476 template <typename ADNumberType>
477 struct ADNumberTraits<
478 ADNumberType,
479 std::enable_if_t<
480 std::is_same_v<ADNumberType,
481 Sacado::Fad::DFad<typename ADNumberType::scalar_type>>>>
482 : NumberTraits<typename ADNumberType::scalar_type,
483 NumberTypes::sacado_dfad>
484 {};
485
486
493 template <typename ADNumberType>
494 struct ADNumberTraits<
495 ADNumberType,
496 std::enable_if_t<
497 std::is_same_v<ADNumberType,
498 std::complex<Sacado::Fad::DFad<
499 typename ADNumberType::value_type::scalar_type>>>>>
500 : NumberTraits<
501 std::complex<typename ADNumberType::value_type::scalar_type>,
502 NumberTypes::sacado_dfad>
503 {};
504
505
510 template <>
511 struct NumberTraits<Sacado::Fad::DFad<float>, NumberTypes::sacado_dfad>
512 : NumberTraits<
513 typename ADNumberTraits<Sacado::Fad::DFad<float>>::scalar_type,
514 NumberTypes::sacado_dfad>
515 {};
516
517
522 template <>
523 struct NumberTraits<std::complex<Sacado::Fad::DFad<float>>,
525 : NumberTraits<typename ADNumberTraits<
526 std::complex<Sacado::Fad::DFad<float>>>::scalar_type,
527 NumberTypes::sacado_dfad>
528 {};
529
530
535 template <>
536 struct NumberTraits<Sacado::Fad::DFad<double>, NumberTypes::sacado_dfad>
537 : NumberTraits<
538 typename ADNumberTraits<Sacado::Fad::DFad<double>>::scalar_type,
539 NumberTypes::sacado_dfad>
540 {};
541
542
547 template <>
548 struct NumberTraits<std::complex<Sacado::Fad::DFad<double>>,
550 : NumberTraits<typename ADNumberTraits<
551 std::complex<Sacado::Fad::DFad<double>>>::scalar_type,
552 NumberTypes::sacado_dfad>
553 {};
554
555
556 /* -------------- NumberTypes::sacado_rad -------------- */
557
558
565 template <typename ADNumberType>
566 struct ADNumberTraits<
567 ADNumberType,
568 std::enable_if_t<std::is_same_v<
569 ADNumberType,
570 Sacado::Rad::ADvar<typename ADNumberType::ADVari::scalar_type>>>>
571 : NumberTraits<typename ADNumberType::ADVari::scalar_type,
572 NumberTypes::sacado_rad>
573 {};
574
575
580 template <>
581 struct NumberTraits<Sacado::Rad::ADvar<float>, NumberTypes::sacado_rad>
582 : NumberTraits<
583 typename ADNumberTraits<Sacado::Rad::ADvar<float>>::scalar_type,
584 NumberTypes::sacado_rad>
585 {};
586
587
592 template <>
593 struct NumberTraits<Sacado::Rad::ADvar<double>, NumberTypes::sacado_rad>
594 : NumberTraits<
595 typename ADNumberTraits<Sacado::Rad::ADvar<double>>::scalar_type,
596 NumberTypes::sacado_rad>
597 {};
598
599
600# ifdef DEAL_II_TRILINOS_CXX_SUPPORTS_SACADO_COMPLEX_RAD
601
602
609 template <typename ADNumberType>
610 struct ADNumberTraits<
611 ADNumberType,
612 std::enable_if_t<std::is_same_v<
613 ADNumberType,
614 std::complex<Sacado::Rad::ADvar<
615 typename ADNumberType::value_type::ADVari::scalar_type>>>>>
616 : NumberTraits<
617 std::complex<typename ADNumberType::value_type::ADVari::scalar_type>,
618 NumberTypes::sacado_rad>
619 {};
620
621
626 template <>
627 struct NumberTraits<std::complex<Sacado::Rad::ADvar<float>>,
629 : NumberTraits<typename ADNumberTraits<
630 std::complex<Sacado::Rad::ADvar<float>>>::scalar_type,
631 NumberTypes::sacado_rad>
632 {};
633
634
639 template <>
640 struct NumberTraits<std::complex<Sacado::Rad::ADvar<double>>,
642 : NumberTraits<typename ADNumberTraits<
643 std::complex<Sacado::Rad::ADvar<double>>>::scalar_type,
644 NumberTypes::sacado_rad>
645 {};
646
647
648# endif
649
650
651 /* -------------- NumberTypes::sacado_dfad_dfad -------------- */
652
660 template <typename ADNumberType>
661 struct ADNumberTraits<
662 ADNumberType,
663 std::enable_if_t<std::is_same_v<ADNumberType,
664 Sacado::Fad::DFad<Sacado::Fad::DFad<
665 typename ADNumberType::scalar_type>>>>>
666 : NumberTraits<typename ADNumberType::scalar_type,
667 NumberTypes::sacado_dfad_dfad>
668 {};
669
670
678 template <typename ADNumberType>
679 struct ADNumberTraits<
680 ADNumberType,
681 std::enable_if_t<
682 std::is_same_v<ADNumberType,
683 std::complex<Sacado::Fad::DFad<Sacado::Fad::DFad<
684 typename ADNumberType::value_type::scalar_type>>>>>>
685 : NumberTraits<
686 std::complex<typename ADNumberType::value_type::scalar_type>,
687 NumberTypes::sacado_dfad_dfad>
688 {};
689
690
695 template <>
696 struct NumberTraits<Sacado::Fad::DFad<Sacado::Fad::DFad<float>>,
698 : NumberTraits<typename ADNumberTraits<Sacado::Fad::DFad<
699 Sacado::Fad::DFad<float>>>::scalar_type,
700 NumberTypes::sacado_dfad_dfad>
701 {};
702
703
708 template <>
709 struct NumberTraits<
710 std::complex<Sacado::Fad::DFad<Sacado::Fad::DFad<float>>>,
712 : NumberTraits<typename ADNumberTraits<std::complex<Sacado::Fad::DFad<
713 Sacado::Fad::DFad<float>>>>::scalar_type,
714 NumberTypes::sacado_dfad_dfad>
715 {};
716
717
722 template <>
723 struct NumberTraits<Sacado::Fad::DFad<Sacado::Fad::DFad<double>>,
725 : NumberTraits<typename ADNumberTraits<Sacado::Fad::DFad<
726 Sacado::Fad::DFad<double>>>::scalar_type,
727 NumberTypes::sacado_dfad_dfad>
728 {};
729
730
735 template <>
736 struct NumberTraits<
737 std::complex<Sacado::Fad::DFad<Sacado::Fad::DFad<double>>>,
739 : NumberTraits<typename ADNumberTraits<std::complex<Sacado::Fad::DFad<
740 Sacado::Fad::DFad<double>>>>::scalar_type,
741 NumberTypes::sacado_dfad_dfad>
742 {};
743
744
745 /* -------------- NumberTypes::sacado_rad_dfad -------------- */
746
754 template <typename ADNumberType>
755 struct ADNumberTraits<ADNumberType,
756 std::enable_if_t<std::is_same_v<
757 ADNumberType,
758 Sacado::Rad::ADvar<Sacado::Fad::DFad<
759 typename ADNumberType::ADVari::scalar_type>>>>>
760 : NumberTraits<typename ADNumberType::ADVari::scalar_type,
761 NumberTypes::sacado_rad_dfad>
762 {};
763
764
769 template <>
770 struct NumberTraits<Sacado::Rad::ADvar<Sacado::Fad::DFad<float>>,
772 : NumberTraits<typename ADNumberTraits<Sacado::Rad::ADvar<
773 Sacado::Fad::DFad<float>>>::scalar_type,
774 NumberTypes::sacado_rad_dfad>
775 {};
776
777
782 template <>
783 struct NumberTraits<Sacado::Rad::ADvar<Sacado::Fad::DFad<double>>,
785 : NumberTraits<typename ADNumberTraits<Sacado::Rad::ADvar<
786 Sacado::Fad::DFad<double>>>::scalar_type,
787 NumberTypes::sacado_rad_dfad>
788 {};
789
790
791# ifdef DEAL_II_TRILINOS_CXX_SUPPORTS_SACADO_COMPLEX_RAD
792
793
801 template <typename ADNumberType>
802 struct ADNumberTraits<
803 ADNumberType,
804 std::enable_if_t<std::is_same_v<
805 ADNumberType,
806 std::complex<Sacado::Rad::ADvar<Sacado::Fad::DFad<
807 typename ADNumberType::value_type::ADVari::scalar_type>>>>>>
808 : NumberTraits<
809 std::complex<typename ADNumberType::value_type::ADVari::scalar_type>,
810 NumberTypes::sacado_rad_dfad>
811 {};
812
813
818 template <>
819 struct NumberTraits<
820 std::complex<Sacado::Rad::ADvar<Sacado::Fad::DFad<float>>>,
822 : NumberTraits<typename ADNumberTraits<std::complex<Sacado::Rad::ADvar<
823 Sacado::Fad::DFad<float>>>>::scalar_type,
824 NumberTypes::sacado_rad_dfad>
825 {};
826
827
832 template <>
833 struct NumberTraits<
834 std::complex<Sacado::Rad::ADvar<Sacado::Fad::DFad<double>>>,
836 : NumberTraits<typename ADNumberTraits<std::complex<Sacado::Rad::ADvar<
837 Sacado::Fad::DFad<double>>>>::scalar_type,
838 NumberTypes::sacado_rad_dfad>
839 {};
840
841
842# endif
843
844
845 /* -------------- Additional type traits -------------- */
846
847
848 template <typename NumberType>
849 struct is_sacado_dfad_number<
850 NumberType,
851 std::enable_if_t<ADNumberTraits<std::decay_t<NumberType>>::type_code ==
852 NumberTypes::sacado_dfad ||
853 ADNumberTraits<std::decay_t<NumberType>>::type_code ==
854 NumberTypes::sacado_dfad_dfad>> : std::true_type
855 {};
856
857
858 template <typename NumberType>
859 struct is_sacado_dfad_number<
860 NumberType,
861 std::enable_if_t<
862 std::is_same_v<NumberType,
863 Sacado::Fad::Expr<typename NumberType::value_type>>>>
864 : std::true_type
865 {};
866
867
868 template <typename NumberType>
869 struct is_sacado_rad_number<
870 NumberType,
871 std::enable_if_t<ADNumberTraits<std::decay_t<NumberType>>::type_code ==
872 NumberTypes::sacado_rad ||
873 ADNumberTraits<std::decay_t<NumberType>>::type_code ==
874 NumberTypes::sacado_rad_dfad>> : std::true_type
875 {};
876
877
878 template <typename NumberType>
879 struct is_sacado_rad_number<
880 NumberType,
881 std::enable_if_t<
882 std::is_same_v<NumberType,
883 Sacado::Rad::ADvari<Sacado::Fad::DFad<
884 typename NumberType::ADVari::scalar_type>>>>>
885 : std::true_type
886 {};
887
888
889 template <typename NumberType>
890 struct is_sacado_number<
891 NumberType,
892 std::enable_if_t<is_sacado_dfad_number<NumberType>::value ||
893 is_sacado_rad_number<NumberType>::value>>
894 : std::true_type
895 {};
896
897 } // namespace AD
898} // namespace Differentiation
899
900
901# endif // DOXYGEN
902
903
904
906
907
908#endif // DEAL_II_TRILINOS_WITH_SACADO
909
910#endif
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
STL namespace.