Vorticity phi with constraint

[malamast]

Work in progress.

• Added an option in the vorticity component to solve phi as a constrained equation.
• Solving phi with Laplacian inversion in the rhs can be very expensive. However, with the constraint option, you can add phi in the sytem of equations and solve a DAE system.
• This might be useful for 2D simulations with drifts.
• It seems its working for now but I need to test it for longer runs.
• I did not include code for the boundaries as is done in the examples for DAE, since we do not evolve the boundaries.
• The non-boussinesq case is only partially implemented. It will require some improvements there to make it non-boussinesq everywhere.