Fingers crossed

When calculating the diffusion source, calculate the strain rate everywhere
before resizing the box for the stress calc. This provides the correct strain
rate on the surface, without making any implicit assumption about dw/dx, dw/dy.

Source/Diffusion/EddyViscosity.H

@@ -42,28 +42,35 @@ ComputeSmnSmn (int& i, int& j, int& k,
                                 + tau12(i+1, j  , k  ) + tau12(i+1, j+1, k  ) );
 
     // NOTES:
-    // - We do not use the strains lying on the bottom boundary with MOST. If
-    //   ERF_EXPLICIT_MOST_STRESS is not used, then these values are corrupted
-    //   with the reflect odd BC used to fill the ghost cells.
-    // - If ERF_EXPLICIT_MOST_STRESS _is_ used, then the assumption in
-    //   MOSTStress.H means that dU/dz at klo and klo+1 are equal.
-    // - In both cases, the strains at klo and klo+1 differ by dW/dx at klo+1.
-    //   This needs to be removed at the surface to be consistent with surface
-    //   conditions.
+    // - If ERF_EXPLICIT_MOST_STRESS is not used, then we do not 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.
+    // - Neglecting the strains on the bottom boundary and using a one-
+    //   sided average implies that the strains are equal at klo and klo+1.
+    //   du/dz and dv/dz can be assumed to be equal through the first cell.
+    //   However, whereas dw/dx = dw/dy = 0 on the surface, this is not true at
+    //   klo+1. Therefore, the strains are not equal in general.
+    // - If ERF_EXPLICIT_MOST_STRESS _is_ used, then strains may be correctly
+    //   evaluated at the surface and no assumptions are needed.
     amrex::Real s13bar;
+#ifndef ERF_EXPLICIT_MOST_STRESS
     if (use_most && k==klo) {
         s13bar = 0.5  * ( tau13(i  , j  , k+1) + tau13(i+1, j  , k+1) );
     }
     else
+#endif
     {
         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;
+#ifndef ERF_EXPLICIT_MOST_STRESS
     if (use_most && k==klo) {
         s23bar = 0.5  * ( tau23(i  , j  , k+1) + tau23(i  , j+1, k+1) );
     }
     else
+#endif
     {
         s23bar = 0.25 * ( tau23(i  , j  , k  ) + tau23(i  , j  , k+1)
                         + tau23(i  , j+1, k  ) + tau23(i  , j+1, k+1) );

Source/TimeIntegration/ERF_advance_dycore.cpp

@@ -118,14 +118,6 @@ void ERF::advance_dycore(int level,
             Box tbxxz = mfi.tilebox(IntVect(1,0,1),IntVect(1,1,0));
             Box tbxyz = mfi.tilebox(IntVect(0,1,1),IntVect(1,1,0));
 
-#ifdef ERF_EXPLICIT_MOST_STRESS
-            if (use_most) {
-                // Don't overwrite modeled total stress value at boundary
-                tbxxz.setSmall(2,1);
-                tbxyz.setSmall(2,1);
-            }
-#endif
-
             const Array4<const Real> & u = xvel_old.array(mfi);
             const Array4<const Real> & v = yvel_old.array(mfi);
             const Array4<const Real> & w = zvel_old.array(mfi);

Source/TimeIntegration/ERF_slow_rhs_pre.cpp

@@ -238,14 +238,6 @@ void erf_slow_rhs_pre (int level, int finest_level,
             Box tbxxz = mfi.tilebox(IntVect(1,0,1));
             Box tbxyz = mfi.tilebox(IntVect(0,1,1));
 
-#ifdef ERF_EXPLICIT_MOST_STRESS
-            if (use_most) {
-                // Don't overwrite modeled total stress value at boundary
-                tbxxz.setSmall(2,1);
-                tbxyz.setSmall(2,1);
-            }
-#endif
-
             // We need a halo cell for terrain
              bxcc.grow(IntVect(1,1,0));
             tbxxy.grow(IntVect(1,1,0));
@@ -336,6 +328,18 @@ void erf_slow_rhs_pre (int level, int finest_level,
                     });
                 }
 
+#ifdef ERF_EXPLICIT_MOST_STRESS
+                // We've updated the strains at all locations including the
+                // surface. This is required to get the correct strain-rate
+                // magnitude. Now, update the stress everywhere but the surface
+                // to retain the values set by MOST.
+                if (use_most) {
+                    // Don't overwrite modeled total stress value at boundary
+                    tbxxz.setSmall(2,1);
+                    tbxyz.setSmall(2,1);
+                }
+#endif
+
                 //-----------------------------------------
                 // Stress tensor compute terrain
                 //-----------------------------------------
@@ -434,6 +438,18 @@ void erf_slow_rhs_pre (int level, int finest_level,
                     });
                 }
 
+#ifdef ERF_EXPLICIT_MOST_STRESS
+                // We've updated the strains at all locations including the
+                // surface. This is required to get the correct strain-rate
+                // magnitude. Now, update the stress everywhere but the surface
+                // to retain the values set by MOST.
+                if (use_most) {
+                    // Don't overwrite modeled total stress value at boundary
+                    tbxxz.setSmall(2,1);
+                    tbxyz.setSmall(2,1);
+                }
+#endif
+
                 //-----------------------------------------
                 // Stress tensor compute no terrain
                 //-----------------------------------------