AMReX bug report: `cell_quartic_interp` is not exact for constants (truncated weights)

Draft of a GitHub issue + PR for AMReX-Codes/amrex (Src/AmrCore/AMReX_Interp_C.H). Prepared 12 Oct 2026 by Sergey Ivliev (sergey@peatland.tech); patch: 0001-AmrCore-exact-weights-in-cell_quartic_interp.patch.

Summary #

cell_quartic_interp_x/y/z (the kernels behind amrex::cell_quartic_interp, CellQuartic) hard-code the five stencil weights as 8-digit decimals:

constexpr Array1D<Real,-2,2> c = {Real(0.01708984), Real(-0.12304688),
                                  Real(0.92285156), Real(0.20507812),
                                  Real(-0.02197266)};

Their sum is 0.99999998, not 1. The exact weights of the 5-point Lagrange interpolation at the fine-cell offsets ±1/4 (nodes −2 … 2) are

All exact values are dyadic rationals, so they are representable exactly in binary floating point and their sum is exactly 1; the decimals truncate each of them at 1e-8 and the truncation errors do not cancel.

Effect #

The interpolation is applied as a tensor product (z, then y, then x), so a constant field is interpolated to (0.99999998)³ = 1 − 6.0e-8 of its value in 3D (1 − 4e-8 in 2D), and every interpolated value carries a relative bias of −6e-8 on top of the usual truncation error. Measured in a 3D run restarted with amr.regrid_on_restart: a newly created level holds exactly 0.99999994000000136 × the coarse data for every component, including fields that are not evolved at all.

Where this matters: the same kernel fills the ghost cells of every fine level at every coarse-fine boundary (FillPatch, FillCoarsePatch, the StateDescriptor interpolater in codes such as GRTeclyn/GRChombo-on-AMReX). In a BSSN/CCZ4 numerical-relativity evolution the right-hand side contains terms like χ ∂_i∂_j α and χ R_ij with the conformal factor χ as a prefactor. At a super-horizon start χ ~ 1e6, and a −6e-8 step between the valid cells and the interpolated ghost cells of the finest level is amplified to a kick of order χ·6e-8/dx² in the traceless extrinsic curvature at the faces of the finest box; in our PBH-formation runs (particle matter, 7–8 levels) this grew into O(1) errors in the conformal metric within 0.2 of the collapse time and was independent of the time step, the dissipation, the matter and the amplitude of the perturbation — only of the finest cell size (∝ dx^−2.5). With the exact weights the artefact disappears (residuals 50× smaller, at the level of the genuine truncation error).

Minimal reproducer: interpolate a constant MultiFab (value 1) with cell_quartic_interp (ratio 2, 3D) and inspect the fine data: 0.99999994… instead of 1. Any test that checks FillCoarsePatch of a constant against the constant to machine precision will fail.

Fix #

Replace the literals with the exact rationals (patch attached; the three kernels change identically):

// exact 5-point Lagrange weights at the fine-cell offsets +-1/4 on the nodes -2..2 (sum exactly 1)
constexpr Array1D<Real,-2,2> c = {Real(35.0/2048.0), Real(-63.0/512.0),
                                  Real(945.0/1024.0), Real(105.0/512.0),
                                  Real(-45.0/2048.0)};

The conservative counterpart ccquartic_interp (CellConservativeQuartic) already uses exact dyadic literals (−0.01171875 = −3/256, 0.0859375 = 11/128, 0.5) and is not affected. It may be worth adding a regression test that every Interpolater reproduces a constant field to round-off.

Environment #

AMReX development at 2ceecb3 ("ParmParse: constexpr bound in queryAsDouble range check"), the literals are unchanged on the current development branch (checked 12 Oct 2026); GCC 13, MPI, 3D, double precision.