Skip to content

Commit

Permalink
Commits for roundoff domain
Browse files Browse the repository at this point in the history 
  * Rework handling of roundoff domain extent (AMReX-Codes#3199)

  * Fix periodic boundary bug in AMReX-Codes#3199 (AMReX-Codes#3243)

  * Fix Roundoff Domain (AMReX-Codes#3247)

  * Add Tests/RoundoffDomain (AMReX-Codes#3248)

  * Roundoff errors in computeRoundoffDomain (AMReX-Codes#3255)
  • Loading branch information
@dpgrote @WeiqunZhang
dpgrote authored and WeiqunZhang committed Sep 28, 2023
1 parent 5631b6e commit 47b2e04
Show file tree
Hide file tree
Showing 9 changed files with 308 additions and 128 deletions.
71 changes: 6 additions & 65 deletions Src/Base/AMReX_Geometry.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,56 +67,6 @@ public:
int coord;
};

namespace detail {
template <typename T>
T bisect_prob_lo (amrex::Real plo, amrex::Real /*phi*/, amrex::Real dxinv, int ilo, int ihi, amrex::Real tol) {
T lo = static_cast<T>(plo + tol);
bool safe;
{
int i = int(std::floor((lo - plo)*dxinv)) + ilo;
safe = i >= ilo && i <= ihi;
}
if (safe) {
return lo;
} else {
// bisect the point at which the cell no longer maps to inside the domain
T hi = static_cast<T>(plo + 0.5_rt/dxinv);
T mid = bisect(lo, hi,
[=] AMREX_GPU_HOST_DEVICE (T x) -> T
{
int i = int(std::floor((x - plo)*dxinv)) + ilo;
bool inside = i >= ilo && i <= ihi;
return static_cast<T>(inside) - T(0.5);
}, static_cast<T>(tol));
return mid - static_cast<T>(tol);
}
}

template <typename T>
T bisect_prob_hi (amrex::Real plo, amrex::Real phi, amrex::Real dxinv, int ilo, int ihi, amrex::Real tol) {
T hi = static_cast<T>(phi - tol);
bool safe;
{
int i = int(std::floor((hi - plo)*dxinv)) + ilo;
safe = i >= ilo && i <= ihi;
}
if (safe) {
return hi;
} else {
// bisect the point at which the cell no longer maps to inside the domain
T lo = static_cast<T>(phi - 0.5_rt/dxinv);
T mid = bisect(lo, hi,
[=] AMREX_GPU_HOST_DEVICE (T x) -> T
{
int i = int(std::floor((x - plo)*dxinv)) + ilo;
bool inside = i >= ilo && i <= ihi;
return static_cast<T>(inside) - T(0.5);
}, static_cast<T>(tol));
return mid - static_cast<T>(tol);
}
}
}

class Geometry
:
public CoordSys
Expand Down Expand Up @@ -242,18 +192,10 @@ public:
}

[[nodiscard]] GpuArray<ParticleReal,AMREX_SPACEDIM> ProbLoArrayInParticleReal () const noexcept {
#ifdef AMREX_SINGLE_PRECISION_PARTICLES
return roundoff_lo_f;
#else
return roundoff_lo_d;
#endif
return roundoff_lo;
}
[[nodiscard]] GpuArray<ParticleReal,AMREX_SPACEDIM> ProbHiArrayInParticleReal () const noexcept {
#ifdef AMREX_SINGLE_PRECISION_PARTICLES
return roundoff_hi_f;
#else
return roundoff_hi_d;
#endif
return roundoff_hi;
}

//! Returns the overall size of the domain by multiplying the ProbLength's together
Expand Down Expand Up @@ -489,11 +431,10 @@ private:
RealBox prob_domain;

// Due to round-off errors, not all floating point numbers for which plo >= x < phi
// will map to a cell that is inside "domain". "roundoff_{lo,hi}_{f,d}" each store
// a position that is very close to that in prob_domain, and for which all doubles and floats less than
// it will map to a cell inside domain.
GpuArray<double, AMREX_SPACEDIM> roundoff_lo_d, roundoff_hi_d;
GpuArray<float , AMREX_SPACEDIM> roundoff_lo_f, roundoff_hi_f;
// will map to a cell that is inside "domain". "roundoff_{lo,hi}" store
// positions that are very close to that in prob_domain, and for which all doubles and floats
// between them will map to a cell inside domain.
GpuArray<ParticleReal , AMREX_SPACEDIM> roundoff_lo, roundoff_hi;

//
Box domain;
Expand Down
179 changes: 157 additions & 22 deletions Src/Base/AMReX_Geometry.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -518,43 +518,178 @@ Geometry::computeRoundoffDomain ()
inv_dx[k] = 1.0_rt/dx[k];
}

constexpr int maxiters = 200;

for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
{
int ilo = Domain().smallEnd(idim);
int ihi = Domain().bigEnd(idim);
int ncells = ihi-ilo+1;
Real plo = ProbLo(idim);
Real phi = ProbHi(idim);
Real dxinv = InvCellSize(idim);
Real deltax = CellSize(idim);

Real ftol = std::max(1.e-4_rt*deltax, 2.e-7_rt*phi);
Real dtol = std::max(1.e-8_rt*deltax, 1.e-14_rt*phi);
// 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];

roundoff_lo_f[idim] = detail::bisect_prob_lo<float> (plo, phi, dxinv, ilo, ihi, ftol);
roundoff_lo_d[idim] = detail::bisect_prob_lo<double>(plo, phi, dxinv, ilo, ihi, dtol);
roundoff_hi_f[idim] = detail::bisect_prob_hi<float> (plo, phi, dxinv, ilo, ihi, ftol);
roundoff_hi_d[idim] = detail::bisect_prob_hi<double>(plo, phi, dxinv, ilo, ihi, dtol);
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);
}
}
}

