Merge pull request #434 from mwarusz/precision

avoiding float to double promotion

Author: pdziekan

libmpdata++/concurr/detail/sharedmem.hpp

@@ -32,9 +32,8 @@ namespace libmpdataxx
 	static_assert(n_dims > 0, "n_dims <= 0");
 	static_assert(n_tlev > 0, "n_tlev <= 0");
 
-        // TODO: T_sumtype (perhaps worh using double even if summing floats?)
         std::unique_ptr<blitz::Array<real_t, 1>> xtmtmp; 
-        std::unique_ptr<blitz::Array<real_t, 1>> sumtmp;
+        std::unique_ptr<blitz::Array<double, 1>> sumtmp;
 
         protected:
 
@@ -93,12 +92,12 @@ namespace libmpdataxx
             throw std::runtime_error("number of subdomains greater than number of gridpoints");
 
           if (n_dims != 1) 
-            sumtmp.reset(new blitz::Array<real_t, 1>(grid_size[0]));
+            sumtmp.reset(new blitz::Array<double, 1>(grid_size[0]));
           xtmtmp.reset(new blitz::Array<real_t, 1>(size));
         }
 
         /// @brief concurrency-aware summation of array elements
-        real_t sum(const arr_t &arr, const idx_t<n_dims> &ijk, const bool sum_khn)
+        double sum(const arr_t &arr, const idx_t<n_dims> &ijk, const bool sum_khn)
         {
 	  // doing a two-step sum to reduce numerical error 
 	  // and make parallel results reproducible
@@ -114,7 +113,7 @@ namespace libmpdataxx
 	      (*sumtmp)(c) = blitz::sum(arr(slice_idx));
           }
           barrier();
-          real_t result;
+          double result;
           if (sum_khn)
             result = blitz::kahan_sum(*sumtmp);
           else
@@ -124,7 +123,7 @@ namespace libmpdataxx
         }
 
         /// @brief concurrency-aware summation of a (element-wise) product of two arrays
-        real_t sum(const arr_t &arr1, const arr_t &arr2, const idx_t<n_dims> &ijk, const bool sum_khn)
+        double sum(const arr_t &arr1, const arr_t &arr2, const idx_t<n_dims> &ijk, const bool sum_khn)
         {
 	  // doing a two-step sum to reduce numerical error 
 	  // and make parallel results reproducible
@@ -140,7 +139,7 @@ namespace libmpdataxx
 	      (*sumtmp)(c) = blitz::sum(arr1(slice_idx) * arr2(slice_idx)); 
           }
           barrier();
-          real_t result;
+          double result;
           if (sum_khn)
             result = blitz::kahan_sum(*sumtmp);
           else

libmpdata++/formulae/common.hpp

