Skip to content

Commit

Permalink
Merge pull request #441 from pdziekan/user_defined_rd_min_max
Browse files Browse the repository at this point in the history 
option to set rd bin edges in opts_init
  • Loading branch information
@pdziekan
pdziekan authored Jan 26, 2023
2 parents 8e98d21 + 3de4959 commit 59c6407
Show file tree
Hide file tree
Showing 2 changed files with 76 additions and 55 deletions.
5 changes: 3 additions & 2 deletions include/libcloudph++/lgrngn/opts_init.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,7 +126,7 @@ namespace libcloudphxx
// SGS mixing length profile [m]
std::vector<real_t> SGS_mix_len;

real_t rd_min; // minimal dry radius of droplets (works only for init from spectrum)
real_t rd_min, rd_max; // min/max dry radius of droplets [m]

bool no_ccn_at_init; // if true, no ccn / SD are put at the start of the simulation

Expand Down Expand Up @@ -230,7 +230,8 @@ namespace libcloudphxx
rlx_timescale(1),
rlx_sd_per_bin(0),
supstp_rlx(1),
rd_min(0.),
rd_min(-1), // negative means that rd_min will be automatically detected
rd_max(-1),
diag_incloud_time(false),
no_ccn_at_init(false),
open_side_walls(false),
Expand Down
126 changes: 73 additions & 53 deletions src/impl/particles_impl_dist_analysis.ipp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,17 +20,10 @@ namespace libcloudphxx
const real_t dt
)
{
// probing the spectrum to find rd_min-rd_max range
// when analysing distro for source, multiplier takes into account that
// the distribution is assumed to represent number of particles created per unit of time!
// TODO: document that

// values to start the search
real_t rd_min = config.rd_min_init, rd_max = config.rd_max_init;

bool found_optimal_range = false;
while (!found_optimal_range)
if(opts_init.rd_min >= 0 && opts_init.rd_max >= 0) // user-defined bin edges
{
real_t rd_min = opts_init.rd_min, rd_max = opts_init.rd_max;

multiplier = log(rd_max / rd_min)
/ sd_conc
* dt
Expand All @@ -41,62 +34,89 @@ namespace libcloudphxx

log_rd_min = log(rd_min);
log_rd_max = log(rd_max);
}
else if (opts_init.rd_min < 0 && opts_init.rd_max < 0) // automatic detection of bin edges
{
// probing the spectrum to find rd_min-rd_max range
// when analysing distro for source, multiplier takes into account that
// the distribution is assumed to represent number of particles created per unit of time!
// TODO: document that

// values to start the search
real_t rd_min = config.rd_min_init, rd_max = config.rd_max_init;

bool found_optimal_range = false;
while (!found_optimal_range)
{
multiplier = log(rd_max / rd_min)
/ sd_conc
* dt
* (n_dims == 0
? dv[0]
: (opts_init.dx * opts_init.dy * opts_init.dz)
);

impl::n_t
n_min = n_of_lnrd_stp(log_rd_min) * multiplier,
n_max = n_of_lnrd_stp(log_rd_max) * multiplier;
log_rd_min = log(rd_min);
log_rd_max = log(rd_max);

if (rd_min == config.rd_min_init && n_min != 0)
throw std::runtime_error(detail::formatter() << "Initial dry radii distribution is non-zero (" << n_min << ") for rd_min_init (" << config.rd_min_init <<")");
if (rd_max == config.rd_max_init && n_max != 0)
throw std::runtime_error(detail::formatter() << "Initial dry radii distribution is non-zero (" << n_max << ") for rd_max_init (" << config.rd_max_init <<")");
impl::n_t
n_min = n_of_lnrd_stp(log_rd_min) * multiplier,
n_max = n_of_lnrd_stp(log_rd_max) * multiplier;

if (n_min == 0) rd_min *= 1.01;
else if (n_max == 0) rd_max /= 1.01;
else found_optimal_range = true;
if (rd_min == config.rd_min_init && n_min != 0)
throw std::runtime_error(detail::formatter() << "Initial dry radii distribution is non-zero (" << n_min << ") for rd_min_init (" << config.rd_min_init <<")");
if (rd_max == config.rd_max_init && n_max != 0)
throw std::runtime_error(detail::formatter() << "Initial dry radii distribution is non-zero (" << n_max << ") for rd_max_init (" << config.rd_max_init <<")");

if (n_min == 0) rd_min *= 1.01;
else if (n_max == 0) rd_max /= 1.01;
else found_optimal_range = true;
}
}
#if !defined(__NVCC__)
using std::max;
#endif
log_rd_min = max(log_rd_min, real_t(log(opts_init.rd_min))); // user-defined lower limit for rd
else assert(false && "opts_init.rd_min * opts_init.rd_max < 0");
};

