From 88311a066f72e0c587aa79d2f721b805d507b85e Mon Sep 17 00:00:00 2001 From: David Wells Date: Wed, 3 Aug 2022 11:39:02 -0400 Subject: [PATCH] Don't include Sacado.hpp in base/numbers.h. Sacado.hpp, when preprocessed, is over 200,000 lines. We can get things working in an equivalent way by using argument dependent lookup instead. --- include/deal.II/base/numbers.h | 17 +-- include/deal.II/base/symmetric_tensor.h | 38 +++--- include/deal.II/base/tensor.h | 8 +- .../deal.II/differentiation/ad/sacado_math.h | 127 +++++++++--------- 4 files changed, 91 insertions(+), 99 deletions(-) diff --git a/include/deal.II/base/numbers.h b/include/deal.II/base/numbers.h index 74d8976ee2..292816d2dd 100644 --- a/include/deal.II/base/numbers.h +++ b/include/deal.II/base/numbers.h @@ -137,15 +137,6 @@ DEAL_II_NAMESPACE_CLOSE # include // Taped double #endif -// Ideally we'd like to #include -// but header indirectly references numbers.h. We therefore simply -// import the whole Sacado header at this point to get the math -// functions imported into the standard namespace. -#ifdef DEAL_II_TRILINOS_WITH_SACADO -DEAL_II_DISABLE_EXTRA_DIAGNOSTICS -# include -DEAL_II_ENABLE_EXTRA_DIAGNOSTICS -#endif namespace std { @@ -608,7 +599,9 @@ namespace numbers typename NumberTraits::real_type NumberTraits::abs(const number &x) { - return std::abs(x); + // Make things work with AD types + using std::abs; + return abs(x); } @@ -626,7 +619,9 @@ namespace numbers typename NumberTraits>::real_type NumberTraits>::abs(const std::complex &x) { - return std::abs(x); + // Make things work with AD types + using std::abs; + return abs(x); } diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h index d69798c9a6..6396f05d58 100644 --- a/include/deal.II/base/symmetric_tensor.h +++ b/include/deal.II/base/symmetric_tensor.h @@ -1563,7 +1563,7 @@ namespace internal // whether we have to fear that the matrix is not regular. Number diagonal_sum = internal::NumberType::value(0.0); for (unsigned int i = 0; i < N; ++i) - diagonal_sum += std::fabs(tmp.data[i][i]); + diagonal_sum += numbers::NumberTraits::abs(tmp.data[i][i]); const Number typical_diagonal_element = diagonal_sum / static_cast(N); (void)typical_diagonal_element; @@ -1576,12 +1576,12 @@ namespace internal { // Pivot search: search that part of the line on and right of the // diagonal for the largest element. - Number max = std::fabs(tmp.data[j][j]); - unsigned int r = j; + Number max = numbers::NumberTraits::abs(tmp.data[j][j]); + unsigned int r = j; for (unsigned int i = j + 1; i < N; ++i) - if (std::fabs(tmp.data[i][j]) > max) + if (numbers::NumberTraits::abs(tmp.data[i][j]) > max) { - max = std::fabs(tmp.data[i][j]); + max = numbers::NumberTraits::abs(tmp.data[i][j]); r = i; } @@ -2343,25 +2343,25 @@ namespace internal compute_norm(const typename SymmetricTensorAccessors:: StorageType<2, dim, Number>::base_tensor_type &data) { + // Make things work with AD types + using std::sqrt; switch (dim) { case 1: return numbers::NumberTraits::abs(data[0]); case 2: - return std::sqrt( - numbers::NumberTraits::abs_square(data[0]) + - numbers::NumberTraits::abs_square(data[1]) + - 2. * numbers::NumberTraits::abs_square(data[2])); + return sqrt(numbers::NumberTraits::abs_square(data[0]) + + numbers::NumberTraits::abs_square(data[1]) + + 2. * numbers::NumberTraits::abs_square(data[2])); case 3: - return std::sqrt( - numbers::NumberTraits::abs_square(data[0]) + - numbers::NumberTraits::abs_square(data[1]) + - numbers::NumberTraits::abs_square(data[2]) + - 2. * numbers::NumberTraits::abs_square(data[3]) + - 2. * numbers::NumberTraits::abs_square(data[4]) + - 2. * numbers::NumberTraits::abs_square(data[5])); + return sqrt(numbers::NumberTraits::abs_square(data[0]) + + numbers::NumberTraits::abs_square(data[1]) + + numbers::NumberTraits::abs_square(data[2]) + + 2. * numbers::NumberTraits::abs_square(data[3]) + + 2. * numbers::NumberTraits::abs_square(data[4]) + + 2. * numbers::NumberTraits::abs_square(data[5])); default: { @@ -2375,7 +2375,7 @@ namespace internal return_value += 2. * numbers::NumberTraits::abs_square(data[d]); - return std::sqrt(return_value); + return sqrt(return_value); } } } @@ -2387,6 +2387,8 @@ namespace internal compute_norm(const typename SymmetricTensorAccessors:: StorageType<4, dim, Number>::base_tensor_type &data) { + // Make things work with AD types + using std::sqrt; switch (dim) { case 1: @@ -2416,7 +2418,7 @@ namespace internal return_value += 4. * numbers::NumberTraits::abs_square(data[i][j]); - return std::sqrt(return_value); + return sqrt(return_value); } } } diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h index 4cde9bde18..b382b2b401 100644 --- a/include/deal.II/base/tensor.h +++ b/include/deal.II/base/tensor.h @@ -1727,7 +1727,9 @@ template inline typename numbers::NumberTraits::real_type Tensor::norm() const { - return std::sqrt(norm_square()); + // Make things work with AD types + using std::sqrt; + return sqrt(norm_square()); } @@ -3009,7 +3011,7 @@ l1_norm(const Tensor<2, dim, Number> &t) { Number sum = internal::NumberType::value(0.0); for (unsigned int i = 0; i < dim; ++i) - sum += std::fabs(t[i][j]); + sum += numbers::NumberTraits::abs(t[i][j]); if (sum > max) max = sum; @@ -3035,7 +3037,7 @@ linfty_norm(const Tensor<2, dim, Number> &t) { Number sum = internal::NumberType::value(0.0); for (unsigned int j = 0; j < dim; ++j) - sum += std::fabs(t[i][j]); + sum += numbers::NumberTraits::abs(t[i][j]); if (sum > max) max = sum; diff --git a/include/deal.II/differentiation/ad/sacado_math.h b/include/deal.II/differentiation/ad/sacado_math.h index 655d281ce7..bd9d8cde93 100644 --- a/include/deal.II/differentiation/ad/sacado_math.h +++ b/include/deal.II/differentiation/ad/sacado_math.h @@ -28,76 +28,69 @@ # ifndef DOXYGEN +DEAL_II_NAMESPACE_OPEN /** - * Define some missing fundamental math functions + * Implementation of the error function for real-valued Sacado numbers. */ -namespace std +template < + typename ADNumberType, + typename = std::enable_if_t< + dealii::Differentiation::AD::is_sacado_number::value && + dealii::Differentiation::AD::is_real_valued_ad_number::value>> +inline ADNumberType +erf(ADNumberType x) { - /** - * Implementation of the error function for real-valued Sacado numbers. - */ - template < - typename ADNumberType, - typename = std::enable_if_t< - dealii::Differentiation::AD::is_sacado_number::value && - dealii::Differentiation::AD::is_real_valued_ad_number< - ADNumberType>::value>> - inline ADNumberType - erf(ADNumberType x) - { - // Reference: - // Handbook of Mathematical Functions: with Formulas, Graphs, and - // Mathematical Tables Abramowitz, M. and Stegun, I. Dover Books on - // Mathematics. 1972. - // - // Current implementation source: - // https://www.johndcook.com/blog/cpp_erf/ - // https://www.johndcook.com/blog/2009/01/19/stand-alone-error-function-erf/ - // - // Note: This implementation has a reported maximum round-off error - // of 1.5e-7. - - // Constants - const double a1 = 0.254829592; - const double a2 = -0.284496736; - const double a3 = 1.421413741; - const double a4 = -1.453152027; - const double a5 = 1.061405429; - const double p = 0.3275911; - - // Save the sign of x - const bool neg_val = - (x < dealii::internal::NumberType::value(0.0) ? true : - false); - x = std::fabs(x); - - // Abramowitz1972a equation 7.1.26 - const ADNumberType t = 1.0 / (1.0 + p * x); - const ADNumberType y = - 1.0 - - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * std::exp(-x * x); - - if (!neg_val) - return y; - else - return -y; - } - - - /** - * Implementation of the complementary error function for Sacado numbers. - */ - template < - typename ADNumberType, - typename = std::enable_if_t< - dealii::Differentiation::AD::is_sacado_number::value>> - inline ADNumberType - erfc(const ADNumberType &x) - { - return 1.0 - std::erf(x); - } - -} // namespace std + // Reference: + // Handbook of Mathematical Functions: with Formulas, Graphs, and + // Mathematical Tables Abramowitz, M. and Stegun, I. Dover Books on + // Mathematics. 1972. + // + // Current implementation source: + // https://www.johndcook.com/blog/cpp_erf/ + // https://www.johndcook.com/blog/2009/01/19/stand-alone-error-function-erf/ + // + // Note: This implementation has a reported maximum round-off error + // of 1.5e-7. + + // Constants + const double a1 = 0.254829592; + const double a2 = -0.284496736; + const double a3 = 1.421413741; + const double a4 = -1.453152027; + const double a5 = 1.061405429; + const double p = 0.3275911; + + // Save the sign of x + const bool neg_val = + (x < dealii::internal::NumberType::value(0.0) ? true : false); + x = std::fabs(x); + + // Abramowitz1972a equation 7.1.26 + const ADNumberType t = 1.0 / (1.0 + p * x); + const ADNumberType y = + 1.0 - + (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * std::exp(-x * x); + + if (!neg_val) + return y; + else + return -y; +} + + +/** + * Implementation of the complementary error function for Sacado numbers. + */ +template ::value>> +inline ADNumberType +erfc(const ADNumberType &x) +{ + return 1.0 - std::erf(x); +} + +DEAL_II_NAMESPACE_CLOSE # endif // DOXYGEN -- 2.39.5