From d0c2d770228bb2ba86ded62a91231374cf39a76b Mon Sep 17 00:00:00 2001 From: Gerard Gorman Date: Mon, 19 Oct 2026 18:26:12 +0000 Subject: [PATCH 1/6] dsl: Compute FD weights at a float x0 exactly A float in x0, e.g. x + 0.5*h_x, made the Taylor weights carry float noise, which only the 9-digit rounding hid. --- devito/finite_differences/tools.py | 4 +++- tests/test_derivatives.py | 10 ++++++---- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/devito/finite_differences/tools.py b/devito/finite_differences/tools.py index 2e2d2847e1..4d816866fb 100644 --- a/devito/finite_differences/tools.py +++ b/devito/finite_differences/tools.py @@ -2,7 +2,7 @@ from itertools import product import numpy as np -from sympy import Expr, Rational, S, cacheit, finite_diff_weights, sympify +from sympy import Expr, Rational, S, cacheit, finite_diff_weights, nsimplify, sympify from devito.logger import warning from devito.tools import Tag, as_tuple @@ -245,6 +245,8 @@ def _numeric_weights(deriv_order, indices, x0): def numeric_weights(function, deriv_order, indices, x0): + # Compute exact weights even for a float `x0`, e.g. `x + 0.5*h_x` + x0 = nsimplify(x0, rational=True, rational_conversion='exact') return _numeric_weights(deriv_order, indices, x0) diff --git a/tests/test_derivatives.py b/tests/test_derivatives.py index 269e5b6e60..f1f2b4ca80 100644 --- a/tests/test_derivatives.py +++ b/tests/test_derivatives.py @@ -603,10 +603,11 @@ def test_fd_adjoint(self, so, ndim, derivative, adjoint_name): b = np.dot(g_deriv.data.reshape(-1), f.data.reshape(-1)) assert np.isclose(1 - a/b, 0, atol=1e-5) + @pytest.mark.parametrize('dtype', [np.float32, np.float64]) @pytest.mark.parametrize('shift, ndim', [(None, 2), (.5, 2), (.5, 3), ((.5, .5, .5), 3)]) - def test_shifted_div(self, shift, ndim): - grid = Grid(tuple([11]*ndim)) + def test_shifted_div(self, shift, ndim, dtype): + grid = Grid(tuple([11]*ndim), dtype=dtype) f = Function(name="f", grid=grid, space_order=4) for order in [None, 2]: df = div(f, shift=shift, order=order).evaluate @@ -617,10 +618,11 @@ def test_shifted_div(self, shift, ndim): ref += getattr(f, f'd{d.name}')(x0=x0, fd_order=order) assert df == ref.evaluate + @pytest.mark.parametrize('dtype', [np.float32, np.float64]) @pytest.mark.parametrize('shift, ndim', [(None, 2), (.5, 2), (.5, 3), ((.5, .5, .5), 3)]) - def test_shifted_grad(self, shift, ndim): - grid = Grid(tuple([11]*ndim)) + def test_shifted_grad(self, shift, ndim, dtype): + grid = Grid(tuple([11]*ndim), dtype=dtype) f = Function(name="f", grid=grid, space_order=4) for order in [None, 2]: g = grad(f, shift=shift, order=order).evaluate From 46106f574d60a9e06a9008a8189d010ae601e42f Mon Sep 17 00:00:00 2001 From: Gerard Gorman Date: Mon, 19 Oct 2026 18:26:12 +0000 Subject: [PATCH 2/6] dsl: Print double-precision FD weights to 20 digits FD weights were rounded to 9 digits whatever the dtype, which floors fp64 accuracy at ~1e-10 relative. Single precision keeps 9 digits. --- devito/finite_differences/finite_difference.py | 11 +++++++++-- tests/test_derivatives.py | 15 +++++++++++++++ 2 files changed, 24 insertions(+), 2 deletions(-) diff --git a/devito/finite_differences/finite_difference.py b/devito/finite_differences/finite_difference.py index 29fa3b4c77..bda2662c56 100644 --- a/devito/finite_differences/finite_difference.py +++ b/devito/finite_differences/finite_difference.py @@ -1,6 +1,7 @@ from collections.abc import Iterable from contextlib import suppress +import numpy as np from sympy import sympify from devito.logger import warning @@ -23,8 +24,10 @@ ] # Number of digits for FD coefficients to avoid roundup errors and non-deterministic -# code generation +# code generation. In double precision, enough digits for the C literals to be +# within one ulp of the exact weights _PRECISION = 9 +_PRECISION_DOUBLE = 20 @check_input @@ -213,7 +216,11 @@ def make_derivative(expr, dim, fd_order, deriv_order, side, matvec, x0, coeffici # Enforce fixed precision FD coefficients to avoid variations in results scale = dim.spacing**(-deriv_order) if scale else 1 - weights = [sympify(scale * w).evalf(_PRECISION) for w in weights] + if expr.dtype in (np.float64, np.complex128): + precision = _PRECISION_DOUBLE + else: + precision = _PRECISION + weights = [sympify(scale * w).evalf(precision) for w in weights] # Transpose the FD, if necessary if matvec == transpose: diff --git a/tests/test_derivatives.py b/tests/test_derivatives.py index f1f2b4ca80..6eef518013 100644 --- a/tests/test_derivatives.py +++ b/tests/test_derivatives.py @@ -943,6 +943,21 @@ def test_unexpand_space_interp_w_saved_timefunc(self): assert_structure(op, ['t,x,y,z', 't,x,y,z,i1', 't,x,y,z,i1,i0']) + @pytest.mark.parametrize('so', [4, 8, 16]) + def test_dx2_exact_double(self, so): + """ + The FD second derivative of `x**2` is exact, so in double precision only + round-off is left (see #2976). + """ + grid = Grid(shape=(31,), extent=(3.,), dtype=np.float64) + u = Function(name='u', grid=grid, space_order=so) + v = Function(name='v', grid=grid, space_order=so) + u.data[:] = np.linspace(0., 3., 31)**2 + + Operator(Eq(v, u.dx2)).apply() + + assert np.allclose(v.data[so//2:-so//2], 2., rtol=0, atol=1e-10) + class TestWeights: From 0dd8ce0f30069f1fd4c8ed6b02e1a835c38ccdb3 Mon Sep 17 00:00:00 2001 From: Gerard Gorman Date: Mon, 19 Oct 2026 20:27:30 +0000 Subject: [PATCH 3/6] dsl: Key the FD weight precision on the dtype --- devito/finite_differences/finite_difference.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/devito/finite_differences/finite_difference.py b/devito/finite_differences/finite_difference.py index bda2662c56..0b96a35e70 100644 --- a/devito/finite_differences/finite_difference.py +++ b/devito/finite_differences/finite_difference.py @@ -27,7 +27,7 @@ # code generation. In double precision, enough digits for the C literals to be # within one ulp of the exact weights _PRECISION = 9 -_PRECISION_DOUBLE = 20 +_PRECISION_DTYPE = {np.float64: 20, np.complex128: 20} @check_input @@ -216,10 +216,7 @@ def make_derivative(expr, dim, fd_order, deriv_order, side, matvec, x0, coeffici # Enforce fixed precision FD coefficients to avoid variations in results scale = dim.spacing**(-deriv_order) if scale else 1 - if expr.dtype in (np.float64, np.complex128): - precision = _PRECISION_DOUBLE - else: - precision = _PRECISION + precision = _PRECISION_DTYPE.get(expr.dtype, _PRECISION) weights = [sympify(scale * w).evalf(precision) for w in weights] # Transpose the FD, if necessary From 09e862c48a17d3a1d7519ba1596b48e0e83430c1 Mon Sep 17 00:00:00 2001 From: Gerard Gorman Date: Tue, 20 Oct 2026 09:51:14 +0000 Subject: [PATCH 4/6] dsl: Turn _PRECISION into a dtype mapper --- devito/finite_differences/finite_difference.py | 11 ++++------- tests/test_derivatives.py | 7 +++---- tests/test_symbolic_coefficients.py | 15 +++++++++------ 3 files changed, 16 insertions(+), 17 deletions(-) diff --git a/devito/finite_differences/finite_difference.py b/devito/finite_differences/finite_difference.py index 0b96a35e70..d105165e66 100644 --- a/devito/finite_differences/finite_difference.py +++ b/devito/finite_differences/finite_difference.py @@ -23,11 +23,9 @@ 'transpose', ] -# Number of digits for FD coefficients to avoid roundup errors and non-deterministic -# code generation. In double precision, enough digits for the C literals to be -# within one ulp of the exact weights -_PRECISION = 9 -_PRECISION_DTYPE = {np.float64: 20, np.complex128: 20} +# Number of digits for FD coefficients: 9 is enough in single precision and keeps code +# generation deterministic; 20 keeps double-precision literals within one ulp +_PRECISION = {np.float32: 9, np.float64: 20, np.complex64: 9, np.complex128: 20} @check_input @@ -216,8 +214,7 @@ def make_derivative(expr, dim, fd_order, deriv_order, side, matvec, x0, coeffici # Enforce fixed precision FD coefficients to avoid variations in results scale = dim.spacing**(-deriv_order) if scale else 1 - precision = _PRECISION_DTYPE.get(expr.dtype, _PRECISION) - weights = [sympify(scale * w).evalf(precision) for w in weights] + weights = [sympify(scale * w).evalf(_PRECISION.get(expr.dtype, 9)) for w in weights] # Transpose the FD, if necessary if matvec == transpose: diff --git a/tests/test_derivatives.py b/tests/test_derivatives.py index 6eef518013..d2be8d1643 100644 --- a/tests/test_derivatives.py +++ b/tests/test_derivatives.py @@ -13,13 +13,12 @@ Add, DiffDerivative, EvalDerivative, IndexDerivative, IndexDerivativeProperty, IndexSum, Weights, interp_for_fd ) +from devito.finite_differences.finite_difference import _PRECISION from devito.symbolics import indexify, retrieve_indexed from devito.types.dimension import StencilDimension from devito.types.grid import DerivedSubDomain from devito.warnings import DevitoWarning -_PRECISION = 9 - def x(grid): return grid.dimensions[0] @@ -204,7 +203,7 @@ def test_derivatives_space(self, derivative, dim, order): else: indices = [(dim + i * dim.spacing) for i in range(-width, width + 1)] - s_expr = u.diff(dim).as_finite_difference(indices).evalf(_PRECISION) + s_expr = u.diff(dim).as_finite_difference(indices).evalf(_PRECISION[u.dtype]) assert(simplify(expr - s_expr) == 0) # Symbolic equality assert type(expr) is EvalDerivative expr1 = s_expr.func(*expr.args) @@ -224,7 +223,7 @@ def test_second_derivatives_space(self, derivative, dim, order): # Establish native sympy derivative expression width = int(order / 2) indices = [(dim + i * dim.spacing) for i in range(-width, width + 1)] - s_expr = u.diff(dim, dim).as_finite_difference(indices).evalf(_PRECISION) + s_expr = u.diff(dim, dim).as_finite_difference(indices).evalf(_PRECISION[u.dtype]) assert(simplify(expr - s_expr) == 0) # Symbolic equality assert type(expr) is EvalDerivative expr1 = s_expr.func(*expr.args) diff --git a/tests/test_symbolic_coefficients.py b/tests/test_symbolic_coefficients.py index 7f4e2e2f29..6a5968f064 100644 --- a/tests/test_symbolic_coefficients.py +++ b/tests/test_symbolic_coefficients.py @@ -23,9 +23,10 @@ class TestSC: ('u.dy2', 4, 2, 1, (0.121, -0.223, 1.648, -2.904, 0), '1.648*u + 0.121*u._subs(y, y - 2*h_y) - 0.223*u._subs(y, y - h_y) \ - 2.904*u._subs(y, y + h_y)')]) - def test_coefficients(self, expr, sorder, dorder, dim, weights, expected): + @pytest.mark.parametrize('dtype', [np.float32, np.float64]) + def test_coefficients(self, expr, sorder, dorder, dim, weights, expected, dtype): """Test that custom coefficients return the expected result""" - grid = Grid(shape=(10, 10)) + grid = Grid(shape=(10, 10), dtype=dtype) u = Function(name='u', grid=grid, space_order=sorder, coefficients='symbolic') x, y = grid.dimensions h_x, h_y = grid.spacing_symbols # noqa @@ -42,7 +43,8 @@ def test_coefficients(self, expr, sorder, dorder, dim, weights, expected): assert isinstance(eq.lhs, Differentiable) s = dim.spacing**(-dorder) - assert sp.simplify(eval(expected).evalf(_PRECISION) * s - eq.evaluate.lhs) == 0 + expected = eval(expected).evalf(_PRECISION[u.dtype]) + assert sp.simplify(expected * s - eq.evaluate.lhs) == 0 def test_function_coefficients(self): """Test that custom function coefficients return the expected result""" @@ -135,9 +137,10 @@ def test_coefficients_expr(self): expr1 = f.dx(w=w).evaluate / x.spacing assert sp.simplify(expr0 - expr1) == 0 - def test_coefficients_w_xreplace(self): + @pytest.mark.parametrize('dtype', [np.float32, np.float64]) + def test_coefficients_w_xreplace(self, dtype): """Test custom coefficients with an xreplace before they are applied""" - grid = Grid(shape=(4, 4)) + grid = Grid(shape=(4, 4), dtype=dtype) u = Function(name='u', grid=grid, space_order=2) x = grid.dimensions[0] h_x = x.spacing @@ -152,7 +155,7 @@ def test_coefficients_w_xreplace(self): s = x.spacing**(-2) expected = (0.1*u - 0.6*u._subs(x, x - h_x) + 0.6*u.subs(x, x + h_x)) * s + 2 - assert sp.simplify(expected.evalf(_PRECISION) - eq.evaluate.lhs) == 0 + assert sp.simplify(expected.evalf(_PRECISION[u.dtype]) - eq.evaluate.lhs) == 0 @pytest.mark.parametrize('order', [2, 4, 6]) def test_staggered_array(self, order): From 9fe210f0b285f58b406bef4a6fa8363b2d8dca1e Mon Sep 17 00:00:00 2001 From: Gerard Gorman Date: Tue, 20 Oct 2026 09:51:14 +0000 Subject: [PATCH 5/6] arch: Turn off FP contraction with safe-math in all GNU languages GCC turns contraction off by default only in ISO C mode (-std=c99), so the C++ backend was fusing multiply-adds under safe-math while the C backend was not. With safe-math on, -ffp-contract=off makes the two match. Fast-math builds are unaffected; bitwise results change for existing safe-math C++ users. --- devito/arch/compiler.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/devito/arch/compiler.py b/devito/arch/compiler.py index 83b57ff51a..9f7d313489 100644 --- a/devito/arch/compiler.py +++ b/devito/arch/compiler.py @@ -494,7 +494,8 @@ def __init_finalize__(self, **kwargs): '-Wno-unused-variable', '-Wno-unused-but-set-variable'] if configuration['safe-math']: - self.cflags.append('-fno-unsafe-math-optimizations') + # GCC disables FMA contraction by default in ISO C, but not in C++ + self.cflags.extend(['-fno-unsafe-math-optimizations', '-ffp-contract=off']) else: self.cflags.append('-ffast-math') From bbb5de2e5afaeb53dcf708d0f6db552369b0b701 Mon Sep 17 00:00:00 2001 From: Gerard Gorman Date: Wed, 21 Oct 2026 08:39:53 +0000 Subject: [PATCH 6/6] dsl: Default the FD weight precision to the single-precision entry Take the fallback digit count from the mapper's `float32` entry instead of an inline literal. An integer or untyped expression has no floating-point precision of its own, so it and any dtype without an entry keep the nine digits their weights had before; half precision gets its own entry. --- devito/finite_differences/finite_difference.py | 13 +++++++++---- tests/test_derivatives.py | 16 ++++++++++++++++ 2 files changed, 25 insertions(+), 4 deletions(-) diff --git a/devito/finite_differences/finite_difference.py b/devito/finite_differences/finite_difference.py index d105165e66..bb81e712a9 100644 --- a/devito/finite_differences/finite_difference.py +++ b/devito/finite_differences/finite_difference.py @@ -23,9 +23,11 @@ 'transpose', ] -# Number of digits for FD coefficients: 9 is enough in single precision and keeps code -# generation deterministic; 20 keeps double-precision literals within one ulp -_PRECISION = {np.float32: 9, np.float64: 20, np.complex64: 9, np.complex128: 20} +# Number of significant digits the FD weights are evaluated to: 9 is enough in single +# (and half) precision and keeps code generation deterministic; 20 keeps +# double-precision literals within one ulp +_PRECISION = {np.float16: 9, np.float32: 9, np.float64: 20, + np.complex64: 9, np.complex128: 20} @check_input @@ -214,7 +216,10 @@ def make_derivative(expr, dim, fd_order, deriv_order, side, matvec, x0, coeffici # Enforce fixed precision FD coefficients to avoid variations in results scale = dim.spacing**(-deriv_order) if scale else 1 - weights = [sympify(scale * w).evalf(_PRECISION.get(expr.dtype, 9)) for w in weights] + # An integer or untyped expression has no floating-point precision of its + # own, so it and any dtype without an entry take the single-precision digits + precision = _PRECISION.get(expr.dtype, _PRECISION[np.float32]) + weights = [sympify(scale * w).evalf(precision) for w in weights] # Transpose the FD, if necessary if matvec == transpose: diff --git a/tests/test_derivatives.py b/tests/test_derivatives.py index d2be8d1643..1e6ca94223 100644 --- a/tests/test_derivatives.py +++ b/tests/test_derivatives.py @@ -229,6 +229,22 @@ def test_second_derivatives_space(self, derivative, dim, order): expr1 = s_expr.func(*expr.args) assert(expr1 == s_expr) # Exact equality + @pytest.mark.parametrize('dtype, precision', [ + (np.float16, _PRECISION[np.float32]), + (np.int32, _PRECISION[np.float32]), + (np.float64, _PRECISION[np.float64]), + ]) + def test_weights_precision_dtype(self, dtype, precision): + """ + The FD weights are evaluated at the digit count of the expression dtype, + with integer expressions taking the single-precision digit count. + """ + x = self.grid.dimensions[0] + u = Function(name='u', grid=self.grid, dtype=dtype, space_order=4) + indices = [(x + i * x.spacing) for i in range(-2, 3)] + s_expr = u.diff(x).as_finite_difference(indices).evalf(precision) + assert simplify(u.dx.evaluate - s_expr) == 0 + @pytest.mark.parametrize('space_order', [2, 4, 6, 8, 10, 12, 14, 16, 18, 20]) # Only test x and t as y and z are the same as x @pytest.mark.parametrize('derivative', ['dx', 'dxl', 'dxr', 'dx2'])