flux reconstruction and transport issues in simplex grid
Summary
I tried running richards+transport and could not get it working.
The cases I tried contain a rain event with solute concentration. DORiE is able to provide the Richards solution, but not the the transport solution and stops in the first time step of the rain event.
EDIT: After Discussion with @lriedel, I could get the simulations running if I allow negative concentrations ([transport.solutionCheck] negativeConcentration = warn). However, the possible issue with the flux reconstruction remains.
I have compiled two examples ('2D_transport' and '3D_transport') that show the issues I had.
It seems like part of the issue comes from the flux reconstruction.
In simplex grids, the reconstructed flux seems to depend on the form of the grid. There are systematic differences between the 'shallow' and the 'deep' part of a grid cell at the surface. This can possibly lead to issues with the transport.
EDIT: I have attached an image showing the issue for the 3D case.
However, I also tried an example with a regular grid ('2d_transport_regular'). This case also fails (but the transport is at least able to make a few time steps). In this case also the resolution is sufficient and should not be an issue. Possibly I am making bad decision in the config file for the transport (however, they are based on the dorie create files). The only way I was able to get this working was by choosing FE order 0.
EDIT: Allowing negative concentrations resolves the regular grid issues. No obvious flux reconstruction issues are visible here.
Steps to reproduce
Please run the provided cases. They all fail. Please have especially a look at the flux reconstruction in '2D_transport' and '3D_transport', which seems rather wrong to me.
What is the current bug behaviour?
Transport fails. Reconstructed flux at surface does not follow boundary condition.
EDIT: Transport does not fail when allowing negative concentrations.
What is the expected correct behaviour?
Transport works. Reconstructed flux is correct and follows boundary condition.
Reproducing input
2D_transport.zip 2D_transport_regular.zip 3D_transport.zip
Ideas how to fix this?
Add them here, if you have any.