I'll go ahead and say that one can roughly classify systems into
- Systems with a unique ground state.
- Systems that have multiple (possibly infinitely many) degenerate ground states, but that will tend to select a unique one by some mechanism.
- Systems with topological ground state degeneracy.
(I won't rule out there existing additional classes, and I'd welcome learning about them in the comments, but this classification seems satisfactory for the current answer.)
Clearly, class 1 is easy to deal with. There might be some system where it's difficult to numerically resolve the gap to the first excited state, but in practice, with finite size effects and all, the gap will usually be larger than machine resolution. An example of such a system would be spins coupled to a polarizing magnetic field.
Class 2 represents includes some of our most well studied systems, such as the antiferromagnetic Ising model $H_I=\sum_{\langle i,j\rangle} S_i^zS_j^z$ on the square lattice, or the Heisenberg model $H_H=\sum_{\langle i,j\rangle} \mathbf{S}_i \cdot \mathbf{S}_j$ on a cubic lattice. Clearly both models allow for a set of degenerate ground states. You're correct that for small systems, and in principle also larger ones, the system could be in a superposition of different several such states. However, in condensed matter we're usually interested in phases of matter, a concept which is only really well-defined in the thermodynamic limit - as the system size tends to infinity.
In that limit we tend to have spontaneous symmetry breaking - that is, for whatever reason (e.g. a small perturbation from the environment) the system picks out a ground state of lower symmetry than the Hamiltonian. If this wasn't the case, a permanent magnet would be a superposition of a state with all spins along $+\hat{z}$ and a state with all spins along $-\hat{z}$, and averaging the magnetization would yield zero - clearly an unphysical result. In other words, the physically observable quantities in macroscopic systems tend to arise from a symmetry broken ground state, not a superposition of degenerate states. Finally, for completeness, I will also mention that in some systems there is also a phenomenon called order by disorder that helps reduce the ground state degeneracy by entropic mechanisms.
For a DMRG calculation (and some other methods), assuming you want to say something about large systems, you will often want to include a small symmetry breaking term to guide the system to the ground state you're interested in. Often one can pick a tiny magnetic field (uniform for a ferromagnetic state, staggered for a Néel state). This picture looks very clean for the kind of magnetic states I've discussed above, one's related by a symmetry of the Hamiltonian, since the symmetry allows us to pick our coordinate system as we like. (That's essentially what we're doing by adding a magnetic field - picking out the natural spin quantization axis.) Now, the picture would be less clear in a system with accidental (or maybe pathological) degeneracies between qualitatively very different states. In that case picking the wrong state likely leads to predictions that disagree with the experiment, so one could just try again.
For smaller systems, where superposition might be a reality, I'd probably try to stick with exact diagonalization techniques myself. People studying systems from class 3 tend to utilize very sophisticated approaches, see e.g. this paper. But based on your comment I think class 2 interests you more.