switch CTU to conservative diffusion

Currently we treat diffusion as an external source and compute the old-time diffusion from S_old, use this in the predictor, and then compute the new-time diffusion after the conservative update, from S_new. The final source is then 1/2(diffusion_old + diffusion_new). However, this is not conservative, and no correction is done at C-F interfaces.

We should instead compute the diffusive flux, -k_th grad T and add this to the energy flux. To be second order, we would still need the old-time diffusion in the predictor to the interface states, and we would also need to evaluate this flux from the Godunov state (including the transport coefficients). This would then automatically be conservative and participate in the flux correction.

This approach will also address #369

Start by tracing the current CTU diffusion path through S_old, S_new, the predictor, and the conservative update. Review how the Godunov state and transport coefficients produce the energy flux, then verify that the diffusive flux participates in C-F interface correction and that the final update is conservative and second order.