Improve Ocean Waves Fill (Exawind#1440)

* set up a set_inflow_sibling operation for ocean waves boundary: the set_velocity was not modifying the right cells
* switch to tight tol for vof checking
* unit test for new capability

amr-wind/equation_systems/icns/icns_advection.cpp

@@ -177,10 +177,8 @@ void MacProjOp::set_inflow_velocity(amrex::Real time)
         velocity.set_inflow_sibling_fields(lev, time, mac_vec);
         if (m_phy_mgr.contains("OceanWaves")) {
             auto& ow = m_phy_mgr.get<amr_wind::ocean_waves::OceanWaves>();
-            for (int dir = 0; dir < ICNS::ndim; ++dir) {
-                ow.ow_bndry().set_velocity(
-                    lev, time, velocity, *mac_vec[dir], 0, dir);
-            }
+            ow.ow_bndry().set_inflow_sibling_velocity(
+                lev, time, velocity, mac_vec);
         }
     }
 }

amr-wind/ocean_waves/boundary_ops/OceanWavesBoundary.H

@@ -42,6 +42,12 @@ public:
         const Field& fld,
         amrex::MultiFab& mfab) const;
 
+    void set_inflow_sibling_velocity(
+        const int lev,
+        const amrex::Real time,
+        const Field& fld,
+        amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM> mfabs) const;
+
     void record_boundary_data_time(const amrex::Real time)
     {
         m_bndry_time = time;

amr-wind/ocean_waves/boundary_ops/OceanWavesBoundary.cpp

@@ -101,7 +101,7 @@ void OceanWavesBoundary::set_velocity(
             amrex::ParallelFor(
                 bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                     for (int n = 0; n < numcomp; n++) {
-                        if (targ_vof(i, j, k) > constants::LOOSE_TOL) {
+                        if (targ_vof(i, j, k) > constants::TIGHT_TOL) {
                             arr(i, j, k, dcomp + n) =
                                 targ_arr(i, j, k, orig_comp + n);
                         }
@@ -212,4 +212,75 @@ void OceanWavesBoundary::set_density(
     }
 }
 
+void OceanWavesBoundary::set_inflow_sibling_velocity(
+    const int lev,
+    const amrex::Real /*time*/,
+    const Field& fld,
+    const amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM> mfabs) const
+{
+
+    if (!m_activate_ow_bndry) {
+        return;
+    }
+
+    BL_PROFILE("amr-wind::OceanWavesBoundary::set_inflow_sibling_velocity");
+
+    const auto& bctype = fld.bc_type();
+    const auto& geom = fld.repo().mesh().Geom(lev);
+
+    for (amrex::OrientationIter oit; oit != nullptr; ++oit) {
+        const auto ori = oit();
+        if ((bctype[ori] != BC::mass_inflow) &&
+            (bctype[ori] != BC::mass_inflow_outflow) &&
+            (bctype[ori] != BC::wave_generation)) {
+            continue;
+        }
+
+        const int idir = ori.coordDir();
+        const auto& domain_box = geom.Domain();
+        for (int fdir = 0; fdir < AMREX_SPACEDIM; ++fdir) {
+
+            // Only face-normal velocities populated here
+            if (idir != fdir) {
+                continue;
+            }
+            const auto& dbx = ori.isLow() ? amrex::bdryLo(domain_box, idir)
+                                          : amrex::bdryHi(domain_box, idir);
+
+            // Shift from valid face index to first cell-centered ghost
+            amrex::IntVect shift_to_cc = {0, 0, 0};
+            if (ori.isLow()) {
+                --shift_to_cc[fdir];
+            }
+
+            auto& mfab = *mfabs[fdir];
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+            for (amrex::MFIter mfi(mfab); mfi.isValid(); ++mfi) {
+                const auto& vbx = mfi.validbox();
+                const auto& bx = vbx & dbx;
+                if (!bx.ok()) {
+                    continue;
+                }
+
+                const auto& targ_vof = m_ow_vof(lev).const_array(mfi);
+                const auto& targ_arr = m_ow_velocity(lev).const_array(mfi);
+                const auto& marr = mfab[mfi].array();
+
+                amrex::ParallelFor(
+                    bx, [=] AMREX_GPU_DEVICE(
+                            const int i, const int j, const int k) noexcept {
+                        amrex::IntVect cc_iv = {i, j, k};
+                        cc_iv += shift_to_cc;
+
+                        if (targ_vof(cc_iv) > constants::TIGHT_TOL) {
+                            marr(i, j, k, 0) = targ_arr(cc_iv, fdir);
+                        }
+                    });
+            }
+        }
+    }
+}
+
 } // namespace amr_wind

unit_tests/ocean_waves/test_relaxation_zones.cpp

@@ -5,6 +5,7 @@
 #include "amr-wind/ocean_waves/OceanWaves.H"
 #include "amr-wind/utilities/constants.H"
 #include "amr-wind/physics/multiphase/MultiPhase.H"
+#include "amr-wind/equation_systems/icns/icns_advection.H"
 
 namespace amr_wind_tests {
 
@@ -203,7 +204,7 @@ void make_target_velocity(
             amrex::ParallelFor(
                 gbx, velocity.num_comp(),
                 [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
-                    if (ow_vof_arr(i, j, k) <= amr_wind::constants::LOOSE_TOL) {
+                    if (ow_vof_arr(i, j, k) <= amr_wind::constants::TIGHT_TOL) {
                         ow_vel_arr(i, j, k, n) = vel_arr(i, j, k, n);
                     }
                 });
@@ -264,6 +265,34 @@ amrex::Real bdy_error(amr_wind::Field& comp, amr_wind::Field& targ)
     return bdy_error(comp, targ, 1);
 }
 
+amrex::Real uface_bdy_error(amr_wind::Field& comp, amr_wind::Field& targ)
+{
+    amrex::Real error_total = 0.0;
+
+    for (int lev = 0; lev < comp.repo().num_active_levels(); ++lev) {
+        error_total += amrex::ReduceSum(
+            comp(lev), targ(lev), 0,
+            [=] AMREX_GPU_HOST_DEVICE(
+                amrex::Box const& bx,
+                amrex::Array4<amrex::Real const> const& comp_arr,
+                amrex::Array4<amrex::Real const> const& targ_arr)
+                -> amrex::Real {
+                amrex::Real error = 0.0;
+
+                amrex::Loop(bx, [=, &error](int i, int j, int k) noexcept {
+                    if (i == 0) {
+                        error += std::abs(
+                            comp_arr(i, j, k, 0) - targ_arr(i - 1, j, k, 0));
+                    }
+                });
+
+                return error;
+            });
+    }
+    amrex::ParallelDescriptor::ReduceRealSum(error_total);
+    return error_total;
+}
+
 } // namespace
 
 TEST_F(OceanWavesOpTest, relaxation_zone)
@@ -427,4 +456,78 @@ TEST_F(OceanWavesOpTest, boundary_fill)
     EXPECT_NEAR(error_total, 0.0, tol);
 }
 
+TEST_F(OceanWavesOpTest, set_inflow_sibling)
+{
+
+    constexpr double tol = 1.0e-3;
+    const amrex::Vector<amrex::Real> gas_vel{{1.0, 0.0, 0.0}};
+
+    populate_parameters();
+    {
+        // Ocean Waves details
+        amrex::ParmParse pp("OceanWaves");
+        pp.add("label", (std::string) "lin_ow");
+        amrex::ParmParse ppow("OceanWaves.lin_ow");
+        ppow.add("type", (std::string) "LinearWaves");
+        ppow.add("wave_height", 0.05);
+        ppow.add("wave_length", 1.0);
+        ppow.add("water_depth", 1.0);
+        // Wave generation and numerical beach
+        ppow.add("relax_zone_gen_length", 2.0);
+        ppow.add("numerical_beach_length", 4.0);
+    }
+    {
+        amrex::ParmParse pp("time");
+        pp.add("fixed_dt", 0.1);
+    }
+    {
+        // Boundary conditions
+        amrex::ParmParse ppxlo("xlo");
+        ppxlo.add("type", (std::string) "wave_generation");
+        ppxlo.addarr("velocity", gas_vel);
+        ppxlo.add("vof", 0.0);
+        ppxlo.add("density", 1.0);
+    }
+
+    initialize_mesh();
+
+    // ICNS must be initialized for MultiPhase physics, which is needed for
+    // OceanWaves
+    auto& pde_mgr = sim().pde_manager();
+    pde_mgr.register_icns();
+    // Initialize physics
+    sim().init_physics();
+    auto& oceanwaves =
+        sim().physics_manager().get<amr_wind::ocean_waves::OceanWaves>();
+    // Initialize fields
+    auto& repo = sim().repo();
+    auto& velocity = repo.get_field("velocity");
+    velocity.setVal(gas_vel, 3);
+    oceanwaves.pre_init_actions();
+    for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
+        oceanwaves.initialize_fields(lev, mesh().Geom(lev));
+    }
+    // Do post-init step, which includes fillpatch calls
+    oceanwaves.post_init_actions();
+
+    // Get MAC velocity in x
+    auto& u_mac = repo.get_field("u_mac");
+    // Set to 0 as a starting point
+    u_mac.setVal(0.0);
+    // Initialize MAC projection operator
+    auto mco = amr_wind::pde::MacProjOp(
+        sim().repo(), sim().physics_manager(), false, false, false, false);
+    // Populate boundary using set inflow
+    mco.set_inflow_velocity(0.0);
+
+    // Get fields for comparison
+    auto& ow_vof = repo.get_field("ow_vof");
+    auto& ow_velocity = repo.get_field("ow_velocity");
+
+    // Check velocity field to confirm not modified
+    make_target_velocity(ow_velocity, velocity, ow_vof);
+    const amrex::Real error_total = uface_bdy_error(u_mac, ow_velocity);
+    EXPECT_NEAR(error_total, 0.0, tol);
+}
+
 } // namespace amr_wind_tests