From ca437affc6fd00490d8b14e244e53bf641207224 Mon Sep 17 00:00:00 2001 From: Michael B Kuhn <31661049+mbkuhn@users.noreply.github.com> Date: Wed, 16 Oct 2024 15:46:05 -0600 Subject: [PATCH] Correctly populate mac vel ghost cells at high boundaries in auxiliary fills (#1298) + new unit test for code moved to new util --- amr-wind/utilities/index_operations.H | 33 ++- amr-wind/wind_energy/ABLBoundaryPlane.cpp | 6 +- amr-wind/wind_energy/ABLModulatedPowerLaw.cpp | 13 +- unit_tests/core/CMakeLists.txt | 1 + unit_tests/core/test_auxiliary_fill.cpp | 261 ++++++++++++++++++ .../projection/test_pressure_offset.cpp | 32 ++- 6 files changed, 318 insertions(+), 28 deletions(-) create mode 100644 unit_tests/core/test_auxiliary_fill.cpp diff --git a/amr-wind/utilities/index_operations.H b/amr-wind/utilities/index_operations.H index 49182fa36f..587480cc47 100644 --- a/amr-wind/utilities/index_operations.H +++ b/amr-wind/utilities/index_operations.H @@ -36,12 +36,43 @@ AMREX_FORCE_INLINE T perpendicular_idx(const int normal) return T{-1, -1}; } +/** Get the intersection with a boundary box while considering if on face or + * cell. Intended for auxiliary boundary fill calls. + * + * \param grown_interior_box box grown from domain interior to overlap with + * boundary + * \param domain_boundary_box box representing domain boundary containing data + * \param ori orientation of current boundary condition being evaluated + * \return The Box instance representing the intersection of the two inputs + * with shifts added to consider the location of the field associated with the + * grown_interior_box + */ +inline amrex::Box face_aware_boundary_box_intersection( + amrex::Box grown_interior_box, + const amrex::Box& domain_boundary_box, + const amrex::Orientation& ori) +{ + const auto& field_location_vector = grown_interior_box.type(); + if (!grown_interior_box.cellCentered()) { + grown_interior_box.enclosedCells(); + } + auto bx = grown_interior_box & domain_boundary_box; + if (bx.isEmpty()) { + return bx; + } + + if (ori.isHigh() && field_location_vector[ori.coordDir()] == 1) { + bx.shift(field_location_vector); + } + return bx; +} + /** Convert a bounding box into amrex::Box index space at a given level * * \param rbx Bounding box as defined in global domain coordinates * \param geom AMReX geometry information for a given level * \return The Box instance that defines the index space equivalent to bounding - * boxt + * box */ amrex::Box realbox_to_box(const amrex::RealBox& rbx, const amrex::Geometry& geom); diff --git a/amr-wind/wind_energy/ABLBoundaryPlane.cpp b/amr-wind/wind_energy/ABLBoundaryPlane.cpp index c9422739ed..5175cf33d9 100644 --- a/amr-wind/wind_energy/ABLBoundaryPlane.cpp +++ b/amr-wind/wind_energy/ABLBoundaryPlane.cpp @@ -911,11 +911,9 @@ void ABLBoundaryPlane::populate_data( ++mfi) { auto sbx = mfi.growntilebox(1); - if (!sbx.cellCentered()) { - sbx.enclosedCells(); - } const auto& src = m_in_data.interpolate_data(ori, lev); - const auto& bx = sbx & src.box(); + const auto& bx = utils::face_aware_boundary_box_intersection( + sbx, src.box(), ori); if (bx.isEmpty()) { continue; } diff --git a/amr-wind/wind_energy/ABLModulatedPowerLaw.cpp b/amr-wind/wind_energy/ABLModulatedPowerLaw.cpp index f26fe89326..fc31e62e6e 100644 --- a/amr-wind/wind_energy/ABLModulatedPowerLaw.cpp +++ b/amr-wind/wind_energy/ABLModulatedPowerLaw.cpp @@ -4,6 +4,7 @@ #include "AMReX_Gpu.H" #include "AMReX_ParmParse.H" #include "amr-wind/utilities/ncutils/nc_interface.H" +#include "amr-wind/utilities/index_operations.H" #include #include @@ -174,10 +175,8 @@ void ABLModulatedPowerLaw::set_velocity( for (amrex::MFIter mfi(mfab); mfi.isValid(); ++mfi) { auto gbx = amrex::grow(mfi.validbox(), nghost); - if (!gbx.cellCentered()) { - gbx.enclosedCells(); - } - const auto& bx = gbx & dbx; + const auto& bx = + utils::face_aware_boundary_box_intersection(gbx, dbx, ori); if (!bx.ok()) { continue; } @@ -251,10 +250,8 @@ void ABLModulatedPowerLaw::set_temperature( for (amrex::MFIter mfi(mfab); mfi.isValid(); ++mfi) { auto gbx = amrex::grow(mfi.validbox(), nghost); - if (!gbx.cellCentered()) { - gbx.enclosedCells(); - } - const auto& bx = gbx & dbx; + const auto& bx = + utils::face_aware_boundary_box_intersection(gbx, dbx, ori); if (!bx.ok()) { continue; } diff --git a/unit_tests/core/CMakeLists.txt b/unit_tests/core/CMakeLists.txt index 6d7173f584..baad472bab 100644 --- a/unit_tests/core/CMakeLists.txt +++ b/unit_tests/core/CMakeLists.txt @@ -7,6 +7,7 @@ target_sources( test_field_ops.cpp test_field_fillpatch_ops.cpp test_physics.cpp + test_auxiliary_fill.cpp ) add_subdirectory(vs) diff --git a/unit_tests/core/test_auxiliary_fill.cpp b/unit_tests/core/test_auxiliary_fill.cpp new file mode 100644 index 0000000000..a5b9c7efd3 --- /dev/null +++ b/unit_tests/core/test_auxiliary_fill.cpp @@ -0,0 +1,261 @@ +#include "aw_test_utils/MeshTest.H" +#include "amr-wind/utilities/index_operations.H" + +namespace amr_wind_tests { + +namespace { + +void auxiliary_fill_boundary(amr_wind::Field& velocity, const int comp = 0) +{ + const int nlevels = velocity.repo().num_active_levels(); + const int ncomp = velocity.num_comp(); + + for (int lev = 0; lev < nlevels; ++lev) { + + for (amrex::OrientationIter oit; oit != nullptr; ++oit) { + auto ori = oit(); + const amrex::Box& minBox = + velocity.repo().mesh().boxArray(lev).minimalBox(); + amrex::IntVect plo(minBox.loVect()); + amrex::IntVect phi(minBox.hiVect()); + const int normal = ori.coordDir(); + plo[normal] = ori.isHigh() ? minBox.hiVect()[normal] + 1 : -1; + phi[normal] = ori.isHigh() ? minBox.hiVect()[normal] + 1 : -1; + const amrex::Box domain_bdy_bx(plo, phi); + + for (amrex::MFIter mfi(velocity(lev), amrex::TilingIfNotGPU()); + mfi.isValid(); ++mfi) { + + auto sbx = mfi.growntilebox(1); + const auto& bx = + amr_wind::utils::face_aware_boundary_box_intersection( + sbx, domain_bdy_bx, ori); + if (bx.isEmpty()) { + continue; + } + + const auto& dest = velocity(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(comp + n + 1); + }); + } + } + } +} + +amrex::Real get_field_err( + amr_wind::Field& field, 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(), + [=] AMREX_GPU_HOST_DEVICE( + amrex::Box const& bx, + amrex::Array4 const& f_arr) -> amrex::Real { + amrex::Real error = 0; + + amrex::Loop(bx, [=, &error](int i, int j, int k) noexcept { + // Do indices manually to check against box operations in code + if (ncomp == 1) { + // Checking a MAC velocity + const bool in_x_bdy = (i == -1 || i == 8); + const bool not_in_yz_bdy = + (j >= 0 && j <= 7) && (k >= 0 && k <= 7); + const bool in_y_bdy = (j == -1 || j == 8); + const bool not_in_xz_bdy = + (i >= 0 && i <= 7) && (k >= 0 && k <= 7); + const bool in_z_bdy = (k == -1 || k == 8); + const bool not_in_xy_bdy = + (j >= 0 && j <= 7) && (i >= 0 && i <= 7); + if (comp == 0) { + // Checking u + if (check_all_ghost) { + // Ghost cells in boundaries, but not corners + const bool in_xf_bdy = (i == -1 || i == 9); + if ((in_xf_bdy && not_in_yz_bdy) || + (in_y_bdy && not_in_xz_bdy) || + (in_z_bdy && not_in_xy_bdy)) { + error += std::abs(f_arr(i, j, k) - 1.0); + } else { + error += std::abs(f_arr(i, j, k) - 0.0); + } + } else { + // Valid face-normal cells + if ((i == 0 || i == 8) && (j >= 0 && j <= 7) && + (k >= 0 && k <= 7)) { + error += std::abs(f_arr(i, j, k) - 1.0); + } else { + error += std::abs(f_arr(i, j, k) - 0.0); + } + } + } else if (comp == 1) { + // Checking v + if (check_all_ghost) { + // Ghost cells in boundaries, but not corners + const bool in_yf_bdy = (j == -1 || j == 9); + if ((in_x_bdy && not_in_yz_bdy) || + (in_yf_bdy && not_in_xz_bdy) || + (in_z_bdy && not_in_xy_bdy)) { + error += std::abs(f_arr(i, j, k) - 2.0); + } else { + error += std::abs(f_arr(i, j, k) - 0.0); + } + } else { + // Valid face-normal cells + if ((j == 0 || j == 8) && (i >= 0 && i <= 7) && + (k >= 0 && k <= 7)) { + error += std::abs(f_arr(i, j, k) - 2.0); + } else { + error += std::abs(f_arr(i, j, k) - 0.0); + } + } + } else { + // Checking w + if (check_all_ghost) { + // Ghost cells in boundaries, but not corners + const bool in_zf_bdy = (k == -1 || k == 9); + if ((in_x_bdy && not_in_yz_bdy) || + (in_y_bdy && not_in_xz_bdy) || + (in_zf_bdy && not_in_xy_bdy)) { + error += std::abs(f_arr(i, j, k) - 3.0); + } else { + error += std::abs(f_arr(i, j, k) - 0.0); + } + } else { + // Valid face-normal cells + if ((k == 0 || k == 8) && (i >= 0 && i <= 7) && + (j >= 0 && j <= 7)) { + error += std::abs(f_arr(i, j, k) - 3.0); + } else { + error += std::abs(f_arr(i, j, k) - 0.0); + } + } + } + } else { + // Checking cell-centered velocity + if (check_all_ghost) { + // Ghost cells in boundaries + if (i <= -1 || j <= -1 || k <= -1 || i >= 8 || j >= 8 || + k >= 8) { + 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); + } else { + error += std::abs(f_arr(i, j, k, 0) - 0.0); + error += std::abs(f_arr(i, j, k, 1) - 0.0); + error += std::abs(f_arr(i, j, k, 2) - 0.0); + } + } else { + // First ghost cells only + const bool is_ibdy = + ((i == -1 || i == 8) && (j >= 0 && j <= 7) && + (k >= 0 && k <= 7)); + const bool is_jbdy = + ((j == -1 || j == 8) && (i >= 0 && i <= 7) && + (k >= 0 && k <= 7)); + const bool is_kbdy = + ((k == -1 || k == 8) && (i >= 0 && i <= 7) && + (j >= 0 && j <= 7)); + if (is_ibdy || is_jbdy || is_kbdy) { + 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); + } else { + error += std::abs(f_arr(i, j, k, 0) - 0.0); + error += std::abs(f_arr(i, j, k, 1) - 0.0); + error += std::abs(f_arr(i, j, k, 2) - 0.0); + } + } + } + }); + + return error; + }); + amrex::ParallelDescriptor::ReduceRealSum(error_total); + return error_total; +} +} // namespace + +class AuxiliaryFillTest : public MeshTest +{ +public: + void set_up_fields_cc() + { + auto& frepo = mesh().field_repo(); + m_vel = &frepo.declare_field("velocity", 3, 1, 1); + + (*m_vel).setVal(0.); + } + + void set_up_fields_face() + { + auto& frepo = mesh().field_repo(); + auto mac_vels = frepo.declare_face_normal_field( + {"u_mac", "v_mac", "w_mac"}, 1, 1, 1); + m_umac = mac_vels[0]; + m_vmac = mac_vels[1]; + m_wmac = mac_vels[2]; + + (*m_umac).setVal(0.); + (*m_vmac).setVal(0.); + (*m_wmac).setVal(0.); + } + + void prep_test(const bool is_cc) + { + // Default dimensions are n_cell = 8 x 8 x 8 + MeshTest::populate_parameters(); + { + amrex::Vector periodic{{0, 0, 0}}; + amrex::ParmParse pp("geometry"); + pp.addarr("is_periodic", periodic); + } + initialize_mesh(); + if (is_cc) { + set_up_fields_cc(); + } else { + set_up_fields_face(); + } + } + + amr_wind::Field* m_vel; + amr_wind::Field* m_umac; + amr_wind::Field* m_vmac; + amr_wind::Field* m_wmac; +}; + +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); + EXPECT_DOUBLE_EQ(err, 0.); +} + +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); + EXPECT_DOUBLE_EQ(u_err, 0.); + + auxiliary_fill_boundary(*m_vmac, 1); + const auto v_err = get_field_err(*m_vmac, 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); + EXPECT_DOUBLE_EQ(w_err, 0.); +} + +} // namespace amr_wind_tests diff --git a/unit_tests/projection/test_pressure_offset.cpp b/unit_tests/projection/test_pressure_offset.cpp index cf6cc692a8..5a5838c942 100644 --- a/unit_tests/projection/test_pressure_offset.cpp +++ b/unit_tests/projection/test_pressure_offset.cpp @@ -10,16 +10,14 @@ void init_vel_z(amr_wind::Field& vel, const amrex::Real w_const) const int nlevels = vel.repo().num_active_levels(); for (int lev = 0; lev < nlevels; ++lev) { - - for (amrex::MFIter mfi(vel(lev)); mfi.isValid(); ++mfi) { - auto gbx = mfi.growntilebox(); - const auto& varr = vel(lev).array(mfi); - - amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { - varr(i, j, k, 2) = w_const; + const auto& varrs = vel(lev).arrays(); + const amrex::IntVect ngs = vel.num_grow(); + amrex::ParallelFor( + vel(lev), ngs, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) { + varrs[nbx](i, j, k, 2) = w_const; }); - } } + amrex::Gpu::synchronize(); } void init_ref_p( @@ -33,21 +31,21 @@ void init_ref_p( for (int lev = 0; lev < nlevels; ++lev) { const auto& dx = geom[lev].CellSizeArray(); const auto& probhi = geom[lev].ProbHiArray(); - for (amrex::MFIter mfi(ref_p(lev)); mfi.isValid(); ++mfi) { - auto nbx = mfi.nodaltilebox(); - const auto& p0_arr = ref_p(lev).array(mfi); - - amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + const auto& p0_arrs = ref_p(lev).arrays(); + const amrex::IntVect ngs(0); + amrex::ParallelFor( + ref_p(lev), ngs, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) { // Height of pressure node const amrex::Real hnode = k * dx[2]; // Integrated density from top const amrex::Real irho = rho_0 * (probhi[2] - hnode); // Multiply with force to get hydrostatic pressure - p0_arr(i, j, k) = -irho * F_g; + p0_arrs[nbx](i, j, k) = -irho * F_g; }); - } } + amrex::Gpu::synchronize(); } amrex::Real get_pbottom(amr_wind::Field& pressure) @@ -85,9 +83,13 @@ void ptest_kernel( my_incflo.init_mesh(); auto& density = my_incflo.sim().repo().get_field("density"); auto& velocity = my_incflo.sim().repo().get_field("velocity"); + auto& gp = my_incflo.sim().repo().get_field("gp"); // Set uniform density density.setVal(rho_0); + // Zero pressure gradient + gp.setVal(0.0); // Set velocity as it would be with gravity forcing + velocity.setVal(0.0); init_vel_z(velocity, w_0); // If requested, form reference_pressure field