diff --git a/libs/cartesiandomain/add_supers_at_domain_top.cpp b/libs/cartesiandomain/add_supers_at_domain_top.cpp index 49fe8da34..29489aef2 100644 --- a/libs/cartesiandomain/add_supers_at_domain_top.cpp +++ b/libs/cartesiandomain/add_supers_at_domain_top.cpp @@ -22,6 +22,25 @@ #include "cartesiandomain/add_supers_at_domain_top.hpp" +/* (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 */ +void move_supers_between_gridboxes_again(const viewd_gbx d_gbxs, const viewd_supers totsupers) { + sort_supers(totsupers); + + const size_t ngbxs(d_gbxs.extent(0)); + Kokkos::parallel_for( + "move_supers_between_gridboxes_again", TeamPolicy(ngbxs, Kokkos::AUTO()), + KOKKOS_LAMBDA(const TeamMember &team_member) { + const int ii = team_member.league_rank(); + + auto &gbx(d_gbxs(ii)); + gbx.supersingbx.set_refs(team_member); + }); +} + /* Call to apply boundary conditions to remove and then add superdroplets to the top of the domain abouve coord3lim. @@ -44,7 +63,7 @@ void AddSupersAtDomainTop::operator()(const CartesianMaps &gbxmaps, viewd_gbx d_ } if (is_supers_added) { // resort totsupers view and set gbx references - move_supers_between_gridboxes(d_gbxs, totsupers); + move_supers_between_gridboxes_again(d_gbxs, totsupers); } } @@ -77,47 +96,22 @@ void AddSupersAtDomainTop::add_superdrops_for_gridbox(const CartesianMaps &gbxma } } -/* (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 */ -void AddSupersAtDomainTop::move_supers_between_gridboxes(const viewd_gbx d_gbxs, - const viewd_supers totsupers) const { - sort_supers(totsupers); - - const size_t ngbxs(d_gbxs.extent(0)); - Kokkos::parallel_for( - "move_supers_between_gridboxes", TeamPolicy(ngbxs, Kokkos::AUTO()), - KOKKOS_CLASS_LAMBDA(const TeamMember &team_member) { - const int ii = team_member.league_rank(); - - auto &gbx(d_gbxs(ii)); - gbx.supersingbx.set_refs(team_member); - }); -} - /* call to create a new superdroplet for gridbox with given gbxindex */ Superdrop CreateSuperdrop::operator()(const CartesianMaps &gbxmaps, const auto gbxindex) const { - URBG urbg{genpool4reset.get_state()}; // thread safe random number generator - const auto sdgbxindex = gbxindex; - const auto coords312 = create_superdrop_coords(gbxmaps, urbg, gbxindex); + const auto coords312 = create_superdrop_coords(gbxmaps, gbxindex); const auto attrs = create_superdrop_attrs(); const auto sdId = sdIdGen->next(); - genpool4reset.free_state(urbg.gen); - return Superdrop(sdgbxindex, coords312[0], coords312[1], coords312[2], attrs, sdId); } /* 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 CreateSuperdrop::create_superdrop_coords(const CartesianMaps &gbxmaps, - URBG &urbg, const auto gbxindex) const { const auto bounds = gbxmaps.coord3bounds(gbxindex); - const auto coord3 = urbg.drand(bounds.first, bounds.second); + const auto coord3 = randgen->drand(bounds.first, bounds.second); const auto coord1 = double{0.0 / dlc::COORD0}; const auto coord2 = double{0.0 / dlc::COORD0}; @@ -126,21 +120,22 @@ std::array CreateSuperdrop::create_superdrop_coords(const CartesianMa /* create attributes for a new super-droplet */ SuperdropAttrs CreateSuperdrop::create_superdrop_attrs() const { - const auto xi_radius = new_xi_radius(); - const auto msol = new_msol(); + const auto [xi, radius] = new_xi_radius(); + const auto msol = new_msol(radius); const auto solute = SoluteProperties{}; - return SuperdropAttrs(solute, xi_radius.first, xi_radius.second, msol); + return SuperdropAttrs(solute, xi, radius, msol); } /* returns radius and xi for a new super-droplet by randomly sampling a distribution. */ std::pair CreateSuperdrop::new_xi_radius() const { - const auto bin = uint64_t{urbg(0, nbins)}; // index of randomly selected log10(r) bin - const auto log10rlow = log10redges(bin); // lower bound of log10(r) - const auto log10rup = log10redges(bin + 1); // upper bound of log10(r) + const auto bin = uint64_t{randgen->urand(0, nbins)}; // index of randomly selected log10(r) bin - const auto frac = urbg.drand(0.0, 1.0); - const auto log10r = double{log10rlow + frac * (log10rup - log10rlow)}; + 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); + const auto log10r = double{log10rlow + frac * log10rwidth}; const auto radius = double{Kokkos::pow(10.0, log10r)}; const auto xi = 100; @@ -149,8 +144,12 @@ std::pair CreateSuperdrop::new_xi_radius() const { } /* returns solute mass for a new super-droplet with a dryradius = 1nano-meter. */ -double CreateSuperdrop::new_msol() const { +double CreateSuperdrop::new_msol(const double radius) const { constexpr double msolconst = 4.0 * Kokkos::numbers::pi * dlc::Rho_sol / 3.0; + if (radius < dryradius) { + throw std::invalid_argument("new radius cannot be < dry radius of droplet"); + } + return msolconst * dryradius * dryradius * dryradius; } diff --git a/libs/cartesiandomain/add_supers_at_domain_top.hpp b/libs/cartesiandomain/add_supers_at_domain_top.hpp index 4db84a4bd..09c13597b 100644 --- a/libs/cartesiandomain/add_supers_at_domain_top.hpp +++ b/libs/cartesiandomain/add_supers_at_domain_top.hpp @@ -24,12 +24,13 @@ #define LIBS_CARTESIANDOMAIN_ADD_SUPERS_AT_DOMAIN_TOP_HPP_ #include +#include #include #include #include -#include #include #include +#include #include #include "../cleoconstants.hpp" @@ -39,11 +40,11 @@ #include "gridboxes/sortsupers.hpp" #include "initialise/optional_config_params.hpp" #include "superdrops/superdrop.hpp" -#include "superdrops/urbg.hpp" struct CreateSuperdrop { private: - GenRandomPool genpool4reset; /**< Kokkos pool for random number generation */ + std::shared_ptr> + randgen; /**< pointer to Kokkos random number generator */ std::shared_ptr sdIdGen; /**< Pointer Superdrop::IDType object for super-droplet ID generation. */ double dryradius; /**< dry radius of new superdrop */ @@ -52,7 +53,7 @@ struct CreateSuperdrop { /* 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 create_superdrop_coords(const CartesianMaps &gbxmaps, URBG &urbg, + std::array create_superdrop_coords(const CartesianMaps &gbxmaps, const auto gbxindex) const; /* create attributes for a new super-droplet */ @@ -62,18 +63,18 @@ struct CreateSuperdrop { std::pair new_xi_radius() const; /* returns solute mass for a new super-droplet with a dryradius = 1nano-meter. */ - double new_msol() const; + double new_msol(const double radius) const; public: /* call to create a new superdroplet for gridbox with given gbxindex */ explicit CreateSuperdrop(const OptionalConfigParams::AddSupersAtDomainTopParams &config) - : genpool4reset(std::random_device {}()), + : randgen(std::make_shared>(std::random_device {}())), sdIdGen(std::make_shared(config.initnsupers)), - dryradius(config.dryradius / dlc::R0), + dryradius(config.DRYRADIUS / dlc::R0), nbins(config.newnsupers), log10redges() { - const auto log10rmin = std::log10(config.minradius / dlc::R0); - const auto log10rmax = std::log10(config.maxradius / dlc::R0); + 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); @@ -96,13 +97,6 @@ struct AddSupersAtDomainTop { void add_superdrops_for_gridbox(const CartesianMaps &gbxmaps, const Gridbox &gbx, const viewd_supers totsupers) const; - /* (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 */ - void move_supers_between_gridboxes(const viewd_gbx d_gbxs, const viewd_supers totsupers) const; - public: /* New super-droplets are added to domain with coord3 >= COORD3LIM [m]. Note generation of * nextsdId assumes it is the only method creating super-droplets - otherwise created sdId may not diff --git a/src/config/config.yaml b/src/config/config.yaml index 0759ee914..6729856ee 100644 --- a/src/config/config.yaml +++ b/src/config/config.yaml @@ -79,5 +79,5 @@ boundary_conditions: COORD3LIM: 805 # SDs added to domain with coord3 >= COORD3LIM [m] newnsupers: 256 # number SDs to add to each gridbox above COORD3LIM DRYRADIUS: 1e-9 # dry radius of new super-droplets (for solute mass) [m] - MINRADIUS: 1e-7 # minimum radius of new super-droplets [m] + MINRADIUS: 1e-6 # minimum radius of new super-droplets [m] MAXRADIUS: 1e-3 # maximum radius of new super-droplets [m]