Ice in super-droplets

CMakeLists.txt

@@ -273,7 +273,7 @@ if (CMAKE_CUDA_COMPILER)
     target_compile_options(cloudphxx_lgrngn PRIVATE  $<$<COMPILE_LANGUAGE:CUDA>: -O3 -use_fast_math -Xcompiler=-Ofast,-march=native,-fopenmp>)
   # Debug mode cuda flags
   elseif(CMAKE_BUILD_TYPE STREQUAL "Debug")
-    target_compile_options(cloudphxx_lgrngn PRIVATE  $<$<COMPILE_LANGUAGE:CUDA>: -g -DTHRUST_DEBUG -lineinfo -Xcompiler=-fopenmp,-g,-Og,-rdynamic>)
+    target_compile_options(cloudphxx_lgrngn PRIVATE  $<$<COMPILE_LANGUAGE:CUDA>: -g -DTHRUST_DEBUG -lineinfo -Xcompiler=-fopenmp,-g,-Og,-rdynamic,-Wfatal-errors>)
   # Release mode cuda flags
   else()
     target_compile_options(cloudphxx_lgrngn PRIVATE  $<$<COMPILE_LANGUAGE:CUDA>: -DNDEBUG -O3 -use_fast_math -Xcompiler=-Ofast,-march=native,-DNDEBUG,-fopenmp>)
@@ -396,7 +396,7 @@ add_subdirectory(include)
 
 if(NOT LIBCLOUDPHXX_DISABLE_BINDINGS) # some tests require bindings, so we disable them too
   enable_testing()
-  add_subdirectory(bindings) 
+  add_subdirectory(bindings)
   add_subdirectory(tests)
 endif()
 

bindings/python/common.hpp

@@ -58,6 +58,12 @@ namespace libcloudphxx
 	return cmn::const_cp::p_vs(T * si::kelvins) / si::pascals;
       }
 
+      template <typename real_t>
+      real_t p_vsi(const real_t &T)
+      {
+        return cmn::const_cp::p_vsi(T * si::kelvins) / si::pascals;
+      }
+
       template <typename real_t>
       real_t p_vs_tet(const real_t &T)
       {
@@ -69,13 +75,31 @@ namespace libcloudphxx
       {
 	return cmn::const_cp::r_vs(T * si::kelvins, p * si::pascals);
       }
+
+      template <typename real_t>
+      real_t r_vsi(const real_t &T, const real_t &p)
+      {
+	return cmn::const_cp::r_vsi(T * si::kelvins, p * si::pascals);
+      }
 
       template <typename real_t>
       real_t l_v(const real_t &T)
       {
 	return cmn::const_cp::l_v(T * si::kelvins) / si::joules * si::kilograms;
       }
 
+      template <typename real_t>
+      real_t l_s(const real_t &T)
+      {
+        return cmn::const_cp::l_s(T * si::kelvins) / si::joules * si::kilograms;
+      }
+
+      template <typename real_t>
+      real_t l_f(const real_t &T)
+      {
+        return cmn::const_cp::l_f(T * si::kelvins) / si::joules * si::kilograms;
+      }
+
       template <typename real_t>
       real_t T(const real_t &th, const real_t &rhod)
       {

bindings/python/lgrngn.hpp

@@ -222,7 +222,7 @@ namespace libcloudphxx
         map_t map = arg->diag_puddle();
         bp::dict dict;
         for(auto& x : map)
-          dict[static_cast<int>(x.first)] = x.second;
+          dict[cmn::output_names.at(x.first)] = x.second;
         return dict;     
       }
 
@@ -234,45 +234,71 @@ namespace libcloudphxx
         return *arg->opts_init;
       }
 
+      // set dry distros from a dict with (kappa, rd_insol) as key
       template <typename real_t>
-      void set_dd( // dry_distro
+      void set_dd(
         lgr::opts_init_t<real_t> *arg,
-        const bp::dict &kappa_func)
+        const bp::dict &kappa_func) // a dict keyed by (kappa, rd_insol)
       {
         arg->dry_distros.clear();
         for (int i = 0; i < len(kappa_func.keys()); ++i)
+        {
+          bp::tuple key = bp::extract<bp::tuple>(kappa_func.keys()[i]);
+          const real_t kappa = bp::extract<real_t>(key[0]);
+          const real_t rd_insol = bp::extract<real_t>(key[1]);
           arg->dry_distros.emplace(
-            bp::extract<real_t>(kappa_func.keys()[i]), 
+            libcloudphxx::lgrngn::kappa_rd_insol_t<real_t>{kappa, rd_insol},
             std::make_shared<detail::pyunary<real_t>>(kappa_func.values()[i])
           );
+        }
       }
 
