Implementing a wall modeled LES channel flow (#836)

Co-authored-by: moprak-nrel <120606615+moprak-nrel@users.noreply.github.com>

amr-wind/boundary_conditions/CMakeLists.txt

@@ -4,3 +4,5 @@ target_sources(${amr_wind_lib_name}
       BCInterface.cpp
       FixedGradientBC.cpp
    )
+
+add_subdirectory(wall_models)

amr-wind/boundary_conditions/wall_models/CMakeLists.txt

@@ -0,0 +1,5 @@
+target_sources(${amr_wind_lib_name}
+   PRIVATE
+      #C++
+      WallFunction.cpp
+   )

amr-wind/boundary_conditions/wall_models/WallFunction.H

@@ -0,0 +1,54 @@
+#ifndef WALLFUNCTION_H
+#define WALLFUNCTION_H
+
+#include "amr-wind/CFDSim.H"
+#include "amr-wind/core/FieldBCOps.H"
+
+namespace amr_wind {
+
+/** Wall-function computations for LES simulations
+ *  \ingroup utilities
+ *
+ *  This class performs the necessary computations at the beginning of
+ *  predictor/corrector steps. The actual BC population in ghost cells is
+ *  performed by VelWallFunc BC interface class.
+ */
+class WallFunction
+{
+public:
+    explicit WallFunction(const CFDSim& sim);
+
+    amrex::Real utau() const { return m_utau; }
+
+    ~WallFunction() = default;
+
+private:
+    const CFDSim& m_sim;
+
+    const amrex::AmrCore& m_mesh;
+
+    amrex::Real m_utau{0.0};
+};
+
+/** Applies a shear-stress value at the domain boundary
+ *  \ingroup field_bc utilities
+ *
+ *  \sa WallFunction
+ */
+class VelWallFunc : public FieldBCIface
+{
+public:
+    VelWallFunc(Field& velocity, const WallFunction& wall_func);
+
+    void operator()(Field& velocity, const FieldState rho_state) override;
+
+    static void wall_model(
+        Field& velocity, const FieldState rho_state, const amrex::Real utau);
+
+private:
+    const WallFunction& m_wall_func;
+    std::string m_wall_shear_stress_type{"constant"};
+};
+} // namespace amr_wind
+
+#endif /* WALLFUNCTION_H */

amr-wind/boundary_conditions/wall_models/WallFunction.cpp

@@ -0,0 +1,125 @@
+#include "amr-wind/boundary_conditions/wall_models/WallFunction.H"
+#include "amr-wind/utilities/tensor_ops.H"
+#include "amr-wind/utilities/trig_ops.H"
+#include "amr-wind/diffusion/diffusion.H"
+
+#include <cmath>
+
+#include "AMReX_ParmParse.H"
+#include "AMReX_Print.H"
+#include "AMReX_ParallelDescriptor.H"
+
+namespace amr_wind {
+
+WallFunction::WallFunction(const CFDSim& sim) : m_sim(sim), m_mesh(m_sim.mesh())
+{
+    {
+        amrex::ParmParse pp("BodyForce");
+        amrex::Vector<amrex::Real> body_force{{0.0, 0.0, 0.0}};
+        pp.getarr("magnitude", body_force);
+        m_utau = std::sqrt(body_force[0]);
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+            std::abs(body_force[1]) < 1e-16,
+            "body force in y should be zero for this wall function");
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+            std::abs(body_force[2]) < 1e-16,
+            "body force in z should be zero for this wall function");
+    }
+}
+
+VelWallFunc::VelWallFunc(Field& /*unused*/, const WallFunction& wall_func)
+    : m_wall_func(wall_func)
+{
+    amrex::ParmParse pp("WallFunction");
+    pp.query("wall_shear_stress_type", m_wall_shear_stress_type);
+    m_wall_shear_stress_type = amrex::toLower(m_wall_shear_stress_type);
+
+    if (m_wall_shear_stress_type == "constant") {
+        amrex::Print() << "Shear Stress model: " << m_wall_shear_stress_type
+                       << std::endl;
+    } else {
+        amrex::Abort("Shear Stress wall model input mistake");
+    }
+}
+
+void VelWallFunc::wall_model(
+    Field& velocity, const FieldState rho_state, const amrex::Real utau)
+{
+    BL_PROFILE("amr-wind::VelWallFunc");
+    constexpr int idim = 2;
+    const auto& repo = velocity.repo();
+    const auto& density = repo.get_field("density", rho_state);
+    const auto& viscosity = repo.get_field("velocity_mueff");
+    const int nlevels = repo.num_active_levels();
+
+    amrex::Orientation zlo(amrex::Direction::z, amrex::Orientation::low);
+    amrex::Orientation zhi(amrex::Direction::z, amrex::Orientation::high);
+    if ((velocity.bc_type()[zlo] != BC::wall_model) &&
+        (velocity.bc_type()[zhi] != BC::wall_model)) {
+        return;
+    }
+
+    for (int lev = 0; lev < nlevels; ++lev) {
+        const auto& geom = repo.mesh().Geom(lev);
+        const auto& domain = geom.Domain();
+        amrex::MFItInfo mfi_info{};
+
+        const auto& rho_lev = density(lev);
+        auto& vel_lev = velocity(lev);
+        const auto& eta_lev = viscosity(lev);
+
+        if (amrex::Gpu::notInLaunchRegion()) {
+            mfi_info.SetDynamic(true);
+        }
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+        for (amrex::MFIter mfi(vel_lev, mfi_info); mfi.isValid(); ++mfi) {
+            const auto& bx = mfi.validbox();
+            auto varr = vel_lev.array(mfi);
+            auto den = rho_lev.array(mfi);
+            auto eta = eta_lev.array(mfi);
+
+            if (bx.smallEnd(idim) == domain.smallEnd(idim) &&
+                velocity.bc_type()[zlo] == BC::wall_model) {
+                amrex::ParallelFor(
+                    amrex::bdryLo(bx, idim),
+                    [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                        const amrex::Real mu = eta(i, j, k);
+
+                        // Dirichlet BC
+                        varr(i, j, k - 1, 2) = 0.0;
+
+                        // Shear stress BC
+                        varr(i, j, k - 1, 0) = utau * utau / mu * den(i, j, k);
+                        varr(i, j, k - 1, 1) = 0.0;
+                    });
+            }
+
+            if (bx.bigEnd(idim) == domain.bigEnd(idim) &&
+                velocity.bc_type()[zhi] == BC::wall_model) {
+                amrex::ParallelFor(
+                    amrex::bdryHi(bx, idim),
+                    [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                        const amrex::Real mu = eta(i, j, k - 1);
+
+                        // Dirichlet BC
+                        varr(i, j, k, 2) = 0.0;
+
+                        // Shear stress BC
+                        varr(i, j, k, 0) = utau * utau / mu * den(i, j, k);
+                        varr(i, j, k, 1) = 0.0;
+                    });
+            }
+        }
+    }
+}
+
+void VelWallFunc::operator()(Field& velocity, const FieldState rho_state)
+{
+    if (m_wall_shear_stress_type == "constant") {
+        wall_model(velocity, rho_state, m_wall_func.utau());
+    }
+}
+
+} // namespace amr_wind

