Fully 3D cosmological simulations of scalar field dark matter with self-interactions, also known as Bose-Einstein condensate dark matter, are performed using a set of effective hydrodynamic equations. These are derived from the non-linear Schrödinger equation by performing a smoothing operation over scales larger than the de Broglie wavelength, but smaller than the self-interaction Jeans’ length. The dynamics on the de Broglie scale become an effective thermal energy in the hydrodynamic approximation, which is assumed to be subdominant in the initial conditions, but become important as structures collapse and the fluid is shock-heated. The halos that form have Navarro-Frenk-White envelopes, while the centers are cored due to the fluid pressures (thermal + self-interaction), confirming the features found by Dawoodbhoy et al. (2021, MNRAS, 506, 2418) using 1D simulations under the assumption of spherical symmetry. The core radii are largely determined by the self-interaction Jeans’ length, even though the effective thermal energy eventually dominates over the self-interaction energy everywhere, a result that is insensitive to the initial ratio of thermal energy to interaction energy, provided it is sufficiently small to not affect the linear and weakly non-linear regimes. Scaling relations for the simulated population of halos are compared to Milky Way dwarf spheroidals and nearby galaxies, assuming a Burkert halo profile, and are found to not match, although they conform better with observations compared to fuzzy dark matter-only simulations.