From f90d44096b4b0853caccff1a25a230278b2cf5ff Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Thu, 3 Jul 2025 16:21:37 -0600 Subject: [PATCH 1/5] improve linear interpolation approach for free surface sampler --- .../utilities/sampling/FreeSurfaceSampler.cpp | 40 ++++++++++++++----- 1 file changed, 30 insertions(+), 10 deletions(-) diff --git a/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp b/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp index e99fe104e5..c6a312a861 100644 --- a/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp +++ b/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp @@ -573,14 +573,17 @@ 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 + // Slopes on either side + const amrex::Real slope_dir_r = + dir == 0 ? dv_xr + : (dir == 1 ? dv_yr : dv_zr); + const amrex::Real slope_dir_l = + dir == 0 ? dv_xl + : (dir == 1 ? dv_yl : dv_zl); + // Central slope 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)); + 0.5 * (slope_dir_r + slope_dir_l); // One-sided slopes for when sign of // distance to interface is known amrex::Real slope_0 = @@ -593,11 +596,26 @@ bool FreeSurfaceSampler::update_sampling_locations() slope_dir *= dxi[dir]; slope_0 *= dxi[gc0]; slope_1 *= dxi[gc1]; - // Trilinear interpolation - ht = (0.5 + slope_dir * xm[dir] - + // Trilinear interpolation: central slope + ht = xm[dir] + + (0.5 - (vof_arr(i, j, k) + slope_0 * dist_0 + slope_1 * dist_1)) / - slope_dir; + (slope_dir + constants::EPS); + // Choose slope based on signed distance + slope_dir = ht - xm[dir] > 0 ? slope_dir_r + : slope_dir_l; + slope_dir *= dxi[dir]; + // Trilinear interpolation: one-sided slope + ht = xm[dir] + + (0.5 - + (vof_arr(i, j, k) + slope_0 * dist_0 + + slope_1 * dist_1)) / + (slope_dir + constants::EPS); + // This way, neighboring cells will agree on + // the reconstruction, avoiding the + // possibility of neither one containing the + // interface } if (calc_flag || calc_flag_diffuse) { @@ -608,9 +626,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); From 5d6c053bb53950e004028dcaa465d65ba526a449 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Thu, 3 Jul 2025 16:22:24 -0600 Subject: [PATCH 2/5] add ability to turn on linear interpolation within numerical beach --- .../utilities/sampling/FreeSurfaceSampler.H | 5 ++ .../utilities/sampling/FreeSurfaceSampler.cpp | 14 ++- unit_tests/utilities/test_free_surface.cpp | 87 +++++++++++++++++++ 3 files changed, 104 insertions(+), 2 deletions(-) diff --git a/amr-wind/utilities/sampling/FreeSurfaceSampler.H b/amr-wind/utilities/sampling/FreeSurfaceSampler.H index 699ddfaf4f..0cde60a1c9 100644 --- a/amr-wind/utilities/sampling/FreeSurfaceSampler.H +++ b/amr-wind/utilities/sampling/FreeSurfaceSampler.H @@ -110,6 +110,11 @@ private: //! Output coordinate amrex::Vector 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 diff --git a/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp b/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp index c6a312a861..e8df264405 100644 --- a/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp +++ b/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp @@ -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(m_start.size()) == AMREX_SPACEDIM); AMREX_ALWAYS_ASSERT(static_cast(m_end.size()) == AMREX_SPACEDIM); AMREX_ALWAYS_ASSERT(static_cast(m_npts_dir.size()) == 2); @@ -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) { @@ -455,6 +461,7 @@ bool FreeSurfaceSampler::update_sampling_locations() geom.InvCellSizeArray(); const amrex::GpuArray plo = geom.ProbLoArray(); + const amrex::Real xhi = geom.ProbHi(0); #ifdef AMREX_USE_OMP #pragma omp parallel if (false) #endif @@ -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 @@ -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( diff --git a/unit_tests/utilities/test_free_surface.cpp b/unit_tests/utilities/test_free_surface.cpp index d0ce299dc3..f1925e23b8 100644 --- a/unit_tests/utilities/test_free_surface.cpp +++ b/unit_tests/utilities/test_free_surface.cpp @@ -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 @@ -667,4 +698,60 @@ 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", 100); + } + + 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); +} + } // namespace amr_wind_tests From 87d8003d59aed62d318d3b8e1faab93f7e4b530c Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Thu, 3 Jul 2025 16:38:02 -0600 Subject: [PATCH 3/5] add test to show behavior when turned off --- unit_tests/utilities/test_free_surface.cpp | 28 ++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/unit_tests/utilities/test_free_surface.cpp b/unit_tests/utilities/test_free_surface.cpp index f1925e23b8..ab2b3621e4 100644 --- a/unit_tests/utilities/test_free_surface.cpp +++ b/unit_tests/utilities/test_free_surface.cpp @@ -754,4 +754,32 @@ TEST_F(FreeSurfaceTest, point_outliers) 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 From d821c9fb0229690a72217bb5063783e4423440e9 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Thu, 3 Jul 2025 16:47:51 -0600 Subject: [PATCH 4/5] tweaking test and check for position --- amr-wind/utilities/sampling/FreeSurfaceSampler.cpp | 2 +- unit_tests/utilities/test_free_surface.cpp | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp b/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp index e8df264405..26fc86937b 100644 --- a/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp +++ b/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp @@ -480,7 +480,7 @@ bool FreeSurfaceSampler::update_sampling_locations() 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); + 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 diff --git a/unit_tests/utilities/test_free_surface.cpp b/unit_tests/utilities/test_free_surface.cpp index ab2b3621e4..99422fc5c7 100644 --- a/unit_tests/utilities/test_free_surface.cpp +++ b/unit_tests/utilities/test_free_surface.cpp @@ -706,7 +706,7 @@ TEST_F(FreeSurfaceTest, point_outliers) setup_grid_0d(1); { amrex::ParmParse pp("freesurface"); - pp.add("linear_interp_extent_from_xhi", 100); + pp.add("linear_interp_extent_from_xhi", 66.1); } amrex::Real water_lev = 61.5; @@ -719,7 +719,8 @@ TEST_F(FreeSurfaceTest, point_outliers) // 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)); + 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); From 949625ed63c4bbd122208ca13501a5e8fe4453d0 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Mon, 21 Jul 2025 13:18:12 -0600 Subject: [PATCH 5/5] instead of picking slopes, check interpolated vof values for intersection --- .../utilities/sampling/FreeSurfaceSampler.cpp | 57 ++++++++++--------- 1 file changed, 30 insertions(+), 27 deletions(-) diff --git a/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp b/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp index 26fc86937b..d09776b124 100644 --- a/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp +++ b/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp @@ -584,16 +584,12 @@ bool FreeSurfaceSampler::update_sampling_locations() const amrex::Real dist_0 = loc0 - xm[gc0]; const amrex::Real dist_1 = loc1 - xm[gc1]; // Slopes on either side - const amrex::Real slope_dir_r = + amrex::Real slope_dir_r = dir == 0 ? dv_xr : (dir == 1 ? dv_yr : dv_zr); - const amrex::Real slope_dir_l = + amrex::Real slope_dir_l = dir == 0 ? dv_xl : (dir == 1 ? dv_yl : dv_zl); - // Central slope for when sign of distance - // to interface is unknown - amrex::Real slope_dir = - 0.5 * (slope_dir_r + slope_dir_l); // One-sided slopes for when sign of // distance to interface is known amrex::Real slope_0 = @@ -603,29 +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: central slope - ht = xm[dir] + - (0.5 - - (vof_arr(i, j, k) + slope_0 * dist_0 + - slope_1 * dist_1)) / - (slope_dir + constants::EPS); - // Choose slope based on signed distance - slope_dir = ht - xm[dir] > 0 ? slope_dir_r - : slope_dir_l; - slope_dir *= dxi[dir]; - // Trilinear interpolation: one-sided slope - ht = xm[dir] + - (0.5 - - (vof_arr(i, j, k) + slope_0 * dist_0 + - slope_1 * dist_1)) / - (slope_dir + constants::EPS); - // This way, neighboring cells will agree on - // the reconstruction, avoiding the - // possibility of neither one containing the - // interface + // 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) {