import numpy and pyJHTDB stuff
In [1]:
import numpy as np
import pyJHTDB
from pyJHTDB.dbinfo import mhd1024, isotropic1024coarse
from pyJHTDB import libJHTDB
now import matplotlib and require plots to be shown inline
In [2]:
%matplotlib inline
import matplotlib.pyplot as plt
Generate points adequate for the isotropic Navier Stokes data set.
In [3]:
info = isotropic1024coarse
nlines = 16
ntimes = 16
spacing = info['dx']
x = np.zeros((nlines, info['nx'], 3), dtype = np.float32)
e = np.random.randn(nlines, 3)
e /= (np.sum(e**2, axis = 1)**.5)[:, None]
x0 = np.random.random(size=(nlines, 3))
x[:, 0, 0] = info['xnodes'][0] + x0[:, 0]*info['lx']
x[:, 0, 1] = info['ynodes'][1] + x0[:, 1]*info['ly']
x[:, 0, 2] = info['znodes'][2] + x0[:, 2]*info['lz']
x[:] = (x[:, 0, None] +
e[:, None]*np.linspace(0,
x.shape[1]*spacing,
num = x.shape[1],
endpoint=False)[None, :, None])
time = np.random.choice(info['time'], size = ntimes)
Get the velocity field at the above points, i.e. on lines. Since the flow is quasistationary, whether we're averaging over space or time shouldn't really matter, therefore we're reshaping the array on the last line to make things easier.
In [4]:
lJHTDB = libJHTDB()
lJHTDB.initialize()
u = []
for t in time:
u.append(lJHTDB.getData(
t,
x,
getFunction = 'getVelocity',
data_set = info['name']))
lJHTDB.finalize()
u = np.array(u).reshape(nlines*ntimes, x.shape[1], -1)
Perform the inverse Fourier transform, and construct an array for the corresponding wavenumbers. While not technically essential, this gives a starting point for the proper treatment of anisotropic grids (such as the channel flow grid).
In [14]:
uk = np.fft.rfft(u, axis = 1) / u.shape[1]
k0 = 2*np.pi / (spacing * x.shape[1])
k = k0*np.linspace(1, uk.shape[1]+1, num = uk.shape[1])
ek = .5*np.average(np.sum(np.abs(uk)**2, axis = 2), axis = 0)
In [17]:
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111)
etaK = ((info['nu']**3)/info['diss'])**.25
ax.plot(k*etaK,
ek / ((info['diss']**(2./3)) * (etaK**(5./3))),
label = '$E(k)\\varepsilon^{-2/3}\\eta_K^{-5/3}$')
ax.plot(k*etaK,
2*(k*etaK)**(-5./3) / 3,
label = '$\\frac{2}{3}(k \\eta_K)^{-5/3}$')
ax.set_xscale('log')
ax.set_yscale('log')
ax.legend(loc = 'best')
Out[17]: