Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix Roundoff Domain #3247

Merged
merged 1 commit into from
Apr 10, 2023
Merged

Fix Roundoff Domain #3247

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
146 changes: 121 additions & 25 deletions Src/Base/AMReX_Geometry.cpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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