Multidimensional viscous fingering is accurately simulated by an extension of the Hartley spectral methods to three dimensions. Two types of initial conditions were used in the rectilinear displacements−white noise and polygonal cells of the wave number of the mode of maximum instability identified by the linear theory. White noise initial conditions demonstrated that the mechanisms of nonlinear interaction of viscous fingers found in two-dimensional (2-D) simulations persist to three dimensions. Further, the long-time rate of advance of viscous fingers remains unchanged from two dimensions, suggesting that 2-D simulations are sufficient to capture essential features of nonlinear viscous fingering. Simulations with polygonal cellular symmetries, specifically square cells and hexagonal cells, illustrated that weak nonlinear theory cannot predict the shape selection of viscous fingers, as the mechanisms of finger interactions that dominate shape selection lack transverse symmetry and symmetry with respect to the displacement front.