Rc adjustment (needed if RH is advected)

bindings/python/lgrngn.hpp

@@ -122,6 +122,18 @@ namespace libcloudphxx
 	);
       }
 
+      // 
+      template <typename real_t>
+      void step_rc_adjust(
+      	lgr::particles_proto_t<real_t> *arg,
+	      const bp::numeric::array &rc_adjust
+      )
+      {
+      	arg->step_rc_adjust(
+      	  np2ai<real_t>(rc_adjust, sz(*arg))
+      	);
+      }
+
       template <typename real_t>
       const lgr::opts_init_t<real_t> get_oi(
         lgr::particles_proto_t<real_t> *arg

bindings/python/lib.cpp

@@ -278,6 +278,9 @@ BOOST_PYTHON_MODULE(libcloudphxx)
         bp::arg("Cz")  = bp::numeric::array(bp::object()),
         bp::arg("ambient_chem") = bp::dict()
       ))
+      .def("step_rc_adjust",    &lgrngn::step_rc_adjust<real_t>, (
+        bp::arg("rc_adjust")  = bp::numeric::array(bp::object())
+      ))
       .def("step_async",   &lgr::particles_proto_t<real_t>::step_async)
       .def("diag_sd_conc", &lgr::particles_proto_t<real_t>::diag_sd_conc)
       .def("diag_all",     &lgr::particles_proto_t<real_t>::diag_all)

include/libcloudph++/lgrngn/particles.hpp

@@ -50,6 +50,13 @@ namespace libcloudphxx
         return 0;
       }  
 
+      // add some cloud water
+      virtual void step_rc_adjust(
+        const arrinfo_t<real_t> rc_adj
+      ) { 
+        assert(false); 
+      }  
+
       // method for accessing super-droplet statistics
       virtual void diag_sd_conc()                                   { assert(false); }
       virtual void diag_all()                                       { assert(false); }
@@ -99,6 +106,10 @@ namespace libcloudphxx
         const opts_t<real_t> &
       );
 
+      void step_rc_adjust(
+        const arrinfo_t<real_t> 
+      ); 
+
       // diagnostic methods
       void diag_sd_conc();
       void diag_dry_rng(

src/particles.tpp

@@ -64,6 +64,7 @@
 #include "particles_impl_fill_outbuf.ipp"
 
 #include "particles_impl_sync.ipp"
+#include "particles_impl_rc_adjust.ipp"
 
 #include "particles_impl_adve.ipp"
 #include "particles_impl_cond.ipp"

src/particles_impl_rc_adjust.ipp

@@ -0,0 +1,101 @@
+// vim:filetype=cpp
+/** @file
+  * @copyright University of Warsaw
+  * @section LICENSE
+  * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/)
+  */
+
+namespace libcloudphxx
+{
+  namespace lgrngn
+  {
+    namespace detail
+    {
+      struct adj_rw2
+      {
+        template<class real_t>
+        BOOST_GPU_ENABLED
+        real_t operator()(real_t rw2, real_t m_added)
+        {
+          // turn into added volume
+          return pow(pow(rw2, real_t(3./2)) + 3. / 4. / 
+#if !defined(__NVCC__)
+            pi<real_t>()
+#else
+            CUDART_PI
+#endif
+            * m_added / (libcloudphxx::common::moist_air::rho_w<real_t>() / si::kilograms * si::cubic_metres),
+            real_t(2./3));
+        }
+      };
+      struct divide
+      {
+        template<class real_t, class n_t>
+        BOOST_GPU_ENABLED
+        real_t operator()(real_t x, n_t y)
+        {
+          return x / real_t(y);
+        }
+      };
+    };
+    template <typename real_t, backend_t device>
+    void particles_t<real_t, device>::impl::rc_adjust(
+      const arrinfo_t<real_t> &_rc_adjust
+    ) 
+    {   
+      // use tmp vector to store adjustment
+      thrust_device::vector<real_t> &rc_adjust(tmp_device_real_cell);
+      if(!l2e.count(&rc_adjust))
+        init_e2l(_rc_adjust, &rc_adjust);
+      // save input adjustment to the temp vector
+      sync(_rc_adjust, rc_adjust);
+
+      hskpng_sort();
+
+      // count number of particles per cell
+      thrust::pair<
+        thrust_device::vector<thrust_size_t>::iterator,
+        thrust_device::vector<n_t>::iterator
+      > np = thrust::reduce_by_key(
+        sorted_ijk.begin(), sorted_ijk.end(),   // input - keys
+        n.begin(),                     // input - values
+        count_ijk.begin(),                      // output - keys
+        count_num.begin()                       // output - values
+      );  
+      count_n = np.first - count_ijk.begin();
+
+      // divide water adjustment by number of drops
+      thrust::transform(
+        thrust::make_permutation_iterator(rc_adjust.begin(), count_ijk.begin()),
+        thrust::make_permutation_iterator(rc_adjust.begin(), count_ijk.end()),
+        thrust::make_permutation_iterator(count_num.begin(), count_ijk.begin()),
+        thrust::make_permutation_iterator(rc_adjust.begin(), count_ijk.begin()),
+        detail::divide()
+      );
+
+      // adjustment is in mass_of_added_water[kg] / mass_of_dry_air[kg]
+      // multiply it now by mass of dry air in the cell
+      thrust::transform(
+        rc_adjust.begin(), rc_adjust.end(),
+        dv.begin(),
+        rc_adjust.begin(),
+        thrust::multiplies<real_t>()
+      );
+      thrust::transform(
+        rc_adjust.begin(), rc_adjust.end(),
+        rhod.begin(),
+        rc_adjust.begin(),
+        thrust::multiplies<real_t>()
+      );
+
+      // apply the adjustment - change rw of drops
+      thrust::transform(
+        rw2.begin(),
+        rw2.end(),
+        thrust::make_permutation_iterator(rc_adjust.begin(), ijk.begin()), // mass to be added to each drop
+        rw2.begin(),
+        detail::adj_rw2()
+      );
+    }
+  };  
+};

src/particles_pimpl_ctor.ipp

@@ -322,6 +322,7 @@ namespace libcloudphxx
       void init_sstp();
       void init_kernel();
       void init_vterm();
+      void rc_adjust(const arrinfo_t<real_t>&);
 
       void fill_outbuf();
 

src/particles_step.ipp

@@ -222,5 +222,12 @@ namespace libcloudphxx
 
       return ret;
     }
