Chemistry

CMakeLists.txt

@@ -55,7 +55,7 @@ target_include_directories(uwlcm_includes
 )
 
 # std::future test
-SET(CMAKE_REQUIRED_FLAGS "-std=c++11") # c++11 needed by the std_future test
+SET(CMAKE_REQUIRED_FLAGS "-std=c++17") 
 include(CheckCXXSourceCompiles)
 check_cxx_source_compiles("#include <future>\n int main() { std::future<void> f; }\n" STD_FUTURE_WORKS)
 if(${STD_FUTURE_WORKS})
@@ -99,6 +99,12 @@ if(UWLCM_DISABLE)
   if(NOT "3D_LGRNGN" IN_LIST UWLCM_DISABLE)
     target_sources(uwlcm PRIVATE src/run_hlpr_3d_lgrngn.cpp)
   endif()
+  if(NOT "2D_LGRNGN_CHEM" IN_LIST UWLCM_DISABLE)
+    target_sources(uwlcm PRIVATE src/run_hlpr_2d_lgrngn_chem.cpp)
+  endif()
+  if(NOT "3D_LGRNGN_CHEM" IN_LIST UWLCM_DISABLE)
+    target_sources(uwlcm PRIVATE src/run_hlpr_3d_lgrngn_chem.cpp)
+  endif()
   if(NOT "2D_NONE" IN_LIST UWLCM_DISABLE)
     target_sources(uwlcm PRIVATE src/run_hlpr_2d_none.cpp)
   endif()
@@ -107,7 +113,7 @@ if(UWLCM_DISABLE)
   endif()
 else()
   # if no options were disabled, add all sources
-  target_sources(uwlcm PRIVATE src/run_hlpr_2d_blk_1m.cpp src/run_hlpr_3d_blk_1m.cpp src/run_hlpr_2d_blk_2m.cpp src/run_hlpr_3d_blk_2m.cpp src/run_hlpr_2d_lgrngn.cpp src/run_hlpr_3d_lgrngn.cpp src/run_hlpr_2d_none.cpp src/run_hlpr_3d_none.cpp)
+  target_sources(uwlcm PRIVATE src/run_hlpr_2d_blk_1m.cpp src/run_hlpr_3d_blk_1m.cpp src/run_hlpr_2d_blk_2m.cpp src/run_hlpr_3d_blk_2m.cpp src/run_hlpr_2d_lgrngn.cpp src/run_hlpr_3d_lgrngn.cpp src/run_hlpr_2d_none.cpp src/run_hlpr_3d_none.cpp src/run_hlpr_2d_lgrngn_chem.cpp src/run_hlpr_3d_lgrngn_chem.cpp )
 endif()
 
 # ensure that uwlcm is not built before current git_revision file is created
@@ -130,7 +136,7 @@ target_link_libraries(uwlcm
 
 # enabling c++14, but not gnu++14
 set_target_properties(uwlcm PROPERTIES CXX_EXTENSIONS OFF)
-target_compile_features(uwlcm PRIVATE cxx_std_14)
+target_compile_features(uwlcm PRIVATE cxx_std_17)
 
 # search for Boost
 find_package(Boost COMPONENTS thread iostreams system timer program_options filesystem REQUIRED)

src/cases/DYCOMS.hpp

@@ -1,6 +1,10 @@
 #pragma once
 #include <random>
+#include <type_traits>
 #include "Anelastic.hpp"
+#include "../detail/blitz_hlpr_fctrs.hpp"
+#include <libcloudph++/common/chem.hpp>
+
 
 namespace cases 
 {
@@ -58,6 +62,31 @@ namespace cases
         rt_above = RF == 1 ? 1.5e-3 : (5. - 3. * (1. - exp((z_i[RF - 1] - z)/500.))) * 1e-3;
       return z < z_i[RF - 1] ? rt_below : rt_above;
     }
+
+    inline quantity<si::dimensionless, real_t> SO2g_dycoms(const real_t &z)
+    {
+      return quantity<si::dimensionless, real_t>(.2e-9);
+    }
+    inline quantity<si::dimensionless, real_t> CO2g_dycoms(const real_t &z)
+    {
+      return quantity<si::dimensionless, real_t>(360e-6);
+    }
+    inline quantity<si::dimensionless, real_t> O3g_dycoms(const real_t &z)
+    {
+      return quantity<si::dimensionless, real_t>(25e-9);
+    }
+    inline quantity<si::dimensionless, real_t> H2O2g_dycoms(const real_t &z)
+    {
+      return quantity<si::dimensionless, real_t>(.4e-9);
+    }
+    inline quantity<si::dimensionless, real_t> NH3g_dycoms(const real_t &z)
+    {
+      return quantity<si::dimensionless, real_t>(.1e-9);
+    }
+    inline quantity<si::dimensionless, real_t> HNO3g_dycoms(const real_t &z)
+    {
+      return quantity<si::dimensionless, real_t>(.1e-9);
+    }
 
     template<class case_ct_params_t, int RF, int n_dims>
     class DycomsCommon : public Anelastic<case_ct_params_t, n_dims>
@@ -99,7 +128,7 @@ namespace cases
         }
         BZ_DECLARE_FUNCTOR(th_std_fctr);
       };
-    
+
       // westerly wind
       struct u_t : hori_vel_t
       {
@@ -124,6 +153,56 @@ namespace cases
         }
         BZ_DECLARE_FUNCTOR(w_LS_fctr);
       };
+
+      // initial ambient chem profiles (mixing ratio?)
+      struct SO2g_fctr
+      {
+        real_t operator()(const real_t &z) const
+        {
+          return SO2g_dycoms(z);
+        }
+        BZ_DECLARE_FUNCTOR(SO2g_fctr);
+      };
+      struct H2O2g_fctr
+      {
+        real_t operator()(const real_t &z) const
+        {
+          return H2O2g_dycoms(z);
+        }
+        BZ_DECLARE_FUNCTOR(H2O2g_fctr);
+      };
+      struct CO2g_fctr
+      {
+        real_t operator()(const real_t &z) const
+        {
+          return CO2g_dycoms(z);
+        }
+        BZ_DECLARE_FUNCTOR(CO2g_fctr);
+      };
+      struct NH3g_fctr
+      {
+        real_t operator()(const real_t &z) const
+        {
+          return NH3g_dycoms(z);
+        }
+        BZ_DECLARE_FUNCTOR(NH3g_fctr);
+      };
+      struct HNO3g_fctr
+      {
+        real_t operator()(const real_t &z) const
+        {
+          return HNO3g_dycoms(z);
+        }
+        BZ_DECLARE_FUNCTOR(HNO3g_fctr);
+      };
+      struct O3g_fctr
+      {
+        real_t operator()(const real_t &z) const
+        {
+          return O3g_dycoms(z);
+        }
+        BZ_DECLARE_FUNCTOR(O3g_fctr);
+      };
 
       // density profile as a function of altitude
       // hydrostatic and assuming constant theta (not used now)
@@ -190,8 +269,11 @@ namespace cases
       }
 
       template <class index_t>
-      void intcond_hlpr(typename parent_t::concurr_any_t &concurr, arr_1D_t &rhod, int rng_seed, index_t index)
+      void intcond_hlpr(typename parent_t::concurr_any_t &concurr, arr_1D_t &rhod, arr_1D_t &p_e, int rng_seed, index_t index)
       {
+        namespace molar_mass  = libcloudphxx::common::molar_mass;
+        using libcloudphxx::common::moist_air::kaBoNA;
+
         int nz = concurr.advectee_global().extent(ix::w);  // ix::w is the index of vertical domension both in 2D and 3D
         real_t dz = (this->Z / si::metres) / (nz-1); 
 
@@ -224,6 +306,27 @@ namespace cases
           this->make_cyclic(th_global);
           concurr.advectee_global_set(th_global, ix::th);
         }
