Merge pull request #388 from pdziekan/cache_mfp

Cache mfp

Author: pdziekan

.travis_scripts/cuda.sh

@@ -6,7 +6,8 @@ cd build
 #if [[ $TRAVIS_OS_NAME == 'linux' && $CXX == 'clang++' ]]; then cmake ../; fi 
 if [[ $TRAVIS_OS_NAME == 'osx' ]]; then cmake .. -DPYTHON_LIBRARY=$PY_LIB -DPYTHON_INCLUDE_DIR=$PY_INC; fi
 cmake -DCMAKE_BUILD_TYPE=Debug ../
-VERBOSE=1 travis_wait 20 make -j2
+#VERBOSE=1 travis_wait 20 make -j2
+VERBOSE=1 make -j2
 cmake -DCMAKE_BUILD_TYPE=Release ../ 
 VERBOSE=1 travis_wait 20 make -j2
 cd ../..

include/libcloudph++/blk_2m/activation_formulae.hpp

@@ -36,8 +36,8 @@ namespace libcloudphxx
         const quantity<si::dimensionless, real_t> &chem_b,
         const quantity<si::dimensionless, real_t> beta = beta_default<real_t>()
       ) {
-        return pow(mean_rd / si::metres, -(1+beta)) * sqrt(4 * pow(common::kelvin::A<real_t>(T) / si::metres, 3) / 27 / chem_b);
-      }           // can't have pow<3/2>(si:: )
+        return std::pow(mean_rd / si::metres, -(1+beta)) * std::sqrt(4 * std::pow(common::kelvin::A<real_t>(T) / si::metres, 3) / 27 / chem_b);
+      }           // can't have std::pow<3/2>(si:: )
 
       // helper for activation formulae (see eq. 10 in Morrison and Grabowski 2007)
       template <typename real_t>
@@ -55,7 +55,7 @@ namespace libcloudphxx
         const quantity<si::dimensionless, real_t> &sdev_rd,
         const quantity<si::dimensionless, real_t> beta = beta_default<real_t>()
       ) {
-        return pow(sdev_rd, 1+beta);
+        return std::pow(sdev_rd, 1+beta);
       }
 
       //helper for activation formulae (see eq. 10 in Morrison and Grabowski 2007)
@@ -73,7 +73,7 @@ namespace libcloudphxx
         return log(
           s_0(T, mean_rd, chem_b) / 
           std::min(real_t(s(p, T, rv)), real_t(RH_max - 1))
-        ) / sqrt(2) / log(sdev_rd_s(sdev_rd));
+        ) / std::sqrt(2) / log(sdev_rd_s(sdev_rd));
       }
 
       // helper for activation formulae (see eq. 10 in Morrison and Grabowski 2007)

include/libcloudph++/blk_2m/common_formulae.hpp

@@ -48,7 +48,8 @@ namespace libcloudphxx
       inline quantity<si::dimensionless, real_t> miu_c(
         const quantity<divide_typeof_helper<si::dimensionless, si::volume>::type, real_t> &n
       ) {
-        auto tmp = real_t(1) / pow<2>(eta(n)) - real_t(1);
+        auto eta_n = eta(n);
+        auto tmp = real_t(1) / (eta_n * eta_n) - real_t(1);
         assert(finite(tmp) && "spectral index n is finite failed");
         return tmp;
       }
