FD weights truncated to 9 significant digits in generated C (_PRECISION=9) — floors fp64 cross-validation at ~1e-10

[ggorman]

Summary

Devito prints finite-difference weights into the generated C as ~9-significant-digit float literals, introducing a ~1.7e-10 relative coefficient error vs the exact rational weights. This is physically harmless for production runs, but it silently floors any fp64 cross-implementation / verification gate at ~1e-10 per step: the achievable Devito-vs-exact-reference agreement is set by the printed-literal truncation, not by the discretization or the port under test.

Observed on 4.8.23.dev101+gcf7c3f30a (and earlier 4.8.x).

Root cause

1. devito/finite_differences/finite_difference.py:27:

   _PRECISION = 9

   and line 187:

   weights = [sympify(scale * w).evalf(_PRECISION) for w in weights]

   — the exact rational Taylor weights (e.g. 8/315, -205/72) are materialized as 9-dps sympy.Floats at construction.
2. devito/ir/cgen/printer.py:328 _print_Float honors the Float's own _prec (dps = prec_to_dps(expr._prec)), so the C literal carries only those ~9 digits (e.g. 2.53968254e-2 for 8/315, -2.847222220 for -205/72).

Minimal reproduction

import numpy as np
from devito import Grid, Function, Eq, Operator

grid = Grid(shape=(51, 51), extent=(50., 50.), dtype=np.float64)
u = Function(name="u", grid=grid, space_order=8)
v = Function(name="v", grid=grid, space_order=8)
op = Operator([Eq(v, u.laplace)])
print("2.53968254e-2 in C:", "2.53968254e-2" in str(op.ccode))   # True — 8/315 truncated

# impulse response: measure the printed tap vs the exact rational
u.data[:] = 0.0; u.data[25, 25] = 1.0
op.apply()
h2 = (50./50.)**2
tap = v.data[25, 29] * h2          # the +-4 tap of the 8th-order Laplacian row
exact = -1.0/560.0                 # exact rational weight
print("relative tap error:", abs(tap - exact) / abs(exact))     # ~1.7e-10, not ~1e-16

Measured downstream: a 10-step fp64 leapfrog cross-check against an exact-coefficient numpy implementation drifts to 2.2e-10 (vs 1.6e-15 when the same stencil is expressed with exact coefficients).

Why it matters

In a verification corpus we swept 332 Devito-comparison gates and found 11 gates in 9 projects whose tolerances had been silently loosened into the 1e-7..1e-9 range specifically to swallow this floor (several with in-code comments like "Devito's generated C rounds the rational stencil coefficients"). Each of those gates is blind to real coefficient defects at the magnitude its tolerance implies. The pattern is invisible-by-construction: the gate passes, so nobody re-reads it.

Verified workarounds (either restores fp64)

• Express the stencil as indexed accesses with sympy.Rational coefficients — the printer's _print_Rational emits exact p.0F/q.0F divisions (residual < 1e-14).
• Raise _PRECISION (e.g. to 30): the emitted literals then round-trip fp64 exactly (byte-match at ~1e-16 measured).

Suggested fix

Either bump _PRECISION from 9 to 17 (the fp64 round-trip minimum; the literal is parsed to the nearest double, so ≥17 sig digits makes the truncation sub-ULP), or keep the weights exact (Rational) through lowering and let the printer decide the target-precision rendering. A 17-digit literal has no runtime cost; the only diff is codegen text.

Happy to PR the one-line _PRECISION = 17 change with a regression test (generated-C literal round-trips to the exact rational within 1 ULP) if that direction is acceptable.