Made AddSupersAtDomainTop GPU compatible

examples/eurec4a1d/src/config/eurec4a1d_config.yaml

@@ -26,7 +26,7 @@
 domain:
   nspacedims : 1                                                 # no. of spatial dimensions to model
   ngbxs : 60                                                     # total number of Gbxs
-  maxnsupers: 15360                                              # maximum number of SDs
+  maxnsupers: 32768                                              # maximum number of SDs
 
 timesteps:
   CONDTSTEP : 0.1                                                # time between SD condensation [s]
@@ -44,6 +44,7 @@ inputfiles:
 initsupers:
   type: frombinary                                               # type of initialisation of super-droplets
   initsupers_filename : ./share/eurec4a1d_ddimlessSDsinit.dat    # binary filename for initialisation of SDs
+  initnsupers: 15360                                             # initial no. of super-droplets to initialise
 
 ### Output Parameters ###
 outputdata:
@@ -69,3 +70,18 @@ coupled_dynamics:
   qvap : ./share/eurec4a1d_dimlessthermo_qvap.dat               # binary filename for vapour mixing ratio
   qcond : ./share/eurec4a1d_dimlessthermo_qcond.dat             # binary filename for liquid mixing ratio
   wvel : ./share/eurec4a1d_dimlessthermo_wvel.dat               # binary filename for vertical (coord3) velocity
+
+### Bounday Conditions Parameters ###
+boundary_conditions:
+  type: addsupersatdomaintop
+  COORD3LIM: 800                                       # SDs added to domain with coord3 >= COORD3LIM [m]
+  newnsupers: 1024                                     # number SDs to add to each gridbox above COORD3LIM
+  DRYRADIUS: 1e-9                                      # dry radius of new super-droplets (for solute mass) [m]
+  MINRADIUS: 1e-8                                      # minimum radius of new super-droplets [m]
+  MAXRADIUS: 1e-4                                      # maximum radius of new super-droplets [m]
+  NUMCONC_a: 2e8                                       # number conc. of 1st droplet lognormal dist [m^-3]
+  GEOMEAN_a: 0.2e-6                                    # geometric mean radius of 1st lognormal dist [m]
+  geosigma_a: 2.3                                      # geometric standard deviation of 1st lognormal dist
+  NUMCONC_b: 4e8                                       # number conc. of 2nd droplet lognormal dist [m^-3]
+  GEOMEAN_b: 3.5e-6                                    # geometric mean radius of 2nd lognormal dist [m]
+  geosigma_b: 2.0                                      # geometric standard deviation of 2nd lognormal dist

examples/eurec4a1d/src/main_eurec4a1D.cpp

@@ -9,7 +9,7 @@
  * Author: Clara Bayley (CB)
  * Additional Contributors:
  * -----
- * Last Modified: Friday 19th April 2024
+ * Last Modified: Saturday 4th May 2024
  * Modified By: CB
  * -----
  * License: BSD 3-Clause "New" or "Revised" License
@@ -29,6 +29,7 @@
 #include <stdexcept>
 #include <string_view>
 
+#include "cartesiandomain/add_supers_at_domain_top.hpp"
 #include "cartesiandomain/cartesianmaps.hpp"
 #include "cartesiandomain/cartesianmotion.hpp"
 #include "cartesiandomain/createcartesianmaps.hpp"
@@ -38,6 +39,7 @@
 #include "gridboxes/gridboxmaps.hpp"
 #include "initialise/config.hpp"
 #include "initialise/init_all_supers_from_binary.hpp"
+#include "initialise/init_supers_from_binary.hpp"
 #include "initialise/initgbxsnull.hpp"
 #include "initialise/initialconditions.hpp"
 #include "initialise/timesteps.hpp"
