Skip to content

Refactor distribution-specific configuration #503

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
pdziekan wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Refactor distribution-specific configuration #503

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
5749604
distribution-specific sd_conc
pdziekan Sep 15, 2026
efe0f15
dry distros-specific const_multi
pdziekan Sep 15, 2026
66fc562
src dry sizes: supstp in key
pdziekan Sep 15, 2026
8f9a130
count_num version of src_dry_distros
pdziekan Sep 15, 2026
abcb427
test tweaks after init distro changes
pdziekan Sep 15, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
51 changes: 26 additions & 25 deletions bindings/python/lgrngn.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -234,46 +234,47 @@ namespace libcloudphxx
return *arg->opts_init;
}

// set dry distros from a dict with (kappa, soluble_fraction) as key
// set dry distros from a dict with (kappa, soluble_fraction, sd_conc, sd_const_multi) as key
template <typename real_t>
void set_dd(
lgr::opts_init_t<real_t> *arg,
const bp::dict &kappa_func) // a dict keyed by (kappa, soluble_fraction)
const bp::dict &kappa_func) // a dict keyed by (kappa, soluble_fraction, sd_conc, sd_const_multi)
{
arg->dry_distros.clear();
for (int i = 0; i < len(kappa_func.keys()); ++i)
{
bp::tuple key = bp::extract<bp::tuple>(kappa_func.keys()[i]);
const real_t kappa = bp::extract<real_t>(key[0]);
const real_t soluble_fraction = bp::extract<real_t>(key[1]);
const unsigned long long sd_conc = bp::extract<unsigned long long>(key[2]);
const unsigned long long sd_const_multi = bp::extract<unsigned long long>(key[3]);
arg->dry_distros.emplace(
libcloudphxx::lgrngn::kappa_soluble_fraction_t<real_t>{kappa, soluble_fraction},
std::make_shared<detail::pyunary<real_t>>(kappa_func.values()[i])
std::make_tuple(kappa, soluble_fraction, sd_conc, sd_const_multi),
std::static_pointer_cast<libcloudphxx::common::unary_function<real_t>>(
std::make_shared<detail::pyunary<real_t>>(kappa_func.values()[i]))
);
}
}

// src_dry_distros moved from opts_init to opts
// src_dry_distros moved from opts_init to opts; key is (kappa, soluble_fraction, sd_conc, sd_const_multi, supstp)
template <typename real_t>
void set_sdd( // src_dry_distro
lgr::opts_t<real_t> *arg,
const bp::dict &kappa_func) // a dict keyed by (kappa, soluble_fraction)
const bp::dict &kappa_func) // a dict keyed by (kappa, soluble_fraction, sd_conc, sd_const_multi, supstp)
{
arg->src_dry_distros.clear();
for (int i = 0; i < len(kappa_func.keys()); ++i)
{
bp::tuple key = bp::extract<bp::tuple>(kappa_func.keys()[i]);
bp::tuple val = bp::extract<bp::tuple>(kappa_func.values()[i]);
const real_t kappa = bp::extract<real_t>(key[0]);
const real_t soluble_fraction = bp::extract<real_t>(key[1]);

const int sd_conc = bp::extract<int>(val[1]);
const int supstp = bp::extract<int>(val[2]);
const int sd_conc = bp::extract<int>(key[2]);
const int sd_const_multi = bp::extract<int>(key[3]);
const int supstp = bp::extract<int>(key[4]);
arg->src_dry_distros.emplace(
libcloudphxx::lgrngn::kappa_soluble_fraction_t<real_t>{kappa, soluble_fraction},
std::make_tuple(std::static_pointer_cast<libcloudphxx::common::unary_function<real_t>>(
std::make_shared<detail::pyunary<real_t>>(val[0])),
sd_conc, supstp)
std::make_tuple(kappa, soluble_fraction, sd_conc, sd_const_multi, supstp),
std::static_pointer_cast<libcloudphxx::common::unary_function<real_t>>(
std::make_shared<detail::pyunary<real_t>>(kappa_func.values()[i]))
);
}
}
Expand Down Expand Up @@ -311,44 +312,44 @@ namespace libcloudphxx
const int count = bp::extract<int> (conc_count_list[1]);
size_conc_map[bp::extract<real_t>(size_conc.keys()[i])] = std::make_pair(conc, count);
}
arg->dry_sizes[libcloudphxx::lgrngn::kappa_soluble_fraction_t<real_t>{kappa, soluble_fraction}] = size_conc_map;
arg->dry_sizes[std::make_tuple(kappa, soluble_fraction)] = size_conc_map;
}
}

// src_dry_sizes moved from opts_init to opts
template <typename real_t>
void set_sds( // src_dry_sizes
lgr::opts_t<real_t> *arg,
const bp::dict &kappa_func // a dict keyed by (kappa, soluble_fraction)
const bp::dict &kappa_func // a dict keyed by (kappa, soluble_fraction, supstp)
)
{
arg->src_dry_sizes.clear();
if(len(kappa_func.keys()) == 0)
return;

// loop over kappas and soluble_fraction
// loop over kappa, soluble_fraction, and source intervals
for (int j = 0; j < len(kappa_func.keys()); ++j)
{
// extract the key tuple (kappa, soluble_fraction)
// extract the key tuple (kappa, soluble_fraction, supstp)
const bp::tuple key = bp::extract<bp::tuple>(kappa_func.keys()[j]);
const real_t kappa = bp::extract<real_t>(key[0]);
const real_t soluble_fraction = bp::extract<real_t>(key[1]);
const int supstp = bp::extract<int>(key[2]);

// extract size : {conc, count, supstp_src} dict for this (kappa, soluble_fraction)
// extract size : {conc, count} dict for this (kappa, soluble_fraction, supstp)
const bp::dict size_conc = bp::extract<bp::dict>(kappa_func.values()[j]);
std::map<double, std::tuple<real_t, int, int>> size_conc_map;
std::map<double, std::pair<real_t, int>> size_conc_map;

// turn the size : {conc, count, supstp_src} dict into a size : {conc, count, supstp_src} map
// turn the size : {conc, count} dict into a size : {conc, count} map
for (int i = 0; i < len(size_conc.keys()); ++i)
{
const bp::list conc_count_list = bp::extract<bp::list>(size_conc.values()[i]);
assert(len(conc_count_list) == 3);
assert(len(conc_count_list) == 2);
const real_t conc = bp::extract<real_t>(conc_count_list[0]);
const int count = bp::extract<int> (conc_count_list[1]);
const int supstp = bp::extract<int> (conc_count_list[2]);
size_conc_map[bp::extract<real_t>(size_conc.keys()[i])] = std::make_tuple(conc, count, supstp);
size_conc_map[bp::extract<real_t>(size_conc.keys()[i])] = std::make_pair(conc, count);
}
arg->src_dry_sizes[libcloudphxx::lgrngn::kappa_soluble_fraction_t<real_t>{kappa, soluble_fraction}] = size_conc_map;
arg->src_dry_sizes[std::make_tuple(kappa, soluble_fraction, supstp)] = size_conc_map;
}
}

Expand Down
2 changes: 0 additions & 2 deletions bindings/python/lib.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -331,10 +331,8 @@ BOOST_PYTHON_MODULE(libcloudphxx)
.def_readwrite("supstp_rlx", &lgr::opts_init_t<real_t>::supstp_rlx)
.def_readwrite("kernel", &lgr::opts_init_t<real_t>::kernel)
.def_readwrite("adve_scheme", &lgr::opts_init_t<real_t>::adve_scheme)
.def_readwrite("sd_conc", &lgr::opts_init_t<real_t>::sd_conc)
.def_readwrite("sd_conc_large_tail", &lgr::opts_init_t<real_t>::sd_conc_large_tail)
.def_readwrite("aerosol_independent_of_rhod", &lgr::opts_init_t<real_t>::aerosol_independent_of_rhod)
.def_readwrite("sd_const_multi", &lgr::opts_init_t<real_t>::sd_const_multi)
.def_readwrite("rlx_bins", &lgr::opts_init_t<real_t>::rlx_bins)
.def_readwrite("rlx_sd_per_bin", &lgr::opts_init_t<real_t>::rlx_sd_per_bin)
.def_readwrite("rlx_timescale", &lgr::opts_init_t<real_t>::rlx_timescale)
Expand Down
4 changes: 1 addition & 3 deletions docs/API.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,6 @@ opts_init.nz = 64; // grid cells in z
opts_init.dx = 100.0; // cell size x [m]
opts_init.dy = 100.0; // cell size y [m]
opts_init.dz = 100.0; // cell size z [m]
opts_init.sd_conc = 64; // super-droplets per cell

