dsl: Print double-precision FD weights to 20 digits

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')
 

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
@@ -22,9 +23,11 @@
     'transpose',
 ]
 
-# Number of digits for FD coefficients to avoid roundup errors and non-deterministic
-# code generation
-_PRECISION = 9
+# 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
@@ -213,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) 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:

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)
 
 

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,12 +223,28 @@ 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)
         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'])
@@ -603,10 +618,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 +633,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
@@ -941,6 +958,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:
 

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):