template <typename real_t, backend_t device>
void particles_t<real_t, device>::impl::dist_analysis_const_multi(
const common::unary_function<real_t> &n_of_lnrd_stp
)
{
// TODO: add same sanity check as above
// how does brent algorithm work for functions with multiple minima??
std::pair<real_t, real_t> init_distr_max; // [ln(position of distribution's maximum), -function value at maximum]
boost::uintmax_t n_iter = config.n_iter;
init_distr_max = boost::math::tools::brent_find_minima(detail::eval_and_mul<real_t>(n_of_lnrd_stp, -1), log(config.rd_min_init), log(config.rd_max_init), 200, n_iter); // bits = 200 to make algorithm choose max precision available
if(opts_init.rd_min >= 0 && opts_init.rd_max >= 0) // user-defined bin edges
{
log_rd_min = log(opts_init.rd_min);
log_rd_max = log(opts_init.rd_max);
}
else if (opts_init.rd_min < 0 && opts_init.rd_max < 0) // automatic detection of bin edges
{
// TODO: add same sanity check as above
// how does brent algorithm work for functions with multiple minima??
std::pair<real_t, real_t> init_distr_max; // [ln(position of distribution's maximum), -function value at maximum]
boost::uintmax_t n_iter = config.n_iter;
init_distr_max = boost::math::tools::brent_find_minima(detail::eval_and_mul<real_t>(n_of_lnrd_stp, -1), log(config.rd_min_init), log(config.rd_max_init), 200, n_iter); // bits = 200 to make algorithm choose max precision available

real_t init_dist_bound_value = -init_distr_max.second / config.threshold; // value of the distribution at which we bind it
n_iter = config.n_iter;
// TODO: it could be written more clearly by creating an object detail::eval_and_oper<real_t>(*n_of_lnrd_stp, -init_dist_bound_value, 1), but for some reason it doesnt give the correct values
log_rd_min =
common::detail::toms748_solve(
detail::eval_and_add<real_t>(n_of_lnrd_stp, -init_dist_bound_value),
real_t(log(config.rd_min_init)), init_distr_max.first,
detail::eval_and_add<real_t>(n_of_lnrd_stp, -init_dist_bound_value)(real_t(log(config.rd_min_init))),
detail::eval_and_add<real_t>(n_of_lnrd_stp, -init_dist_bound_value)(init_distr_max.first),
config.eps_tolerance, n_iter
);
real_t init_dist_bound_value = -init_distr_max.second / config.threshold; // value of the distribution at which we bind it
n_iter = config.n_iter;
// TODO: it could be written more clearly by creating an object detail::eval_and_oper<real_t>(*n_of_lnrd_stp, -init_dist_bound_value, 1), but for some reason it doesnt give the correct values
log_rd_min =
common::detail::toms748_solve(
detail::eval_and_add<real_t>(n_of_lnrd_stp, -init_dist_bound_value),
real_t(log(config.rd_min_init)), init_distr_max.first,
detail::eval_and_add<real_t>(n_of_lnrd_stp, -init_dist_bound_value)(real_t(log(config.rd_min_init))),
detail::eval_and_add<real_t>(n_of_lnrd_stp, -init_dist_bound_value)(init_distr_max.first),
config.eps_tolerance, n_iter
);

n_iter = config.n_iter;
log_rd_max =
common::detail::toms748_solve(
detail::eval_and_add<real_t>(n_of_lnrd_stp, -init_dist_bound_value),
init_distr_max.first, real_t(log(config.rd_max_init)),
detail::eval_and_add<real_t>(n_of_lnrd_stp, -init_dist_bound_value)(init_distr_max.first),
detail::eval_and_add<real_t>(n_of_lnrd_stp, -init_dist_bound_value)(real_t(log(config.rd_max_init))),
config.eps_tolerance, n_iter
);
#if !defined(__NVCC__)
using std::max;
#endif
log_rd_min = max(log_rd_min, real_t(log(opts_init.rd_min))); // user-defined lower limit for rd
n_iter = config.n_iter;
log_rd_max =
common::detail::toms748_solve(
detail::eval_and_add<real_t>(n_of_lnrd_stp, -init_dist_bound_value),
init_distr_max.first, real_t(log(config.rd_max_init)),
detail::eval_and_add<real_t>(n_of_lnrd_stp, -init_dist_bound_value)(init_distr_max.first),
detail::eval_and_add<real_t>(n_of_lnrd_stp, -init_dist_bound_value)(real_t(log(config.rd_max_init))),
config.eps_tolerance, n_iter
);
}
else assert(false && "opts_init.rd_min * opts_init.rd_max < 0");
}
};
};

0 comments on commit 59c6407

Please sign in to comment.