Fixing bug in DiffusionIP implementation

Issue/feature addressed

The Interior Penalty diffusion implementation did not give the same results in serial and parallel. It was found that the issue was related to the "discontinuity penalty factor", which is proportional to 1/h, where h is the average length of the two elements that share an interface (i.e., an edge in 2D or face in 3D). The problem occurred because the length was not computed consistently in serial and in parallel. In addition to this, the length was not computed consistently at periodic boundaries. The reason for this is that the length should be multiplied by a factor of 0.5 at all boundaries except periodic boundaries. In the old implementation, however, the length was multiplied by a factor of 0.5 at parallel interfaces (where MPI communication should be done) and at periodic boundaries. These two issues have now been resolved.

Proposed solution

The solution to the above problems was to not compute the length of the "Bwd" trace inside ExpList::GetElmtNormalLength(). As a result, the lengthBwd array will be empty on all exterior boundaries, all parallel boundaries, and all periodic boundaries. After GetElmtNormalLength() has been called, we fix the missing values in lengthBwd by first copying the lengthFwd into lengthBwd on periodic boundaries, and then do a parallel exchange. After this, only the normal boundaries remain. On these boundaries, which we find by checking the magnitude of lengthBwd, we simply add the 0.5 factor and copy lengthFwd into lengthBwd.

Implementation

For this simple fix, the code should be quite self-explanatory.

Tests

Notes

The current implementation is a a bit "hacky". A nicer solution would be to put GetElmtNormalLength inside the DisContField class instead of the ExpList class, and then add a virtual function inside ExpList which we override inside DisContField, that contains the implementation of GetElmtNormalLength(). We could then raise an error if the implementation of the virtual function in the base class (ExpList) is called, to avoid unintended usage of this function. After all, the current implementation only makes sense if a DG discretisation is used.

Checklist

• Functions and classes, or changes to them, are documented.
• User guide/documentation is updated.
• Changelog is updated.
• Suitable tests added for new functionality.
• Contributed code is correctly formatted. (See the contributing guidelines).
• License added to any new files.
• No extraneous files have been added (e.g. compiler output or test data files).

Warning

On the 19.07 the code formatting (code style) was standardised using clang-format, over the whole Nektar++ code. This means changes in your branch will conflict with formatting changes on the master branch. To resolve these conflicts , see #295 (closed)