From 99ca119ee7adb608fb2b683ab9ea24dec64c24d2 Mon Sep 17 00:00:00 2001
From: Maciej Waruszewski <mwarusz@igf.fuw.edu.pl>
Date: Tue, 12 Mar 2019 17:04:50 +0100
Subject: [PATCH 1/4] initial work on avoiding float to double promotion

---
 libmpdata++/formulae/common.hpp               |   6 +
 .../formulae/mpdata/formulae_mpdata_1d.hpp    |   6 +-
 .../formulae/mpdata/formulae_mpdata_2d.hpp    |  10 +-
 .../formulae/mpdata/formulae_mpdata_3d.hpp    |  10 +-
 .../mpdata/formulae_mpdata_common.hpp         |   4 +
 .../mpdata/formulae_mpdata_dfl_1d.hpp         |   6 +-
 .../mpdata/formulae_mpdata_dfl_2d.hpp         |   4 +-
 .../mpdata/formulae_mpdata_dfl_3d.hpp         |   4 +-
 .../mpdata/formulae_mpdata_hot_1d.hpp         |   2 +-
 .../mpdata/formulae_mpdata_hot_2d.hpp         |  22 +--
 .../mpdata/formulae_mpdata_hot_3d.hpp         |   8 +-
 .../mpdata/formulae_mpdata_psi_3d.hpp         |   2 +-
 libmpdata++/formulae/nabla_formulae.hpp       |  16 +-
 libmpdata++/formulae/stress_formulae.hpp      | 140 +++++++++---------
 .../solvers/detail/boussinesq_expl.hpp        |  12 +-
 .../solvers/detail/boussinesq_impl.hpp        |  20 +--
 .../solvers/detail/boussinesq_sgs_common.hpp  |  12 +-
 .../solvers/detail/mpdata_rhs_vip_common.hpp  |  22 +--
 .../detail/mpdata_rhs_vip_prs_common.hpp      |   4 +-
 .../detail/mpdata_rhs_vip_prs_sgs_common.hpp  |   9 +-
 .../detail/mpdata_rhs_vip_prs_sgs_smg.hpp     |   9 +-
 libmpdata++/solvers/detail/solver_1d.hpp      |   2 +-
 libmpdata++/solvers/detail/solver_2d.hpp      |   2 +-
 libmpdata++/solvers/detail/solver_3d.hpp      |   2 +-
 libmpdata++/solvers/detail/solver_common.hpp  |   2 +-
 libmpdata++/solvers/mpdata_rhs_vip.hpp        |  14 +-
 26 files changed, 185 insertions(+), 165 deletions(-)

diff --git a/libmpdata++/formulae/common.hpp b/libmpdata++/formulae/common.hpp
index 66622dca..3b9846b6 100644
--- a/libmpdata++/formulae/common.hpp
+++ b/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)
diff --git a/libmpdata++/formulae/mpdata/formulae_mpdata_1d.hpp b/libmpdata++/formulae/mpdata/formulae_mpdata_1d.hpp
index 02956ae2..c7ef5073 100644
--- a/libmpdata++/formulae/mpdata/formulae_mpdata_1d.hpp
+++ b/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)
diff --git a/libmpdata++/formulae/mpdata/formulae_mpdata_2d.hpp b/libmpdata++/formulae/mpdata/formulae_mpdata_2d.hpp
index 8a8f2b12..1b7ce8b0 100644
--- a/libmpdata++/formulae/mpdata/formulae_mpdata_2d.hpp
+++ b/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)
diff --git a/libmpdata++/formulae/mpdata/formulae_mpdata_3d.hpp b/libmpdata++/formulae/mpdata/formulae_mpdata_3d.hpp
index d6b90213..2581b85b 100644
--- a/libmpdata++/formulae/mpdata/formulae_mpdata_3d.hpp
+++ b/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)
diff --git a/libmpdata++/formulae/mpdata/formulae_mpdata_common.hpp b/libmpdata++/formulae/mpdata/formulae_mpdata_common.hpp
index ecf0d583..09d14544 100644
--- a/libmpdata++/formulae/mpdata/formulae_mpdata_common.hpp
+++ b/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;
 
diff --git a/libmpdata++/formulae/mpdata/formulae_mpdata_dfl_1d.hpp b/libmpdata++/formulae/mpdata/formulae_mpdata_dfl_1d.hpp
index a942e5e4..e262ac98 100644
--- a/libmpdata++/formulae/mpdata/formulae_mpdata_dfl_1d.hpp
+++ b/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
diff --git a/libmpdata++/formulae/mpdata/formulae_mpdata_dfl_2d.hpp b/libmpdata++/formulae/mpdata/formulae_mpdata_dfl_2d.hpp
index 0b6dc7cc..68a3347c 100644
--- a/libmpdata++/formulae/mpdata/formulae_mpdata_dfl_2d.hpp
+++ b/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) 
           * 
diff --git a/libmpdata++/formulae/mpdata/formulae_mpdata_dfl_3d.hpp b/libmpdata++/formulae/mpdata/formulae_mpdata_dfl_3d.hpp
index acb0ad08..cba75024 100644
--- a/libmpdata++/formulae/mpdata/formulae_mpdata_dfl_3d.hpp
+++ b/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)
           * 
diff --git a/libmpdata++/formulae/mpdata/formulae_mpdata_hot_1d.hpp b/libmpdata++/formulae/mpdata/formulae_mpdata_hot_1d.hpp
index 3ab12b60..94b70891 100644
--- a/libmpdata++/formulae/mpdata/formulae_mpdata_hot_1d.hpp
+++ b/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
         );
diff --git a/libmpdata++/formulae/mpdata/formulae_mpdata_hot_2d.hpp b/libmpdata++/formulae/mpdata/formulae_mpdata_hot_2d.hpp
index c755c2f5..16d40160 100644
--- a/libmpdata++/formulae/mpdata/formulae_mpdata_hot_2d.hpp
+++ b/libmpdata++/formulae/mpdata/formulae_mpdata_hot_2d.hpp
@@ -29,7 +29,7 @@ namespace libmpdataxx
         return return_helper<ix_t>(
           (
             3 * GC(pi<dim>(i+h, j)) * abs(GC(pi<dim>(i+h, j))) / G_bar_x<opts, dim>(G, i, j)
-            - 2 * pow(GC(pi<dim>(i+h, j)), 3) / pow(G_bar_x<opts, dim>(G, i, j), 2)
+            - 2 * pow3(GC(pi<dim>(i+h, j))) / pow2(G_bar_x<opts, dim>(G, i, j))
             - GC(pi<dim>(i+h, j))
           ) / 6
         );
@@ -44,7 +44,7 @@ namespace libmpdataxx
       )
       {
         return return_helper<ix_t>(
-          (abs(GC[dim](pi<dim>(i+h, j))) - 2 * pow(GC[dim](pi<dim>(i+h, j)), 2) / G_bar_x<opts, dim>(G, i, j)) 
+          (abs(GC[dim](pi<dim>(i+h, j))) - 2 * pow2(GC[dim](pi<dim>(i+h, j))) / G_bar_x<opts, dim>(G, i, j)) 
            * GC1_bar_xy<dim>(GC[dim+1], i, j) / (2 * G_bar_x<opts, dim>(G, i, j))
         );
       }
@@ -90,9 +90,9 @@ namespace libmpdataxx
       {
         return return_helper<ix_t>(
           (
-            6 * abs(GC(pi<dim>(i+h, j))) * pow(GC(pi<dim>(i+h, j)), 2) / pow(G_bar_x<opts, dim>(G, i, j), 2)
-            - 3 * pow(GC(pi<dim>(i+h, j)), 2) / G_bar_x<opts, dim>(G, i, j)
-            - 3 * pow(GC(pi<dim>(i+h, j)), 4) / pow(G_bar_x<opts, dim>(G, i, j), 3)
+            6 * abs(GC(pi<dim>(i+h, j))) * pow2(GC(pi<dim>(i+h, j))) / pow2(G_bar_x<opts, dim>(G, i, j))
+            - 3 * pow2(GC(pi<dim>(i+h, j))) / G_bar_x<opts, dim>(G, i, j)
+            - 3 * pow4(GC(pi<dim>(i+h, j))) / pow3(G_bar_x<opts, dim>(G, i, j))
           ) / 2
         );
       }
@@ -107,9 +107,9 @@ namespace libmpdataxx
       {
         return return_helper<ix_t>(
           3 * 
-          (abs(GC[dim](pi<dim>(i+h, j))) - pow(GC[dim](pi<dim>(i+h, j)), 2) / G_bar_x<opts, dim>(G, i, j)) 
+          (abs(GC[dim](pi<dim>(i+h, j))) - pow2(GC[dim](pi<dim>(i+h, j))) / G_bar_x<opts, dim>(G, i, j)) 
            * GC[dim](pi<dim>(i+h, j)) * GC1_bar_xy<dim>(GC[dim-1], i, j)
-           / pow(G_bar_x<opts, dim>(G, i, j), 2)
+           / pow2(G_bar_x<opts, dim>(G, i, j))
         );
       }
       