+      // src_dry_distros moved from opts_init to opts
       template <typename real_t>
       void set_sdd( // src_dry_distro
-        lgr::opts_init_t<real_t> *arg,
-        const bp::dict &kappa_func)
+        lgr::opts_t<real_t> *arg,
+        const bp::dict &kappa_func) // a dict keyed by (kappa, rd_insol)
       {
         arg->src_dry_distros.clear();
         for (int i = 0; i < len(kappa_func.keys()); ++i)
+        {
+          bp::tuple key = bp::extract<bp::tuple>(kappa_func.keys()[i]);
+          bp::tuple val = bp::extract<bp::tuple>(kappa_func.values()[i]);
+          const real_t kappa = bp::extract<real_t>(key[0]);
+          const real_t rd_insol = bp::extract<real_t>(key[1]);
+
+          const int sd_conc   = bp::extract<int>(val[1]);
+          const int supstp    = bp::extract<int>(val[2]);
           arg->src_dry_distros.emplace(
-            bp::extract<real_t>(kappa_func.keys()[i]), 
-            std::make_shared<detail::pyunary<real_t>>(kappa_func.values()[i])
+            libcloudphxx::lgrngn::kappa_rd_insol_t<real_t>{kappa, rd_insol},
+            std::make_tuple(std::static_pointer_cast<libcloudphxx::common::unary_function<real_t>>(
+          std::make_shared<detail::pyunary<real_t>>(val[0])),
+              sd_conc, supstp)
           );
+        }
       }
 
+      // set dry sizes from a dict with (kappa, rd_insol) as key
       template <typename real_t>
-      void set_ds( // dry_sizes
+      void set_ds(
         lgr::opts_init_t<real_t> *arg,
-        const bp::dict &kappa_func
+        const bp::dict &kappa_func // a dict keyed by (kappa, rd_insol)
       )
       {
         arg->dry_sizes.clear();
         if(len(kappa_func.keys()) == 0)
           return;
 
-        // loop over kappas
+        // loop over kappas and rd_insol
         for (int j = 0; j < len(kappa_func.keys()); ++j)
         {
+
+          // extract the key tuple (kappa, rd_insol)
+          const bp::tuple key = bp::extract<bp::tuple>(kappa_func.keys()[j]);
+          const real_t kappa    = bp::extract<real_t>(key[0]);
+          const real_t rd_insol = bp::extract<real_t>(key[1]);
+
+          // extract size : {conc, count} dict for this (kappa, rd_insol)
           const bp::dict size_conc = bp::extract<bp::dict>(kappa_func.values()[j]);
           std::map<real_t, std::pair<real_t, int>> size_conc_map;
 
@@ -283,43 +309,50 @@ namespace libcloudphxx
             assert(len(conc_count_list) == 2);
             const real_t conc = bp::extract<real_t>(conc_count_list[0]);
             const int count   = bp::extract<int>   (conc_count_list[1]);
-            size_conc_map[bp::extract<real_t>(size_conc.keys()[i])] = std::make_pair(conc, count);
+            size_conc_map[bp::extract<real_t>(size_conc.keys()[i])] = std::make_pair(conc, count); 
           }
-          const real_t kappa = bp::extract<real_t>(kappa_func.keys()[j]);
-          arg->dry_sizes[kappa] = size_conc_map;
+          arg->dry_sizes[libcloudphxx::lgrngn::kappa_rd_insol_t<real_t>{kappa, rd_insol}] = size_conc_map;
         }
       }
 
