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