Terrain blanking informed by ocean waves

amr-wind/equation_systems/icns/source_terms/DragForcing.H

@@ -34,6 +34,7 @@ private:
     const CFDSim& m_sim;
     const amrex::AmrCore& m_mesh;
     const Field& m_velocity;
+    const Field* m_target_vel{nullptr};
     amrex::Gpu::DeviceVector<amrex::Real> m_device_vel_ht;
     amrex::Gpu::DeviceVector<amrex::Real> m_device_vel_vals;
     amrex::Real m_drag_coefficient{10.0};
@@ -48,6 +49,7 @@ private:
     int m_sponge_south{0};
     int m_sponge_north{1};
     bool m_is_laminar{false};
+    bool m_terrain_is_waves{false};
 };
 
 } // namespace amr_wind::pde::icns

amr-wind/equation_systems/icns/source_terms/DragForcing.cpp

@@ -4,6 +4,7 @@
 #include "AMReX_Gpu.H"
 #include "AMReX_Random.H"
 #include "amr-wind/wind_energy/ABL.H"
+#include "amr-wind/physics/TerrainDrag.H"
 #include "amr-wind/utilities/linear_interpolation.H"
 
 namespace amr_wind::pde::icns {
@@ -42,6 +43,13 @@ DragForcing::DragForcing(const CFDSim& sim)
     } else {
         m_sponge_strength = 0.0;
     }
+    if (phy_mgr.contains("OceanWaves")) {
+        const auto terrain_phys =
+            m_sim.physics_manager().get<amr_wind::terraindrag::TerrainDrag>();
+        const auto target_vel_name = terrain_phys.wave_velocity_field_name();
+        m_target_vel = &sim.repo().get_field(target_vel_name);
+        m_terrain_is_waves = true;
+    }
 }
 
 DragForcing::~DragForcing() = default;
@@ -68,6 +76,12 @@ void DragForcing::operator()(
     const auto& drag = (*m_terrain_drag)(lev).const_array(mfi);
     auto* const m_terrainz0 = &this->m_sim.repo().get_field("terrainz0");
     const auto& terrainz0 = (*m_terrainz0)(lev).const_array(mfi);
+
+    const bool is_waves = m_terrain_is_waves;
+    const auto& target_vel_arr = is_waves
+                                     ? (*m_target_vel)(lev).const_array(mfi)
+                                     : amrex::Array4<amrex::Real>();
+
     const auto& geom = m_mesh.Geom(lev);
     const auto& dx = geom.CellSizeArray();
     const auto& prob_lo = geom.ProbLoArray();
@@ -158,20 +172,29 @@ void DragForcing::operator()(
             Dyz = -ustar * ustar * uy1 /
                   (tiny + std::sqrt(ux1 * ux1 + uy1 * uy1)) / dx[2];
         }
+        amrex::Real target_u = 0.;
+        amrex::Real target_v = 0.;
+        amrex::Real target_w = 0.;
+        if (is_waves) {
+            target_u = target_vel_arr(i, j, k, 0);
+            target_v = target_vel_arr(i, j, k, 1);
+            target_w = target_vel_arr(i, j, k, 2);
+        }
+
         const amrex::Real CdM =
             std::min(Cd / (m + tiny), cd_max / scale_factor);
         src_term(i, j, k, 0) -=
-            (CdM * m * ux1 * blank(i, j, k) + Dxz * drag(i, j, k) +
+            (CdM * m * (ux1 - target_u) * blank(i, j, k) + Dxz * drag(i, j, k) +
              bc_forcing_x * drag(i, j, k) +
              (xstart_damping + xend_damping + ystart_damping + yend_damping) *
                  (ux1 - sponge_density * spongeVelX));
         src_term(i, j, k, 1) -=
-            (CdM * m * uy1 * blank(i, j, k) + Dyz * drag(i, j, k) +
+            (CdM * m * (uy1 - target_v) * blank(i, j, k) + Dyz * drag(i, j, k) +
              bc_forcing_y * drag(i, j, k) +
              (xstart_damping + xend_damping + ystart_damping + yend_damping) *
                  (uy1 - sponge_density * spongeVelY));
         src_term(i, j, k, 2) -=
-            (CdM * m * uz1 * blank(i, j, k) +
+            (CdM * m * (uz1 - target_w) * blank(i, j, k) +
              (xstart_damping + xend_damping + ystart_damping + yend_damping) *
                  (uz1 - sponge_density * spongeVelZ));
     });

amr-wind/ocean_waves/OceanWaves.H

@@ -52,6 +52,8 @@ protected:
 private:
     void relaxation_zones();
 
+    void update_target_volume_fraction();
+
     CFDSim& m_sim;
 
     //! Unique pointer to the ocean waves model
@@ -68,6 +70,11 @@ private:
 
     //! Ocean waves pressure gradient
     // Field& m_ow_pressure;
+
+    //! Multiphase mode forces the velocity and vof in the domain
+    //! When off, single-phase mode simply provides the target variables, and
+    //! forcing is handled in other parts of the code
+    bool m_multiphase_mode{true};
 };
 
 } // namespace ocean_waves

amr-wind/ocean_waves/OceanWaves.cpp

