add implementaiton for a log-law based wall model (#847)

* add implementaiton for a log-law based wall model

* make clang-format happy

* make LogLaw header only and add gpu fix

* linter fix

* Add unit test and remove print for gpu

* remove .h from cmakelist and fix gpu

* Add a Schumann type model too

* change default constructor for loglaw, and fix bug in utau_mean initialization

* add a base amd model that can work without ABL physics to test the log-law wall-model for channel flow

* Remove default constructors

* Add the newly added wall-models to the doc, along with documentation for the base AMD model and Smagorisnky model

* allow the specification of a reference index for the log-law instead of just the first cell

* Change to a while loop, and abort if u_tau hasn't converged

* make clang happy

* formatting change

* fix breaking changes from pr #842

* gpu fix for for the notherm model, needed to use CellSizeArray, and a
GPUArray for grid spacing

---------

Co-authored-by: Marc T. Henry de Frahan <marc.henrydefrahan@nrel.gov>

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

amr-wind/boundary_conditions/wall_models/WallFunction.cpp

@@ -1,4 +1,5 @@
 #include "amr-wind/boundary_conditions/wall_models/WallFunction.H"
+#include "amr-wind/boundary_conditions/wall_models/ShearStressSimple.H"
 #include "amr-wind/utilities/tensor_ops.H"
 #include "amr-wind/utilities/trig_ops.H"
 #include "amr-wind/diffusion/diffusion.H"
@@ -11,37 +12,142 @@
 
 namespace amr_wind {
 
-WallFunction::WallFunction(const CFDSim& sim) : m_sim(sim), m_mesh(m_sim.mesh())
+WallFunction::WallFunction(CFDSim& sim)
+    : m_sim(sim), m_mesh(m_sim.mesh()), m_pa_vel(sim, m_direction)
 {
     {
         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");
+        m_log_law.utau_mean = std::sqrt(std::sqrt(
+            body_force[0] * body_force[0] + body_force[1] * body_force[1]));
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
             std::abs(body_force[2]) < 1e-16,
             "body force in z should be zero for this wall function");
     }
+    {
+        amrex::ParmParse pp("transport");
+        pp.query("viscosity", m_log_law.nu);
+    }
+    {
+        amrex::ParmParse pp("WallFunction");
+        // Reference height for log-law
+        if (pp.contains("log_law_ref_index")) {
+            pp.get("log_law_ref_index", m_log_law.ref_index);
+        }
+        const auto& geom = m_mesh.Geom(0);
+        m_log_law.zref =
+            (geom.ProbLo(m_direction) +
+             (m_log_law.ref_index + 0.5) * geom.CellSize(m_direction));
+    }
 }
 
-VelWallFunc::VelWallFunc(Field& /*unused*/, const WallFunction& wall_func)
+VelWallFunc::VelWallFunc(Field& /*unused*/, 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") {
+    if (m_wall_shear_stress_type == "constant" ||
+        m_wall_shear_stress_type == "log_law" ||
+        m_wall_shear_stress_type == "schumann") {
         amrex::Print() << "Shear Stress model: " << m_wall_shear_stress_type
                        << std::endl;
     } else {
         amrex::Abort("Shear Stress wall model input mistake");
     }
 }
 
+template <typename ShearStressSimple>
+void VelWallFunc::wall_model(
+    Field& velocity, const FieldState rho_state, const ShearStressSimple& tau)
+{
+    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();
+    const auto idx_offset = tau.m_ll.ref_index;
+
+    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);
+        auto& vold_lev = velocity.state(FieldState::Old)(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 vold_arr = vold_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);
+                        const amrex::Real uu =
+                            vold_arr(i, j, k + idx_offset, 0);
+                        const amrex::Real vv =
+                            vold_arr(i, j, k + idx_offset, 1);
+                        const amrex::Real wspd = std::sqrt(uu * uu + vv * vv);
+                        // Dirichlet BC
+                        varr(i, j, k - 1, 2) = 0.0;
+
+                        // Shear stress BC
+                        varr(i, j, k - 1, 0) =
+                            tau.get_shear(uu, wspd) / mu * den(i, j, k);
+                        varr(i, j, k - 1, 1) =
+                            tau.get_shear(vv, wspd) / mu * den(i, j, k);
+                    });
+            }
+
+            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);
+                        const amrex::Real uu =
+                            vold_arr(i, j, k - 1 - idx_offset, 0);
+                        const amrex::Real vv =
+                            vold_arr(i, j, k - 1 - idx_offset, 1);
+                        const amrex::Real wspd = std::sqrt(uu * uu + vv * vv);
+                        // Dirichlet BC
+                        varr(i, j, k, 2) = 0.0;
+
+                        // Shear stress BC
+                        varr(i, j, k, 0) =
+                            -tau.get_shear(uu, wspd) / mu * den(i, j, k);
+                        varr(i, j, k, 1) =
+                            -tau.get_shear(vv, wspd) / mu * den(i, j, k);
+                    });
+            }
+        }
+    }
+}
+
 void VelWallFunc::wall_model(
     Field& velocity, const FieldState rho_state, const amrex::Real utau)
 {
@@ -119,7 +225,24 @@ 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());
+    } else if (m_wall_shear_stress_type == "log_law") {
+        m_wall_func.update_umean();
+        m_wall_func.update_utau_mean();
+        auto tau = SimpleShearLogLaw(m_wall_func.log_law());
+        wall_model(velocity, rho_state, tau);
+    } else if (m_wall_shear_stress_type == "schumann") {
+        m_wall_func.update_umean();
+        auto tau = SimpleShearSchumann(m_wall_func.log_law());
+        wall_model(velocity, rho_state, tau);
     }
 }
 
+void WallFunction::update_umean()
+{
+    m_pa_vel();
+    m_log_law.wspd_mean =
+        m_pa_vel.line_hvelmag_average_interpolated(m_log_law.zref);
+}
+
+void WallFunction::update_utau_mean() { m_log_law.update_utau_mean(); }
 } // namespace amr_wind