From f329979aa2bcabdcc1bf27223a59d0b2bb75f4f3 Mon Sep 17 00:00:00 2001
From: Marc Henry de Frahan <marc.henrydefrahan@nrel.gov>
Date: Mon, 25 Nov 2024 19:13:13 -0700
Subject: [PATCH] beta in turb models

---
 amr-wind/turbulence/LES/AMD.cpp       | 4 ----
 amr-wind/turbulence/LES/OneEqKsgs.cpp | 6 +++---
 amr-wind/turbulence/RANS/KLAxell.cpp  | 9 +++++----
 3 files changed, 8 insertions(+), 11 deletions(-)

diff --git a/amr-wind/turbulence/LES/AMD.cpp b/amr-wind/turbulence/LES/AMD.cpp
index 9999d4b5e1..b2be6229a4 100644
--- a/amr-wind/turbulence/LES/AMD.cpp
+++ b/amr-wind/turbulence/LES/AMD.cpp
@@ -25,10 +25,6 @@ AMD<Transport>::AMD(CFDSim& sim)
 {
     auto& phy_mgr = this->m_sim.physics_manager();
     if (phy_mgr.contains("ABL")) {
-        {
-            amrex::ParmParse pp("ABL");
-            pp.get("reference_temperature", m_ref_theta);
-        }
         {
             amrex::ParmParse pp("incflo");
             pp.queryarr("gravity", m_gravity);
diff --git a/amr-wind/turbulence/LES/OneEqKsgs.cpp b/amr-wind/turbulence/LES/OneEqKsgs.cpp
index da1c6653e4..393186988c 100644
--- a/amr-wind/turbulence/LES/OneEqKsgs.cpp
+++ b/amr-wind/turbulence/LES/OneEqKsgs.cpp
@@ -43,7 +43,6 @@ OneEqKsgsM84<Transport>::OneEqKsgsM84(CFDSim& sim)
 
     {
         amrex::ParmParse pp("ABL");
-        pp.get("reference_temperature", m_ref_theta);
         pp.query("enable_hybrid_rl_mode", m_hybrid_rl);
     }
 
@@ -101,7 +100,7 @@ void OneEqKsgsM84<Transport>::update_turbulent_viscosity(
 
     const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> gravity{
         m_gravity[0], m_gravity[1], m_gravity[2]};
-    const amrex::Real beta = 1.0 / m_ref_theta;
+    const auto beta = (this->m_transport).beta();
 
     auto& mu_turb = this->mu_turb();
     const amrex::Real Ce = this->m_Ce;
@@ -127,6 +126,7 @@ void OneEqKsgsM84<Transport>::update_turbulent_viscosity(
             const auto& tke_arr = (*this->m_tke)(lev).array(mfi);
             const auto& buoy_prod_arr = (this->m_buoy_prod)(lev).array(mfi);
             const auto& shear_prod_arr = (this->m_shear_prod)(lev).array(mfi);
+            const auto& beta_arr = (*beta)(lev).array(mfi);
 
             amrex::ParallelFor(
                 bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
@@ -134,7 +134,7 @@ void OneEqKsgsM84<Transport>::update_turbulent_viscosity(
                         -(gradT_arr(i, j, k, 0) * gravity[0] +
                           gradT_arr(i, j, k, 1) * gravity[1] +
                           gradT_arr(i, j, k, 2) * gravity[2]) *
-                        beta;
+                      beta_arr(i,j,k);
                     if (stratification > 1e-10) {
                         tlscale_arr(i, j, k) = amrex::min<amrex::Real>(
                             ds, 0.76 * std::sqrt(
diff --git a/amr-wind/turbulence/RANS/KLAxell.cpp b/amr-wind/turbulence/RANS/KLAxell.cpp
index edf25a8ed2..c74a72c2ec 100644
--- a/amr-wind/turbulence/RANS/KLAxell.cpp
+++ b/amr-wind/turbulence/RANS/KLAxell.cpp
@@ -31,7 +31,6 @@ KLAxell<Transport>::KLAxell(CFDSim& sim)
     }
     {
         amrex::ParmParse pp("ABL");
-        pp.get("reference_temperature", m_ref_theta);
         pp.get("surface_temp_flux", m_surf_flux);
         pp.query("length_scale_switch", m_lengthscale_switch);
     }
@@ -84,7 +83,7 @@ void KLAxell<Transport>::update_turbulent_viscosity(
 
     const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> gravity{
         m_gravity[0], m_gravity[1], m_gravity[2]};
-    const amrex::Real beta = 1.0 / m_ref_theta;
+    const auto beta = (this->m_transport).beta();
     const amrex::Real Cmu = m_Cmu;
     const amrex::Real Cb_stable = m_Cb_stable;
     const amrex::Real Cb_unstable = m_Cb_unstable;
@@ -113,6 +112,8 @@ void KLAxell<Transport>::update_turbulent_viscosity(
             const auto& tke_arr = (*this->m_tke)(lev).array(mfi);
             const auto& buoy_prod_arr = (this->m_buoy_prod)(lev).array(mfi);
             const auto& shear_prod_arr = (this->m_shear_prod)(lev).array(mfi);
+            const auto& beta_arr = (*beta)(lev).array(mfi);
+
             //! Add terrain components
             const bool has_terrain =
                 this->m_sim.repo().int_field_exists("terrain_blank");
@@ -130,7 +131,7 @@ void KLAxell<Transport>::update_turbulent_viscosity(
                             -(gradT_arr(i, j, k, 0) * gravity[0] +
                               gradT_arr(i, j, k, 1) * gravity[1] +
                               gradT_arr(i, j, k, 2) * gravity[2]) *
-                            beta;
+                          beta_arr(i,j,k);
                         const amrex::Real z = std::max(
                             problo[2] + (k + 0.5) * dz - ht_arr(i, j, k),
                             0.5 * dz);
@@ -201,7 +202,7 @@ void KLAxell<Transport>::update_turbulent_viscosity(
                             -(gradT_arr(i, j, k, 0) * gravity[0] +
                               gradT_arr(i, j, k, 1) * gravity[1] +
                               gradT_arr(i, j, k, 2) * gravity[2]) *
-                            beta;
+                          beta_arr(i,j,k);
                         const amrex::Real z = problo[2] + (k + 0.5) * dz;
                         const amrex::Real lscale_s =
                             (lambda * kappa * z) / (lambda + kappa * z);