Merge pull request #415 from pdziekan/pos_rng_seed

pdziekan · web-flow · commit 904b6adb1384 · 2020-12-21T15:35:00.000+01:00
rng seed for initialization
diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp
@@ -305,6 +305,7 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .def_readwrite("chem_rho", &lgr::opts_init_t<real_t>::chem_rho)
       .def_readwrite("RH_max", &lgr::opts_init_t<real_t>::RH_max)
       .def_readwrite("rng_seed", &lgr::opts_init_t<real_t>::rng_seed)
+      .def_readwrite("rng_seed_init", &lgr::opts_init_t<real_t>::rng_seed_init)
       .def_readwrite("diag_incloud_time", &lgr::opts_init_t<real_t>::diag_incloud_time)
       .add_property("w_LS", &lgrngn::get_w_LS<real_t>, &lgrngn::set_w_LS<real_t>)
       .add_property("SGS_mix_len", &lgrngn::get_SGS_mix_len<real_t>, &lgrngn::set_SGS_mix_len<real_t>)
diff --git a/include/libcloudph++/lgrngn/opts_init.hpp b/include/libcloudph++/lgrngn/opts_init.hpp
@@ -100,17 +100,16 @@ namespace libcloudphxx
       // coalescence kernel parameters
       std::vector<real_t> kernel_parameters;
 
-      // chem
       bool chem_switch,  // if false no chemical reactions throughout the whole simulation (no memory allocation)
            coal_switch,  // if false no coalescence throughout the whole simulation
            sedi_switch,  // if false no sedimentation throughout the whole simulation
            subs_switch,  // if false no subsidence throughout the whole simulation
            src_switch,   // if false no source throughout the whole simulation
-           exact_sstp_cond, // if true, use per-particle sstp_cond logic, if false, use per-cell
            turb_adve_switch,   // if true, turbulent motion of SDs is modeled
            turb_cond_switch,   // if true, turbulent condensation of SDs is modeled
-           turb_coal_switch;   // if true, turbulent coalescence kernels can be used
-
+           turb_coal_switch,   // if true, turbulent coalescence kernels can be used
+           exact_sstp_cond;    // if true, use per-particle sstp_cond logic, if false, use per-cell
+           
       int sstp_chem;
       real_t chem_rho;
 
@@ -121,7 +120,8 @@ namespace libcloudphxx
       real_t RH_max;
 
       // rng seed
-      int rng_seed;
+      int rng_seed,
+          rng_seed_init; // seed used to init SD (positions and dry sizes)
 
       // no of GPUs per MPI node to use, 0 for all available
       int dev_count; 
@@ -168,6 +168,7 @@ namespace libcloudphxx
         RH_max(.95), // value seggested in Lebo and Seinfeld 2011
         chem_rho(0), // dry particle density  //TODO add checking if the user gave a different value (np w init)  (was 1.8e-3)
         rng_seed(44),
+        rng_seed_init(rng_seed),
         terminal_velocity(vt_t::undefined),
         kernel(kernel_t::undefined),
         adve_scheme(as_t::implicit),