+      // src_dry_sizes moved from opts_init to opts
       template <typename real_t>
       void set_sds( // src_dry_sizes
-        lgr::opts_init_t<real_t> *arg,
-        const bp::dict &kappa_func
+        lgr::opts_t<real_t> *arg,
+        const bp::dict &kappa_func // a dict keyed by (kappa, rd_insol)
       )
       {
         arg->src_dry_sizes.clear();
         if(len(kappa_func.keys()) == 0)
           return;
 
-        // loop over kappas
+        // loop over kappas and rd_insol
         for (int j = 0; j < len(kappa_func.keys()); ++j)
         {
+          // extract the key tuple (kappa, rd_insol)
+          const bp::tuple key = bp::extract<bp::tuple>(kappa_func.keys()[j]);
+          const real_t kappa    = bp::extract<real_t>(key[0]);
+          const real_t rd_insol = bp::extract<real_t>(key[1]);
+
+          // extract size : {conc, count, supstp_src} dict for this (kappa, rd_insol)
           const bp::dict size_conc = bp::extract<bp::dict>(kappa_func.values()[j]);
-          std::map<real_t, std::pair<real_t, int>> size_conc_map;
+          std::map<double, std::tuple<real_t, int, int>> size_conc_map;
 
-          // turn the size : {conc, count} dict into a size : {conc, count} map
+          // turn the size : {conc, count, supstp_src} dict into a size : {conc, count, supstp_src} map
           for (int i = 0; i < len(size_conc.keys()); ++i)
           {
             const bp::list conc_count_list = bp::extract<bp::list>(size_conc.values()[i]);
-            assert(len(conc_count_list) == 2);
+            assert(len(conc_count_list) == 3);
             const real_t conc = bp::extract<real_t>(conc_count_list[0]);
             const int count   = bp::extract<int>   (conc_count_list[1]);
-            size_conc_map[bp::extract<real_t>(size_conc.keys()[i])] = std::make_pair(conc, count);
+            const int supstp  = bp::extract<int>   (conc_count_list[2]);
+            size_conc_map[bp::extract<real_t>(size_conc.keys()[i])] = std::make_tuple(conc, count, supstp);
           }
-          const real_t kappa = bp::extract<real_t>(kappa_func.keys()[j]);
-          arg->src_dry_sizes[kappa] = size_conc_map;
+          arg->src_dry_sizes[libcloudphxx::lgrngn::kappa_rd_insol_t<real_t>{kappa, rd_insol}] = size_conc_map;
         }
       }
 
