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 e99fe104e5..d09776b124 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( @@ -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 = @@ -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) { @@ -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); diff --git a/unit_tests/utilities/test_free_surface.cpp b/unit_tests/utilities/test_free_surface.cpp index 5d681b8509..2dd4e907e7 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,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