GAP has an undocumented attribute
DirectFactorsOfGroup which is used by
StructureDescription. With its help, I've calculated that there are 3079 groups for which
DirectFactorsOfGroup can not find factorisation into the direct product, hence the remaining 3459 groups of order 2016 are direct products.
Just in case of
DirectFactorsOfGroup works under some assumptions and may miss some factorisations (after all, it's not documented!), I think it is safe to say that there are AT LEAST 3459 groups of order 2016 which are direct products of smaller groups.
I guess that your calculation of the upper bound returned $7466$ groups:
gap> n:=Filtered([2..46],i -> IsInt(2016/i));
[ 2, 3, 4, 6, 7, 8, 9, 12, 14, 16, 18, 21, 24, 28, 32, 36, 42 ]
which is more than twice more. This is not surprising, since some groups may be represented as direct products in more than one way.
It should be still feasible to generate these 7466 direct products and then classify them up to isomorphism. You need to start to break their set into clusters using some easily computed invariants (e.g. the order of the centre, the number of conjugacy classes etc.) and then refine them using more and more invariants, and apply hard tests only to the refined clusters.