Skip to content

FreeSurface sampler: avoid issues when sampling within numerical beach #1678

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

Draft
wants to merge 6 commits into
base: main
Choose a base branch
from
Draft

FreeSurface sampler: avoid issues when sampling within numerical beach #1678

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
f90d440
improve linear interpolation approach for free surface sampler
mbkuhn Jul 3, 2025
5d6c053
add ability to turn on linear interpolation within numerical beach
mbkuhn Jul 3, 2025
87d8003
add test to show behavior when turned off
mbkuhn Jul 3, 2025
d821c9f
tweaking test and check for position
mbkuhn Jul 3, 2025
385e4ec
Merge branch 'main' into free_surface_skip_outliers
mbkuhn Jul 21, 2025
949625e
instead of picking slopes, check interpolated vof values for intersec…
mbkuhn Jul 21, 2025
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
5 changes: 5 additions & 0 deletions amr-wind/utilities/sampling/FreeSurfaceSampler.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,6 +110,11 @@ private:
//! Output coordinate
amrex::Vector<amrex::Real> m_out;

//! Flag to use linear interpolation for surface finding in part of domain
bool m_use_linear{false};
//! Parameter for optional linear interpolation
amrex::Real m_lx_linear{0.};

//! Max number of sample points found in a single cell
int m_ncomp{1};
//! Max number of sample points allowed in a single cell
Expand Down
67 changes: 50 additions & 17 deletions amr-wind/utilities/sampling/FreeSurfaceSampler.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,10 @@ void FreeSurfaceSampler::initialize(const std::string& key)
pp.query("search_direction", m_coorddir);
pp.query("num_instances", m_ninst);
pp.query("max_sample_points_per_cell", m_ncmax);
m_use_linear = pp.contains("linear_interp_extent_from_xhi");
if (m_use_linear) {
pp.get("linear_interp_extent_from_xhi", m_lx_linear);
}
AMREX_ALWAYS_ASSERT(static_cast<int>(m_start.size()) == AMREX_SPACEDIM);
AMREX_ALWAYS_ASSERT(static_cast<int>(m_end.size()) == AMREX_SPACEDIM);
AMREX_ALWAYS_ASSERT(static_cast<int>(m_npts_dir.size()) == 2);
Expand Down Expand Up @@ -436,6 +440,8 @@ bool FreeSurfaceSampler::update_sampling_locations()
const int gc1 = m_gc1;
const int ncomp = m_ncomp;