@@ -124,11 +124,11 @@ namespace libmpdataxx
         return return_helper<ix_t>(
           3 *
           (
-              3 * abs(GC1_bar_xy<dim>(GC[dim-1], i, j)) * pow(GC[dim](pi<dim>(i+h, j)), 2) / G_bar_x<opts, dim>(G, i, j)
-            + 3 * pow(GC1_bar_xy<dim>(GC[dim-1], i, j), 2) * abs(GC[dim](pi<dim>(i+h, j))) / G_bar_x<opts, dim>(G, i, j)
+              3 * abs(GC1_bar_xy<dim>(GC[dim-1], i, j)) * pow2(GC[dim](pi<dim>(i+h, j))) / G_bar_x<opts, dim>(G, i, j)
+            + 3 * pow2(GC1_bar_xy<dim>(GC[dim-1], i, j)) * abs(GC[dim](pi<dim>(i+h, j))) / G_bar_x<opts, dim>(G, i, j)
             - 2 * abs(GC1_bar_xy<dim>(GC[dim-1], i, j)) * abs(GC[dim](pi<dim>(i+h, j)))
-            - 9 * pow(GC1_bar_xy<dim>(GC[dim-1], i, j), 2) * pow(GC[dim](pi<dim>(i+h, j)), 2)
-                / pow(G_bar_x<opts, dim>(G, i, j), 2)
+            - 9 * pow2(GC1_bar_xy<dim>(GC[dim-1], i, j)) * pow2(GC[dim](pi<dim>(i+h, j)))
+                / pow2(G_bar_x<opts, dim>(G, i, j))
           ) / G_bar_x<opts, dim>(G, i, j) / 4
         );
       }
diff --git a/libmpdata++/formulae/mpdata/formulae_mpdata_hot_3d.hpp b/libmpdata++/formulae/mpdata/formulae_mpdata_hot_3d.hpp
index 69a73854..fc2afd3f 100644
--- a/libmpdata++/formulae/mpdata/formulae_mpdata_hot_3d.hpp
+++ b/libmpdata++/formulae/mpdata/formulae_mpdata_hot_3d.hpp
@@ -30,7 +30,7 @@ namespace libmpdataxx
         return return_helper<ix_t>(
           (
             3 * GC(pi<dim>(i+h, j, k)) * abs(GC(pi<dim>(i+h, j, k))) / G_bar_x<opts, dim>(G, i, j, k)
-          - 2 * pow(GC(pi<dim>(i+h, j, k)), 3) / pow(G_bar_x<opts, dim>(G, i, j, k), 2)
+          - 2 * pow3(GC(pi<dim>(i+h, j, k))) / pow2(G_bar_x<opts, dim>(G, i, j, k))
           - GC(pi<dim>(i+h, j, k))
           ) / 6
         );
@@ -48,7 +48,7 @@ namespace libmpdataxx
         return return_helper<ix_t>(
           (
             abs(GC[dim](pi<dim>(i+h, j, k)))
-          - 2 * pow(GC[dim](pi<dim>(i+h, j, k)), 2) / G_bar_x<opts, dim>(G, i, j, k)
+          - 2 * pow2(GC[dim](pi<dim>(i+h, j, k))) / G_bar_x<opts, dim>(G, i, j, k)
           ) * GC1_bar_xy<dim>(GC[dim+1], i, j, k) / (2 * G_bar_x<opts, dim>(G, i, j, k))
         );
       }
@@ -65,7 +65,7 @@ namespace libmpdataxx
         return return_helper<ix_t>(
           (
             abs(GC[dim](pi<dim>(i+h, j, k)))
-          - 2 * pow(GC[dim](pi<dim>(i+h, j, k)), 2) / G_bar_x<opts, dim>(G, i, j, k)
+          - 2 * pow2(GC[dim](pi<dim>(i+h, j, k))) / G_bar_x<opts, dim>(G, i, j, k)
           ) * GC2_bar_xz<dim>(GC[dim-1], i, j, k) / (2 * G_bar_x<opts, dim>(G, i, j, k))
         );
       }
@@ -81,7 +81,7 @@ namespace libmpdataxx
       {
         return return_helper<ix_t>(
           -2 * GC[dim](pi<dim>(i+h, j, k)) * GC1_bar_xy<dim>(GC[dim+1], i, j, k) * GC2_bar_xz<dim>(GC[dim-1], i, j, k) / 3
-             / pow(G_bar_x<opts, dim>(G, i, j, k), 2)
+             / pow2(G_bar_x<opts, dim>(G, i, j, k))
         );
       }
 
diff --git a/libmpdata++/formulae/mpdata/formulae_mpdata_psi_3d.hpp b/libmpdata++/formulae/mpdata/formulae_mpdata_psi_3d.hpp
index cc061009..2f0c30c6 100644
--- a/libmpdata++/formulae/mpdata/formulae_mpdata_psi_3d.hpp
+++ b/libmpdata++/formulae/mpdata/formulae_mpdata_psi_3d.hpp
@@ -933,7 +933,7 @@ namespace libmpdataxx
       {
         static_assert(!opts::isset(opts, opts::abs), "abs & iga options are mutually exclusive");
         return return_helper<ix_t>(
-          0.5 * (psi_np1(pi<d>(i+1, j, k)) - psi_n(pi<d>(i+1, j, k)) + psi_np1(pi<d>(i, j, k)) - psi_n(pi<d>(i, j, k)))
+          fconst<arr_3d_t>(0.5) * (psi_np1(pi<d>(i+1, j, k)) - psi_n(pi<d>(i+1, j, k)) + psi_np1(pi<d>(i, j, k)) - psi_n(pi<d>(i, j, k)))
         );
       }
     } // namespace mpdata
diff --git a/libmpdata++/formulae/nabla_formulae.hpp b/libmpdata++/formulae/nabla_formulae.hpp
index 07a16f1d..6ff594b1 100644
--- a/libmpdata++/formulae/nabla_formulae.hpp
+++ b/libmpdata++/formulae/nabla_formulae.hpp
@@ -30,7 +30,7 @@ namespace libmpdataxx
   	  (
 	    x(i+1) - 
 	    x(i-1)
-	  ) / dx / 2.
+	  ) / dx / 2
         );
       }
 
@@ -47,7 +47,7 @@ namespace libmpdataxx
   	  (
   	    x(pi<d>(i+1, j)) - 
   	    x(pi<d>(i-1, j))
-  	  ) / dx / 2.
+  	  ) / dx / 2
         );
       }
       
@@ -65,7 +65,7 @@ namespace libmpdataxx
   	  (
 	    x(pi<d>(i+1, j, k)) - 
 	    x(pi<d>(i-1, j, k))
-	  ) / dx / 2.
+	  ) / dx / 2
         );
       }
       
@@ -174,9 +174,9 @@ namespace libmpdataxx
       ) 
       {
         return blitz::safeToReturn(
-	  (v[0](ijk[0]+1, ijk[1]) - v[0](ijk[0]-1, ijk[1])) / dijk[0] / 2.
+	  (v[0](ijk[0]+1, ijk[1]) - v[0](ijk[0]-1, ijk[1])) / dijk[0] / 2
 	  +
-	  (v[1](ijk[0], ijk[1]+1) - v[1](ijk[0], ijk[1]-1)) / dijk[1] / 2.
+	  (v[1](ijk[0], ijk[1]+1) - v[1](ijk[0], ijk[1]-1)) / dijk[1] / 2
         );
       }
       
@@ -190,11 +190,11 @@ namespace libmpdataxx
       ) 
       {
         return blitz::safeToReturn(
-	  (v[0](ijk[0]+1, ijk[1], ijk[2]) - v[0](ijk[0]-1, ijk[1], ijk[2])) / dijk[0] / 2.
+	  (v[0](ijk[0]+1, ijk[1], ijk[2]) - v[0](ijk[0]-1, ijk[1], ijk[2])) / dijk[0] / 2
 	  +
-	  (v[1](ijk[0], ijk[1]+1, ijk[2]) - v[1](ijk[0], ijk[1]-1, ijk[2])) / dijk[1] / 2.
+	  (v[1](ijk[0], ijk[1]+1, ijk[2]) - v[1](ijk[0], ijk[1]-1, ijk[2])) / dijk[1] / 2
 	  +
-	  (v[2](ijk[0], ijk[1], ijk[2]+1) - v[2](ijk[0], ijk[1], ijk[2]-1)) / dijk[2] / 2.
+	  (v[2](ijk[0], ijk[1], ijk[2]+1) - v[2](ijk[0], ijk[1], ijk[2]-1)) / dijk[2] / 2
         );
       }
     } // namespace nabla_op
diff --git a/libmpdata++/formulae/stress_formulae.hpp b/libmpdata++/formulae/stress_formulae.hpp
index 71d6efdf..df3e117c 100644
--- a/libmpdata++/formulae/stress_formulae.hpp
+++ b/libmpdata++/formulae/stress_formulae.hpp
@@ -264,20 +264,20 @@ namespace libmpdataxx
         =
         2 * (
           (v[0](ijkm[0] + 1, ijk[1]) - v[0](ijkm[0], ijk[1])) / dijk[0]
-          - 0.25 / 2 * ( div_v(ijkm[0], ijk[1] + h) + div_v(ijkm[0], ijk[1] - h) 
-                       + div_v(ijkm[0] + 1, ijk[1] + h) + div_v(ijkm[0] + 1, ijk[1] - h))
+          - fconst<arr_t>(0.25 / 2) * ( div_v(ijkm[0], ijk[1] + h) + div_v(ijkm[0], ijk[1] - h) 
+                                      + div_v(ijkm[0] + 1, ijk[1] + h) + div_v(ijkm[0] + 1, ijk[1] - h))
         );
         
         tau[1](ijk[0], ijkm[1] + h)
         =
         2 * (
           (v[1](ijk[0], ijkm[1] + 1) - v[1](ijk[0], ijkm[1])) / dijk[1]
-          - 1.0 / 2 * div_v(ijk[0], ijkm[1] + h)
+          - fconst<arr_t>(1.0 / 2) * div_v(ijk[0], ijkm[1] + h)
         );
 
         tau[2](ijkm[0] + h, ijkm[1] + h)
         =
