In [1]:
# Written for JHTDB by German G Saltar Rivera (2019)
# To use K3D capabilities, use Firefox or Chrome browser.
# Safari has trouble with K3D generated objects
#
import pyJHTDB
from pyJHTDB import libJHTDB
import time as tt
import numpy as np
import k3d #https://github.com/K3D-tools/K3D-jupyter
import ipywidgets as widgets
from ipywidgets import interact, interactive, fixed
import math
from numpy import sin,cos,pi
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)
#Set domain to be queried
time = 0
FD4Lag4 = 44
deltax = 0.01
deltay = 0.008
deltaz = 0.008
nx=100
ny=100
nz=100
xmin, xmax = 0, deltax*nx
ymin, ymax = -1, -1+deltay*ny
zmin, zmax = 0, deltaz*nz
In [2]:
#Creates query points and arranges their coordinates into the required (n,3)-type array
points = np.zeros((nx*ny*nz,3),dtype='float32')
x=np.linspace(xmin,xmax,nx,dtype='float32')
y=np.linspace(ymin,ymax,ny,dtype='float32')
z=np.linspace(zmin,zmax,nz,dtype='float32')
count = 0
for ii in range(np.size(x)):
for jj in range(np.size(y)):
for kk in range(np.size(z)):
points[count,0] = x[ii]
points[count,1] = y[jj]
points[count,2] = z[kk]
count = count + 1
print(np.shape(points))
In [3]:
#Queries the velocity gradient from JHTDB
start = tt.time()
velgrad = lJHTDB.getData(
time, point_coords=points,sinterp = FD4Lag4,
data_set = 'channel',
getFunction = 'getVelocityGradient')
lJHTDB.finalize()
end = tt.time()
print(end - start)
print(velgrad.shape)
In [4]:
#Calculates the q-criterion
qc = np.zeros((np.size(velgrad[:,0]),1))
qc[:,0] = -0.5*(velgrad[:,0]**2+velgrad[:,4]**2+velgrad[:,8]**2+2*(velgrad[:,1]*velgrad[:,3]+velgrad[:,2]*velgrad[:,6]+velgrad[:,5]*velgrad[:,7]))
count2 = 0
qcriterion = np.zeros((nx,ny,nz))
for ii in range(np.size(x)):
for jj in range(np.size(y)):
for kk in range(np.size(z)):
qcriterion[ii,jj,kk] = qc[count2]
count2 = count2 + 1
print(qcriterion.shape)
In [ ]:
#Creates a K3D-volume rendering object
#In order to export the plot as html object, in the controls panel, click on "Snapshot"
vol = k3d.volume(qcriterion, color_range=[2,100], color_map=np.array(k3d.basic_color_maps.Jet,dtype=np.float32),
bounds=[xmin,xmax,ymin,ymax,zmin,zmax],
alpha_coef=100,name="Channel Flow Q vizualization")
plt = k3d.plot()
plt.camera_auto_fit = False
plt.camera = [1.5,0.2,1.5,0,-1,-0.5,0,1,0]
plt += vol
plt.axes = ['x','y','z']
plt.display()
In [ ]:
#Set domain to be queried
#Generates a 3D plot of Q iso-surface with overlayed kinetic energy volume
#rendering in a [0,0.5]^3 subcube in isotropic turbulence
time1 = 3.00
nx1=80
ny1=80
nz1=80
xmin1, xmax1 = 0, 0.5
ymin1, ymax1 = 0, 0.5
zmin1, zmax1 = 0, 0.5
#Creates query points and arranges their coordinates into the required (n,3)-type array
points1 = np.zeros((nx1*ny1*nz1,3),dtype='float32')
x1=np.linspace(xmin1,xmax1,nx1,dtype='float32')
y1=np.linspace(ymin1,ymax1,ny1,dtype='float32')
z1=np.linspace(zmin1,zmax1,nz1,dtype='float32')
count = 0
for ii in range(np.size(x1)):
for jj in range(np.size(y1)):
for kk in range(np.size(z1)):
points1[count,0] = x1[ii]
points1[count,1] = y1[jj]
points1[count,2] = z1[kk]
count = count + 1
In [ ]:
#Queries the velocity from JHTDB
lJHTDB.initialize()
start = tt.time()
Lag4 = 4
vel1 = lJHTDB.getData(
time1, point_coords=points1,sinterp = Lag4,
data_set = 'isotropic1024coarse',
getFunction = 'getVelocity')
end = tt.time()
print(end - start)
lJHTDB.finalize()
print(vel1.shape)
In [ ]:
#Queries the velocity gradient from JHTDB
lJHTDB.initialize()
start = tt.time()
grad1 = lJHTDB.getData(
time1, point_coords=points1,sinterp = FD4Lag4,
data_set = 'isotropic1024coarse',
getFunction = 'getVelocityGradient')
end = tt.time()
print(end - start)
lJHTDB.finalize()
print(grad1.shape)
In [ ]:
#Calculates the q-criterion
q1 = np.zeros((np.size(grad1[:,0]),1))
q1[:,0] = -0.5*(grad1[:,0]**2+grad1[:,4]**2+grad1[:,8]**2+2*(grad1[:,1]*grad1[:,3]+grad1[:,2]*grad1[:,6]+grad1[:,5]*grad1[:,7]))
#Calculates the kinetic energy
e1 = np.zeros((np.size(vel1[:,0]),1))
e1[:,0] = vel1[:,0]**2 + vel1[:,1]**2 + vel1[:,2]**2
In [ ]:
#Arrange 1D arrays into 3D arrays
qcrit = np.zeros((nx1,ny1,nz1))
energ = np.zeros((nx1,ny1,nz1))
count2 = 0
for ii in range(nx1):
for jj in range(ny1):
for kk in range(nz1):
qcrit[ii,jj,kk] = q1[count2]
energ[ii,jj,kk] = e1[count2]
count2 += 1
In [ ]:
#Creates a K3D isosurface object
isosurface = k3d.marching_cubes(qcrit,xmin=xmin1,xmax=xmax1,ymin=ymin1,ymax=ymax1, zmin=zmin1, zmax=zmax1,
level=250, color = 0xf4ea0e,name = 'isotropic: Q Isosurface')
#Creates a K3D volume rendering object
volume = k3d.volume(energ, color_range=[0.1*np.max(energ),0.8*np.max(energ)], color_map=np.array(k3d.basic_color_maps.Jet,dtype=np.float32),
bounds=[xmin1,xmax1,ymin1,ymax1,zmin1,zmax1]
,alpha_coef=15,name="isotropic: kinetic energy")
plot = k3d.plot()
plot.camera_auto_fit = True
plot += volume
plot += isosurface
plot.axes = ['x','y','z']
plot.display()
In [ ]:
In [ ]: