Hi Philipp,

this is an interesting problem, and one that I observed at another occasion when using triangular FEs — when you simulate density driven flow and you have, say, a high density (high-salinity) fluid at the top boundary and a low-density (low-salinity) fluid everywhere else, you observe different element velocities in the FEs at the top boundary, depending if they share one or two FEs with the top boundary, which impacts how density is averaged across the FE. And I don't think you are doing something wrong with the PDE operators (or they are wrong per se for this instance).

Now, your mesh is obviously very coarse and the error should decrease as you refine the mesh. Your problem is essentially 1D and you could use the 1D analytical solution for the diffusion equation to compare your numerical and analytical results at each x-y coordinate, calculate the different between the anayltical and numerical result (in $L^2$ space, I can show you how) and plot it as a function of, say, the average FE size $h$. If plotted in log-log space, you should see how the error decreases, I cannot remember of the top of my head if it is with $\mathcal{O}(h)$ or $\mathcal{O}(h^2)$, i.e. with linearly with $h$ or quadratically with $h$. If you use quadratic FEs, you should see how your error decreases faster (i.e., quadratically or cubically). In the end, it is probably something one needs to live with as any numerical scheme is just an approximation to the PDE and only a small enough $h$ and $\Delta t$ will approximate the equation accurately. But I will check with people at mathematics if there are tricks that can be used in such a case.

Cheers, Sebastian