Reference documentation for deal.II version GIT 276cfaa7ee 2022-08-16 14:00: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\}}\)
sacado_math.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2018 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_differentiation_ad_sacado_math_h
17 #define dealii_differentiation_ad_sacado_math_h
18 
19 #include <deal.II/base/config.h>
20 
21 #ifdef DEAL_II_TRILINOS_WITH_SACADO
22 
23 # include <deal.II/base/numbers.h>
24 
27 
28 
29 # ifndef DOXYGEN
30 
35 template <
36  typename ADNumberType,
37  typename = std::enable_if_t<
40 inline ADNumberType
41 erf(ADNumberType x)
42 {
43  // Reference:
44  // Handbook of Mathematical Functions: with Formulas, Graphs, and
45  // Mathematical Tables Abramowitz, M. and Stegun, I. Dover Books on
46  // Mathematics. 1972.
47  //
48  // Current implementation source:
49  // https://www.johndcook.com/blog/cpp_erf/
50  // https://www.johndcook.com/blog/2009/01/19/stand-alone-error-function-erf/
51  //
52  // Note: This implementation has a reported maximum round-off error
53  // of 1.5e-7.
54 
55  // Constants
56  const double a1 = 0.254829592;
57  const double a2 = -0.284496736;
58  const double a3 = 1.421413741;
59  const double a4 = -1.453152027;
60  const double a5 = 1.061405429;
61  const double p = 0.3275911;
62 
63  // Save the sign of x
64  const bool neg_val =
65  (x < ::internal::NumberType<ADNumberType>::value(0.0) ? true : false);
66  x = std::fabs(x);
67 
68  // Abramowitz1972a equation 7.1.26
69  const ADNumberType t = 1.0 / (1.0 + p * x);
70  const ADNumberType y =
71  1.0 -
72  (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * std::exp(-x * x);
73 
74  if (!neg_val)
75  return y;
76  else
77  return -y;
78 }
79 
80 
84 template <typename ADNumberType,
85  typename = std::enable_if_t<
87 inline ADNumberType
88 erfc(const ADNumberType &x)
89 {
90  return 1.0 - std::erf(x);
91 }
92 
94 
95 # endif // DOXYGEN
96 
97 #endif // DEAL_II_TRILINOS_WITH_SACADO
98 
99 #endif
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
Expression fabs(const Expression &x)
Expression erfc(const Expression &x)
Expression erf(const Expression &x)