Incorrect representation of fluxes in the output
Since fluxes are quite important for the transport model, in
dune-modelling they are calculated directly from the local operator, later reconstructed and stored using Raviart-Thomas elements, meanwhile, in
DORiE they are calculated only for postprocessing purposes using
GridFunctionInterface and later adapted to the vtk output with
VTKGridFunctionAdapter. Details of both are quite diferent, but in the end, I would say they must representing the very same flux!
Steps to reproduce
Enforce a strong Neumann boundary condition.
What is the current bug behaviour?
Here I'm posting the differences I have found between both for a 1D case with infiltration. (left
DORiE VTKAdapter, right
dune-modelling flux reconstruction)
DORiE, the flux is being calculated point-wise, which is wrong because in the DG model they have only meaning in the interface of elements. Moreover, the evaluation of the fluxes in the interface is also wrong because they are being computed without taking into account that the head is a discontinuous quantity. This can be seen in the figure: Although the influx is 5.5e-8, only the "inner" flux of the element is computed in the
DORiE adapter. The missing part of the flux is given by the jump between the boundary and the element face. In this case, with a big jump due to a strong boundary condition, the output already differs by three orders of magnitude.
What is the expected correct behaviour?
If the printing is continuous, then, fluxes must be reconstructed. Another option would be to enforce printing only vertex values in the interfaces and restrict somehow the interpolation in paraview such that is clear it doesn't have a clear continuous representation.