Messages for Asserts in `computeRoundoffDomain`

[ax3l]

I [...] find the assert in computeRoundoffDomain particularly annoying, as it has a very nondescript message, and I think it even sometimes triggers for perfectly valid geometries if one of the sides is at exactly zero. Annoyingly, the roundoff domain is not even used in our code.

Originally posted by @AlexanderSinn in #5219 (comment)

Currently, those throw something that is not actionable for users. Looking into the code it already looks like a work-around/hack what is in there. We need to revisit those, add messages, etc.

[ax3l]

amrex/Src/Base/AMReX_Geometry.cpp

Lines 535 to 682 in 2fd316b

	// Check that the grid is well formed and that deltax > roundoff
	AMREX_ASSERT((plo + ihi*CellSize(idim)) < (plo + (ihi + 1)*CellSize(idim)));
	// roundoff_lo will be the lowest value that will be inside the domain
	// roundoff_hi will be the highest value that will be inside the domain
	roundoff_lo[idim] = static_cast<ParticleReal>(plo);
	roundoff_hi[idim] = static_cast<ParticleReal>(phi);
	auto& rlo = roundoff_lo[idim];
	auto& rhi = roundoff_hi[idim];
	auto is_outside = [=] (auto x) -> bool
	{
	auto idx = int(std::floor((x - plo)*dxinv));
	return (idx < 0) || (idx >= ncells);
	};
	auto is_inside = [=] (auto x) -> bool
	{
	return !is_outside(x);
	};
	ParticleReal rlo_out;
	if (is_inside(rlo))
	{
	auto eps = std::numeric_limits<ParticleReal>::epsilon() * (rhi-rlo);
	int iters = 0;
	rlo_out = rlo - eps;
	while (is_inside(rlo_out) && iters < maxiters) {
	eps *= ParticleReal(2.);
	rlo_out = rlo - eps;
	++iters;
	}
	// The assertion on rlo_out makes sure the compiler cannot optimize it away.
	AMREX_ALWAYS_ASSERT(rlo_out > std::numeric_limits<ParticleReal>::lowest()
	&& iters < maxiters);
	}
	else
	{
	auto eps = std::numeric_limits<ParticleReal>::epsilon() * (rhi-rlo);
	int iters = 0;
	auto rtmp = rlo + eps;
	while (is_outside(rtmp) && iters < maxiters) {
	eps *= ParticleReal(2.);
	rtmp = rlo + eps;
	++iters;
	}
	rlo_out = rlo;
	rlo = rtmp;
	// The assertion on rtmp makes sure the compiler cannot optimize it away.
	AMREX_ALWAYS_ASSERT(rtmp > std::numeric_limits<ParticleReal>::lowest()
	&& iters < maxiters);
	}
	{
	// rlo_out is outside
	// rlo is inside
	int iters = 0;
	auto epsilon = std::numeric_limits<ParticleReal>::epsilon()
	* std::max(ParticleReal(CellSize(idim)),std::abs(rlo))
	* ParticleReal(2.0);
	auto rlo_minus = rlo-epsilon;
	bool rlo_minus_is_inside = is_inside(rlo_minus);
	while (rlo_minus_is_inside && iters < maxiters) {
	auto rmid = ParticleReal(0.5)*(rlo_out+rlo);
	if (rmid == rlo_out || rmid == rlo) {
	break;
	}
	if (is_inside(rmid)) {
	rlo = rmid;
	epsilon = std::numeric_limits<ParticleReal>::epsilon()
	* std::max(ParticleReal(CellSize(idim)),std::abs(rlo))
	* ParticleReal(2.0);
	rlo_minus = rlo - epsilon;
	rlo_minus_is_inside = is_inside(rlo_minus);
	} else {
	rlo_out = rmid;
	}
	++iters;
	}
	// The assertion on rlo_minus makes sure the compiler cannot optimize it away.
	AMREX_ALWAYS_ASSERT(rlo_minus > std::numeric_limits<ParticleReal>::lowest()
	&& iters < maxiters);
	}
	ParticleReal rhi_out;
	if (is_inside(rhi))
	{
	auto eps = std::numeric_limits<ParticleReal>::epsilon() * (rhi-rlo);
	int iters = 0;
	rhi_out = rhi + eps;
	while (is_inside(rhi_out) && iters < maxiters) {
	eps *= ParticleReal(2.);
	rhi_out = rhi + eps;
	++iters;
	}
	// The assertion on rhi_out makes sure the compiler cannot optimize it away.
	AMREX_ALWAYS_ASSERT(rhi_out > std::numeric_limits<ParticleReal>::lowest()
	&& iters < maxiters);
	}
	else
	{
	auto eps = std::numeric_limits<ParticleReal>::epsilon() * (rhi-rlo);
	int iters = 0;
	// Yes, we have to write it this way for Intel compiler.
	// is_outside(rhi-eps) could be different from is_outside(rtmp),
	// where rtmp = rhs-eps.
	auto rtmp = rhi - eps;
	while (is_outside(rtmp) && iters < maxiters) {
	eps *= ParticleReal(2.);
	rtmp = rhi - eps;
	++iters;
	}
	rhi_out = rhi;
	rhi = rtmp;
	// The assertion on rtmp makes sure the compiler cannot optimize it away.
	AMREX_ALWAYS_ASSERT(rtmp > std::numeric_limits<ParticleReal>::lowest()
	&& iters < maxiters);
	}
	{
	// rhi is inside
	// rhi_out is outside
	int iters = 0;
	auto epsilon = std::numeric_limits<ParticleReal>::epsilon()
	* std::max(ParticleReal(CellSize(idim)),std::abs(rhi))
	* ParticleReal(2.0);
	auto rhi_plus = rhi+epsilon;
	bool rhi_plus_is_inside = is_inside(rhi_plus);
	while (rhi_plus_is_inside && iters < maxiters) {
	auto rmid = ParticleReal(0.5)*(rhi+rhi_out);
	if (rmid == rhi || rmid == rhi_out) {
	break;
	}
	if (is_inside(rmid)) {
	rhi = rmid;
	epsilon = std::numeric_limits<ParticleReal>::epsilon()
	* std::max(ParticleReal(CellSize(idim)),std::abs(rhi))
	* ParticleReal(2.0);
	rhi_plus = rhi + epsilon;
	rhi_plus_is_inside = is_inside(rhi_plus);
	} else {
	rhi_out = rmid;
	}
	++iters;
	}
	// The assertion on rhi_plus makes sure the compiler cannot optimize it away.
	AMREX_ALWAYS_ASSERT(rhi_plus > std::numeric_limits<ParticleReal>::lowest()
	&& iters < maxiters);

[AlexanderSinn]

Also see Hi-PACE/hipace#1364 and BLAST-ImpactX/impactx#1348

[WeiqunZhang]

Let's try to figure out why the hipace++ case failed. Could you add a line to the beginning of computeRoundoffDomain?

amrex::Print().SetPrecision(17) << "computeRoundoffDomain: " << *this << std::endl;

What compiler was used?

[AlexanderSinn]

In that specific issue, the domain box was not ok() with (lo[2] == hi[2] + 1), but I also remember it failing with a valid box in another case. I would have to experiment myself again to find a case where it fails.

[WeiqunZhang]

For the invalid input case, we could add an assert at the beginning of the computeRoundoffDomain with an error message.

[ax3l]

Improved via #5224