P2 value-only fit of an affine field is about 100 times less accurate than FDI and PLI (misses 1e-8 in every frame) #399

Summary

The P2 interpolator reproduces an affine field from value constraints about 100 times less accurately than the FDI and PLI interpolators. The error is the same in an unrotated frame at the origin, so it is not a coordinate-frame defect. It is a solver or regularisation accuracy limit.

Found while testing spec 005 (coordinate and fault-frame correctness). The spec's tolerance for this case (SC-002) is 1e-8*F for values and 1e-8*max(1, |gradient|) for gradients. FDI and PLI meet it. P2 does not. Spec 005 decided (C-017) that this is out of its scope and that a separate issue tracks it. The 005 tests use the exact nodal field for P2 for this reason.

Reproduction

No rotation, no translation. The field is f = 2 + 0.2 x - 0.3 y + 0.4 z on the box (0,0,0) to (8,10,12), with 27 value observations, a mesh of 6 x 6 x 6 cells and the default regularisation.

import numpy as np
from loop_common.supports import P2TetMesh, StructuredGrid, TetMesh
from loop_interpolation import FiniteDifferenceInterpolator, P1Interpolator, P2Interpolator

gradient = np.array([0.2, -0.3, 0.4])
xs, ys, zs = np.meshgrid([1, 4, 7], [1, 5, 9], [1, 6, 11], indexing="ij")
points = np.c_[xs.ravel(), ys.ravel(), zs.ravel()].astype(float)
rows = np.c_[points, 2 + points @ gradient, np.ones(len(points))]
query = np.array([[2.2, 3.7, 4.3], [6.1, 8.2, 9.1]])
step = np.array([8, 10, 12]) / 6

engines = {
    "FDI": FiniteDifferenceInterpolator(
        StructuredGrid(origin=np.zeros(3), nsteps_cells=[6, 6, 6], step_vector=step)),
    "PLI": P1Interpolator(
        TetMesh(origin=np.zeros(3), nsteps_cells=[6, 6, 6], step_vector=step)),
    "P2": P2Interpolator(P2TetMesh(origin=np.zeros(3), nsteps=[6, 6, 6], step_vector=step)),
}
for name, interpolator in engines.items():
    interpolator.set_value_constraints(rows)
    interpolator.setup()
    interpolator.solve_system(
        solver="cg", solver_kwargs={"rtol": 1e-10, "atol": 1e-12, "maxiter": 20000})
    value_error = np.abs(interpolator.evaluate_value(query) - (2 + query @ gradient)).max()
    gradient_error = np.abs(interpolator.evaluate_gradient(query) - gradient).max()
    print(f"{name}: max value error {value_error:.2e}, max gradient error {gradient_error:.2e}")

Result:

The limits of the spec are 7.5e-8 for the value (1e-8 times the largest observation, 7.5) and 1e-8 for the gradient. P2 misses the gradient limit by a factor of about 39 and the value limit by about 4 percent. With a 30 degree rotation of the local frame the P2 value error is 2e-7 to 8e-7 and the gradient error is 1.2e-7 to 2.1e-7; FDI and PLI stay below 4e-9 there.

Notes

• The field is affine, so P2 can represent it exactly. The error comes from the fit, not from the basis.
• Hypothesis, not checked: the P2 mesh has about 1935 unknowns for 27 observations, so the regularisation terms (ridge and gradient regularisation) dominate. FDI and PLI have 343 unknowns.
• Possible work: check the regularisation weights of P2Interpolator, a tighter solve, or a documented accuracy level for P2 fits.
• Spec 005, plan point OP7, and the tests test_rebuild_in_new_frame_preserves_world_evaluation (PR fix: saved interpolator frames, copies and stale-frame check, stage S3 (spec 005) #397) refer to this issue.