diff --git a/examples/eurec4a1d/src/config/eurec4a1d_config.yaml b/examples/eurec4a1d/src/config/eurec4a1d_config.yaml index 8ca132de4..05795f1bc 100644 --- a/examples/eurec4a1d/src/config/eurec4a1d_config.yaml +++ b/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 diff --git a/examples/eurec4a1d/src/main_eurec4a1D.cpp b/examples/eurec4a1d/src/main_eurec4a1D.cpp index c5bce5084..5f9ab82be 100644 --- a/examples/eurec4a1d/src/main_eurec4a1D.cpp +++ b/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 #include +#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 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); } diff --git a/libs/cartesiandomain/add_supers_at_domain_top.cpp b/libs/cartesiandomain/add_supers_at_domain_top.cpp index de3a6b87c..3679cf2ab 100644 --- a/libs/cartesiandomain/add_supers_at_domain_top.cpp +++ b/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 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 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 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("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 gbxindexes, + const size_t newnsupers_pergbx) { + auto newnsupers_total = size_t{0}; + Kokkos::parallel_reduce( + "newnsupers_total", Kokkos::RangePolicy(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 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(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 hostcopy_coord3bounds(const CartesianMaps &gbxmaps, + const unsigned int gbxindex) { + const Kokkos::View[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::random_device {}())), + sdIdGen(std::make_shared(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 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(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 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(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(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)}; diff --git a/libs/cartesiandomain/add_supers_at_domain_top.hpp b/libs/cartesiandomain/add_supers_at_domain_top.hpp index 9fa460e60..76e385a13 100644 --- a/libs/cartesiandomain/add_supers_at_domain_top.hpp +++ b/libs/cartesiandomain/add_supers_at_domain_top.hpp @@ -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 @@ -24,13 +24,13 @@ #define LIBS_CARTESIANDOMAIN_ADD_SUPERS_AT_DOMAIN_TOP_HPP_ #include -#include #include #include #include #include #include #include +#include #include #include @@ -74,8 +74,7 @@ class TwoLognormalsDistribution { struct CreateSuperdrop { private: - std::shared_ptr> - randgen; /**< pointer to Kokkos random number generator */ + std::shared_ptr randgen; /**< pointer to random number generator */ std::shared_ptr sdIdGen; /**< Pointer Superdrop::IDType object for super-droplet ID generation. */ size_t nbins; /**< number of bins for sampling superdroplet radius */ @@ -83,10 +82,10 @@ struct CreateSuperdrop { double dryradius; /**< dry radius of new superdrop */ TwoLognormalsDistribution dist; /**< distribution for creating superdroplet xi */ - /* create spatial coordinates for super-droplet by setting coord1 = coord2 = 0.0 and coord3 to a - random value within the gridbox's bounds */ + /* 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, - const auto gbxindex) const; + const unsigned int gbxindex) const; /* create attributes for a new super-droplet */ SuperdropAttrs create_superdrop_attrs(const double gbxvolume) const; @@ -99,22 +98,9 @@ struct CreateSuperdrop { public: /* call to create a new superdroplet for gridbox with given gbxindex */ - explicit CreateSuperdrop(const OptionalConfigParams::AddSupersAtDomainTopParams &config) - : randgen(std::make_shared>(std::random_device {}())), - sdIdGen(std::make_shared(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); - } - } - - Superdrop operator()(const CartesianMaps &gbxmaps, const auto gbxindex) const; + explicit CreateSuperdrop(const OptionalConfigParams::AddSupersAtDomainTopParams &config); + + Superdrop operator()(const CartesianMaps &gbxmaps, const unsigned int gbxindex) const; }; struct AddSupersAtDomainTop { @@ -123,13 +109,6 @@ struct AddSupersAtDomainTop { double coord3lim; /**< gridboxes with upper bound > coord3lim get new super-droplets */ CreateSuperdrop create_superdrop; /**< methods to create a new superdrop */ - /* set super-droplet sdgbxindex to out of bounds value if superdrop coord3 > coord3lim */ - void remove_superdrops_from_gridbox(const Gridbox &gbx) const; - - /* create 'newnsupers' number of new superdroplets from the create_superdrop function */ - void add_superdrops_for_gridbox(const CartesianMaps &gbxmaps, const Gridbox &gbx, - 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/libs/gridboxes/movesupersindomain.hpp b/libs/gridboxes/movesupersindomain.hpp index f75baab4a..6a45cf556 100644 --- a/libs/gridboxes/movesupersindomain.hpp +++ b/libs/gridboxes/movesupersindomain.hpp @@ -8,7 +8,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 @@ -34,113 +34,101 @@ #include "superdrops/motion.hpp" #include "superdrops/superdrop.hpp" -/* struct for functionality to move superdroplets throughtout +/* +struct for functionality to move superdroplets throughtout the domain by updating their spatial coordinates (according to some type of Motion) and then moving them between gridboxes -after updating their gridbox indexes concordantly */ +after updating their gridbox indexes concordantly +*/ template M, typename BoundaryConditions> struct MoveSupersInDomain { - M motion; - BoundaryConditions apply_domain_boundary_conditions; - - /* enact steps (1) and (2) movement of superdroplets for 1 gridbox: - (1) update their spatial coords according to type of motion. (device) - (1b) optional detect precipitation (device) - (2) update their sdgbxindex accordingly (device). - Kokkos::parallel_for([...]) is equivalent to: - for (size_t kk(0); kk < supers.extent(0); ++kk) {[...]} - when in serial */ - KOKKOS_INLINE_FUNCTION - void move_supers_in_gbx(const TeamMember &team_member, const unsigned int gbxindex, - const GbxMaps &gbxmaps, const State &state, - const subviewd_supers supers) const { - const size_t nsupers(supers.extent(0)); - Kokkos::parallel_for(Kokkos::TeamThreadRange(team_member, nsupers), [&, this](const size_t kk) { - /* step (1) */ - motion.superdrop_coords(gbxindex, gbxmaps, state, supers(kk)); - - /* optional step (1b) */ - // gbx.detectors -> detect_precipitation(area, drop); // TODO(CB) detectors - - /* step (2) */ - motion.superdrop_gbx(gbxindex, gbxmaps, supers(kk)); - }); - } - - /* enact steps (1) and (2) movement of superdroplets - throughout domain (i.e. for all gridboxes): - (1) update their spatial coords according to type of motion. (device) - (1b) optional detect precipitation (device) - (2) update their sdgbxindex accordingly (device). - Kokkos::parallel_for([...]) is equivalent to: - for (size_t ii(0); ii < ngbxs; ++ii) {[...]} - when in serial */ - void move_supers_in_gridboxes(const GbxMaps &gbxmaps, const viewd_gbx d_gbxs) const { - const size_t ngbxs(d_gbxs.extent(0)); - - Kokkos::parallel_for( - "move_supers_in_gridboxes", TeamPolicy(ngbxs, Kokkos::AUTO()), - KOKKOS_CLASS_LAMBDA(const TeamMember &team_member) { - const int ii = team_member.league_rank(); - - auto &gbx(d_gbxs(ii)); - move_supers_in_gbx(team_member, gbx.get_gbxindex(), gbxmaps, gbx.state, - gbx.supersingbx()); - }); - } - - /* (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 - _Note:_ totsupers is view of all superdrops (both in and out of bounds of domain). - */ - void 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); - }); - - // /* optional (expensive!) test to raise error if - // superdrops' gbxindex doesn't match gridbox's gbxindex */ - // for (size_t ii(0); ii < ngbxs; ++ii) - // { - // d_gbxs(ii).supersingbx.iscorrect(); - // } - } - - /* enact movement of superdroplets throughout domain in three stages: - (1) update their spatial coords according to type of motion. (device) - (1b) optional detect precipitation (device) - (2) update their sdgbxindex accordingly (device) - (3) move superdroplets between gridboxes (host) - (4) (optional) apply domain boundary conditions (host and device) - _Note:_ totsupers is view of all superdrops (both in and out of bounds of domain). - // TODO(all) use tasking to convert all 3 team policy - // loops from first two function calls into 1 loop? + /* + EnactMotion struct encapsulates motion so that parallel loops with KOKKOS_CLASS_LAMBDA + (ie. [=] on CPUs) functors only captures motion and not other members of MoveSupersInDomain + coincidentally (which may not be GPU compatible). */ - void move_superdrops_in_domain(const unsigned int t_sdm, const GbxMaps &gbxmaps, viewd_gbx d_gbxs, - const viewd_supers totsupers) const { - /* steps (1 - 2) */ - move_supers_in_gridboxes(gbxmaps, d_gbxs); + struct EnactMotion { + M motion; + + /* enact steps (1) and (2) movement of superdroplets for 1 gridbox: + (1) update their spatial coords according to type of motion. (device) + (1b) optional detect precipitation (device) + (2) update their sdgbxindex accordingly (device). + Kokkos::parallel_for([...]) is equivalent to: + for (size_t kk(0); kk < supers.extent(0); ++kk) {[...]} + when in serial */ + KOKKOS_INLINE_FUNCTION + void move_supers_in_gbx(const TeamMember &team_member, const unsigned int gbxindex, + const GbxMaps &gbxmaps, const State &state, + const subviewd_supers supers) const { + const size_t nsupers(supers.extent(0)); + Kokkos::parallel_for(Kokkos::TeamThreadRange(team_member, nsupers), + [&, this](const size_t kk) { + /* step (1) */ + motion.superdrop_coords(gbxindex, gbxmaps, state, supers(kk)); + + /* optional step (1b) */ + // gbx.detectors -> detect_precipitation(area, drop); // TODO(CB) + // detectors + + /* step (2) */ + motion.superdrop_gbx(gbxindex, gbxmaps, supers(kk)); + }); + } - /* step (3) */ - move_supers_between_gridboxes(d_gbxs, totsupers); + /* enact steps (1) and (2) movement of superdroplets + throughout domain (i.e. for all gridboxes): + (1) update their spatial coords according to type of motion. (device) + (1b) optional detect precipitation (device) + (2) update their sdgbxindex accordingly (device). + Kokkos::parallel_for([...]) is equivalent to: + for (size_t ii(0); ii < ngbxs; ++ii) {[...]} + when in serial */ + void move_supers_in_gridboxes(const GbxMaps &gbxmaps, const viewd_gbx d_gbxs) const { + const size_t ngbxs(d_gbxs.extent(0)); + + Kokkos::parallel_for( + "move_supers_in_gridboxes", TeamPolicy(ngbxs, Kokkos::AUTO()), + KOKKOS_CLASS_LAMBDA(const TeamMember &team_member) { + const int ii = team_member.league_rank(); + + auto &gbx(d_gbxs(ii)); + move_supers_in_gbx(team_member, gbx.get_gbxindex(), gbxmaps, gbx.state, + gbx.supersingbx()); + }); + } - /* step (4) */ - apply_domain_boundary_conditions(gbxmaps, d_gbxs, totsupers); - } + /* (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 + _Note:_ totsupers is view of all superdrops (both in and out of bounds of domain). + */ + void 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); + }); + + // /* optional (expensive!) test to raise error if + // superdrops' gbxindex doesn't match gridbox's gbxindex */ + // for (size_t ii(0); ii < ngbxs; ++ii) + // { + // d_gbxs(ii).supersingbx.iscorrect(); + // } + } + } enactmotion; MoveSupersInDomain(const M mtn, const BoundaryConditions boundary_conditions) - : motion(mtn), apply_domain_boundary_conditions(boundary_conditions) {} + : enactmotion({mtn}), apply_domain_boundary_conditions(boundary_conditions) {} /* extra constructor useful to help when compiler cannot deduce type of GBxMaps */ MoveSupersInDomain(const GbxMaps &gbxmaps, const M mtn, @@ -150,7 +138,9 @@ struct MoveSupersInDomain { /* returns time when superdroplet motion is next due to occur given current time, t_sdm */ KOKKOS_INLINE_FUNCTION - unsigned int next_step(const unsigned int t_sdm) const { return motion.next_step(t_sdm); } + unsigned int next_step(const unsigned int t_sdm) const { + return enactmotion.motion.next_step(t_sdm); + } /* * if current time, t_sdm, is time when superdrop motion should occur, enact movement of @@ -161,10 +151,35 @@ struct MoveSupersInDomain { */ void run_step(const unsigned int t_sdm, const GbxMaps &gbxmaps, viewd_gbx d_gbxs, const viewd_supers totsupers) const { - if (motion.on_step(t_sdm)) { + if (enactmotion.motion.on_step(t_sdm)) { move_superdrops_in_domain(t_sdm, gbxmaps, d_gbxs, totsupers); } } + + private: + BoundaryConditions apply_domain_boundary_conditions; + + /* enact movement of superdroplets throughout domain in three stages: + (1) update their spatial coords according to type of motion. (device) + (1b) optional detect precipitation (device) + (2) update their sdgbxindex accordingly (device) + (3) move superdroplets between gridboxes (host) + (4) (optional) apply domain boundary conditions (host and device) + _Note:_ totsupers is view of all superdrops (both in and out of bounds of domain). + // TODO(all) use tasking to convert all 3 team policy + // loops from first two function calls into 1 loop? + */ + void move_superdrops_in_domain(const unsigned int t_sdm, const GbxMaps &gbxmaps, viewd_gbx d_gbxs, + const viewd_supers totsupers) const { + /* steps (1 - 2) */ + enactmotion.move_supers_in_gridboxes(gbxmaps, d_gbxs); + + /* step (3) */ + enactmotion.move_supers_between_gridboxes(d_gbxs, totsupers); + + /* step (4) */ + apply_domain_boundary_conditions(gbxmaps, d_gbxs, totsupers); + } }; #endif // LIBS_GRIDBOXES_MOVESUPERSINDOMAIN_HPP_ diff --git a/src/config/config.yaml b/src/config/config.yaml index 5a9e481ce..317e833d8 100644 --- a/src/config/config.yaml +++ b/src/config/config.yaml @@ -32,9 +32,9 @@ timesteps: CONDTSTEP : 0.1 # time between SD condensation [s] COLLTSTEP : 1 # time between SD collision [s] MOTIONTSTEP : 2 # time between SDM motion [s] - COUPLTSTEP : 4 # time between dynamic couplings [s] + COUPLTSTEP : 8 # time between dynamic couplings [s] OBSTSTEP : 2 # time between SDM observations [s] - T_END : 4 # time span of integration from 0s to T_END [s] + T_END : 16 # time span of integration from 0s to T_END [s] ### Initialisation Parameters ### inputfiles: diff --git a/src/main_impl.hpp b/src/main_impl.hpp index 0f2bfad0b..f922744d6 100644 --- a/src/main_impl.hpp +++ b/src/main_impl.hpp @@ -9,7 +9,7 @@ * Author: Clara Bayley (CB) * Additional Contributors: * ----- - * Last Modified: Friday 3rd May 2024 + * Last Modified: Saturday 4th May 2024 * Modified By: CB * ----- * License: BSD 3-Clause "New" or "Revised" License @@ -112,8 +112,8 @@ inline Motion auto create_motion(const unsigned int motionstep) { } inline auto create_boundary_conditions(const Config &config) { - // return AddSupersAtDomainTop(config.get_addsupersatdomaintop()); - return NullBoundaryConditions{}; + return AddSupersAtDomainTop(config.get_addsupersatdomaintop()); + // return NullBoundaryConditions{}; } template