Full domain init

include/libcloudph++/lgrngn/opts_init.hpp

@@ -53,7 +53,7 @@ namespace libcloudphxx
       // Lagrangian domain extents
       real_t x0, y0, z0, x1, y1, z1;
 
-      // no. of super-droplets per cell
+      // no. of super-droplets per cell (or in the entire domain if domain_sd_init==true)
       unsigned long long sd_conc; 
 
       // should more SDs be added to better represent large tail of the distribution
@@ -128,9 +128,9 @@ namespace libcloudphxx
 
       real_t rd_min, rd_max; // min/max dry radius of droplets [m] 
 
-      bool no_ccn_at_init; // if true, no ccn / SD are put at the start of the simulation
-
-      bool open_side_walls,       // if true, side walls are "open", i.e. SD are removed at contact. Periodic otherwise.
+      bool no_ccn_at_init,        // if true, no ccn / SD are put at the start of the simulation
+           domain_sd_init,        // if true, SDs are initialized in the entire domain as if it was a single cell (experimental, could mess with SD sources and other things...)
+           open_side_walls,       // if true, side walls are "open", i.e. SD are removed at contact. Periodic otherwise.
            periodic_topbot_walls; // if true, top and bot walls are periodic. Open otherwise
 
 
@@ -234,6 +234,7 @@ namespace libcloudphxx
         rd_max(-1),
         diag_incloud_time(false),
         no_ccn_at_init(false),
+        domain_sd_init(false),
         open_side_walls(false),
         periodic_topbot_walls(false),
         variable_dt_switch(false),

src/impl/particles_impl_dist_analysis.ipp

@@ -28,8 +28,9 @@ namespace libcloudphxx
           / sd_conc
           * dt
           * (n_dims == 0
-            ? dv[0]
-            : (opts_init.dx * opts_init.dy * opts_init.dz)
+            ? dv[0]                                                                          // 0D (box/parcel)
+            : (opts_init.dx * opts_init.dy * opts_init.dz)                                   // 1D 2D 3D
+              * (opts_init.domain_sd_init ? opts_init.nx * opts_init.ny * opts_init.nz : 1)  // if initializing in the entire domain, not per-cell
           );
 
         log_rd_min = log(rd_min);
@@ -52,8 +53,9 @@ namespace libcloudphxx
             / sd_conc
             * dt
             * (n_dims == 0
-              ? dv[0]
-              : (opts_init.dx * opts_init.dy * opts_init.dz)
+              ? dv[0]                                                                          // 0D (box/parcel)
+              : (opts_init.dx * opts_init.dy * opts_init.dz)                                   // 1D 2D 3D
+                * (opts_init.domain_sd_init ? opts_init.nx * opts_init.ny * opts_init.nz : 1)  // if initializing in the entire domain, not per-cell
             );
 
           log_rd_min = log(rd_min);

src/impl/particles_impl_init_SD_with_distros.ipp

@@ -52,6 +52,17 @@ namespace libcloudphxx
     template <typename real_t, backend_t device>
     void particles_t<real_t, device>::impl::init_SD_with_distros_finalize(const real_t &kappa, const bool unravel_ijk_switch)
     {
+      // ijk -> i, j, k
+      if(unravel_ijk_switch)
+        unravel_ijk(n_part_old);
+
+      // initialising particle positions
+      init_xyz();
+
+      // get correct ijk if domain-initializing
+      if(opts_init.domain_sd_init)
+        hskpng_ijk();
+
       // init kappa
       init_kappa(kappa);
 
@@ -77,13 +88,6 @@ namespace libcloudphxx
       if (opts_init.chem_switch){
         chem_vol_ante();
       }
-
-      // ijk -> i, j, k
-      if(unravel_ijk_switch)
-        unravel_ijk(n_part_old);
-
-      // initialising particle positions
-      init_xyz();
 
       if(opts_init.diag_incloud_time)
         init_incloud_time();

src/impl/particles_impl_init_count_num.ipp

@@ -31,7 +31,13 @@ namespace libcloudphxx
     template <typename real_t, backend_t device>
     void particles_t<real_t, device>::impl::init_count_num_sd_conc(const real_t &ratio)
     {
-      thrust::fill(count_num.begin(), count_num.end(), ratio * opts_init.sd_conc);
+      if(!opts_init.domain_sd_init) // usual per-cell init
+        thrust::fill(count_num.begin(), count_num.end(), ratio * opts_init.sd_conc);
+      else // if initializing in the entire domain, pretend that the first cell is the entire domain
+      {
+        thrust::fill(count_num.begin(), count_num.end(), 0);
+        *(count_num.begin()) = ratio * opts_init.sd_conc;
+      }
     }
 
     // calculate number of droplets in a cell from concentration [1/m^3], taking into account cell volume and air density
@@ -44,22 +50,49 @@ namespace libcloudphxx
       namespace arg = thrust::placeholders;
       using common::earth::rho_stp;
 
-      // cell volume
-      thrust::transform(
-        dv.begin(), dv.end(), 
-        arr.begin(), 
-        arr.begin(),
-        arg::_2 * arg::_1 
-      );
-
-      // correct for density with respect to STP
-      if(!opts_init.aerosol_independent_of_rhod)
+      if(!opts_init.domain_sd_init || n_dims==0) // usual per-cell init
+      {
+        // cell volume
         thrust::transform(
-          rhod.begin(), rhod.end(), 
+          dv.begin(), dv.end(), 
           arr.begin(), 
           arr.begin(),
-          arg::_1 / real_t(rho_stp<real_t>() / si::kilograms * si::cubic_metres) * arg::_2  
+          arg::_2 * arg::_1 
         );
+
+        // correct for density with respect to STP
+        if(!opts_init.aerosol_independent_of_rhod)
+          thrust::transform(
+            rhod.begin(), rhod.end(), 
+            arr.begin(), 
+            arr.begin(),
+            arg::_1 / real_t(rho_stp<real_t>() / si::kilograms * si::cubic_metres) * arg::_2  
+          );
+      }
+      else // full-domain init 1D 2D 3D
+      {
+        // all SDs are initialized as if they were in 0-th cell, so set zero number in other cells
+        thrust::fill(arr.begin()+1, arr.end(), 0);
+
+        // domain volume
+        thrust::transform(
+          arr.begin(), arr.begin()+1, 
+          arr.begin(), 
+          arg::_1
+            * (opts_init.x1 - opts_init.x0)
+            * (n_dims >= 2 ? opts_init.z1 - opts_init.z0 : 1)
+            * (n_dims == 3 ? opts_init.y1 - opts_init.y0 : 1)
+        );
+
+        // correct for density with respect to STP; TODO: 0-th cell density is used, although later these SDs are spread around entire domain... 
+        if(!opts_init.aerosol_independent_of_rhod)
+          thrust::transform(
+            rhod.begin(), rhod.begin()+1, 
+            arr.begin(), 
+            arr.begin(),
+            arg::_1 / real_t(rho_stp<real_t>() / si::kilograms * si::cubic_metres) * arg::_2  
+          );
+      }
     }
 
     template <typename real_t, backend_t device>

