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

Jump to bottom

Add roundoff_lo corresponding to roundoff_hi for domains that don't start at 0 #2950

Merged
merged 5 commits into from
Sep 16, 2022
Merged

Add roundoff_lo corresponding to roundoff_hi for domains that don't start at 0 #2950

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
1c8abb0
Lay groundwork for roundoff_lo
PhilMiller Sep 15, 2022
60a1a5c
Add dummy implementation of roundoff_lo computation
PhilMiller Sep 15, 2022
58ff22f
implement bisect_prob_lo
atmyers Sep 15, 2022
9f161cb
change idx -> dxinv
atmyers Sep 15, 2022
4ee6f83
use rlo instead of plo in locateParticle
atmyers Sep 15, 2022
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
47 changes: 39 additions & 8 deletions Src/Base/AMReX_Geometry.H
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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 @@ -1242,6 +1243,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 @@ -1263,7 +1265,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
Open in desktop
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