diff --git a/cmake/modules/DumuxMacros.cmake b/cmake/modules/DumuxMacros.cmake index d89484b592551353e431229e845e437df7bef7b5..45e5dd894d2581f2b8e5b2b46b2d0bfc65399944 100644 --- a/cmake/modules/DumuxMacros.cmake +++ b/cmake/modules/DumuxMacros.cmake @@ -12,3 +12,4 @@ find_package(NLOPT) find_package(PTScotch) find_package(PVPython) find_package(Valgrind) +find_package(Quadmath) diff --git a/cmake/modules/FindQuadmath.cmake b/cmake/modules/FindQuadmath.cmake new file mode 100644 index 0000000000000000000000000000000000000000..874d8efe67677f69994e862de2665ec6fd6ae339 --- /dev/null +++ b/cmake/modules/FindQuadmath.cmake @@ -0,0 +1,48 @@ +# Module that checks whether the compiler supports the +# quadruple precision floating point math +# +# Sets the following variables: +# HAVE_QUAD +# QUADMATH_LIBRARIES +# +# perform tests +include(CheckCSourceCompiles) +include(CheckCXXSourceCompiles) +include(CMakePushCheckState) +include(CheckCXXCompilerFlag) + +if(NOT DEFINED USE_QUADMATH OR USE_QUADMATH) + if(NOT DEFINED HAVE_EXTENDED_NUMERIC_LITERALS) + check_cxx_compiler_flag("-Werror -fext-numeric-literals" HAVE_EXTENDED_NUMERIC_LITERALS) + endif() + + if (HAVE_EXTENDED_NUMERIC_LITERALS) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fext-numeric-literals") + endif() + + cmake_push_check_state(RESET) + list(APPEND CMAKE_REQUIRED_LIBRARIES "quadmath") + CHECK_CXX_SOURCE_COMPILES(" +#include <quadmath.h> +int main(void){ + __float128 foo = sqrtq(123.456); + foo = FLT128_MIN; +}" QUADMATH_FOUND) + cmake_pop_check_state() + + if (QUADMATH_FOUND) + set(QUADMATH_LIBRARIES "quadmath") + set(HAVE_QUAD "${QUADMATH_FOUND}") + endif() +endif() + +if (USE_QUADMATH AND NOT QUADMATH_FOUND) + message(FATAL_ERROR "Quadruple precision math support was explicitly requested but is unavailable!") +endif() + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(Quadmath + DEFAULT_MSG + QUADMATH_LIBRARIES + HAVE_QUAD + ) diff --git a/config.h.cmake b/config.h.cmake index 859fa68708b447269b138d5ee9e3f825a91c5db6..bad26002e566f3a0c89e128936dbb9401235cc23 100644 --- a/config.h.cmake +++ b/config.h.cmake @@ -64,6 +64,9 @@ /* Define the path to pvpython */ #define PVPYTHON_EXECUTABLE "${PVPYTHON_EXECUTABLE}" +/* Define to 1 if quadmath was found */ +#cmakedefine HAVE_QUAD 1 + /* end dumux Everything below here will be overwritten */ diff --git a/dumux/common/quad.hh b/dumux/common/quad.hh index f0a50fc500b38f4081e39d112071e92137b4bfe3..00118db02a753fcf96c51522473e6fc6d5ce7df9 100644 --- a/dumux/common/quad.hh +++ b/dumux/common/quad.hh @@ -33,34 +33,52 @@ extern "C" { #include <quadmath.h> } -using quad = __float128; +#include <dune/common/typetraits.hh> -namespace std -{ +namespace Dumux { + +using Quad = __float128; + +} // namespace Dumux + +// Dune's desired way of enabling their algebraic operations for +// extended-precision data types, see dune/common/typetraits.hh. +namespace Dune { + +template <> +struct IsNumber<Dumux::Quad> : std::true_type {}; + +} // namespace Dune + +namespace std { // provide the numeric limits for the quad precision type template <> -class numeric_limits<quad> +class numeric_limits<Dumux::Quad> { public: static constexpr bool is_specialized = true; - static constexpr quad min() throw() + static constexpr Dumux::Quad min() noexcept { return FLT128_MIN; } - static constexpr quad max() throw() + static constexpr Dumux::Quad max() noexcept { return FLT128_MAX; } + static constexpr Dumux::Quad lowest() noexcept + { return -FLT128_MAX; } // number of bits in mantissa - static constexpr int digits = FLT128_MANT_DIG; + static constexpr int digits = FLT128_MANT_DIG; // number of decimal digits - static constexpr int digits10 = FLT128_DIG; + static constexpr int digits10 = FLT128_DIG; + static constexpr int max_digits10 = FLT128_MANT_DIG; + static constexpr bool is_signed = true; static constexpr bool is_integer = false; static constexpr bool is_exact = false; static constexpr int radix = 0; - static constexpr quad epsilon() throw() + static constexpr Dumux::Quad epsilon() noexcept { return FLT128_EPSILON; } - static constexpr quad round_error() throw() + static constexpr Dumux::Quad round_error() noexcept { return 0.5; } static constexpr int min_exponent = FLT128_MIN_EXP; @@ -73,13 +91,13 @@ public: static constexpr bool has_signaling_NaN = true; static constexpr float_denorm_style has_denorm = denorm_present; static constexpr bool has_denorm_loss = false; - static constexpr quad infinity() throw() + static constexpr Dumux::Quad infinity() noexcept { return __builtin_huge_valq(); }; - static constexpr quad quiet_NaN() throw() + static constexpr Dumux::Quad quiet_NaN() noexcept { return __builtin_nan(""); } - static constexpr quad signaling_NaN() throw() + static constexpr Dumux::Quad signaling_NaN() noexcept { return __builtin_nans(""); } - static constexpr quad denorm_min() throw() + static constexpr Dumux::Quad denorm_min() noexcept { return FLT128_DENORM_MIN; } static constexpr bool is_iec559 = true; @@ -91,131 +109,81 @@ public: static constexpr float_round_style round_style = round_to_nearest; }; -inline std::ostream& operator<<(std::ostream& os, const quad &val) +} // namespace std + +// Putting the following definitions in the namespace std yields undefined +// behavior. Unfortunately, ADL doesn't work for aliases (Quad) because the underlying type +// is used for the lookup. Putting these functions into any other namespace yields compiler errors. +namespace std { + +inline ostream& operator<<(ostream& os, const Dumux::Quad &val) { return (os << double(val)); } -inline std::istream& operator>>(std::istream& is, quad &val) +inline istream& operator>>(istream& is, Dumux::Quad &val) { double tmp; - std::istream &ret = (is >> tmp); + istream &ret = (is >> tmp); val = tmp; return ret; } -inline quad abs(quad val) +inline Dumux::Quad abs(Dumux::Quad val) { return (val < 0)?-val:val; } -inline quad floor(quad val) +inline Dumux::Quad floor(Dumux::Quad val) { return floorq(val); } -inline quad ceil(quad val) +inline Dumux::Quad ceil(Dumux::Quad val) { return ceilq(val); } -inline quad max(quad a, quad b) +inline Dumux::Quad max(Dumux::Quad a, Dumux::Quad b) { return (a>b)?a:b; } -inline quad min(quad a, quad b) +inline Dumux::Quad min(Dumux::Quad a, Dumux::Quad b) { return (a<b)?a:b; } -inline quad sqrt(quad val) +inline Dumux::Quad sqrt(Dumux::Quad val) { return sqrtq(val); } template <class ExpType> -inline quad pow(quad base, ExpType exp) +inline Dumux::Quad pow(Dumux::Quad base, ExpType exp) { return powq(base, exp); } -inline quad exp(quad val) +inline Dumux::Quad exp(Dumux::Quad val) { return expq(val); } -inline quad log(quad val) +inline Dumux::Quad log(Dumux::Quad val) { return logq(val); } -inline quad sin(quad val) +inline Dumux::Quad sin(Dumux::Quad val) { return sinq(val); } -inline quad cos(quad val) +inline Dumux::Quad cos(Dumux::Quad val) { return cosq(val); } -inline quad tan(quad val) +inline Dumux::Quad tan(Dumux::Quad val) { return tanq(val); } -inline quad atan(quad val) +inline Dumux::Quad atan(Dumux::Quad val) { return atanq(val); } -inline quad atan2(quad a, quad b) +inline Dumux::Quad atan2(Dumux::Quad a, Dumux::Quad b) { return atan2q(a, b); } -inline bool isfinite(quad val) +inline bool signbit(Dumux::Quad val) +{ return signbitq(val); } + +inline bool isfinite(Dumux::Quad val) { return !isnanq(val) && !isinfq(val); } -inline bool isnan(quad val) +inline bool isnan(Dumux::Quad val) { return isnanq(val); } -inline bool isinf(quad val) +inline bool isinf(Dumux::Quad val) { return isinfq(val); } } // namespace std -// this is a hack so that Dune's classname.hh does not get -// included. this is required because GCC 4.6 and earlier does not -// generate a type_info structure for __float128 which causes the -// linker to fail. -#define DUNE_CLASSNAME_HH - -#include <cstdlib> -#include <string> -#include <typeinfo> - -#if defined(__GNUC__) && ! defined(__clang__) -#include <cxxabi.h> -#endif // #ifdef __GNUC__ - -namespace Dune -{ -template <class T> -class ClassNameHelper_ -{ public: - static std::string name() - { - std::string className = typeid( T ).name(); -#if defined(__GNUC__) && ! defined(__clang__) - int status; - char *demangled = abi::__cxa_demangle( className.c_str(), 0, 0, &status ); - if( demangled ) - { - className = demangled; - std::free( demangled ); - } -#endif // #ifdef __GNUC__ - return className; - } -}; - -#if HAVE_QUAD -// specialize for quad precision floating point values to avoid -// needing a type_info structure -template <> -class ClassNameHelper_<__float128> -{ public: - static std::string name() - { return "quad"; } -}; -#endif - -template <class T> -std::string className() -{ - return ClassNameHelper_<T>::name(); -} - -template <class T> -std::string className(const T &t) -{ - return ClassNameHelper_<T>::name(); -} - -} - #endif // DUMUX_QUAD_HH diff --git a/test/porousmediumflow/1p/implicit/incompressible/CMakeLists.txt b/test/porousmediumflow/1p/implicit/incompressible/CMakeLists.txt index 457e1b19802f4152b596505a7f731b34a153d475..3fd0f237001dd1fb1c1d8d8c76def2b138e7fdba 100644 --- a/test/porousmediumflow/1p/implicit/incompressible/CMakeLists.txt +++ b/test/porousmediumflow/1p/implicit/incompressible/CMakeLists.txt @@ -10,6 +10,18 @@ dune_add_test(NAME test_1p_incompressible_tpfa ${CMAKE_CURRENT_BINARY_DIR}/1ptestcctpfa-00001.vtu --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_incompressible_tpfa test_1p.input -Problem.Name 1ptestcctpfa") +# using tpfa, analytical Jacobian and quad precision +dune_add_test(NAME test_1p_incompressible_tpfa_quad + SOURCES test_1pfv.cc + CMAKE_GUARD HAVE_QUAD + COMPILE_DEFINITIONS TYPETAG=OnePIncompressibleTpfaQuad NUMDIFFMETHOD=DiffMethod::analytic + LINK_LIBRARIES ${QUADMATH_LIBRARIES} + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/1ptestcc-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/1ptestcctpfaquad-00001.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_incompressible_tpfa_quad test_1p.input -Problem.Name 1ptestcctpfaquad") + # using mpfa and analytical Jacobian dune_add_test(NAME test_1p_incompressible_mpfa SOURCES test_1pfv.cc diff --git a/test/porousmediumflow/1p/implicit/incompressible/problem.hh b/test/porousmediumflow/1p/implicit/incompressible/problem.hh index 03bd696b0bdd523e71c0b78970d78b5d4f200beb..f7d379dc559ff4cd2e00900ba744275fdefa8349 100644 --- a/test/porousmediumflow/1p/implicit/incompressible/problem.hh +++ b/test/porousmediumflow/1p/implicit/incompressible/problem.hh @@ -74,6 +74,11 @@ SET_BOOL_PROP(OnePIncompressible, EnableGridVolumeVariablesCache, false); SET_BOOL_PROP(OnePIncompressible, EnableGridFluxVariablesCache, false); SET_BOOL_PROP(OnePIncompressible, EnableFVGridGeometryCache, false); +// define a TypeTag for a quad precision test +#if HAVE_QUAD +NEW_TYPE_TAG(OnePIncompressibleTpfaQuad, INHERITS_FROM(OnePIncompressibleTpfa)); +SET_TYPE_PROP(OnePIncompressibleTpfaQuad, Scalar, Quad); +#endif } // end namespace Properties /*! * \ingroup OnePTests diff --git a/test/porousmediumflow/1p/implicit/incompressible/test_1pfv.cc b/test/porousmediumflow/1p/implicit/incompressible/test_1pfv.cc index 5380c2f0fa8c2ef2e62f062ee217b09048e3e304..35eaadff20c7091c66f78d5df38238d1ee9610e1 100644 --- a/test/porousmediumflow/1p/implicit/incompressible/test_1pfv.cc +++ b/test/porousmediumflow/1p/implicit/incompressible/test_1pfv.cc @@ -26,6 +26,9 @@ #include <ctime> #include <iostream> +// Support for quad precision has to be included before any other Dune module: +#include <dumux/common/quad.hh> + #include <dune/common/parallel/mpihelper.hh> #include <dune/common/timer.hh> #include <dune/grid/io/file/dgfparser/dgfexception.hh>