@@ -18,7 +18,7 @@ OceanWaves::OceanWaves(CFDSim& sim)
           sim.repo().declare_field("ow_velocity", AMREX_SPACEDIM, 3, 1))
 {
     if (!sim.physics_manager().contains("MultiPhase")) {
-        amrex::Abort("OceanWaves requires Multiphase physics to be active");
+        m_multiphase_mode = false;
     }
     m_ow_levelset.set_default_fillpatch_bc(sim.time());
     m_ow_vof.set_default_fillpatch_bc(sim.time());
@@ -32,6 +32,13 @@ void OceanWaves::pre_init_actions()
     BL_PROFILE("amr-wind::ocean_waves::OceanWaves::pre_init_actions");
     amrex::ParmParse pp(identifier());
 
+    if (!(m_multiphase_mode ||
+          m_sim.physics_manager().contains("TerrainDrag"))) {
+        amrex::Abort(
+            "OceanWaves requires MultiPhase or TerrainDrag physics to be "
+            "active");
+    }
+
     std::string label;
     pp.query("label", label);
     const std::string& tname = label;
@@ -54,13 +61,18 @@ void OceanWaves::pre_init_actions()
 void OceanWaves::initialize_fields(int level, const amrex::Geometry& geom)
 {
     BL_PROFILE("amr-wind::ocean_waves::OceanWaves::initialize_fields");
-    m_owm->init_waves(level, geom);
+    m_owm->init_waves(level, geom, m_multiphase_mode);
 }
 
 void OceanWaves::post_init_actions()
 {
     BL_PROFILE("amr-wind::ocean_waves::OceanWaves::post_init_actions");
-    relaxation_zones();
+    if (m_multiphase_mode) {
+        relaxation_zones();
+    } else {
+        m_owm->update_target_fields(true);
+        m_owm->update_target_volume_fraction();
+    }
 }
 
 void OceanWaves::post_regrid_actions()
@@ -72,25 +84,38 @@ void OceanWaves::post_regrid_actions()
 void OceanWaves::pre_advance_work()
 {
     BL_PROFILE("amr-wind::ocean_waves::OceanWaves::pre_advance_work");
+    if (!m_multiphase_mode) {
+        m_owm->update_target_fields(true);
+        m_owm->update_target_volume_fraction();
+    }
 }
 
 void OceanWaves::post_advance_work()
 {
     BL_PROFILE("amr-wind::ocean_waves::OceanWaves::post_init_actions");
-    relaxation_zones();
+    if (m_multiphase_mode) {
+        relaxation_zones();
+    }
 }
 
 /** Update ocean waves relaxation zones
  *
  */
 void OceanWaves::relaxation_zones()
 {
-    BL_PROFILE("amr-wind::ocean_waves::OceanWaves::update_relaxation_zones");
-    m_owm->update_relax_zones();
+    BL_PROFILE("amr-wind::ocean_waves::OceanWaves::relaxation_zones");
+    m_owm->update_target_fields(false);
     m_owm->apply_relax_zones();
     m_owm->reset_regrid_flag();
 }
 
+void OceanWaves::update_target_volume_fraction()
+{
+    BL_PROFILE(
+        "amr-wind::ocean_waves::OceanWaves::update_target_volume_fraction");
+    m_owm->update_target_volume_fraction();
+}
+
 void OceanWaves::prepare_outputs()
 {
     const std::string post_dir = m_sim.io_manager().post_processing_directory();

amr-wind/ocean_waves/OceanWavesModel.H

@@ -32,12 +32,14 @@ public:
 
     virtual void read_inputs(const ::amr_wind::utils::MultiParser&) = 0;
 
-    virtual void init_waves(int, const amrex::Geometry&) = 0;
+    virtual void init_waves(int, const amrex::Geometry&, bool) = 0;
 
-    virtual void update_relax_zones() = 0;
+    virtual void update_target_fields(bool) = 0;
 
     virtual void apply_relax_zones() = 0;
 
+    virtual void update_target_volume_fraction() = 0;
+
     virtual void prepare_outputs(const std::string&) = 0;
 
     virtual void write_outputs() = 0;
@@ -91,16 +93,21 @@ public:
         m_out_op.read_io_options(pp);
     }
 
-    void update_relax_zones() override
+    void update_target_fields(bool for_forcing_term) override
     {
-        ops::UpdateRelaxZonesOp<WaveTheoryTrait>()(m_data);
+        ops::UpdateTargetFieldsOp<WaveTheoryTrait>()(m_data, for_forcing_term);
     }
 
     void apply_relax_zones() override
     {
         ops::ApplyRelaxZonesOp<WaveTheoryTrait>()(m_data);
     }
 
+    void update_target_volume_fraction() override
+    {
+        ops::UpdateTargetVolumeFractionOp<WaveTheoryTrait>()(m_data);
+    }
+
     void prepare_outputs(const std::string& out_dir) override
     {
         m_out_op.prepare_outputs(out_dir);
@@ -112,9 +119,11 @@ public:
 
     void write_outputs() override { m_out_op.write_outputs(); }
 
-    void init_waves(int level, const amrex::Geometry& geom) override
+    void init_waves(
+        int level, const amrex::Geometry& geom, bool multiphase_mode) override
     {
-        ops::InitDataOp<WaveTheoryTrait>()(m_data, level, geom);
+        ops::InitDataOp<WaveTheoryTrait>()(
+            m_data, level, geom, multiphase_mode);
     }
 };
 

amr-wind/ocean_waves/OceanWavesOps.H

@@ -32,11 +32,14 @@ template <typename WaveTheoryTrait, typename = void>
 struct InitDataOp;
 
 template <typename WaveTheoryTrait, typename = void>
-struct UpdateRelaxZonesOp;
+struct UpdateTargetFieldsOp;
 
 template <typename WaveTheoryTrait, typename = void>
 struct ApplyRelaxZonesOp;
 
+template <typename WaveTheoryTrait, typename = void>
+struct UpdateTargetVolumeFractionOp;
+
 template <typename WaveTheoryTrait, typename = void>
 struct ProcessOutputsOp;