diff --git a/Src/Base/AMReX_Geometry.cpp b/Src/Base/AMReX_Geometry.cpp index 5be69baa478..a240d62788b 100644 --- a/Src/Base/AMReX_Geometry.cpp +++ b/Src/Base/AMReX_Geometry.cpp @@ -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); @@ -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(plo); roundoff_hi[idim] = static_cast(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::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::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::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::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::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::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::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::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; } diff --git a/Src/Particle/AMReX_ParticleContainerI.H b/Src/Particle/AMReX_ParticleContainerI.H index db70ab6c963..52215eb71c4 100644 --- a/Src/Particle/AMReX_ParticleContainerI.H +++ b/Src/Particle/AMReX_ParticleContainerI.H @@ -306,45 +306,18 @@ void ParticleContainer::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 rhi = Geom(0).ProbHiArrayInParticleReal(); - GpuArray 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); diff --git a/Src/Particle/AMReX_ParticleUtil.H b/Src/Particle/AMReX_ParticleUtil.H index 0c4a5c74050..ec7b92ca9fd 100644 --- a/Src/Particle/AMReX_ParticleUtil.H +++ b/Src/Particle/AMReX_ParticleUtil.H @@ -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(phi[idim] - plo[idim]); } // clamp to avoid precision issues; @@ -584,12 +584,12 @@ bool enforcePeriodic (P& p, p.pos(idim) += static_cast(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;