Skip to content

Commit

Permalink
Underwater terrain-aware wave forcing (#1475)
Browse files Browse the repository at this point in the history 
* turn off apply_relaxation_zones where terrain blanking (or terrain drag forcing) happens

* modify ow boundary for underwater terrain

* add regression test
  • Loading branch information
@mbkuhn
mbkuhn authored Feb 3, 2025
1 parent 4f11d0b commit 38d1b9f
Show file tree
Hide file tree
Showing 6 changed files with 2,291 additions and 8 deletions.
8 changes: 8 additions & 0 deletions amr-wind/ocean_waves/boundary_ops/OceanWavesBoundary.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
#define OCEANWAVESBOUNDARY_H

#include "amr-wind/core/Field.H"
#include "amr-wind/core/IntField.H"
#include "amr-wind/CFDSim.H"
#include "AMReX_Gpu.H"
#include <AMReX_BndryRegister.H>
Expand Down Expand Up @@ -61,10 +62,17 @@ private:
const amrex::AmrCore& m_mesh;
Field& m_ow_velocity;
Field& m_ow_vof;
IntField* m_terrain_blank_ptr{nullptr};

// Should be active unless boundary planes are used
bool m_activate_ow_bndry{true};

// Check for single-phase simulation
bool m_vof_exists{true};

// Check for terrain
bool m_terrain_exists{false};

// Store time corresponding to current boundary data
amrex::Real m_bndry_time{0.0};

Expand Down
55 changes: 49 additions & 6 deletions amr-wind/ocean_waves/boundary_ops/OceanWavesBoundary.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,14 +42,20 @@ void OceanWavesBoundary::post_init_actions()
m_repo.get_field("velocity")
.register_fill_patch_op<OceanWavesFillInflow>(
m_mesh, m_time, *this);
if (m_repo.field_exists("vof")) {
m_vof_exists = m_repo.field_exists("vof");
if (m_vof_exists) {
m_repo.get_field("vof")
.register_fill_patch_op<OceanWavesFillInflow>(
m_mesh, m_time, *this);
m_repo.get_field("density")
.register_fill_patch_op<OceanWavesFillInflow>(
m_mesh, m_time, *this);
}

m_terrain_exists = m_repo.int_field_exists("terrain_blank");
if (m_terrain_exists) {
m_terrain_blank_ptr = &m_repo.get_int_field("terrain_blank");
}
}
}