diff --git a/src/detail/urand.hpp b/src/detail/urand.hpp
@@ -20,150 +20,162 @@ namespace libcloudphxx
       class rng
       {
 #if !defined(__NVCC__)
-	// serial version using C++11's <random>
-	using engine_t = std::mt19937;
+        // serial version using C++11's <random>
+        using engine_t = std::mt19937;
         using dist_u01_t = std::uniform_real_distribution<real_t>;
         using dist_normal01_t = std::normal_distribution<real_t>;
         using dist_un_t = std::uniform_int_distribution<unsigned int>;
-	engine_t engine;
-	dist_u01_t dist_u01;
-	dist_normal01_t dist_normal01;
-	dist_un_t dist_un;
+        engine_t engine;
+        dist_u01_t dist_u01;
+        dist_normal01_t dist_normal01;
+        dist_un_t dist_un;
 
-	struct fnctr_u01
-	{
+        struct fnctr_u01
+        {
           engine_t &engine;
           dist_u01_t &dist_u01;
-	  real_t operator()() { return dist_u01(engine); }
-	};
+          real_t operator()() { return dist_u01(engine); }
+        };
 
-	struct fnctr_normal01
-	{
+        struct fnctr_normal01
+        {
           engine_t &engine;
           dist_normal01_t &dist_normal01;
-	  real_t operator()() { return dist_normal01(engine); }
-	};
+          real_t operator()() { return dist_normal01(engine); }
+        };
 
-	struct fnctr_un
-	{
+        struct fnctr_un
+        {
           engine_t &engine;
           dist_un_t &dist_un;
-	  real_t operator()() { return dist_un(engine); }
-	};
+          real_t operator()() { return dist_un(engine); }
+        };
 
-	public:
+        public:
 
         // ctor
         rng(int seed) : engine(seed), dist_u01(0,1), dist_normal01(0,1), dist_un(0, std::numeric_limits<unsigned int>::max()) {}
 
-	void generate_n(
-	  thrust_device::vector<real_t> &u01, 
-	  const thrust_size_t n
-	) {
+        void reseed(int seed)
+        {
+          engine.seed(seed);
+        }
+
+        void generate_n(
+          thrust_device::vector<real_t> &u01, 
+          const thrust_size_t n
+        ) {
           // note: generate_n copies the third argument!!!
-	  std::generate_n(u01.begin(), n, fnctr_u01({engine, dist_u01})); 
-	}
+          std::generate_n(u01.begin(), n, fnctr_u01({engine, dist_u01})); 
+        }
 
-	void generate_normal_n(
-	  thrust_device::vector<real_t> &normal01, 
-	  const thrust_size_t n
-	) {
+        void generate_normal_n(
+          thrust_device::vector<real_t> &normal01, 
+          const thrust_size_t n
+        ) {
           // note: generate_n copies the third argument!!!
-	  std::generate_n(normal01.begin(), n, fnctr_normal01({engine, dist_normal01})); 
-	}
+          std::generate_n(normal01.begin(), n, fnctr_normal01({engine, dist_normal01})); 
+        }
 
-	void generate_n(
-	  thrust_device::vector<unsigned int> &un, 
-	  const thrust_size_t n
-	) {
+        void generate_n(
+          thrust_device::vector<unsigned int> &un, 
+          const thrust_size_t n
+        ) {
           // note: generate_n copies the third argument!!!
-	  std::generate_n(un.begin(), n, fnctr_un({engine, dist_un})); 
-	}
+          std::generate_n(un.begin(), n, fnctr_un({engine, dist_un})); 
+        }
 #endif
       };
  
       template <typename real_t>
       class rng<real_t, CUDA>
       {
 #if defined(__NVCC__)
-	// CUDA parallel version using curand
+        // CUDA parallel version using curand
 
-	// private member fields
-	curandGenerator_t gen;
-	
-	public:
+        // private member fields
+        curandGenerator_t gen;
+        
+        public:
 
-	rng(int seed)
-	{
+        rng(int seed)
+        {
           {
-	    int status = curandCreateGenerator(&gen, CURAND_RNG_PSEUDO_MTGP32);
-	    assert(status == CURAND_STATUS_SUCCESS /* && "curandCreateGenerator failed"*/);
-	    _unused(status);
+            int status = curandCreateGenerator(&gen, CURAND_RNG_PSEUDO_MTGP32);
+            assert(status == CURAND_STATUS_SUCCESS /* && "curandCreateGenerator failed"*/);
+            _unused(status);
           }
           {
-	    int status = curandSetPseudoRandomGeneratorSeed(gen, seed);
-	    assert(status == CURAND_STATUS_SUCCESS /* && "curandSetPseudoRandomGeneratorSeed failed"*/);
+            int status = curandSetPseudoRandomGeneratorSeed(gen, seed);
+            assert(status == CURAND_STATUS_SUCCESS /* && "curandSetPseudoRandomGeneratorSeed failed"*/);
             _unused(status);
-	  }
+          }
+        }
+
+        void reseed(int seed)
+        {
+          int status = curandSetPseudoRandomGeneratorSeed(gen, seed);
+          assert(status == CURAND_STATUS_SUCCESS /* && "curandSetPseudoRandomGeneratorSeed failed"*/);
+          _unused(status);
         }
 
-	~rng()
-	{
-	  int status = curandDestroyGenerator(gen); 
-	  assert(status == CURAND_STATUS_SUCCESS /* && "curandDestroyGenerator failed"*/);
+        ~rng()
+        {
+          int status = curandDestroyGenerator(gen); 
+          assert(status == CURAND_STATUS_SUCCESS /* && "curandDestroyGenerator failed"*/);
           _unused(status);
-	}
-
-	void generate_n(
-	  thrust_device::vector<float> &v, 
-	  const thrust_size_t n
-	)
-	{
-	  int status = curandGenerateUniform(gen, thrust::raw_pointer_cast(v.data()), n);
+        }
+
+        void generate_n(
+          thrust_device::vector<float> &v, 
+          const thrust_size_t n
+        )
+        {
+          int status = curandGenerateUniform(gen, thrust::raw_pointer_cast(v.data()), n);
           assert(status == CURAND_STATUS_SUCCESS /* && "curandGenerateUniform failed"*/);
           _unused(status);
 
-	}
+        }
 
-	void generate_n(
-	  thrust_device::vector<double> &v, 
-	  const thrust_size_t n
-	)
-	{
-	  int status = curandGenerateUniformDouble(gen, thrust::raw_pointer_cast(v.data()), n);
-	  assert(status == CURAND_STATUS_SUCCESS /* && "curandGenerateUniform failed"*/);
+        void generate_n(
+          thrust_device::vector<double> &v, 
+          const thrust_size_t n
+        )
+        {
+          int status = curandGenerateUniformDouble(gen, thrust::raw_pointer_cast(v.data()), n);
+          assert(status == CURAND_STATUS_SUCCESS /* && "curandGenerateUniform failed"*/);
           _unused(status);
-	}
-
-	void generate_normal_n(
-	  thrust_device::vector<float> &v, 
-	  const thrust_size_t n
-	)
-	{
-	  int status = curandGenerateNormal(gen, thrust::raw_pointer_cast(v.data()), n, float(0), float(1));
-	  assert(status == CURAND_STATUS_SUCCESS /* && "curandGenerateUniform failed"*/);
+        }
+
+        void generate_normal_n(
+          thrust_device::vector<float> &v, 
+          const thrust_size_t n
+        )
+        {
+          int status = curandGenerateNormal(gen, thrust::raw_pointer_cast(v.data()), n, float(0), float(1));
+          assert(status == CURAND_STATUS_SUCCESS /* && "curandGenerateUniform failed"*/);
           _unused(status);
-	}
-
-	void generate_normal_n(
-	  thrust_device::vector<double> &v, 
-	  const thrust_size_t n
-	)
-	{
-	  int status = curandGenerateNormalDouble(gen, thrust::raw_pointer_cast(v.data()), n, double(0), double(1));
-	  assert(status == CURAND_STATUS_SUCCESS /* && "curandGenerateUniform failed"*/);
+        }
+
+        void generate_normal_n(
+          thrust_device::vector<double> &v, 
+          const thrust_size_t n
+        )
+        {
+          int status = curandGenerateNormalDouble(gen, thrust::raw_pointer_cast(v.data()), n, double(0), double(1));
+          assert(status == CURAND_STATUS_SUCCESS /* && "curandGenerateUniform failed"*/);
           _unused(status);
-	}
-
-	void generate_n(
-	  thrust_device::vector<unsigned int> &v, 
-	  const thrust_size_t n
-	)
-	{
-	  int status = curandGenerate(gen, thrust::raw_pointer_cast(v.data()), n);
-	  assert(status == CURAND_STATUS_SUCCESS /* && "curandGenerateUniform failed"*/);
+        }
+
+        void generate_n(
+          thrust_device::vector<unsigned int> &v, 
+          const thrust_size_t n
+        )
+        {
+          int status = curandGenerate(gen, thrust::raw_pointer_cast(v.data()), n);
+          assert(status == CURAND_STATUS_SUCCESS /* && "curandGenerateUniform failed"*/);
           _unused(status);
-	}
+        }
 #endif
       };
     };
