Testing and visualizing large index arrays

Question

I will be implementing nodal discontinuous Galerkin method soon, and having done this before I know the basic indexing arrays I will need to compute, given a mesh and polynomial data.

The problem I ran into in previous code were subtle mistakes I made in computing things like interior/exterior trace indexing. Problems which didn't arise on simpler test cases would arise on larger meshes, and usually this yields an unstable scheme since boundary conditions are not properly imposed (so no chance of just watching the simulation every 10 steps or so and seeing a localized problem).

I'm hoping some more experienced folks here know good tests to run on the index arrays to get some confidence that they are right. Quadratures, derivatives and the like are very easy to test, but other things I can't figure out.

Some tests I have done in the past is adding interior normals to exterior normals, which should yield 0 or +-2. Being able to quickly see the result of some code change is helpful, but I can't think of a meaningful way to do this with indexing.

I should also mention that these are going in for quads and hexes, with the potential for curvilinear elements. Not much existing code is there to compare against a working library already.

Bonus points if there are good unit tests that I can write which wouldn't rely on an existing correct answer to compare against. I'll settle for a lot of good heuristics.

Answer 1

The obvious answer would of course be to not implement it all yourself but to use what others have done before. For example, deal.II (http://www.dealii.org -- disclaimer: I am one of the authors of this library) already has DG elements on quads, hexes, and curved elements, tutorial programs that show how to use them, etc. All this has been tested for a decade already and is tested via 2,700 regression tests after every change to the repository. There is much work to be saved for you.

But to come back to the original question: What I often do for these sort of things when writing test cases for deal.II is to do something for which I know the exact solution and verify that I indeed get it. Example: If you have an interpolating (nodal) element, then on one cell $K$, assign to every degree of freedom the value of a function $f(x)$ at the nodes $x_i$. Do this by hand, using your knowledge of how degrees of freedom are ordered. This yields a function $u_h(x)=\sum_i U_i \varphi_i(x)$ where $U_i=f(x_i)$. Then use a quadrature formula to compute something like $\|u_h-f\|^2 = \int_K |u_h(x)-f(x)|^2 \; dx$. If $f$ was a polynomial of the same order as your ansatz space, then the result should be zero. If you mixed up indices, either when assigning values to degrees of freedom by hand or in your program when evaluating $u_h$ at quadrature points, you will get something that is nonzero.

The point of all this is that in your test, you use two independent methods to deal with indexing (your knowledge, when assigning values to DoFs by hand, and what your code thinks is correct when evaluating $u_h$) and they need to agree. If they don't, at least one of them is incorrect.

Of course, once you have this running, you need to save the testcase and run it periodically to ensure that you are not accidentally breaking this kind of functionality again.