Merge pull request #450 from pdziekan/aerosol_profile

aerosol conc factor option

.github/workflows/test_libclouphxx.yml

@@ -23,7 +23,7 @@ jobs:
     - uses: actions/checkout@v2
 
     - name: build libcloudph++
-      uses: igfuw/libcloudphxx_build@v0.2-beta
+      uses: igfuw/libcloudphxx_build@v0.3
       with:
         disable_cuda: ${{matrix.disable_cuda}}
         build_type: ${{matrix.build_type}}
@@ -53,7 +53,7 @@ jobs:
     - uses: actions/checkout@v2
 
     - name: build libcloudph++
-      uses: igfuw/libcloudphxx_build@v0.2-beta
+      uses: igfuw/libcloudphxx_build@v0.3
       with:
         disable_cuda: ${{matrix.disable_cuda}}
         build_type: ${{matrix.build_type}}
@@ -100,7 +100,7 @@ jobs:
       run: tar -xvf build.tar
 
     - name: load UWLCM Singularity image
-      uses: igfuw/load_UWLCM_singularity_image@v0.1
+      uses: igfuw/load_UWLCM_singularity_image@v0.1b
       with:
         path: ${{ github.workspace }}/singularity_images
         tag: ${{ matrix.tag }}
@@ -145,7 +145,7 @@ jobs:
         path: libmpdataxx
 
     - name: Install libmpdata++
-      uses: igfuw/libmpdataxx_install@v0.4
+      uses: igfuw/libmpdataxx_install@v0.5
       with:
         build_type: ${{matrix.build_type}}
         threads: 4
@@ -204,7 +204,7 @@ jobs:
       run: tar -xvf build.tar
 
     - name: load UWLCM Singularity image
-      uses: igfuw/load_UWLCM_singularity_image@v0.1
+      uses: igfuw/load_UWLCM_singularity_image@v0.1b
       with:
         path: ${{ github.workspace }}/singularity_images
 
@@ -263,7 +263,7 @@ jobs:
     - run: grep diag_accr include/libcloudph++/lgrngn/particles.hpp
 
     - name: build libcloudph++
-      uses: igfuw/libcloudphxx_build@v0.2-beta
+      uses: igfuw/libcloudphxx_build@v0.3
       with:
         disable_cuda: ${{matrix.disable_cuda}}
         build_type: ${{matrix.build_type}}

include/libcloudph++/lgrngn/opts_init.hpp

@@ -60,6 +60,7 @@ namespace libcloudphxx
       bool sd_conc_large_tail;
 
       // should aerosol concentration init be independent of rhod (assumed to be in cm^{-3} and not at STP)
+      // NOTE: it also affects aerosol sources (not sure if all types of sources...)
       bool aerosol_independent_of_rhod;
 
       // is it allowed to change dt during simulation through opts.dt
@@ -126,6 +127,10 @@ namespace libcloudphxx
       // SGS mixing length profile [m]
       std::vector<real_t> SGS_mix_len;
 
+      // profile of aerosol concentration factor
+      // NOTE: it also affects aerosol sources (not sure if all types of sources...)
+      std::vector<real_t> aerosol_conc_factor;
+
       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

src/impl/particles_impl.ipp

@@ -151,6 +151,7 @@ namespace libcloudphxx
 
       thrust_device::vector<real_t> w_LS; // large-scale subsidence velocity profile
       thrust_device::vector<real_t> SGS_mix_len; // SGS mixing length profile
+      thrust_device::vector<real_t> aerosol_conc_factor; // profile of aerosol concentration factor
 
       // time steps to be used, considering that opts.dt can override opts_init.dt
       real_t dt;
@@ -335,6 +336,7 @@ namespace libcloudphxx
         ),
         w_LS(_opts_init.w_LS),
         SGS_mix_len(_opts_init.SGS_mix_len),
+        aerosol_conc_factor(_opts_init.aerosol_conc_factor),
         adve_scheme(_opts_init.adve_scheme),
         allow_sstp_cond(_opts_init.sstp_cond > 1 || _opts_init.variable_dt_switch),
         allow_sstp_chem(_opts_init.sstp_chem > 1 || _opts_init.variable_dt_switch),

src/impl/particles_impl_init_count_num.ipp

@@ -60,6 +60,18 @@ namespace libcloudphxx
           arr.begin(),
           arg::_1 / real_t(rho_stp<real_t>() / si::kilograms * si::cubic_metres) * arg::_2  
         );
+
+      // accounting for the aerosol concentration profile
+      if(opts_init.aerosol_conc_factor.size()>0)
+        thrust::transform(
+          arr.begin(), arr.end(),            // input - 1st arg
+          thrust::make_permutation_iterator( // input - 2nd arg
+            aerosol_conc_factor.begin(),
+            thrust::make_transform_iterator(thrust::make_counting_iterator<thrust_size_t>(0), arg::_1 % opts_init.nz) // k index
+          ),
+          arr.begin(),                       // output
+          arg::_1 * arg::_2
+        );
     }
 
     template <typename real_t, backend_t device>

src/impl/particles_impl_init_n.ipp

@@ -91,6 +91,18 @@ namespace libcloudphxx
             arg::_1 * arg::_2 / real_t(rho_stp<real_t>() / si::kilograms * si::cubic_metres)
           );
 
+        // accounting for the aerosol concentration profile
+        if(opts_init.aerosol_conc_factor.size()>0)
+          thrust::transform(
+            tmp_real.begin(), tmp_real.end(),            // input - 1st arg
+            thrust::make_permutation_iterator( // input - 2nd arg
+              opts_init.aerosol_conc_factor.begin(), 
+              thrust::make_transform_iterator(tmp_ijk.begin(), arg::_1 % opts_init.nz) // k index
+            ),
+            tmp_real.begin(),                       // output
+            arg::_1 * arg::_2
+          );
+
 
         // adjust to cell volume
         if(n_dims > 0)

src/impl/particles_impl_init_sanity_check.ipp

@@ -125,6 +125,15 @@ namespace libcloudphxx
         throw std::runtime_error("libcloudph++: at least one of opts_init.turb_adve_switch, opts_init.turb_cond_switch is true, but SGS mixing length profile size != nz");
       if(opts_init.SGS_mix_len.size() > 0 && *std::min(opts_init.SGS_mix_len.begin(), opts_init.SGS_mix_len.end()) <= 0)
         throw std::runtime_error("libcloudph++: SGS_mix_len <= 0");
+      if (!opts_init.aerosol_conc_factor.empty() && n_dims<2)
+        throw std::runtime_error("libcloudph++: aerosol_conc_factor can only be used in 2D and 3D");
+      if (!opts_init.aerosol_conc_factor.empty() && opts_init.nz != opts_init.aerosol_conc_factor.size())
+        throw std::runtime_error("libcloudph++: aerosol_conc_factor size needs to be either 0 or nz");
+      if (!opts_init.aerosol_conc_factor.empty() && opts_init.aerosol_independent_of_rhod==false)
+      {
+        std::cerr << "aerosol conc factor size: " << opts_init.aerosol_conc_factor.size() << std::endl;
+        throw std::runtime_error("libcloudph++: aerosol_conc_factor can only be used if aerosol_independent_of_rhod==true");
+        }
       #if defined(USE_MPI)
         if(opts_init.rlx_switch)
           std::cerr << "libcloudph++ WARNING: relaxation is not fully supported in MPI runs. Mean calculation and addition of SD will be done locally on each node." << std::endl;