Skip to content

Commit

Permalink
Rework handling of roundoff domain extent (AMReX-Codes#3199)
Browse files Browse the repository at this point in the history 
## Summary
This reworks previous changes that created the round off domain which
was done to avoid round off errors when locating particles in the
domain, specifically when handling periodic boundaries. This PR changes
how this is done so that the roundoff domain makes as small a change as
possible.

## Additional background
Previous PRs AMReX-Codes#2839 and AMReX-Codes#2950 created a round off domain that has a
slightly smaller extent that the actual domain. Any particles that ended
up outside the round off domain but inside the actual domain were
shifted to be on the boundary of the round off domain. The main purpose
was to handle problems with periodic boundary with mixed precision
(single precision particles and double precision grid). In that case, it
is possible for a particle to have `z < lo` and `z+(hi-lo) > hi`. It is
also possible for a particle to have `z < hi` but `(z-lo)*dzinv > ihi`
putting it in a grid cell beyond the domain. Unfortunately, the roundoff
domains was being set further inside that was needed and would affect
valid particles that would not have those issues. While the problem is
small, the effects can be noticeable when examining simulation
properties, for example charge conservation in WarpX.

This PR changes the roundoff domain so that only problematic particles
are affected. The code now uses `std::nextafter` to find the exact place
where the problems occurs. At the lower end, it ensures that the casted
position will be greater than ProbLo. At the upper end, it ensures that
`(z-lo)*dzinv <= ihi`.

One comment is that this assumes that particle grid coordinates will
always be calculated at the higher precision, with expressions like `iz
= (int)(z - lo)`, without `lo` being cast to `ParticleReal`.

Co-authored-by: atmyers <[email protected]>
  • Loading branch information
@dpgrote @atmyers @guj
2 people authored and guj committed Jul 13, 2023
1 parent b67e1b6 commit 7ef97e2
Show file tree
Hide file tree
Showing 2 changed files with 41 additions and 87 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
57 changes: 35 additions & 22 deletions Src/Base/AMReX_Geometry.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -525,36 +525,49 @@ Geometry::computeRoundoffDomain ()
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 lowest value that will be outside the domain
roundoff_lo[idim] = static_cast<ParticleReal>(plo);
roundoff_hi[idim] = static_cast<ParticleReal>(phi);

// 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++;
}
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;
}
iters++;
}

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);
// 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);
}
}

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

0 comments on commit 7ef97e2

Please sign in to comment.