Ice

CMakeLists.txt

@@ -241,7 +241,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>)
@@ -259,6 +259,8 @@ if (CMAKE_CUDA_COMPILER)
 endif()
 
 ############################################################################################
+# find Thrust using libcloudphxx/cmake/obsolete/FindThrust.cmake
+
 # Thrust, location of Thrust can be hinted by setting THRUST_INCLUDE_DIR
 find_package(Thrust)
 
@@ -275,6 +277,13 @@ endif()
 # include thrust found here instead of the one shipped with cuda
 target_include_directories(cloudphxx_lgrngn PRIVATE ${THRUST_INCLUDE_DIR})
 
+## alternative based on libcloudphxx/cmake/alternative/thrust-config.cmake
+## TODO: doesn't work
+#SET(Thrust_DIR ${PROJECT_SOURCE_DIR}/cmake)
+#find_package(Thrust REQUIRED)
+#thrust_create_target(Thrust)
+#target_link_libraries(cloudphxx_lgrngn Thrust)
+
 ############################################################################################
 # Boost libraries
 find_package(Boost COMPONENTS) # we only need header-only boost libs

bindings/python/common.hpp

@@ -69,6 +69,12 @@ 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)

bindings/python/lgrngn.hpp

@@ -242,11 +242,13 @@ namespace libcloudphxx
         arg->dry_distros.clear();
         for (int i = 0; i < len(kappa_func.keys()); ++i)
           arg->dry_distros.emplace(
-            bp::extract<real_t>(kappa_func.keys()[i]), 
+            std::make_pair(real_t(bp::extract<real_t>(kappa_func.keys()[i])), 0),  // assume ice=0
             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,
@@ -255,10 +257,11 @@ namespace libcloudphxx
         arg->src_dry_distros.clear();
         for (int i = 0; i < len(kappa_func.keys()); ++i)
           arg->src_dry_distros.emplace(
-            bp::extract<real_t>(kappa_func.keys()[i]), 
+            std::make_pair(bp::extract<real_t>(kappa_func.keys()[i]), 0),  // assume ice=0
             std::make_shared<detail::pyunary<real_t>>(kappa_func.values()[i])
           );
       }
+      */
 
       template <typename real_t>
       void set_ds( // dry_sizes
@@ -283,13 +286,15 @@ 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[std::make_pair(kappa, 0)] = size_conc_map; // assume ice=0
         }
       }
 
+// 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,
@@ -319,6 +324,7 @@ namespace libcloudphxx
           arg->src_dry_sizes[kappa] = size_conc_map;
         }
       }
