Add guards to prevent divide by zero #481
Conversation
mjschmidt271
left a comment
mjschmidt271
left a comment
There was a problem hiding this comment.
I am strongly against merging this in its current form and would argue that major questions need to be addressed before including anything resembling these changes.
Also, note that gpu tests are failing as well as clang-format check.
| (bad_boltzmann * temperature)) * | ||
| Kcoll_dust_a3 * icnlx; | ||
|
|
||
| // ---- dt guards and small-dt finite-difference limits ---- |
There was a problem hiding this comment.
Why are we doing this? This raises some concerns for me:
- Assuming
dt < 1e-12is the chosen setting for a run, does it seem appropriate to effectively reducedtto zero for these calculations?- This isn't yet to the point of losing precision, and how small would we expect
dtto ever be?
- This isn't yet to the point of losing precision, and how small would we expect
- This change introduces a spurious infinity in order to give the same result as a comparison we've already done.
- This is an inappropriate way of accomplishing something much simpler, so can we outline what the goal is?
- Also, this leads to guaranteed unnecessary calculations:
- introducing 3
ifs to assign the desired values ofrate_over_dt_cnt<i>when the result is already determined by the initialif (dt < 1e-12). - introducing 3 calls to
min()when this is also predetermined - Note that this all assumes that the
Jvalues are strictly positive, but that seems to be the case because things would look straaaange otherwise
- introducing 3
- Now, assuming setting
dt = 0is an appropriate choice, nearly all of the new code is unnecessary.- It effectively defines division by zero to be equal to infinity (it is not) and then appears to take great pains to not explicitly divide by zero.
Here, for the quantities q<i> = {frzbccnt, frzducnt} and in the case that their corresponding do_<xyz> flag is true, this appears to lead to the same result as
if (dt <= 1e-12) {
// this will never be smaller than x2, so why do the comparison?
x1 = z * infinity;
x2 = z * J<i>;
q<i> += min(x1, x2);
} else {
// both of these should be strictly positive, and if I'm not mistaken,
// for fixed J > 0, we always have x2 >= x1
// so, again, introducing all of the comparisons is unnecessary
x1 = z * ((1.0 - haero::exp(-Jcnt_bc * deltat)) / deltat);
x2 = z / deltat;
q<i> += min(x1, x2);
}| // dust_a1 | ||
| Real coat_ratio1 = vol_shell[0] * (r_bc * 2.0) * fac_volsfc_pcarbon; | ||
|
|
||
| const Real n_so4_monolayers_dust = 1.0; |
There was a problem hiding this comment.
I strongly disagree with this approach. Where do the choices for these in-line constants come from?
Our goal has been to identify and subsequently get rid of all in-line constants in MAM4xx (work in progress--maybe one day). So, adding more is not something I can agree with.
| Real rhw_ss_clamped = haero::max(rhw_ss, 1.0 + 1.0e-12); | ||
| Real denom = haero::max(bad_boltzmann * temperature * haero::log(rhw_ss_clamped),1.0e-30); | ||
| const Real rgdep = 2.0 * vwice * sigma_iv / denom; | ||
| //(bad_boltzmann * temperature * haero::log(rhw_ss_clamped)); |
There was a problem hiding this comment.
please remove leftover comments
| } | ||
|
|
||
| KOKKOS_INLINE_FUNCTION | ||
| Real get_Aimm(const Real vwice, const Real rgimm, const Real temperature, |
There was a problem hiding this comment.
This is causing test failure for gpu. No idea why cpu still passes
Also, please delete leftover comments
| // ---- dt guards and small-dt finite-difference limits ---- | ||
| constexpr Real eps_dt = 1e-12; | ||
| const bool small_dt = haero::abs(deltat) <= eps_dt; | ||
| const Real inv_dt = small_dt ? std::numeric_limits<Real>::infinity() |
There was a problem hiding this comment.
@jeff-cohere @bartgol, can we use std::numeric_limits<Real>::infinity() in GPUs (CUDA)? What I recall is that we cannot use standard features inside GPU kernels, but I do not know if this is okay nowadays.
There was a problem hiding this comment.
It looks like Kokkos supports these now: https://kokkos.org/kokkos-core-wiki/API/core/numerics/numeric-traits.html
There was a problem hiding this comment.
Yeah, you can use the kokkos numeric traits.
| haero::min(Hetfrz::limfacbc * uncoated_aer_num[Hetfrz::id_bc] / deltat, | ||
| uncoated_aer_num[Hetfrz::id_bc] / deltat * | ||
| (1.0 - haero::exp(-Jcnt_bc * deltat))); | ||
| std::fmin(Hetfrz::limfacbc * uncoated_aer_num[Hetfrz::id_bc] * inv_dt, |
There was a problem hiding this comment.
Should we replace std::fmin with the Kokkos version of min?
| frzducnt + haero::min(uncoated_aer_num[Hetfrz::id_dst1] / deltat, | ||
| uncoated_aer_num[Hetfrz::id_dst1] / deltat * | ||
| (1.0 - haero::exp(-Jcnt_dust_a1 * deltat))); | ||
| frzducnt + std::fmin(uncoated_aer_num[Hetfrz::id_dst1] * inv_dt, |
| frzducnt + haero::min(uncoated_aer_num[Hetfrz::id_dst3] / deltat, | ||
| uncoated_aer_num[Hetfrz::id_dst3] / deltat * | ||
| (1.0 - haero::exp(-Jcnt_dust_a3 * deltat))); | ||
| frzducnt + std::fmin(uncoated_aer_num[Hetfrz::id_dst3] * inv_dt, |
please summarize your changes