Fix non-Newtonian viscosity bugs: mu_max sentinel, shear-rate clamp, D_xz/D_yz boundary fallback, weno_Re_flux validation

Author: sbryngelson

src/simulation/m_checker.fpp

@@ -113,7 +113,7 @@ contains
                 @:PROHIBIT(fluid_pp(i)%nn <= 0._wp, "Non-Newtonian fluid flow behavior index nn must be > 0")
                 @:PROHIBIT(fluid_pp(i)%tau0 < 0._wp, "Non-Newtonian fluid yield stress tau0 must be >= 0")
                 @:PROHIBIT(fluid_pp(i)%mu_min < 0._wp, "Non-Newtonian fluid mu_min must be >= 0")
-                @:PROHIBIT(fluid_pp(i)%mu_max < dflt_real .and. fluid_pp(i)%mu_max <= fluid_pp(i)%mu_min, &
+                @:PROHIBIT(.not. f_is_default(fluid_pp(i)%mu_max) .and. fluid_pp(i)%mu_max <= fluid_pp(i)%mu_min, &
                            & "Non-Newtonian fluid mu_max must be > mu_min when set")
                 @:PROHIBIT(fluid_pp(i)%hb_m <= 0._wp, "Non-Newtonian Papanastasiou parameter hb_m must be > 0")
             end if

src/simulation/m_hb_function.fpp

@@ -42,7 +42,7 @@ contains
 
         mu = yield_term + power_law_term
         mu = max(mu, mu_min_val)
-        if (mu_max_val > 0._wp) mu = min(mu, mu_max_val)
+        if (mu_max_val > dflt_real) mu = min(mu, mu_max_val)
 
     end function f_compute_hb_viscosity
 
@@ -65,9 +65,6 @@ contains
         ! 2*D_ij*D_ij = 2*(D_xx^2+D_yy^2+D_zz^2+2*(D_xy^2+D_xz^2+D_yz^2))
         shear_rate = sqrt(2._wp*(D_xx*D_xx + D_yy*D_yy + D_zz*D_zz + 2._wp*(D_xy*D_xy + D_xz*D_xz + D_yz*D_yz)))
 
-        ! Clamp for numerical safety
-        shear_rate = min(max(shear_rate, 1.0e-2_wp), 1.0e5_wp)
-
     end function f_compute_shear_rate_from_components
 
 end module m_hb_function

src/simulation/m_re_visc.fpp

@@ -130,6 +130,14 @@ contains
                                & l) - q_prim_vf(eqn_idx%mom%beg + 2)%sf(j - 1, k, l))/(x_cc(j + 1) - x_cc(j - 1)))
             else
                 D_xz = 0._wp
+                if (l - 1 >= l_lo .and. l + 1 <= l_hi) then
+                    D_xz = 0.5_wp*(q_prim_vf(eqn_idx%mom%beg)%sf(j, k, l + 1) - q_prim_vf(eqn_idx%mom%beg)%sf(j, k, &
+                                   & l - 1))/(z_cc(l + 1) - z_cc(l - 1))
+                end if
+                if (j - 1 >= j_lo .and. j + 1 <= j_hi) then
+                    D_xz = D_xz + 0.5_wp*(q_prim_vf(eqn_idx%mom%beg + 2)%sf(j + 1, k, &
+                                          & l) - q_prim_vf(eqn_idx%mom%beg + 2)%sf(j - 1, k, l))/(x_cc(j + 1) - x_cc(j - 1))
+                end if
             end if
         else
             D_xz = 0._wp
@@ -143,6 +151,14 @@ contains
                                & l) - q_prim_vf(eqn_idx%mom%beg + 2)%sf(j, k - 1, l))/(y_cc(k + 1) - y_cc(k - 1)))
             else
                 D_yz = 0._wp
+                if (l - 1 >= l_lo .and. l + 1 <= l_hi) then
+                    D_yz = 0.5_wp*(q_prim_vf(eqn_idx%mom%beg + 1)%sf(j, k, l + 1) - q_prim_vf(eqn_idx%mom%beg + 1)%sf(j, k, &
+                                   & l - 1))/(z_cc(l + 1) - z_cc(l - 1))
+                end if
+                if (k - 1 >= k_lo .and. k + 1 <= k_hi) then
+                    D_yz = D_yz + 0.5_wp*(q_prim_vf(eqn_idx%mom%beg + 2)%sf(j, k + 1, l) - q_prim_vf(eqn_idx%mom%beg + 2)%sf(j, &
+                                          & k - 1, l))/(y_cc(k + 1) - y_cc(k - 1))
+                end if
             end if
         else
             D_yz = 0._wp

toolchain/mfc/case_validator.py

@@ -877,6 +877,9 @@ def check_viscosity(self):
             # Check Re(1) requirement
             self.prohibit(Re1 is None and viscous, f"viscous is set to true, but fluid_pp({i})%Re(1) is not specified")
 
+        weno_Re_flux = self.get("weno_Re_flux", "F") == "T"
+        self.prohibit(weno_Re_flux and not viscous, "weno_Re_flux requires viscous to be enabled")
+
     def check_non_newtonian(self):
         """Checks constraints on non-Newtonian fluid parameters"""
         viscous = self.get("viscous", "F") == "T"