Skip to content

Commit

Permalink
Fix global variable read on GPU error (Exawind#621)
Browse files Browse the repository at this point in the history
  • Loading branch information
@jrood-nrel
jrood-nrel authored May 9, 2022
1 parent 89d7ec0 commit 2b176ca
Show file tree
Hide file tree
Showing 2 changed files with 69 additions and 54 deletions.
4 changes: 2 additions & 2 deletions amr-wind/utilities/sampling/FreeSurface.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,8 +109,8 @@ private:
//! (this also determines the meaning of start and end points)
int m_coorddir{2};
//! Grid coordinates, determined as a function of m_coorddir
int gc1 = 0;
int gc2 = 1;
int m_gc1 = 0;
int m_gc2 = 1;

//! Parameters to set up plane
amrex::Vector<amrex::Real> m_start, m_end;
Expand Down
119 changes: 67 additions & 52 deletions amr-wind/utilities/sampling/FreeSurface.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,18 +39,18 @@ void FreeSurface::initialize()

switch (m_coorddir) {
case 0: {
gc1 = 1;
gc2 = 2;
m_gc1 = 1;
m_gc2 = 2;
break;
}
case 1: {
gc1 = 0;
gc2 = 2;
m_gc1 = 0;
m_gc2 = 2;
break;
}
case 2: {
gc1 = 0;
gc2 = 1;
m_gc1 = 0;
m_gc2 = 1;
break;
}
default: {
Expand All @@ -71,8 +71,8 @@ void FreeSurface::initialize()

// Get size of sample grid spacing
amrex::Vector<amrex::Real> dx = {0.0, 0.0};
dx[0] = (m_end[gc1] - m_start[gc1]) / amrex::max(m_npts_dir[0] - 1, 1);
dx[1] = (m_end[gc2] - m_start[gc2]) / amrex::max(m_npts_dir[1] - 1, 1);
dx[0] = (m_end[m_gc1] - m_start[m_gc1]) / amrex::max(m_npts_dir[0] - 1, 1);
dx[1] = (m_end[m_gc2] - m_start[m_gc2]) / amrex::max(m_npts_dir[1] - 1, 1);

// Store locations
int idx = 0;
Expand All @@ -83,9 +83,9 @@ void FreeSurface::initialize()
m_out[idx * m_ninst + ni] = m_start[m_coorddir];
}
// Grid direction 1
m_locs[idx][0] = m_start[gc1] + dx[0] * i;
m_locs[idx][0] = m_start[m_gc1] + dx[0] * i;
// Grid direction 2
m_locs[idx][1] = m_start[gc2] + dx[1] * j;
m_locs[idx][1] = m_start[m_gc2] + dx[1] * j;

++idx;
}
Expand Down Expand Up @@ -114,7 +114,7 @@ void FreeSurface::post_advance_work()
// Set up device vector of outputs, initialize to above phi0
const auto& plo0 = m_sim.mesh().Geom(0).ProbLoArray();
const auto& phi0 = m_sim.mesh().Geom(0).ProbHiArray();
amrex::Gpu::DeviceVector<double> dout(m_npts, phi0[2] + 1.0);
amrex::Gpu::DeviceVector<amrex::Real> dout(m_npts, phi0[2] + 1.0);
auto* dout_ptr = dout.data();

// Sum of interface locations at each point (assumes one interface only)
Expand Down Expand Up @@ -147,6 +147,10 @@ void FreeSurface::post_advance_work()
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> plo =
geom.ProbLoArray();

const int captured_m_coorddir = m_coorddir;
const int captured_m_gc1 = m_gc1;
const int captured_m_gc2 = m_gc2;

// Loop points in 2D grid
for (int n = 0; n < m_npts; ++n) {
amrex::GpuArray<amrex::Real, 2> loc;
Expand All @@ -169,21 +173,24 @@ void FreeSurface::post_advance_work()
bx,
[=, &height_fab](int i, int j, int k) noexcept {
// Initialize height measurement
amrex::Real ht = plo[m_coorddir];
amrex::Real ht = plo[captured_m_coorddir];
// Cell location
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>
xm;
xm[0] = plo[0] + (i + 0.5) * dx[0];
xm[1] = plo[1] + (j + 0.5) * dx[1];
xm[2] = plo[2] + (k + 0.5) * dx[2];
int ip = (m_coorddir == 0 ? 1 : 0);
int im = i - ip;
int ip = static_cast<int>(
captured_m_coorddir == 0);
const int im = i - ip;
ip = i + ip;
int jp = (m_coorddir == 1 ? 1 : 0);
int jm = j - jp;
int jp = static_cast<int>(
captured_m_coorddir == 1);
const int jm = j - jp;
jp = j + jp;
int kp = (m_coorddir == 2 ? 1 : 0);
int km = k - kp;
int kp = static_cast<int>(
captured_m_coorddir == 2);
const int km = k - kp;
kp = k + kp;
// (1) Check that cell height is below
// previous instance. (2) Check if cell
Expand All @@ -194,18 +201,22 @@ void FreeSurface::post_advance_work()
// multiphase, then check if cell might have
// interface at top or bottom
if ((dout_ptr[n] >
xm[m_coorddir] +
0.5 * dx[m_coorddir]) &&
(((plo[gc1] == loc[0] &&
xm[gc1] - loc[0] == 0.5 * dx[gc1]) ||
(xm[gc1] - loc[0] < 0.5 * dx[gc1] &&
loc[0] - xm[gc1] <=
0.5 * dx[gc1])) &&
((plo[gc2] == loc[1] &&
xm[gc2] - loc[1] == 0.5 * dx[gc2]) ||
(xm[gc2] - loc[1] < 0.5 * dx[gc2] &&
loc[1] - xm[gc2] <=
0.5 * dx[gc2]))) &&
xm[captured_m_coorddir] +
0.5 * dx[captured_m_coorddir]) &&
(((plo[captured_m_gc1] == loc[0] &&
xm[captured_m_gc1] - loc[0] ==
0.5 * dx[captured_m_gc1]) ||
(xm[captured_m_gc1] - loc[0] <
0.5 * dx[captured_m_gc1] &&
loc[0] - xm[captured_m_gc1] <=
0.5 * dx[captured_m_gc1])) &&
((plo[captured_m_gc2] == loc[1] &&
xm[captured_m_gc2] - loc[1] ==
0.5 * dx[captured_m_gc2]) ||
(xm[captured_m_gc2] - loc[1] <
0.5 * dx[captured_m_gc2] &&
loc[1] - xm[captured_m_gc2] <=
0.5 * dx[captured_m_gc2]))) &&
((vof_arr(i, j, k) < (1.0 - 1e-12) &&
vof_arr(i, j, k) > 1e-12) ||
(vof_arr(i, j, k) < 1e-12 &&
Expand All @@ -227,7 +238,7 @@ void FreeSurface::post_advance_work()
int kdn = k;

// Determine which cells to use for grid
if (m_coorddir != 0) {
if (captured_m_coorddir != 0) {
// If x is a grid coord, it is
// always first (e.g., xy or xz)
int li = 0;
Expand All @@ -244,9 +255,10 @@ void FreeSurface::post_advance_work()
(loc[li] - xm[0]) * dxi[0];
}
}
if (m_coorddir != 1) {
if (captured_m_coorddir != 1) {
// y can be first or second (xy, yz)
int li = (gc1 == 1 ? 0 : 1);
int li =
(captured_m_gc1 == 1 ? 0 : 1);
if (loc[li] < xm[1]) {
jup = j;
jdn = j - 1;
Expand All @@ -260,7 +272,7 @@ void FreeSurface::post_advance_work()
(loc[li] - xm[1]) * dxi[1];
}
}
if (m_coorddir != 2) {
if (captured_m_coorddir != 2) {
// If z is a grid coord, it is
// always second (e.g., yz or xz)
int li = 1;
Expand All @@ -277,15 +289,15 @@ void FreeSurface::post_advance_work()
(loc[li] - xm[2]) * dxi[2];
}
}
amrex::Real wx_lo = 1.0 - wx_hi;
amrex::Real wy_lo = 1.0 - wy_hi;
amrex::Real wz_lo = 1.0 - wz_hi;
const amrex::Real wx_lo = 1.0 - wx_hi;
const amrex::Real wy_lo = 1.0 - wy_hi;
const amrex::Real wz_lo = 1.0 - wz_hi;

amrex::Real vof_above = 0.0;
amrex::Real vof_below = 0.0;
amrex::Real vof_here = 0.0;

if (m_coorddir == 0) {
if (captured_m_coorddir == 0) {
vof_above =
wz_lo * wy_lo *
vof_arr(i + 1, jdn, kdn) +
Expand Down Expand Up @@ -314,7 +326,7 @@ void FreeSurface::post_advance_work()
wz_hi * wy_hi *
vof_arr(i - 1, jup, kup);
}
if (m_coorddir == 1) {
if (captured_m_coorddir == 1) {
vof_above =
wx_lo * wz_lo *
vof_arr(idn, j + 1, kdn) +
Expand Down Expand Up @@ -343,7 +355,7 @@ void FreeSurface::post_advance_work()
wx_hi * wz_hi *
vof_arr(iup, j - 1, kup);
}
if (m_coorddir == 2) {
if (captured_m_coorddir == 2) {
vof_above =
wx_lo * wy_lo *
vof_arr(idn, jdn, k + 1) +
Expand Down Expand Up @@ -374,26 +386,28 @@ void FreeSurface::post_advance_work()
}
// Determine which cell to
// interpolate with
bool above = (vof_above - 0.5) *
(vof_here - 0.5) <=
0.0;
bool below = (vof_below - 0.5) *
(vof_here - 0.5) <=
0.0;
const bool above =
(vof_above - 0.5) *
(vof_here - 0.5) <=
0.0;
const bool below =
(vof_below - 0.5) *
(vof_here - 0.5) <=
0.0;
if (above) {
// Interpolate positive
// direction
ht = xm[m_coorddir] +
(dx[m_coorddir]) /
ht = xm[captured_m_coorddir] +
(dx[captured_m_coorddir]) /
(vof_above - vof_here) *
(0.5 - vof_here);
} else {
if (below) {
// Interpolate negative
// direction
ht =
xm[m_coorddir] -
(dx[m_coorddir]) /
xm[captured_m_coorddir] -
(dx[captured_m_coorddir]) /
(vof_below - vof_here) *
(0.5 - vof_here);
}
Expand All @@ -405,8 +419,9 @@ void FreeSurface::post_advance_work()
// Offset by removing lo and contribute to
// whole
height_fab = amrex::max(
height_fab, mask_arr(i, j, k) *
(ht - plo[m_coorddir]));
height_fab,
mask_arr(i, j, k) *
(ht - plo[captured_m_coorddir]));
});
return height_fab;
}));
Expand Down

0 comments on commit 2b176ca

Please sign in to comment.