Not only that, but I did not apply the Neumann condition to the same boundary in the two simulations… This is the Fenics solution (complete weak form with cross terms) when the influx boundary condition is applied to the same boundary as in the Elmer simulation:
Seems like things are looking pretty good. Results are very similar. Thank you for bearing with me guys, I think this is solved!