-        0.5 *
+        fconst<arr_t>(0.5) *
         (
           ( v[0](ijkm[0], ijkm[1] + 1) - v[0](ijkm[0], ijkm[1]) 
           + v[0](ijkm[0] + 1, ijkm[1] + 1) - v[0](ijkm[0] + 1, ijkm[1])
@@ -303,28 +303,28 @@ namespace libmpdataxx
         =
         2 * (
           (v[0](ijkm[0] + 1, ijk[1], ijk[2]) - v[0](ijkm[0], ijk[1], ijk[2])) / dijk[0]
-          - 0.25 / 3 * ( div_v(ijkm[0], ijk[1], ijk[2] + h) + div_v(ijkm[0], ijk[1], ijk[2] - h)
-                       + div_v(ijkm[0] + 1, ijk[1], ijk[2] + h) + div_v(ijkm[0] + 1, ijk[1], ijk[2] - h))
+          - fconst<arr_t>(0.25 / 3) * ( div_v(ijkm[0], ijk[1], ijk[2] + h) + div_v(ijkm[0], ijk[1], ijk[2] - h)
+                                      + div_v(ijkm[0] + 1, ijk[1], ijk[2] + h) + div_v(ijkm[0] + 1, ijk[1], ijk[2] - h))
         );
         
         tau[1](ijk[0], ijkm[1] + h, ijk[2])
         =
         2 * (
           (v[1](ijk[0], ijkm[1] + 1, ijk[2]) - v[1](ijk[0], ijkm[1], ijk[2])) / dijk[1]
-          - 0.25 / 3 * ( div_v(ijk[0], ijkm[1], ijk[2] + h) + div_v(ijk[0], ijkm[1], ijk[2] - h)
-                       + div_v(ijk[0], ijkm[1] + 1, ijk[2] + h) + div_v(ijk[0], ijkm[1] + 1, ijk[2] - h))
+          - fconst<arr_t>(0.25 / 3) * ( div_v(ijk[0], ijkm[1], ijk[2] + h) + div_v(ijk[0], ijkm[1], ijk[2] - h)
+                                    + div_v(ijk[0], ijkm[1] + 1, ijk[2] + h) + div_v(ijk[0], ijkm[1] + 1, ijk[2] - h))
         );
         
         tau[2](ijk[0], ijk[1], ijkm[2] + h)
         =
         2 * (
           (v[2](ijk[0], ijk[1], ijkm[2] + 1) - v[2](ijk[0], ijk[1], ijkm[2])) / dijk[2]
-          - 1.0 / 3  * div_v(ijk[0], ijk[1], ijkm[2] + h) 
+          - fconst<arr_t>(1.0 / 3)  * div_v(ijk[0], ijk[1], ijkm[2] + h) 
         );
 
         tau[3](ijkm[0] + h, ijkm[1] + h, ijk[2])
         =
-        0.5 *
+        fconst<arr_t>(0.5) *
         (
           ( v[0](ijkm[0], ijkm[1] + 1, ijk[2]) - v[0](ijkm[0], ijkm[1], ijk[2]) 
           + v[0](ijkm[0] + 1, ijkm[1] + 1, ijk[2]) - v[0](ijkm[0] + 1, ijkm[1], ijk[2])
@@ -337,7 +337,7 @@ namespace libmpdataxx
         
         tau[4](ijkm[0] + h, ijk[1], ijkm[2] + h)
         =
-        0.5 *
+        fconst<arr_t>(0.5) *
         (
           ( v[0](ijkm[0], ijk[1], ijkm[2] + 1) - v[0](ijkm[0], ijk[1], ijkm[2])
           + v[0](ijkm[0] + 1, ijk[1], ijkm[2] + 1) - v[0](ijkm[0] + 1, ijk[1], ijkm[2])
@@ -350,7 +350,7 @@ namespace libmpdataxx
         
         tau[5](ijk[0], ijkm[1] + h, ijkm[2] + h)
         =
-        0.5 *
+        fconst<arr_t>(0.5) *
         (
           ( v[1](ijk[0], ijkm[1], ijkm[2] + 1) - v[1](ijk[0], ijkm[1], ijkm[2])
           + v[1](ijk[0], ijkm[1] + 1, ijkm[2] + 1) - v[1](ijk[0], ijkm[1] + 1, ijkm[2])
@@ -462,12 +462,12 @@ namespace libmpdataxx
                                       typename std::enable_if<nd == 2>::type* = 0)
       {
         av[0](ijk[0] + h, ijk[1]) *= coeff * 
-                           0.5 * (k_m(ijk[0] + 1, ijk[1]) + k_m(ijk[0], ijk[1])) *
-                           0.5 * (G<opts, 0>(rho, ijk[0] + 1, ijk[1]) + G<opts, 0>(rho, ijk[0], ijk[1]));
+                           real_t(0.5) * (k_m(ijk[0] + 1, ijk[1]) + k_m(ijk[0], ijk[1])) *
+                           real_t(0.5) * (G<opts, 0>(rho, ijk[0] + 1, ijk[1]) + G<opts, 0>(rho, ijk[0], ijk[1]));
         
         av[1](ijk[0], ijk[1] + h) *= coeff *
-                           0.5 * (k_m(ijk[0], ijk[1] + 1) + k_m(ijk[0], ijk[1])) * 
-                           0.5 * (G<opts, 0>(rho, ijk[0], ijk[1] + 1) + G<opts, 0>(rho, ijk[0], ijk[1]));
+                           real_t(0.5) * (k_m(ijk[0], ijk[1] + 1) + k_m(ijk[0], ijk[1])) * 
+                           real_t(0.5) * (G<opts, 0>(rho, ijk[0], ijk[1] + 1) + G<opts, 0>(rho, ijk[0], ijk[1]));
       }
 
       // 3D version
@@ -480,16 +480,16 @@ namespace libmpdataxx
                                       typename std::enable_if<nd == 3>::type* = 0)
       {
         av[0](ijk[0] + h, ijk[1], ijk[2]) *= coeff *
-                           0.5 * (k_m(ijk[0] + 1, ijk[1], ijk[2]) + k_m(ijk[0], ijk[1], ijk[2])) *
-                           0.5 * (G<opts, 0>(rho, ijk[0] + 1, ijk[1], ijk[2]) + G<opts, 0>(rho,ijk[0], ijk[1], ijk[2]));
+                           real_t(0.5) * (k_m(ijk[0] + 1, ijk[1], ijk[2]) + k_m(ijk[0], ijk[1], ijk[2])) *
+                           real_t(0.5) * (G<opts, 0>(rho, ijk[0] + 1, ijk[1], ijk[2]) + G<opts, 0>(rho,ijk[0], ijk[1], ijk[2]));
         
         av[1](ijk[0], ijk[1] + h, ijk[2]) *= coeff *
-                           0.5 * (k_m(ijk[0], ijk[1] + 1, ijk[2]) + k_m(ijk[0], ijk[1], ijk[2])) *
-                           0.5 * (G<opts, 0>(rho, ijk[0], ijk[1] + 1, ijk[2]) + G<opts, 0>(rho, ijk[0], ijk[1], ijk[2]));
+                           real_t(0.5) * (k_m(ijk[0], ijk[1] + 1, ijk[2]) + k_m(ijk[0], ijk[1], ijk[2])) *
+                           real_t(0.5) * (G<opts, 0>(rho, ijk[0], ijk[1] + 1, ijk[2]) + G<opts, 0>(rho, ijk[0], ijk[1], ijk[2]));
         
         av[2](ijk[0], ijk[1], ijk[2] + h) *= coeff *
-                           0.5 * (k_m(ijk[0], ijk[1], ijk[2] + 1) + k_m(ijk[0], ijk[1], ijk[2])) *
-                           0.5 * (G<opts, 0>(rho, ijk[0], ijk[1], ijk[2] + 1) + G<opts, 0>(rho, ijk[0], ijk[1], ijk[2]));
+                           real_t(0.5) * (k_m(ijk[0], ijk[1], ijk[2] + 1) + k_m(ijk[0], ijk[1], ijk[2])) *
+                           real_t(0.5) * (G<opts, 0>(rho, ijk[0], ijk[1], ijk[2] + 1) + G<opts, 0>(rho, ijk[0], ijk[1], ijk[2]));
       }
       
       // multiplication of compact tensor components by constant molecular viscosity
@@ -528,16 +528,16 @@ namespace libmpdataxx
       {
         multiply_vctr_cmpct<nd, opts>(av, coeff, k_m, rho, ijk);
         av[2](ijk[0] + h, ijk[1] + h) *= coeff * 
-                                         0.25 * ( k_m(ijk[0] + 1, ijk[1]    )
-                                                + k_m(ijk[0]    , ijk[1]    )
-                                                + k_m(ijk[0] + 1, ijk[1] + 1)
-                                                + k_m(ijk[0]    , ijk[1] + 1)
-                                                ) *
-                                         0.25 * ( G<opts, 0>(rho, ijk[0] + 1, ijk[1]    )
-                                                + G<opts, 0>(rho, ijk[0]    , ijk[1]    )
-                                                + G<opts, 0>(rho, ijk[0] + 1, ijk[1] + 1)
-                                                + G<opts, 0>(rho, ijk[0]    , ijk[1] + 1)
-                                                );
+                                         real_t(0.25) * ( k_m(ijk[0] + 1, ijk[1]    )
+                                                        + k_m(ijk[0]    , ijk[1]    )
+                                                        + k_m(ijk[0] + 1, ijk[1] + 1)
+                                                        + k_m(ijk[0]    , ijk[1] + 1)
+                                                        ) *
+                                         real_t(0.25) * ( G<opts, 0>(rho, ijk[0] + 1, ijk[1]    )
+                                                        + G<opts, 0>(rho, ijk[0]    , ijk[1]    )
+                                                        + G<opts, 0>(rho, ijk[0] + 1, ijk[1] + 1)
+                                                        + G<opts, 0>(rho, ijk[0]    , ijk[1] + 1)
+                                                        );
       }
 
       // 3D version
@@ -551,40 +551,40 @@ namespace libmpdataxx
       {
         multiply_vctr_cmpct<nd, opts>(av, coeff, k_m, rho, ijk);
         av[3](ijk[0] + h, ijk[1] + h, ijk[2]) *= coeff *
-                                                 0.25 * ( k_m(ijk[0] + 1, ijk[1]    , ijk[2])
-                                                        + k_m(ijk[0]    , ijk[1]    , ijk[2])
-                                                        + k_m(ijk[0] + 1, ijk[1] + 1, ijk[2])
-                                                        + k_m(ijk[0]    , ijk[1] + 1, ijk[2])
-                                                        ) *
-                                                 0.25 * ( G<opts, 0>(rho, ijk[0] + 1, ijk[1]    , ijk[2])
-                                                        + G<opts, 0>(rho, ijk[0]    , ijk[1]    , ijk[2])
-                                                        + G<opts, 0>(rho, ijk[0] + 1, ijk[1] + 1, ijk[2])
-                                                        + G<opts, 0>(rho, ijk[0]    , ijk[1] + 1, ijk[2])
-                                                        );
+                                                 real_t(0.25) * ( k_m(ijk[0] + 1, ijk[1]    , ijk[2])
+                                                                + k_m(ijk[0]    , ijk[1]    , ijk[2])
+                                                                + k_m(ijk[0] + 1, ijk[1] + 1, ijk[2])
+                                                                + k_m(ijk[0]    , ijk[1] + 1, ijk[2])
+                                                                ) *
+                                                 real_t(0.25) * ( G<opts, 0>(rho, ijk[0] + 1, ijk[1]    , ijk[2])
+                                                                + G<opts, 0>(rho, ijk[0]    , ijk[1]    , ijk[2])
+                                                                + G<opts, 0>(rho, ijk[0] + 1, ijk[1] + 1, ijk[2])
+                                                                + G<opts, 0>(rho, ijk[0]    , ijk[1] + 1, ijk[2])
+                                                                );
         
         av[4](ijk[0] + h, ijk[1], ijk[2] + h) *= coeff *
-                                                 0.25 * ( k_m(ijk[0] + 1, ijk[1], ijk[2]    )
-                                                        + k_m(ijk[0]    , ijk[1], ijk[2]    )
-                                                        + k_m(ijk[0] + 1, ijk[1], ijk[2] + 1)
-                                                        + k_m(ijk[0]    , ijk[1], ijk[2] + 1)
-                                                        ) *
-                                                 0.25 * ( G<opts, 0>(rho, ijk[0] + 1, ijk[1], ijk[2]    )
-                                                        + G<opts, 0>(rho, ijk[0]    , ijk[1], ijk[2]    )
-                                                        + G<opts, 0>(rho, ijk[0] + 1, ijk[1], ijk[2] + 1)
-                                                        + G<opts, 0>(rho, ijk[0]    , ijk[1], ijk[2] + 1)
-                                                        );
+                                                 real_t(0.25) * ( k_m(ijk[0] + 1, ijk[1], ijk[2]    )
+                                                                + k_m(ijk[0]    , ijk[1], ijk[2]    )
+                                                                + k_m(ijk[0] + 1, ijk[1], ijk[2] + 1)
+                                                                + k_m(ijk[0]    , ijk[1], ijk[2] + 1)
+                                                                ) *
+                                                 real_t(0.25) * ( G<opts, 0>(rho, ijk[0] + 1, ijk[1], ijk[2]    )
+                                                                + G<opts, 0>(rho, ijk[0]    , ijk[1], ijk[2]    )
+                                                                + G<opts, 0>(rho, ijk[0] + 1, ijk[1], ijk[2] + 1)
+                                                                + G<opts, 0>(rho, ijk[0]    , ijk[1], ijk[2] + 1)
+                                                                );
 
         av[5](ijk[0], ijk[1] + h, ijk[2] + h) *= coeff *
-                                                 0.25 * ( k_m(ijk[0], ijk[1]    , ijk[2] + 1)
-                                                        + k_m(ijk[0], ijk[1]    , ijk[2]    )
-                                                        + k_m(ijk[0], ijk[1] + 1, ijk[2] + 1)
-                                                        + k_m(ijk[0], ijk[1] + 1, ijk[2]    )
-                                                        ) *
-                                                 0.25 * ( G<opts, 0>(rho, ijk[0], ijk[1]    , ijk[2] + 1)
-                                                        + G<opts, 0>(rho, ijk[0], ijk[1]    , ijk[2]    )
-                                                        + G<opts, 0>(rho, ijk[0], ijk[1] + 1, ijk[2] + 1)
-                                                        + G<opts, 0>(rho, ijk[0], ijk[1] + 1, ijk[2]    )
-                                                        );
+                                                 real_t(0.25) * ( k_m(ijk[0], ijk[1]    , ijk[2] + 1)
+                                                                + k_m(ijk[0], ijk[1]    , ijk[2]    )
+                                                                + k_m(ijk[0], ijk[1] + 1, ijk[2] + 1)
+                                                                + k_m(ijk[0], ijk[1] + 1, ijk[2]    )
+                                                                ) *
+                                                 real_t(0.25) * ( G<opts, 0>(rho, ijk[0], ijk[1]    , ijk[2] + 1)
+                                                                + G<opts, 0>(rho, ijk[0], ijk[1]    , ijk[2]    )
+                                                                + G<opts, 0>(rho, ijk[0], ijk[1] + 1, ijk[2] + 1)
+                                                                + G<opts, 0>(rho, ijk[0], ijk[1] + 1, ijk[2]    )
+                                                                );
       }
 
       // flux divergence
@@ -642,7 +642,7 @@ namespace libmpdataxx
         coeff *
         (
           (tau[0](ijk[0] + h, ijk[1]) - tau[0](ijk[0] - h, ijk[1])) / dijk[0]
-          + 0.5 * 
+          + real_t(0.5) * 
           ( tau[2](ijk[0] + h, ijk[1] + h) - tau[2](ijk[0] + h, ijk[1] - h)
           + tau[2](ijk[0] - h, ijk[1] + h) - tau[2](ijk[0] - h, ijk[1] - h)
           ) / dijk[1]
@@ -652,7 +652,7 @@ namespace libmpdataxx
         +=
         coeff *
         (
-          0.5 * 
+          real_t(0.5) * 
           ( tau[2](ijk[0] + h, ijk[1] + h) - tau[2](ijk[0] - h, ijk[1] + h)
           + tau[2](ijk[0] + h, ijk[1] - h) - tau[2](ijk[0] - h, ijk[1] - h)
           ) / dijk[0]
@@ -676,11 +676,11 @@ namespace libmpdataxx
         coeff *
         (
           (tau[0](ijk[0] + h, ijk[1], ijk[2]) - tau[0](ijk[0] - h, ijk[1], ijk[2])) / dijk[0]
-          + 0.5 * 
+          + real_t(0.5) * 
           ( tau[3](ijk[0] + h, ijk[1] + h, ijk[2]) - tau[3](ijk[0] + h, ijk[1] - h, ijk[2])
           + tau[3](ijk[0] - h, ijk[1] + h, ijk[2]) - tau[3](ijk[0] - h, ijk[1] - h, ijk[2])
           ) / dijk[1]
-          + 0.5 * 
+          + real_t(0.5) * 
           ( tau[4](ijk[0] + h, ijk[1], ijk[2] + h) - tau[4](ijk[0] + h, ijk[1], ijk[2] - h)
           + tau[4](ijk[0] - h, ijk[1], ijk[2] + h) - tau[4](ijk[0] - h, ijk[1], ijk[2] - h)
           ) / dijk[2]
@@ -690,13 +690,13 @@ namespace libmpdataxx
         +=
         coeff *
         (
-          0.5 * 
+          real_t(0.5) * 
           ( tau[3](ijk[0] + h, ijk[1] + h, ijk[2]) - tau[3](ijk[0] - h, ijk[1] + h, ijk[2])
           + tau[3](ijk[0] + h, ijk[1] - h, ijk[2]) - tau[3](ijk[0] - h, ijk[1] - h, ijk[2])
           ) / dijk[0]
           +
           (tau[1](ijk[0], ijk[1] + h, ijk[2]) - tau[1](ijk[0], ijk[1] - h, ijk[2])) / dijk[1]
-          + 0.5 * 
+          + real_t(0.5) * 
           ( tau[5](ijk[0], ijk[1] + h, ijk[2] + h) - tau[5](ijk[0], ijk[1] + h, ijk[2] - h)
           + tau[5](ijk[0], ijk[1] - h, ijk[2] + h) - tau[5](ijk[0], ijk[1] - h, ijk[2] - h)
           ) / dijk[2]
@@ -706,11 +706,11 @@ namespace libmpdataxx
         +=
         coeff *
         (
-          0.5 * 
+          real_t(0.5) * 
           ( tau[4](ijk[0] + h, ijk[1], ijk[2] + h) - tau[4](ijk[0] - h, ijk[1], ijk[2] + h)
           + tau[4](ijk[0] + h, ijk[1], ijk[2] - h) - tau[4](ijk[0] - h, ijk[1], ijk[2] - h)
           ) / dijk[0]
-          + 0.5 * 
+          + real_t(0.5) * 
           ( tau[5](ijk[0], ijk[1] + h, ijk[2] + h) - tau[5](ijk[0], ijk[1] - h, ijk[2] + h)
           + tau[5](ijk[0], ijk[1] + h, ijk[2] - h) - tau[5](ijk[0], ijk[1] - h, ijk[2] - h)
           ) / dijk[1]
diff --git a/libmpdata++/solvers/detail/boussinesq_expl.hpp b/libmpdata++/solvers/detail/boussinesq_expl.hpp
index a9d7a38f..48b68966 100644
--- a/libmpdata++/solvers/detail/boussinesq_expl.hpp
+++ b/libmpdata++/solvers/detail/boussinesq_expl.hpp
@@ -36,7 +36,7 @@ namespace libmpdataxx
         {
           const auto &i(this->i), &j(this->j);
           this->xchng_sclr(tmp1, this->ijk);
-          tmp2(i, j) = 0.25 * (tmp1(i, j + 1) + 2 * tmp1(i, j) + tmp1(i, j - 1));
+          tmp2(i, j) = real_t(0.25) * (tmp1(i, j + 1) + 2 * tmp1(i, j) + tmp1(i, j - 1));
         }
 
         template <int nd = ct_params_t::n_dims> 
@@ -44,7 +44,7 @@ namespace libmpdataxx
         {
           const auto &i(this->i), &j(this->j), &k(this->k);
           this->xchng_sclr(tmp1, this->ijk);
-          tmp2(i, j, k) = 0.25 * (tmp1(i, j, k + 1) + 2 * tmp1(i, j, k) + tmp1(i, j, k - 1));
+          tmp2(i, j, k) = real_t(0.25) * (tmp1(i, j, k + 1) + 2 * tmp1(i, j, k) + tmp1(i, j, k - 1));
         }
         
         // helpers for buoyancy forces
@@ -62,8 +62,8 @@ namespace libmpdataxx
           return return_helper<rng_t>(
             this->g * (
                 (  this->state(ix::tht)(ijk)
-                 + 0.5 * this->dt * this->hflux_frc(ijk) + 0.5 * this->dt * this->tht_abs(ijk) * this->tht_e(ijk))
-                / (1 + 0.5 * this->dt * this->tht_abs(ijk))
+                 + real_t(0.5) * this->dt * this->hflux_frc(ijk) + real_t(0.5) * this->dt * this->tht_abs(ijk) * this->tht_e(ijk))
+                / (1 + real_t(0.5) * this->dt * this->tht_abs(ijk))
                 - this->tht_e(ijk)
               ) / this->Tht_ref
           );
@@ -110,8 +110,8 @@ namespace libmpdataxx
             case (1):
             {
               rhs.at(ix::tht)(ijk) += this->hflux_frc(ijk) - this->tht_abs(ijk) * (
-                (tht(ijk) + 0.5 * this->dt * this->hflux_frc(ijk) + 0.5 * this->dt * this->tht_abs(ijk) * this->tht_e(ijk))
-                / (1 + 0.5 * this->dt * this->tht_abs(ijk))
+                (tht(ijk) + real_t(0.5) * this->dt * this->hflux_frc(ijk) + real_t(0.5) * this->dt * this->tht_abs(ijk) * this->tht_e(ijk))
+                / (1 + real_t(0.5) * this->dt * this->tht_abs(ijk))
                 - this->tht_e(ijk)
               );
 
diff --git a/libmpdata++/solvers/detail/boussinesq_impl.hpp b/libmpdata++/solvers/detail/boussinesq_impl.hpp
index eb5d0ad5..be1ccb5b 100644
--- a/libmpdata++/solvers/detail/boussinesq_impl.hpp
+++ b/libmpdata++/solvers/detail/boussinesq_impl.hpp
@@ -61,18 +61,18 @@ namespace libmpdataxx
           {
             for (int d = 0; d < ct_params_t::n_dims - 1; ++d)
             {
-              v[d](this->ijk) /= (1 + 0.5 * this->dt * (*this->mem->vab_coeff)(this->ijk));
+              v[d](this->ijk) /= (1 + real_t(0.5) * this->dt * (*this->mem->vab_coeff)(this->ijk));
             }
             v[ct_params_t::n_dims - 1](this->ijk) /=
-            (1 + 0.5 * this->dt * (*this->mem->vab_coeff)(this->ijk)
-               + 0.25 * this->dt * this->dt * this->g / this->Tht_ref * this->dtht_e(this->ijk)
-                 / (1 + 0.5 * this->dt * this->tht_abs(this->ijk)));
+            (1 + real_t(0.5) * this->dt * (*this->mem->vab_coeff)(this->ijk)
+               + real_t(0.25) * this->dt * this->dt * this->g / this->Tht_ref * this->dtht_e(this->ijk)
+                 / (1 + real_t(0.5) * this->dt * this->tht_abs(this->ijk)));
           }
           else
           {
             v[ct_params_t::n_dims - 1](this->ijk) /=
-            (1 + 0.25 * this->dt * this->dt * this->g / this->Tht_ref * this->dtht_e(this->ijk)
-                 / (1 + 0.5 * this->dt * this->tht_abs(this->ijk)));
+            (1 + real_t(0.25) * this->dt * this->dt * this->g / this->Tht_ref * this->dtht_e(this->ijk)
+                 / (1 + real_t(0.5) * this->dt * this->tht_abs(this->ijk)));
           }
         }
         
@@ -104,8 +104,8 @@ namespace libmpdataxx
             {
               rhs.at(ix::tht)(ijk) += this->hflux_frc(ijk);
 
-              rhs.at(ix_w)(ijk) += this->g * (tht(ijk) + 0.5 * this->dt * rhs.at(ix::tht)(ijk))
-                                   / (this->Tht_ref * (1 + 0.5 * this->dt * this->tht_abs(ijk)));
+              rhs.at(ix_w)(ijk) += this->g * (tht(ijk) + real_t(0.5) * this->dt * rhs.at(ix::tht)(ijk))
+                                   / (this->Tht_ref * (1 + real_t(0.5) * this->dt * this->tht_abs(ijk)));
               break;
             }
           }
@@ -117,8 +117,8 @@ namespace libmpdataxx
           
           const auto &w = this->vips()[ct_params_t::n_dims - 1];
           this->state(ix::tht)(this->ijk) = ( this->state(ix::tht)(this->ijk) 
-                                            - 0.5 * this->dt * w(this->ijk) * this->dtht_e(this->ijk))
-                                            / (1 + 0.5 * this->dt * this->tht_abs(this->ijk));
+                                            - real_t(0.5) * this->dt * w(this->ijk) * this->dtht_e(this->ijk))
+                                            / (1 + real_t(0.5) * this->dt * this->tht_abs(this->ijk));
           this->rhs.at(ix::tht)(this->ijk) += -w(this->ijk) * this->dtht_e(this->ijk)
                                               -this->tht_abs(this->ijk) * this->state(ix::tht)(this->ijk);
         }
diff --git a/libmpdata++/solvers/detail/boussinesq_sgs_common.hpp b/libmpdata++/solvers/detail/boussinesq_sgs_common.hpp
index e25be230..f5dc2449 100644
--- a/libmpdata++/solvers/detail/boussinesq_sgs_common.hpp
+++ b/libmpdata++/solvers/detail/boussinesq_sgs_common.hpp
@@ -30,7 +30,7 @@ namespace libmpdataxx
         template <int nd = ct_params_t::n_dims> 
         void calc_rcdsn_num(typename std::enable_if<nd == 2>::type* = 0)
         {
-          rcdsn_num(this->ijk) = this->g * 0.5 * (
+          rcdsn_num(this->ijk) = this->g * real_t(0.5) * (
                                              grad_tht[ct_params_t::n_dims - 1](this->i, this->j - h)
                                            + grad_tht[ct_params_t::n_dims - 1](this->i, this->j + h)
                                            ) / (this->Tht_ref * tdef_sq(this->ijk));
@@ -39,7 +39,7 @@ namespace libmpdataxx
         template <int nd = ct_params_t::n_dims> 
         void calc_rcdsn_num(typename std::enable_if<nd == 3>::type* = 0)
         {
-          rcdsn_num(this->ijk) = this->g * 0.5 * (
+          rcdsn_num(this->ijk) = this->g * real_t(0.5) * (
                                              grad_tht[ct_params_t::n_dims - 1](this->i, this->j, this->k - h)
                                            + grad_tht[ct_params_t::n_dims - 1](this->i, this->j, this->k + h)
                                            ) / (this->Tht_ref * tdef_sq(this->ijk));
@@ -96,12 +96,16 @@ namespace libmpdataxx
           if (this->rank > 0)
             ijkm_aux[0] = this->ijk[0];
 
-          formulae::stress::multiply_tnsr_cmpct<ct_params_t::n_dims, ct_params_t::opts>(this->tau, 1.0, this->k_m, *this->mem->G, ijkm_aux);
+          formulae::stress::multiply_tnsr_cmpct<ct_params_t::n_dims, ct_params_t::opts>(this->tau,
+                                                                                        real_t(1.0),
+                                                                                        this->k_m,
+                                                                                        *this->mem->G,
+                                                                                        ijkm_aux);
 
           this->xchng_sgs_tnsr_offdiag(this->tau, this->tau_srfc, this->ijk, this->ijkm);
 
           formulae::stress::multiply_vctr_cmpct<ct_params_t::n_dims, ct_params_t::opts>(grad_tht,
-                                                                                        1.0 / prandtl_num,
+                                                                                        real_t(1.0) / prandtl_num,
                                                                                         this->k_m,
                                                                                         *this->mem->G,
                                                                                         this->ijk);
diff --git a/libmpdata++/solvers/detail/mpdata_rhs_vip_common.hpp b/libmpdata++/solvers/detail/mpdata_rhs_vip_common.hpp
index 4387f8f1..7f8c4fbf 100644
--- a/libmpdata++/solvers/detail/mpdata_rhs_vip_common.hpp
+++ b/libmpdata++/solvers/detail/mpdata_rhs_vip_common.hpp
@@ -43,16 +43,16 @@ namespace libmpdataxx
       template <class ct_params_t, int minhalo> 
       class mpdata_rhs_vip_common : public mpdata_rhs<ct_params_vip_default_t<ct_params_t>, minhalo>
       {
+        public:
         using ix = typename ct_params_t::ix;
-
-	protected:
-	
         using parent_t = mpdata_rhs<ct_params_vip_default_t<ct_params_t>, minhalo>;
+        using real_t = typename ct_params_t::real_t;
 
+	protected:
 	// member fields
         std::array<int, parent_t::n_dims> vip_ixs;
 	arrvec_t<typename parent_t::arr_t> &stash, &vip_rhs;
-        typename parent_t::real_t eps;
+        real_t eps;
 
         arrvec_t<typename parent_t::arr_t>& vip_stash(const int t_lev)
         {
@@ -189,7 +189,7 @@ namespace libmpdataxx
             }
             else
             {
-              vip_rhs[d](this->ijk) = 0.0;
+              vip_rhs[d](this->ijk) = 0;
             }
           }
         }
@@ -220,7 +220,7 @@ namespace libmpdataxx
             for (int d = 0; d < parent_t::n_dims; ++d)
             {
               vip_rhs[d](this->ijk) += vips()[d](this->ijk);
-              vip_rhs[d](this->ijk) /= (0.5 * this->dt);
+              vip_rhs[d](this->ijk) /= (real_t(0.5) * this->dt);
             }
           }
         }
@@ -229,7 +229,7 @@ namespace libmpdataxx
         {    
           for (int d = 0; d < parent_t::n_dims; ++d)
           {
-            vips()[d](this->ijk) += 0.5 * this->dt * vip_rhs[d](this->ijk);
+            vips()[d](this->ijk) += real_t(0.5) * this->dt * vip_rhs[d](this->ijk);
             vip_rhs[d](this->ijk) = 0;
           }
         }
@@ -240,7 +240,7 @@ namespace libmpdataxx
           {
             for (int d = 0; d < parent_t::n_dims; ++d)
             {
-              v[d](this->ijk) /= (1.0 + 0.5 * this->dt * (*this->mem->vab_coeff)(this->ijk));
+              v[d](this->ijk) /= (1 + real_t(0.5) * this->dt * (*this->mem->vab_coeff)(this->ijk));
             }
           }
         }
@@ -250,7 +250,7 @@ namespace libmpdataxx
           for (int d = 0; d < parent_t::n_dims; ++d)
           {
             this->vips()[d](this->ijk) +=
-              0.5 * this->dt * (*this->mem->vab_coeff)(this->ijk) * this->mem->vab_relax[d](this->ijk);
+              real_t(0.5) * this->dt * (*this->mem->vab_coeff)(this->ijk) * this->mem->vab_relax[d](this->ijk);
           }
         }
 
@@ -313,7 +313,7 @@ namespace libmpdataxx
                                                 - (this->dt_stash[1] + this->dt_stash[0]) * vip_stash(-1)[d](this->ijk)
                                                 + this->dt_stash[0] * vip_stash(-2)[d](this->ijk)
                                              ) / ( this->dt_stash[0] * this->dt_stash[1]
-                                                 * 0.5 * (this->dt_stash[0] + this->dt_stash[1])
+                                                 * real_t(0.5) * (this->dt_stash[0] + this->dt_stash[1])
                                              );
                 this->xchng_pres(this->vip_stash(0)[d], this->ijk, ex);
               }
@@ -349,7 +349,7 @@ namespace libmpdataxx
 	
         struct rt_params_t : parent_t::rt_params_t
         {
-          typename parent_t::real_t vip_eps = 0;
+          real_t vip_eps = 0;
         };
 
 	static void alloc(
diff --git a/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_common.hpp b/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_common.hpp
index 984f258e..ffbddcd6 100644
--- a/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_common.hpp
+++ b/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_common.hpp
@@ -87,7 +87,7 @@ namespace libmpdataxx
           int npoints = 1;
           for (int d = 0; d < parent_t::n_dims; ++d)
           {
-	    Phi(this->ijk) -= 0.5 * pow2(this->vips()[d](this->ijk));
+	    Phi(this->ijk) -= real_t(0.5) * pow2(this->vips()[d](this->ijk));
             npoints *= (this->mem->grid_size[d].last() + 1);
           }
           
@@ -184,7 +184,7 @@ namespace libmpdataxx
           for (int d = 0; d < parent_t::n_dims; ++d)
           {
             this->vip_rhs[d](this->ijk) += this->vips()[d](this->ijk);
-            this->vip_rhs[d](this->ijk) /= (0.5 * this->dt);
+            this->vip_rhs[d](this->ijk) /= (real_t(0.5) * this->dt);
           }
         }
 
diff --git a/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_common.hpp b/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_common.hpp
index 9104d333..aa061507 100644
--- a/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_common.hpp
+++ b/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_common.hpp
@@ -36,6 +36,7 @@ namespace libmpdataxx
       class mpdata_rhs_vip_prs_sgs_common : public mpdata_rhs_vip_prs<ct_params_t, minhalo>
       {
         using parent_t = mpdata_rhs_vip_prs<ct_params_t, minhalo>;
+        using real_t = typename ct_params_t::real_t;
 
         protected:
 
@@ -47,7 +48,7 @@ namespace libmpdataxx
         arrvec_t<typename parent_t::arr_t> &wrk;
 
         std::array<rng_t, ct_params_t::n_dims> ijkm;
-        typename parent_t::real_t cdrag;
+        real_t cdrag;
 
         virtual void multiply_sgs_visc() = 0;
 
@@ -61,7 +62,7 @@ namespace libmpdataxx
             formulae::stress::pade_dispatch<ct_params_t::n_dims>(wrk, this->ijk, d);
             // finish calculation of wrk[1] before modyfying wrk[0]
             this->mem->barrier();
-            wrk[0](this->ijk) += (drv[d](this->ijk) - 0.25 * wrk[1](this->ijk));
+            wrk[0](this->ijk) += (drv[d](this->ijk) - real_t(0.25) * wrk[1](this->ijk));
           }
           drv[d](this->ijk) = wrk[0](this->ijk);
           // needed because otherwise other threads could start calculating pade correction
@@ -115,7 +116,7 @@ namespace libmpdataxx
                                                                                             *this->mem->G,
                                                                                             this->ijk,
                                                                                             this->dijk,
-                                                                                            2.0);
+                                                                                            real_t(2.0));
           }
           else
           {
@@ -151,7 +152,7 @@ namespace libmpdataxx
             }
 
             // update forces
-            formulae::stress::calc_stress_rhs<ct_params_t::n_dims>(this->vip_rhs, drv, this->ijk, 2.0);
+            formulae::stress::calc_stress_rhs<ct_params_t::n_dims>(this->vip_rhs, drv, this->ijk, real_t(2.0));
           }
         }
 
