Skip to content

Commit 8c1a6a3

Browse files
mbkuhnmbkuhn
authored
Correct field norms calculation for variables that are not cell-centered (#1703)
* set up test that fails with current field norms not at cell center * give realistic tolerances * get new test to pass: - use a real mask instead of int - keep box cell centered (using mask) - extend box based on face or node index - keep mask indices bounded * make test harder to pass * test with different refinement, increase field value * remove unused variable
1 parent 66a70ff commit 8c1a6a3

File tree

2 files changed

+119
-14
lines changed

2 files changed

+119
-14
lines changed

‎amr-wind/utilities/sampling/FieldNorms.cpp‎

Lines changed: 38 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -53,37 +53,65 @@ amrex::Real FieldNorms::l2_norm(
5353
geom[lev].CellSize()[1] *
5454
geom[lev].CellSize()[2];
5555

56-
amrex::iMultiFab level_mask;
56+
auto index_type = field(lev).boxArray().ixType();
57+
int it_sum = 0;
58+
int node_dir = 0;
59+
for (int ix = 0; ix < AMREX_SPACEDIM; ++ix) {
60+
it_sum += index_type[ix];
61+
node_dir = index_type[ix] == 1 ? ix : node_dir;
62+
}
63+
64+
amrex::MultiFab level_mask;
5765
if (use_mask) {
5866
if (lev < finest_level) {
5967
level_mask = makeFineMask(
6068
mesh.boxArray(lev), mesh.DistributionMap(lev),
61-
mesh.boxArray(lev + 1), mesh.refRatio(lev), 1, 0);
69+
mesh.boxArray(lev + 1), mesh.refRatio(lev), 1., 0.);
6270
} else {
6371
level_mask.define(
6472
mesh.boxArray(lev), mesh.DistributionMap(lev), 1, 0,
6573
amrex::MFInfo());
66-
level_mask.setVal(1);
74+
level_mask.setVal(1.);
6775
}
6876
} else {
6977
// Always on
7078
level_mask.define(
7179
mesh.boxArray(lev), mesh.DistributionMap(lev), 1, 0,
7280
amrex::MFInfo());
73-
level_mask.setVal(1);
81+
level_mask.setVal(1.);
7482
}
7583

7684
nrm += amrex::ReduceSum(
77-
field(lev), level_mask, nghost,
85+
level_mask, field(lev), nghost,
7886
[=] AMREX_GPU_HOST_DEVICE(
7987
amrex::Box const& bx,
80-
amrex::Array4<amrex::Real const> const& field_arr,
81-
amrex::Array4<int const> const& mask_arr) -> amrex::Real {
88+
amrex::Array4<amrex::Real const> const& mask_arr,
89+
amrex::Array4<amrex::Real const> const& field_arr)
90+
-> amrex::Real {
8291
amrex::Real nrm_fab = 0.0;
8392

84-
amrex::Loop(bx, [=, &nrm_fab](int i, int j, int k) noexcept {
85-
nrm_fab += cell_vol * field_arr(i, j, k, comp) *
86-
field_arr(i, j, k, comp) * mask_arr(i, j, k);
93+
const auto& fbx =
94+
it_sum == 1
95+
? amrex::surroundingNodes(bx, node_dir)
96+
: (it_sum > 1 ? amrex::surroundingNodes(bx) : bx);
97+
98+
amrex::Loop(fbx, [=, &nrm_fab](int i, int j, int k) noexcept {
99+
const amrex::IntVect iv{i, j, k};
100+
amrex::IntVect iv_cc = iv;
101+
// adjust volume for different data locations
102+
amrex::Real data_vol = cell_vol;
103+
for (int nn = 0; nn < AMREX_SPACEDIM; ++nn) {
104+
const bool at_node_dir =
105+
index_type[nn] == 1 && (iv[nn] == fbx.bigEnd(nn) ||
106+
iv[nn] == fbx.smallEnd(nn));
107+
data_vol *= at_node_dir ? 0.5 : 1.0;
108+
// limit mask array indices to cell-centered
109+
iv_cc[nn] = amrex::min(
110+
bx.bigEnd(nn),
111+
amrex::max(bx.smallEnd(nn), iv_cc[nn]));
112+
}
113+
nrm_fab += data_vol * field_arr(iv, comp) *
114+
field_arr(iv, comp) * mask_arr(iv_cc);
87115
});
88116
return nrm_fab;
89117
});

‎unit_tests/utilities/test_field_norms.cpp‎

Lines changed: 81 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -101,6 +101,7 @@ class FieldNormsImpl : public amr_wind::field_norms::FieldNorms
101101
{}
102102
void check_output(
103103
amrex::Real check_val0, amrex::Real check_val1, amrex::Real check_val2);
104+
void check_output(amrex::Real check_val);
104105

105106
protected:
106107
// No file output during test
@@ -117,10 +118,24 @@ void FieldNormsImpl::check_output(
117118
EXPECT_EQ(var_names()[1], (std::string) "velocityy");
118119
EXPECT_EQ(var_names()[2], (std::string) "velocityz");
119120
// Loop through norm values and check them
120-
const amrex::Real tiny = std::numeric_limits<amrex::Real>::epsilon();
121-
EXPECT_NEAR(field_norms()[0], check_val0, 500 * tiny);
122-
EXPECT_NEAR(field_norms()[1], check_val1, 500 * tiny);
123-
EXPECT_NEAR(field_norms()[2], check_val2, 500 * tiny);
121+
const amrex::Real tol = amr_wind::constants::TIGHT_TOL;
122+
EXPECT_NEAR(field_norms()[0], check_val0, tol);
123+
EXPECT_NEAR(field_norms()[1], check_val1, tol);
124+
EXPECT_NEAR(field_norms()[2], check_val2, tol);
125+
}
126+
127+
void FieldNormsImpl::check_output(amrex::Real check_val)
128+
{
129+
// Get variable names and check
130+
EXPECT_EQ(var_names()[0], (std::string) "pressure");
131+
EXPECT_EQ(var_names()[1], (std::string) "u_mac");
132+
EXPECT_EQ(var_names()[2], (std::string) "v_mac");
133+
EXPECT_EQ(var_names()[3], (std::string) "w_mac");
134+
// Loop through norm values and check them
135+
const amrex::Real tol = amr_wind::constants::TIGHT_TOL;
136+
for (int n = 0; n < 4; ++n) {
137+
EXPECT_NEAR(field_norms()[n], check_val, tol * check_val * 0.1);
138+
}
124139
}
125140

126141
} // namespace
@@ -473,4 +488,66 @@ TEST_F(FieldNormsTest, levelmask_on_int_refinement)
473488
tool.check_output(unorm, vnorm, wnorm);
474489
}
475490

491+
TEST_F(FieldNormsTest, levelmask_not_cc)
492+
{
493+
bool levelmask = true;
494+
// Set up parameters for domain
495+
populate_parameters();
496+
// Set up parameters for refinement
497+
setup_fieldrefinement();
498+
// Set up parameters for sampler
499+
setup_fnorm(levelmask);
500+
// Create mesh and initialize
501+
reset_prob_domain();
502+
auto rmesh = FNRefinemesh();
503+
rmesh.initialize_mesh(0.0);
504+
505+
// Repo and fields
506+
auto& repo = rmesh.field_repo();
507+
auto& pressure = repo.declare_nd_field("pressure", 1, 1);
508+
auto& u_mac = repo.declare_xf_field("u_mac", 1, 1);
509+
auto& v_mac = repo.declare_yf_field("v_mac", 1, 1);
510+
auto& w_mac = repo.declare_zf_field("w_mac", 1, 1);
511+
auto& flag = repo.declare_field(m_fname, 1, 1);
512+
513+
// Set up scalar for determining refinement - all fine level
514+
flag.setVal(2.0 * m_fref_val);
515+
516+
// Initialize mesh refiner and remesh
517+
rmesh.init_refiner();
518+
rmesh.remesh();
519+
520+
// Initialize pressure distribution and access sim
521+
const amrex::Real fval = 10000.3;
522+
pressure.setVal(fval);
523+
u_mac.setVal(fval);
524+
v_mac.setVal(fval);
525+
w_mac.setVal(fval);
526+
auto& rsim = rmesh.sim();
527+
528+
// Initialize IOManager because FieldNorms relies on it
529+
auto& io_mgr = rsim.io_manager();
530+
// Set up output (plot) variables
531+
io_mgr.register_output_var("pressure");
532+
io_mgr.register_output_var("u_mac");
533+
io_mgr.register_output_var("v_mac");
534+
io_mgr.register_output_var("w_mac");
535+
io_mgr.initialize_io();
536+
537+
// Initialize sampler and check result on initial mesh
538+
FieldNormsImpl tool(rsim, "fieldnorm");
539+
tool.initialize();
540+
tool.output_actions();
541+
// Only highest level will be counted
542+
tool.check_output(fval);
543+
544+
// Turn off refinements
545+
flag.setVal(0.0);
546+
rmesh.remesh();
547+
548+
// Check again
549+
tool.output_actions();
550+
tool.check_output(fval);
551+
}
552+
476553
} // namespace amr_wind_tests

0 commit comments

Comments
 (0)