Merge branch 'Exawind:main' into mesoscale_zone

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

@@ -34,6 +34,8 @@ private:
     const CFDSim& m_sim;
     const amrex::AmrCore& m_mesh;
     const Field& m_velocity;
+    const Field* m_target_vel{nullptr};
+    const Field* m_target_levelset{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 +50,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

@@ -1,11 +1,80 @@
 #include "amr-wind/equation_systems/icns/source_terms/DragForcing.H"
+#include "amr-wind/equation_systems/vof/volume_fractions.H"
 #include "amr-wind/utilities/IOManager.H"
 
 #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 {
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real viscous_drag_calculations(
+    amrex::Real& Dxz,
+    amrex::Real& Dyz,
+    const amrex::Real ux1r,
+    const amrex::Real uy1r,
+    const amrex::Real ux2r,
+    const amrex::Real uy2r,
+    const amrex::Real z0,
+    const amrex::Real dz,
+    const amrex::Real kappa,
+    const amrex::Real tiny)
+{
+    const amrex::Real m2 = std::sqrt(ux2r * ux2r + uy2r * uy2r);
+    const amrex::Real ustar = m2 * kappa / std::log(1.5 * dz / z0);
+    Dxz += -ustar * ustar * ux1r /
+           (tiny + std::sqrt(ux1r * ux1r + uy1r * uy1r)) / dz;
+    Dyz += -ustar * ustar * uy1r /
+           (tiny + std::sqrt(ux1r * ux1r + uy1r * uy1r)) / dz;
+    return ustar;
+}
+
+// Implementation comes from MOSD approach in boundary_conditions/wall_models/
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void form_drag_calculations(
+    amrex::Real& Dxz,
+    amrex::Real& Dyz,
+    const int i,
+    const int j,
+    const int k,
+    amrex::Array4<amrex::Real const> const& phi,
+    const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx,
+    const amrex::Real ux1r,
+    const amrex::Real uy1r)
+{
+    // phi = eta - z, so eta derivatives in x and y can be calculate with phi
+    amrex::Real n_x, n_y, n_z;
+    amr_wind::multiphase::youngs_finite_difference_normal(
+        i, j, k, phi, n_x, n_y, n_z);
+    // factor of 32 has to do with finite differences, number of points used
+    // negative to make normals point away from waves into air
+    const amrex::Real dx_eta_wave = -n_x / 32. / dx[0];
+    const amrex::Real dy_eta_wave = -n_y / 32. / dx[1];
+    const amrex::Real grad_eta_wave =
+        std::sqrt(dx_eta_wave * dx_eta_wave + dy_eta_wave * dy_eta_wave);
+    n_x = dx_eta_wave / grad_eta_wave;
+    n_y = dy_eta_wave / grad_eta_wave;
+
+    // Relative velocity while considering interface normal
+    const amrex::Real ur_mag =
+        std::sqrt(ux1r * ux1r * n_x * n_x + uy1r * uy1r * n_y * n_y);
+    // Heaviside function changes behavior for velocity surplus/deficit
+    const amrex::Real Heavi_arg = (ux1r * dx_eta_wave + uy1r * dy_eta_wave);
+    const amrex::Real Heavi =
+        (Heavi_arg + std::abs(Heavi_arg)) / (2 * Heavi_arg);
+
+    // Stress in each direction
+    const amrex::Real tau_xz = (1 / M_PI) * ur_mag * ur_mag * grad_eta_wave *
+                               grad_eta_wave * Heavi * n_x;
+    const amrex::Real tau_yz = (1 / M_PI) * ur_mag * ur_mag * grad_eta_wave *
+                               grad_eta_wave * Heavi * n_y;
+
+    // Drag terms from waves added to log law
+    Dxz += -tau_xz / dx[2];
+    Dyz += -tau_yz / dx[2];
+}
+} // namespace
+
 namespace amr_wind::pde::icns {
 
 DragForcing::DragForcing(const CFDSim& sim)
@@ -42,6 +111,16 @@ 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);
+        const auto target_levelset_name =
+            terrain_phys.wave_negative_elevation_name();
+        m_target_levelset = &sim.repo().get_field(target_levelset_name);
+        m_terrain_is_waves = true;
+    }
 }
 
 DragForcing::~DragForcing() = default;
