Merge pull request #146 from pdziekan/rng_seed_init_switch_merge_mast…

…er_rd_min_max

Rng seed init switch merge master rd min max

src/detail/user_params.hpp

@@ -8,16 +8,15 @@
 struct user_params_t
 {
   int nt, outfreq, spinup, rng_seed, rng_seed_init;
-  bool relax_th_rv;
   setup::real_t X, Y, Z, dt;
   std::string outdir, model_case;
-//  bool th_src, rv_src, rc_src, rr_src, nc_src, nr_src, uv_src, w_src, relax_ccn, relax_th_rv;
   setup::real_t sgs_delta;
   quantity<si::length, setup::real_t> mean_rd1, mean_rd2;		
   quantity<si::dimensionless, setup::real_t> sdev_rd1, sdev_rd2;		
   quantity<power_typeof_helper<si::length, static_rational<-3>>::type, setup::real_t> n1_stp, n2_stp;		
   quantity<si::dimensionless, setup::real_t> kappa1, kappa2;
   quantity<si::dimensionless, setup::real_t> case_n_stp_multiplier;
 
-  bool relax_ccn = false; // relevant only for lgrngn micro, hence needs a default value as otherwise it might be undefined in blk_1m/blk_2m
+  bool relax_th_rv,
+       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

@@ -79,6 +79,8 @@ void setopts_micro(
     ("ReL", po::value<setup::real_t>()->default_value(100) , "taylor-microscale reynolds number (onishi kernel)")
     ("out_dry_spec", po::value<bool>()->default_value(false), "enable output for plotting dry spectrum")
     ("out_wet_spec", po::value<bool>()->default_value(false), "enable output for plotting wet spectrum")
+    ("rd_min", po::value<setup::real_t>()->default_value(rt_params.cloudph_opts_init.rd_min), "minimum dry radius of initialized droplets [m] (negative means automatic detection)")
+    ("rd_max", po::value<setup::real_t>()->default_value(rt_params.cloudph_opts_init.rd_max), "maximum dry radius of initialized droplets [m] (negative means automatic detection); sd_conc_large_tail==true may result in initialization of even larger droplets")
     ("relax_ccn", po::value<bool>()->default_value(false) , "add CCN if per-level mean of CCN concentration is lower than (case-specific) desired concentration")
     // TODO: MAC, HAC, vent_coef
   ;
@@ -99,6 +101,9 @@ void setopts_micro(
   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>();
 
+  rt_params.cloudph_opts_init.rd_min = vm["rd_min"].as<setup::real_t>();
+  rt_params.cloudph_opts_init.rd_max = vm["rd_max"].as<setup::real_t>();
+
   std::string adve_scheme_str = vm["adve_scheme"].as<std::string>();
   if (adve_scheme_str == "euler") rt_params.cloudph_opts_init.adve_scheme = libcloudphxx::lgrngn::as_t::euler;
   else if (adve_scheme_str == "implicit") rt_params.cloudph_opts_init.adve_scheme = libcloudphxx::lgrngn::as_t::implicit;
@@ -397,6 +402,7 @@ void setopts_micro(
 
   rt_params.cloudph_opts_init.rng_seed = user_params.rng_seed;
   rt_params.cloudph_opts_init.rng_seed_init = user_params.rng_seed_init;
+  rt_params.cloudph_opts_init.rng_seed_init_switch = true;
 
   // coalescence kernel choice
   if(!vm["turb_coal"].as<bool>())

src/solvers/lgrngn/hook_ante_loop_lgrngn.hpp

@@ -153,6 +153,7 @@ void slvr_lgrngn<ct_params_t>::hook_ante_loop(int nt)
     this->record_aux_const("diag_incloud_time", "lgrngn", params.cloudph_opts_init.diag_incloud_time);  
     this->record_aux_const("rng_seed", "lgrngn", params.cloudph_opts_init.rng_seed);  
     this->record_aux_const("rng_seed_init", "lgrngn", params.cloudph_opts_init.rng_seed_init);  
+    this->record_aux_const("rng_seed_init_switch", "lgrngn", params.cloudph_opts_init.rng_seed_init_switch);  
     this->record_aux_const("async", "lgrngn", params.async);  
     this->record_aux_const("adve", "lgrngn", params.cloudph_opts.adve);  
     this->record_aux_const("sedi", "lgrngn", params.cloudph_opts.sedi);  
@@ -192,6 +193,8 @@ void slvr_lgrngn<ct_params_t>::hook_ante_loop(int nt)
     this->record_aux_const("src_type", "lgrngn", libcloudphxx::lgrngn::src_name.at(params.cloudph_opts_init.src_type));  
     this->record_aux_const("terminal_velocity", "lgrngn", libcloudphxx::lgrngn::vt_name.at(params.cloudph_opts_init.terminal_velocity));  
     this->record_aux_const("RH_formula", "lgrngn", libcloudphxx::lgrngn::RH_formula_name.at(params.cloudph_opts_init.RH_formula));  
+    this->record_aux_const("rd_min", "lgrngn", params.cloudph_opts_init.rd_min);  
+    this->record_aux_const("rd_max", "lgrngn", params.cloudph_opts_init.rd_max);  
     this->record_aux_const("relax_ccn", "user_params", params.user_params.relax_ccn);  
   }
   this->mem->barrier();

tests/unit/api_test.cpp

@@ -18,7 +18,7 @@ int main(int ac, char** av)
   bool run_piggy = std::stoi(av[2]);
 
   string opts_common = 
-    "--outfreq=1000 --nt=2 --spinup=1 --dt=1 --serial=true --prs_tol=1e-3 --rng_seed=44 --case_n_stp_multiplier=1e-8"; 
+    "--outfreq=1000 --nt=2 --spinup=1 --dt=1 --serial=true --prs_tol=1e-3 --case_n_stp_multiplier=1e-8"; 
   vector<string> opts_dim({
     "--nx=4 --nz=4",
     "--nx=4 --nz=4 --X=1000 --Z=-1",
@@ -32,7 +32,8 @@ int main(int ac, char** av)
     "--async=false --micro=lgrngn --backend=serial --sd_conc=8",
     "--async=false --micro=lgrngn --backend=serial --sd_conc=8 --gccn=1",
     "--async=false --micro=lgrngn --backend=serial --sd_conc=8 --relax_ccn=1",
-    "--async=false --micro=lgrngn --backend=serial --sd_conc=8 --gccn=1 --relax_ccn=1"
+    "--async=false --micro=lgrngn --backend=serial --sd_conc=8 --gccn=1 --relax_ccn=1",
+    "--async=false --micro=lgrngn --backend=serial --sd_conc=8 --rd_min=1e-9 --rd_max=1e-6"
   });
   vector<string> opts_case({
     "--case=moist_thermal",
@@ -49,6 +50,10 @@ int main(int ac, char** av)
     "--piggy=0",
     "--piggy=0 --save_vel=1"
   });
+  vector<string> opts_rng({
+    "--rng_seed=44",
+    "--rng_seed=44 --rng_seed_init=55"
+  });
   // run the piggybacker, if required
   if(run_piggy)
     opts_piggy.push_back("--piggy=1 --vel_in=velocity_out.dat");  // take vel file from blk, cause it's ran first
@@ -61,46 +66,47 @@ int main(int ac, char** av)
   for (auto &opts_d : opts_dim)
     for (auto &opts_m : opts_micro)
       for (auto &opts_c : opts_case)
-        for (auto &opts_p : opts_piggy) // piggy has to be last to prevent overwriting of vel_out
-        {
-          if((opts_c == opts_case[1]) && (opts_d == opts_dim[2] || opts_d == opts_dim[3]))
-          {
-            std::cout << "skipping 3d dry thermal tests" << std::endl;
-            continue;
-          }
-          if((opts_c == opts_case[1]) && opts_m != opts_micro[0])
+        for (auto &opts_r : opts_rng)
+          for (auto &opts_p : opts_piggy) // piggy has to be last to prevent overwriting of vel_out
           {
-            std::cout << "skipping dry thermal tests with microphysics other than 'none'" << std::endl;
-            continue; 
-          }
-          if((opts_c == opts_case[7]) && opts_m != opts_micro[0])
-          {
-            std::cout << "skipping dry pbl tests with microphysics other than 'none'" << std::endl;
-            continue; 
-          }
+            if((opts_c == opts_case[1]) && (opts_d == opts_dim[2] || opts_d == opts_dim[3]))
+            {
+              std::cout << "skipping 3d dry thermal tests" << std::endl;
+              continue;
+            }
+            if((opts_c == opts_case[1]) && opts_m != opts_micro[0])
+            {
+              std::cout << "skipping dry thermal tests with microphysics other than 'none'" << std::endl;
+              continue; 
+            }
+            if((opts_c == opts_case[7]) && opts_m != opts_micro[0])
+            {
+              std::cout << "skipping dry pbl tests with microphysics other than 'none'" << std::endl;
+              continue; 
+            }
 
-          ostringstream cmd, opts;
-          opts << opts_common << " " << opts_m << " " << opts_d << " " << opts_c << " " << opts_p << " " << opts_additional;
-          // we want outdir=opts.str(), but that gives a long outdir that h5diff has trouble with reading and gives error when comparing const.h5 with refdata
-          // hence we hash opts to get outdir
-          auto outdir = std::hash<std::string>{}(opts.str());
-          ofdict << outdir << " : " << opts.str() << std::endl;
+            ostringstream cmd, opts;
+            opts << opts_common << " " << opts_m << " " << opts_d << " " << opts_c << " " << opts_r << " " << opts_p << " "<< opts_additional;
+            // we want outdir=opts.str(), but that gives a long outdir that h5diff has trouble with reading and gives error when comparing const.h5 with refdata
+            // hence we hash opts to get outdir
+            auto outdir = std::hash<std::string>{}(opts.str());
+            ofdict << outdir << " : " << opts.str() << std::endl;
 
-          cmd << av[1] <<  "/../../build/uwlcm " << opts.str() << " --outdir=\"output/" << outdir << "\"";
+            cmd << av[1] <<  "/../../build/uwlcm " << opts.str() << " --outdir=\"output/" << outdir << "\"";
 
-          cerr << endl << "=========" << endl;
-          notice_macro("about to call: " << cmd.str())
+            cerr << endl << "=========" << endl;
+            notice_macro("about to call: " << cmd.str())
 
-          if (EXIT_SUCCESS != system(cmd.str().c_str()))
-            error_macro("model run failed: " << cmd.str())
+            if (EXIT_SUCCESS != system(cmd.str().c_str()))
+              error_macro("model run failed: " << cmd.str())
 
-          // copy the stored velocity for the next run
-          if(opts_p == opts_piggy[1])
-          {
-            ostringstream cpcmd;
-            cpcmd << "cp \"output/" << outdir << "/velocity_out.dat\" .";
-            notice_macro("about to call: " << cpcmd.str())
-            system(cpcmd.str().c_str());
+            // copy the stored velocity for the next run
+            if(opts_p == opts_piggy[1])
+            {
+              ostringstream cpcmd;
+              cpcmd << "cp \"output/" << outdir << "/velocity_out.dat\" .";
+              notice_macro("about to call: " << cpcmd.str())
+              system(cpcmd.str().c_str());
+            }
           }
-        }
 }