+
+    template <typename real_t, backend_t device>
+    void particles_t<real_t, device>::step_rc_adjust(const arrinfo_t<real_t> rc_adjust)
+    {
+      // apply the adjustment
+      pimpl->rc_adjust(rc_adjust);
+    }
   };
 };

tests/python/unit/CMakeLists.txt

@@ -1,5 +1,5 @@
 # non-pytest tests
-foreach(test api_blk_1m api_blk_2m api_lgrngn api_common segfault_20150216 test_col_kernels test_terminal_velocities test_SD_removal test_uniform_init test_source test_chem_coal)
+foreach(test api_blk_1m api_blk_2m api_lgrngn api_common segfault_20150216 test_col_kernels test_terminal_velocities test_SD_removal test_uniform_init test_source test_chem_coal test_rc_adjustment)
   #TODO: indicate that tests depend on the lib
   add_test(
     NAME ${test}

tests/python/unit/test_rc_adjustment.py

@@ -0,0 +1,69 @@
+#test if after initialization we have approximately the same water content in each cell
+
+import sys 
+sys.path.insert(0, "../../bindings/python/")
+
+from libcloudphxx import lgrngn
+import numpy as np
+import time
+from math import pi
+
+# initial exponential distribution in droplet volume
+# as a function of ln(r)
+def expvolumelnr(lnr):
+  r=np.exp(lnr)
+  return n_zero * 3.*np.power(r,3)/np.power(r_zero,3)*np.exp(- np.power((r/r_zero),3));
+
+#initial conditions, ca. 2g / m^3
+r_zero = 15e-6
+n_zero = 1.42e8
+
+opts_init = lgrngn.opts_init_t()
+opts_init.coal_switch=0
+opts_init.sedi_switch=0
+opts_init.cond_switch=0
+opts_init.adve_switch=0
+opts_init.dt = 1
+opts_init.dx = 1
+opts_init.dz = 3
+opts_init.dy = 2
+opts_init.nx = 2 
+opts_init.nz = 3 
+opts_init.ny = 4 
+opts_init.x1 = opts_init.dx * opts_init.nx
+opts_init.z1 = opts_init.dz * opts_init.nz
+opts_init.y1 = opts_init.dy * opts_init.ny
+opts_init.rng_seed = int(time.time())
+
+th   = 300 * np.ones((opts_init.nx, opts_init.ny, opts_init.nz))  
+rv   = 0.01 * np.ones((opts_init.nx, opts_init.ny, opts_init.nz))  
+rc_adj   = 0.001 * np.ones((opts_init.nx, opts_init.ny, opts_init.nz))  
+rhod = 1. * np.ones((opts_init.nx, opts_init.ny, opts_init.nz)) + .1 * np.mgrid[1:1+opts_init.nx, 1:1+opts_init.ny, 1:1+opts_init.nz][1] # different densities, hence different water content
+
+kappa = 1e-6
+
+opts_init.dry_distros = {kappa:expvolumelnr}
+
+opts_init.sd_conc = 64
+opts_init.n_sd_max = opts_init.sd_conc * opts_init.nx * opts_init.ny * opts_init.nz
+
+try:
+  prtcls = lgrngn.factory(lgrngn.backend_t.OpenMP, opts_init)
+except:
+  prtcls = lgrngn.factory(lgrngn.backend_t.serial, opts_init)
+
+prtcls.init(th, rv, rhod)
+
+prtcls.diag_all()
+prtcls.diag_wet_mom(3) # gives specific moment (divided by rhod)
+water_pre =1e3* 4/3*pi* np.frombuffer(prtcls.outbuf())#.mean() # dropping a constant
+print water_pre
+
+prtcls.step_rc_adjust(rc_adj)
+
+prtcls.diag_wet_mom(3) # gives specific moment (divided by rhod)
+water_post = 1e3*4/3*pi*np.frombuffer(prtcls.outbuf())#.mean() # dropping a constant
+water_post -= rc_adj[0][0][0] * np.ones((opts_init.nx * opts_init.ny * opts_init.nz)) # 
+print water_post
+if(not np.allclose(water_pre, water_post,rtol=0,atol=1e-8)):
+  raise Exception("Error in adding rc adjustment")