Merge branch 'main' into ci-fix

.github/workflows/ci.yml

@@ -178,11 +178,11 @@ jobs:
       - name: Prepare SyCL environment
         run: |
           export DEBIAN_FRONTEND=noninteractive
-          sudo wget https://apt.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS-2023.PUB
-          sudo apt-key add GPG-PUB-KEY-INTEL-SW-PRODUCTS-2023.PUB
+          sudo wget https://apt.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS.PUB
+          sudo apt-key add GPG-PUB-KEY-INTEL-SW-PRODUCTS.PUB
           echo "deb https://apt.repos.intel.com/oneapi all main" | sudo tee /etc/apt/sources.list.d/oneAPI.list
           sudo apt-get update
-          sudo apt-get install -y --no-install-recommends ninja-build intel-oneapi-dpcpp-cpp-compiler intel-oneapi-mkl-devel
+          sudo apt-get install -y --no-install-recommends ninja-build intel-oneapi-compiler-dpcpp-cpp intel-oneapi-mkl-devel
       - name: Configure and build
         run: |
           set +e

amr-wind/CFDSim.cpp

@@ -23,18 +23,19 @@ CFDSim::~CFDSim() = default;
 
 void CFDSim::create_turbulence_model()
 {
-    std::string transport_model = "ConstTransport";
-    std::string turbulence_model = "Laminar";
+    std::string transport_model_name = "ConstTransport";
+    std::string turbulence_model_name = "Laminar";
     {
         amrex::ParmParse pp("transport");
-        pp.query("model", transport_model);
+        pp.query("model", transport_model_name);
     }
     {
         amrex::ParmParse pp("turbulence");
-        pp.query("model", turbulence_model);
+        pp.query("model", turbulence_model_name);
     }
 
-    const std::string identifier = turbulence_model + "-" + transport_model;
+    const std::string identifier =
+        turbulence_model_name + "-" + transport_model_name;
     m_turbulence = turbulence::TurbulenceModel::create(identifier, *this);
     m_turbulence->parse_model_coeffs();
 }
@@ -45,7 +46,7 @@ void CFDSim::init_physics()
     amrex::Vector<std::string> phys_names;
     pp.queryarr("physics", phys_names);
 
-    for (auto& phy : phys_names) {
+    for (const auto& phy : phys_names) {
         m_physics_mgr.create(phy, *this);
     }
 }

amr-wind/boundary_conditions/BCInterface.cpp

@@ -131,6 +131,7 @@ std::pair<const std::string, const std::string> BCIface::get_dirichlet_udfs()
                         "faces");
                 } else {
                     inflow_udf = val;
+                    has_inflow_udf = true;
                 }
             }
         }
@@ -146,6 +147,7 @@ std::pair<const std::string, const std::string> BCIface::get_dirichlet_udfs()
                         "faces");
                 } else {
                     wall_udf = val;
+                    has_wall_udf = true;
                 }
             }
         }

amr-wind/boundary_conditions/wall_models/LogLaw.H

@@ -0,0 +1,59 @@
+#ifndef LogLaw_H
+#define LogLaw_H
+
+#include "AMReX_AmrCore.H"
+#include <AMReX.H>
+namespace amr_wind {
+struct LogLaw
+{
+    /*
+     * A simple wall model that sets the wall-shear stress
+     * based on computing u_tau given the horizontal velocity
+     * magnitude at a zref. This is akin to an explicit non-linear
+     * Robin boundary condition at the wall.
+     */
+
+    // Log law constants from Lee & Moser 2015
+    // https://doi.org/10.1017/jfm.2015.268.
+    amrex::Real B{4.27};
+    amrex::Real kappa{0.384};
+    int max_iters = 25; // Max iterations for u_tau Newton-Raphson solve
+    // Reference height for log law
+    amrex::Real zref;
+    int ref_index{0};
+    amrex::Real nu; // molecular viscosity
+    // u_tau state variable, gets updated in update_utau depending on
+    // the type of wall model used
+    amrex::Real utau_mean{1.0};
+    amrex::Real wspd_mean; // mean horizontal velocity magnitude
+
+    void update_utau_mean() { utau_mean = get_utau(wspd_mean); }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
+    get_utau(amrex::Real wspd) const
+    {
+        amrex::Real utau_iter = -1;
+        amrex::Real wspd_pred;
+        amrex::Real wspd_deriv;
+        amrex::Real zplus;
+        amrex::Real utau = utau_mean;
+        int iter = 0;
+        while ((std::abs(utau_iter - utau) > 1e-5) && iter <= max_iters) {
+            utau_iter = utau;
+            zplus = zref * utau / nu;
+            // Get wspd for a given utau from log-law
+            wspd_pred = utau * (std::log(zplus) / kappa + B);
+            wspd_deriv = (1 + std::log(zplus)) / kappa + B; // d(wspd)/d(utau)
+            utau =
+                utau - (wspd_pred - wspd) / wspd_deriv; // Newton-Raphson update
+            ++iter;
+        }
+        if (iter == max_iters) {
+            amrex::Abort();
+        }
+        return utau;
+    }
+};
+} /* namespace amr_wind */
+
+#endif /* LogLaw_H */