bool
Geometry::outsideRoundoffDomain (AMREX_D_DECL(ParticleReal x, ParticleReal y, ParticleReal z)) const
{
#ifdef AMREX_SINGLE_PRECISION_PARTICLES
bool outside = AMREX_D_TERM(x < roundoff_lo_f[0]
|| x >= roundoff_hi_f[0],
|| y < roundoff_lo_f[1]
|| y >= roundoff_hi_f[1],
|| z < roundoff_lo_f[2]
|| z >= roundoff_hi_f[2]);
#else
bool outside = AMREX_D_TERM(x < roundoff_lo_d[0]
|| x >= roundoff_hi_d[0],
|| y < roundoff_lo_d[1]
|| y >= roundoff_hi_d[1],
|| z < roundoff_lo_d[2]
|| z >= roundoff_hi_d[2]);
#endif
bool outside = AMREX_D_TERM(x < roundoff_lo[0]
|| x > roundoff_hi[0],
|| y < roundoff_lo[1]
|| y > roundoff_hi[1],
|| z < roundoff_lo[2]
|| z > roundoff_hi[2]);
return outside;
}

Expand Down
43 changes: 8 additions & 35 deletions Src/Particle/AMReX_ParticleContainerI.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -306,45 +306,18 @@ void
ParticleContainer<NStructReal, NStructInt, NArrayReal, NArrayInt, Allocator>::locateParticle (ParticleType& p, ParticleLocData& pld,
int lev_min, int lev_max, int nGrow, int local_grid) const
{
bool outside = AMREX_D_TERM(p.pos(0) < Geom(0).ProbLo(0)
|| p.pos(0) >= Geom(0).ProbHi(0),
|| p.pos(1) < Geom(0).ProbLo(1)
|| p.pos(1) >= Geom(0).ProbHi(1),
|| p.pos(2) < Geom(0).ProbLo(2)
|| p.pos(2) >= Geom(0).ProbHi(2));

if (! outside)
bool success;
if (Geom(0).outsideRoundoffDomain(AMREX_D_DECL(p.pos(0), p.pos(1), p.pos(2))))
{
if (Geom(0).outsideRoundoffDomain(AMREX_D_DECL(p.pos(0), p.pos(1), p.pos(2))))
// Note that EnforcePeriodicWhere may shift the particle if it is successful.
success = EnforcePeriodicWhere(p, pld, lev_min, lev_max, local_grid);
if (!success && lev_min == 0)
{
GpuArray<ParticleReal, AMREX_SPACEDIM> rhi = Geom(0).ProbHiArrayInParticleReal();
GpuArray<ParticleReal, AMREX_SPACEDIM> rlo = Geom(0).ProbLoArrayInParticleReal();
for (int idim=0; idim < AMREX_SPACEDIM; ++idim)
{
if (p.pos(idim) <= rlo[idim]) {
p.pos(idim) = std::nextafter(rlo[idim], rhi[idim]);
}
if (p.pos(idim) >= rhi[idim]) {
p.pos(idim) = std::nextafter(rhi[idim], rlo[idim]);
}
}

AMREX_ASSERT(! Geom(0).outsideRoundoffDomain(AMREX_D_DECL(p.pos(0), p.pos(1), p.pos(2))));
// The particle has left the domain; invalidate it.
p.id() = -p.id();
success = true;
}
}

bool success;
if (outside)
{
// Note that EnforcePeriodicWhere may shift the particle if it is successful.
success = EnforcePeriodicWhere(p, pld, lev_min, lev_max, local_grid);
if (!success && lev_min == 0)
{
// The particle has left the domain; invalidate it.
p.id() = -p.id();
success = true;
}
}
else
{
success = Where(p, pld, lev_min, lev_max, 0, local_grid);
Expand Down
10 changes: 5 additions & 5 deletions Src/Particle/AMReX_ParticleUtil.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -569,8 +569,8 @@ bool enforcePeriodic (P& p,
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
{
if (! is_per[idim]) continue;
if (p.pos(idim) >= phi[idim]) {
while (p.pos(idim) >= phi[idim]) {
if (p.pos(idim) > rhi[idim]) {
while (p.pos(idim) > rhi[idim]) {
p.pos(idim) -= static_cast<ParticleReal>(phi[idim] - plo[idim]);
}
// clamp to avoid precision issues;
Expand All @@ -579,8 +579,8 @@ bool enforcePeriodic (P& p,
}
shifted = true;
}
else if (p.pos(idim) < plo[idim]) {
while (p.pos(idim) < plo[idim]) {
else if (p.pos(idim) < rlo[idim]) {
while (p.pos(idim) < rlo[idim]) {
p.pos(idim) += static_cast<ParticleReal>(phi[idim] - plo[idim]);
}
// clamp to avoid precision issues;
Expand All @@ -589,7 +589,7 @@ bool enforcePeriodic (P& p,
}
shifted = true;
}
AMREX_ASSERT( (p.pos(idim) >= plo[idim] ) && ( p.pos(idim) < phi[idim] ));
AMREX_ASSERT( (p.pos(idim) >= rlo[idim] ) && ( p.pos(idim) <= rhi[idim] ));
}

return shifted;
Expand Down
2 changes: 1 addition & 1 deletion Tests/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
#
# List of subdirectories to search for CMakeLists.
#
set( AMREX_TESTS_SUBDIRS AsyncOut MultiBlock Reinit Amr CLZ Parser CTOParFor)
set( AMREX_TESTS_SUBDIRS AsyncOut MultiBlock Reinit Amr CLZ Parser CTOParFor RoundoffDomain)

if (AMReX_PARTICLES)
list(APPEND AMREX_TESTS_SUBDIRS Particles)
Expand Down
8 changes: 8 additions & 0 deletions Tests/RoundoffDomain/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,8 @@
set(_sources main.cpp)
set(_input_files)

setup_test(_sources _input_files)

unset(_sources)
unset(_input_files)

Loading

0 comments on commit 47b2e04

Please sign in to comment.