Skip to content

Commit

Permalink
feat: generate initial radius from sampling of bins in logr using sin…
Browse files Browse the repository at this point in the history 
…gle random number generator
  • Loading branch information
@yoctoyotta1024
yoctoyotta1024 committed Apr 29, 2024
1 parent ba22698 commit ee74ee9
Show file tree
Hide file tree
Showing 3 changed files with 47 additions and 54 deletions.
73 changes: 36 additions & 37 deletions libs/cartesiandomain/add_supers_at_domain_top.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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);
}
}

Expand Down Expand Up @@ -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<ExecSpace> 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<double, 3> CreateSuperdrop::create_superdrop_coords(const CartesianMaps &gbxmaps,
URBG<ExecSpace> &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};

Expand All @@ -126,21 +120,22 @@ std::array<double, 3> 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<size_t, double> 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;
Expand All @@ -149,8 +144,12 @@ std::pair<size_t, double> 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;
}
26 changes: 10 additions & 16 deletions libs/cartesiandomain/add_supers_at_domain_top.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,12 +24,13 @@
#define LIBS_CARTESIANDOMAIN_ADD_SUPERS_AT_DOMAIN_TOP_HPP_

#include <Kokkos_Core.hpp>
#include <Kokkos_Random.hpp>
#include <array>
#include <cmath>
#include <memory>
#include <pair>
#include <random>
#include <stdexcept>
#include <utility>
#include <vector>

#include "../cleoconstants.hpp"
Expand All @@ -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<Kokkos::Random_XorShift64<HostSpace>>
randgen; /**< pointer to Kokkos random number generator */
std::shared_ptr<Superdrop::IDType::Gen>
sdIdGen; /**< Pointer Superdrop::IDType object for super-droplet ID generation. */
double dryradius; /**< dry radius of new superdrop */
Expand All @@ -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<double, 3> create_superdrop_coords(const CartesianMaps &gbxmaps, URBG<ExecSpace> &urbg,
std::array<double, 3> create_superdrop_coords(const CartesianMaps &gbxmaps,
const auto gbxindex) const;

/* create attributes for a new super-droplet */
Expand All @@ -62,18 +63,18 @@ struct CreateSuperdrop {
std::pair<size_t, double> 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<Kokkos::Random_XorShift64<HostSpace>>(std::random_device {}())),
sdIdGen(std::make_shared<Superdrop::IDType::Gen>(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);
Expand All @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/config/config.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]

0 comments on commit ee74ee9

Please sign in to comment.