@@ -60,7 +61,7 @@ namespace libcloudphxx
          const quantity<si::dimensionless, real_t> &rc,
         const quantity<si::mass_density, real_t> &rhod
       ) {
-        auto tmp = pow(
+        auto tmp = std::pow(
 	  c_md<real_t>() * n * std::tgamma(miu_c(n*rhod) + d_md<real_t>() + real_t(1))
           /
 	  (rc * std::tgamma(miu_c(n*rhod) + real_t(1))) * si::cubic_metres
@@ -82,7 +83,7 @@ namespace libcloudphxx
          const quantity<si::mass_density, real_t> &rhod
       ) {
         auto tmp = nc
-          * pow(lambda_c(nc, rc, rhod) * si::metres, miu_c(nc * rhod) + real_t(1))
+          * std::pow(lambda_c(nc, rc, rhod) * si::metres, miu_c(nc * rhod) + real_t(1))
           / std::tgamma(miu_c(nc * rhod) + real_t(1))
           / si::metres;
         assert(finite(tmp * si::metres * si::cubic_metres) && "intercept N0_c is finite failed");
@@ -97,7 +98,7 @@ namespace libcloudphxx
          const quantity<divide_typeof_helper<si::dimensionless, si::mass>::type, real_t> &nr,
          const quantity<si::dimensionless, real_t> &rr
       ) {
-        auto tmp = pow(
+        auto tmp = std::pow(
 	  c_md<real_t>() * nr * std::tgamma(d_md<real_t>() + real_t(1)) / rr * si::cubic_metres
           ,
 	  real_t(1) / d_md<real_t>()

include/libcloudph++/blk_2m/rhs_cellwise.hpp

@@ -215,7 +215,8 @@ namespace libcloudphxx
 
               // sink of N for cloud droplets is combined with the sink due to accretion
               // source of N for drizzle assumes that all the drops have the same radius
-              local_dot_nr += tmp / (real_t(4)/3 * pi<real_t>() * rho_w<real_t>() * pow<3>(drizzle_radius<real_t>()))
+              local_dot_nr += tmp / (real_t(4)/3 * pi<real_t>() * rho_w<real_t>() *
+                drizzle_radius<real_t>() * drizzle_radius<real_t>() * drizzle_radius<real_t>())
                 * si::kilograms * si::seconds; // to make it dimensionless
 
               if(tmp == rc / (dt * si::seconds)) // all cloud water turned into rain

include/libcloudph++/blk_2m/terminal_vel_formulae.hpp

@@ -60,8 +60,8 @@ namespace libcloudphxx
           const quantity<divide_typeof_helper<si::dimensionless, si::length>::type, real_t> &lbd, //slope of assumed exponential size distribution
           const quantity<si::length, real_t> &D // diameter
         ) {
-          auto tmp =  - pow(lbd * si::metres, -6) * exp(-lbd * D) *
-            (pow(real_t(lbd * D), 5) + 5 * pow(real_t(lbd * D), 4) * 20 * pow(real_t(lbd * D), 3) + 60 * pow(real_t(lbd * D), 2) + 120 * real_t(lbd * D) + 120);
+          auto tmp =  - std::pow(lbd * si::metres, -6) * exp(-lbd * D) *
+            (pow(real_t(lbd * D), 5) + 5 * std::pow(real_t(lbd * D), 4) * 20 * std::pow(real_t(lbd * D), 3) + 60 * std::pow(real_t(lbd * D), 2) + 120 * real_t(lbd * D) + 120);
 
           assert(finite(tmp) && "mint_1 is finite failed");
           return tmp;
@@ -72,8 +72,8 @@ namespace libcloudphxx
           const quantity<divide_typeof_helper<si::dimensionless, si::length>::type, real_t> &lbd, //slope of assumed exponential size distribution
           const quantity<si::length, real_t> &D // diameter
         ) {
-          auto tmp = - pow(real_t(lbd * si::metres), -5) * exp(-lbd * D) *
-            (pow(real_t(lbd * D), 4) + 4 * pow(real_t(lbd * D), 3) * 12 * pow(real_t(lbd * D), 2) + 24 * real_t(lbd * D) + 24);
+          auto tmp = - std::pow(real_t(lbd * si::metres), -5) * exp(-lbd * D) *
+            (pow(real_t(lbd * D), 4) + 4 * std::pow(real_t(lbd * D), 3) * 12 * std::pow(real_t(lbd * D), 2) + 24 * real_t(lbd * D) + 24);
 
           assert(finite(tmp) && "mint_2 is finite failed");
           return tmp;
@@ -84,9 +84,9 @@ namespace libcloudphxx
           const quantity<divide_typeof_helper<si::dimensionless, si::length>::type, real_t> &lbd, // slope of assumed exponential size distribution
           const quantity<si::length, real_t> &D // diameter
         ) {
-          auto tmp = real_t(1./16) / pow(real_t(lbd * si::metres), real_t(9./2))
-                   * (105 * sqrt(pi<real_t>()) * std::erf(sqrt(lbd * D))
-                      - 2 * sqrt(lbd * D) * exp(-lbd * D) * (8 * pow(real_t(lbd * D), 3) + 28 * pow(real_t(lbd * D), 2) + 70 * real_t(lbd * D) + 105)
+          auto tmp = real_t(1./16) / std::pow(real_t(lbd * si::metres), real_t(9./2))
+                   * (105 * std::sqrt(pi<real_t>()) * std::erf(sqrt(lbd * D))
+                      - 2 * std::sqrt(lbd * D) * exp(-lbd * D) * (8 * std::pow(real_t(lbd * D), 3) + 28 * std::pow(real_t(lbd * D), 2) + 70 * real_t(lbd * D) + 105)
                    );
 
           assert(finite(tmp) && "mint_3 is finite failed");
@@ -99,7 +99,7 @@ namespace libcloudphxx
           const quantity<divide_typeof_helper<si::dimensionless, si::length>::type, real_t> &lbd, //slope of assumed exponential size distribution
           const quantity<si::length, real_t> &D // diameter
         ) {
-          auto tmp = - pow(real_t(lbd * si::metres), -1) * exp(-lbd * D);
+          auto tmp = - std::pow(real_t(lbd * si::metres), -1) * exp(-lbd * D);
 
           assert(finite(tmp) && "int_4 is finite failed");
           return tmp;
@@ -112,7 +112,7 @@ namespace libcloudphxx
          const quantity<divide_typeof_helper<si::dimensionless, si::length>::type, real_t> &lbd, // slope of assumed exponential size distribution
          const quantity<si::length, real_t> &D // diameter
        ) {
-         auto tmp = pow(real_t(lbd * si::metres), -3) * exp(-lbd * D) * (- real_t(lbd * D) * (real_t(lbd * D) + 2) - 2);
+         auto tmp = std::pow(real_t(lbd * si::metres), -3) * exp(-lbd * D) * (- real_t(lbd * D) * (real_t(lbd * D) + 2) - 2);
 
          assert(finite(tmp) && "nint_1 is finite failed");
          return tmp;
@@ -124,7 +124,7 @@ namespace libcloudphxx
          const quantity<divide_typeof_helper<si::dimensionless, si::length>::type, real_t> &lbd, //slope of assumed exponential size distribution
          const quantity<si::length, real_t> &D // diameter
        ) {
-         auto tmp =  - pow(real_t(lbd * si::metres), -2) * exp(-lbd * D) * (lbd * D + 1);
+         auto tmp =  - std::pow(real_t(lbd * si::metres), -2) * exp(-lbd * D) * (lbd * D + 1);
 
          assert(finite(tmp) && "nint_2 is finite failed");
          return tmp;
@@ -136,8 +136,8 @@ namespace libcloudphxx
          const quantity<divide_typeof_helper<si::dimensionless, si::length>::type, real_t> &lbd, //slope of assumed exponential size distribution
          const quantity<si::length, real_t> &D //diameter
        ) {
-         auto tmp = sqrt(pi<real_t>()) * std::erf(sqrt(lbd * D)) / 2 / pow(real_t(lbd * si::metres), 3./2)
-           - sqrt(D / si::metres) * exp(-lbd*D) / (lbd * si::metres);
+         auto tmp = std::sqrt(pi<real_t>()) * std::erf(sqrt(lbd * D)) / 2 / std::pow(real_t(lbd * si::metres), 3./2)
+           - std::sqrt(D / si::metres) * exp(-lbd*D) / (lbd * si::metres);
 
          assert(finite(tmp) && "nint_3 is finite failed");
          return tmp;
@@ -161,15 +161,15 @@ namespace libcloudphxx
           * real_t(1000) // mass of the drop in grams
           * (
 	    alpha_fall(d1<real_t>() / real_t(2))
-            * pow(c_md<real_t>() * si::cubic_metres / si::kilograms * real_t(1000), beta_fall(d1<real_t>() / real_t(2)))
+            * std::pow(c_md<real_t>() * si::cubic_metres / si::kilograms * real_t(1000), beta_fall(d1<real_t>() / real_t(2)))
 	    * (mint_1(lbd, d1<real_t>()) - mint_1(lbd, real_t(0) * si::metres))
 	    +
 	    alpha_fall(d1<real_t>() + d2<real_t>() / real_t(2))
-            * pow(c_md<real_t>() * si::cubic_metres / si::kilograms * real_t(1000), beta_fall(d1<real_t>() + d2<real_t>() / real_t(2)))
+            * std::pow(c_md<real_t>() * si::cubic_metres / si::kilograms * real_t(1000), beta_fall(d1<real_t>() + d2<real_t>() / real_t(2)))
 	    * (mint_2(lbd, d2<real_t>()) - mint_2(lbd, d1<real_t>()))
 	    +
 	    alpha_fall(d2<real_t>() + d3<real_t>() / real_t(2))
-            * pow(c_md<real_t>() * si::cubic_metres / si::kilograms * real_t(1000), beta_fall(d2<real_t>() + d3<real_t>() / real_t(2)))
+            * std::pow(c_md<real_t>() * si::cubic_metres / si::kilograms * real_t(1000), beta_fall(d2<real_t>() + d3<real_t>() / real_t(2)))
 	    * (mint_3(lbd, d3<real_t>()) - mint_3(lbd, d2<real_t>()))
 	    +
 	    alpha_fall(real_t(2) * d3<real_t>())
@@ -197,15 +197,15 @@ namespace libcloudphxx
         auto tmp = rho_stp<real_t>() / rhod
           * (
 	   alpha_fall(d1<real_t>() / real_t(2))
-           * pow(c_md<real_t>() * si::cubic_metres / si::kilograms * real_t(1000), beta_fall(d1<real_t>() / real_t(2)))
+           * std::pow(c_md<real_t>() * si::cubic_metres / si::kilograms * real_t(1000), beta_fall(d1<real_t>() / real_t(2)))
            * (nint_1(lbd, d1<real_t>()) - nint_1(lbd, real_t(0) * si::metres))
 	   +
            alpha_fall(d1<real_t>() + d2<real_t>() / real_t(2))
-           * pow(c_md<real_t>() * si::cubic_metres / si::kilograms * real_t(1000), beta_fall(d1<real_t>() + d2<real_t>() / real_t(2)))
+           * std::pow(c_md<real_t>() * si::cubic_metres / si::kilograms * real_t(1000), beta_fall(d1<real_t>() + d2<real_t>() / real_t(2)))
            * (nint_2(lbd, d2<real_t>()) - nint_2(lbd, d1<real_t>()))
            +
            alpha_fall(d2<real_t>() + d3<real_t>() / real_t(2))
-           * pow(c_md<real_t>() * si::cubic_metres / si::kilograms * real_t(1000), beta_fall(d2<real_t>() + d3<real_t>() / real_t(2)))
+           * std::pow(c_md<real_t>() * si::cubic_metres / si::kilograms * real_t(1000), beta_fall(d2<real_t>() + d3<real_t>() / real_t(2)))
            * (nint_3(lbd, d3<real_t>()) - nint_3(lbd, d2<real_t>()))
            +
            alpha_fall(real_t(2) * d3<real_t>()) * (real_t(0) - int_4(lbd, d3<real_t>()))

include/libcloudph++/common/GA17_turbulence.hpp

@@ -7,6 +7,13 @@ namespace libcloudphxx
     // Grabowski and Abade 2017
     namespace GA17_turbulence
     {
+
+#if !defined(__NVCC__)
+      using std::pow;
+      using std::cbrt;
+      using std::sqrt;
+#endif
+
       libcloudphxx_const_derived(si::dimensionless, C_E, real_t(0.845));
       libcloudphxx_const_derived(si::dimensionless, C_tau, real_t(1.5));
       libcloudphxx_const_derived(si::dimensionless, cube_root_of_two_pi, pow(real_t(2) *
@@ -18,11 +25,6 @@ namespace libcloudphxx
         , real_t(1./3.))
       );
 
-#if !defined(__NVCC__)
-      using std::pow;
-      using std::sqrt;
-#endif
-
       typedef divide_typeof_helper<
         si::dimensionless,  
         si::length
@@ -62,7 +64,8 @@ namespace libcloudphxx
         const quantity<si::length, real_t> &L // characteristic length-scale
       )
       {
-        return quantity<area_over_time_squared, real_t>(pow((L * diss_rate) / si::cubic_metres * si::seconds * si::seconds * si::seconds / C_E<real_t>(), real_t(2./3.)) * si::metres * si::metres / si::seconds / si::seconds);
+        const real_t ret = cbrt( (L * diss_rate) / si::cubic_metres * si::seconds * si::seconds * si::seconds / C_E<real_t>());
+        return quantity<area_over_time_squared, real_t>(ret * ret * si::metres * si::metres / si::seconds / si::seconds);
       };
 
       template <typename real_t>

include/libcloudph++/common/SGS_length_scale.hpp

@@ -48,7 +48,7 @@ namespace libcloudphxx
         quantity<si::length, real_t> dy,
         quantity<si::length, real_t> dz
       )
-      { return quantity<si::length, real_t>(pow((dx/si::metres)*(dy/si::metres)*(dz/si::metres), real_t(1./3.)) * si::metres);}
+      { return quantity<si::length, real_t>(cbrt((dx/si::metres)*(dy/si::metres)*(dz/si::metres)) * si::metres);}
 
       // lambda = (dx+dy+dz) / 3
       template <typename real_t>

include/libcloudph++/common/const_cp.hpp

@@ -33,29 +33,29 @@ namespace libcloudphxx
       ) 
 //</listing-1>
       {
-	return p_tri<real_t>() * exp(
-	  (l_tri<real_t>() + (c_pw<real_t>() - c_pv<real_t>()) * T_tri<real_t>()) / R_v<real_t>() * (real_t(1) / T_tri<real_t>() - real_t(1) / T)
-	  - (c_pw<real_t>() - c_pv<real_t>()) / R_v<real_t>() * std::log(T / T_tri<real_t>())
-	);
+        return p_tri<real_t>() * exp(
+          (l_tri<real_t>() + (c_pw<real_t>() - c_pv<real_t>()) * T_tri<real_t>()) / R_v<real_t>() * (real_t(1) / T_tri<real_t>() - real_t(1) / T)
+          - (c_pw<real_t>() - c_pv<real_t>()) / R_v<real_t>() * std::log(T / T_tri<real_t>())
+        );
       }
 
       // saturation vapour mixing ratio for water as a function of pressure and temperature
       template <typename real_t>
       BOOST_GPU_ENABLED
       quantity<si::dimensionless, real_t> r_vs(
-	const quantity<si::temperature, real_t> &T,
-	const quantity<si::pressure, real_t> &p
+        const quantity<si::temperature, real_t> &T,
+        const quantity<si::pressure, real_t> &p
       ) {
-	return eps<real_t>() / (p / p_vs<real_t>(T) - 1);
+        return eps<real_t>() / (p / p_vs<real_t>(T) - 1);
       }
 
       // latent heat for constant c_p
       template <typename real_t>
       BOOST_GPU_ENABLED
       quantity<divide_typeof_helper<si::energy, si::mass>::type , real_t> l_v(
-	const quantity<si::temperature, real_t> &T
+        const quantity<si::temperature, real_t> &T
       ) {
-	return l_tri<real_t>() + (c_pv<real_t>() - c_pw<real_t>()) * (T - T_tri<real_t>());
+        return l_tri<real_t>() + (c_pv<real_t>() - c_pw<real_t>()) * (T - T_tri<real_t>());
       }
     };
   };