bool use_linear = m_use_linear;
const amrex::Real lx_linear = m_lx_linear;
bool has_overset = m_sim.has_overset();
amr_wind::IntField* iblank_ptr{nullptr};
if (has_overset) {
Expand All @@ -455,6 +461,7 @@ bool FreeSurfaceSampler::update_sampling_locations()
geom.InvCellSizeArray();
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> plo =
geom.ProbLoArray();
const amrex::Real xhi = geom.ProbHi(0);
#ifdef AMREX_USE_OMP
#pragma omp parallel if (false)
#endif
Expand All @@ -472,6 +479,8 @@ bool FreeSurfaceSampler::update_sampling_locations()
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];
bool linear_on =
use_linear && (xm[0] >= xhi - lx_linear);
// Loop number of components
for (int n = 0; n < ncomp; ++n) {
// Get index of current component and cell
Expand Down Expand Up @@ -515,8 +524,9 @@ bool FreeSurfaceSampler::update_sampling_locations()
// Multiphase cell case
if (vof_arr(i, j, k) < (1.0 - 1e-12) &&
vof_arr(i, j, k) > 1e-12) {
if (!has_overset ||
ibl_arr(i, j, k) != -1) {
if ((!has_overset ||
ibl_arr(i, j, k) != -1) &&
!linear_on) {
// Planar reconstruction
calc_flag = true;
multiphase::fit_plane(
Expand Down Expand Up @@ -573,14 +583,13 @@ bool FreeSurfaceSampler::update_sampling_locations()
// Distances from cell center to probe loc
const amrex::Real dist_0 = loc0 - xm[gc0];
const amrex::Real dist_1 = loc1 - xm[gc1];
// Central slopes for when sign of distance
// to interface is unknown
amrex::Real slope_dir =
dir == 0
? 0.5 * (dv_xl + dv_xr)
: (dir == 1
? 0.5 * (dv_yl + dv_yr)
: 0.5 * (dv_zl + dv_zr));
// Slopes on either side
amrex::Real slope_dir_r =
dir == 0 ? dv_xr
: (dir == 1 ? dv_yr : dv_zr);
amrex::Real slope_dir_l =
dir == 0 ? dv_xl
: (dir == 1 ? dv_yl : dv_zl);
// One-sided slopes for when sign of
// distance to interface is known
amrex::Real slope_0 =
Expand All @@ -590,14 +599,36 @@ bool FreeSurfaceSampler::update_sampling_locations()
dist_1 > 0 ? (dir == 2 ? dv_yr : dv_zr)
: (dir == 2 ? dv_yl : dv_zl);
// Turn finite differences into true slopes
slope_dir *= dxi[dir];
slope_dir_r *= dxi[dir];
slope_dir_l *= dxi[dir];
slope_0 *= dxi[gc0];
slope_1 *= dxi[gc1];
// Trilinear interpolation
ht = (0.5 + slope_dir * xm[dir] -
(vof_arr(i, j, k) + slope_0 * dist_0 +
slope_1 * dist_1)) /
slope_dir;
// Trilinear interpolation for vof
const amrex::Real vof_c = vof_arr(i, j, k) +
slope_0 * dist_0 +
slope_1 * dist_1;
const amrex::Real vof_r =
vof_c + 0.5 * slope_dir_r * dx[dir];
const amrex::Real vof_l =
vof_c - 0.5 * slope_dir_l * dx[dir];
// Check for intersect with 0.5
if ((vof_c - 0.5) * (vof_r - 0.5) <= 0.) {
ht = xm[dir] +
(0.5 - (vof_arr(i, j, k) +
slope_0 * dist_0 +
slope_1 * dist_1)) /
(slope_dir_r + constants::EPS);
} else if (
(vof_c - 0.5) * (vof_l - 0.5) <= 0.) {
ht = xm[dir] +
(0.5 - (vof_arr(i, j, k) +
slope_0 * dist_0 +
slope_1 * dist_1)) /
(slope_dir_l + constants::EPS);
} else {
// Skip if no intersection
ht = plo[dir];
}
}

if (calc_flag || calc_flag_diffuse) {
Expand All @@ -608,9 +639,11 @@ bool FreeSurfaceSampler::update_sampling_locations()
}
// If interface is above upper
// bound, limit it
// Ignore it in the diffuse case
if (ht > xm[dir] +
0.5 * dx[dir] * (1.0 + 1e-8)) {
ht = xm[dir] + 0.5 * dx[dir];
ht = calc_flag ? xm[dir] + 0.5 * dx[dir]
: plo[dir];
}
// Save interface location by atomic max
amrex::Gpu::Atomic::Max(&dout_ptr[idx], ht);
Expand Down
116 changes: 116 additions & 0 deletions unit_tests/utilities/test_free_surface.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -147,6 +147,37 @@ void init_vof_diffuse(amr_wind::Field& vof_fld, amrex::Real water_level)
amrex::Gpu::streamSynchronize();
}

void init_vof_outliers(
amr_wind::Field& vof_fld, amrex::Real water_level, bool distort_above)
{
const auto& mesh = vof_fld.repo().mesh();
const int nlevels = vof_fld.repo().num_active_levels();

// Since VOF is cell centered
amrex::Real offset = 0.5;

for (int lev = 0; lev < nlevels; ++lev) {
const auto& dx = mesh.Geom(lev).CellSizeArray();
const auto& problo = mesh.Geom(lev).ProbLoArray();
const auto& farrs = vof_fld(lev).arrays();

amrex::ParallelFor(
vof_fld(lev), vof_fld.num_grow(),
[=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
const amrex::Real z = problo[2] + (k + offset) * dx[2];
if (std::abs(water_level - z) < 0.5 * dx[2]) {
farrs[nbx](i, j, k) =
(water_level - (z - 0.5 * dx[2])) / dx[2];
} else if (z > water_level) {
farrs[nbx](i, j, k) = distort_above ? 0.1 : 0.0;
} else {
farrs[nbx](i, j, k) = !distort_above ? 0.9 : 1.0;
}
});
}
amrex::Gpu::streamSynchronize();
}

//! Custom mesh class to be able to refine like a simulation would
// - combination of AmrTestMesh and incflo classes
// - with ability to initialize the refiner and regrid
Expand Down Expand Up @@ -667,4 +698,89 @@ TEST_F(FreeSurfaceTest, point_diffuse)
ASSERT_EQ(nout, 1);
}

TEST_F(FreeSurfaceTest, point_outliers)
{
initialize_mesh();
auto& repo = sim().repo();
auto& vof = repo.declare_field("vof", 1, 2);
setup_grid_0d(1);
{
amrex::ParmParse pp("freesurface");
pp.add("linear_interp_extent_from_xhi", 66.1);
}

amrex::Real water_lev = 61.5;
// Sets up field with multiphase cells above interface
init_vof_outliers(vof, water_lev, true);
auto& m_sim = sim();
FreeSurfaceImpl tool(m_sim);
tool.initialize("freesurface");
tool.update_sampling_locations();

// Linear interpolation will get different value
amrex::Real local_vof = (water_lev - 60.) / 2.0;
amrex::Real water_lev_linear =
61. + (0.5 - local_vof) * (2. / (0.1 - local_vof));

// Check output value
int nout = tool.check_output("~", water_lev_linear);
ASSERT_EQ(nout, 1);
// Check sampling locations
int nsloc = tool.check_sloc("~");
ASSERT_EQ(nsloc, 1);

// Now distort VOF field below interface
init_vof_outliers(vof, water_lev, false);
tool.update_sampling_locations();
water_lev_linear = 61. + (0.5 - local_vof) * (2. / (0. - local_vof));
nout = tool.check_output("~", water_lev_linear);
ASSERT_EQ(nout, 1);
nsloc = tool.check_sloc("~");
ASSERT_EQ(nsloc, 1);

// Repeat with target cell having 0 < vof < 0.5
water_lev = 60.5;
init_vof_outliers(vof, water_lev, true);
tool.update_sampling_locations();
local_vof = (water_lev - 60.) / 2.0;
water_lev_linear = 61. - (0.5 - local_vof) * (2. / (1.0 - local_vof));
nout = tool.check_output("~", water_lev_linear);
ASSERT_EQ(nout, 1);
init_vof_outliers(vof, water_lev, false);
tool.update_sampling_locations();
water_lev_linear = 61. - (0.5 - local_vof) * (2. / (0.9 - local_vof));
nout = tool.check_output("~", water_lev_linear);
ASSERT_EQ(nout, 1);
nsloc = tool.check_sloc("~");
ASSERT_EQ(nsloc, 1);
}

TEST_F(FreeSurfaceTest, point_outliers_off)
{
initialize_mesh();
auto& repo = sim().repo();
auto& vof = repo.declare_field("vof", 1, 2);
setup_grid_0d(1);
{
amrex::ParmParse pp("freesurface");
pp.add("linear_interp_extent_from_xhi", 64);
}

amrex::Real water_lev = 61.5;
// Sets up field with multiphase cells above interface
init_vof_outliers(vof, water_lev, true);
auto& m_sim = sim();
FreeSurfaceImpl tool(m_sim);
tool.initialize("freesurface");
tool.update_sampling_locations();

// Because linear interpolation is off and there are multiphase cells all
// the way above the interface, it will get a bad value, middle of top cell
const amrex::Real water_lev_geom = 123.;

// Check output value
int nout = tool.check_output("~", water_lev_geom);
ASSERT_EQ(nout, 1);
}

} // namespace amr_wind_tests
Loading