It has been recognized for a while now that the usual, velocity-pressure FEM
can result in \(\mathcal {O}(1)\)  errors in mass conservation, \(||\nabla \cdot u||=\mathcal {O}(1)\) , e.g., [1]} and Akbas, Rebholz and Tone
[2]}. This \(||\nabla \cdot u||=\mathcal {O}(1)\)  is clearly evident in
the \(\gamma =\alpha =0\)  tests in Section 5. The cure for this is added grad-div
stabilization, another simple idea with strong positive consequences. Its
origin seems to be in SUPG type local residual stabilization methods, Brooks
and Hughes [3]}, and the idea of adding an operator positive definite
on the constraint set in optimization. Detailed analysis of the discretization
including an added grad-div term can be found in many papers including Case
and Ervin [4]}, Olshanskii and Reusken [5]}, Olshanskii
[6]}, Jenkins, John, Linke and Rebholz [7]}, Braack,
Bürman, John and Lube [8]}, Layton, Manica, Neda, Olshanskii and
Rebholz [9]}, Galvin, Linke, Rebholz and Wilson [10]} and
Connors, Jenkins and Rebholz [11]}. Preselection of the grad-div
parameter \(\gamma \)  is treated in many places such as Heavner [12]} and
self-adaptive selection in Xie [13]}.