diff --git a/src/particles_init.ipp b/src/particles_init.ipp
@@ -27,6 +27,9 @@ namespace libcloudphxx
     {
       pimpl->init_sanity_check(th, rv, rhod, p, courant_x, courant_y, courant_z, ambient_chem);
 
+      // set rng seed to be used during init
+      pimpl->rng.reseed(pimpl->opts_init.rng_seed_init);
+
       // is a constant pressure profile used?
       pimpl->const_p = !p.is_null();
       // if pressure comes from a profile, sstp_tmp_p also needs to be copied between distributed memories
@@ -120,6 +123,9 @@ namespace libcloudphxx
 
       // init count_num and count_ijk
       pimpl->hskpng_count();
+
+      // set rng seed to be used after init
+      pimpl->rng.reseed(pimpl->opts_init.rng_seed);
     }
   };
 };
diff --git a/tests/python/physics/coalescence_golovin.py b/tests/python/physics/coalescence_golovin.py
@@ -136,7 +136,7 @@ def calc_golovin(res,t,n0,v0,b):
   
   if(i==0):
     print('sd_conc RMSD = ' + str(rmsd));
-    limit = 1e-5;
+    limit = 1.1e-5;
   else:
     print('const_multi RMSD = ' + str(rmsd));
     limit = 3e-5; # constant multiplicity doesn't represent tails of the distribution so well and mass densty function depends on large tail?
diff --git a/tests/python/unit/api_lgrngn.py b/tests/python/unit/api_lgrngn.py
@@ -33,6 +33,7 @@ def lognormal(lnr):
 opts_init.sd_conc = 64
 opts_init.n_sd_max = int(1e6) # some space for tail SDs
 opts_init.rng_seed = 396
+opts_init.rng_seed_init = 456
 opts_init.src_dry_distros = {kappa1:lognormal}
 opts_init.src_sd_conc = 64
 opts_init.src_z1 = opts_init.dz
@@ -54,6 +55,7 @@ def lognormal(lnr):
 print("sd_conc = ", opts_init.sd_conc)
 print("RH_max = ", opts_init.RH_max)
 print("rng_seed = ", opts_init.rng_seed)
+print("rng_seed_init = ", opts_init.rng_seed_init)
 print("kernel =", opts_init.kernel) 
 print("sd_conc =", opts_init.sd_conc)
 print("terminal_velocity =", opts_init.terminal_velocity)
