Skip to content

Commit

Permalink
Fix Roundoff Domain
Browse files Browse the repository at this point in the history 
There was a flaw in amrex::Geometry::computeRoundoffDomain. If probhi is
close to zero, std::nextafter towards negative infinity will be a very small
number (e.g., 1.e-300). So the function will be able to find the lowest
value of roundoff_hi that's outside the domain within a reasonable number of
iterations.

In the new implementation, we use bisection to find the lowest value of
roundoff_lo that's inside the domain and the highest value fo roundoff_hi
that's inside the domain, up to a tolerance.
  • Loading branch information
@WeiqunZhang
WeiqunZhang committed Apr 10, 2023
1 parent 8cae802 commit 5a1ca61
Show file tree
Hide file tree
Showing 3 changed files with 134 additions and 65 deletions.
146 changes: 121 additions & 25 deletions Src/Base/AMReX_Geometry.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -518,10 +518,14 @@ Geometry::computeRoundoffDomain ()
inv_dx[k] = 1.0_rt/dx[k];
}

constexpr int maxiters = 200;
int numiters = 0;

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);
Expand All @@ -530,44 +534,136 @@ Geometry::computeRoundoffDomain ()
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 lowest value that will be outside 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);
};

// Hopefully, the loops should never take more than 1 or 2 iterations.
// In obscure cases, more may be needed, for example if (phi - plo) is near round off
// compared to (phi + plo).
int iters = 0;
while (roundoff_lo[idim] < plo && iters < 20) {
roundoff_lo[idim] = std::nextafter(roundoff_lo[idim], roundoff_hi[idim]);
iters++;
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;
while (is_inside(rlo-eps) && iters < maxiters) {
eps *= ParticleReal(2.);
++iters;
}
rlo_out = rlo - eps;
numiters = std::max(numiters,iters);
}
ParticleReal rhi = roundoff_hi[idim];
while (iters < 20) {
rhi = std::nextafter(rhi, roundoff_lo[idim]);
if (int(std::floor((rhi - plo)*dxinv)) >= ihi + 1 - ilo) {
roundoff_hi[idim] = rhi;
} else {
break;
else
{
auto eps = std::numeric_limits<ParticleReal>::epsilon() * (rhi-rlo);
int iters = 0;
while (is_outside(rlo+eps) && iters < maxiters) {
eps *= ParticleReal(2.);
++iters;
}
rlo_out = rlo;
rlo += eps;
numiters = std::max(numiters,iters);
}

{
// 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);
while (is_inside(rlo-epsilon) && 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);
} else {
rlo_out = rmid;
}
++iters;
}
iters++;
numiters = std::max(numiters,iters);
}

// If iters is at the limit, then the grid is ill formed and reasonable values could not
// be found for the round off domain extent
AMREX_ASSERT(iters < 20);
ParticleReal rhi_out;
if (is_inside(rhi))
{
auto eps = std::numeric_limits<ParticleReal>::epsilon() * (rhi-rlo);
int iters = 0;
while (is_inside(rhi+eps) && iters < maxiters) {
eps *= ParticleReal(2.);
++iters;
}
rhi_out = rhi + eps;
numiters = std::max(numiters,iters);
}
else
{
auto eps = std::numeric_limits<ParticleReal>::epsilon() * (rhi-rlo);
int iters = 0;
while (is_outside(rhi-eps) && iters < maxiters) {
eps *= ParticleReal(2.);
++iters;
}
rhi_out = rhi;
rhi -= eps;
numiters = std::max(numiters,iters);
}

{
// 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);
while (is_inside(rhi+epsilon) && 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);
} else {
rhi_out = rmid;
}
++iters;
}
numiters = std::max(numiters,iters);
}
}

AMREX_ALWAYS_ASSERT(numiters < maxiters);
}

bool
Geometry::outsideRoundoffDomain (AMREX_D_DECL(ParticleReal x, ParticleReal y, ParticleReal z)) const
{
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]);
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) >= rhi[idim]) {
while (p.pos(idim) >= rhi[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 @@ -584,12 +584,12 @@ bool enforcePeriodic (P& p,
p.pos(idim) += static_cast<ParticleReal>(phi[idim] - plo[idim]);
}
// clamp to avoid precision issues;
if (p.pos(idim) >= rhi[idim]) {
p.pos(idim) = std::nextafter(p.pos(idim), rlo[idim]);
if (p.pos(idim) > rhi[idim]) {
p.pos(idim) = rhi[idim];
}
shifted = true;
}
AMREX_ASSERT( (p.pos(idim) >= plo[idim] ) && ( p.pos(idim) < phi[idim] ));
AMREX_ASSERT( (p.pos(idim) >= plo[idim] ) && ( p.pos(idim) <= phi[idim] ));
}

return shifted;
Expand Down

0 comments on commit 5a1ca61

Please sign in to comment.