@@ -13,6 +13,12 @@ namespace libmpdataxx
 {
   namespace formulae
   {
+    // helper to cast floating literals to correct precision based on the blitz array underlaying type
+    template<class arr_t>
+    constexpr auto fconst(const double v)
+    {
+      return static_cast<typename arr_t::T_numtype>(v);
+    }
     // overloads of abs/min/max/where that pick out the correct version based on ix_t
     template<class ix_t, class arg_t>
     forceinline_macro auto abs(const arg_t &a, typename std::enable_if<std::is_same<ix_t, int>::value>::type* = 0)

libmpdata++/formulae/mpdata/formulae_mpdata_1d.hpp

@@ -130,7 +130,7 @@ namespace libmpdataxx
       )
       {
         return return_helper<ix_t>(
-          - 1.0 / 24 *
+          - fconst<arr_1d_t>(1.0 / 24) *
           (
               4 * GC[0](i+h) * ndxx_psi<opts>(psi, i)
             + 2 * ndx_psi<opts>(psi, i) * ndx_GC0(GC[0], i)
@@ -170,9 +170,9 @@ namespace libmpdataxx
           // spatial terms
           + div_3rd_spatial<opts, sptl_intrp>(psi, GC, G, i)
           // mixed terms
-          + 0.5 * abs(GC[0](i+h)) * ndx_fdiv<opts>(psi, GC, G, i)
+          + fconst<arr_1d_t>(0.5) * abs(GC[0](i+h)) * ndx_fdiv<opts>(psi, GC, G, i)
           // temporal terms
-          + 1.0 / 24 *
+          + fconst<arr_1d_t>(1.0 / 24) *
           (
               - 8 * GC[0](i+h) *  nfdiv_fdiv<opts>(psi, GC, G, i)
               + div_3rd_temporal<opts, tmprl_extrp>(psi, ndtt_GC, i)

libmpdata++/formulae/mpdata/formulae_mpdata_2d.hpp

@@ -145,7 +145,7 @@ namespace libmpdataxx
       )
       {
         return return_helper<ix_t>(
-          - 1.0 / 24 *
+          - fconst<arr_2d_t>(1.0 / 24) *
           (
               4 * GC[dim](pi<dim>(i+h, j)) * ndxx_psi<opts, dim>(psi, i, j)
             + 2 * ndx_psi<opts, dim>(psi, i, j) * ndx_GC0<dim>(GC[dim], i, j)
@@ -193,9 +193,9 @@ namespace libmpdataxx
           // spatial terms
           + div_3rd_spatial<opts, dim, sptl_intrp>(psi_np1, GC, G, i, j)
           // mixed terms
-          + 0.5 * abs(GC[dim](pi<dim>(i+h, j))) * ndx_fdiv<opts, dim>(psi_np1, GC, G, i, j)
+          + fconst<arr_2d_t>(0.5) * abs(GC[dim](pi<dim>(i+h, j))) * ndx_fdiv<opts, dim>(psi_np1, GC, G, i, j)
           // temporal terms
-          + 1.0 / 24 *
+          + fconst<arr_2d_t>(1.0 / 24) *
           (
               - 8 * GC[dim](pi<dim>(i+h, j)) *  nfdiv_fdiv<opts, dim>(psi_np1, GC, G, i, j)
               + div_3rd_temporal<opts, dim, tmprl_extrp>(psi_np1, ndtt_GC, i, j)
@@ -226,9 +226,9 @@ namespace libmpdataxx
           // spatial terms
           + div_3rd_spatial<opts, dim, sptl_intrp>(psi_np1, GC, G, i, j)
           // mixed terms
-          - 0.5 * abs(GC[dim](pi<dim>(i+h, j))) * ndtx_psi<opts, dim>(psi_np1, psi_n, i, j)
+          - fconst<arr_2d_t>(0.5) * abs(GC[dim](pi<dim>(i+h, j))) * ndtx_psi<opts, dim>(psi_np1, psi_n, i, j)
           // temporal terms
-          + 1.0 / 24 *
+          + fconst<arr_2d_t>(1.0 / 24) *
           (
               + 8 * GC[dim](pi<dim>(i+h, j)) *  nfdiv_dt<opts, dim>(psi_np1, psi_n, GC, G, i, j)
               + 1 * ndtt_GC0<opts, dim>(psi_np1, ndtt_GC[dim], i, j)

libmpdata++/formulae/mpdata/formulae_mpdata_3d.hpp

@@ -154,7 +154,7 @@ namespace libmpdataxx
       )
       {
         return return_helper<ix_t>(
-          - 1.0 / 24 *
+          - fconst<arr_3d_t>(1.0 / 24) *
           (
               4 * GC[dim](pi<dim>(i+h, j, k)) * ndxx_psi<opts, dim>(psi, i, j, k)
             + 2 * ndx_psi<opts, dim>(psi, i, j, k) * ndx_GC0<dim>(GC[dim], i, j, k)
@@ -204,9 +204,9 @@ namespace libmpdataxx
           // spatial terms
           + div_3rd_spatial<opts, dim, sptl_intrp>(psi_np1, GC, G, i, j, k)
           // mixed terms
-          + 0.5 * abs(GC[dim](pi<dim>(i+h, j, k))) * ndx_fdiv<opts, dim>(psi_np1, GC, G, i, j, k)
+          + fconst<arr_3d_t>(0.5) * abs(GC[dim](pi<dim>(i+h, j, k))) * ndx_fdiv<opts, dim>(psi_np1, GC, G, i, j, k)
           // temporal terms
-          + 1.0 / 24 *
+          + fconst<arr_3d_t>(1.0 / 24) *
           (
               - 8 * GC[dim](pi<dim>(i+h, j, k)) *  nfdiv_fdiv<opts, dim>(psi_np1, GC, G, i, j, k)
               + div_3rd_temporal<opts, dim, tmprl_extrp>(psi_np1, ndtt_GC, i, j, k)
@@ -238,9 +238,9 @@ namespace libmpdataxx
           // spatial terms
           + div_3rd_spatial<opts, dim, sptl_intrp>(psi_np1, GC, G, i, j, k)
           // mixed terms
-          - 0.5 * abs(GC[dim](pi<dim>(i+h, j, k))) * ndtx_psi<opts, dim>(psi_np1, psi_n, i, j, k)
+          - fconst<arr_3d_t>(0.5) * abs(GC[dim](pi<dim>(i+h, j, k))) * ndtx_psi<opts, dim>(psi_np1, psi_n, i, j, k)
           // temporal terms
-          + 1.0 / 24 *
+          + fconst<arr_3d_t>(1.0 / 24) *
           (
               + 8 * GC[dim](pi<dim>(i+h, j, k)) *  nfdiv_dt<opts, dim>(psi_np1, psi_n, GC, G, i, j, k)
               + div_3rd_temporal<opts, dim, tmprl_extrp>(psi_np1, ndtt_GC, i, j, k)

libmpdata++/formulae/mpdata/formulae_mpdata_common.hpp

@@ -40,6 +40,10 @@ namespace libmpdataxx
       using idxperm::pi;
       using opts::opts_t;
       using std::abs;
+      
+      using blitz::pow2;
+      using blitz::pow3;
+      using blitz::pow4;
 
       const int n_tlev = 2;
 

libmpdata++/formulae/mpdata/formulae_mpdata_dfl_1d.hpp

@@ -38,7 +38,7 @@ namespace libmpdataxx
       )
       {
         return return_helper<ix_t>(
-          - 0.5 * GC(i+h) 
+          - fconst<arr_1d_t>(0.5) * GC(i+h) 
           / 
           (formulae::G<opts>(G, i+1) + formulae::G<opts>(G, i)) 
           * 
@@ -56,13 +56,13 @@ namespace libmpdataxx
       )
       {
         return return_helper<ix_t>(
-          - 0.5 * GC(i+h) 
+          - fconst<arr_1d_t>(0.5) * GC(i+h) 
           / 
           (formulae::G<opts>(G, i+1) + formulae::G<opts>(G, i)) 
           * 
           (GC((i+1)+h) - GC(i-h))
           *
-          0.5 *  (psi(i+1) + psi(i)) //to be compatible with iga formulation
+          fconst<arr_1d_t>(0.5) *  (psi(i+1) + psi(i)) //to be compatible with iga formulation
         );
       }
     } // namespace mpdata

libmpdata++/formulae/mpdata/formulae_mpdata_dfl_2d.hpp

@@ -42,7 +42,7 @@ namespace libmpdataxx
       )
       {
         return return_helper<ix_t>(
-          - 0.25 * GC[dim](pi<dim>(i+h, j)) 
+          - fconst<arr_2d_t>(0.25) * GC[dim](pi<dim>(i+h, j)) 
           /
           G_bar_x<opts, dim>(G, i, j) 
           * 
@@ -74,7 +74,7 @@ namespace libmpdataxx
       )
       {
         return return_helper<ix_t>(
-          - 0.25 * GC[dim](pi<dim>(i+h, j)) 
+          - fconst<arr_2d_t>(0.25) * GC[dim](pi<dim>(i+h, j)) 
           /
           G_bar_x<opts, dim>(G, i, j) 
           * 

libmpdata++/formulae/mpdata/formulae_mpdata_dfl_3d.hpp

@@ -44,7 +44,7 @@ namespace libmpdataxx
       )
       {
         return return_helper<ix_t>(
-          - 0.25 * GC[dim](pi<dim>(i+h, j, k)) 
+          - fconst<arr_3d_t>(0.25) * GC[dim](pi<dim>(i+h, j, k)) 
           /
           G_bar_x<opts, dim>(G, i, j, k)
           * 
@@ -84,7 +84,7 @@ namespace libmpdataxx
       )
       {
         return return_helper<ix_t>(
-          - 0.25 * GC[dim](pi<dim>(i+h, j, k)) 
+          - fconst<arr_3d_t>(0.25) * GC[dim](pi<dim>(i+h, j, k)) 
           /
           G_bar_x<opts, dim>(G, i, j, k)
           * 

libmpdata++/formulae/mpdata/formulae_mpdata_hot_1d.hpp

@@ -26,7 +26,7 @@ namespace libmpdataxx
         return return_helper<ix_t>(
           (
             3 * GC(i+h) * abs(GC(i+h)) / G_bar_x<opts>(G, i)
-            - 2 * pow(GC(i+h), 3) / pow(G_bar_x<opts>(G, i), 2)  
+            - 2 * pow3(GC(i+h)) / pow2(G_bar_x<opts>(G, i))
             - GC(i+h)
           ) / 6
         );