diff --git a/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_smg.hpp b/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_smg.hpp
index 29451c6e..56d102ac 100644
--- a/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_smg.hpp
+++ b/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_smg.hpp
@@ -18,22 +18,23 @@ namespace libmpdataxx
       class mpdata_rhs_vip_prs_sgs_smg : public detail::mpdata_rhs_vip_prs_sgs_common<ct_params_t, minhalo>
       {
 	using parent_t = detail::mpdata_rhs_vip_prs_sgs_common<ct_params_t, minhalo>;
+        using real_t = typename ct_params_t::real_t;
 
         protected:
         
-	typename parent_t::real_t smg_c, c_m;
+	real_t smg_c, c_m;
         typename parent_t::arr_t &k_m;
 
         void multiply_sgs_visc()
         {
-          const auto dlta = std::accumulate(this->dijk.begin(), this->dijk.end(), 0.) / 3;
+          const auto dlta = std::accumulate(this->dijk.begin(), this->dijk.end(), real_t(0.)) / 3;
 
           if (static_cast<stress_diff_t>(ct_params_t::stress_diff) == compact)
           {
             k_m(this->ijk) = pow(smg_c * dlta, 2) * formulae::stress::calc_tdef_sq_cmpct<ct_params_t::n_dims>(this->tau, this->ijk);
 
             formulae::stress::multiply_tnsr_cmpct<ct_params_t::n_dims, ct_params_t::opts>(this->tau,
-                                                                                          1.0,
+                                                                                          real_t(1.0),
                                                                                           this->k_m,
                                                                                           *this->mem->G,
                                                                                           this->ijk);
@@ -56,7 +57,7 @@ namespace libmpdataxx
 
         struct rt_params_t : parent_t::rt_params_t
         {
-          typename ct_params_t::real_t smg_c, c_m;
+          real_t smg_c, c_m;
         };
 
 	// ctor
diff --git a/libmpdata++/solvers/detail/solver_1d.hpp b/libmpdata++/solvers/detail/solver_1d.hpp
index 0c2bef04..9ca9939f 100644
--- a/libmpdata++/solvers/detail/solver_1d.hpp
+++ b/libmpdata++/solvers/detail/solver_1d.hpp
@@ -67,7 +67,7 @@ namespace libmpdataxx
 
         typename parent_t::real_t courant_number(const arrvec_t<typename parent_t::arr_t> &arrvec) final
         {
-          stat_field(this->ijk) = 0.5 * (abs(arrvec[0](i+h) + arrvec[0](i-h)));
+          stat_field(this->ijk) = typename parent_t::real_t(0.5) * (abs(arrvec[0](i+h) + arrvec[0](i-h)));
           return this->mem->max(this->rank, stat_field(this->ijk));
         }
         
diff --git a/libmpdata++/solvers/detail/solver_2d.hpp b/libmpdata++/solvers/detail/solver_2d.hpp
index ff385014..2f63a84a 100644
--- a/libmpdata++/solvers/detail/solver_2d.hpp
+++ b/libmpdata++/solvers/detail/solver_2d.hpp
@@ -200,7 +200,7 @@ namespace libmpdataxx
 
         typename parent_t::real_t courant_number(const arrvec_t<typename parent_t::arr_t> &arrvec) final
         {
-          stat_field(this->ijk) = 0.5 * (
+          stat_field(this->ijk) = typename parent_t::real_t(0.5) * (
                                            abs(arrvec[0](i+h, j) + arrvec[0](i-h, j))
                                          + abs(arrvec[1](i, j+h) + arrvec[1](i, j-h))
                                         ) / formulae::G<ct_params_t::opts, 0>(*this->mem->G, i, j);
diff --git a/libmpdata++/solvers/detail/solver_3d.hpp b/libmpdata++/solvers/detail/solver_3d.hpp
index 7b0dc051..895211f6 100644
--- a/libmpdata++/solvers/detail/solver_3d.hpp
+++ b/libmpdata++/solvers/detail/solver_3d.hpp
@@ -241,7 +241,7 @@ namespace libmpdataxx
 
         typename parent_t::real_t courant_number(const arrvec_t<typename parent_t::arr_t> &arrvec) final
         {
-          stat_field(this->ijk) =  0.5 * (
+          stat_field(this->ijk) =  typename parent_t::real_t(0.5) * (
                                             abs(arrvec[0](i+h, j, k) + arrvec[0](i-h, j, k))
                                           + abs(arrvec[1](i, j+h, k) + arrvec[1](i, j-h, k))
                                           + abs(arrvec[2](i, j, k+h) + arrvec[2](i, j, k-h))
diff --git a/libmpdata++/solvers/detail/solver_common.hpp b/libmpdata++/solvers/detail/solver_common.hpp
index e7224263..5c404a32 100644
--- a/libmpdata++/solvers/detail/solver_common.hpp
+++ b/libmpdata++/solvers/detail/solver_common.hpp
@@ -209,7 +209,7 @@ namespace libmpdataxx
         struct rt_params_t 
         {
           std::array<int, n_dims> grid_size;
-          real_t dt=0, max_abs_div_eps = blitz::epsilon(real_t(44)), max_courant = 0.5;
+          real_t dt=0, max_abs_div_eps = blitz::epsilon(real_t(44)), max_courant = real_t(0.5);
         };
 
 	// ctor
diff --git a/libmpdata++/solvers/mpdata_rhs_vip.hpp b/libmpdata++/solvers/mpdata_rhs_vip.hpp
index 9ea85205..e363a775 100644
--- a/libmpdata++/solvers/mpdata_rhs_vip.hpp
+++ b/libmpdata++/solvers/mpdata_rhs_vip.hpp
@@ -40,9 +40,11 @@ namespace libmpdataxx
                                 const arrvec_t<typename parent_t::arr_t> &src) final
       {
         using namespace libmpdataxx::arakawa_c;
+        using real_t = typename ct_params_t::real_t;
+
 	if (!this->mem->G)
 	{
-	  dst[0](im + h) = this->dt / this->di * .5 * (
+	  dst[0](im + h) = this->dt / this->di * real_t(.5) * (
 	    src[0](im    ) + 
 	    src[0](im + 1)
 	  );
@@ -104,17 +106,18 @@ namespace libmpdataxx
       {   
 	using idxperm::pi;
 	using namespace arakawa_c;
+        using real_t = typename ct_params_t::real_t;
   
 	if (!this->mem->G)
 	{
-	  dst(pi<d>(i+h,j)) = this->dt / di * .5 * (
+	  dst(pi<d>(i+h,j)) = this->dt / di * real_t(.5) * (
 	    src(pi<d>(i,    j)) + 
 	    src(pi<d>(i + 1,j))
 	  );
 	} 
 	else
 	{ 
-	  dst(pi<d>(i+h,j)) = this->dt / di * .5 * (
+	  dst(pi<d>(i+h,j)) = this->dt / di * real_t(.5) * (
 	    (*this->mem->G)(pi<d>(i,    j)) * src(pi<d>(i,    j)) + 
 	    (*this->mem->G)(pi<d>(i + 1,j)) * src(pi<d>(i + 1,j))
 	  );
@@ -214,17 +217,18 @@ namespace libmpdataxx
       {   
 	using idxperm::pi;
 	using namespace arakawa_c;
+        using real_t = typename ct_params_t::real_t;
   
 	if (!this->mem->G)
 	{
-	  dst(pi<d>(i+h, j, k)) = this->dt / di * .5 * (
+	  dst(pi<d>(i+h, j, k)) = this->dt / di * real_t(.5) * (
 	    src(pi<d>(i,     j, k)) + 
 	    src(pi<d>(i + 1, j, k))
 	  );
 	} 
 	else
 	{ 
-	  dst(pi<d>(i+h, j, k)) = this->dt / di * .5 * (
+	  dst(pi<d>(i+h, j, k)) = this->dt / di * real_t(.5) * (
 	    (*this->mem->G)(pi<d>(i  , j, k)) * src(pi<d>(i,   j, k)) + 
 	    (*this->mem->G)(pi<d>(i+1, j, k)) * src(pi<d>(i+1, j, k))
 	  );

From 94173c77a9f391078a2f95f455b835b4f5deb9d9 Mon Sep 17 00:00:00 2001
From: Maciej Waruszewski <mwarusz@igf.fuw.edu.pl>
Date: Wed, 13 Mar 2019 12:16:28 +0100
Subject: [PATCH 2/4] use a different typedef for real_t

---
 libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_common.hpp | 2 +-
 libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_dns.hpp    | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_common.hpp b/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_common.hpp
index aa061507..ce1395e2 100644
--- a/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_common.hpp
+++ b/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_common.hpp
@@ -160,7 +160,7 @@ namespace libmpdataxx
 
         struct rt_params_t : parent_t::rt_params_t 
         { 
-          typename parent_t::real_t cdrag = 0; 
+          real_t cdrag = 0; 
         };
 
         // ctor
diff --git a/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_dns.hpp b/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_dns.hpp
index f22e5dd7..021efcae 100644
--- a/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_dns.hpp
+++ b/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_dns.hpp
@@ -21,7 +21,7 @@ namespace libmpdataxx
 
         protected:
         
-	typename parent_t::real_t eta;
+	typename ct_params_t::real_t eta;
 
         void multiply_sgs_visc() final
         {

From da38808867080008542f8c005122a35f42a145b3 Mon Sep 17 00:00:00 2001
From: Maciej Waruszewski <mwarusz@igf.fuw.edu.pl>
Date: Wed, 13 Mar 2019 12:31:12 +0100
Subject: [PATCH 3/4] always use double in sums

---
 libmpdata++/concurr/detail/sharedmem.hpp | 13 ++++++-------
 1 file changed, 6 insertions(+), 7 deletions(-)

diff --git a/libmpdata++/concurr/detail/sharedmem.hpp b/libmpdata++/concurr/detail/sharedmem.hpp
index 20be0319..d90a6784 100644
--- a/libmpdata++/concurr/detail/sharedmem.hpp
+++ b/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

From b79e7f3053e63a2bdd29c123d2b2a2449e88e648 Mon Sep 17 00:00:00 2001
From: Maciej Waruszewski <mwarusz@igf.fuw.edu.pl>
Date: Wed, 13 Mar 2019 15:35:45 +0100
Subject: [PATCH 4/4] make real_t typedefs public

---
 .../solvers/detail/mpdata_rhs_vip_common.hpp  |  3 ++-
 .../detail/mpdata_rhs_vip_prs_sgs_common.hpp  |  3 +++
 .../detail/mpdata_rhs_vip_prs_sgs_smg.hpp     |  3 +++
 libmpdata++/solvers/detail/solver_1d.hpp      | 18 ++++++++++-------
 libmpdata++/solvers/detail/solver_2d.hpp      | 20 +++++++++++--------
 libmpdata++/solvers/detail/solver_3d.hpp      | 20 +++++++++++--------
 libmpdata++/solvers/detail/solver_common.hpp  |  2 +-
 7 files changed, 44 insertions(+), 25 deletions(-)

diff --git a/libmpdata++/solvers/detail/mpdata_rhs_vip_common.hpp b/libmpdata++/solvers/detail/mpdata_rhs_vip_common.hpp
index 7f8c4fbf..b0898474 100644
--- a/libmpdata++/solvers/detail/mpdata_rhs_vip_common.hpp
+++ b/libmpdata++/solvers/detail/mpdata_rhs_vip_common.hpp
@@ -43,9 +43,10 @@ namespace libmpdataxx
       template <class ct_params_t, int minhalo> 
       class mpdata_rhs_vip_common : public mpdata_rhs<ct_params_vip_default_t<ct_params_t>, minhalo>
       {
+        using parent_t = mpdata_rhs<ct_params_vip_default_t<ct_params_t>, minhalo>;
+
         public:
         using ix = typename ct_params_t::ix;
-        using parent_t = mpdata_rhs<ct_params_vip_default_t<ct_params_t>, minhalo>;
         using real_t = typename ct_params_t::real_t;
 
 	protected:
diff --git a/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_common.hpp b/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_common.hpp
index ce1395e2..fd53692a 100644
--- a/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_common.hpp
+++ b/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_common.hpp
@@ -36,6 +36,9 @@ namespace libmpdataxx
       class mpdata_rhs_vip_prs_sgs_common : public mpdata_rhs_vip_prs<ct_params_t, minhalo>
       {
         using parent_t = mpdata_rhs_vip_prs<ct_params_t, minhalo>;
+
+        public:
+
         using real_t = typename ct_params_t::real_t;
 
         protected:
diff --git a/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_smg.hpp b/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_smg.hpp
index 56d102ac..aff23fba 100644
--- a/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_smg.hpp
+++ b/libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_smg.hpp
@@ -18,6 +18,9 @@ namespace libmpdataxx
       class mpdata_rhs_vip_prs_sgs_smg : public detail::mpdata_rhs_vip_prs_sgs_common<ct_params_t, minhalo>
       {
 	using parent_t = detail::mpdata_rhs_vip_prs_sgs_common<ct_params_t, minhalo>;
+
+        public:
+
         using real_t = typename ct_params_t::real_t;
 
         protected:
diff --git a/libmpdata++/solvers/detail/solver_1d.hpp b/libmpdata++/solvers/detail/solver_1d.hpp
index 9ca9939f..f579a6bf 100644
--- a/libmpdata++/solvers/detail/solver_1d.hpp
+++ b/libmpdata++/solvers/detail/solver_1d.hpp
@@ -25,6 +25,10 @@ namespace libmpdataxx
       {
 	using parent_t = solver_common<ct_params_t, n_tlev, minhalo>;
 
+        public:
+
+	using real_t = typename ct_params_t::real_t;
+
 	protected:
 
 	const rng_t i; //TODO: to be removed
@@ -65,21 +69,21 @@ namespace libmpdataxx
           this->mem->barrier();
         }
 
-        typename parent_t::real_t courant_number(const arrvec_t<typename parent_t::arr_t> &arrvec) final
+        real_t courant_number(const arrvec_t<typename parent_t::arr_t> &arrvec) final
         {
-          stat_field(this->ijk) = typename parent_t::real_t(0.5) * (abs(arrvec[0](i+h) + arrvec[0](i-h)));
+          stat_field(this->ijk) = real_t(0.5) * (abs(arrvec[0](i+h) + arrvec[0](i-h)));
           return this->mem->max(this->rank, stat_field(this->ijk));
         }
         
-        typename parent_t::real_t max_abs_vctr_div(const arrvec_t<typename parent_t::arr_t> &arrvec) final
+        real_t max_abs_vctr_div(const arrvec_t<typename parent_t::arr_t> &arrvec) final
         {
           stat_field(this->ijk) = abs((arrvec[0](i+h) - arrvec[0](i-h)));
           return this->mem->max(this->rank, stat_field(this->ijk));
         }
 
-        void scale_gc(const typename parent_t::real_t time,
-                      const typename parent_t::real_t cur_dt,
-                      const typename parent_t::real_t old_dt) final
+        void scale_gc(const real_t time,
+                      const real_t cur_dt,
+                      const real_t old_dt) final
         {
           this->mem->GC[0](rng_t(i.first(), i.last()-1)^h) *= cur_dt / old_dt;
           this->xchng_vctr_alng(this->mem->GC);
@@ -99,7 +103,7 @@ namespace libmpdataxx
 
         struct rt_params_t : parent_t::rt_params_t
         {
-          typename parent_t::real_t di = 0;
+          real_t di = 0;
         };
 
         protected:
diff --git a/libmpdata++/solvers/detail/solver_2d.hpp b/libmpdata++/solvers/detail/solver_2d.hpp
index 2f63a84a..82c216a3 100644
--- a/libmpdata++/solvers/detail/solver_2d.hpp
+++ b/libmpdata++/solvers/detail/solver_2d.hpp
@@ -26,6 +26,10 @@ namespace libmpdataxx
       {
 	using parent_t = solver_common<ct_params_t, n_tlev, minhalo>;
 
+        public:
+
+	using real_t = typename ct_params_t::real_t;
+
 	protected:
       
 	const rng_t i, j; // TODO: to be removed
@@ -191,23 +195,23 @@ namespace libmpdataxx
           // TODO: same in 1D
           if (!opts::isset(ct_params_t::opts, opts::dfl))
           {
-            typename ct_params_t::real_t max_abs_div = max_abs_vctr_div(this->mem->GC);
+            real_t max_abs_div = max_abs_vctr_div(this->mem->GC);
 
 	    if (max_abs_div > this->max_abs_div_eps) 
 	      throw std::runtime_error("initial advector field is divergent");
           }
         }
 
-        typename parent_t::real_t courant_number(const arrvec_t<typename parent_t::arr_t> &arrvec) final
+        real_t courant_number(const arrvec_t<typename parent_t::arr_t> &arrvec) final
         {
-          stat_field(this->ijk) = typename parent_t::real_t(0.5) * (
+          stat_field(this->ijk) = real_t(0.5) * (
                                            abs(arrvec[0](i+h, j) + arrvec[0](i-h, j))
                                          + abs(arrvec[1](i, j+h) + arrvec[1](i, j-h))
                                         ) / formulae::G<ct_params_t::opts, 0>(*this->mem->G, i, j);
           return this->mem->max(this->rank, stat_field(this->ijk));
         }
         
-        typename parent_t::real_t max_abs_vctr_div(const arrvec_t<typename parent_t::arr_t> &arrvec) final
+        real_t max_abs_vctr_div(const arrvec_t<typename parent_t::arr_t> &arrvec) final
         {
           stat_field(this->ijk) = abs(
                                         (arrvec[0](i+h, j) - arrvec[0](i-h, j))
@@ -216,9 +220,9 @@ namespace libmpdataxx
           return this->mem->max(this->rank, stat_field(this->ijk));
         }
         
-        void scale_gc(const typename parent_t::real_t time,
-                      const typename parent_t::real_t cur_dt,
-                      const typename parent_t::real_t old_dt) final
+        void scale_gc(const real_t time,
+                      const real_t cur_dt,
+                      const real_t old_dt) final
         {
           this->mem->GC[0](rng_t(i.first(), i.last()-1)^h, j) *= cur_dt / old_dt;
           this->mem->GC[1](i, rng_t(j.first(), j.last()-1)^h) *= cur_dt / old_dt;
@@ -241,7 +245,7 @@ namespace libmpdataxx
         
         struct rt_params_t : parent_t::rt_params_t
         {
-          typename parent_t::real_t di = 0, dj = 0;
+          real_t di = 0, dj = 0;
         };
 
         protected:
diff --git a/libmpdata++/solvers/detail/solver_3d.hpp b/libmpdata++/solvers/detail/solver_3d.hpp
index 895211f6..595bf8e8 100644
--- a/libmpdata++/solvers/detail/solver_3d.hpp
+++ b/libmpdata++/solvers/detail/solver_3d.hpp
@@ -26,6 +26,10 @@ namespace libmpdataxx
       {
 	using parent_t = solver_common<ct_params_t, n_tlev, minhalo>;
 
+        public:
+
+	using real_t = typename ct_params_t::real_t;
+
 	protected:
 
 	const rng_t i, j, k; // TODO: we have ijk in solver_common - could it be removed?
@@ -232,16 +236,16 @@ namespace libmpdataxx
           // TODO: same in 1D
           if (!opts::isset(ct_params_t::opts, opts::dfl))
           {
-            typename ct_params_t::real_t max_abs_div = max_abs_vctr_div(this->mem->GC);
+            real_t max_abs_div = max_abs_vctr_div(this->mem->GC);
 
 	    if (max_abs_div > this->max_abs_div_eps) 
 	      throw std::runtime_error("initial advector field is divergent");
           }
         }
 
-        typename parent_t::real_t courant_number(const arrvec_t<typename parent_t::arr_t> &arrvec) final
+        real_t courant_number(const arrvec_t<typename parent_t::arr_t> &arrvec) final
         {
-          stat_field(this->ijk) =  typename parent_t::real_t(0.5) * (
+          stat_field(this->ijk) = real_t(0.5) * (
                                             abs(arrvec[0](i+h, j, k) + arrvec[0](i-h, j, k))
                                           + abs(arrvec[1](i, j+h, k) + arrvec[1](i, j-h, k))
                                           + abs(arrvec[2](i, j, k+h) + arrvec[2](i, j, k-h))
@@ -249,7 +253,7 @@ namespace libmpdataxx
           return this->mem->max(this->rank, stat_field(this->ijk));
         }
         
-        typename parent_t::real_t max_abs_vctr_div(const arrvec_t<typename parent_t::arr_t> &arrvec) final
+        real_t max_abs_vctr_div(const arrvec_t<typename parent_t::arr_t> &arrvec) final
         {
           stat_field(this->ijk) =  abs(
                                          (arrvec[0](i+h, j, k) - arrvec[0](i-h, j, k))
@@ -260,9 +264,9 @@ namespace libmpdataxx
           return this->mem->max(this->rank, stat_field(this->ijk));
         }
 
-        void scale_gc(const typename parent_t::real_t time,
-                      const typename parent_t::real_t cur_dt,
-                      const typename parent_t::real_t old_dt) final
+        void scale_gc(const real_t time,
+                      const real_t cur_dt,
+                      const real_t old_dt) final
         {
           this->mem->GC[0](rng_t(i.first(), i.last()-1)^h, j, k) *= cur_dt / old_dt;
           this->mem->GC[1](i, rng_t(j.first(), j.last()-1)^h, k) *= cur_dt / old_dt;
@@ -289,7 +293,7 @@ namespace libmpdataxx
 
         struct rt_params_t : parent_t::rt_params_t
         {
-          typename parent_t::real_t di = 0, dj = 0, dk = 0;
+          real_t di = 0, dj = 0, dk = 0;
         };
 
         protected:
diff --git a/libmpdata++/solvers/detail/solver_common.hpp b/libmpdata++/solvers/detail/solver_common.hpp
index 5c404a32..0dcc8074 100644
--- a/libmpdata++/solvers/detail/solver_common.hpp
+++ b/libmpdata++/solvers/detail/solver_common.hpp
@@ -40,7 +40,7 @@ namespace libmpdataxx
         enum { n_tlev = n_tlev_ };
         
         using ct_params_t_ = ct_params_t; // propagate ct_params_t mainly for output purposes
-        typedef typename ct_params_t::real_t real_t;
+        using real_t = typename ct_params_t::real_t;
         typedef blitz::Array<real_t, n_dims> arr_t;
         using bcp_t = std::unique_ptr<bcond::detail::bcond_common<real_t, halo>>;