+
       template <typename real_t>
       void set_rdd( // rlx_dry_distros
         lgr::opts_init_t<real_t> *arg,

bindings/python/lib.cpp

@@ -63,6 +63,7 @@ BOOST_PYTHON_MODULE(libcloudphxx)
     bp::scope().attr("eps") = (real_t) (cmn::moist_air::eps<real_t>());
     bp::scope().attr("rho_stp") = (real_t) (cmn::earth::rho_stp<real_t>() * si::cubic_metres / si::kilograms);
     bp::scope().attr("rho_w") = (real_t) (cmn::moist_air::rho_w<real_t>() * si::cubic_metres / si::kilograms);
+    bp::scope().attr("rho_i") = (real_t) (cmn::moist_air::rho_i<real_t>() * si::cubic_metres / si::kilograms);
     //molar mass of trace gases
     bp::scope().attr("M_SO2")   = (real_t) (cmn::molar_mass::M_SO2<real_t>()   * si::moles / si::kilograms);
     bp::scope().attr("M_H2O2")  = (real_t) (cmn::molar_mass::M_H2O2<real_t>()  * si::moles / si::kilograms);
@@ -131,6 +132,7 @@ BOOST_PYTHON_MODULE(libcloudphxx)
     bp::def("p_v", &common::p_v<real_t>);
     bp::def("p_vs", &common::p_vs<real_t>);
     bp::def("r_vs", &common::r_vs<real_t>);
+    bp::def("r_vsi", &common::r_vsi<real_t>);
     bp::def("p_vs_tet", &common::p_vs_tet<real_t>);
     bp::def("l_v", &common::l_v<real_t>);
     bp::def("T", &common::T<real_t>);
@@ -278,13 +280,14 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .def_readwrite("turb_adve", &lgr::opts_t<real_t>::turb_adve)
       .def_readwrite("turb_cond", &lgr::opts_t<real_t>::turb_cond)
       .def_readwrite("turb_coal", &lgr::opts_t<real_t>::turb_coal)
+      .def_readwrite("ice_nucl", &lgr::opts_t<real_t>::ice_nucl)
       .def_readwrite("dt", &lgr::opts_t<real_t>::dt)
+      .add_property("src_dry_distros", &lgrngn::get_sdd<real_t>, &lgrngn::set_sdd<real_t>)
+      .add_property("src_dry_sizes", &lgrngn::get_ds<real_t>, &lgrngn::set_sds<real_t>)
     ;
     bp::class_<lgr::opts_init_t<real_t>>("opts_init_t")
       .add_property("dry_distros", &lgrngn::get_dd<real_t>, &lgrngn::set_dd<real_t>)
       .add_property("dry_sizes", &lgrngn::get_ds<real_t>, &lgrngn::set_ds<real_t>)
-      .add_property("src_dry_distros", &lgrngn::get_sdd<real_t>, &lgrngn::set_sdd<real_t>)
-      .add_property("src_dry_sizes", &lgrngn::get_ds<real_t>, &lgrngn::set_sds<real_t>)
       .add_property("rlx_dry_distros", &lgrngn::get_rdd<real_t>, &lgrngn::set_rdd<real_t>)
       .def_readwrite("nx", &lgr::opts_init_t<real_t>::nx)
       .def_readwrite("ny", &lgr::opts_init_t<real_t>::ny)
@@ -321,15 +324,13 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .def_readwrite("sstp_cond", &lgr::opts_init_t<real_t>::sstp_cond)
       .def_readwrite("sstp_coal", &lgr::opts_init_t<real_t>::sstp_coal)
       .def_readwrite("sstp_chem", &lgr::opts_init_t<real_t>::sstp_chem)
-      .def_readwrite("supstp_src", &lgr::opts_init_t<real_t>::supstp_src)
       .def_readwrite("supstp_rlx", &lgr::opts_init_t<real_t>::supstp_rlx)
       .def_readwrite("kernel", &lgr::opts_init_t<real_t>::kernel)
       .def_readwrite("adve_scheme", &lgr::opts_init_t<real_t>::adve_scheme)
       .def_readwrite("sd_conc", &lgr::opts_init_t<real_t>::sd_conc)
       .def_readwrite("sd_conc_large_tail", &lgr::opts_init_t<real_t>::sd_conc_large_tail)
       .def_readwrite("aerosol_independent_of_rhod", &lgr::opts_init_t<real_t>::aerosol_independent_of_rhod)
       .def_readwrite("sd_const_multi", &lgr::opts_init_t<real_t>::sd_const_multi)
-      .def_readwrite("src_sd_conc", &lgr::opts_init_t<real_t>::src_sd_conc)
       .def_readwrite("rlx_bins", &lgr::opts_init_t<real_t>::rlx_bins)
       .def_readwrite("rlx_sd_per_bin", &lgr::opts_init_t<real_t>::rlx_sd_per_bin)
       .def_readwrite("rlx_timescale", &lgr::opts_init_t<real_t>::rlx_timescale)
@@ -347,6 +348,8 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .add_property("SGS_mix_len", &lgrngn::get_SGS_mix_len<real_t>, &lgrngn::set_SGS_mix_len<real_t>)
       .add_property("kernel_parameters", &lgrngn::get_kp<real_t>, &lgrngn::set_kp<real_t>)
       .def_readwrite("variable_dt_switch", &lgr::opts_init_t<real_t>::variable_dt_switch)
+      .def_readwrite("ice_switch", &lgr::opts_init_t<real_t>::ice_switch)
+      .def_readwrite("time_dep_ice_nucl", &lgr::opts_init_t<real_t>::time_dep_ice_nucl)
       .def_readwrite("const_p", &lgr::opts_init_t<real_t>::const_p)
       .def_readwrite("th_dry", &lgr::opts_init_t<real_t>::th_dry)
     ;
@@ -409,6 +412,13 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .def("diag_chem",    &lgr::particles_proto_t<real_t>::diag_chem)
       .def("diag_precip_rate",    &lgr::particles_proto_t<real_t>::diag_precip_rate)
       .def("diag_puddle",    &lgrngn::diag_puddle<real_t>)
+      .def("diag_ice",    &lgr::particles_proto_t<real_t>::diag_ice)
+      .def("diag_water",    &lgr::particles_proto_t<real_t>::diag_water)
+      .def("diag_ice_cons",    &lgr::particles_proto_t<real_t>::diag_ice_cons)
+      .def("diag_water_cons",    &lgr::particles_proto_t<real_t>::diag_water_cons)
+      .def("diag_ice_a_mom",    &lgr::particles_proto_t<real_t>::diag_ice_a_mom)
+      .def("diag_ice_c_mom",    &lgr::particles_proto_t<real_t>::diag_ice_c_mom)
+      .def("diag_ice_mass",    &lgr::particles_proto_t<real_t>::diag_ice_mass)
       .def("outbuf",       &lgrngn::outbuf<real_t>)
       .def("get_attr",    &lgr::particles_proto_t<real_t>::get_attr)
     ;