src/impl/particles_impl_init_n.ipp

@@ -93,7 +93,7 @@ namespace libcloudphxx
 
 
         // adjust to cell volume
-        if(n_dims > 0)
+        if(n_dims > 0 && !opts_init.domain_sd_init)
         {
           // rhod not needed anymore, reuse it to store dv on host
           thrust::copy(

src/impl/particles_impl_init_sanity_check.ipp

@@ -141,6 +141,11 @@ namespace libcloudphxx
         throw std::runtime_error("libcloudph++: rlx_timescale <= 0");        
       if(opts_init.rlx_switch && opts_init.chem_switch)
         throw std::runtime_error("libcloudph++: CCN relaxation does not work with chemistry");
+
+      if(opts_init.domain_sd_init && opts_init.dry_sizes.size() > 0)
+        throw std::runtime_error("libcloudph++: entire domain initialization (domain_sd_init==true) does not work with initializing from a dry sizes dictionary (dry_sizes.size > 0)");
+      if(opts_init.domain_sd_init && opts_init.sd_const_multi > 0)
+        throw std::runtime_error("libcloudph++: entire domain initialization (domain_sd_init==true) does not work with the constant-mutliplicity init (sd_const_multi > 0)");
     }
   };
 };

src/impl/particles_impl_init_wet.ipp

@@ -44,7 +44,7 @@ namespace libcloudphxx
     {
       // initialising values of rw2
       {
-	// calculating rw_eq
+        // calculating rw_eq
         {
           typedef thrust::permutation_iterator<
             typename thrust_device::vector<real_t>::iterator,
@@ -67,11 +67,11 @@ namespace libcloudphxx
             pi_t(T.begin(),  ijk.begin() + n_part_old)
           ));
 
-	  thrust::transform(
-	    zip_it, zip_it + n_part_to_init, // input
-	    rw2.begin() + n_part_old, // output
+          thrust::transform(
+            zip_it, zip_it + n_part_to_init, // input
+            rw2.begin() + n_part_old, // output
             detail::rw2_eq<real_t>(opts_init.RH_max) 
-	  );
+          );
         }
       }
     }

src/impl/particles_impl_init_xyz.ipp

@@ -29,7 +29,7 @@ namespace libcloudphxx
           using std::min;
           using std::max;
 #endif
-        	
+                
           return u01 * min(p1, (ii+1) * dp) + (1. - u01) * max(p0, ii * dp); 
         }
       };
@@ -55,17 +55,27 @@ namespace libcloudphxx
         // tossing random numbers
         rand_u01(n_part_to_init);
 
-	// shifting from [0,1] to random position within respective cell 
-  // TODO: now the rand range is [0,1), include this here
+        // shifting from [0,1] to random position within respective cell 
+        // TODO: now the rand range is [0,1), include this here
+        if(!opts_init.domain_sd_init)
         {
-          namespace arg = thrust::placeholders;
-	  thrust::transform(
-	    u01.begin(), 
-	    u01.begin() + n_part_to_init,
+          thrust::transform(
+            u01.begin(), 
+            u01.begin() + n_part_to_init,
             ii[ix]->begin() + n_part_old, 
-	    v[ix]->begin() + n_part_old, 
+            v[ix]->begin() + n_part_old, 
             detail::pos_lgrngn_domain<real_t>(a[ix], b[ix], d[ix])
-	  );
+          );
+        }
+        else // random position in the entire domain
+        {
+          namespace arg = thrust::placeholders;
+          thrust::transform(
+            u01.begin(), 
+            u01.begin() + n_part_to_init,
+            v[ix]->begin() + n_part_old, 
+            a[ix] + arg::_1 * (b[ix] - a[ix])
+          );
         }
       }
     }