diff --git a/src/detail/negtozero.hpp b/src/detail/negtozero.hpp new file mode 100644 index 000000000..f159529b6 --- /dev/null +++ b/src/detail/negtozero.hpp @@ -0,0 +1,34 @@ +#pragma once + +// functor for testing not-positive values +struct isneg { + template + BOOST_GPU_ENABLED + bool operator()(const real_t a) const { + return a <= real_t(0); + } +}; + +#define _LIBCLOUDPHXX_SMALL_POSITIVE_VALUE 1e-10 + +// make not-positive values slightly positive... +#ifdef NDEBUG +#define negtozero(arr, name) {\ + thrust::replace_if(\ + arr.begin(), arr.end(), \ + isneg(), \ + real_t(_LIBCLOUDPHXX_SMALL_POSITIVE_VALUE) \ + );} +#else +#define negtozero(arr, name) {\ + auto min_position = thrust::min_element(arr.begin(), arr.end());\ + if(*min_position <= 0.)\ + {\ + printf("A non-positive number (%g) detected in: " name " at position %g\n", *min_position, min_position);\ + printf("Cheating in libcloudphxx: replacing non-positive values with small positive ones\n");\ + thrust::replace_if(\ + arr.begin(), arr.end(), \ + isneg(), \ + real_t(_LIBCLOUDPHXX_SMALL_POSITIVE_VALUE) \ + );}} +#endif diff --git a/src/impl/particles_impl_sstp.ipp b/src/impl/particles_impl_sstp.ipp index e69d66044..2fe961eed 100644 --- a/src/impl/particles_impl_sstp.ipp +++ b/src/impl/particles_impl_sstp.ipp @@ -88,6 +88,7 @@ namespace libcloudphxx arg::_1 + arg::_2 / sstp // op: rv = rv + drv_adv / sstp ); } + negtozero((*scl[ix]), "Eulerian field after per-cell sstp_step"); } } @@ -138,6 +139,7 @@ namespace libcloudphxx arg::_1 + arg::_2 // op: rv = rv + drv_adv ); } + negtozero((*tmp[ix]), "Eulerian field after per-particle sstp_step"); } } }; diff --git a/src/impl/particles_impl_update_th_rv.ipp b/src/impl/particles_impl_update_th_rv.ipp index f27fba0d3..5d06813aa 100644 --- a/src/impl/particles_impl_update_th_rv.ipp +++ b/src/impl/particles_impl_update_th_rv.ipp @@ -60,7 +60,7 @@ namespace libcloudphxx rv.begin(), // output thrust::plus() ); - assert(*thrust::min_element(rv.begin(), rv.end()) >= 0); + negtozero(rv, "rv after post-cond updating rv"); nancheck(rv, "update_th_rv: rv after update"); // updating th diff --git a/src/particles.tpp b/src/particles.tpp index 4ddd40cb8..7b15f7f92 100644 --- a/src/particles.tpp +++ b/src/particles.tpp @@ -18,6 +18,7 @@ #include "detail/wang_collision_enhancement.hpp" #include "detail/kernel_onishi_nograv.hpp" #include "detail/checknan.hpp" +#include "detail/negtozero.hpp" #include "detail/formatter.cpp" #include "detail/tpl_calc_wrapper.hpp" #include "detail/kernels.hpp"