Merge pull request #1348 from baperry2/pbl-terrain

Add terrain/stretched grid capability for PBL models

Source/Diffusion/ComputeTurbulentViscosity.cpp

@@ -16,7 +16,8 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
                               const TurbChoice& turbChoice,
                               std::unique_ptr<ABLMost>& most,
                               const amrex::BCRec* bc_ptr,
-                              bool /*vert_only*/);
+                              bool /*vert_only*/,
+                              const std::unique_ptr<amrex::MultiFab>& z_phys_nd);
 
 /**
  * Function for computing the turbulent viscosity with LES.
@@ -368,6 +369,7 @@ void ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::Multi
                                 amrex::MultiFab& Hfx1, amrex::MultiFab& Hfx2, amrex::MultiFab& Hfx3, amrex::MultiFab& Diss,
                                 const amrex::Geometry& geom,
                                 const amrex::MultiFab& mapfac_u, const amrex::MultiFab& mapfac_v,
+                                const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
                                 const TurbChoice& turbChoice, const Real const_grav,
                                 std::unique_ptr<ABLMost>& most,
                                 const amrex::BCRec* bc_ptr,
@@ -407,6 +409,6 @@ void ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::Multi
     if (turbChoice.pbl_type != PBLType::None) {
         // NOTE: state_new is passed in for Cons_old (due to ptr swap in advance)
         ComputeTurbulentViscosityPBL(xvel, yvel, cons_in, eddyViscosity,
-                                     geom, turbChoice, most, bc_ptr, vert_only);
+                                     geom, turbChoice, most, bc_ptr, vert_only, z_phys_nd);
     }
 }

Source/Diffusion/DiffusionSrcForState_N.cpp

@@ -1,6 +1,6 @@
 #include <Diffusion.H>
 #include <EddyViscosity.H>
-#include <ComputeQKESourceTerm.H>
+#include <PBLModels.H>
 
 using namespace amrex;
 

Source/Diffusion/DiffusionSrcForState_T.cpp

@@ -1,7 +1,7 @@
 #include <Diffusion.H>
 #include <EddyViscosity.H>
 #include <TerrainMetrics.H>
-#include <ComputeQKESourceTerm.H>
+#include <PBLModels.H>
 
 using namespace amrex;
 
@@ -565,11 +565,13 @@ DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain,
         bool v_ext_dir_on_zhi = ( (bc_ptr[BCVars::yvel_bc].lo(5) == ERFBCType::ext_dir) );
         amrex::ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
+            const Real met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd);
             cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim,
                                                                   mu_turb,dxInv,domain,pbl_B1_l,tm_arr(i,j,0),
                                                                   c_ext_dir_on_zlo, c_ext_dir_on_zhi,
                                                                   u_ext_dir_on_zlo, u_ext_dir_on_zhi,
-                                                                  v_ext_dir_on_zlo, v_ext_dir_on_zhi);
+                                                                  v_ext_dir_on_zlo, v_ext_dir_on_zhi,
+                                                                  met_h_zeta);
         });
     }
 }

Source/Diffusion/EddyViscosity.H

@@ -16,6 +16,7 @@ ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::MultiFab&
                            amrex::MultiFab& Hfx1, amrex::MultiFab& Hfx2, amrex::MultiFab& Hfx3, amrex::MultiFab& Diss,
                            const amrex::Geometry& geom,
                            const amrex::MultiFab& mapfac_u, const amrex::MultiFab& mapfac_v,
+                           const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
                            const TurbChoice& turbChoice, const amrex::Real const_grav,
                            std::unique_ptr<ABLMost>& most,
                            const amrex::BCRec* bc_ptr,

Source/Diffusion/Make.package

@@ -3,7 +3,7 @@ CEXE_sources += DiffusionSrcForMom_T.cpp
 
 CEXE_sources += DiffusionSrcForState_N.cpp
 CEXE_sources += DiffusionSrcForState_T.cpp
-    
+
 CEXE_sources += ComputeStress_N.cpp
 CEXE_sources += ComputeStress_T.cpp
 
@@ -13,8 +13,8 @@ CEXE_sources += ComputeStrain_T.cpp
 CEXE_sources += PBLModels.cpp
 CEXE_sources += NumericalDiffusion.cpp
 CEXE_sources += ComputeTurbulentViscosity.cpp
-        
+
 CEXE_headers += Diffusion.H
 CEXE_headers += EddyViscosity.H
 CEXE_headers += NumericalDiffusion.H
-CEXE_headers += ComputeQKESourceTerm.H
+CEXE_headers += PBLModels.H

Source/Diffusion/ComputeQKESourceTerm.H → Source/Diffusion/PBLModels.H

@@ -2,44 +2,32 @@
 #define _PBLMODELS_H_
 
 /**
- * Function for computing the QKE source terms.
+ * Function for computing vertical derivatives for use in PBL model
  *
  * @param[in] u velocity in x-dir
  * @param[in] v velocity in y-dir
  * @param[in] cell_data conserved cell center vars
- * @param[in] cell_prim primitive cell center vars
- * @param[in] K_turb turbulent viscosity
- * @param[in] cellSizeInv inverse cell size array
- * @param[in] domain box of the whole domain
- * @param[in] pbl_B1_l a parameter
- * @param[in] theta_mean average theta
  */
 AMREX_GPU_DEVICE
 AMREX_FORCE_INLINE
-amrex::Real
-ComputeQKESourceTerms (int i, int j, int k,
-                       const amrex::Array4<const amrex::Real>& uvel,
-                       const amrex::Array4<const amrex::Real>& vvel,
-                       const amrex::Array4<const amrex::Real>& cell_data,
-                       const amrex::Array4<const amrex::Real>& cell_prim,
-                       const amrex::Array4<const amrex::Real>& K_turb,
-                       const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
-                       const amrex::Box& domain,
-                       amrex::Real pbl_B1_l,
-                       const amrex::Real theta_mean,
-                       bool c_ext_dir_on_zlo,
-                       bool c_ext_dir_on_zhi,
-                       bool u_ext_dir_on_zlo,
-                       bool u_ext_dir_on_zhi,
-                       bool v_ext_dir_on_zlo,
-                       bool v_ext_dir_on_zhi)
+void
+ComputeVerticalDerivativesPBL(int i, int j, int k,
+                              const amrex::Array4<const amrex::Real>& uvel,
+                              const amrex::Array4<const amrex::Real>& vvel,
+                              const amrex::Array4<const amrex::Real>& cell_data,
+                              const int izmin,
+                              const int izmax,
+                              const amrex::Real& dz_inv,
+                              const bool c_ext_dir_on_zlo,
+                              const bool c_ext_dir_on_zhi,
+                              const bool u_ext_dir_on_zlo,
+                              const bool u_ext_dir_on_zhi,
+                              const bool v_ext_dir_on_zlo,
+                              const bool v_ext_dir_on_zhi,
+                              amrex::Real& dthetadz,
+                              amrex::Real& dudz,
+                              amrex::Real& dvdz)
 {
-    // Compute some relevant derivatives
-    amrex::Real dz_inv = cellSizeInv[2];
-    int izmin = domain.smallEnd(2);
-    int izmax = domain.bigEnd(2);
-    amrex::Real dthetadz, dudz, dvdz;
-    amrex::Real source_term = 0.0;
     if ( k==izmax && c_ext_dir_on_zhi ) {
         dthetadz = (1.0/3.0)*(-cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp)
                               - 3.0 * cell_data(i,j,k  ,RhoTheta_comp)/cell_data(i,j,k  ,Rho_comp)
@@ -72,6 +60,56 @@ ComputeQKESourceTerms (int i, int j, int k,
     } else {
         dvdz = 0.25*( vvel(i,j,k+1) - vvel(i,j,k-1) + vvel(i,j+1,k+1) - vvel(i,j+1,k-1) )*dz_inv;
     }
+}
+
+/**
+ * Function for computing the QKE source terms.
+ *
+ * @param[in] u velocity in x-dir
+ * @param[in] v velocity in y-dir
+ * @param[in] cell_data conserved cell center vars
+ * @param[in] cell_prim primitive cell center vars
+ * @param[in] K_turb turbulent viscosity
+ * @param[in] cellSizeInv inverse cell size array
+ * @param[in] domain box of the whole domain
+ * @param[in] pbl_B1_l a parameter
+ * @param[in] theta_mean average theta
+ */
+AMREX_GPU_DEVICE
+AMREX_FORCE_INLINE
+amrex::Real
+ComputeQKESourceTerms (int i, int j, int k,
+                       const amrex::Array4<const amrex::Real>& uvel,
+                       const amrex::Array4<const amrex::Real>& vvel,
+                       const amrex::Array4<const amrex::Real>& cell_data,
+                       const amrex::Array4<const amrex::Real>& cell_prim,
+                       const amrex::Array4<const amrex::Real>& K_turb,
+                       const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
+                       const amrex::Box& domain,
+                       amrex::Real pbl_B1_l,
+                       const amrex::Real theta_mean,
+                       bool c_ext_dir_on_zlo,
+                       bool c_ext_dir_on_zhi,
+                       bool u_ext_dir_on_zlo,
+                       bool u_ext_dir_on_zhi,
+                       bool v_ext_dir_on_zlo,
+                       bool v_ext_dir_on_zhi,
+                       const amrex::Real met_h_zeta=1.0)
+{
+    // Compute some relevant derivatives
+    amrex::Real dthetadz, dudz, dvdz;
+    amrex::Real source_term = 0.0;
+
+    amrex::Real dz_inv = cellSizeInv[2];
+    int izmin = domain.smallEnd(2);
+    int izmax = domain.bigEnd(2);
+
+    ComputeVerticalDerivativesPBL(i, j, k,
+                                  uvel, vvel, cell_data, izmin, izmax, dz_inv/met_h_zeta,
+                                  c_ext_dir_on_zlo, c_ext_dir_on_zhi,
+                                  u_ext_dir_on_zlo, u_ext_dir_on_zhi,
+                                  v_ext_dir_on_zlo, v_ext_dir_on_zhi,
+                                  dthetadz, dudz, dvdz);
 
     // Production (We store mu_turb, which is 0.5*K_turb)
     source_term += 4.0*K_turb(i,j,k,EddyDiff::Mom_v) * (dudz*dudz + dvdz*dvdz);

Source/Diffusion/PBLModels.cpp

@@ -3,6 +3,7 @@
 #include "Diffusion.H"
 #include "ERF_Constants.H"
 #include "TurbStruct.H"
+#include "PBLModels.H"
 
 /**
  * Function to compute turbulent viscosity with PBL.
@@ -24,8 +25,11 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
                               const TurbChoice& turbChoice,
                               std::unique_ptr<ABLMost>& most,
                               const amrex::BCRec* bc_ptr,
-                              bool /*vert_only*/)
+                              bool /*vert_only*/,
+                              const std::unique_ptr<amrex::MultiFab>& z_phys_nd)
 {
+    const bool use_terrain = (z_phys_nd != nullptr);
+
     // MYNN Level 2.5 PBL Model
     if (turbChoice.pbl_type == PBLType::MYNN25) {
 
@@ -79,22 +83,45 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
             const amrex::Array4<amrex::Real> qvel = qturb.array();
             const amrex::Array4<amrex::Real> qvel_old = qturb_old.array();
 
-            amrex::ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), bx,
-                               [=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept
-            {
-                qvel(i,j,k)     = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp));
-                qvel_old(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp) + eps);
-                AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > 0.0, "QKE must have a positive value");
-                AMREX_ASSERT_WITH_MESSAGE(qvel_old(i,j,k) > 0.0, "Old QKE must have a positive value");
-
-                const amrex::Real Zval = gdata.ProbLo(2) + (k + 0.5)*gdata.CellSize(2);
-                if (sbx.contains(i,j,k)) {
-                    amrex::Gpu::deviceReduceSum(&qint(i,j,0,0), Zval*qvel(i,j,k), handler);
-                    amrex::Gpu::deviceReduceSum(&qint(i,j,0,1), qvel(i,j,k), handler);
-                }
-            });
+            // vertical integrals to compute lengthscale
+            if (use_terrain) {
+                const amrex::Array4<amrex::Real const> &z_nd = z_phys_nd->array(mfi);
+                const auto invCellSize = geom.InvCellSizeArray();
+                amrex::ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), bx,
+                                   [=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept
+                {
+                    qvel(i,j,k)     = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp));
+                    qvel_old(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp) + eps);
+                    AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > 0.0, "QKE must have a positive value");
+                    AMREX_ASSERT_WITH_MESSAGE(qvel_old(i,j,k) > 0.0, "Old QKE must have a positive value");
+
+                    if (sbx.contains(i,j,k)) {
+                        const amrex::Real Zval = Compute_Zrel_AtCellCenter(i,j,k,z_nd);
+                        const amrex::Real dz = Compute_h_zeta_AtCellCenter(i,j,k,invCellSize,z_nd);
+                        amrex::Gpu::deviceReduceSum(&qint(i,j,0,0), Zval*qvel(i,j,k)*dz, handler);
+                        amrex::Gpu::deviceReduceSum(&qint(i,j,0,1), qvel(i,j,k)*dz, handler);
+                    }
+                });
+            } else {
+                amrex::ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), bx,
+                                   [=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept
+                {
+                    qvel(i,j,k)     = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp));
+                    qvel_old(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp) + eps);
+                    AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > 0.0, "QKE must have a positive value");
+                    AMREX_ASSERT_WITH_MESSAGE(qvel_old(i,j,k) > 0.0, "Old QKE must have a positive value");
+
+                    // Not multiplying by dz: its constant and would fall out when we divide qint0/qint1 anyway
+                    if (sbx.contains(i,j,k)) {
+                        const amrex::Real Zval = gdata.ProbLo(2) + (k + 0.5)*gdata.CellSize(2);
+                        amrex::Gpu::deviceReduceSum(&qint(i,j,0,0), Zval*qvel(i,j,k), handler);
+                        amrex::Gpu::deviceReduceSum(&qint(i,j,0,1), qvel(i,j,k), handler);
+                    }
+                });
+            }
 
             amrex::Real dz_inv = geom.InvCellSize(2);
+            const auto& dxInv = geom.InvCellSizeArray();
             int izmin = geom.Domain().smallEnd(2);
             int izmax = geom.Domain().bigEnd(2);
 
@@ -110,43 +137,23 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
             const auto& u_star_arr = u_star_mf->array(mfi);
             const auto& t_star_arr = t_star_mf->array(mfi);
 
+            const amrex::Array4<amrex::Real const> z_nd;
+            if (use_terrain) {
+                const auto& z_nd = z_phys_nd->array(mfi);
+            }
+
             amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 // Compute some partial derivatives that we will need (second order)
                 // U and V derivatives are interpolated to account for staggered grid
+                const amrex::Real met_h_zeta = use_terrain ? Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd) : 1.0;
                 amrex::Real dthetadz, dudz, dvdz;
-                if ( k==izmax && c_ext_dir_on_zhi ) {
-                    dthetadz = (1.0/3.0)*(-cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp)
-                                          - 3.0 * cell_data(i,j,k  ,RhoTheta_comp)/cell_data(i,j,k  ,Rho_comp)
-                                          + 4.0 * cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) )*dz_inv;
-                } else if ( k==izmin && c_ext_dir_on_zlo ) {
-                    dthetadz = (1.0/3.0)*( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp)
-                                           + 3.0 * cell_data(i,j,k  ,RhoTheta_comp)/cell_data(i,j,k  ,Rho_comp)
-                                           - 4.0 * cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) )*dz_inv;
-                } else {
-                    dthetadz = 0.5*(cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp)
-                                    - cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp))*dz_inv;
-                }
-
-                if ( k==izmax && u_ext_dir_on_zhi ) {
-                    dudz = (1.0/6.0)*( (-uvel(i  ,j,k-1) - 3.0 * uvel(i  ,j,k  ) + 4.0 * uvel(i  ,j,k+1))
-                                       + (-uvel(i+1,j,k-1) - 3.0 * uvel(i+1,j,k  ) + 4.0 * uvel(i+1,j,k+1)) )*dz_inv;
-                } else if ( k==izmin && u_ext_dir_on_zlo ) {
-                    dudz = (1.0/6.0)*( (uvel(i  ,j,k+1) + 3.0 * uvel(i  ,j,k  ) - 4.0 * uvel(i  ,j,k-1))
-                                       + (uvel(i+1,j,k+1) + 3.0 * uvel(i+1,j,k  ) - 4.0 * uvel(i+1,j,k-1)) )*dz_inv;
-                } else {
-                    dudz = 0.25*(uvel(i,j,k+1) - uvel(i,j,k-1) + uvel(i+1,j,k+1) - uvel(i+1,j,k-1))*dz_inv;
-                }
-
-                if ( k==izmax && v_ext_dir_on_zhi ) {
-                    dvdz = (1.0/6.0)*( (-vvel(i,j  ,k-1) - 3.0 * vvel(i,j  ,k  ) + 4.0 * vvel(i,j  ,k+1))
-                                       + (-vvel(i,j+1,k-1) - 3.0 * vvel(i,j+1,k  ) + 4.0 * vvel(i,j+1,k+1)) )*dz_inv;
-                } else if ( k==izmin && v_ext_dir_on_zlo ) {
-                    dvdz = (1.0/6.0)*( (vvel(i,j  ,k+1) + 3.0 * vvel(i,j  ,k  ) - 4.0 * vvel(i,j  ,k-1))
-                                       + (vvel(i,j+1,k+1) + 3.0 * vvel(i,j+1,k  ) - 4.0 * vvel(i,j+1,k-1)) )*dz_inv;
-                } else {
-                    dvdz = 0.25*(vvel(i,j,k+1) - vvel(i,j,k-1) + vvel(i,j+1,k+1) - vvel(i,j+1,k-1))*dz_inv;
-                }
+                ComputeVerticalDerivativesPBL(i, j, k,
+                                              uvel, vvel, cell_data, izmin, izmax, dz_inv/met_h_zeta,
+                                              c_ext_dir_on_zlo, c_ext_dir_on_zhi,
+                                              u_ext_dir_on_zlo, u_ext_dir_on_zhi,
+                                              v_ext_dir_on_zlo, v_ext_dir_on_zhi,
+                                              dthetadz, dudz, dvdz);
 
                 // Spatially varying MOST
                 amrex::Real surface_heat_flux = -u_star_arr(i,j,0) * t_star_arr(i,j,0);

Source/TimeIntegration/ERF_advance_dycore.cpp

@@ -220,6 +220,7 @@ void ERF::advance_dycore(int level,
                                   state_old[IntVar::cons],
                                   *eddyDiffs, *Hfx1, *Hfx2, *Hfx3, *Diss, // to be updated
                                   fine_geom, *mapfac_u[level], *mapfac_v[level],
+                                  z_phys_nd[level],
                                   tc, solverChoice.gravity, m_most, bc_ptr_d);
     }
 

Source/Utils/TerrainMetrics.H

@@ -337,6 +337,26 @@ Compute_h_eta_AtEdgeCenterI (const int &i, const int &j, const int &k,
     return met_h_eta;
 }
 
+// Relative height above terrain surface at cell center from z_nd (nodal absolute height)
+AMREX_GPU_DEVICE
+AMREX_FORCE_INLINE
+amrex::Real
+Compute_Zrel_AtCellCenter (const int &i, const int &j, const int &k,
+                           const amrex::Array4<const amrex::Real>& z_nd)
+{
+    const amrex::Real z_cc = 0.125*(  z_nd(i  ,j  ,k  ) + z_nd(i  ,j  ,k+1) +
+                                      + z_nd(i+1,j  ,k  ) + z_nd(i  ,j  ,k+1)
+                                      + z_nd(i  ,j+1,k  ) + z_nd(i  ,j+1,k+1)
+                                      + z_nd(i+1,j+1,k  ) + z_nd(i  ,j+1,k+1));
+
+    // Note: we assume the z_nd array spans from the bottom to top of the domain
+    // i.e. no domain decomposition across processors in vertical direction
+    const amrex::Real z0_cc = 0.25*(  z_nd(i  ,j  ,0) + z_nd(i  ,j+1,0)
+                                      + z_nd(i+1,j  ,0) + z_nd(i+1,j+1,0));
+
+    return (z_cc - z0_cc);
+}
+
 /**
  * Define omega given u,v and w
  */