tgranularTemperature.py - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
(HTM) git clone git://src.adamsgaard.dk/sphere
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tgranularTemperature.py (2989B)
---
1 #!/usr/bin/env python
2 # Written by Indraneel Kasmalkar <ineel@stanford.edu>
3 import numpy
4 import sphere
5 """Function to compute granular temperature, namely, velocity fluctuations.
6 @param sim: Sphere object.
7 @param step1: int. step number from which we start estimating velocity.
8 @param step2: int. step number at which we end estimating velocity.
9 @param project2D: bool. If true, the values are averaged over the x-axis to create a 2D array.
10 @returns None.
11 @post: creates a Numpy array sim.Theta which has shape equal to sim.num, and
12 which stores spatial velocity fluctuation of the average velocity from step1 to step2.
13 """
14 def computeGranularTemperature(sim, step1, step2, project2D = True):
15 #Create arrays to store grain velocities.
16 sim.s_vel_x = numpy.zeros(sim.num, dtype = float)
17 sim.s_vel_y = numpy.zeros(sim.num, dtype = float)
18 sim.s_vel_z = numpy.zeros(sim.num, dtype = float)
19 #Array for square of magnitude of velocity
20 sim.m_velsqr = numpy.zeros(sim.num, dtype = float)
21
22 #Read information at steps 1 and 2.
23 sim.readstep(step1)
24 t1 = sim.time_current
25 xyzsum1 = sim.xyzsum
26
27
28 sim.readstep(step2)
29 t2 = sim.time_current
30 xyzsum2 = sim.xyzsum
31
32 #Compute total displacement between t1 and t2.
33 xyzsum = xyzsum2 - xyzsum1
34
35 #Create array for counting grains in each cell, and a map indicating which cell each grain is in.
36 sim.count = numpy.zeros(sim.num, dtype = int)
37 sim.cellInd = sim.x//(sim.L[0]/sim.num[0])
38 sim.cellInd=sim.cellInd.astype(int)
39
40 for i in range(sim.np):
41 #Increase count of the grain for the respective cell location
42 sim.count[sim.cellInd[i,0],sim.cellInd[i,1],sim.cellInd[i,2]] += 1
43
44 #Compute velocities of grains. Will be divided by cell count at the end.
45 sim.s_vel_x[sim.cellInd[i,0],sim.cellInd[i,1],sim.cellInd[i,2]] += xyzsum[i,0]/(t2-t1)
46 sim.s_vel_y[sim.cellInd[i,0],sim.cellInd[i,1],sim.cellInd[i,2]] += xyzsum[i,1]/(t2-t1)
47 sim.s_vel_z[sim.cellInd[i,0],sim.cellInd[i,1],sim.cellInd[i,2]] += xyzsum[i,2]/(t2-t1)
48 sim.m_velsqr[sim.cellInd[i,0],sim.cellInd[i,1],sim.cellInd[i,2]] += (xyzsum[i,0]**2+xyzsum[i,1]**2+xyzsum[i,2]**2)/((t2-t1)**2)
49
50 #To avoid division by zero.
51 sim.count[sim.count==0] = 1
52
53 #Compute average velocities in each direction for each cell.
54 #If doing calculations averaged along x:
55 if project2D:
56 sim.s_vel_x = sim.s_vel_x.sum(axis = 0)/sim.count.sum(axis = 0)
57 sim.s_vel_y = sim.s_vel_y.sum(axis = 0)/sim.count.sum(axis = 0)
58 sim.s_vel_z = sim.s_vel_z.sum(axis = 0)/sim.count.sum(axis = 0)
59 sim.m_velsqr = sim.m_velsqr.sum(axis = 0)/sim.count(axis = 0)
60
61 else:
62 sim.s_vel_x = sim.s_vel_x/sim.count
63 sim.s_vel_y = sim.s_vel_y/sim.count
64 sim.s_vel_z = sim.s_vel_z/sim.count
65 sim.m_velsqr = sim.m_velsqr/sim.count
66
67 #Create sim.Theta.
68 sim.Theta = numpy.sqrt(sim.m_velsqr - sim.s_vel_x**2 - sim.s_vel_y**2 - sim.s_vel_z**2)