FD weights: 17-digit C literals + canonicalized Taylor weights (fixes #2976)

devito/finite_differences/finite_difference.py

@@ -1,6 +1,7 @@
 from collections.abc import Iterable
 from contextlib import suppress
 
+import sympy
 from sympy import sympify
 
 from devito.logger import warning
@@ -22,9 +23,38 @@
     'transpose',
 ]
 
-# Number of digits for FD coefficients to avoid roundup errors and non-deterministic
-# code generation
-_PRECISION = 9
+# Number of digits for FD coefficients. 17 significant digits is the minimum that
+# round-trips an IEEE-754 double exactly, so the C literals the printer emits parse
+# back to the exact rational weight (sub-ULP); evalf is deterministic, so code
+# generation stays deterministic. (The previous value, 9, truncated the weights at
+# ~1.7e-10 relative -- see issue #2976.)
+_PRECISION = 17
+
+
+def _canonicalize_weight(w):
+    """Render a weight's numeric atoms as canonical fixed-precision Floats.
+
+    Taylor weights are mathematically rational, but depending on the route a
+    stencil is built through (e.g. a float ``x0`` shift) they can arrive as
+    binary-float approximations that differ in the last ULP between routes,
+    which would make otherwise-identical stencils compare unequal. Recover the
+    nearby rational when one exists within 5e-14 relative (float-route Fornberg
+    recursions accumulate ~1e-14 relative error; an arbitrary float far from a
+    small rational maps to a rational of the same double value, so
+    non-rational user-supplied coefficients are value-preserved), then render
+    at ``_PRECISION`` so the emitted C literal round-trips the intended double
+    exactly.
+    """
+    def _canon(a):
+        try:
+            r = sympy.Rational(a).limit_denominator(10**12)
+        except (TypeError, ValueError, ZeroDivisionError):
+            return a
+        if a == 0 or abs(float(r) - float(a)) <= 5e-14 * abs(float(a)):
+            return r
+        return a
+
+    return w.replace(lambda a: a.is_Float, _canon).evalf(_PRECISION)
 
 
 @check_input
@@ -168,23 +198,32 @@ def make_derivative(expr, dim, fd_order, deriv_order, side, matvec, x0, coeffici
                                    x0=x0, nweights=nweights)
     # Finite difference weights corresponding to the indices. Computed via the
     # `coefficients` method (`taylor` or `symbolic`)
+    computed_weights = False
     if weights is None:
         weights = fd_weights_registry[coefficients](expr, deriv_order, indices, x0)
         _, wdim, _ = process_weights(weights, expr, dim)
+        computed_weights = coefficients == 'taylor'
     elif isinstance(weights, Iterable) and len(weights) != len(indices):
         warning(f"Number of weights ({len(weights)}) does not match "
                 f"number of indices ({len(indices)}), reverting to Taylor")
         scale = False
         wdim = None
         weights = fd_weights_registry['taylor'](expr, deriv_order, indices, x0)
+        computed_weights = True
 
     # Did fd_weights_registry return a new Function/Expression instead of a values?
     if wdim is not None:
         weights = [weights._subs(wdim, i) for i in range(len(indices))]
 
-    # Enforce fixed precision FD coefficients to avoid variations in results
+    # Enforce fixed precision FD coefficients to avoid variations in results.
+    # Taylor-computed weights are additionally canonicalized (rational
+    # recovery) so route-dependent float error cannot make identical stencils
+    # compare unequal; user-supplied weights are rendered as-is.
     scale = dim.spacing**(-deriv_order) if scale else 1
-    weights = [sympify(scale * w).evalf(_PRECISION) for w in weights]
+    if computed_weights:
+        weights = [_canonicalize_weight(sympify(scale * w)) for w in weights]
+    else:
+        weights = [sympify(scale * w).evalf(_PRECISION) for w in weights]
 
     # Transpose the FD, if necessary
     if matvec == transpose:

tests/test_derivatives.py

@@ -1,5 +1,6 @@
 import numpy as np
 import pytest
+import sympy
 from sympy import Float, Symbol, diff, simplify, sympify
 
 from conftest import assert_structure
@@ -16,7 +17,7 @@
 from devito.types.dimension import StencilDimension
 from devito.warnings import DevitoWarning
 
-_PRECISION = 9
+from devito.finite_differences.finite_difference import _PRECISION  # noqa
 
 
 def x(grid):
@@ -423,15 +424,15 @@ def test_fd_new_side(self):
 
     @pytest.mark.parametrize('so, expected', [
         (2, 'u(x)/h_x - u(x - h_x)/h_x'),
-        (4, '1.125*u(x)/h_x + 0.0416666667*u(x - 2*h_x)/h_x - '
-            '1.125*u(x - h_x)/h_x - 0.0416666667*u(x + h_x)/h_x'),
+        (4, '1.125*u(x)/h_x + 0.041666666666666667*u(x - 2*h_x)/h_x - '
+            '1.125*u(x - h_x)/h_x - 0.041666666666666667*u(x + h_x)/h_x'),
         (6, '1.171875*u(x)/h_x - 0.0046875*u(x - 3*h_x)/h_x + '
-            '0.0651041667*u(x - 2*h_x)/h_x - 1.171875*u(x - h_x)/h_x - '
-            '0.0651041667*u(x + h_x)/h_x + 0.0046875*u(x + 2*h_x)/h_x'),
-        (8, '1.19628906*u(x)/h_x + 0.000697544643*u(x - 4*h_x)/h_x - '
-            '0.0095703125*u(x - 3*h_x)/h_x + 0.0797526042*u(x - 2*h_x)/h_x - '
-            '1.19628906*u(x - h_x)/h_x - 0.0797526042*u(x + h_x)/h_x + '
-            '0.0095703125*u(x + 2*h_x)/h_x - 0.000697544643*u(x + 3*h_x)/h_x')])
+            '0.065104166666666667*u(x - 2*h_x)/h_x - 1.171875*u(x - h_x)/h_x - '
+            '0.065104166666666667*u(x + h_x)/h_x + 0.0046875*u(x + 2*h_x)/h_x'),
+        (8, '1.1962890625*u(x)/h_x + 0.00069754464285714286*u(x - 4*h_x)/h_x - '
+            '0.0095703125*u(x - 3*h_x)/h_x + 0.079752604166666667*u(x - 2*h_x)/h_x - '
+            '1.1962890625*u(x - h_x)/h_x - 0.079752604166666667*u(x + h_x)/h_x + '
+            '0.0095703125*u(x + 2*h_x)/h_x - 0.00069754464285714286*u(x + 3*h_x)/h_x')])
     def test_fd_new_x0(self, so, expected):
         grid = Grid((10,))
         x = grid.dimensions[0]
@@ -560,10 +561,10 @@ def test_transpose_simple(self):
 
         f = TimeFunction(name='f', grid=grid, space_order=4)
 
-        assert str(f.dx.T.evaluate) == ("-0.0833333333*f(t, x - 2*h_x, y)/h_x "
-                                        "+ 0.666666667*f(t, x - h_x, y)/h_x "
-                                        "- 0.666666667*f(t, x + h_x, y)/h_x "
-                                        "+ 0.0833333333*f(t, x + 2*h_x, y)/h_x")
+        assert str(f.dx.T.evaluate) == ("-0.083333333333333333*f(t, x - 2*h_x, y)/h_x "
+                                        "+ 0.66666666666666667*f(t, x - h_x, y)/h_x "
+                                        "- 0.66666666666666667*f(t, x + h_x, y)/h_x "
+                                        "+ 0.083333333333333333*f(t, x + 2*h_x, y)/h_x")
 
     @pytest.mark.parametrize('so', [2, 4, 8, 12])
     @pytest.mark.parametrize('ndim', [1, 2])
@@ -1461,3 +1462,39 @@ def test_unevaluated(self):
         assert Derivative(self.x, self.t)
         assert Derivative(self.x, self.y, self.t)
         assert Derivative(self.x, (self.x, 0))
+
+
+class TestCoefficientPrecision:
+
+    """
+    _PRECISION must round-trip every FD weight to its exact rational value as an
+    IEEE-754 double: the C literal the printer emits is parsed by the compiler to
+    the nearest double, so >= 17 significant digits makes the truncation sub-ULP.
+    With the previous _PRECISION = 9 the order-8 weight 8/315 was emitted as
+    2.53968254e-2 -- a ~1.7e-10 relative coefficient error that silently floored
+    fp64 cross-implementation checks (issue #2976).
+    """
+
+    @pytest.mark.parametrize('p, q', [
+        (-205, 72), (8, 5), (-1, 5), (8, 315), (-1, 560),  # 8th-order d2 weights
+        (-49, 18), (3, 2), (-3, 20), (1, 90),              # 6th-order d2 weights
+        (9, 8), (-1, 24),                                  # 4th-order staggered d1
+    ])
+    def test_fd_weight_literal_roundtrips_fp64(self, p, q):
+        w = sympy.Rational(p, q)
+        emitted = float(sympy.sympify(w).evalf(_PRECISION))
+        exact = p / q
+        assert emitted == exact, (
+            f"{p}/{q}: evalf(_PRECISION={_PRECISION}) -> {emitted!r} != {exact!r} "
+            f"(rel err {abs(emitted - exact) / abs(exact):.3e})")
+
+    def test_evaluated_derivative_floats_roundtrip_fp64(self):
+        grid = Grid(shape=(11, 11))
+        u = Function(name='u_prec', grid=grid, space_order=8)
+        floats = u.dx2.evaluate.atoms(sympy.Float)
+        assert floats, "expected Float weights in the evaluated derivative"
+        for f in floats:
+            r = sympy.nsimplify(f, rational=True, tolerance=1e-8, full=True)
+            assert float(f) == float(r), (
+                f"weight {f} does not round-trip to its rational {r} "
+                f"(rel err {abs(float(f) - float(r)) / abs(float(r)):.3e})")