@@ -68,6 +147,15 @@ 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& target_lvs_arr =
+        is_waves ? (*m_target_levelset)(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();
@@ -134,43 +222,75 @@ void DragForcing::operator()(
         amrex::Real Dyz = 0.0;
         amrex::Real bc_forcing_x = 0;
         amrex::Real bc_forcing_y = 0;
+
         const amrex::Real ux1 = vel(i, j, k, 0);
         const amrex::Real uy1 = vel(i, j, k, 1);
         const amrex::Real uz1 = vel(i, j, k, 2);
         const amrex::Real m = std::sqrt(ux1 * ux1 + uy1 * uy1 + uz1 * uz1);
         if (drag(i, j, k) == 1 && (!is_laminar)) {
-            const amrex::Real ux2 = vel(i, j, k + 1, 0);
-            const amrex::Real uy2 = vel(i, j, k + 1, 1);
-            const amrex::Real m2 = std::sqrt(ux2 * ux2 + uy2 * uy2);
+            // Check if close enough to interface to use current cell or below
+            int k_off = -1;
+            if (is_waves) {
+                const amrex::Real cell_length_2D =
+                    std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]);
+                if (target_lvs_arr(i, j, k) + cell_length_2D >= 0) {
+                    // Current cell will be used for wave velocity
+                    k_off = 0;
+                }
+                // Cell below will be used if not (default of -1)
+            }
+            // Establish wall velocity
+            // - estimate wave velocity using target velocity in cells below
+            const amrex::Real wall_u =
+                !is_waves ? 0.0 : target_vel_arr(i, j, k + k_off, 0);
+            const amrex::Real wall_v =
+                !is_waves ? 0.0 : target_vel_arr(i, j, k + k_off, 1);
+            // Relative velocities for calculating shear
+            const amrex::Real ux1r = ux1 - wall_u;
+            const amrex::Real uy1r = uy1 - wall_v;
+            const amrex::Real ux2r = vel(i, j, k + 1, 0) - wall_u;
+            const amrex::Real uy2r = vel(i, j, k + 1, 1) - wall_v;
             const amrex::Real z0 = std::max(terrainz0(i, j, k), z0_min);
-            const amrex::Real ustar = m2 * kappa / std::log(1.5 * dx[2] / z0);
+            const amrex::Real ustar = viscous_drag_calculations(
+                Dxz, Dyz, ux1r, uy1r, ux2r, uy2r, z0, dx[2], kappa, tiny);
+            if (is_waves) {
+                form_drag_calculations(
+                    Dxz, Dyz, i, j, k, target_lvs_arr, dx, ux1r, uy1r);
+            }
             const amrex::Real uTarget =
                 ustar / kappa * std::log(0.5 * dx[2] / z0);
             const amrex::Real uxTarget =
-                uTarget * ux2 / (tiny + std::sqrt(ux2 * ux2 + uy2 * uy2));
+                uTarget * ux2r / (tiny + std::sqrt(ux2r * ux2r + uy2r * uy2r));
             const amrex::Real uyTarget =
-                uTarget * uy2 / (tiny + std::sqrt(ux2 * ux2 + uy2 * uy2));
+                uTarget * uy2r / (tiny + std::sqrt(ux2r * ux2r + uy2r * uy2r));
+            // BC forcing pushes nonrelative velocity toward target velocity
             bc_forcing_x = -(uxTarget - ux1) / dt;
             bc_forcing_y = -(uyTarget - uy1) / dt;
-            Dxz = -ustar * ustar * ux1 /
-                  (tiny + std::sqrt(ux1 * ux1 + uy1 * uy1)) / dx[2];
-            Dyz = -ustar * ustar * uy1 /
-                  (tiny + std::sqrt(ux1 * ux1 + uy1 * uy1)) / dx[2];
         }
+        // Target velocity intended for within terrain
+        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

@@ -72,8 +72,10 @@ private:
     //! Ocean waves target velocity
     Field& m_ow_velocity;
 
-    //! 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());
@@ -34,6 +34,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;
@@ -56,16 +63,19 @@ 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");
     m_ow_bndry->post_init_actions();
-    m_owm->update_relax_zones(m_sim.time().current_time());
+    m_owm->update_target_fields(m_sim.time().current_time());
+    m_owm->update_target_volume_fraction();
     m_ow_bndry->record_boundary_data_time(m_sim.time().current_time());
-    m_owm->apply_relax_zones();
+    if (m_multiphase_mode) {
+        m_owm->apply_relax_zones();
+    }
     m_owm->reset_regrid_flag();
 }
 
@@ -81,7 +91,8 @@ void OceanWaves::pre_advance_work()
     // Update ow values for advection boundaries
     const amrex::Real adv_bdy_time =
         0.5 * (m_sim.time().current_time() + m_sim.time().new_time());
-    m_owm->update_relax_zones(adv_bdy_time);
+    m_owm->update_target_fields(adv_bdy_time);
+    m_owm->update_target_volume_fraction();
     m_ow_bndry->record_boundary_data_time(adv_bdy_time);
 }
 
@@ -90,15 +101,17 @@ void OceanWaves::pre_predictor_work()
     BL_PROFILE("amr-wind::ocean_waves::OceanWaves::pre_predictor_work");
     // Update ow values for boundary fills at new time
     const amrex::Real bdy_fill_time = m_sim.time().new_time();
-    m_owm->update_relax_zones(bdy_fill_time);
+    m_owm->update_target_fields(bdy_fill_time);
+    m_owm->update_target_volume_fraction();
     m_ow_bndry->record_boundary_data_time(bdy_fill_time);
 }
 
 void OceanWaves::post_advance_work()
 {
     BL_PROFILE("amr-wind::ocean_waves::OceanWaves::post_init_actions");
-    // Relax zones are up-to-date because of pre-predictor
-    m_owm->apply_relax_zones();
+    if (m_multiphase_mode) {
+        m_owm->apply_relax_zones();
+    }
     m_owm->reset_regrid_flag();
 }
 

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(const amrex::Real) = 0;
+    virtual void update_target_fields(const amrex::Real) = 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(const amrex::Real time) override
+    void update_target_fields(const amrex::Real time) override
     {
-        ops::UpdateRelaxZonesOp<WaveTheoryTrait>()(m_data, time);
+        ops::UpdateTargetFieldsOp<WaveTheoryTrait>()(m_data, time);
     }
 
     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

@@ -14,11 +14,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;