+
+        arr_1D_t mixr_hlpr(nz);
+        mixr_hlpr = p_e / // pressure
+                    real_t(kaBoNA<real_t>() / si::joules * si::kelvin * si::mole) // Blotzman * Avogadro
+                    / (th_std_fctr()(blitz::firstIndex() * dz) * calc_exner()(p_e))  / rhod(blitz::Range(0,nz-1));
+//        std::cerr << "mixr_hlpr: " << mixr_hlpr << std::endl;
+
+        // chemical composition
+        if constexpr(has_SO2g<ix>)
+          concurr.advectee(ix::SO2g)  = mixr_hlpr(index) * SO2g_fctr() (index * dz) * real_t(molar_mass::M_SO2<real_t>()  * si::moles / si::kilograms);
+        if constexpr(has_O3g<ix>)
+          concurr.advectee(ix::O3g)  = mixr_hlpr(index) * O3g_fctr() (index * dz) * real_t(molar_mass::M_O3<real_t>()  * si::moles / si::kilograms);
+        if constexpr(has_H2O2g<ix>)
+          concurr.advectee(ix::H2O2g)  = mixr_hlpr(index) * H2O2g_fctr() (index * dz) * real_t(molar_mass::M_H2O2<real_t>()  * si::moles / si::kilograms);
+        if constexpr(has_CO2g<ix>)
+          concurr.advectee(ix::CO2g)  = mixr_hlpr(index) * CO2g_fctr() (index * dz) * real_t(molar_mass::M_CO2<real_t>()  * si::moles / si::kilograms);
+        if constexpr(has_NH3g<ix>)
+          concurr.advectee(ix::NH3g)  = mixr_hlpr(index) * NH3g_fctr() (index * dz) * real_t(molar_mass::M_NH3<real_t>()  * si::moles / si::kilograms);
+        if constexpr(has_HNO3g<ix>)
+          concurr.advectee(ix::HNO3g)  = mixr_hlpr(index) * HNO3g_fctr() (index * dz) * real_t(molar_mass::M_HNO3<real_t>()  * si::moles / si::kilograms);
+
       }
 
       // calculate the initial environmental theta and rv profiles
@@ -353,7 +456,7 @@ namespace cases
                    arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, arr_1D_t &rl_e, arr_1D_t &p_e, int rng_seed)
       {
         blitz::secondIndex k;
-        this->intcond_hlpr(concurr, rhod, rng_seed, k);
+        this->intcond_hlpr(concurr, rhod, p_e, rng_seed, k);
       };
 
       // ctor
@@ -395,7 +498,7 @@ namespace cases
                    arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, arr_1D_t &rl_e, arr_1D_t &p_e, int rng_seed)
       {
         blitz::thirdIndex k;
-        this->intcond_hlpr(concurr, rhod, rng_seed, k);
+        this->intcond_hlpr(concurr, rhod, p_e, rng_seed, k);
 
         int nz = concurr.advectee_global().extent(ix::w);
         real_t dz = (this->Z / si::metres) / (nz-1); 

src/detail/ct_params.hpp

@@ -53,6 +53,38 @@ struct ct_params_3D_lgrngn : ct_params_common
   enum { delayed_step = libmpdataxx::opts::bit(ix::th) | libmpdataxx::opts::bit(ix::rv) };
 };
 
+struct ct_params_2D_lgrngn_chem : ct_params_common
+{
+  enum { n_dims = 2 };
+  enum { n_eqns = 10 };
+  struct ix { enum {
+    u, w, th, rv, 
+    SO2g, O3g, H2O2g, CO2g, NH3g, HNO3g,
+    vip_i=u, vip_j=w, vip_den=-1
+  }; };
+  enum { delayed_step = libmpdataxx::opts::bit(ix::th) | libmpdataxx::opts::bit(ix::rv) 
+                      | libmpdataxx::opts::bit(ix::SO2g) | libmpdataxx::opts::bit(ix::O3g)
+                      | libmpdataxx::opts::bit(ix::H2O2g) | libmpdataxx::opts::bit(ix::CO2g)
+                      | libmpdataxx::opts::bit(ix::NH3g) | libmpdataxx::opts::bit(ix::HNO3g)
+       };
+};
+
+struct ct_params_3D_lgrngn_chem : ct_params_common
+{
+  enum { n_dims = 3 };
+  enum { n_eqns = 11 };
+  struct ix { enum {
+    u, v, w, th, rv, 
+    SO2g, O3g, H2O2g, CO2g, NH3g, HNO3g,
+    vip_i=u, vip_j=v, vip_k=w, vip_den=-1
+  }; };
+  enum { delayed_step = libmpdataxx::opts::bit(ix::th) | libmpdataxx::opts::bit(ix::rv) 
+                      | libmpdataxx::opts::bit(ix::SO2g) | libmpdataxx::opts::bit(ix::O3g)
+                      | libmpdataxx::opts::bit(ix::H2O2g) | libmpdataxx::opts::bit(ix::CO2g)
+                      | libmpdataxx::opts::bit(ix::NH3g) | libmpdataxx::opts::bit(ix::HNO3g)
+       };
+};
+
 struct ct_params_2D_blk_2m : ct_params_common
 {
   enum { n_dims = 2 };
@@ -112,3 +144,34 @@ struct ct_params_3D_dry : ct_params_common
     vip_i=u, vip_j=v, vip_k=w, vip_den=-1
   }; };
 };