amr-wind/boundary_conditions/wall_models/ShearStressSimple.H

@@ -0,0 +1,40 @@
+#ifndef SHEARSTRESSSIMPLE_H
+#define SHEARSTRESSSIMPLE_H
+
+#include "amr-wind/boundary_conditions/wall_models/LogLaw.H"
+#include "amr-wind/wind_energy/ShearStress.H"
+
+namespace amr_wind {
+
+struct SimpleShearSchumann
+{
+    explicit SimpleShearSchumann(const amr_wind::LogLaw& ll)
+        : utau2(ll.utau_mean * ll.utau_mean), wspd_mean(ll.wspd_mean), m_ll(ll)
+    {}
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
+    get_shear(amrex::Real u, amrex::Real /* wspd */) const
+    {
+        return u / wspd_mean * utau2;
+    };
+
+    amrex::Real utau2;
+    amrex::Real wspd_mean;
+    const amr_wind::LogLaw m_ll;
+};
+struct SimpleShearLogLaw
+{
+    explicit SimpleShearLogLaw(const amr_wind::LogLaw& ll) : m_ll(ll) {}
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
+    get_shear(amrex::Real u, amrex::Real wspd) const
+    {
+        amrex::Real utau = m_ll.get_utau(wspd);
+        return utau * utau * u / wspd;
+    };
+
+    const amr_wind::LogLaw m_ll;
+};
+} // namespace amr_wind
+
+#endif /* SHEARSTRESSSIMPLE_H */

amr-wind/boundary_conditions/wall_models/WallFunction.H

@@ -3,6 +3,8 @@
 
 #include "amr-wind/CFDSim.H"
 #include "amr-wind/core/FieldBCOps.H"
+#include "amr-wind/utilities/FieldPlaneAveragingFine.H"
+#include "amr-wind/boundary_conditions/wall_models/LogLaw.H"
 
 namespace amr_wind {
 
@@ -16,9 +18,14 @@ namespace amr_wind {
 class WallFunction
 {
 public:
-    explicit WallFunction(const CFDSim& sim);
+    explicit WallFunction(CFDSim& sim);
 
-    amrex::Real utau() const { return m_utau; }
+    amrex::Real utau() const { return m_log_law.utau_mean; }
+    LogLaw log_law() const { return m_log_law; }
+
+    //! Update the mean velocity at a given timestep
+    void update_umean();
+    void update_utau_mean();
 
     ~WallFunction() = default;
 
@@ -27,7 +34,10 @@ private:
 
     const amrex::AmrCore& m_mesh;
 
-    amrex::Real m_utau{0.0};
+    //! LogLaw instance
+    LogLaw m_log_law;
+    int m_direction{2}; ///< Direction normal to wall, hardcoded to z
+    VelPlaneAveragingFine m_pa_vel;
 };
 
 /** Applies a shear-stress value at the domain boundary
@@ -38,15 +48,21 @@ private:
 class VelWallFunc : public FieldBCIface
 {
 public:
-    VelWallFunc(Field& velocity, const WallFunction& wall_func);
+    VelWallFunc(Field& velocity, 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);
 
+    template <typename ShearStressSimple>
+    static void wall_model(
+        Field& velocity,
+        const FieldState rho_state,
+        const ShearStressSimple& tau);
+
 private:
-    const WallFunction& m_wall_func;
+    WallFunction& m_wall_func;
     std::string m_wall_shear_stress_type{"constant"};
 };
 } // namespace amr_wind