When you perform a DNS, no matter what the discretization is (FD, FV, FE, SM) you need to compute the solution over a grid enough fine that the Kolmogorov scale is resolved (or at least the Taylor micro-scale is covered).

In terms of computational cost, that means you must fulfill the condition Re_h =O(1).

To give you a quantitative idea, in 2009 the largest Reynolds number achieved with DNS demanded a mesh with 4096^3 nodes. That is more than 68 billion mesh nodes. (Have a look at:

http://www.annualreviews.org/doi/full/10.1146/annurev.fluid.010908.165203)

The main problem ist that the turbulent structures do not scale up with an enlarged region considered, or, the ratio of the viscous sublayer thickness to the boundary-layer thickness gets smaller with increasing Re number. That means that you have to cope with a multiscale problem, and the ratio of the largest to the smallest scales increases with Re number - i.e. the resolution requirements go up dramatically (~Re**9/4).

Some hot fix is to use a not fully resolved, stabilised gross scale "DNS" which merges, depending on the approach, with a LES. Of course this is not "the DNS" anymore but can deliver better results with a local expert and a high-order code than any other approach.

The type of discretisation concerns only the local expert.

In order to carry out simulations at very high Reynolds number, spectral-element methods are rather common to use. For example, this code:

https://nek5000.mcs.anl.gov/

is commonly used for parallel simulations of turbulent flows.

I recommend using spectral method than other algorithms, if the boundary is not complex. Because this method has higher computational efficiency and needs a smaller resource.