Underwater terrain-aware wave forcing (Exawind#1475)

* turn off apply_relaxation_zones where terrain blanking (or terrain drag forcing) happens

* modify ow boundary for underwater terrain

* add regression test

Author: mbkuhn

amr-wind/ocean_waves/boundary_ops/OceanWavesBoundary.H

@@ -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>
@@ -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};
 

amr-wind/ocean_waves/boundary_ops/OceanWavesBoundary.cpp

@@ -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");
+        }
     }
 }
 
@@ -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) &&
@@ -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;
             }
@@ -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;
+                            }
                         }
                     }
                 });
@@ -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);
 
@@ -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
@@ -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 {
@@ -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;
+                            }
+                        }
                     });
             }
         }

amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp

@@ -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();
@@ -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;
@@ -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
@@ -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);

test/CMakeLists.txt

@@ -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)