From: Jean-Paul Pelteret Date: Thu, 16 Nov 2017 07:38:25 +0000 (+0100) Subject: Add Sacado math header X-Git-Tag: v9.0.0-rc1~764^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d0e9d0b454c6f59f2ae39313ac1b47a14c811b2d;p=dealii.git Add Sacado math header This adds some missing fundamental math functions for the supported Sacado numbers. --- diff --git a/doc/doxygen/headers/automatic_and_symbolic_differentiation.h b/doc/doxygen/headers/automatic_and_symbolic_differentiation.h index 9b6195edb8..e54972b15d 100644 --- a/doc/doxygen/headers/automatic_and_symbolic_differentiation.h +++ b/doc/doxygen/headers/automatic_and_symbolic_differentiation.h @@ -73,6 +73,8 @@ * use Adol-C numbers. * - adolc_product_types.h: Defines some product and scalar types that allow the use of * Adol-C numbers in conjunction with the Tensor and SymmetricTensor classes. + * - sacado_math.h: Extension of the sacado math operations that allow these numbers to be used + * consistently throughout the library. * - sacado_number_types.h: Implementation of the internal classes that define how we * use the supported Sacado numbers. * - sacado_product_types.h: Defines some product and scalar types that allow the use of diff --git a/include/deal.II/differentiation/ad.h b/include/deal.II/differentiation/ad.h index ccfbc517b6..74a9421407 100644 --- a/include/deal.II/differentiation/ad.h +++ b/include/deal.II/differentiation/ad.h @@ -27,6 +27,7 @@ #include #include +#include #include #include diff --git a/include/deal.II/differentiation/ad/sacado_math.h b/include/deal.II/differentiation/ad/sacado_math.h new file mode 100644 index 0000000000..3b2598713f --- /dev/null +++ b/include/deal.II/differentiation/ad/sacado_math.h @@ -0,0 +1,96 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii_differentiation_ad_sacado_math_h +#define dealii_differentiation_ad_sacado_math_h + +#include + +#ifdef DEAL_II_WITH_TRILINOS + +#include +#include +#include + + +#ifndef DOXYGEN + +/** + * Define some missing fundamental math functions + */ +namespace std +{ + + /** + * Implementation of the error function for real-valued Sacado numbers. + */ + template::value && + dealii::Differentiation::AD::is_real_valued_ad_number::value + >::type> + 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::value + >::type> + inline ADNumberType erfc(const ADNumberType &x) + { + return 1.0 - std::erf(x); + } + +} // namespace std + +#endif // DOXYGEN + +#endif // DEAL_II_WITH_TRILINOS + +#endif