Expand All @@ -73,6 +79,8 @@ void OceanWavesBoundary::set_velocity(
const int nghost = 1;
const auto& domain = geom.growPeriodicDomain(nghost);

const bool terrain_and_vof_exist = m_terrain_exists && m_vof_exists;

for (amrex::OrientationIter oit; oit != nullptr; ++oit) {
auto ori = oit();
if ((bctype[ori] != BC::mass_inflow) &&
Expand All @@ -85,10 +93,14 @@ void OceanWavesBoundary::set_velocity(
const auto& dbx = ori.isLow() ? amrex::adjCellLo(domain, idir, nghost)
: amrex::adjCellHi(domain, idir, nghost);

auto shift_to_interior =
amrex::IntVect::TheDimensionVector(idir) * (ori.isLow() ? 1 : -1);

for (amrex::MFIter mfi(mfab); mfi.isValid(); ++mfi) {
auto gbx = amrex::grow(mfi.validbox(), nghost);
const auto& bx =
utils::face_aware_boundary_box_intersection(gbx, dbx, ori);
amrex::IntVect shift_to_cc = {0, 0, 0};
const auto& bx = utils::face_aware_boundary_box_intersection(
shift_to_cc, gbx, dbx, ori);
if (!bx.ok()) {
continue;
}
Expand All @@ -98,12 +110,25 @@ void OceanWavesBoundary::set_velocity(
const auto& arr = mfab[mfi].array();
const int numcomp = mfab.nComp();

const auto terrain_blank_flags =
terrain_and_vof_exist
? (*m_terrain_blank_ptr)(lev).const_array(mfi)
: amrex::Array4<int const>();

amrex::ParallelFor(
bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
const amrex::IntVect iv{i, j, k};
for (int n = 0; n < numcomp; n++) {
if (targ_vof(i, j, k) > constants::TIGHT_TOL) {
arr(i, j, k, dcomp + n) =
targ_arr(i, j, k, orig_comp + n);
if (targ_vof(iv) > constants::TIGHT_TOL) {
arr(iv, dcomp + n) = targ_arr(iv, orig_comp + n);
}
if (terrain_and_vof_exist) {
// Terrain-blanked adjacent cell means 0 velocity
if (terrain_blank_flags(
iv + shift_to_cc + shift_to_interior) ==
1) {
arr(iv, dcomp + n) = 0.0;
}
}
}
});
Expand Down Expand Up @@ -225,6 +250,7 @@ void OceanWavesBoundary::set_inflow_sibling_velocity(

BL_PROFILE("amr-wind::OceanWavesBoundary::set_inflow_sibling_velocity");

const bool terrain_and_vof_exist = m_terrain_exists && m_vof_exists;
const auto& bctype = fld.bc_type();
const auto& geom = fld.repo().mesh().Geom(lev);

Expand All @@ -238,6 +264,10 @@ void OceanWavesBoundary::set_inflow_sibling_velocity(

const int idir = ori.coordDir();
const auto& domain_box = geom.Domain();

auto shift_to_interior =
amrex::IntVect::TheDimensionVector(idir) * (ori.isLow() ? 1 : -1);

for (int fdir = 0; fdir < AMREX_SPACEDIM; ++fdir) {

// Only face-normal velocities populated here
Expand Down Expand Up @@ -268,6 +298,11 @@ void OceanWavesBoundary::set_inflow_sibling_velocity(
const auto& targ_arr = m_ow_velocity(lev).const_array(mfi);
const auto& marr = mfab[mfi].array();

const auto terrain_blank_flags =
terrain_and_vof_exist
? (*m_terrain_blank_ptr)(lev).const_array(mfi)
: amrex::Array4<int const>();

amrex::ParallelFor(
bx, [=] AMREX_GPU_DEVICE(
const int i, const int j, const int k) noexcept {
Expand All @@ -277,6 +312,14 @@ void OceanWavesBoundary::set_inflow_sibling_velocity(
if (targ_vof(cc_iv) > constants::TIGHT_TOL) {
marr(i, j, k, 0) = targ_arr(cc_iv, fdir);
}
if (terrain_and_vof_exist) {
// Terrain-blanked boundary-adjacent cell should set
// boundary velocity to 0
if (terrain_blank_flags(
cc_iv + shift_to_interior) == 1) {
marr(i, j, k, 0) = 0.0;
}
}
});
}
}
Expand Down
27 changes: 25 additions & 2 deletions amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,15 @@ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata)
auto& velocity = sim.repo().get_field("velocity");
auto& density = sim.repo().get_field("density");

amr_wind::IntField* terrain_blank_ptr{nullptr};
amr_wind::IntField* terrain_drag_ptr{nullptr};
const bool terrain_exists = sim.repo().int_field_exists("terrain_blank");
// Get fields to prevent forcing in or near underwater terrain
if (terrain_exists) {
terrain_blank_ptr = &sim.repo().get_int_field("terrain_blank");
terrain_drag_ptr = &sim.repo().get_int_field("terrain_drag");
}

for (int lev = 0; lev < nlevels; ++lev) {
const auto& dx = geom[lev].CellSizeArray();
const auto& problo = geom[lev].ProbLoArray();
Expand All @@ -101,6 +110,13 @@ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata)
const auto target_volfrac_arrs = ow_vof(lev).const_arrays();
const auto target_vel_arrs = ow_vel(lev).const_arrays();

const auto terrain_blank_flags =
terrain_exists ? (*terrain_blank_ptr)(lev).const_arrays()
: amrex::MultiArray4<int const>();
const auto terrain_drag_flags =
terrain_exists ? (*terrain_drag_ptr)(lev).const_arrays()
: amrex::MultiArray4<int const>();

const amrex::Real gen_length = wdata.gen_length;
const amrex::Real beach_length = wdata.beach_length;
const amrex::Real beach_length_factor = wdata.beach_length_factor;
Expand All @@ -125,8 +141,15 @@ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata)
const auto target_volfrac = target_volfrac_arrs[nbx];
const auto target_vel = target_vel_arrs[nbx];

bool in_or_near_terrain{false};
if (terrain_exists) {
in_or_near_terrain =
(terrain_blank_flags[nbx](i, j, k) == 1 ||
terrain_drag_flags[nbx](i, j, k) == 1);
}

// Generation region
if (x <= problo[0] + gen_length) {
if (x <= problo[0] + gen_length && !in_or_near_terrain) {
const amrex::Real Gamma =
utils::gamma_generate(x - problo[0], gen_length);
// Get bounded new vof, incorporate with increment
Expand Down Expand Up @@ -164,7 +187,7 @@ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata)
}
}
// Outlet region
if (x + beach_length >= probhi[0]) {
if (x + beach_length >= probhi[0] && !in_or_near_terrain) {
const amrex::Real Gamma = utils::gamma_absorb(
x - (probhi[0] - beach_length), beach_length,
beach_length_factor);
Expand Down
1 change: 1 addition & 0 deletions test/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -254,6 +254,7 @@ add_test_re(ow_linear)
add_test_re(ow_linear_init_waves)
add_test_re(ow_stokes)
add_test_re(ow_current_wing)
add_test_re(ow_current_bathymetry)
add_test_re(scalar_advection_uniform)
add_test_re(scalar_advection_refined)
add_test_re(freestream_bds)
Expand Down
Loading

0 comments on commit 38d1b9f

Please sign in to comment.