Reference documentation for deal.II version GIT 7deb6c54a6 2023-06-09 18: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\}}\)
numbers.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2006 - 2022 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 // compiling CUDA code, i.e, when Kokkos was configured with CUDA support.
41 // This is problematic when calling non-constexpr functions in constexpr
42 // functions. Hence, we need a way to annotate functions explicitly as
43 // host-only.
44 #if defined(__clang__) && defined(__CUDA__)
45 # define DEAL_II_HOST __host__
46 #else
47 # define DEAL_II_HOST
48 #endif
49 
50 // Forward-declare the automatic differentiation types so we can add prototypes
51 // for our own wrappers.
52 #ifdef DEAL_II_WITH_ADOLC
53 class adouble;
54 namespace adtl
55 {
56  class adouble;
57 }
58 #endif
59 
61 
62 namespace internal
63 {
80  template <typename Number>
82  {
86  constexpr static unsigned int max_width = 1;
87  };
88 
95  template <>
97  {
101  constexpr static unsigned int max_width =
102 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512
103  8;
104 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256
105  4;
106 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128
107  2;
108 #else
109  1;
110 #endif
111  };
112 
119  template <>
121  {
125  constexpr static unsigned int max_width =
126 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ALTIVEC__)
127  4;
128 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512 && defined(__AVX512F__)
129  16;
130 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256 && defined(__AVX__)
131  8;
132 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__SSE2__)
133  4;
134 #else
135  1;
136 #endif
137  };
138 
139 
140 } // namespace internal
141 
142 // forward declarations to support abs or sqrt operations on VectorizedArray
143 #ifndef DOXYGEN
144 template <typename Number,
145  std::size_t width =
147 class VectorizedArray;
148 template <typename T>
149 struct EnableIfScalar;
150 #endif
151 
152 #ifdef DEAL_II_WITH_ADOLC
153 # ifndef DOXYGEN
154 // Prototype some inline functions present in adolc_math.h for use in
155 // NumberTraits.
156 //
157 // ADOL-C uses fabs(), but for genericity we want to use abs(). Simultaneously,
158 // though, we don't want to include ADOL-C headers in this header since
159 // numbers.h is in everything. To get around this: use C++ rules which permit
160 // the use of forward-declared classes in function prototypes to declare some
161 // functions which are defined in adolc_math.h. This permits us to write "using
162 // ::abs;" in NumberTraits which will allow us to select the correct
163 // overload (the one in ::) when instantiating NumberTraits for ADOL-C
164 // types.
165 
166 adouble
167 abs(const adouble &x);
168 
169 adtl::adouble
170 abs(const adtl::adouble &x);
171 # endif
172 #endif
173 
175 
176 namespace std
177 {
178  template <typename Number, std::size_t width>
179  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
180  sqrt(const ::VectorizedArray<Number, width> &);
181  template <typename Number, std::size_t width>
182  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
183  abs(const ::VectorizedArray<Number, width> &);
184  template <typename Number, std::size_t width>
185  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
186  max(const ::VectorizedArray<Number, width> &,
187  const ::VectorizedArray<Number, width> &);
188  template <typename Number, std::size_t width>
189  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
190  min(const ::VectorizedArray<Number, width> &,
191  const ::VectorizedArray<Number, width> &);
192  template <typename Number, size_t width>
194  pow(const ::VectorizedArray<Number, width> &, const Number p);
195  template <typename Number, size_t width>
197  sin(const ::VectorizedArray<Number, width> &);
198  template <typename Number, size_t width>
200  cos(const ::VectorizedArray<Number, width> &);
201  template <typename Number, size_t width>
203  tan(const ::VectorizedArray<Number, width> &);
204  template <typename Number, size_t width>
206  exp(const ::VectorizedArray<Number, width> &);
207  template <typename Number, size_t width>
209  log(const ::VectorizedArray<Number, width> &);
210 } // namespace std
211 
213 
229 namespace numbers
230 {
234  static constexpr double E = 2.7182818284590452354;
235 
239  static constexpr double LOG2E = 1.4426950408889634074;
240 
244  static constexpr double LOG10E = 0.43429448190325182765;
245 
249  static constexpr double LN2 = 0.69314718055994530942;
250 
254  static constexpr double LN10 = 2.30258509299404568402;
255 
259  static constexpr double PI = 3.14159265358979323846;
260 
264  static constexpr double PI_2 = 1.57079632679489661923;
265 
269  static constexpr double PI_4 = 0.78539816339744830962;
270 
274  static constexpr double SQRT2 = 1.41421356237309504880;
275 
279  static constexpr double SQRT1_2 = 0.70710678118654752440;
280 
290  bool
291  is_finite(const double x);
292 
297  bool
298  is_finite(const std::complex<double> &x);
299 
304  bool
305  is_finite(const std::complex<float> &x);
306 
315  bool
316  is_finite(const std::complex<long double> &x);
317 
328  template <typename Number1, typename Number2>
329  constexpr DEAL_II_HOST_DEVICE bool
330  values_are_equal(const Number1 &value_1, const Number2 &value_2);
331 
342  template <typename Number1, typename Number2>
343  bool
344  values_are_not_equal(const Number1 &value_1, const Number2 &value_2);
345 
353  template <typename Number>
354  constexpr DEAL_II_HOST_DEVICE bool
355  value_is_zero(const Number &value);
356 
367  template <typename Number1, typename Number2>
368  bool
369  value_is_less_than(const Number1 &value_1, const Number2 &value_2);
370 
381  template <typename Number1, typename Number2>
382  bool
383  value_is_less_than_or_equal_to(const Number1 &value_1,
384  const Number2 &value_2);
385 
386 
387 
398  template <typename Number1, typename Number2>
399  bool
400  value_is_greater_than(const Number1 &value_1, const Number2 &value_2);
401 
412  template <typename Number1, typename Number2>
413  bool
414  value_is_greater_than_or_equal_to(const Number1 &value_1,
415  const Number2 &value_2);
416 
425  template <typename number>
427  {
433  static constexpr bool is_complex = false;
434 
441  using real_type = number;
442 
447 
455  static constexpr DEAL_II_HOST_DEVICE const number &
456  conjugate(const number &x);
457 
466  static constexpr DEAL_II_HOST_DEVICE real_type
467  abs_square(const number &x);
468 
472  static real_type
473  abs(const number &x);
474  };
475 
476 
481  template <typename number>
482  struct NumberTraits<std::complex<number>>
483  {
489  static constexpr bool is_complex = true;
490 
497  using real_type = number;
498 
502  using double_type = std::complex<double>;
503 
507  static constexpr std::complex<number>
508  conjugate(const std::complex<number> &x);
509 
516  static constexpr real_type
517  abs_square(const std::complex<number> &x);
518 
519 
523  static real_type
524  abs(const std::complex<number> &x);
525  };
526 
527  // --------------- inline and template functions ---------------- //
528 
529  inline bool
530  is_nan(const double x)
531  {
532  return std::isnan(x);
533  }
534 
535 
536 
537  inline bool
538  is_finite(const double x)
539  {
540  return std::isfinite(x);
541  }
542 
543 
544 
545  inline bool
546  is_finite(const std::complex<double> &x)
547  {
548  // Check complex numbers for infinity
549  // by testing real and imaginary part
550  return (is_finite(x.real()) && is_finite(x.imag()));
551  }
552 
553 
554 
555  inline bool
556  is_finite(const std::complex<float> &x)
557  {
558  // Check complex numbers for infinity
559  // by testing real and imaginary part
560  return (is_finite(x.real()) && is_finite(x.imag()));
561  }
562 
563 
564 
565  inline bool
566  is_finite(const std::complex<long double> &x)
567  {
568  // Same for std::complex<long double>
569  return (is_finite(x.real()) && is_finite(x.imag()));
570  }
571 
572 
573  template <typename number>
574  constexpr DEAL_II_HOST_DEVICE const number &
576  {
577  return x;
578  }
579 
580 
581 
582  template <typename number>
585  {
586  return x * x;
587  }
588 
589 
590 
591  template <typename number>
593  NumberTraits<number>::abs(const number &x)
594  {
595  // Make things work with AD types
596  using std::abs;
597 #ifdef DEAL_II_WITH_ADOLC
598  // This one is a little tricky - we have our own abs function in ::,
599  // prototyped with forward-declared types in this file, but it only exists
600  // if we have ADOL-C: hence we only add this using statement in that
601  // situation
602  using ::abs;
603 #endif
604  return abs(x);
605  }
606 
607 
608 
609  template <typename number>
610  constexpr std::complex<number>
611  NumberTraits<std::complex<number>>::conjugate(const std::complex<number> &x)
612  {
613  return std::conj(x);
614  }
615 
616 
617 
618  template <typename number>
619  typename NumberTraits<std::complex<number>>::real_type
620  NumberTraits<std::complex<number>>::abs(const std::complex<number> &x)
621  {
622  // Make things work with AD types
623  using std::abs;
624 #ifdef DEAL_II_WITH_ADOLC
625  // Same comment as the non-complex case holds here
626  using ::abs;
627 #endif
628  return abs(x);
629  }
630 
631 
632 
633  template <typename number>
634  constexpr typename NumberTraits<std::complex<number>>::real_type
635  NumberTraits<std::complex<number>>::abs_square(const std::complex<number> &x)
636  {
637  return std::norm(x);
638  }
639 
640 } // namespace numbers
641 
642 
643 // Forward declarations
645 {
646  namespace AD
647  {
648  namespace internal
649  {
650  // Defined in differentiation/ad/ad_number_traits.h
651  template <typename T>
652  struct NumberType;
653  } // namespace internal
654 
655  // Defined in differentiation/ad/ad_number_traits.h
656  template <typename NumberType>
657  struct is_ad_number;
658  } // namespace AD
659 } // namespace Differentiation
660 
661 
662 namespace internal
663 {
668  template <typename From, typename To>
670  {
671  // Source: https://stackoverflow.com/a/16944130
672  private:
673  template <typename T>
674  static void f(T);
675 
676  template <typename F, typename T>
677  static constexpr auto
678  test(int) -> decltype(f(static_cast<T>(std::declval<F>())), true)
679  {
680  return true;
681  }
682 
683  template <typename F, typename T>
684  static constexpr auto
685  test(...) -> bool
686  {
687  return false;
688  }
689 
690  public:
691  static bool const value = test<From, To>(0);
692  };
693 
694  /*
695  * The structs below are needed to convert between some special number types.
696  * Also see tensor.h for another specialization.
697  */
698  template <typename T>
699  struct NumberType
700  {
701  static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T &
702  value(const T &t)
703  {
704  return t;
705  }
706 
707  // Below are generic functions that allows an overload for any
708  // type U that is transformable to type T. This is particularly
709  // useful when needing to cast exotic number types
710  // (e.g. auto-differentiable or symbolic numbers) to a floating
711  // point one, such as might happen when converting between tensor
712  // types.
713 
714  // Type T is constructible from F.
715  template <typename F>
716  static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T
717  value(const F &f,
718  std::enable_if_t<!std::is_same<typename std::decay<T>::type,
719  typename std::decay<F>::type>::value &&
720  std::is_constructible<T, F>::value> * = 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<typename std::decay<T>::type,
730  typename std::decay<F>::type>::value &&
731  !std::is_constructible<T, F>::value &&
733  nullptr)
734  {
735  return static_cast<T>(f);
736  }
737 
738  // Sacado doesn't provide any conversion operators, so we have
739  // to extract the value and perform further conversions from there.
740  // To be safe, we extend this to other possible AD numbers that
741  // might fall into the same category.
742  template <typename F>
743  static T
745  const F &f,
746  std::enable_if_t<!std::is_same<typename std::decay<T>::type,
747  typename std::decay<F>::type>::value &&
748  !std::is_constructible<T, F>::value &&
751  {
753  }
754  };
755 
756  template <typename T>
757  struct NumberType<std::complex<T>>
758  {
759  static constexpr const std::complex<T> &
760  value(const std::complex<T> &t)
761  {
762  return t;
763  }
764 
765  static constexpr std::complex<T>
766  value(const T &t)
767  {
768  return std::complex<T>(t);
769  }
770 
771  // Facilitate cast from complex<double> to complex<float>
772  template <typename U>
773  static constexpr std::complex<T>
774  value(const std::complex<U> &t)
775  {
776  return std::complex<T>(NumberType<T>::value(t.real()),
777  NumberType<T>::value(t.imag()));
778  }
779  };
780 
781 #ifdef DEAL_II_WITH_CUDA
782  template <>
783  struct NumberType<cuComplex>
784  {
785  static cuComplex
786  value(const float t)
787  {
788  return make_cuComplex(t, 0.f);
789  }
790  };
791 
792  template <>
793  struct NumberType<cuDoubleComplex>
794  {
795  static cuDoubleComplex
796  value(const double t)
797  {
798  return make_cuDoubleComplex(t, 0.);
799  }
800  };
801 #endif
802 } // namespace internal
803 
804 namespace numbers
805 {
806 #ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
807 
818  // Defined in differentiation/ad/adolc_number_types.cc
819  bool
820  values_are_equal(const adouble &value_1, const adouble &value_2);
821 
822 
833  template <typename Number>
834  bool
835  values_are_equal(const adouble &value_1, const Number &value_2)
836  {
837  // Use the specialized definition for two ADOL-C taped types
838  return values_are_equal(
839  value_1, ::internal::NumberType<adouble>::value(value_2));
840  }
841 
842 
853  template <typename Number>
854  bool
855  values_are_equal(const Number &value_1, const adouble &value_2)
856  {
857  // Use the above definition
858  return values_are_equal(value_2, value_1);
859  }
860 
872  // Defined in differentiation/ad/adolc_number_types.cc
873  bool
874  value_is_less_than(const adouble &value_1, const adouble &value_2);
875 
876 
888  template <typename Number>
889  bool
890  value_is_less_than(const adouble &value_1, const Number &value_2)
891  {
892  // Use the specialized definition for two ADOL-C taped types
893  return value_is_less_than(
894  value_1, ::internal::NumberType<adouble>::value(value_2));
895  }
896 
897 
909  template <typename Number>
910  bool
911  value_is_less_than(const Number &value_1, const adouble &value_2)
912  {
913  // Use the specialized definition for two ADOL-C taped types
914  return value_is_less_than(
915  ::internal::NumberType<adouble>::value(value_1), value_2);
916  }
917 
918 #endif
919 
920 
921  template <typename Number1, typename Number2>
922  constexpr DEAL_II_HOST_DEVICE bool
923  values_are_equal(const Number1 &value_1, const Number2 &value_2)
924  {
925  return (value_1 == ::internal::NumberType<Number1>::value(value_2));
926  }
927 
928 
929  template <typename Number1, typename Number2>
930  inline bool
931  values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
932  {
933  return !(values_are_equal(value_1, value_2));
934  }
935 
936 
937  template <typename Number>
938  constexpr DEAL_II_HOST_DEVICE bool
939  value_is_zero(const Number &value)
940  {
941  return values_are_equal(value, 0.0);
942  }
943 
944 
945  template <typename Number1, typename Number2>
946  inline bool
947  value_is_less_than(const Number1 &value_1, const Number2 &value_2)
948  {
949  return (value_1 < ::internal::NumberType<Number1>::value(value_2));
950  }
951 
952 
953  template <typename Number1, typename Number2>
954  inline bool
955  value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
956  {
957  return (value_is_less_than(value_1, value_2) ||
958  values_are_equal(value_1, value_2));
959  }
960 
961 
962  template <typename Number1, typename Number2>
963  bool
964  value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
965  {
966  return !(value_is_less_than_or_equal_to(value_1, value_2));
967  }
968 
969 
970  template <typename Number1, typename Number2>
971  inline bool
972  value_is_greater_than_or_equal_to(const Number1 &value_1,
973  const Number2 &value_2)
974  {
975  return !(value_is_less_than(value_1, value_2));
976  }
977 } // namespace numbers
978 
980 
981 #endif
VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &x)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:475
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:476
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:55
static constexpr double LOG10E
Definition: numbers.h:244
static constexpr double PI_2
Definition: numbers.h:264
bool value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:955
static constexpr double E
Definition: numbers.h:234
static constexpr double PI
Definition: numbers.h:259
bool value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:964
static constexpr double SQRT2
Definition: numbers.h:274
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:939
bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:931
static constexpr double SQRT1_2
Definition: numbers.h:279
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:923
static constexpr double PI_4
Definition: numbers.h:269
static constexpr double LN10
Definition: numbers.h:254
static constexpr double LN2
Definition: numbers.h:249
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:947
bool is_finite(const double x)
Definition: numbers.h:538
static constexpr double LOG2E
Definition: numbers.h:239
bool value_is_greater_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:972
bool is_nan(const double x)
Definition: numbers.h:530
#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:786
static cuDoubleComplex value(const double t)
Definition: numbers.h:796
static constexpr std::complex< T > value(const std::complex< U > &t)
Definition: numbers.h:774
static constexpr std::complex< T > value(const T &t)
Definition: numbers.h:766
static constexpr const std::complex< T > & value(const std::complex< T > &t)
Definition: numbers.h:760
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T value(const F &f, std::enable_if_t<!std::is_same< typename std::decay< T >::type, typename std::decay< F >::type >::value &&!std::is_constructible< T, F >::value &&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< typename std::decay< T >::type, typename std::decay< F >::type >::value &&!std::is_constructible< T, F >::value &&!is_explicitly_convertible< const F, T >::value &&Differentiation::AD::is_ad_number< F >::value > *=nullptr)
Definition: numbers.h:744
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T value(const F &f, std::enable_if_t<!std::is_same< typename std::decay< T >::type, typename std::decay< F >::type >::value &&std::is_constructible< T, F >::value > *=nullptr)
Definition: numbers.h:717
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition: numbers.h:702
constexpr static unsigned int max_width
Definition: numbers.h:86
static constexpr auto test(...) -> bool
Definition: numbers.h:685
static constexpr auto test(int) -> decltype(f(static_cast< T >(std::declval< F >())), true)
Definition: numbers.h:678
static constexpr const number & conjugate(const number &x)
Definition: numbers.h:575
static constexpr bool is_complex
Definition: numbers.h:433
static real_type abs(const number &x)
Definition: numbers.h:593
static constexpr real_type abs_square(const number &x)
Definition: numbers.h:584