"""Standalone pipe-root checks; only a C compiler and Python stdlib are needed.""" import ctypes import math import os from pathlib import Path import shlex import shutil import subprocess import tempfile import unittest ROOT = Path(__file__).resolve().parents[1] LAMINAR_END = 89.96829989 def reference_friction(reynolds, roughness): """Independent evaluation of the retained Darcy blend, without its slope.""" laminar = 64 / reynolds if reynolds <= LAMINAR_END: return laminar smooth = (-1.8 * math.log10(6.9 / reynolds)) ** -2 turbulent = smooth if roughness: fully_rough = (-2 * math.log10(roughness / 3.7)) ** -2 weight = 1 / (1 + (180 / (reynolds * roughness)) ** 2) turbulent = (1 - weight) * smooth + weight * fully_rough transition = ((reynolds - LAMINAR_END) / 2741.96700831) ** 8.37293695 return (laminar + transition * turbulent) / (1 + transition) def reference_root(constant, roughness): if constant == 0: return 0.0 low, high = 0.0, 1.0 while high * high * reference_friction(high, roughness) < constant: high *= 2 for _ in range(120): middle = low + (high - low) / 2 if middle == low or middle == high: break if middle * middle * reference_friction(middle, roughness) < constant: low = middle else: high = middle return low + (high - low) / 2 class PipeStatus(ctypes.Structure): _fields_ = [ ('converged', ctypes.c_int), ('iterations', ctypes.c_int), ('bisections', ctypes.c_int), ('relative_residual', ctypes.c_double), ] class NativePipeSolverTests(unittest.TestCase): @classmethod def setUpClass(cls): command = shlex.split(os.environ.get('CC', '')) if not command: compiler = shutil.which('gcc') or shutil.which('clang') if not compiler: raise unittest.SkipTest('A native C compiler is required') command = [compiler] cls.directory = tempfile.TemporaryDirectory(prefix='native-pipe-solver-') cls.addClassCleanup(cls.directory.cleanup) cls.compiler = command source = (ROOT / 'native/components/kernels.c').read_text() cls.library = cls.build_library(source, 'ordinary') # Fault injection only in the temporary test translation unit. The # resistance equation and public production ABI have no test switches. prepared = 'static double pipe_friction_prepared(' derivative = 'static double pipe_friction_derivative(' assert source.count(prepared) == source.count(derivative) == 1 injected = 'int test_pipe_fault_mode=0;\n' + source.replace( prepared, 'static double pipe_friction_prepared_original(', 1) injected = injected.replace('static double pipe_friction(double', r''' static double pipe_friction_prepared(double re,double rr,double rough) { if(test_pipe_fault_mode==5)return 0; /* No upper sign change. */ if(test_pipe_fault_mode==6)return re<2000 ? .01 : 100; return pipe_friction_prepared_original(re,rr,rough); } static double pipe_friction(double''', 1) injected = injected.replace( derivative, 'static double pipe_friction_derivative_original(', 1) injected = injected.replace('static double pipe_checked_solution(', r''' static double pipe_friction_derivative(double re,double rr,double rough,double *df) { double f=pipe_friction_derivative_original(re,rr,rough,df); double slope=2*re*f+re*re*(*df); if(test_pipe_fault_mode==1)*df=(1e6*slope-2*re*f)/(re*re); if(test_pipe_fault_mode==2)*df=NAN; if(test_pipe_fault_mode==3)*df=(-slope-2*re*f)/(re*re); if(test_pipe_fault_mode==4)return NAN; if(test_pipe_fault_mode==6){*df=0;return re<2000 ? .01 : 100;} return f; } static double pipe_checked_solution(''', 1) cls.fault_library = cls.build_library(injected, 'faults') cls.fault_mode = ctypes.c_int.in_dll(cls.fault_library, 'test_pipe_fault_mode') @classmethod def build_library(cls, source, name): directory = Path(cls.directory.name) source_path = directory / (name + '.c') library_path = directory / (name + ('.dll' if os.name == 'nt' else '.so')) source_path.write_text(source) command = cls.compiler + [ '-std=c11', '-O3', '-Wall', '-Wextra', '-Werror', '-ffp-contract=off', '-fno-fast-math', '-shared', ] if os.name != 'nt': command.append('-fPIC') command += ['-I', str(ROOT / 'native/include'), str(source_path), '-lm', '-o', str(library_path)] run = subprocess.run(command, capture_output=True, text=True, timeout=60) if run.returncode: raise AssertionError(run.stderr) library = ctypes.CDLL(str(library_path)) if os.name == 'nt': import _ctypes cls.addClassCleanup(_ctypes.FreeLibrary, library._handle) library.native_pipe_resistance.argtypes = [ ctypes.c_double, ctypes.c_double, ctypes.c_double, ctypes.POINTER(PipeStatus), ] library.native_pipe_resistance.restype = ctypes.c_double return library def solve(self, constant, roughness=0, flow_scale=1, library=None): # Nonzero sentinel fields ensure every call resets a reused status. status = PipeStatus(1, 999, 999, -1) result = (library or self.library).native_pipe_resistance( constant, roughness, flow_scale, ctypes.byref(status)) return result, status def assert_root(self, constant, roughness, result, status): self.assertTrue(status.converged) self.assertTrue(math.isfinite(result)) self.assertLessEqual(status.relative_residual, 1e-9) expected = reference_root(constant, roughness) self.assertLessEqual(abs(result - expected), 1e-12 + expected * 1e-9) residual = abs(result * result * reference_friction(result, roughness) / constant - 1) self.assertLessEqual(residual, 1e-9) def test_independent_bisection_across_reynolds_roughness_and_flow_scales(self): fallbacks = 0 maximum_iterations = 0 for exponent in range(151): reynolds = 10 ** (-6 + exponent * .1) for roughness in (0, 1e-5, 1e-4, 1e-3, 1e-2, .1): constant = reynolds ** 2 * reference_friction(reynolds, roughness) for scale in (1e-12, 1e-6, 1): with self.subTest(reynolds=reynolds, roughness=roughness, scale=scale): result, status = self.solve(constant, roughness, scale) self.assert_root(constant, roughness, result, status) self.assertLess(status.iterations, 20) fallbacks += status.bisections maximum_iterations = max(maximum_iterations, status.iterations) self.assertGreater(fallbacks, 0) self.assertGreater(maximum_iterations, 0) def test_laminar_transition_neighbors_and_extreme_flow_scales(self): for center in (LAMINAR_END, 1000, 2300, LAMINAR_END + 2741.96700831, 4000): for reynolds in (math.nextafter(center, 0), center, math.nextafter(center, math.inf)): for roughness in (0, 1e-5, .1): constant = reynolds ** 2 * reference_friction(reynolds, roughness) for scale in (1e-300, 1e-12, 1, 1e300): with self.subTest(reynolds=reynolds, roughness=roughness, scale=scale): result, status = self.solve(constant, roughness, scale) self.assert_root(constant, roughness, result, status) def test_zero_invalid_inputs_and_unrepresentable_values_fail_explicitly(self): result, status = self.solve(0) self.assertEqual(result, 0) self.assertTrue(status.converged) self.assertEqual(status.relative_residual, 0) result, status = self.solve(64, 0, 1.7e308) self.assertEqual(result, 1) self.assertTrue(status.converged) for constant, roughness, scale in ( (-1, 0, 1), (math.nan, 0, 1), (math.inf, 0, 1), (1, -1, 1), (1, math.nan, 1), (1, math.inf, 1), (1, 0, 0), (1, 0, -1), (1, 0, math.nan), (1, 0, math.inf), (math.ulp(0.0), 0, 1), # A positive root rounds to zero. (1e8, 3.7, 1), # Non-finite roughness limit. (1e100, 0, 1), # Non-finite friction in bracket expansion. (1e308, 0, 1), # Non-finite starting bound. (128, 0, 1.7e308), (1e8, 0, 1.7e308), ): with self.subTest(constant=constant, roughness=roughness, scale=scale): result, status = self.solve(constant, roughness, scale) self.assertTrue(math.isnan(result)) self.assertFalse(status.converged) self.assertLessEqual(status.iterations, 128) self.assertNotEqual(status.relative_residual, -1) def test_inaccurate_or_unusable_newton_slopes_trigger_convergent_bisection(self): try: for mode in (1, 2, 3): self.fault_mode.value = mode for reynolds in (1000, 2800, 1e5, 1e9): for roughness in (0, .1): constant = reynolds ** 2 * reference_friction(reynolds, roughness) with self.subTest(mode=mode, reynolds=reynolds, roughness=roughness): result, status = self.solve(constant, roughness, library=self.fault_library) self.assert_root(constant, roughness, result, status) self.assertGreater(status.bisections, 0) self.assertLess(status.iterations, 128) finally: self.fault_mode.value = 0 def test_nonfinite_equation_missing_bracket_and_float_stagnation_fail(self): try: for mode in (4, 5, 6): self.fault_mode.value = mode with self.subTest(mode=mode): result, status = self.solve(1e6, library=self.fault_library) self.assertTrue(math.isnan(result)) self.assertFalse(status.converged) if mode == 4: self.assertEqual(status.iterations, 1) elif mode == 5: self.assertEqual(status.iterations, 0) else: # A discontinuous equation has no valid root even when # the bracket narrows to adjacent representable values. self.assertGreater(status.bisections, 0) self.assertLess(status.iterations, 128) finally: self.fault_mode.value = 0 if __name__ == '__main__': unittest.main()