diff --git a/amr-wind/equation_systems/icns/icns_advection.cpp b/amr-wind/equation_systems/icns/icns_advection.cpp index 417c401551..67d06d7c69 100644 --- a/amr-wind/equation_systems/icns/icns_advection.cpp +++ b/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(); - 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); } } } diff --git a/amr-wind/ocean_waves/boundary_ops/OceanWavesBoundary.H b/amr-wind/ocean_waves/boundary_ops/OceanWavesBoundary.H index f64f93ea0c..c97d218c1d 100644 --- a/amr-wind/ocean_waves/boundary_ops/OceanWavesBoundary.H +++ b/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 mfabs) const; + void record_boundary_data_time(const amrex::Real time) { m_bndry_time = time; diff --git a/amr-wind/ocean_waves/boundary_ops/OceanWavesBoundary.cpp b/amr-wind/ocean_waves/boundary_ops/OceanWavesBoundary.cpp index 4570fa5e1c..aeb85d0a93 100644 --- a/amr-wind/ocean_waves/boundary_ops/OceanWavesBoundary.cpp +++ b/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 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 diff --git a/unit_tests/ocean_waves/test_relaxation_zones.cpp b/unit_tests/ocean_waves/test_relaxation_zones.cpp index 5b148796e2..69663198a4 100644 --- a/unit_tests/ocean_waves/test_relaxation_zones.cpp +++ b/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 const& comp_arr, + amrex::Array4 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 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(); + // 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