+
+// ct_params type traits
+template< typename, typename = void >
+constexpr bool has_SO2g{};
+template< typename T >
+constexpr bool has_SO2g<T, std::void_t<decltype(T::SO2g)>> = true;
+
+template< typename, typename = void >
+constexpr bool has_HNO3g{};
+template< typename T >
+constexpr bool has_HNO3g<T, std::void_t<decltype(T::HNO3g)>> = true;
+
+template< typename, typename = void >
+constexpr bool has_NH3g{};
+template< typename T >
+constexpr bool has_NH3g<T, std::void_t<decltype(T::NH3g)>> = true;
+
+template< typename, typename = void >
+constexpr bool has_CO2g{};
+template< typename T >
+constexpr bool has_CO2g<T, std::void_t<decltype(T::CO2g)>> = true;
+
+template< typename, typename = void >
+constexpr bool has_H2O2g{};
+template< typename T >
+constexpr bool has_H2O2g<T, std::void_t<decltype(T::H2O2g)>> = true;
+
+template< typename, typename = void >
+constexpr bool has_O3g{};
+template< typename T >
+constexpr bool has_O3g<T, std::void_t<decltype(T::O3g)>> = true;

src/detail/user_params.hpp

@@ -19,5 +19,6 @@ struct user_params_t
 
   bool relax_th_rv,
        window,
+       chem,
        relax_ccn = false; // relevant only for lgrngn micro, hence needs a default value as otherwise it might be undefined in blk_1m/blk_2m
 };

src/opts/opts_lgrngn.hpp

