Adding k-exact properties tests to PyFR operators

General
1

Dear PyFR developers,

In our in-house Spectral-Differences (SD) code, we have recently found a small bug related to the gradient computation. This bug was discovered by testing the k-exact properties of the gradient operator.
I first implemented this test in PyFR to better understand if our error was related to the implementation or the formulation of the method.
This allowed to verify that the error we observed was indeed a bug and not a consequence of the method formulation.
I think that PyFR could benefit from the addition of such tests to help debug future issues, to ease the implementation of new features (ie non-linear cells), to assess very bad-quality meshes, etc.
Let me first briefly describe what a mean by k-exactness of the operators.

Suppose that the solution points of all mesh elements are initialized with a polynomial of the type P(x, y, z) = (x + y + z)^p, where (x, y, z) represents the solution points position in physical space.
A FR scheme which presents local p+1 degree interpolation errors, should be capable of reconstructing the value of this polynomial at any point of the domain from the value of the polynomial at the solution points.
This also implies that the gradient \nabla P should also be exactly reconstructed at any point.
One can validate the numerical implementation of the operators by using the previously described properties.
In particular the following checks seem interesting to me:

  1. Computation of the gradient at the solution points (operator M4)
  2. Interpolation from solution points to flux points (operator M0)
  3. Computation of the divergence at solution points from the values at the flux points (operator M3). Here one should use a vector of polynomials.
  4. Testing the connectivity of the face elements by assessing if the left and right reconstruction states of connected flux points are the same up to machine precision (there is a small caveat for periodic boundaries which I can further describe in future comments). This allows to validate the connectivity of the local and MPI flux points.

Would you be interested in adding these kind of tests to PyFR? I have already implemented some of them. However, I’m sure that there a lot of implementation details to be discussed before including these tests in the code suite.

Regards,

Gonzalo

2

Hi,

So we do currently have some limited tests in PyFR. These can be found here:

https://github.com/PyFR/PyFR/blob/develop/pyfr/tests/test_ele_mats.py

and serve to check the operator matrix construction (this being done for a p = 3 hexahedral element). We could certainly look to expand the scope of these tests to include more elements and polynomial orders (although due to the generic means by which the operators are constructed, if one works the others very likely do, too).

With regards to the k-exact test, I think that this under-stresses the capabilities of some elements. For example (x + y + z)^p does not contain any terms of the form x^py^pz^p which is a mode which can be represented inside of a hexahedral element. Thus, some more thought might be warranted to ensure that we suitably stress the polynomial spaces inside of each element.

Regards, Freddie.

3

Hi Freddie,

Thank you for your answer. I have just made a pull-request with the tests that I described in my first post. I added some discussions on the results in the pull-request answering some of your concerns.

Regards,

Gonzalo