  The dynamics of self-gravitating fluid bodies is described by the
Euler-Einstein system of partial differential equations. The break-down of
well-posedness on the fluid-vacuum interface remains a challenging open
problem, which is manifested in simulations of oscillating or inspiraling
binary neutron-stars. We formulate and implement a well-posed canonical
hydrodynamic scheme, suitable for neutron-star simulations in numerical general
relativity. The scheme uses a variational principle by Carter-Lichnerowicz
stating that barotropic fluid motions are conformally geodesic and Helmholtz's
third theorem stating that initially irrotational flows remain irrotational. We
apply this scheme in 3+1 numerical general relativity to evolve the canonical
momentum of a fluid element via the Hamilton-Jacobi equation. We explore a
regularization scheme for the Euler equations, that uses a fiducial atmosphere
in hydrostatic equilibrium and allows the pressure to vanish, while preserving
strong hyperbolicity on the vacuum boundary. The new regularization scheme
resolves a larger number of radial oscillation modes compared to standard,
non-equilibrium atmosphere treatments.