@@ -58,6 +58,8 @@ void setopts_micro(
     ("chem_dsl", po::value<bool>()->default_value(rt_params.cloudph_opts.chem_dsl) , "dissolving trace gases (1=on, 0=off)")
     ("chem_dsc", po::value<bool>()->default_value(rt_params.cloudph_opts.chem_dsc) , "dissociation           (1=on, 0=off)")
     ("chem_rct", po::value<bool>()->default_value(rt_params.cloudph_opts.chem_rct) , "aqueous chemistry      (1=on, 0=off)")
+    ("chem_rho", po::value<setup::real_t>()->default_value(1.8e3),    "dry particle density (used in chemistry) [kg/m3]")
+
     ("dev_count", po::value<int>()->default_value(0), "no. of CUDA devices")
     ("dev_id", po::value<int>()->default_value(-1), "CUDA backend - id of device to be used")
     // free parameters
@@ -70,6 +72,8 @@ void setopts_micro(
     // 
     ("out_dry", po::value<std::string>()->default_value(""),  "dry radius ranges and moment numbers (r1:r2|n1,n2...;...)")
     ("out_wet", po::value<std::string>()->default_value(""),  "wet radius ranges and moment numbers (r1:r2|n1,n2...;...)")
+    ("out_wet_pH", po::value<std::string>()->default_value("0:1|0"), "wet radius ranges for output of H+ and S_VI)")
+    ("out_chem", po::value<std::string>()->default_value("0:1|0"),   "dry radius ranges for which chem mass is outputted")
     ("gccn", po::value<setup::real_t>()->default_value(0) , "concentration of giant aerosols = gccn * VOCALS observations")
 //    ("unit_test", po::value<bool>()->default_value(false) , "very low number concentration for unit tests")
     ("adve_scheme", po::value<std::string>()->default_value("euler") , "one of: euler, implicit, pred_corr")
@@ -98,6 +102,7 @@ void setopts_micro(
 //  bool unit_test = vm["unit_test"].as<bool>();
   setup::real_t ReL = vm["ReL"].as<setup::real_t>();
 
+  rt_params.cloudph_opts_init.chem_switch = user_params.chem;
   rt_params.cloudph_opts_init.sd_conc = vm["sd_conc"].as<unsigned long long>();
   rt_params.cloudph_opts_init.sd_const_multi = vm["sd_const_multi"].as<double>();
 
@@ -389,6 +394,7 @@ void setopts_micro(
   rt_params.cloudph_opts.chem_dsl = vm["chem_dsl"].as<bool>();
   rt_params.cloudph_opts.chem_dsc = vm["chem_dsc"].as<bool>();
   rt_params.cloudph_opts.chem_rct = vm["chem_rct"].as<bool>();
+  rt_params.cloudph_opts_init.chem_rho = vm["chem_rho"].as<thrust_real_t>();
 
   rt_params.cloudph_opts_init.dev_count = vm["dev_count"].as<int>();
   rt_params.cloudph_opts_init.dev_id = vm["dev_id"].as<int>();
@@ -430,27 +436,28 @@ void setopts_micro(
   rt_params.cloudph_opts_init.subs_switch = rt_params.subsidence;
   rt_params.cloudph_opts.subs = rt_params.subsidence;
 
-  // parsing --out_dry and --out_wet options values
+  // parsing binned output options
   // the format is: "rmin:rmax|0,1,2;rmin:rmax|3;..."
-  for (auto &opt : std::set<std::string>({"out_dry", "out_wet"}))
+  for (auto &opt : std::set<std::string>({"out_dry", "out_wet", "out_chem", "out_wet_pH"}))
   {
     namespace qi = boost::spirit::qi;
     namespace phoenix = boost::phoenix;
 
+    if(!user_params.chem && (opt == "out_chem" || opt == "out_wet_pH")) 
+      continue;
+
+    outmom_t<thrust_real_t> moms;
+
     std::string val = vm[opt].as<std::string>();
     auto first = val.begin();
     auto last  = val.end();
 
     std::vector<std::pair<std::string, std::string>> min_maxnum;
-    outmom_t<thrust_real_t> &moms = 
-      opt == "out_dry"
-        ? rt_params.out_dry
-        : rt_params.out_wet;
 
     const bool result = qi::phrase_parse(first, last, 
       *(
-	*(qi::char_-":")  >>  qi::lit(":") >>  
-	*(qi::char_-";")  >> -qi::lit(";") 
+        *(qi::char_-":")  >>  qi::lit(":") >>  
+        *(qi::char_-";")  >> -qi::lit(";") 
       ),
       boost::spirit::ascii::space, min_maxnum
     );    
@@ -479,14 +486,27 @@ void setopts_micro(
       const bool result = qi::phrase_parse(
         nums_first, 
         nums_last, 
-	(
-	  qi::int_[phoenix::push_back(phoenix::ref(moms.back().second), qi::_1)]
-	      >> *(',' >> qi::int_[phoenix::push_back(phoenix::ref(moms.back().second), qi::_1)])
-	),
-	boost::spirit::ascii::space
+        (
+          qi::int_[phoenix::push_back(phoenix::ref(moms.back().second), qi::_1)]
+              >> *(',' >> qi::int_[phoenix::push_back(phoenix::ref(moms.back().second), qi::_1)])
+        ),
+        boost::spirit::ascii::space
       );    
+
+      if(opt == "out_dry")
+        rt_params.out_dry = moms;
+      if(opt == "out_wet")
+        rt_params.out_wet = moms;
+      if constexpr(has_SO2g<typename solver_t::ix>) // with chemistry
+      {
+        if(opt == "out_chem")
+          rt_params.out_chem = moms;
+        if(opt == "out_wet_pH")
+          rt_params.out_wet_pH = moms;
+      }
+
       if (!result || nums_first != nums_last) BOOST_THROW_EXCEPTION(po::validation_error(
-	  po::validation_error::invalid_option_value, opt, val // TODO: report only the relevant part?
+          po::validation_error::invalid_option_value, opt, val // TODO: report only the relevant part?
       ));  
     }
   }