Skip to content

Commit

Permalink
Diagnostics bug fix + de-hardcoding refRatio in makeFineMask calls (#…
Browse files Browse the repository at this point in the history 
…1050)

* unit test for multilevel diagnostics
   - this test detects out-of-bounds when compiled in debug
* avoid out-of-bounds and fix makeFineMesh command in PrintMaxMACVelLocations
* replacing hard-coded refratio IntVect(2) in makeFineMesh calls
  • Loading branch information
@mbkuhn
mbkuhn authored May 10, 2024
1 parent 083f42b commit 3edc314
Show file tree
Hide file tree
Showing 16 changed files with 175 additions and 28 deletions.
2 changes: 1 addition & 1 deletion amr-wind/physics/BurggrafFlow.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,7 +121,7 @@ amrex::Real BurggrafFlow::compute_error(const Field& field)
if (lev < nlevels - 1) {
level_mask = makeFineMask(
m_mesh.boxArray(lev), m_mesh.DistributionMap(lev),
m_mesh.boxArray(lev + 1), amrex::IntVect(2), 1, 0);
m_mesh.boxArray(lev + 1), m_mesh.refRatio(lev), 1, 0);
} else {
level_mask.define(
m_mesh.boxArray(lev), m_mesh.DistributionMap(lev), 1, 0,
Expand Down
4 changes: 2 additions & 2 deletions amr-wind/physics/ChannelFlow.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -271,7 +271,7 @@ amrex::Real ChannelFlow::compute_error(const IndexSelector& idxOp)
if (lev < nlevels - 1) {
level_mask = makeFineMask(
m_mesh.boxArray(lev), m_mesh.DistributionMap(lev),
m_mesh.boxArray(lev + 1), amrex::IntVect(2), 1, 0);
m_mesh.boxArray(lev + 1), m_mesh.refRatio(lev), 1, 0);
} else {
level_mask.define(
m_mesh.boxArray(lev), m_mesh.DistributionMap(lev), 1, 0,
Expand Down Expand Up @@ -352,7 +352,7 @@ amrex::Real ChannelFlow::compute_analytical_smagorinsky_error()
if (lev < nlevels - 1) {
level_mask = makeFineMask(
m_mesh.boxArray(lev), m_mesh.DistributionMap(lev),
m_mesh.boxArray(lev + 1), amrex::IntVect(2), 1, 0);
m_mesh.boxArray(lev + 1), m_mesh.refRatio(lev), 1, 0);
} else {
level_mask.define(
m_mesh.boxArray(lev), m_mesh.DistributionMap(lev), 1, 0,
Expand Down
2 changes: 1 addition & 1 deletion amr-wind/physics/ConvectingTaylorVortex.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -231,7 +231,7 @@ amrex::Real ConvectingTaylorVortex::compute_error(const Field& field)
if (lev < nlevels - 1) {
level_mask = makeFineMask(
m_mesh.boxArray(lev), m_mesh.DistributionMap(lev),
m_mesh.boxArray(lev + 1), amrex::IntVect(2), 1, 0);
m_mesh.boxArray(lev + 1), m_mesh.refRatio(lev), 1, 0);
} else {
level_mask.define(
m_mesh.boxArray(lev), m_mesh.DistributionMap(lev), 1, 0,
Expand Down
2 changes: 1 addition & 1 deletion amr-wind/physics/EkmanSpiral.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,7 +146,7 @@ amrex::Real EkmanSpiral::compute_error(const Field& field)
if (lev < nlevels - 1) {
level_mask = makeFineMask(
m_mesh.boxArray(lev), m_mesh.DistributionMap(lev),
m_mesh.boxArray(lev + 1), amrex::IntVect(2), 1, 0);
m_mesh.boxArray(lev + 1), m_mesh.refRatio(lev), 1, 0);
} else {
level_mask.define(
m_mesh.boxArray(lev), m_mesh.DistributionMap(lev), 1, 0,
Expand Down
3 changes: 2 additions & 1 deletion amr-wind/physics/ScalarAdvection.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -257,7 +257,8 @@ ScalarAdvection::compute_error(const Shape& scalar_function)
level_mask = makeFineMask(
m_repo.mesh().boxArray(level),
m_repo.mesh().DistributionMap(level),
m_repo.mesh().boxArray(level + 1), amrex::IntVect(2), 1, 0);
m_repo.mesh().boxArray(level + 1),
m_repo.mesh().refRatio(level), 1, 0);
} else {
level_mask.define(
m_repo.mesh().boxArray(level),
Expand Down
2 changes: 1 addition & 1 deletion amr-wind/physics/mms/MMS.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,7 @@ MMS::compute_error(const int comp, const Field& field, amr_wind::mms::FuncDef f)
if (lev < nlevels - 1) {
level_mask = makeFineMask(
m_mesh.boxArray(lev), m_mesh.DistributionMap(lev),
m_mesh.boxArray(lev + 1), amrex::IntVect(2), 1, 0);
m_mesh.boxArray(lev + 1), m_mesh.refRatio(lev), 1, 0);
} else {
level_mask.define(
m_mesh.boxArray(lev), m_mesh.DistributionMap(lev), 1, 0,
Expand Down
4 changes: 2 additions & 2 deletions amr-wind/physics/multiphase/MultiPhase.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -214,7 +214,7 @@ amrex::Real MultiPhase::volume_fraction_sum()
if (lev < nlevels - 1) {
level_mask = makeFineMask(
mesh.boxArray(lev), mesh.DistributionMap(lev),
mesh.boxArray(lev + 1), amrex::IntVect(2), 1, 0);
mesh.boxArray(lev + 1), mesh.refRatio(lev), 1, 0);
} else {
level_mask.define(
mesh.boxArray(lev), mesh.DistributionMap(lev), 1, 0,
Expand Down Expand Up @@ -261,7 +261,7 @@ amrex::Real MultiPhase::momentum_sum(int n)
if (lev < nlevels - 1) {
level_mask = makeFineMask(
mesh.boxArray(lev), mesh.DistributionMap(lev),
mesh.boxArray(lev + 1), amrex::IntVect(2), 1, 0);
mesh.boxArray(lev + 1), mesh.refRatio(lev), 1, 0);
} else {
level_mask.define(
mesh.boxArray(lev), mesh.DistributionMap(lev), 1, 0,
Expand Down
3 changes: 2 additions & 1 deletion amr-wind/physics/multiphase/ZalesakDiskScalarVel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -225,7 +225,8 @@ amrex::Real ZalesakDiskScalarVel::compute_error(const Field& field)
if (lev < nlevels - 1) {
level_mask = makeFineMask(
m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
m_sim.mesh().boxArray(lev + 1), amrex::IntVect(2), 1, 0);
m_sim.mesh().boxArray(lev + 1), m_sim.mesh().refRatio(lev), 1,
0);
} else {
level_mask.define(
m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
Expand Down
4 changes: 2 additions & 2 deletions amr-wind/utilities/FieldPlaneAveragingFine.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -219,7 +219,7 @@ void FPlaneAveragingFine<FType>::compute_averages(const IndexSelector& idxOp)
if (lev < finestLevel) {
level_mask = makeFineMask(
mesh.boxArray(lev), mesh.DistributionMap(lev),
mesh.boxArray(lev + 1), amrex::IntVect(2), 1, 0);
mesh.boxArray(lev + 1), mesh.refRatio(lev), 1, 0);
} else {
level_mask.define(
mesh.boxArray(lev), mesh.DistributionMap(lev), 1, 0,
Expand Down Expand Up @@ -407,7 +407,7 @@ void VelPlaneAveragingFine::compute_hvelmag_averages(const IndexSelector& idxOp)
if (lev < finestLevel) {
level_mask = makeFineMask(
mesh.boxArray(lev), mesh.DistributionMap(lev),
mesh.boxArray(lev + 1), amrex::IntVect(2), 1, 0);
mesh.boxArray(lev + 1), mesh.refRatio(lev), 1, 0);
} else {
level_mask.define(
mesh.boxArray(lev), mesh.DistributionMap(lev), 1, 0,
Expand Down
40 changes: 32 additions & 8 deletions amr-wind/utilities/diagnostics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -167,7 +167,7 @@ amrex::Array<amrex::Real, 24> amr_wind::diagnostics::PrintMaxVelLocations(
if (lev < finest_level) {
level_mask = makeFineMask(
repo.mesh().boxArray(lev), repo.mesh().DistributionMap(lev),
repo.mesh().boxArray(lev + 1), amrex::IntVect(2), 1, 0);
repo.mesh().boxArray(lev + 1), repo.mesh().refRatio(lev), 1, 0);
} else {
level_mask.define(
repo.mesh().boxArray(lev), repo.mesh().DistributionMap(lev), 1,
Expand Down Expand Up @@ -221,7 +221,7 @@ amrex::Array<amrex::Real, 24> amr_wind::diagnostics::PrintMaxVelLocations(
if (lev < finest_level) {
level_mask = makeFineMask(
repo.mesh().boxArray(lev), repo.mesh().DistributionMap(lev),
repo.mesh().boxArray(lev + 1), amrex::IntVect(2), 1, 0);
repo.mesh().boxArray(lev + 1), repo.mesh().refRatio(lev), 1, 0);
} else {
level_mask.define(
repo.mesh().boxArray(lev), repo.mesh().DistributionMap(lev), 1,
Expand Down Expand Up @@ -346,11 +346,23 @@ amrex::Array<amrex::Real, 24> amr_wind::diagnostics::PrintMaxMACVelLocations(
amrex::Real uMAC_min{-1e8}, vMAC_min{-1e8}, wMAC_min{-1e8};
for (int lev = 0; lev <= finest_level; lev++) {
// Use level_mask to only count finest level present
// Do it with a ghost cell for the sake of checking faces
amrex::iMultiFab level_mask;
if (lev < finest_level) {
level_mask = makeFineMask(
repo.mesh().boxArray(lev), repo.mesh().DistributionMap(lev),
repo.mesh().boxArray(lev + 1), amrex::IntVect(2), 1, 1);
// MultiFab with ghost cell
level_mask.define(
repo.mesh().boxArray(lev), repo.mesh().DistributionMap(lev), 1,
1, amrex::MFInfo());
// Set to 0 (ghost cells will be 0)
level_mask.setVal(0);
// Populate non-ghost cells with fine mask
amrex::iMultiFab::Copy(
level_mask,
makeFineMask(
repo.mesh().boxArray(lev), repo.mesh().DistributionMap(lev),
repo.mesh().boxArray(lev + 1), repo.mesh().refRatio(lev), 1,
0),
0, 0, 1, amrex::IntVect(0));
} else {
level_mask.define(
repo.mesh().boxArray(lev), repo.mesh().DistributionMap(lev), 1,
Expand Down Expand Up @@ -412,11 +424,23 @@ amrex::Array<amrex::Real, 24> amr_wind::diagnostics::PrintMaxMACVelLocations(
problo[0], problo[1], problo[2]};
for (int lev = 0; lev <= finest_level; lev++) {
// Use level_mask to only count finest level present
// Do it with a ghost cell for the sake of checking faces
amrex::iMultiFab level_mask;
if (lev < finest_level) {
level_mask = makeFineMask(
repo.mesh().boxArray(lev), repo.mesh().DistributionMap(lev),
repo.mesh().boxArray(lev + 1), amrex::IntVect(2), 1, 1);
// MultiFab with ghost cell
level_mask.define(
repo.mesh().boxArray(lev), repo.mesh().DistributionMap(lev), 1,
1, amrex::MFInfo());
// Set to 0 (ghost cells will be 0)
level_mask.setVal(0);
// Populate non-ghost cells with fine mask
amrex::iMultiFab::Copy(
level_mask,
makeFineMask(
repo.mesh().boxArray(lev), repo.mesh().DistributionMap(lev),
repo.mesh().boxArray(lev + 1), repo.mesh().refRatio(lev), 1,
0),
0, 0, 1, amrex::IntVect(0));
} else {
level_mask.define(
repo.mesh().boxArray(lev), repo.mesh().DistributionMap(lev), 1,
Expand Down
3 changes: 2 additions & 1 deletion amr-wind/utilities/sampling/Enstrophy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,8 @@ amrex::Real Enstrophy::calculate_enstrophy()
if (lev < finest_level) {
level_mask = makeFineMask(
m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
m_sim.mesh().boxArray(lev + 1), amrex::IntVect(2), 1, 0);
m_sim.mesh().boxArray(lev + 1), m_sim.mesh().refRatio(lev), 1,
0);
} else {
level_mask.define(
m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
Expand Down
2 changes: 1 addition & 1 deletion amr-wind/utilities/sampling/FieldNorms.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ FieldNorms::L2_norm(amr_wind::Field& field, const int comp, const bool use_mask)
if (lev < finest_level) {
level_mask = makeFineMask(
mesh.boxArray(lev), mesh.DistributionMap(lev),
mesh.boxArray(lev + 1), amrex::IntVect(2), 1, 0);
mesh.boxArray(lev + 1), mesh.refRatio(lev), 1, 0);
} else {
level_mask.define(
mesh.boxArray(lev), mesh.DistributionMap(lev), 1, 0,
Expand Down
9 changes: 6 additions & 3 deletions amr-wind/utilities/sampling/FreeSurface.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,8 @@ void FreeSurface::initialize()
if (lev < finest_level) {
level_mask = makeFineMask(
m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
m_sim.mesh().boxArray(lev + 1), amrex::IntVect(2), 1, 0);
m_sim.mesh().boxArray(lev + 1), m_sim.mesh().refRatio(lev), 1,
0);
} else {
level_mask.define(
m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
Expand Down Expand Up @@ -221,7 +222,8 @@ void FreeSurface::initialize()
if (lev < finest_level) {
level_mask = makeFineMask(
m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
m_sim.mesh().boxArray(lev + 1), amrex::IntVect(2), 1, 0);
m_sim.mesh().boxArray(lev + 1), m_sim.mesh().refRatio(lev), 1,
0);
} else {
level_mask.define(
m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
Expand Down Expand Up @@ -548,7 +550,8 @@ void FreeSurface::post_regrid_actions()
if (lev < finest_level) {
level_mask = makeFineMask(
m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
m_sim.mesh().boxArray(lev + 1), amrex::IntVect(2), 1, 0);
m_sim.mesh().boxArray(lev + 1), m_sim.mesh().refRatio(lev), 1,
0);
} else {
level_mask.define(
m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
Expand Down
3 changes: 2 additions & 1 deletion amr-wind/utilities/sampling/KineticEnergy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,8 @@ amrex::Real KineticEnergy::calculate_kinetic_energy()
if (lev < finest_level) {
level_mask = makeFineMask(
m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
m_sim.mesh().boxArray(lev + 1), amrex::IntVect(2), 1, 0);
m_sim.mesh().boxArray(lev + 1), m_sim.mesh().refRatio(lev), 1,
0);
} else {
level_mask.define(
m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
Expand Down
6 changes: 4 additions & 2 deletions amr-wind/utilities/sampling/WaveEnergy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,8 @@ amrex::Real WaveEnergy::calculate_kinetic_energy()
if (lev < finest_level) {
level_mask = makeFineMask(
m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
m_sim.mesh().boxArray(lev + 1), amrex::IntVect(2), 1, 0);
m_sim.mesh().boxArray(lev + 1), m_sim.mesh().refRatio(lev), 1,
0);
} else {
level_mask.define(
m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
Expand Down Expand Up @@ -117,7 +118,8 @@ amrex::Real WaveEnergy::calculate_potential_energy()
if (lev < finest_level) {
level_mask = makeFineMask(
m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
m_sim.mesh().boxArray(lev + 1), amrex::IntVect(2), 1, 0);
m_sim.mesh().boxArray(lev + 1), m_sim.mesh().refRatio(lev), 1,
0);
} else {
level_mask.define(
m_sim.mesh().boxArray(lev), m_sim.mesh().DistributionMap(lev),
Expand Down
114 changes: 114 additions & 0 deletions unit_tests/utilities/test_diagnostics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,13 @@ void init_velocity(amr_wind::Field& velocity)
farr(i, j, k, 0) = 1.0 - std::pow(xc, 2.0);
farr(i, j, k, 1) = -1.0 + std::pow(zc, 2.0);
farr(i, j, k, 2) = 5.0 * std::cos(yc);

if (lev == 0 && nlevels > 1) {
// Set base level to large values to detect masking errors
farr(i, j, k, 0) = 1e5;
farr(i, j, k, 1) = -1e5;
farr(i, j, k, 2) = 1e5;
}
});
}
}
Expand Down Expand Up @@ -58,6 +65,13 @@ void init_mac_velocity(
uarr(i, j, k) = 1.0 - std::pow(x, 2.0);
varr(i, j, k) = -1.0 + std::pow(zc, 2.0);
warr(i, j, k) = -3.0 * std::cos(yc);

if (lev == 0 && nlevels > 1) {
// Set base level to large values to detect masking errors
uarr(i, j, k) = 1e5;
varr(i, j, k) = -1e5;
warr(i, j, k) = 1e5;
}
});
}
}
Expand Down Expand Up @@ -155,4 +169,104 @@ TEST_F(DiagnosticsTest, Max_MACvel)
EXPECT_NEAR(std::abs(fc_results[22]), 0.5 * 10.0 / 24.0, tol);
}

TEST_F(DiagnosticsTest, Max_Vel_MultiLevel)
{
populate_parameters();
{
amrex::ParmParse pp("amr");
pp.add("max_level", 1);
}
// Create the refinement input file
// Cover the whole domain for easier testing
std::stringstream ss;
ss << "1 // Number of levels" << std::endl;
ss << "1 // Number of boxes at this level" << std::endl;
ss << "-5 -5 -2 5 5 2" << std::endl;
create_mesh_instance<RefineMesh>();
std::unique_ptr<amr_wind::CartBoxRefinement> box_refine(
new amr_wind::CartBoxRefinement(sim()));
box_refine->read_inputs(mesh(), ss);
mesh<RefineMesh>()->refine_criteria_vec().push_back(std::move(box_refine));
initialize_mesh();

auto& repo = sim().repo();
auto& velocity = repo.declare_field("velocity", 3, 0);
init_velocity(velocity);

auto cc_results =
amr_wind::diagnostics::PrintMaxVelLocations(repo, "cell-centered");

// Check max's and min's, according to profiles
const amrex::Real tol = 1.0e-10;
// max(u)
EXPECT_NEAR(cc_results[0], 1.0 - std::pow(0.5 * 10.0 / 48.0, 2.0), tol);
// min(u)
EXPECT_NEAR(cc_results[4], 1.0 - std::pow(23.5 * 10.0 / 48.0, 2.0), tol);
// max(v)
EXPECT_NEAR(cc_results[8], -1.0 + std::pow(7.5 * 4.0 / 16.0, 2.0), tol);
// min(v)
EXPECT_NEAR(cc_results[12], -1.0 + std::pow(0.5 * 4.0 / 16.0, 2.0), tol);
// max(w)
EXPECT_NEAR(cc_results[16], 5.0 * std::cos(0.5 * 10.0 / 48.0), tol);

// Check locations (abs due to symmetry)
EXPECT_NEAR(std::abs(cc_results[1]), 0.5 * 10.0 / 48.0, tol);
EXPECT_NEAR(std::abs(cc_results[5]), 23.5 * 10.0 / 48.0, tol);
EXPECT_NEAR(std::abs(cc_results[11]), 7.5 * 4.0 / 16.0, tol);
EXPECT_NEAR(std::abs(cc_results[15]), 0.5 * 4.0 / 16.0, tol);
EXPECT_NEAR(std::abs(cc_results[18]), 0.5 * 10.0 / 48.0, tol);
}

TEST_F(DiagnosticsTest, Max_MACvel_MultiLevel)
{
populate_parameters();
{
amrex::ParmParse pp("amr");
pp.add("max_level", 1);
}
// Create the refinement input file
// Cover the whole domain for easier testing
std::stringstream ss;
ss << "1 // Number of levels" << std::endl;
ss << "1 // Number of boxes at this level" << std::endl;
ss << "-5 -5 -2 5 5 2" << std::endl;
create_mesh_instance<RefineMesh>();
std::unique_ptr<amr_wind::CartBoxRefinement> box_refine(
new amr_wind::CartBoxRefinement(sim()));
box_refine->read_inputs(mesh(), ss);
mesh<RefineMesh>()->refine_criteria_vec().push_back(std::move(box_refine));
initialize_mesh();

auto& repo = sim().repo();
repo.declare_face_normal_field({"u_mac", "v_mac", "w_mac"}, 1, 1, 1);
auto& umac = repo.get_field("u_mac");
auto& vmac = repo.get_field("v_mac");
auto& wmac = repo.get_field("w_mac");
auto& cc = repo.declare_field("cc", 1, 0);
init_mac_velocity(cc, umac, vmac, wmac);

auto fc_results =
amr_wind::diagnostics::PrintMaxMACVelLocations(repo, "face-centered");

// Check max's and min's, according to profiles
const amrex::Real tol = 1.0e-10;
// max(umac)
EXPECT_NEAR(fc_results[0], 1.0 - std::pow(0.0 * 10.0 / 48.0, 2.0), tol);
// min(umac)
EXPECT_NEAR(fc_results[4], 1.0 - std::pow(24 * 10.0 / 48.0, 2.0), tol);
// max(vmac)
EXPECT_NEAR(fc_results[8], -1.0 + std::pow(7.5 * 4.0 / 16.0, 2.0), tol);
// min(vmac)
EXPECT_NEAR(fc_results[12], -1.0 + std::pow(0.5 * 4.0 / 16.0, 2.0), tol);
// min(wmac)
EXPECT_NEAR(fc_results[20], -3.0 * std::cos(0.5 * 10.0 / 48.0), tol);

// Check locations
EXPECT_NEAR(fc_results[1], 0.0 * 10.0 / 24.0, tol);
EXPECT_NEAR(std::abs(fc_results[5]), 24 * 10.0 / 48.0, tol);
EXPECT_NEAR(std::abs(fc_results[11]), 7.5 * 4.0 / 16.0, tol);
EXPECT_NEAR(std::abs(fc_results[15]), 0.5 * 4.0 / 16.0, tol);
EXPECT_NEAR(std::abs(fc_results[22]), 0.5 * 10.0 / 48.0, tol);
}

} // namespace amr_wind_tests

0 comments on commit 3edc314

Please sign in to comment.