fix: stable log_cosh difference and prohibit int_comp with viscous/surface_tension

f_thinc_integral_1d used log_cosh(a+h) - log_cosh(a-h) which suffers
catastrophic cancellation in single precision when h (= beta*n_transverse/2)
is small (interface nearly axis-aligned). Replace with the identity
2*atanh(tanh(a)*tanh(h)), which is algebraically equivalent and numerically
stable at all precisions. Fixes the single-precision CI failure on the
2D and 3D MTHINC WENO tests (5126B21F, 4F3722DB).

Also add validator rules prohibiting int_comp > 0 with viscous=T or
surface_tension=T: the split reconstruction paths in m_rhs prevent
s_thinc_compression from ever running in those cases, so permitting the
combination would silently produce unsharpened results.

Author: sbryngelson

src/simulation/m_thinc.fpp

@@ -56,6 +56,18 @@ contains
 
     end function f_log_cosh
 
+    !> @brief Stable difference: log_cosh(a+h) - log_cosh(a-h) = 2*atanh(tanh(a)*tanh(h)). Avoids catastrophic cancellation when h
+    !! is small relative to a.
+    function f_log_cosh_diff(a, h) result(res)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+        real(wp), intent(in) :: a, h
+        real(wp)             :: res
+
+        res = 2._wp*atanh(tanh(a)*tanh(h))
+
+    end function f_log_cosh_diff
+
     !> @brief Analytical 1-D integral of the THINC function
     function f_thinc_integral_1d(a, b) result(res)
 
@@ -66,7 +78,7 @@ contains
         if (abs(b) < verysmall) then
             res = 5e-1_wp*(1._wp + tanh(a))
         else
-            res = 5e-1_wp + (f_log_cosh(a + 5e-1_wp*b) - f_log_cosh(a - 5e-1_wp*b))/(2._wp*b)
+            res = 5e-1_wp + f_log_cosh_diff(a, 5e-1_wp*b)/(2._wp*b)
         end if
 
     end function f_thinc_integral_1d

toolchain/mfc/case_validator.py

@@ -391,6 +391,10 @@ def check_interface_compression(self):
         self.prohibit(int_comp not in [0, 1, 2], "int_comp must be 0 (off), 1 (THINC), or 2 (MTHINC)")
         self.prohibit(int_comp == 2 and n == 0, "int_comp = 2 (MTHINC) requires at least 2D (n > 0)")
         self.prohibit(int_comp != 0 and num_fluids != 2, "int_comp > 0 requires num_fluids = 2")
+        viscous = self.get("viscous", "F") == "T"
+        surface_tension = self.get("surface_tension", "F") == "T"
+        self.prohibit(int_comp != 0 and viscous, "int_comp > 0 is not supported with viscosity (reconstruction is split; compression would be silently skipped)")
+        self.prohibit(int_comp != 0 and surface_tension, "int_comp > 0 is not supported with surface tension (reconstruction is split; compression would be silently skipped)")
 
         recon_type = self.get("recon_type", 1)
         if recon_type == 1:  # WENO