+      */
 
       template <typename real_t>
       void set_rdd( // rlx_dry_distros

bindings/python/lib.cpp

@@ -131,6 +131,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>);
@@ -264,8 +265,8 @@ BOOST_PYTHON_MODULE(libcloudphxx)
     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("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)

include/libcloudph++/common/const_cp.hpp

@@ -11,6 +11,7 @@ namespace libcloudphxx
     namespace const_cp
     {
       using moist_air::c_pw;
+      using moist_air::c_pi;
       using moist_air::c_pv;
       using moist_air::R_v;
       using moist_air::eps;
@@ -20,7 +21,8 @@ namespace libcloudphxx
       // water triple point parameters
       libcloudphxx_const(si::pressure, p_tri, 611.73, si::pascals) // pressure
       libcloudphxx_const(si::temperature, T_tri, 273.16, si::kelvins) // temperature
-      libcloudphxx_const(energy_over_mass, l_tri, 2.5e6, si::joules / si::kilograms) // latent heat of evaporation
+      libcloudphxx_const(energy_over_mass, lv_tri, 2.5e6, si::joules / si::kilograms)   // latent heat of evaporation
+      libcloudphxx_const(energy_over_mass, ls_tri, 2.834e6, si::joules / si::kilograms) // latent heat of sublimation
 
       // saturation vapour pressure for water assuming constant c_p_v and c_p_w
       // with constants taken at triple point
@@ -34,11 +36,26 @@ 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)
+          (lv_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 pressure for ice assuming constant c_p_v and c_p_i
+      // with constants taken at triple point
+      // (solution to the Clausius-Clapeyron equation assuming rho_vapour << rho_solid)
+      template <typename real_t>
+      BOOST_GPU_ENABLED
+      quantity<si::pressure, real_t> p_vsi(
+        const quantity<si::temperature, real_t> &T
+      ) 
+      {
+        return p_tri<real_t>() * exp(
+          (ls_tri<real_t>() + (c_pi<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_pi<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
@@ -49,13 +66,32 @@ namespace libcloudphxx
         return eps<real_t>() / (p / p_vs<real_t>(T) - 1);
       }
 
-      // latent heat for constant c_p
+      // saturation vapour mixing ratio for ice as a function of pressure and temperature
+      template <typename real_t>
+      BOOST_GPU_ENABLED
+      quantity<si::dimensionless, real_t> r_vsi(
+        const quantity<si::temperature, real_t> &T,
+        const quantity<si::pressure, real_t> &p
+      ) {
+        return eps<real_t>() / (p / p_vsi<real_t>(T) - 1);
+      }
+
+      // latent heat of evaporation 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
       ) {
-        return l_tri<real_t>() + (c_pv<real_t>() - c_pw<real_t>()) * (T - T_tri<real_t>());
+        return lv_tri<real_t>() + (c_pv<real_t>() - c_pw<real_t>()) * (T - T_tri<real_t>());
+      }
+
+      // latent heat sublimation for constant c_p
+      template <typename real_t>
+      BOOST_GPU_ENABLED
+      quantity<divide_typeof_helper<si::energy, si::mass>::type , real_t> l_s(
+        const quantity<si::temperature, real_t> &T
+      ) {
+        return ls_tri<real_t>() + (c_pv<real_t>() - c_pi<real_t>()) * (T - T_tri<real_t>());
       }
     };
   };

include/libcloudph++/common/detail/toms748.hpp

@@ -301,9 +301,9 @@ T toms748_solve(F f, const T& ax, const T& bx, const T& fax, const T& fbx, Tol t
    fb = fbx;
 
 #if !defined(NDEBUG)
-   if(a >= b)
+   if(a > b)
    {
-     printf("toms a >= b; a = %g b = %g fa = %g fb = %g\n", a, b, fa, fb);
+     printf("toms a > b; a = %g b = %g fa = %g fb = %g\n", a, b, fa, fb);
      assert(0);
    }
 #endif

include/libcloudph++/common/kappa_koehler.hpp

@@ -36,7 +36,7 @@ namespace libcloudphxx
       {   
         assert(RH > .03); // kappa-koehler assumes well dissolved matter
         assert(RH < 1.); // no equilibrium over RH=100%
-        assert(kappa > 0); // pure-water case left out
+        assert(kappa >= 0); 
         return rd3 * (1 - RH * (1 - kappa)) / (1 - RH);
       }   
 
@@ -49,7 +49,7 @@ namespace libcloudphxx
         quantity<si::dimensionless, real_t> kappa
       )
       {
-        assert(kappa > 0); // pure-water case left out
+        assert(kappa > 0);
         return (rw3 - rd3) / (rw3 - rd3 * (real_t(1) - kappa));
       }
 
@@ -134,7 +134,9 @@ namespace libcloudphxx
       )   
       {   
         assert(RH < 1); // no equilibrium over RH=100%
-        assert(kappa > 0); // pure-water case left out
+        assert(kappa >= 0);
+
+        if(kappa == 0) return rd3;
 
         return common::detail::toms748_solve(
           detail::rw3_eq_minfun<real_t>(RH, rd3, kappa, T), // the above-defined functor
@@ -154,7 +156,7 @@ namespace libcloudphxx
         quantity<si::temperature, real_t> T
       )   
       {   
-        assert(kappa > 0); // pure-water case left out
+        assert(kappa > 0);
 
         quantity<si::volume, double> rd3_dbl = static_cast<quantity<si::volume, double> >(rd3);
         quantity<si::temperature, double> T_dbl = static_cast<quantity<si::temperature, double> >(T);
@@ -175,7 +177,7 @@ namespace libcloudphxx
         quantity<si::temperature, real_t> T
       )
       {
-        assert(kappa > 0); // pure-water case left out
+        assert(kappa > 0);
 #if !defined(__NVCC__)
         using std::pow;
         using std::cbrt;

include/libcloudph++/common/maxwell-mason.hpp

@@ -15,34 +15,71 @@ namespace libcloudphxx
       quantity<divide_typeof_helper<si::area, si::time>::type, real_t> rdrdt(
         const quantity<diffusivity, real_t> D,            // D 
         const quantity<thermal_conductivity, real_t> K,   // K
-	const quantity<si::mass_density, real_t> rho_v,   // ambient water vapour density
-	const quantity<si::temperature, real_t> T,        // ambient temperature
-	const quantity<si::pressure, real_t> p,           // ambient pressure
-	const quantity<si::dimensionless, real_t> RH,     // p_v/p_vs = relative humidity
-	const quantity<si::dimensionless, real_t> a_w,    // water activity
-	const quantity<si::dimensionless, real_t> klvntrm // the Kelvin term
+        const quantity<si::mass_density, real_t> rho_v,   // ambient water vapour density
+        const quantity<si::temperature, real_t> T,        // ambient temperature
+        const quantity<si::pressure, real_t> p,           // ambient pressure
+        const quantity<si::dimensionless, real_t> RH,     // p_v/p_vs = relative humidity
+        const quantity<si::dimensionless, real_t> a_w,    // water activity
+        const quantity<si::dimensionless, real_t> klvntrm // the Kelvin term
       )
       {
         using moist_air::rho_w;
+        using moist_air::rho_i;
         using moist_air::R_v;
 
-	quantity<divide_typeof_helper<si::energy, si::mass>::type, real_t> 
+        quantity<divide_typeof_helper<si::energy, si::mass>::type, real_t> 
           l_v = const_cp::l_v<real_t>(T);
 
-	return (real_t(1) - a_w * klvntrm / RH)
-	  / rho_w<real_t>() 
-	  / ( 
-	    real_t(1) 
+        return (real_t(1) - a_w * klvntrm / RH)
+          / rho_w<real_t>() 
+          / ( 
+            real_t(1) 
               / D
               / rho_v 
-	    +
-	    l_v 
+            +
+            l_v 
               / K
               / RH  
               / T 
               * (l_v / R_v<real_t>() / T - real_t(1))
-	  )   
-	;   
+          )   
+        ;   
+      }
+
+      // for ice
+      // mass accommodation coeff = 1
+      // no solute nor curvature effects
+      template <typename real_t>
+      BOOST_GPU_ENABLED
+      quantity<divide_typeof_helper<si::area, si::time>::type, real_t> rdrdt_i(
+        const quantity<diffusivity, real_t> D,            // D 
+        const quantity<thermal_conductivity, real_t> K,   // K
+        const quantity<si::mass_density, real_t> rho_v,   // ambient water vapour density
+        const quantity<si::temperature, real_t> T,        // ambient temperature
+        const quantity<si::pressure, real_t> p,           // ambient pressure
+        const quantity<si::dimensionless, real_t> RH_i    // p_v/p_vsi = relative humidity w.r.t. ice
+      )
+      {
+        using moist_air::rho_i;
+        using moist_air::R_v;
+
+        quantity<divide_typeof_helper<si::energy, si::mass>::type, real_t> 
+          l_s = const_cp::l_s<real_t>(T);
+
+        return (real_t(1) - real_t(1) / RH_i)
+          / rho_i<real_t>() 
+          / ( 
+            real_t(1) 
+              / D
+              / rho_v 
+            +
+            l_s 
+              / K
+              / RH_i
+              / T 
+              * (l_s / R_v<real_t>() / T - real_t(1))
+          )   
+        ;   
       }
     };
   };

include/libcloudph++/common/moist_air.hpp

@@ -26,6 +26,7 @@ namespace libcloudphxx
       libcloudphxx_const(energy_over_temperature_mass, c_pd, 1005, si::joules / si::kilograms / si::kelvins) // dry air
       libcloudphxx_const(energy_over_temperature_mass, c_pv, 1850, si::joules / si::kilograms / si::kelvins) // water vapour
       libcloudphxx_const(energy_over_temperature_mass, c_pw, 4218, si::joules / si::kilograms / si::kelvins) // liquid water
+      libcloudphxx_const(energy_over_temperature_mass, c_pi, 2114, si::joules / si::kilograms / si::kelvins) // ice
 
       // molar masses
       libcloudphxx_const(mass_over_amount, M_d, 0.02897, si::kilograms / si::moles) // dry air (Curry & Webster / Seinfeld & Pandis)
@@ -47,6 +48,8 @@ namespace libcloudphxx
 
       // water density
       libcloudphxx_const(si::mass_density, rho_w, 1e3, si::kilograms / si::cubic_metres)
+      // ice density
+      libcloudphxx_const(si::mass_density, rho_i, 910, si::kilograms / si::cubic_metres)
 
       // mixing rule for extensive quantitites (i.e. using mass mixing ratio)
       template <typename real_t, typename quant>

include/libcloudph++/common/output.hpp

@@ -17,7 +17,9 @@ namespace libcloudphxx
       outH = chem::H,                          
       outliq_vol,
       outdry_vol,
-      outprtcl_num
+      outice_vol,
+      outliq_num,
+      outice_num
     };
 
     const std::map<output_t, std::string> output_names
@@ -30,9 +32,11 @@ namespace libcloudphxx
       {outO3     , "O3"},
       {outS_VI   , "S_VI"},                        
       {outH      , "H"},                          
-      {outliq_vol, "liquid_volume"},
+      {outliq_vol, "water_volume"},
       {outdry_vol, "dry_volume"},
-      {outprtcl_num, "particle_number"}
+      {outice_vol, "ice_volume"},
+      {outliq_num, "water_number"},
+      {outice_num, "ice_number"}
     };
   };
 };