In [1]:
import numpy as np
import pyJHTDB
from pyJHTDB.dbinfo import isotropic1024coarse as info
import time as tt
npoints = 2
nparticles = 2**5
nsteps = 2**7
x = np.zeros(shape = (npoints, nparticles, 3), dtype = np.float32)
x[..., 0] = info['lx']*np.random.random(size = (npoints,))[:, None]
# this value is adequate for channel flow
x[..., 1] = info['ynodes'][info['ynodes'].shape[0]//2]
x[..., 2] = info['lz']*np.random.random(size = (npoints,))[:, None]
In [2]:
from pyJHTDB import libJHTDB
from pyJHTDB.dbinfo import interpolation_code
lJHTDB = libJHTDB()
lJHTDB.initialize()
#Add token
auth_token = "edu.jhu.pha.turbulence.testing-201311" #Replace with your own token here
lJHTDB.add_token(auth_token)
t = info['time'][-1] #final time
dt = info['time'][1] - info['time'][0] # this may be too big
xfull = np.zeros(shape = (nsteps+1, npoints, nparticles, 3),
dtype = np.float32)
xfull[0] = x
kappa = (2*info['nu'])**.5
for tindex in range(nsteps):
print('step {0}'.format(tindex))
# get velocity
u = lJHTDB.getData(
t,
xfull[tindex],
sinterp = interpolation_code['M2Q8'],
tinterp = interpolation_code['NoTInt'],
data_set = info['name'],
getFunction = 'getVelocity')
# Euler Maruyama
dW = np.random.randn(*xfull.shape[1:])*(dt**.5)
xfull[tindex+1] = xfull[tindex] - u*dt + kappa*dW
t -= dt
lJHTDB.finalize()
In [4]:
import numpy as np
import ipyvolume as ipv
fig = ipv.figure()
#ax = fig.add_subplot(111, projection = '3d')
for traj in range(xfull.shape[2]):
p=ipv.pylab.plot(xfull[:, 0, traj, 0],
xfull[:, 0, traj, 1],
xfull[:, 0, traj, 2],
color = 'red')
#print(xfull[:, 0, :, 0].min())
ipv.xlim(xfull[:, 0, :, 0].min(),xfull[:, 0, :, 0].max())
ipv.ylim(xfull[:, 0, :, 1].min(),xfull[:, 0, :, 1].max())
ipv.zlim(xfull[:, 0, :, 2].min(),xfull[:, 0, :, 2].max())
ipv.show()
In [5]:
lJHTDB = libJHTDB()
lJHTDB.initialize()
#Add token
lTDB.add_token(auth_token)
start = tt.time()
nx=128
ny=512
nz=128
result = lJHTDB.getRawData(
9,
start = np.array([0, 0, 0], dtype = np.int),
size = np.array([nx, ny, nz], dtype = np.int),
data_set = 'channel',
getFunction = 'Velocity')
end = tt.time()
print(end - start)
lJHTDB.finalize()
print(result.shape)
In [6]:
#fig = ipv.figure()
#p=ipv.pylab.volshow(result[:,:,:,0])
#ipv.show()
ipv.quickvolshow(result[:,:,:,0], level=[0.2, 0.5, 1.0], opacity=[0.1, 0.1, 0.1], extent=[[0, 1.5],[-1,1],[0,0.78]])
ipv.xlim(0, 1.5)
ipv.ylim(-1, 1)
ipv.zlim(0, 0.78)
ipv.show()
In [7]:
fig = ipv.figure()
ipv.pylab.plot_isosurface(result[:,:,:,0], level=0.5, extent=[[0, 1.5],[-1,1],[0,0.78]])
#ipv.quickvolshow(result[:,:,:,0], level=[0.2, 0.5, 1.0], opacity=[0.1, 0.1, 0.1])
ipv.xlim(0, 1.5)
ipv.ylim(-1, 1)
ipv.zlim(0, 0.78)
ipv.show()
In [9]:
import ipyvolume as ipv
skip=20
lengthx=int(np.ceil(nx/skip))
lengthy=int(np.ceil(ny/skip))
lengthz=int(np.ceil(nz/skip))
x=np.linspace(0, 8*np.pi, num=2049)
x=x[0:nx:skip]
y=pyJHTDB.dbinfo.channel['ynodes']
y=y[0:ny:skip]
z=np.linspace(0, 3*np.pi, num=1537)
z=z[0:nz:skip]
print(z)
print(lengthy)
y,x,z=np.meshgrid(y,x,z)
x=np.reshape(x,lengthx*lengthy*lengthz)
y=np.reshape(y,lengthx*lengthy*lengthz)
z=np.reshape(z,lengthx*lengthy*lengthz)
u=np.reshape(result[0:nx:skip,0:ny:skip,0:nz:skip,0],lengthx*lengthy*lengthz)
v=np.reshape(result[0:nx:skip,0:ny:skip,0:nz:skip,1],lengthx*lengthy*lengthz)
w=np.reshape(result[0:nx:skip,0:ny:skip,0:nz:skip,2],lengthx*lengthy*lengthz)
fig = ipv.figure()
quiver = ipv.pylab.quiver(x,y,z,u,v,w,size=5**u)
ipv.xlim(0, 1.5)
ipv.ylim(-1, 1)
ipv.zlim(0, 0.78)
#ipv.show()
from ipywidgets import FloatSlider, ColorPicker, VBox, jslink
size = FloatSlider(min=0, max=30, step=0.1)
size_selected = FloatSlider(min=0, max=30, step=0.1)
color = ColorPicker()
color_selected = ColorPicker()
jslink((quiver, 'size'), (size, 'value'))
jslink((quiver, 'size_selected'), (size_selected, 'value'))
jslink((quiver, 'color'), (color, 'value'))
jslink((quiver, 'color_selected'), (color_selected, 'value'))
VBox([ipv.gcc(), size, size_selected, color, color_selected])
In [ ]: