Make rF code more robust

Per Peter, QUICK is not TVD after rF=5 and SOU fails beyond rF=2.
Consequently, a coefficient of 1e6 is pretty ludicrous. Cutting to
1e3 seems to make nonlinear convergence quite a bit better for
a lid problem with symmetry in x-velocity about the vertical center-line.

Refs advection sheme test case introduced in idaholab#20504

framework/include/utils/MathFVUtils.h

@@ -430,10 +430,16 @@ template <typename Scalar, typename Vector>
 Scalar
 rF(const Scalar & phiC, const Scalar & phiD, const Vector & gradC, const RealVectorValue & dCD)
 {
-  if ((phiD - phiC) == 0)
-    return 1e6 * MathUtils::sign(gradC * dCD) + 0 * (gradC * dCD + phiD + phiC);
+  const auto fgrad = phiD - phiC;
+  const auto fgradC = gradC * dCD;
+  constexpr Real coeff = 1e3;
+  static const zero_vec = RealVectorValue(0);
 
-  return 2. * gradC * dCD / (phiD - phiC) - 1.;
+  if (std::abs(fgradC) >= coeff * std::abs(fgrad))
+    // Second term (multiplied by zero) is to keep same sparsity pattern as else branch
+    return 2 * coeff * MathUtils::sign(fgradC) * MathUtils::sign(fgrad) + zero_vec * gradC;
+
+  return 2. * fgradC / fgrad - 1.;
 }
 
 /**