Skip to content

Commit 5a79177

Browse files
mbkuhnmbkuhn
authored
FreeSurface sampler: avoid issues when sampling within numerical beach (#1678)
* improve linear interpolation approach for free surface sampler * add ability to turn on linear interpolation within numerical beach * add test to show behavior when turned off * tweaking test and check for position * instead of picking slopes, check interpolated vof values for intersection * test that fails with current implementation * correct test value in unit test * implementation that passes test! * refactor a bit * add factor of safety * add to documentation
1 parent f3e85e6 commit 5a79177

File tree

4 files changed

+267
-23
lines changed

4 files changed

+267
-23
lines changed

‎amr-wind/utilities/sampling/FreeSurfaceSampler.H‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -110,6 +110,11 @@ private:
110110
//! Output coordinate
111111
amrex::Vector<amrex::Real> m_out;
112112

113+
//! Flag to use linear interpolation for surface finding in part of domain
114+
bool m_use_linear{false};
115+
//! Parameter for optional linear interpolation
116+
amrex::Real m_lx_linear{0.};
117+
113118
//! Max number of sample points found in a single cell
114119
int m_ncomp{1};
115120
//! Max number of sample points allowed in a single cell

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

Lines changed: 83 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,10 @@ void FreeSurfaceSampler::initialize(const std::string& key)
2727
pp.query("search_direction", m_coorddir);
2828
pp.query("num_instances", m_ninst);
2929
pp.query("max_sample_points_per_cell", m_ncmax);
30+
m_use_linear = pp.contains("linear_interp_extent_from_xhi");
31+
if (m_use_linear) {
32+
pp.get("linear_interp_extent_from_xhi", m_lx_linear);
33+
}
3034
AMREX_ALWAYS_ASSERT(static_cast<int>(m_start.size()) == AMREX_SPACEDIM);
3135
AMREX_ALWAYS_ASSERT(static_cast<int>(m_end.size()) == AMREX_SPACEDIM);
3236
AMREX_ALWAYS_ASSERT(static_cast<int>(m_npts_dir.size()) == 2);
@@ -436,6 +440,8 @@ bool FreeSurfaceSampler::update_sampling_locations()
436440
const int gc1 = m_gc1;
437441
const int ncomp = m_ncomp;
438442

443+
bool use_linear = m_use_linear;
444+
const amrex::Real lx_linear = m_lx_linear;
439445
bool has_overset = m_sim.has_overset();
440446
amr_wind::IntField* iblank_ptr{nullptr};
441447
if (has_overset) {
@@ -455,6 +461,7 @@ bool FreeSurfaceSampler::update_sampling_locations()
455461
geom.InvCellSizeArray();
456462
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> plo =
457463
geom.ProbLoArray();
464+
const amrex::Real xhi = geom.ProbHi(0);
458465
#ifdef AMREX_USE_OMP
459466
#pragma omp parallel if (false)
460467
#endif
@@ -472,6 +479,8 @@ bool FreeSurfaceSampler::update_sampling_locations()
472479
xm[0] = plo[0] + (i + 0.5) * dx[0];
473480
xm[1] = plo[1] + (j + 0.5) * dx[1];
474481
xm[2] = plo[2] + (k + 0.5) * dx[2];
482+
bool linear_on =
483+
use_linear && (xm[0] >= xhi - lx_linear);
475484
// Loop number of components
476485
for (int n = 0; n < ncomp; ++n) {
477486
// Get index of current component and cell
@@ -498,10 +507,19 @@ bool FreeSurfaceSampler::update_sampling_locations()
498507
// is not single-phase)
499508
bool calc_flag = false;
500509
bool calc_flag_diffuse = false;
501-
if ((ni % 2 == 0 &&
510+
const bool single_phase_below_interface =
511+
(ni % 2 == 0 &&
502512
vof_arr(i, j, k) >= 1.0 - 1e-12) ||
503-
(ni % 2 != 0 &&
504-
vof_arr(i, j, k) <= 1e-12)) {
513+
(ni % 2 != 0 && vof_arr(i, j, k) <= 1e-12);
514+
const bool multiphase =
515+
vof_arr(i, j, k) < (1.0 - 1e-12) &&
516+
vof_arr(i, j, k) > 1e-12;
517+
const bool use_linear_interp =
518+
(has_overset && ibl_arr(i, j, k) == -1) ||
519+
linear_on;
520+
const bool single_phase_above_interface = !(
521+
multiphase || single_phase_below_interface);
522+
if (single_phase_below_interface) {
505523
// put bdy at top
506524
alpha = 1.0;
507525
if (ni % 2 != 0) {
@@ -513,10 +531,8 @@ bool FreeSurfaceSampler::update_sampling_locations()
513531
calc_flag = true;
514532
}
515533
// Multiphase cell case
516-
if (vof_arr(i, j, k) < (1.0 - 1e-12) &&
517-
vof_arr(i, j, k) > 1e-12) {
518-
if (!has_overset ||
519-
ibl_arr(i, j, k) != -1) {
534+
if (multiphase) {
535+
if (!use_linear_interp) {
520536
// Planar reconstruction
521537
calc_flag = true;
522538
multiphase::fit_plane(
@@ -527,6 +543,31 @@ bool FreeSurfaceSampler::update_sampling_locations()
527543
calc_flag_diffuse = true;
528544
}
529545
}
546+
// Single-phase cell bordering other phase, does
547+
// not fit into other categories
548+
if (single_phase_above_interface &&
549+
use_linear_interp) {
550+
// With linear interp, interface can show up
551+
// in single-phase cell
552+
const amrex::IntVect iv{i, j, k};
553+
amrex::IntVect iv_p = iv;
554+
amrex::IntVect iv_m = iv;
555+
iv_p[dir] += 1;
556+
iv_m[dir] -= 1;
557+
// Check for 0.5 intersect
558+
const bool intersect_above =
559+
(vof_arr(iv_p) - 0.5) *
560+
(vof_arr(iv) - 0.5) <=
561+
0;
562+
const bool intersect_below =
563+
(vof_arr(iv) - 0.5) *
564+
(vof_arr(iv_m) - 0.5) <=
565+
0;
566+
calc_flag_diffuse =
567+
calc_flag_diffuse || intersect_above;
568+
calc_flag_diffuse =
569+
calc_flag_diffuse || intersect_below;
570+
}
530571

531572
// Initialize height measurement
532573
amrex::Real ht = plo[dir];
@@ -573,14 +614,13 @@ bool FreeSurfaceSampler::update_sampling_locations()
573614
// Distances from cell center to probe loc
574615
const amrex::Real dist_0 = loc0 - xm[gc0];
575616
const amrex::Real dist_1 = loc1 - xm[gc1];
576-
// Central slopes for when sign of distance
577-
// to interface is unknown
578-
amrex::Real slope_dir =
579-
dir == 0
580-
? 0.5 * (dv_xl + dv_xr)
581-
: (dir == 1
582-
? 0.5 * (dv_yl + dv_yr)
583-
: 0.5 * (dv_zl + dv_zr));
617+
// Slopes on either side
618+
amrex::Real slope_dir_r =
619+
dir == 0 ? dv_xr
620+
: (dir == 1 ? dv_yr : dv_zr);
621+
amrex::Real slope_dir_l =
622+
dir == 0 ? dv_xl
623+
: (dir == 1 ? dv_yl : dv_zl);
584624
// One-sided slopes for when sign of
585625
// distance to interface is known
586626
amrex::Real slope_0 =
@@ -590,14 +630,36 @@ bool FreeSurfaceSampler::update_sampling_locations()
590630
dist_1 > 0 ? (dir == 2 ? dv_yr : dv_zr)
591631
: (dir == 2 ? dv_yl : dv_zl);
592632
// Turn finite differences into true slopes
593-
slope_dir *= dxi[dir];
633+
slope_dir_r *= dxi[dir];
634+
slope_dir_l *= dxi[dir];
594635
slope_0 *= dxi[gc0];
595636
slope_1 *= dxi[gc1];
596-
// Trilinear interpolation
597-
ht = (0.5 + slope_dir * xm[dir] -
598-
(vof_arr(i, j, k) + slope_0 * dist_0 +
599-
slope_1 * dist_1)) /
600-
slope_dir;
637+
// Trilinear interpolation for vof
638+
const amrex::Real vof_c = vof_arr(i, j, k) +
639+
slope_0 * dist_0 +
640+
slope_1 * dist_1;
641+
// Extrapolate to cell edge (0.5) plus a
642+
// factor of safety (0.1) because
643+
// reconstruction is not identical in
644+
// neighboring cells
645+
const amrex::Real vof_r =
646+
vof_c + 0.6 * slope_dir_r * dx[dir];
647+
const amrex::Real vof_l =
648+
vof_c - 0.6 * slope_dir_l * dx[dir];
649+
// Check for intersect with 0.5
650+
if ((vof_c - 0.5) * (vof_r - 0.5) <= 0.) {
651+
ht = xm[dir] +
652+
(0.5 - vof_c) /
653+
(slope_dir_r + constants::EPS);
654+
} else if (
655+
(vof_c - 0.5) * (vof_l - 0.5) <= 0.) {
656+
ht = xm[dir] +
657+
(0.5 - vof_c) /
658+
(slope_dir_l + constants::EPS);
659+
} else {
660+
// Skip if no intersection
661+
ht = plo[dir];
662+
}
601663
}
602664

603665
if (calc_flag || calc_flag_diffuse) {

‎docs/sphinx/user/inputs_Sampling.rst‎

Lines changed: 14 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -199,15 +199,27 @@ direction, specified by ``search_direction``, and the other coordinates are
199199
determined by the ``plane_`` parameters. The default search direction parameter
200200
is 2, indicating the samplers will search for the interface along the z-direction.
201201
Due to this design, it is best to specify a plane that is normal to the intended
202-
search direction. Another optional parameter is ``num_instances``, which is available
202+
search direction.
203+
204+
Another optional parameter is ``num_instances``, which is available
203205
for cases where the interface location is multi-valued along the search direction,
204206
such as during wave breaking. This parameter defaults to 1, and the sampler will
205207
automatically select the highest position along the search direction when the interface
206208
location is multi-valued.
207209

210+
The free surface location is calculated with
211+
a geometric approach using the reconstruction of the interface in a computational
212+
cell. However, within the numerical beach of a wave simulation, the volume fraction distribution
213+
can become noisy, and the geometric approach produce noisy results. To avoid this,
214+
there is an option to use a linear interpolation approach instead within the beach,
215+
which helps to reduce the noise. This can be turned on using the argument
216+
``linear_interp_extent_from_xhi``, which specifies the extent from the upper domain boundary (in x)
217+
where linear interpolation should be used instead of the standard geometric approach. This
218+
input parameter should be set to the length of the numerical beach.
219+
208220
Example::
209221

210222
sampling.fs1.type = FreeSurfaceSampler
211223
sampling.fs1.plane_start = 4.0 -1.0 0.0
212224
sampling.fs1.plane_end = 0.0 1.0 0.0
213-
sampling.fs1.plane_num_points = 20 10
225+
sampling.fs1.plane_num_points = 20 10

‎unit_tests/utilities/test_free_surface.cpp‎

Lines changed: 165 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -147,6 +147,37 @@ void init_vof_diffuse(amr_wind::Field& vof_fld, amrex::Real water_level)
147147
amrex::Gpu::streamSynchronize();
148148
}
149149

150+
void init_vof_outliers(
151+
amr_wind::Field& vof_fld, amrex::Real water_level, bool distort_above)
152+
{
153+
const auto& mesh = vof_fld.repo().mesh();
154+
const int nlevels = vof_fld.repo().num_active_levels();
155+
156+
// Since VOF is cell centered
157+
amrex::Real offset = 0.5;
158+
159+
for (int lev = 0; lev < nlevels; ++lev) {
160+
const auto& dx = mesh.Geom(lev).CellSizeArray();
161+
const auto& problo = mesh.Geom(lev).ProbLoArray();
162+
const auto& farrs = vof_fld(lev).arrays();
163+
164+
amrex::ParallelFor(
165+
vof_fld(lev), vof_fld.num_grow(),
166+
[=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
167+
const amrex::Real z = problo[2] + (k + offset) * dx[2];
168+
if (std::abs(water_level - z) < 0.5 * dx[2]) {
169+
farrs[nbx](i, j, k) =
170+
(water_level - (z - 0.5 * dx[2])) / dx[2];
171+
} else if (z > water_level) {
172+
farrs[nbx](i, j, k) = distort_above ? 0.1 : 0.0;
173+
} else {
174+
farrs[nbx](i, j, k) = !distort_above ? 0.9 : 1.0;
175+
}
176+
});
177+
}
178+
amrex::Gpu::streamSynchronize();
179+
}
180+
150181
//! Custom mesh class to be able to refine like a simulation would
151182
// - combination of AmrTestMesh and incflo classes
152183
// - with ability to initialize the refiner and regrid
@@ -380,6 +411,17 @@ class FreeSurfaceTest : public MeshTest
380411
pp.addarr("plane_end", m_pt_coord);
381412
}
382413
void setup_grid_0d(int ninst) { setup_grid_0d(ninst, "freesurface"); }
414+
void setup_grid_0d(int ninst, amrex::Real xy_sample)
415+
{
416+
amrex::ParmParse pp("freesurface");
417+
pp.add("output_interval", 1);
418+
pp.add("num_instances", ninst);
419+
pp.addarr("plane_num_points", amrex::Vector<int>{1, 1});
420+
const amrex::Vector<amrex::Real> pt_coord{
421+
xy_sample, xy_sample, m_pt_coord[2]};
422+
pp.addarr("plane_start", pt_coord);
423+
pp.addarr("plane_end", pt_coord);
424+
}
383425
void setup_grid_2d(int ninst)
384426
{
385427
amrex::ParmParse pp("freesurface");
@@ -667,4 +709,127 @@ TEST_F(FreeSurfaceTest, point_diffuse)
667709
ASSERT_EQ(nout, 1);
668710
}
669711

712+
TEST_F(FreeSurfaceTest, point_outliers)
713+
{
714+
initialize_mesh();
715+
auto& repo = sim().repo();
716+
auto& vof = repo.declare_field("vof", 1, 2);
717+
setup_grid_0d(1);
718+
{
719+
amrex::ParmParse pp("freesurface");
720+
pp.add("linear_interp_extent_from_xhi", 66.1);
721+
}
722+
723+
amrex::Real water_lev = 61.5;
724+
// Sets up field with multiphase cells above interface
725+
init_vof_outliers(vof, water_lev, true);
726+
auto& m_sim = sim();
727+
FreeSurfaceImpl tool(m_sim);
728+
tool.initialize("freesurface");
729+
tool.update_sampling_locations();
730+
731+
// Linear interpolation will get different value
732+
amrex::Real local_vof = (water_lev - 60.) / 2.0;
733+
amrex::Real water_lev_linear =
734+
61. + (0.5 - local_vof) * (2. / (0.1 - local_vof));
735+
736+
// Check output value
737+
int nout = tool.check_output("~", water_lev_linear);
738+
ASSERT_EQ(nout, 1);
739+
// Check sampling locations
740+
int nsloc = tool.check_sloc("~");
741+
ASSERT_EQ(nsloc, 1);
742+
743+
// Now distort VOF field below interface
744+
init_vof_outliers(vof, water_lev, false);
745+
tool.update_sampling_locations();
746+
water_lev_linear = 61. + (0.5 - local_vof) * (2. / (0. - local_vof));
747+
nout = tool.check_output("~", water_lev_linear);
748+
ASSERT_EQ(nout, 1);
749+
nsloc = tool.check_sloc("~");
750+
ASSERT_EQ(nsloc, 1);
751+
752+
// Repeat with target cell having 0 < vof < 0.5
753+
water_lev = 60.5;
754+
init_vof_outliers(vof, water_lev, true);
755+
tool.update_sampling_locations();
756+
local_vof = (water_lev - 60.) / 2.0;
757+
water_lev_linear = 61. - (0.5 - local_vof) * (2. / (1.0 - local_vof));
758+
nout = tool.check_output("~", water_lev_linear);
759+
ASSERT_EQ(nout, 1);
760+
init_vof_outliers(vof, water_lev, false);
761+
tool.update_sampling_locations();
762+
water_lev_linear = 61. - (0.5 - local_vof) * (2. / (0.9 - local_vof));
763+
nout = tool.check_output("~", water_lev_linear);
764+
ASSERT_EQ(nout, 1);
765+
nsloc = tool.check_sloc("~");
766+
ASSERT_EQ(nsloc, 1);
767+
}
768+
769+
TEST_F(FreeSurfaceTest, point_outliers_off)
770+
{
771+
initialize_mesh();
772+
auto& repo = sim().repo();
773+
auto& vof = repo.declare_field("vof", 1, 2);
774+
setup_grid_0d(1);
775+
{
776+
amrex::ParmParse pp("freesurface");
777+
pp.add("linear_interp_extent_from_xhi", 64);
778+
}
779+
780+
amrex::Real water_lev = 61.5;
781+
// Sets up field with multiphase cells above interface
782+
init_vof_outliers(vof, water_lev, true);
783+
auto& m_sim = sim();
784+
FreeSurfaceImpl tool(m_sim);
785+
tool.initialize("freesurface");
786+
tool.update_sampling_locations();
787+
788+
// Because linear interpolation is off and there are multiphase cells all
789+
// the way above the interface, it will get a bad value, middle of top cell
790+
const amrex::Real water_lev_geom = 123.;
791+
792+
// Check output value
793+
int nout = tool.check_output("~", water_lev_geom);
794+
ASSERT_EQ(nout, 1);
795+
}
796+
797+
TEST_F(FreeSurfaceTest, point_diffuse_in_single_phase)
798+
{
799+
initialize_mesh();
800+
sim().activate_overset();
801+
auto& repo = sim().repo();
802+
auto& vof = repo.declare_field("vof", 1, 2);
803+
auto& iblank = repo.get_int_field("iblank_cell");
804+
iblank.setVal(-1);
805+
// Set up sampler to be on lateral edge of cell
806+
setup_grid_0d(1, 64. - 0.1);
807+
808+
// Set up water level to be just below bottom edge of cell
809+
const amrex::Real water_lev_diffuse = 60. - 0.1;
810+
const amrex::Real vof_slope = 0.1;
811+
// Set up vof with nonzero slope but 0 in current cell
812+
init_vof_slope(vof, water_lev_diffuse, vof_slope, 124.);
813+
auto& m_sim = sim();
814+
FreeSurfaceImpl tool(m_sim);
815+
tool.initialize("freesurface");
816+
tool.update_sampling_locations();
817+
818+
// Check output value
819+
const amrex::Real vof_cell = 0.;
820+
const amrex::Real vof_mz = (water_lev_diffuse - 58.) / 2.;
821+
const amrex::Real vof_pxy = (water_lev_diffuse + 4. * vof_slope - 60.) / 2.;
822+
const amrex::Real vof_c =
823+
vof_cell + 2. * (vof_pxy - vof_cell) / 4.0 * (2.0 - 0.1);
824+
const amrex::Real slope_z = (vof_cell - vof_mz) / 2.;
825+
const amrex::Real ht =
826+
61. + (0.5 - vof_c) / (slope_z + amr_wind::constants::EPS);
827+
const amrex::Real ht_est =
828+
water_lev_diffuse + 2.0 * vof_slope * (2.0 - 0.1);
829+
int nout = tool.check_output("~", ht);
830+
ASSERT_EQ(nout, 1);
831+
// The linear interpolation answer should be close to the naive calc
832+
EXPECT_NEAR(ht, ht_est, 5e-2);
833+
}
834+
670835
} // namespace amr_wind_tests

0 commit comments

Comments
 (0)