From e79c8e156989da4d40c7d0bb7c57c541971b5432 Mon Sep 17 00:00:00 2001
From: Eliot Quon <eliot.quon@nrel.gov>
Date: Mon, 13 Nov 2023 08:55:12 -0700
Subject: [PATCH] Update GABLS1 problem for LES (#1290)

* Init QKE only if state array has that component

* Perturb theta, not rho*theta, if T_0_Pert_Mag is given

* Make pert ref height ~ the height of the divergence-free perturbation layer
---
 Exec/RegTests/GABLS1/prob.cpp | 38 ++++++++++++++++++++++++++++-------
 1 file changed, 31 insertions(+), 7 deletions(-)

diff --git a/Exec/RegTests/GABLS1/prob.cpp b/Exec/RegTests/GABLS1/prob.cpp
index c9bffb30f..314897e03 100644
--- a/Exec/RegTests/GABLS1/prob.cpp
+++ b/Exec/RegTests/GABLS1/prob.cpp
@@ -70,6 +70,28 @@ Problem::init_custom_pert(
     amrex::Array4<amrex::Real const> const& /*mf_v*/,
     const SolverChoice& /*sc*/)
 {
+  if (state.nComp() > RhoQKE_comp) {
+    amrex::Print() << "Initializing QKE" << std::endl;
+  }
+
+  if (parms.pert_ref_height > 0) {
+    if (parms.pert_deltaU != 0.0) {
+      amrex::Print() << "Adding divergence-free x-velocity perturbations" << std::endl;
+    }
+    if (parms.pert_deltaV != 0.0) {
+      amrex::Print() << "Adding divergence-free y-velocity perturbations" << std::endl;
+    }
+    if (parms.T_0_Pert_Mag != 0.0) {
+      amrex::Print() << "Adding random theta perturbations" << std::endl;
+    }
+    if (parms.U_0_Pert_Mag != 0.0) {
+      amrex::Print() << "Adding random x-velocity perturbations" << std::endl;
+    }
+    if (parms.V_0_Pert_Mag != 0.0) {
+      amrex::Print() << "Adding random y-velocity perturbations" << std::endl;
+    }
+  }
+
   ParallelForRNG(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept {
     // Geometry
     const Real* prob_lo = geomdata.ProbLo();
@@ -89,18 +111,20 @@ Problem::init_custom_pert(
     // Add temperature perturbations
     if ((z <= parms.pert_ref_height) && (parms.T_0_Pert_Mag != 0.0)) {
         Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0
-        state(i, j, k, RhoTheta_comp) = (rand_double*2.0 - 1.0)*parms.T_0_Pert_Mag;
+        state(i, j, k, RhoTheta_comp) = r_hse(i,j,k) * (rand_double*2.0 - 1.0)*parms.T_0_Pert_Mag;
     }
 
     // Set scalar = A_0*exp(-10r^2), where r is distance from center of domain
     state(i, j, k, RhoScalar_comp) = parms.A_0 * exp(-10.*r*r);
 
     // Set an initial value for QKE = 2*e
-    if (z < 250) {
-        // Following Cuxart et al. 2006
-        state(i, j, k, RhoQKE_comp) = r_hse(i,j,k) * 0.8 * std::pow(1 - z/250, 3);
-    } else {
-        state(i, j, k, RhoQKE_comp) = r_hse(i,j,k) * parms.QKE_0;
+    if (state.nComp() > RhoQKE_comp) {
+        if (z < 250) {
+            // Following Cuxart et al. 2006
+            state(i, j, k, RhoQKE_comp) = r_hse(i,j,k) * 0.8 * std::pow(1 - z/250, 3);
+        } else {
+            state(i, j, k, RhoQKE_comp) = r_hse(i,j,k) * parms.QKE_0;
+        }
     }
 
 #if defined(ERF_USE_MOISTURE)
@@ -130,7 +154,7 @@ Problem::init_custom_pert(
     if (parms.pert_deltaU != 0.0)
     {
         const amrex::Real yl = y - prob_lo[1];
-        const amrex::Real zl = z / parms.pert_ref_height;
+        const amrex::Real zl = z / (parms.pert_ref_height/2);
         const amrex::Real damp = std::exp(-0.5 * zl * zl);
         x_vel(i, j, k) += parms.ufac * damp * z * std::cos(parms.aval * yl);
     }