Reference documentation for deal.II version GIT e22cb6a53e 2023-09-21 14:00: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\}}\)
numbers.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2006 - 2023 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_numbers_h
17 #define dealii_numbers_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/types.h>
23 
24 #ifdef DEAL_II_WITH_CUDA
25 # include <cuComplex.h>
26 #endif
27 
28 #include <Kokkos_Macros.hpp>
29 
30 #include <cmath>
31 #include <complex>
32 #include <cstddef>
33 #include <type_traits>
34 
35 #define DEAL_II_HOST_DEVICE KOKKOS_FUNCTION
36 #define DEAL_II_CUDA_HOST_DEV DEAL_II_HOST_DEVICE
37 #define DEAL_II_HOST_DEVICE_ALWAYS_INLINE KOKKOS_FORCEINLINE_FUNCTION
38 
39 // clang++ assumes that all constexpr functions are __host__ __device__ when
40 // Kokkos was configured with CUDA or HIP support. This is problematic
41 // when calling non-constexpr functions in constexpr functions. Hence, we
42 // need a way to annotate functions explicitly as host-only.
43 #if (defined(__clang__) && defined(__CUDA__)) || defined(KOKKOS_ENABLE_HIP)
44 # define DEAL_II_HOST __host__
45 #else
46 # define DEAL_II_HOST
47 #endif
48 
49 // Forward-declare the automatic differentiation types so we can add prototypes
50 // for our own wrappers.
51 #ifdef DEAL_II_WITH_ADOLC
52 class adouble;
53 namespace adtl
54 {
55  class adouble;
56 }
57 #endif
58 
60 
61 namespace internal
62 {
79  template <typename Number>
81  {
85  constexpr static unsigned int max_width = 1;
86  };
87 
94  template <>
96  {
100  constexpr static unsigned int max_width =
101 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512
102  8;
103 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256
104  4;
105 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128
106  2;
107 #else
108  1;
109 #endif
110  };
111 
118  template <>
120  {
124  constexpr static unsigned int max_width =
125 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ALTIVEC__)
126  4;
127 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512 && defined(__AVX512F__)
128  16;
129 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256 && defined(__AVX__)
130  8;
131 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__SSE2__)
132  4;
133 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ARM_NEON)
134  4;
135 #else
136  1;
137 #endif
138  };
139 
140 
141 } // namespace internal
142 
143 // forward declarations to support abs or sqrt operations on VectorizedArray
144 #ifndef DOXYGEN
145 template <typename Number,
146  std::size_t width =
148 class VectorizedArray;
149 template <typename T>
150 struct EnableIfScalar;
151 #endif
152 
153 #ifdef DEAL_II_WITH_ADOLC
154 # ifndef DOXYGEN
155 // Prototype some inline functions present in adolc_math.h for use in
156 // NumberTraits.
157 //
158 // ADOL-C uses fabs(), but for genericity we want to use abs(). Simultaneously,
159 // though, we don't want to include ADOL-C headers in this header since
160 // numbers.h is in everything. To get around this: use C++ rules which permit
161 // the use of forward-declared classes in function prototypes to declare some
162 // functions which are defined in adolc_math.h. This permits us to write "using
163 // ::abs;" in NumberTraits which will allow us to select the correct
164 // overload (the one in ::) when instantiating NumberTraits for ADOL-C
165 // types.
166 
167 adouble
168 abs(const adouble &x);
169 
170 adtl::adouble
171 abs(const adtl::adouble &x);
172 # endif
173 #endif
174 
176 
177 namespace std
178 {
179  template <typename Number, std::size_t width>
180  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
181  sqrt(const ::VectorizedArray<Number, width> &);
182  template <typename Number, std::size_t width>
183  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
184  abs(const ::VectorizedArray<Number, width> &);
185  template <typename Number, std::size_t width>
186  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
187  max(const ::VectorizedArray<Number, width> &,
188  const ::VectorizedArray<Number, width> &);
189  template <typename Number, std::size_t width>
190  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
191  min(const ::VectorizedArray<Number, width> &,
192  const ::VectorizedArray<Number, width> &);
193  template <typename Number, size_t width>
195  pow(const ::VectorizedArray<Number, width> &, const Number p);
196  template <typename Number, size_t width>
198  sin(const ::VectorizedArray<Number, width> &);
199  template <typename Number, size_t width>
201  cos(const ::VectorizedArray<Number, width> &);
202  template <typename Number, size_t width>
204  tan(const ::VectorizedArray<Number, width> &);
205  template <typename Number, size_t width>
207  exp(const ::VectorizedArray<Number, width> &);
208  template <typename Number, size_t width>
210  log(const ::VectorizedArray<Number, width> &);
211 } // namespace std
212 
214 
230 namespace numbers
231 {
235  static constexpr double E = 2.7182818284590452354;
236 
240  static constexpr double LOG2E = 1.4426950408889634074;
241 
245  static constexpr double LOG10E = 0.43429448190325182765;
246 
250  static constexpr double LN2 = 0.69314718055994530942;
251 
255  static constexpr double LN10 = 2.30258509299404568402;
256 
260  static constexpr double PI = 3.14159265358979323846;
261 
265  static constexpr double PI_2 = 1.57079632679489661923;
266 
270  static constexpr double PI_4 = 0.78539816339744830962;
271 
275  static constexpr double SQRT2 = 1.41421356237309504880;
276 
280  static constexpr double SQRT1_2 = 0.70710678118654752440;
281 
291  bool
292  is_finite(const double x);
293 
298  bool
299  is_finite(const std::complex<double> &x);
300 
305  bool
306  is_finite(const std::complex<float> &x);
307 
316  bool
317  is_finite(const std::complex<long double> &x);
318 
329  template <typename Number1, typename Number2>
330  constexpr DEAL_II_HOST_DEVICE bool
331  values_are_equal(const Number1 &value_1, const Number2 &value_2);
332 
343  template <typename Number1, typename Number2>
344  bool
345  values_are_not_equal(const Number1 &value_1, const Number2 &value_2);
346 
354  template <typename Number>
355  constexpr DEAL_II_HOST_DEVICE bool
356  value_is_zero(const Number &value);
357 
368  template <typename Number1, typename Number2>
369  bool
370  value_is_less_than(const Number1 &value_1, const Number2 &value_2);
371 
382  template <typename Number1, typename Number2>
383  bool
384  value_is_less_than_or_equal_to(const Number1 &value_1,
385  const Number2 &value_2);
386 
387 
388 
399  template <typename Number1, typename Number2>
400  bool
401  value_is_greater_than(const Number1 &value_1, const Number2 &value_2);
402 
413  template <typename Number1, typename Number2>
414  bool
415  value_is_greater_than_or_equal_to(const Number1 &value_1,
416  const Number2 &value_2);
417 
426  template <typename number>
428  {
434  static constexpr bool is_complex = false;
435 
442  using real_type = number;
443 
448 
456  static constexpr DEAL_II_HOST_DEVICE const number &
457  conjugate(const number &x);
458 
467  static constexpr DEAL_II_HOST_DEVICE real_type
468  abs_square(const number &x);
469 
473  static real_type
474  abs(const number &x);
475  };
476 
477 
482  template <typename number>
483  struct NumberTraits<std::complex<number>>
484  {
490  static constexpr bool is_complex = true;
491 
498  using real_type = number;
499 
503  using double_type = std::complex<double>;
504 
508  static constexpr std::complex<number>
509  conjugate(const std::complex<number> &x);
510 
517  static constexpr real_type
518  abs_square(const std::complex<number> &x);
519 
520 
524  static real_type
525  abs(const std::complex<number> &x);
526  };
527 
528  // --------------- inline and template functions ---------------- //
529 
530  inline bool
531  is_nan(const double x)
532  {
533  return std::isnan(x);
534  }
535 
536 
537 
538  inline bool
539  is_finite(const double x)
540  {
541  return std::isfinite(x);
542  }
543 
544 
545 
546  inline bool
547  is_finite(const std::complex<double> &x)
548  {
549  // Check complex numbers for infinity
550  // by testing real and imaginary part
551  return (is_finite(x.real()) && is_finite(x.imag()));
552  }
553 
554 
555 
556  inline bool
557  is_finite(const std::complex<float> &x)
558  {
559  // Check complex numbers for infinity
560  // by testing real and imaginary part
561  return (is_finite(x.real()) && is_finite(x.imag()));
562  }
563 
564 
565 
566  inline bool
567  is_finite(const std::complex<long double> &x)
568  {
569  // Same for std::complex<long double>
570  return (is_finite(x.real()) && is_finite(x.imag()));
571  }
572 
573 
574  template <typename number>
575  constexpr DEAL_II_HOST_DEVICE const number &
577  {
578  return x;
579  }
580 
581 
582 
583  template <typename number>
586  {
587  return x * x;
588  }
589 
590 
591 
592  template <typename number>
594  NumberTraits<number>::abs(const number &x)
595  {
596  // Make things work with AD types
597  using std::abs;
598 #ifdef DEAL_II_WITH_ADOLC
599  // This one is a little tricky - we have our own abs function in ::,
600  // prototyped with forward-declared types in this file, but it only exists
601  // if we have ADOL-C: hence we only add this using statement in that
602  // situation
603  using ::abs;
604 #endif
605  return abs(x);
606  }
607 
608 
609 
610  template <typename number>
611  constexpr std::complex<number>
612  NumberTraits<std::complex<number>>::conjugate(const std::complex<number> &x)
613  {
614  return std::conj(x);
615  }
616 
617 
618 
619  template <typename number>
620  typename NumberTraits<std::complex<number>>::real_type
621  NumberTraits<std::complex<number>>::abs(const std::complex<number> &x)
622  {
623  // Make things work with AD types
624  using std::abs;
625 #ifdef DEAL_II_WITH_ADOLC
626  // Same comment as the non-complex case holds here
627  using ::abs;
628 #endif
629  return abs(x);
630  }
631 
632 
633 
634  template <typename number>
635  constexpr typename NumberTraits<std::complex<number>>::real_type
636  NumberTraits<std::complex<number>>::abs_square(const std::complex<number> &x)
637  {
638  return std::norm(x);
639  }
640 
641 } // namespace numbers
642 
643 
644 // Forward declarations
646 {
647  namespace AD
648  {
649  namespace internal
650  {
651  // Defined in differentiation/ad/ad_number_traits.h
652  template <typename T>
653  struct NumberType;
654  } // namespace internal
655 
656  // Defined in differentiation/ad/ad_number_traits.h
657  template <typename NumberType>
658  struct is_ad_number;
659  } // namespace AD
660 } // namespace Differentiation
661 
662 
663 namespace internal
664 {
669  template <typename From, typename To>
671  {
672  // Source: https://stackoverflow.com/a/16944130
673  private:
674  template <typename T>
675  static void f(T);
676 
677  template <typename F, typename T>
678  static constexpr auto
679  test(int) -> decltype(f(static_cast<T>(std::declval<F>())), true)
680  {
681  return true;
682  }
683 
684  template <typename F, typename T>
685  static constexpr auto
686  test(...) -> bool
687  {
688  return false;
689  }
690 
691  public:
692  static const bool value = test<From, To>(0);
693  };
694 
695  /*
696  * The structs below are needed to convert between some special number types.
697  * Also see tensor.h for another specialization.
698  */
699  template <typename T>
700  struct NumberType
701  {
702  static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T &
703  value(const T &t)
704  {
705  return t;
706  }
707 
708  // Below are generic functions that allows an overload for any
709  // type U that is transformable to type T. This is particularly
710  // useful when needing to cast exotic number types
711  // (e.g. auto-differentiable or symbolic numbers) to a floating
712  // point one, such as might happen when converting between tensor
713  // types.
714 
715  // Type T is constructible from F.
716  template <typename F>
717  static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T
718  value(const F &f,
719  std::enable_if_t<!std::is_same_v<std::decay_t<T>, std::decay_t<F>> &&
720  std::is_constructible_v<T, F>> * = nullptr)
721  {
722  return T(f);
723  }
724 
725  // Type T is explicitly convertible (but not constructible) from F.
726  template <typename F>
727  static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T
728  value(const F &f,
729  std::enable_if_t<!std::is_same_v<std::decay_t<T>, std::decay_t<F>> &&
730  !std::is_constructible_v<T, F> &&
732  nullptr)
733  {
734  return static_cast<T>(f);
735  }
736 
737  // Sacado doesn't provide any conversion operators, so we have
738  // to extract the value and perform further conversions from there.
739  // To be safe, we extend this to other possible AD numbers that
740  // might fall into the same category.
741  template <typename F>
742  static T
744  const F &f,
745  std::enable_if_t<!std::is_same_v<std::decay_t<T>, std::decay_t<F>> &&
746  !std::is_constructible_v<T, F> &&
749  {
751  }
752  };
753 
754  template <typename T>
755  struct NumberType<std::complex<T>>
756  {
757  static constexpr const std::complex<T> &
758  value(const std::complex<T> &t)
759  {
760  return t;
761  }
762 
763  static constexpr std::complex<T>
764  value(const T &t)
765  {
766  return std::complex<T>(t);
767  }
768 
769  // Facilitate cast from complex<double> to complex<float>
770  template <typename U>
771  static constexpr std::complex<T>
772  value(const std::complex<U> &t)
773  {
774  return std::complex<T>(NumberType<T>::value(t.real()),
775  NumberType<T>::value(t.imag()));
776  }
777  };
778 
779 #ifdef DEAL_II_WITH_CUDA
780  template <>
781  struct NumberType<cuComplex>
782  {
783  static cuComplex
784  value(const float t)
785  {
786  return make_cuComplex(t, 0.f);
787  }
788  };
789 
790  template <>
791  struct NumberType<cuDoubleComplex>
792  {
793  static cuDoubleComplex
794  value(const double t)
795  {
796  return make_cuDoubleComplex(t, 0.);
797  }
798  };
799 #endif
800 } // namespace internal
801 
802 namespace numbers
803 {
804 #ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
805 
816  // Defined in differentiation/ad/adolc_number_types.cc
817  bool
818  values_are_equal(const adouble &value_1, const adouble &value_2);
819 
820 
831  template <typename Number>
832  bool
833  values_are_equal(const adouble &value_1, const Number &value_2)
834  {
835  // Use the specialized definition for two ADOL-C taped types
836  return values_are_equal(
837  value_1, ::internal::NumberType<adouble>::value(value_2));
838  }
839 
840 
851  template <typename Number>
852  bool
853  values_are_equal(const Number &value_1, const adouble &value_2)
854  {
855  // Use the above definition
856  return values_are_equal(value_2, value_1);
857  }
858 
870  // Defined in differentiation/ad/adolc_number_types.cc
871  bool
872  value_is_less_than(const adouble &value_1, const adouble &value_2);
873 
874 
886  template <typename Number>
887  bool
888  value_is_less_than(const adouble &value_1, const Number &value_2)
889  {
890  // Use the specialized definition for two ADOL-C taped types
891  return value_is_less_than(
892  value_1, ::internal::NumberType<adouble>::value(value_2));
893  }
894 
895 
907  template <typename Number>
908  bool
909  value_is_less_than(const Number &value_1, const adouble &value_2)
910  {
911  // Use the specialized definition for two ADOL-C taped types
912  return value_is_less_than(
913  ::internal::NumberType<adouble>::value(value_1), value_2);
914  }
915 
916 #endif
917 
918 
919  template <typename Number1, typename Number2>
920  constexpr DEAL_II_HOST_DEVICE bool
921  values_are_equal(const Number1 &value_1, const Number2 &value_2)
922  {
923  return (value_1 == ::internal::NumberType<Number1>::value(value_2));
924  }
925 
926 
927  template <typename Number1, typename Number2>
928  inline bool
929  values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
930  {
931  return !(values_are_equal(value_1, value_2));
932  }
933 
934 
935  template <typename Number>
936  constexpr DEAL_II_HOST_DEVICE bool
937  value_is_zero(const Number &value)
938  {
939  return values_are_equal(value, 0.0);
940  }
941 
942 
943  template <typename Number1, typename Number2>
944  inline bool
945  value_is_less_than(const Number1 &value_1, const Number2 &value_2)
946  {
947  return (value_1 < ::internal::NumberType<Number1>::value(value_2));
948  }
949 
950 
951  template <typename Number1, typename Number2>
952  inline bool
953  value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
954  {
955  return (value_is_less_than(value_1, value_2) ||
956  values_are_equal(value_1, value_2));
957  }
958 
959 
960  template <typename Number1, typename Number2>
961  bool
962  value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
963  {
964  return !(value_is_less_than_or_equal_to(value_1, value_2));
965  }
966 
967 
968  template <typename Number1, typename Number2>
969  inline bool
970  value_is_greater_than_or_equal_to(const Number1 &value_1,
971  const Number2 &value_2)
972  {
973  return !(value_is_less_than(value_1, value_2));
974  }
975 } // namespace numbers
976 
978 
979 #endif
VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &x)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static const char T
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
Definition: numbers.h:54
static constexpr double LOG10E
Definition: numbers.h:245
static constexpr double PI_2
Definition: numbers.h:265
bool value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:953
static constexpr double E
Definition: numbers.h:235
static constexpr double PI
Definition: numbers.h:260
bool value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:962
static constexpr double SQRT2
Definition: numbers.h:275
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:937
bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:929
static constexpr double SQRT1_2
Definition: numbers.h:280
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:921
static constexpr double PI_4
Definition: numbers.h:270
static constexpr double LN10
Definition: numbers.h:255
static constexpr double LN2
Definition: numbers.h:250
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:945
bool is_finite(const double x)
Definition: numbers.h:539
static constexpr double LOG2E
Definition: numbers.h:240
bool value_is_greater_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:970
bool is_nan(const double x)
Definition: numbers.h:531
#define DEAL_II_HOST_DEVICE
Definition: numbers.h:35
#define DEAL_II_HOST_DEVICE_ALWAYS_INLINE
Definition: numbers.h:37
static cuComplex value(const float t)
Definition: numbers.h:784
static cuDoubleComplex value(const double t)
Definition: numbers.h:794
static constexpr std::complex< T > value(const std::complex< U > &t)
Definition: numbers.h:772
static constexpr std::complex< T > value(const T &t)
Definition: numbers.h:764
static constexpr const std::complex< T > & value(const std::complex< T > &t)
Definition: numbers.h:758
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T value(const F &f, std::enable_if_t<!std::is_same_v< std::decay_t< T >, std::decay_t< F >> &&!std::is_constructible_v< T, F > &&is_explicitly_convertible< const F, T >::value > *=nullptr)
Definition: numbers.h:728
static T value(const F &f, std::enable_if_t<!std::is_same_v< std::decay_t< T >, std::decay_t< F >> &&!std::is_constructible_v< T, F > &&!is_explicitly_convertible< const F, T >::value &&Differentiation::AD::is_ad_number< F >::value > *=nullptr)
Definition: numbers.h:743
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T value(const F &f, std::enable_if_t<!std::is_same_v< std::decay_t< T >, std::decay_t< F >> &&std::is_constructible_v< T, F >> *=nullptr)
Definition: numbers.h:718
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition: numbers.h:703
constexpr static unsigned int max_width
Definition: numbers.h:85
static constexpr auto test(...) -> bool
Definition: numbers.h:686
static constexpr auto test(int) -> decltype(f(static_cast< T >(std::declval< F >())), true)
Definition: numbers.h:679
static constexpr const number & conjugate(const number &x)
Definition: numbers.h:576
static constexpr bool is_complex
Definition: numbers.h:434
static real_type abs(const number &x)
Definition: numbers.h:594
static constexpr real_type abs_square(const number &x)
Definition: numbers.h:585