Use one-sided diff for smag with MOST

Source/Diffusion/ComputeTurbulentViscosity.cpp

@@ -47,11 +47,12 @@ void ComputeTurbulentViscosityLES (const amrex::MultiFab& Tau11, const amrex::Mu
                                    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)
+                                   const TurbChoice& turbChoice, const Real const_grav, std::unique_ptr<ABLMost>& most)
 {
     const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> cellSizeInv = geom.InvCellSizeArray();
     const Box& domain = geom.Domain();
-
+    const int& klo    = domain.smallEnd(2);
+    const bool use_most = (most != nullptr);
     const bool use_terrain = (z_phys_nd != nullptr);
 
     // SMAGORINSKY: Fill Kturb for momentum in horizontal and vertical
@@ -86,7 +87,7 @@ void ComputeTurbulentViscosityLES (const amrex::MultiFab& Tau11, const amrex::Mu
 
           ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
           {
-              Real SmnSmn = ComputeSmnSmn(i,j,k,tau11,tau22,tau33,tau12,tau13,tau23);
+              Real SmnSmn = ComputeSmnSmn(i,j,k,tau11,tau22,tau33,tau12,tau13,tau23,klo,use_most);
               Real dxInv = cellSizeInv[0];
               Real dyInv = cellSizeInv[1];
               Real dzInv = cellSizeInv[2];
@@ -426,7 +427,7 @@ void ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::Multi
                                      Hfx1, Hfx2, Hfx3, Diss,
                                      geom, mapfac_u, mapfac_v,
                                      z_phys_nd,
-                                     turbChoice, const_grav);
+                                     turbChoice, const_grav, most);
     }
 
     if (turbChoice.pbl_type != PBLType::None) {

Source/Diffusion/EddyViscosity.H

@@ -31,17 +31,33 @@ ComputeSmnSmn (int& i, int& j, int& k,
                const amrex::Array4<amrex::Real const>& tau33,
                const amrex::Array4<amrex::Real const>& tau12,
                const amrex::Array4<amrex::Real const>& tau13,
-               const amrex::Array4<amrex::Real const>& tau23)
+               const amrex::Array4<amrex::Real const>& tau23,
+               const int& klo,
+               const bool& use_most)
 {
     amrex::Real s11bar = tau11(i,j,k);
     amrex::Real s22bar = tau22(i,j,k);
     amrex::Real s33bar = tau33(i,j,k);
     amrex::Real s12bar = 0.25 * ( tau12(i  , j  , k  ) + tau12(i  , j+1, k  )
                                 + tau12(i+1, j  , k  ) + tau12(i+1, j+1, k  ) );
-    amrex::Real s13bar = 0.25 * ( tau13(i  , j  , k  ) + tau13(i  , j  , k+1)
-                                + tau13(i+1, j  , k  ) + tau13(i+1, j  , k+1) );
-    amrex::Real s23bar = 0.25 * ( tau23(i  , j  , k  ) + tau23(i  , j  , k+1)
-                                + tau23(i  , j+1, k  ) + tau23(i  , j+1, k+1) );
+
+    // NOTE: We cannot use the strains lying on the bottom boundary with MOST.
+    //       These values are corrupted with the reflect odd BC used to fill
+    //       the ghost cells. We use one-sided for this case.
+    amrex::Real s13bar;
+    if (use_most && k==klo) {
+        s13bar = 0.5  * ( tau13(i  , j  , k+1) + tau13(i+1, j  , k+1) );
+    } else {
+        s13bar = 0.25 * ( tau13(i  , j  , k  ) + tau13(i  , j  , k+1)
+                        + tau13(i+1, j  , k  ) + tau13(i+1, j  , k+1) );
+    }
+    amrex::Real s23bar;
+    if (use_most && k==klo) {
+        s23bar = 0.5  * ( tau23(i  , j  , k+1) + tau23(i  , j+1, k+1) );
+    } else {
+        s23bar = 0.25 * ( tau23(i  , j  , k  ) + tau23(i  , j  , k+1)
+                        + tau23(i  , j+1, k  ) + tau23(i  , j+1, k+1) );
+    }
     amrex::Real SmnSmn = s11bar*s11bar + s22bar*s22bar + s33bar*s33bar
                        + 2.0*s12bar*s12bar + 2.0*s13bar*s13bar + 2.0*s23bar*s23bar;
 

Source/TimeIntegration/ERF_slow_rhs_pre.cpp

@@ -140,12 +140,14 @@ void erf_slow_rhs_pre (int level, int finest_level,
                                     tc.pbl_type == PBLType::YSU );
 
     const bool use_moisture = (solverChoice.moisture_type != MoistureType::None);
+    const bool use_most     = (most != nullptr);
 
     const amrex::BCRec* bc_ptr   = domain_bcs_type_d.data();
     const amrex::BCRec* bc_ptr_h = domain_bcs_type.data();
 
     const Box& domain = geom.Domain();
-    const int domhi_z = domain.bigEnd()[2];
+    const int domhi_z = domain.bigEnd(2);
+    const int domlo_z = domain.smallEnd(2);
 
     const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
     const Real* dx = geom.CellSize();
@@ -319,7 +321,7 @@ void erf_slow_rhs_pre (int level, int finest_level,
                     SmnSmn_a = SmnSmn->array(mfi);
                     ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                     {
-                        SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23);
+                        SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23,domlo_z,use_most);
                     });
                 }
 
@@ -417,7 +419,7 @@ void erf_slow_rhs_pre (int level, int finest_level,
                     SmnSmn_a = SmnSmn->array(mfi);
                     ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                     {
-                        SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23);
+                        SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23,domlo_z,use_most);
                     });
                 }