// Create particle system
auto particles = libcloudphxx::lgrngn::factory<double>(
Expand Down Expand Up @@ -439,7 +438,7 @@ arrinfo_t<double> info_3d_v(arr_3d, str);

See [USER_OPTIONS.md](USER_OPTIONS.md#lagrangian-lgrngn-initialization-options) for complete documentation of initialization options including:
- Domain configuration (nx, ny, nz, dx, dy, dz)
- Super-droplet configuration (sd_conc, sd_conc_mean)
- Super-droplet configuration (per-distribution sd_conc, sd_const_multi)
- Aerosol size distributions (dry_distros, dry_sizes)
- Chemistry options (chem_switch, chem_rho)
- GPU settings (dev_count, dev_id)
Expand Down Expand Up @@ -954,7 +953,6 @@ int main() {
opts_init.dt = dt;
opts_init.nx = nx; opts_init.ny = ny; opts_init.nz = nz;
opts_init.dx = dx; opts_init.dy = dy; opts_init.dz = dz;
opts_init.sd_conc = 64;
opts_init.dry_distros = {{0.04e-6, 60e6}}; // 40 nm mode
opts_init.kappa = 0.61;

Expand Down
16 changes: 6 additions & 10 deletions docs/USER_OPTIONS.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -135,8 +135,7 @@ These options are set once at initialization and cannot be changed during the si

| Option | Type | Default | Description |
|--------|------|---------|-------------|
| `sd_conc` | `unsigned long long` | `0` | Number of super-droplets per cell |
| `sd_const_multi` | `unsigned long long` | `0` | Alternative to `sd_conc`: constant multiplicity for all SDs |
| `dry_distros[].sd_const_multi` | `unsigned long long` | `0` | Constant multiplicity for the distribution; cannot be used with that distribution's `sd_conc` |
| `n_sd_max` | `unsigned long long` | `0` | Maximum number of super-droplets in the system (should account for sources) |
| `sd_conc_large_tail` | `bool` | `false` | Add more SDs to better represent large tail of the distribution |
| `rd_min`, `rd_max` | `real_t` | `-1` | Min/max dry radius of droplets [m]; negative = auto-detect |
Expand All @@ -147,8 +146,8 @@ Two methods are available:

**1. Distribution-based (recommended):**
```cpp
typedef std::unordered_map<
kappa_soluble_fraction_t<real_t>, // (kappa - hygroscopicity parameter, soluble_fraction - volume fraction of soluble part)
typedef std::map<
std::tuple<real_t, real_t, unsigned long long, unsigned long long>, // kappa, soluble_fraction, sd_conc, sd_const_multi
std::shared_ptr<unary_function<real_t>> // n(ln(rd)) @ STP
> dry_distros_t;
dry_distros_t dry_distros;
Expand All @@ -157,7 +156,7 @@ dry_distros_t dry_distros;
**2. Size-number pairs:**
```cpp
typedef std::map<
kappa_soluble_fraction_t<real_t> // (kappa - hygroscopicity parameter, soluble_fraction - volume fraction of soluble part)
std::tuple<real_t, real_t>, // kappa, soluble_fraction
std::map<real_t, // radius [m]
std::pair<real_t, int> // STP concentration [1/m^3], number of SDs
>
Expand Down Expand Up @@ -275,13 +274,11 @@ dry_sizes_t dry_sizes;
| Option | Type | Default | Description |
|--------|------|---------|-------------|
| `src_type` | `src_t` | `off` | Type of CCN source (`off`, `constant`, etc.) |
| `src_dry_distros` | `dry_distros_t` | - | Source distribution per unit time |
| `src_dry_distros` | `src_dry_distros_t` | - | Source distribution per unit time; map key is `(kappa, soluble_fraction, sd_conc, sd_const_multi, supstp)` and simple sources require exactly one of `sd_conc` or `sd_const_multi` to be nonzero |
| `src_dry_sizes` | `dry_sizes_t` | - | Alternative source specification using size-number pairs |
| `src_sd_conc` | `unsigned long long` | `0` | Number of SDs created per cell per source iteration |
| `src_x0`, `src_x1` | `real_t` | `0` | Source box x-boundaries [m] (rounded to cell boundaries) |
| `src_y0`, `src_y1` | `real_t` | `0` | Source box y-boundaries [m] |
| `src_z0`, `src_z1` | `real_t` | `0` | Source box z-boundaries [m] |
| `supstp_src` | `int` | `1` | Timestep interval for applying source |

#### Aerosol Relaxation Options

Expand Down Expand Up @@ -385,7 +382,6 @@ opts_init.ny = 100;
opts_init.nz = 100;
opts_init.dx = opts_init.dy = opts_init.dz = 10; // 10 m grid spacing
opts_init.dt = 1.0; // 1 s timestep
opts_init.sd_conc = 64; // 64 SDs per cell
opts_init.sstp_cond = 4; // 4 condensation substeps
opts_init.adaptive_sstp_cond = true; // Enable adaptive substepping
opts_init.exact_sstp_cond = true; // Per-particle substepping
Expand All @@ -400,7 +396,7 @@ auto lognormal = [](double lnr) {
return n_tot * exp(-pow((lnr - log(mean_r)), 2) / 2 / pow(log(stdev), 2))
/ log(stdev) / sqrt(2 * M_PI);
};
opts_init.dry_distros[0.61] = std::make_shared<decltype(lognormal)>(lognormal);
opts_init.dry_distros[{0.61, 1., 64, 0}] = std::make_shared<decltype(lognormal)>(lognormal); // kappa, soluble_fraction, sd_conc, sd_const_multi

// Runtime options (can change each step)
libcloudphxx::lgrngn::opts_t<double> opts;
Expand Down
28 changes: 6 additions & 22 deletions include/libcloudph++/lgrngn/distro_t.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,35 +6,19 @@ namespace libcloudphxx
{
using common::unary_function;

template <typename real_t>
struct kappa_soluble_fraction_t {
real_t kappa;
real_t soluble_fraction; // volume fraction of the soluble part (0-1)

kappa_soluble_fraction_t(real_t kappa_, real_t soluble_fraction_)
: kappa(kappa_), soluble_fraction(soluble_fraction_)
{}

bool operator<(const kappa_soluble_fraction_t &other) const
{
if (kappa != other.kappa) return kappa < other.kappa;
return soluble_fraction < other.soluble_fraction;
}
};

// initial dry sizes of aerosol
// defined with a distribution
// uses shared_ptr to make opts_init copyable
template<typename real_t>
using dry_distros_t = std::map<
kappa_soluble_fraction_t<real_t>, // (kappa, soluble_fraction); dry_distros defines total dry radius
std::tuple<real_t, real_t, unsigned long long, unsigned long long>, // kappa, soluble_fraction, sd_conc, sd_const_multi
std::shared_ptr<unary_function<real_t>> // n(ln(rd)) @ STP; alternatively it's n(ln(rd)) independent of rhod if aerosol_independent_of_rhod=true
>;

// defined with a size-number pair
template<typename real_t>
using dry_sizes_t = std::map<
kappa_soluble_fraction_t<real_t>, // (kappa, soluble_fraction); dry_sizes defines total dry radius
std::tuple<real_t, real_t>, // (kappa, soluble_fraction); dry_sizes defines total dry radius
std::map<real_t, // radius [m]
std::pair<real_t, int> // STP_concentration [1/m^3], number of SD that represent this radius kappa and concentration
>
Expand All @@ -43,16 +27,16 @@ namespace libcloudphxx
// similar, but for sources of aerosols after initialization
template<typename real_t>
using src_dry_distros_t = std::map<
kappa_soluble_fraction_t<real_t>, // (kappa, soluble_fraction); src_dry_distros defines total dry radius
std::tuple<std::shared_ptr<unary_function<real_t>>, int, int> // 1st: n(ln(rd)) @ STP created per second; alternatively it's n(ln(rd)) independent of rhod if aerosol_independent_of_rhod=true; 2nd: sd_conc for this distribution ; 3rd: supstp for this aerosol (interval in timesteps beween addition of these aerosols)
std::tuple<real_t, real_t, int, int, int>, // kappa, soluble_fraction, sd_conc, sd_const_multi, supstp
std::shared_ptr<unary_function<real_t>> // n(ln(rd)) @ STP created per second; alternatively it's n(ln(rd)) independent of rhod if aerosol_independent_of_rhod=true
>;

// defined with a size-number pair
template<typename real_t>
using src_dry_sizes_t = std::map<
kappa_soluble_fraction_t<real_t>, // (kappa, soluble_fraction); src_dry_sizes defines total dry radius
std::tuple<real_t, real_t, int>, // kappa, soluble_fraction, supstp
std::map<real_t, // radius [m]
std::tuple<real_t, int, int> // STP_concentration [1/m^3] created per second, number of SD that represent this radius kappa and concentration, supstp
std::pair<real_t, int> // STP_concentration [1/m^3] created per second, number of SD that represent this radius kappa and concentration
>
>;

Expand Down
8 changes: 0 additions & 8 deletions include/libcloudph++/lgrngn/opts_init.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,9 +52,6 @@ namespace libcloudphxx
// Lagrangian domain extents
real_t x0, y0, z0, x1, y1, z1;

// no. of super-droplets per cell
unsigned long long sd_conc;

// should more SDs be added to better represent large tail of the distribution
bool sd_conc_large_tail;

Expand All @@ -65,9 +62,6 @@ namespace libcloudphxx
// is it allowed to change dt during simulation through opts.dt
bool variable_dt_switch;

// or, alternatively to sd_conc_mean, multiplicity of all SDs = const
unsigned long long sd_const_multi;

// max no. of super-droplets in the system
// should be enough to store particles from sources
unsigned long long n_sd_max;
Expand Down Expand Up @@ -187,10 +181,8 @@ namespace libcloudphxx
dx(1), dy(1), dz(1),
x0(0), y0(0), z0(0),
x1(1), y1(1), z1(1),
sd_conc(0),
sd_conc_large_tail(false),
aerosol_independent_of_rhod(false),
sd_const_multi(0),
dt(0),
sstp_cond(1), sstp_coal(1), sstp_chem(1), sstp_cond_act(1),
chem_switch(false), // chemical reactions turned off by default
Expand Down
31 changes: 11 additions & 20 deletions src/impl/initialization/particles_impl_init_SD_with_distros.ipp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,52 +14,43 @@ namespace libcloudphxx
template <typename real_t, backend_t device>
void particles_t<real_t, device>::impl::init_SD_with_distros()
{
// calc sum of ln(rd) ranges of all distributions
real_t tot_lnrd_rng = 0.;
if(opts_init.sd_conc > 0)
for (auto ddi = opts_init.dry_distros.cbegin(); ddi != opts_init.dry_distros.cend(); ++ddi)
{
init_dist_analysis_sd_conc(
*(ddi->second),
opts_init.sd_conc
);
tot_lnrd_rng += log_rd_max - log_rd_min;
}

// initialize SDs of each kappa-type
for (auto ddi = opts_init.dry_distros.cbegin(); ddi != opts_init.dry_distros.cend(); ++ddi)
{
if(opts_init.sd_conc > 0)
const auto &distro = ddi->second;
const auto sd_conc = std::get<2>(ddi->first);
const auto sd_const_multi = std::get<3>(ddi->first);
if(sd_conc > 0)
{
init_SD_with_distros_sd_conc(*(ddi->second), tot_lnrd_rng);
init_SD_with_distros_sd_conc(*distro, sd_conc);
init_SD_with_distros_finalize(ddi->first);

if(opts_init.sd_conc_large_tail)
{
init_SD_with_distros_tail(*(ddi->second), log_rd_max);
init_SD_with_distros_tail(*distro, log_rd_max);
init_SD_with_distros_finalize(ddi->first);
}
}
if(opts_init.sd_const_multi > 0)
if(sd_const_multi > 0)
{
init_SD_with_distros_const_multi(*(ddi->second));
init_SD_with_distros_const_multi(*distro, sd_const_multi);
init_SD_with_distros_finalize(ddi->first);
}
}
}

// final inits common for tail/sd_conc/const_multi
template <typename real_t, backend_t device>
void particles_t<real_t, device>::impl::init_SD_with_distros_finalize(const kappa_soluble_fraction_t<real_t> &kpa_sol_frac, const bool unravel_ijk_switch)
void particles_t<real_t, device>::impl::init_SD_with_distros_finalize(const std::tuple<real_t, real_t, unsigned long long, unsigned long long> &kpa_sol_frac, const bool unravel_ijk_switch)
{
// dry_distros defines total dry radius; insoluble part determined by soluble_fraction

// init kappa
init_kappa(kpa_sol_frac.kappa, kpa_sol_frac.soluble_fraction);
init_kappa(std::get<0>(kpa_sol_frac), std::get<1>(kpa_sol_frac));

if (opts_init.ice_switch)
{
init_insol(kpa_sol_frac.soluble_fraction);
init_insol(std::get<1>(kpa_sol_frac));

init_a_c_rho_ice();
if (! opts_init.time_dep_ice_nucl)
Expand Down
Loading
Loading