  In this paper we present a numerical method for hydrodynamic models that
arise from time dependent density functional theories of freezing. The models
take the form of compressible Navier-Stokes equations whose pressure is
determined by the variational derivative of a free energy, which is a
functional of the density field. We present unconditionally energy stable and
mass conserving implicit finite difference methods for the models. The methods
are based on a convex splitting of the free energy and that ensures that a
discrete energy is non-increasing for any choice of time and space step. The
methods are applicable to a large class of models, including both local and
non-local free energy functionals. The theoretical basis for the numerical
method is presented in a general context. The method is applied to problems
using two specific free energy functionals: one local and one non-local
functional. A nonlinear multigrid method is used to solve the numerical method,
which is nonlinear at the implicit time step. The non-local functional, which
is a convolution operator, is approximated using the Discrete Fourier
Transform. Numerical simulations that confirm the stability and accuracy of the
numerical method are presented.
