Skip to content

Commit

Permalink
Add roundoff_lo corresponding to roundoff_hi for domains that don't s…
Browse files Browse the repository at this point in the history 
…tart at 0 (AMReX-Codes#2950)

* Lay groundwork for roundoff_lo

* Add dummy implementation of roundoff_lo computation

* implement bisect_prob_lo

* change idx -> dxinv

* use rlo instead of plo in locateParticle

Co-authored-by: atmyers <atmyers2@gmail.com>
  • Loading branch information
@PhilMiller @atmyers
PhilMiller and atmyers authored Sep 16, 2022
1 parent 6a5a056 commit 826cd37
Show file tree
Hide file tree
Showing 4 changed files with 67 additions and 28 deletions.
47 changes: 39 additions & 8 deletions Src/Base/AMReX_Geometry.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,22 +69,46 @@ public:

namespace detail {
template <typename T>
T bisect_prob_hi (amrex::Real plo, amrex::Real phi, amrex::Real idx, int ilo, int ihi, amrex::Real tol) {
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(Math::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(Math::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(Math::floor((hi - plo)*idx)) + ilo;
int i = int(Math::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/idx);
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(Math::floor((x - plo)*idx)) + ilo;
int i = int(Math::floor((x - plo)*dxinv)) + ilo;
bool inside = i >= ilo && i <= ihi;
return static_cast<T>(inside) - T(0.5);
}, static_cast<T>(tol));
Expand Down Expand Up @@ -217,6 +241,13 @@ public:
return {{AMREX_D_DECL(prob_domain.hi(0),prob_domain.hi(1),prob_domain.hi(2))}};
}

GpuArray<ParticleReal,AMREX_SPACEDIM> ProbLoArrayInParticleReal () const noexcept {
#ifdef AMREX_SINGLE_PRECISION_PARTICLES
return roundoff_lo_f;
#else
return roundoff_lo_d;
#endif
}
GpuArray<ParticleReal,AMREX_SPACEDIM> ProbHiArrayInParticleReal () const noexcept {
#ifdef AMREX_SINGLE_PRECISION_PARTICLES
return roundoff_hi_f;
Expand Down Expand Up @@ -454,11 +485,11 @@ 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_hi_d" and "roundoff_hi_f" each store
// a phi that is very close to that in prob_domain, and for which all doubles and floats less than
// 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_hi_d;
GpuArray<float , AMREX_SPACEDIM> roundoff_hi_f;
GpuArray<double, AMREX_SPACEDIM> roundoff_lo_d, roundoff_hi_d;
GpuArray<float , AMREX_SPACEDIM> roundoff_lo_f, roundoff_hi_f;

//
Box domain;
Expand Down
20 changes: 11 additions & 9 deletions Src/Base/AMReX_Geometry.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -512,33 +512,35 @@ Geometry::computeRoundoffDomain ()
int ihi = Domain().bigEnd(idim);
Real plo = ProbLo(idim);
Real phi = ProbHi(idim);
Real idx = InvCellSize(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);

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

bool
Geometry::outsideRoundoffDomain (AMREX_D_DECL(ParticleReal x, ParticleReal y, ParticleReal z)) const
{
#ifdef AMREX_SINGLE_PRECISION_PARTICLES
bool outside = AMREX_D_TERM(x < prob_domain.lo(0)
bool outside = AMREX_D_TERM(x < roundoff_lo_f[0]
|| x >= roundoff_hi_f[0],
|| y < prob_domain.lo(1)
|| y < roundoff_lo_f[1]
|| y >= roundoff_hi_f[1],
|| z < prob_domain.lo(2)
|| z < roundoff_lo_f[2]
|| z >= roundoff_hi_f[2]);
#else
bool outside = AMREX_D_TERM(x < prob_domain.lo(0)
bool outside = AMREX_D_TERM(x < roundoff_lo_d[0]
|| x >= roundoff_hi_d[0],
|| y < prob_domain.lo(1)
|| y < roundoff_lo_d[1]
|| y >= roundoff_hi_d[1],
|| z < prob_domain.lo(2)
|| z < roundoff_lo_d[2]
|| z >= roundoff_hi_d[2]);
#endif
return outside;
Expand Down
18 changes: 10 additions & 8 deletions Src/Particle/AMReX_ParticleContainerI.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -239,10 +239,11 @@ ParticleContainer<NStructReal, NStructInt, NArrayReal, NArrayInt, Allocator>
const auto& geom = Geom(0);
const auto plo = geom.ProbLoArray();
const auto phi = geom.ProbHiArray();
const auto rlo = geom.ProbLoArrayInParticleReal();
const auto rhi = geom.ProbHiArrayInParticleReal();
const auto is_per = geom.isPeriodicArray();

return enforcePeriodic(p, plo, phi, rhi, is_per);
return enforcePeriodic(p, plo, phi, rlo, rhi, is_per);
}

template <int NStructReal, int NStructInt, int NArrayReal, int NArrayInt,
Expand Down Expand Up @@ -316,15 +317,15 @@ ParticleContainer<NStructReal, NStructInt, NArrayReal, NArrayInt, Allocator>::lo
{
if (Geom(0).outsideRoundoffDomain(AMREX_D_DECL(p.pos(0), p.pos(1), p.pos(2))))
{
RealBox prob_domain = Geom(0).ProbDomain();
GpuArray<ParticleReal, AMREX_SPACEDIM> phi = Geom(0).ProbHiArrayInParticleReal();
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) <= prob_domain.lo(idim)) {
p.pos(idim) = std::nextafter((ParticleReal) prob_domain.lo(idim), phi[idim]);
if (p.pos(idim) <= rlo[idim]) {
p.pos(idim) = std::nextafter(rlo[idim], rhi[idim]);
}
if (p.pos(idim) >= phi[idim]) {
p.pos(idim) = std::nextafter(phi[idim], (ParticleReal) prob_domain.lo(idim));
if (p.pos(idim) >= rhi[idim]) {
p.pos(idim) = std::nextafter(rhi[idim], rlo[idim]);
}
}

Expand Down Expand Up @@ -1250,6 +1251,7 @@ ParticleContainer<NStructReal, NStructInt, NArrayReal, NArrayInt, Allocator>
Vector<std::map<int, int> > new_sizes(num_levels);
const auto plo = Geom(0).ProbLoArray();
const auto phi = Geom(0).ProbHiArray();
const auto rlo = Geom(0).ProbLoArrayInParticleReal();
const auto rhi = Geom(0).ProbHiArrayInParticleReal();
const auto is_per = Geom(0).isPeriodicArray();
for (int lev = lev_min; lev <= finest_lev_particles; ++lev)
Expand All @@ -1271,7 +1273,7 @@ ParticleContainer<NStructReal, NStructInt, NArrayReal, NArrayInt, Allocator>
"perhaps particles have not been initialized correctly?");

int num_stay = partitionParticlesByDest(src_tile, assign_grid, BufferMap(),
plo, phi, rhi, is_per, lev, gid, tid,
plo, phi, rlo, rhi, is_per, lev, gid, tid,
lev_min, lev_max, nGrow, remove_negative);

int num_move = np - num_stay;
Expand Down
10 changes: 7 additions & 3 deletions Src/Particle/AMReX_ParticleUtil.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -517,6 +517,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
bool enforcePeriodic (P& p,
amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& phi,
amrex::GpuArray<amrex::ParticleReal,AMREX_SPACEDIM> const& rlo,
amrex::GpuArray<amrex::ParticleReal,AMREX_SPACEDIM> const& rhi,
amrex::GpuArray<int,AMREX_SPACEDIM> const& is_per) noexcept
{
Expand All @@ -529,7 +530,9 @@ bool enforcePeriodic (P& p,
p.pos(idim) -= static_cast<ParticleReal>(phi[idim] - plo[idim]);
}
// clamp to avoid precision issues;
if (p.pos(idim) < plo[idim]) p.pos(idim) = static_cast<ParticleReal>(plo[idim]);
if (p.pos(idim) < rlo[idim]) {
p.pos(idim) = rlo[idim];
}
shifted = true;
}
else if (p.pos(idim) < plo[idim]) {
Expand All @@ -538,7 +541,7 @@ bool enforcePeriodic (P& p,
}
// clamp to avoid precision issues;
if (p.pos(idim) > rhi[idim]) {
p.pos(idim) = static_cast<ParticleReal>(rhi[idim]);
p.pos(idim) = rhi[idim];
}
shifted = true;
}
Expand All @@ -555,6 +558,7 @@ int
partitionParticlesByDest (PTile& ptile, const PLocator& ploc, const ParticleBufferMap& pmap,
const GpuArray<Real,AMREX_SPACEDIM>& plo,
const GpuArray<Real,AMREX_SPACEDIM>& phi,
const GpuArray<ParticleReal,AMREX_SPACEDIM>& rlo,
const GpuArray<ParticleReal,AMREX_SPACEDIM>& rhi,
const GpuArray<int ,AMREX_SPACEDIM>& is_per,
int lev, int gid, int /*tid*/,
Expand Down Expand Up @@ -602,7 +606,7 @@ partitionParticlesByDest (PTile& ptile, const PLocator& ploc, const ParticleBuff
else
{
auto p_prime = p;
enforcePeriodic(p_prime, plo, phi, rhi, is_per);
enforcePeriodic(p_prime, plo, phi, rlo, rhi, is_per);
auto tup_prime = ploc(p_prime, lev_min, lev_max, nGrow);
assigned_grid = amrex::get<0>(tup_prime);
assigned_lev = amrex::get<1>(tup_prime);
Expand Down

0 comments on commit 826cd37

Please sign in to comment.