Dear Firedrakers, Attached is a code for a skew symmetric FEM scheme, that should conserve the discrete energy of the system. I have two issues: 1) When I run the code in serial, it conserves energy as expected but when I run it in parallel there seems to be an error. Attached a graph of energy drifts from running on 1/2 cores. 2) To make the code skew symmetric, I have to project physical variables ( density, velocity) into their variational derivatives. Currently this is done by extending the system from (vel, dens) to (vel, dens, vel_deriv, den_deriv) and I include the projection as part of the mixed system. Obviously this is rather expensive as I double the amount of variables used, and I'm looking to use this for more complicated systems. I believe Colin mentioned to Onno, there may be a better way of including this projection step without extending the mixed system. Thanks Will