@@ -76,7 +78,8 @@ inline CoupledDynamics auto create_coupldyn(const Config &config, const Cartesia
 }
 
 inline InitialConditions auto create_initconds(const Config &config) {
-  const InitAllSupersFromBinary initsupers(config.get_initsupersfrombinary());
+  // const InitAllSupersFromBinary initsupers(config.get_initsupersfrombinary());
+  const InitSupersFromBinary initsupers(config.get_initsupersfrombinary());
   const InitGbxsNull initgbxs(config.get_ngbxs());
 
   return InitConds(initsupers, initgbxs);
@@ -94,7 +97,8 @@ inline auto create_movement(const Config &config, const Timesteps &tsteps,
   const Motion<CartesianMaps> auto motion =
       CartesianMotion(tsteps.get_motionstep(), &step2dimlesstime, terminalv);
 
-  const auto boundary_conditions = NullBoundaryConditions{};
+  // const auto boundary_conditions = NullBoundaryConditions{};
+  const auto boundary_conditions = AddSupersAtDomainTop(config.get_addsupersatdomaintop());
 
   return MoveSupersInDomain(gbxmaps, motion, boundary_conditions);
 }

libs/cartesiandomain/add_supers_at_domain_top.cpp

@@ -9,7 +9,7 @@
  * Author: Clara Bayley (CB)
  * Additional Contributors:
  * -----
- * Last Modified: Wednesday 1st May 2024
+ * Last Modified: Saturday 4th May 2024
  * Modified By: CB
  * -----
  * License: BSD 3-Clause "New" or "Revised" License
@@ -22,11 +22,41 @@
 
 #include "cartesiandomain/add_supers_at_domain_top.hpp"
 
+Kokkos::View<unsigned int *> remove_superdrops_from_gridboxes(const CartesianMaps &gbxmaps,
+                                                              const viewd_gbx d_gbxs,
+                                                              const double coord3lim);
+viewd_supers create_newsupers_for_gridboxes(const CartesianMaps &gbxmaps,
+                                            const CreateSuperdrop &create_superdrop,
+                                            Kokkos::View<unsigned int *> gbxindexes,
+                                            const size_t newnsupers_pergbx);
+void add_superdrops_for_gridboxes(const viewd_supers totsupers, const viewd_constgbx d_gbxs,
+                                  const viewd_constsupers newsupers);
+void move_supers_between_gridboxes_again(const viewd_gbx d_gbxs, const viewd_supers totsupers);
+
+/*
+Call to apply boundary conditions to remove and then add superdroplets to the top of the domain
+above coord3lim.
+
+_Note:_ totsupers is view of all superdrops (both in and out of bounds of domain).
+*/
+void AddSupersAtDomainTop::operator()(const CartesianMaps &gbxmaps, viewd_gbx d_gbxs,
+                                      const viewd_supers totsupers) const {
+  const auto gbxindexes_for_newsupers =
+      remove_superdrops_from_gridboxes(gbxmaps, d_gbxs, coord3lim);
+
+  auto newsupers_for_gridboxes = create_newsupers_for_gridboxes(
+      gbxmaps, create_superdrop, gbxindexes_for_newsupers, newnsupers);
+
+  add_superdrops_for_gridboxes(totsupers, d_gbxs, newsupers_for_gridboxes);
+
+  move_supers_between_gridboxes_again(d_gbxs, totsupers);  // resort totsupers view and set gbx refs
+}
+
 /* (re)sorting supers based on their gbxindexes and then updating the span for each
 gridbox accordingly.
-Kokkos::parallel_for([...]) (on host) is equivalent to:
-for (size_t ii(0); ii < ngbxs; ++ii){[...]}
-when in serial */
+Kokkos::parallel_for([...]) is equivalent in serial to:
+for (size_t ii(0); ii < d_gbxs.extent(0); ++ii){[...]}.
+*/
 void move_supers_between_gridboxes_again(const viewd_gbx d_gbxs, const viewd_supers totsupers) {
   sort_supers(totsupers);
 
@@ -41,63 +71,179 @@ void move_supers_between_gridboxes_again(const viewd_gbx d_gbxs, const viewd_sup
       });
 }
 
-/*
-Call to apply boundary conditions to remove and then add superdroplets to the top of the domain
-abouve coord3lim.
+/* set super-droplet sdgbxindex to out of bounds value if superdrop coord3 > coord3lim.
+Kokkos::parallel_for([...]) is equivalent in serial to:
+for (size_t kk(0); kk < supers.extent(0); ++kk){[...]}.
+*/
+KOKKOS_FUNCTION
+void remove_superdrop_above_coord3lim(const TeamMember &team_member, const Gridbox &gbx,
+                                      const double coord3lim) {
+  const auto supers = gbx.supersingbx();
+  Kokkos::parallel_for(
+      Kokkos::TeamThreadRange(team_member, supers.extent(0)), [supers, coord3lim](const size_t kk) {
+        if (supers(kk).get_coord3() >= coord3lim) {
+          supers(kk).set_sdgbxindex(outofbounds_gbxindex());  // remove super-droplet from domain
+        }
+      });
+}
 
-_Note:_ totsupers is view of all superdrops (both in and out of bounds of domain).
+/* for gridboxes with coordinates above coord3lim, set super-droplet sdgbxindex to out of bounds
+value if superdrop coord3 > coord3lim.
+
+Kokkos::parallel_for([...]) is equivalent in serial to:
+for (size_t ii(0); ii < d_gbxs.extent(0); ++ii){[...]}.
+
+Function returns view of all the gridboxes indexes in d_gbxs where the value of the gridbox index
+has been replaced by outofbounds_gbxindex unless superdrops were removed from that gridbox.
 */
-void AddSupersAtDomainTop::operator()(const CartesianMaps &gbxmaps, viewd_gbx d_gbxs,
-                                      const viewd_supers totsupers) const {
+Kokkos::View<unsigned int *> remove_superdrops_from_gridboxes(const CartesianMaps &gbxmaps,
+                                                              const viewd_gbx d_gbxs,
+                                                              const double coord3lim) {
   const size_t ngbxs(d_gbxs.extent(0));
+  auto gbxindexes_of_removedsupers = Kokkos::View<unsigned int *>("gbxs_of_removedsupers", ngbxs);
+  Kokkos::parallel_for(
+      "remove_superdrops", TeamPolicy(ngbxs, Kokkos::AUTO()),
+      KOKKOS_LAMBDA(const TeamMember &team_member) {
+        const int ii = team_member.league_rank();
 
-  bool is_supers_added = false;
-  for (size_t ii(0); ii < ngbxs; ++ii) {  // TODO(CB) parallelise?
-    auto &gbx(d_gbxs(ii));
-    const auto bounds = gbxmaps.coord3bounds(gbx.get_gbxindex());
-    if (bounds.second > coord3lim) {
-      remove_superdrops_from_gridbox(gbx);
-      add_superdrops_for_gridbox(gbxmaps, gbx, totsupers);
-      is_supers_added = true;
-    }
-  }
+        const auto ubound = gbxmaps.coord3bounds(d_gbxs(ii).get_gbxindex()).second;
+        if (ubound > coord3lim) {
+          remove_superdrop_above_coord3lim(team_member, d_gbxs(ii), coord3lim);
+          gbxindexes_of_removedsupers(ii) = d_gbxs(ii).get_gbxindex();  // add newsupers
+        } else {
+          gbxindexes_of_removedsupers(ii) = outofbounds_gbxindex();  // don't add newsupers
+        }
+      });
 
-  if (is_supers_added) {  // resort totsupers view and set gbx references
-    move_supers_between_gridboxes_again(d_gbxs, totsupers);
-  }
+  return gbxindexes_of_removedsupers;
 }
 
-/* set super-droplet sdgbxindex to out of bounds value if superdrop coord3 > coord3lim */
-void AddSupersAtDomainTop::remove_superdrops_from_gridbox(const Gridbox &gbx) const {
-  const auto supers = gbx.supersingbx();
-  for (size_t kk(0); kk < supers.extent(0); ++kk) {
-    if (supers(kk).get_coord3() >= coord3lim) {
-      supers(kk).set_sdgbxindex(outofbounds_gbxindex());  // remove super-droplet from domain
+/* Given a view of gridboxes where the value of the gridbox index has been replaced by
+outofbounds_gbxindex unless superdrops should be added to that gridbox,
+count the total number of new superdroplets to create.
+
+Kokkos::parallel_reduce([...]) is equivalent in serial to suming up newnsupers_total in loop:
+for (size_t ii(0); ii < d_gbxs.extent(0); ++ii){[...]}.
+*/
+size_t total_newnsupers_to_create(Kokkos::View<unsigned int *> gbxindexes,
+                                  const size_t newnsupers_pergbx) {
+  auto newnsupers_total = size_t{0};
+  Kokkos::parallel_reduce(
+      "newnsupers_total", Kokkos::RangePolicy<ExecSpace>(0, gbxindexes.extent(0)),
+      KOKKOS_LAMBDA(const size_t ii, size_t &nsupers) {
+        if (gbxindexes(ii) != outofbounds_gbxindex()) {
+          nsupers = newnsupers_pergbx;
+        } else {
+          nsupers = 0;
+        }
+      },
+      newnsupers_total);
+
+  return newnsupers_total;
+}
+
+/* Given a view of gridboxes where the value of the gridbox index has been replaced by
+outofbounds_gbxindex unless superdrops should be added to that gridbox, function create 'newnsupers'
+new superdroplets per gridbox by calling the create_superdrop function on host and then
+copies resultant view to device memory.
+*/
+viewd_supers create_newsupers_for_gridboxes(const CartesianMaps &gbxmaps,
+                                            const CreateSuperdrop &create_superdrop,
+                                            Kokkos::View<unsigned int *> gbxindexes,
+                                            const size_t newnsupers_pergbx) {
+  viewd_supers newsupers("newsupers", total_newnsupers_to_create(gbxindexes, newnsupers_pergbx));
+  auto h_newsupers = Kokkos::create_mirror_view(newsupers);
+
+  auto h_gbxindexes = Kokkos::create_mirror_view(gbxindexes);
+  Kokkos::deep_copy(h_gbxindexes, gbxindexes);
+
+  auto nn = size_t{0};  // number of super_droplets created
+  for (size_t ii(0); ii < h_gbxindexes.extent(0); ++ii) {
+    if (h_gbxindexes(ii) != outofbounds_gbxindex()) {
+      for (size_t kk(0); kk < newnsupers_pergbx; ++kk) {
+        h_newsupers(nn) = create_superdrop(gbxmaps, h_gbxindexes(ii));
+        ++nn;
+      }
     }
   }
+  Kokkos::deep_copy(newsupers, h_newsupers);
+
+  assert((newsupers.extent(0) == nn) &&
+         "total number of superdrops created must equal newsupers view size");
+
+  return newsupers;
 }
 
-/* create 'newnsupers' number of new superdroplets from the create_superdrop function */
-void AddSupersAtDomainTop::add_superdrops_for_gridbox(const CartesianMaps &gbxmaps,
-                                                      const Gridbox &gbx,
-                                                      const viewd_supers totsupers) const {
-  const auto gbxindex = gbx.get_gbxindex();
-  const size_t start = gbx.domain_totnsupers();
+/* returns copy of 1 gridbox within a view on host memory of the ii'th gridbox in a
+device view 'd_gbxs' */
+viewh_constgbx hostcopy_one_gridbox(const viewd_constgbx d_gbxs, const size_t ii) {
+  const auto d_gbx = viewd_gbx("gbx", 1);
+  Kokkos::parallel_for(
+      "copy_gbx", 1, KOKKOS_LAMBDA(const unsigned int i) { d_gbx(0) = d_gbxs(ii); });
+
+  auto h_gbx = Kokkos::create_mirror_view(d_gbx);
+  Kokkos::deep_copy(h_gbx, d_gbx);
+  return h_gbx;
+}
 
-  if (start + newnsupers > totsupers.extent(0)) {
+/* throws error if the size of newnsupers + oldnsupers > total space in totsupers view */
+size_t check_space_in_totsupers(const viewd_constsupers totsupers, const viewd_constgbx d_gbxs,
+                                const viewd_constsupers newsupers) {
+  auto h_gbx = hostcopy_one_gridbox(d_gbxs, 0);
+  const auto oldnsupers = h_gbx(0).domain_totnsupers();
+  if (oldnsupers + newsupers.extent(0) > totsupers.extent(0)) {
     const auto err = std::string(
         "UNDEFINED BEHAVIOUR: Number of super-droplets in the domain cannot become larger than the "
         "size of the super-droplets' view");
     throw std::invalid_argument(err);
   }
 
-  for (size_t kk(start); kk < start + newnsupers; ++kk) {
-    totsupers(kk) = create_superdrop(gbxmaps, gbxindex);
+  return oldnsupers;
+}
+
+/* check there is space in totsupers for newsupers, then append superdrops in newsupers to end of
+totsupers view */
+void add_superdrops_for_gridboxes(const viewd_supers totsupers, const viewd_constgbx d_gbxs,
+                                  const viewd_constsupers newsupers) {
+  auto og_totnsupers = check_space_in_totsupers(totsupers, d_gbxs, newsupers);
+
+  Kokkos::parallel_for(
+      "add_superdrops", Kokkos::RangePolicy<ExecSpace>(0, newsupers.extent(0)),
+      KOKKOS_LAMBDA(const size_t kk) { totsupers(kk + og_totnsupers) = newsupers(kk); });
+}
+
+/* returns host copy of {lower, upper} coord3 boundaries from gbxmaps for 'gbxindex' on device */
+Kokkos::pair<double, double> hostcopy_coord3bounds(const CartesianMaps &gbxmaps,
+                                                   const unsigned int gbxindex) {
+  const Kokkos::View<Kokkos::pair<double, double>[1]> d_bound("d_bound");
+  Kokkos::parallel_for(
+      "copy_coord3bound", 1,
+      KOKKOS_LAMBDA(const unsigned int i) { d_bound(0) = gbxmaps.coord3bounds(gbxindex); });
+
+  auto h_bound = Kokkos::create_mirror_view(d_bound);
+  Kokkos::deep_copy(h_bound, d_bound);
+  return h_bound(0);
+}
+
+/* call to create a new superdroplet for gridbox with given gbxindex */
+CreateSuperdrop::CreateSuperdrop(const OptionalConfigParams::AddSupersAtDomainTopParams &config)
+    : randgen(std::make_shared<std::mt19937>(std::random_device {}())),
+      sdIdGen(std::make_shared<Superdrop::IDType::Gen>(config.initnsupers)),
+      nbins(config.newnsupers),
+      log10redges(),
+      dryradius(config.DRYRADIUS / dlc::R0),
+      dist(config) {
+  const auto log10rmin = std::log10(config.MINRADIUS / dlc::R0);
+  const auto log10rmax = std::log10(config.MAXRADIUS / dlc::R0);
+  const auto log10deltar = double{(log10rmax - log10rmin) / nbins};
+  for (size_t nn(0); nn < nbins + 1; ++nn) {
+    log10redges.push_back(log10rmin + nn * log10deltar);
   }
 }
 
 /* call to create a new superdroplet for gridbox with given gbxindex */
-Superdrop CreateSuperdrop::operator()(const CartesianMaps &gbxmaps, const auto gbxindex) const {
+Superdrop CreateSuperdrop::operator()(const CartesianMaps &gbxmaps,
+                                      const unsigned int gbxindex) const {
   const auto sdgbxindex = gbxindex;
   const auto coords312 = create_superdrop_coords(gbxmaps, gbxindex);
   const auto attrs = create_superdrop_attrs(gbxmaps.get_gbxvolume(gbxindex));
@@ -109,9 +255,11 @@ Superdrop CreateSuperdrop::operator()(const CartesianMaps &gbxmaps, const auto g
 /* create spatial coordinates for super-droplet by setting coord1 = coord2 = 0.0 and coord3 to a
 random value within the gridbox's bounds */
 std::array<double, 3> CreateSuperdrop::create_superdrop_coords(const CartesianMaps &gbxmaps,
-                                                               const auto gbxindex) const {
-  const auto bounds = gbxmaps.coord3bounds(gbxindex);
-  const auto coord3 = randgen->drand(bounds.first, bounds.second);
+                                                               const unsigned int gbxindex) const {
+  const auto bounds = hostcopy_coord3bounds(gbxmaps, gbxindex);
+  auto dist = std::uniform_real_distribution<double>(bounds.first, bounds.second);
+  const double coord3 = dist(*randgen);
+
   const auto coord1 = double{0.0 / dlc::COORD0};
   const auto coord2 = double{0.0 / dlc::COORD0};
 
@@ -129,12 +277,15 @@ SuperdropAttrs CreateSuperdrop::create_superdrop_attrs(const double gbxvolume) c
 
 /* returns radius and xi for a new super-droplet by randomly sampling a distribution. */
 std::pair<size_t, double> CreateSuperdrop::new_xi_radius(const double gbxvolume) const {
-  const auto bin = uint64_t{randgen->urand(0, nbins)};  // index of randomly selected log10(r) bin
+  auto uintdist = std::uniform_int_distribution<uint64_t>(0, nbins - 1);
+  const uint64_t bin = uintdist(*randgen);  // index of randomly selected log10(r) bin
 
   const auto log10rlow = log10redges.at(bin);     // lower bound of log10(r)
   const auto log10rup = log10redges.at(bin + 1);  // upper bound of log10(r)
   const auto log10rwidth = (log10rup - log10rlow);
-  const auto frac = randgen->drand(0.0, 1.0);
+  auto dbldist = std::uniform_real_distribution<double>(0.0, 1.0);
+  const auto frac = dbldist(*randgen);
+
   const auto log10r = double{log10rlow + frac * log10rwidth};
   const auto radius = double{std::pow(10.0, log10r)};