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

Avoid OOB for hi boundary planes #1307

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
+104 −17
Open

Avoid OOB for hi boundary planes #1307

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2bad64f
include shift for hi boundaries, abl planes
mbkuhn Oct 21, 2024
9c20fe0
unit test shift indices of mac velocities
mbkuhn Oct 21, 2024
e2a1504
unit test indices for cell centers and ensure we start with bad values.
mbkuhn Oct 21, 2024
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
16 changes: 16 additions & 0 deletions amr-wind/utilities/index_operations.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,8 @@ AMREX_FORCE_INLINE T perpendicular_idx(const int normal)
/** Get the intersection with a boundary box while considering if on face or
* cell. Intended for auxiliary boundary fill calls.
*
* \param shift_to_cc output index vector to shift from output box to
* a cell-centered representation.
* \param grown_interior_box box grown from domain interior to overlap with
* boundary
* \param domain_boundary_box box representing domain boundary containing data
Expand All @@ -48,10 +50,12 @@ AMREX_FORCE_INLINE T perpendicular_idx(const int normal)
* grown_interior_box
*/
inline amrex::Box face_aware_boundary_box_intersection(
amrex::IntVect& shift_to_cc,
amrex::Box grown_interior_box,
const amrex::Box& domain_boundary_box,
const amrex::Orientation& ori)
{
shift_to_cc = {0, 0, 0};
const auto& field_location_vector = grown_interior_box.type();
if (!grown_interior_box.cellCentered()) {
grown_interior_box.enclosedCells();
Expand All @@ -63,10 +67,22 @@ inline amrex::Box face_aware_boundary_box_intersection(

if (ori.isHigh() && field_location_vector[ori.coordDir()] == 1) {
bx.shift(field_location_vector);
shift_to_cc = -field_location_vector;
}
return bx;
}

// This version omits the shift to cell-centered data, often not needed
inline amrex::Box face_aware_boundary_box_intersection(
amrex::Box grown_interior_box,
const amrex::Box& domain_boundary_box,
const amrex::Orientation& ori)
{
amrex::IntVect shift_to_cc = {0, 0, 0};
return face_aware_boundary_box_intersection(
shift_to_cc, grown_interior_box, domain_boundary_box, ori);
}

/** Convert a bounding box into amrex::Box index space at a given level
*
* \param rbx Bounding box as defined in global domain coordinates
Expand Down
8 changes: 5 additions & 3 deletions amr-wind/wind_energy/ABLBoundaryPlane.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -922,8 +922,9 @@ void ABLBoundaryPlane::populate_data(

auto sbx = mfi.growntilebox(1);
const auto& src = m_in_data.interpolate_data(ori, lev);
auto shift_to_cc = amrex::IntVect(0);
const auto& bx = utils::face_aware_boundary_box_intersection(
sbx, src.box(), ori);
shift_to_cc, sbx, src.box(), ori);
if (bx.isEmpty()) {
continue;
}
Expand All @@ -934,8 +935,9 @@ void ABLBoundaryPlane::populate_data(
amrex::ParallelFor(
bx, nc,
[=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
dest(i, j, k, n + dcomp) =
src_arr(i, j, k, n + nstart + orig_comp);
dest(i, j, k, n + dcomp) = src_arr(
i + shift_to_cc[0], j + shift_to_cc[1],
k + shift_to_cc[2], n + nstart + orig_comp);
});
}
}
Expand Down
97 changes: 83 additions & 14 deletions unit_tests/core/test_auxiliary_fill.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,8 @@ namespace amr_wind_tests {

namespace {

void auxiliary_fill_boundary(amr_wind::Field& velocity, const int comp = 0)
void auxiliary_fill_boundary(
amr_wind::Field& velocity, amr_wind::IntField& indices, const int comp = 0)
{
const int nlevels = velocity.repo().num_active_levels();
const int ncomp = velocity.num_comp();
Expand All @@ -27,37 +28,51 @@ void auxiliary_fill_boundary(amr_wind::Field& velocity, const int comp = 0)
mfi.isValid(); ++mfi) {

auto sbx = mfi.growntilebox(1);
auto shift_to_cc = amrex::IntVect(0);
const auto& bx =
amr_wind::utils::face_aware_boundary_box_intersection(
sbx, domain_bdy_bx, ori);
shift_to_cc, sbx, domain_bdy_bx, ori);
if (bx.isEmpty()) {
continue;
}

const auto& dest = velocity(lev).array(mfi);
const auto& idx = indices(lev).array(mfi);
amrex::ParallelFor(
bx, ncomp,
[=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
dest(i, j, k, n) =
static_cast<amrex::Real>(comp + n + 1);
idx(i + shift_to_cc[0], j + shift_to_cc[1],
k + shift_to_cc[2], 3 * comp + 3 * n) =
shift_to_cc[0];
idx(i + shift_to_cc[0], j + shift_to_cc[1],
k + shift_to_cc[2], 3 * comp + 3 * n + 1) =
shift_to_cc[1];
idx(i + shift_to_cc[0], j + shift_to_cc[1],
k + shift_to_cc[2], 3 * comp + 3 * n + 2) =
shift_to_cc[2];
});
}
}
}
}

amrex::Real get_field_err(
amr_wind::Field& field, const bool check_all_ghost, const int comp = 0)
amr_wind::Field& field,
amr_wind::IntField& indices,
const bool check_all_ghost,
const int comp = 0)
{
const int lev = 0;
amrex::Real error_total = 0;
const int ncomp = field.num_comp();

error_total += amrex::ReduceSum(
field(lev), field(lev).nGrow(),
field(lev), indices(lev), field(lev).nGrow(),
[=] AMREX_GPU_HOST_DEVICE(
amrex::Box const& bx,
amrex::Array4<amrex::Real const> const& f_arr) -> amrex::Real {
amrex::Box const& bx, amrex::Array4<amrex::Real const> const& f_arr,
amrex::Array4<int const> const& i_arr) -> amrex::Real {
amrex::Real error = 0;

amrex::Loop(bx, [=, &error](int i, int j, int k) noexcept {
Expand All @@ -82,6 +97,21 @@ amrex::Real get_field_err(
(in_y_bdy && not_in_xz_bdy) ||
(in_z_bdy && not_in_xy_bdy)) {
error += std::abs(f_arr(i, j, k) - 1.0);
if (i == 9) {
error += std::abs((amrex::Real)(
i_arr(i - 1, j, k, 0) + 1));
error += std::abs((amrex::Real)(
i_arr(i - 1, j, k, 1) - 0));
error += std::abs((amrex::Real)(
i_arr(i - 1, j, k, 2) - 0));
} else if (i == -1) {
error += std::abs(
(amrex::Real)(i_arr(i, j, k, 0) - 0));
error += std::abs(
(amrex::Real)(i_arr(i, j, k, 1) - 0));
error += std::abs(
(amrex::Real)(i_arr(i, j, k, 2) - 0));
}
} else {
error += std::abs(f_arr(i, j, k) - 0.0);
}
Expand All @@ -103,6 +133,21 @@ amrex::Real get_field_err(
(in_yf_bdy && not_in_xz_bdy) ||
(in_z_bdy && not_in_xy_bdy)) {
error += std::abs(f_arr(i, j, k) - 2.0);
if (j == 9) {
error += std::abs((amrex::Real)(
i_arr(i, j - 1, k, 3) - 0));
error += std::abs((amrex::Real)(
i_arr(i, j - 1, k, 4) + 1));
error += std::abs((amrex::Real)(
i_arr(i, j - 1, k, 5) - 0));
} else if (j == -1) {
error += std::abs(
(amrex::Real)(i_arr(i, j, k, 3) - 0));
error += std::abs(
(amrex::Real)(i_arr(i, j, k, 4) - 0));
error += std::abs(
(amrex::Real)(i_arr(i, j, k, 5) - 0));
}
} else {
error += std::abs(f_arr(i, j, k) - 0.0);
}
Expand All @@ -124,6 +169,21 @@ amrex::Real get_field_err(
(in_y_bdy && not_in_xz_bdy) ||
(in_zf_bdy && not_in_xy_bdy)) {
error += std::abs(f_arr(i, j, k) - 3.0);
if (k == 9) {
error += std::abs((amrex::Real)(
i_arr(i, j, k - 1, 6) - 0));
error += std::abs((amrex::Real)(
i_arr(i, j, k - 1, 7) - 0));
error += std::abs((amrex::Real)(
i_arr(i, j, k - 1, 8) + 1));
} else if (k == -1) {
error += std::abs(
(amrex::Real)(i_arr(i, j, k, 6) - 0));
error += std::abs(
(amrex::Real)(i_arr(i, j, k, 7) - 0));
error += std::abs(
(amrex::Real)(i_arr(i, j, k, 8) - 0));
}
} else {
error += std::abs(f_arr(i, j, k) - 0.0);
}
Expand Down Expand Up @@ -166,6 +226,10 @@ amrex::Real get_field_err(
error += std::abs(f_arr(i, j, k, 0) - 1.0);
error += std::abs(f_arr(i, j, k, 1) - 2.0);
error += std::abs(f_arr(i, j, k, 2) - 3.0);
for (int ni = 0; ni < 9; ++ni) {
error += std::abs(
(amrex::Real)(i_arr(i, j, k, ni) - 0));
}
} else {
error += std::abs(f_arr(i, j, k, 0) - 0.0);
error += std::abs(f_arr(i, j, k, 1) - 0.0);
Expand All @@ -189,8 +253,10 @@ class AuxiliaryFillTest : public MeshTest
{
auto& frepo = mesh().field_repo();
m_vel = &frepo.declare_field("velocity", 3, 1, 1);
m_ind = &frepo.declare_int_field("indices", 9, 1);

(*m_vel).setVal(0.);
(*m_ind).setVal(5);
}

void set_up_fields_face()
Expand All @@ -201,10 +267,12 @@ class AuxiliaryFillTest : public MeshTest
m_umac = mac_vels[0];
m_vmac = mac_vels[1];
m_wmac = mac_vels[2];
m_ind = &frepo.declare_int_field("indices", 9, 1);

(*m_umac).setVal(0.);
(*m_vmac).setVal(0.);
(*m_wmac).setVal(0.);
(*m_ind).setVal(5);
}

void prep_test(const bool is_cc)
Expand All @@ -228,15 +296,16 @@ class AuxiliaryFillTest : public MeshTest
amr_wind::Field* m_umac;
amr_wind::Field* m_vmac;
amr_wind::Field* m_wmac;
amr_wind::IntField* m_ind;
};

TEST_F(AuxiliaryFillTest, velocity_cc)
{
prep_test(true);

// Do fill and check ghost cells
auxiliary_fill_boundary(*m_vel);
const auto err = get_field_err(*m_vel, false);
auxiliary_fill_boundary(*m_vel, *m_ind);
const auto err = get_field_err(*m_vel, *m_ind, false);
EXPECT_DOUBLE_EQ(err, 0.);
}

Expand All @@ -245,16 +314,16 @@ TEST_F(AuxiliaryFillTest, velocity_face)
prep_test(false);

// Do fill and check ghost cells
auxiliary_fill_boundary(*m_umac, 0);
const auto u_err = get_field_err(*m_umac, true, 0);
auxiliary_fill_boundary(*m_umac, *m_ind, 0);
const auto u_err = get_field_err(*m_umac, *m_ind, true, 0);
EXPECT_DOUBLE_EQ(u_err, 0.);

auxiliary_fill_boundary(*m_vmac, 1);
const auto v_err = get_field_err(*m_vmac, true, 1);
auxiliary_fill_boundary(*m_vmac, *m_ind, 1);
const auto v_err = get_field_err(*m_vmac, *m_ind, true, 1);
EXPECT_DOUBLE_EQ(v_err, 0.);

auxiliary_fill_boundary(*m_wmac, 2);
const auto w_err = get_field_err(*m_wmac, true, 2);
auxiliary_fill_boundary(*m_wmac, *m_ind, 2);
const auto w_err = get_field_err(*m_wmac, *m_ind, true, 2);
EXPECT_DOUBLE_EQ(w_err, 0.);
}

Expand Down