amr-wind/diffusion/diffusion.H

@@ -7,23 +7,6 @@
 
 namespace diffusion {
 
-void wall_model_bc_moeng(
-    amr_wind::Field& velocity,
-    const amrex::Real utau,
-    const amr_wind::FieldState fstate,
-    const amrex::FArrayBox& instplanar);
-
-void wall_model_bc(
-    amr_wind::Field& velocity,
-    const amrex::Real utau,
-    const amrex::Real umag,
-    const amr_wind::FieldState fstate);
-
-void temp_wall_model_bc(
-    amr_wind::Field& temperature,
-    const amrex::FArrayBox& instplanar,
-    const amr_wind::FieldState fstate);
-
 amrex::Vector<amrex::Array<amrex::LinOpBCType, AMREX_SPACEDIM>>
 get_diffuse_tensor_bc(
     amr_wind::Field& velocity, amrex::Orientation::Side side) noexcept;

amr-wind/physics/ChannelFlow.H

@@ -4,8 +4,36 @@
 #include "amr-wind/core/Physics.H"
 #include "amr-wind/core/Field.H"
 #include "amr-wind/CFDSim.H"
+#include "amr-wind/boundary_conditions/wall_models/WallFunction.H"
+#include "amr-wind/turbulence/TurbulenceModel.H"
 #include <AMReX_REAL.H>
 
+namespace {
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
+analytical_smagorinsky_profile(
+    const amrex::Real h,
+    const amrex::Real Cs,
+    const amrex::Real dx,
+    const amrex::Real rho,
+    const amrex::Real mu,
+    const amrex::Real dpdx,
+    const amrex::Real C0,
+    const amrex::Real C1)
+{
+    const amrex::Real Cs2 = Cs * Cs;
+    const amrex::Real dx2 = dx * dx;
+    const amrex::Real Cs2dx2 = Cs2 * dx2;
+    return -mu / (rho * Cs2dx2) * h +
+           Cs2dx2 / (3 * dpdx) *
+               std::pow(
+                   (2.0 / Cs2dx2 * dpdx * h +
+                    (mu / (rho * Cs2dx2)) * (mu / (rho * Cs2dx2)) + C1),
+                   3.0 / 2.0) +
+           C0;
+}
+
+} // namespace
+
 namespace amr_wind::channel_flow {
 
 /** Channel Flow physics
@@ -33,6 +61,8 @@ public:
     template <typename IndexSelector>
     amrex::Real compute_error(const IndexSelector& idxOp);
 
+    amrex::Real compute_analytical_smagorinsky_error();
+
     void output_error();
 
     void post_init_actions() override;
@@ -44,10 +74,15 @@ public:
     void post_advance_work() override;
 
 private:
+    //! CFD simulation controller instance
+    CFDSim& m_sim;
+
     const amr_wind::SimTime& m_time;
-    const FieldRepo& m_repo;
+    FieldRepo& m_repo;
     const amrex::AmrCore& m_mesh;
 
+    WallFunction m_wall_func;
+
     //! Wall normal directon - Default y direction
     int m_norm_dir{1};
 
@@ -72,9 +107,24 @@ private:
     //! initial sdr value
     amrex::Real m_sdr0{1000.0};
 
+    //! flag for perturbations of the initial condition
+    bool m_perturb_vel{false};
+
+    //! perturbation period (y-direction)
+    amrex::Real m_perturb_y_period{1.0};
+
+    //! perturbation period (z-direction)
+    amrex::Real m_perturb_z_period{1.0};
+
+    //! perturbation factor (fraction of utau)
+    amrex::Real m_perturb_fac{0.1};
+
     //! Von-Karman constant
     amrex::Real m_kappa{0.41};
 
+    //! Turbulence model
+    std::string m_turbulence_model;
+
     bool m_laminar{false};
 
     bool m_mesh_mapping{false};
@@ -83,6 +133,18 @@ private:
 
     bool m_half{false};
 
+    //! flag for analytical smagorinsky test
+    bool m_analytical_smagorinsky_test{false};
+
+    //! Analytical Smagorinsky first coefficient
+    amrex::Real m_C0{0.0};
+
+    //! Analytical Smagorinsky second coefficient
+    amrex::Real m_C1{0.0};
+
+    //! Body forcing (x direction)
+    amrex::Real m_dpdx{0.0};
+
     //! direction of mean velocity
     int m_mean_vel_dir{0};