Merge pull request #399 from pdziekan/sgs_tke

add diags for up vp wp moms

Author: pdziekan

include/libcloudph++/lgrngn/particles.hpp

@@ -102,6 +102,9 @@ namespace libcloudphxx
       virtual void diag_chem(const enum common::chem::chem_species_t&)          { assert(false); }
       virtual void diag_precip_rate()                                           { assert(false); }
       virtual void diag_kappa_mom(const int&)                                   { assert(false); }
+      virtual void diag_up_mom(const int&)                                      { assert(false); }
+      virtual void diag_vp_mom(const int&)                                      { assert(false); }
+      virtual void diag_wp_mom(const int&)                                      { assert(false); }
       virtual void diag_incloud_time_mom(const int&)                            { assert(false); } // requires opts_init.diag_incloud_time==true
       virtual void diag_max_rw()                                                { assert(false); }
       virtual void diag_vel_div()                                               { assert(false); }
@@ -181,6 +184,9 @@ namespace libcloudphxx
       void diag_dry_mom(const int &k);
       void diag_wet_mom(const int &k);
       void diag_kappa_mom(const int &k);
+      void diag_up_mom(const int&);
+      void diag_vp_mom(const int&);
+      void diag_wp_mom(const int&);
       void diag_incloud_time_mom(const int &k);
       void diag_wet_mass_dens(const real_t&, const real_t&);
 
@@ -274,6 +280,9 @@ namespace libcloudphxx
       void diag_dry_mom(const int &k);
       void diag_wet_mom(const int &k);
       void diag_kappa_mom(const int&);
+      void diag_up_mom(const int&);
+      void diag_vp_mom(const int&);
+      void diag_wp_mom(const int&);
       void diag_incloud_time_mom(const int&);
       void diag_wet_mass_dens(const real_t&, const real_t&);
       real_t *outbuf();

src/impl/particles_impl_moms.ipp

@@ -155,14 +155,22 @@ namespace libcloudphxx
           const real_t n = thrust::get<0>(tpl);
           const real_t x = thrust::get<1>(tpl);
 #if !defined(NDEBUG)
-          real_t res = n * pow(x, xp); // TODO: check if xp=0 is optimised
+          real_t res;
+          if(x >= 0)
+            res = n * pow(x, xp);
+          else
+          {
+            int ixp = xp;
+            assert(ixp == xp);
+            res = n * pow(x, ixp);
+          }
           if(isnaninf()(res))
           {
             printf("nan/inf res in moment counter, n = %g x = %g res = %g xp = %g\n", n, x, res, xp);
           }
           return res;
 #else
-          return n * pow(x, xp); // TODO: check if xp=0 is optimised
+          return x >= 0 ? n * pow(x, xp) : n * pow(x, int(xp)); // for negative x and non-integer xp, pow = NaN
 #endif
         }
       };

src/particles_diag.ipp

@@ -318,6 +318,27 @@ namespace libcloudphxx
       pimpl->moms_calc(pimpl->kpa.begin(), n);
     }   
 
+    // compute n-th moment of SGS velocity in x for selected particles
+    template <typename real_t, backend_t device>
+    void particles_t<real_t, device>::diag_up_mom(const int &n)
+    {   
+      pimpl->moms_calc(pimpl->up.begin(), n);
+    }   
+
+    // compute n-th moment of SGS velocity in y for selected particles
+    template <typename real_t, backend_t device>
+    void particles_t<real_t, device>::diag_vp_mom(const int &n)
+    {   
+      pimpl->moms_calc(pimpl->vp.begin(), n);
+    }   
+
+    // compute n-th moment of SGS velocity in z for selected particles
+    template <typename real_t, backend_t device>
+    void particles_t<real_t, device>::diag_wp_mom(const int &n)
+    {   
+      pimpl->moms_calc(pimpl->wp.begin(), n);
+    }   
+
     // compute n-th moment of incloud_time for selected particles
     template <typename real_t, backend_t device>
     void particles_t<real_t, device>::diag_incloud_time_mom(const int &n)

src/particles_multi_gpu_diag.ipp

@@ -121,6 +121,24 @@ namespace libcloudphxx
       pimpl->mcuda_run(&particles_t<real_t, CUDA>::diag_kappa_mom, k);
     }   
 
+    template <typename real_t>
+    void particles_t<real_t, multi_CUDA>::diag_up_mom(const int &k)
+    {
+      pimpl->mcuda_run(&particles_t<real_t, CUDA>::diag_up_mom, k);
+    }   
+ 
+    template <typename real_t>
+    void particles_t<real_t, multi_CUDA>::diag_vp_mom(const int &k)
+    {
+      pimpl->mcuda_run(&particles_t<real_t, CUDA>::diag_vp_mom, k);
+    }   
+ 
+    template <typename real_t>
+    void particles_t<real_t, multi_CUDA>::diag_wp_mom(const int &k)
+    {
+      pimpl->mcuda_run(&particles_t<real_t, CUDA>::diag_wp_mom, k);
+    }   
+ 
     template <typename real_t>
     void particles_t<real_t, multi_CUDA>::diag_incloud_time_mom(const int &k)
     {