Precision of the Interpolation and Gradient operators of Flux Reconstruction methods

[Gonzalo Sáez]

Dear PyFR developers

I have a question regarding regarding the precision of the interpolation and gradient operators of the FR/SD methods in hexahedral meshes.

Let's say that I initialize a variable in the solution points using a general polynomial expression $f$:

$f(x, y, z) = (x + y + z)^p$

Where $p$ is the order of the polynomial and $x, y, z$ is the position in the physical space of a given solution point. Then, I use the interpolation operators (Legendre/Lagrange interpolation basis based on the elements' solution points) to compute the numerical interpolated values of $f$ in the flux/Gauss points $f_g$. I have observed that, if the flux reconstruction scheme is of order $m$ then the interpolation operators will interpolate exactly a function $f$ of order $p = m - 1$ in any type of hexahedral mesh, am I right? Moreover, regarding the approximation of the gradients $\nabla f$ at the flux/gauss points $\nabla f_g$ I have observed that this extrapolated gradient will be numerically exact for $p = m - 1$ in cartesian hexahedral meshes. However, for non-cartesian hexahedral meshes I observed that I can only compute numerically exact the gradient of a function $f$ of order $p = m - 2$. This implies 0th truncation errors of the gradient computation for the FR2/SD2 method in general non-cartesian hexahedral meshes...

Do you know of any references which discuss about the topics that I have described beforehand?

Thank you very much for your help and insight,

Gonzalo