Gizmo cooling

This MR makes cooling work (maybe) with Gizmo.

Cooling models generally set du/dt, some of them also occasionally set u directly. We chose to keep that for now, although Gizmo does not use u or du/dt internally.

When du/dt is set, we compute dE=(du/dt - du/dt_old)* dt * m and add that to the energy flux. We use the dt stored in flux.dt, which contains the current time step of the particle. Since the cooling is only applied to active particles, this dt will match that particle's time step. The dE will almost immediately be applied to the conserved energy in kick1. du/dt_old is computed from the gradients and primitive variables using the Euler equations:

du/dt = -vec{v}.nabla(u) + P/rho nabla.vec{v} = 1/(gamma-1) vec{v}.[1/rho^2 nabla(rho) - 1/rho nabla(P)] + P/rho nabla.vec{v}

This is also the value returned to the cooling scheme through hydro_get_comoving_internal_energy_dt().

When u is set directly, we directly update E = m * u + Ekin, where Ekin is computed from the current conserved variables. We also update P, so that the change in u is also reflected in the primitive variables.

Testing is in progress, but results so far seem consistent with SPHENIX for all cooling implementations I managed to run. I still need to test this with comoving integration.

@mivkov @yuyttenhove Please break things!