tgranularTemperature.py: script to compute granular kinetics - 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
---
(DIR) commit 18390320045b78f12e079e53f9534ffbd26879fd
(DIR) parent 7da84647e15a6ce230d9c2f0fcbbe1111fea6bf8
(HTM) Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 19 Aug 2024 16:27:29 +0200
granularTemperature.py: script to compute granular kinetics
Diffstat:
A python/granularTemperature.py | 68 +++++++++++++++++++++++++++++++
1 file changed, 68 insertions(+), 0 deletions(-)
---
(DIR) diff --git a/python/granularTemperature.py b/python/granularTemperature.py
t@@ -0,0 +1,68 @@
+#!/usr/bin/env python
+# Written by Indraneel Kasmalkar <ineel@stanford.edu>
+import numpy
+import sphere
+"""Function to compute granular temperature, namely, velocity fluctuations.
+@param sim: Sphere object.
+@param step1: int. step number from which we start estimating velocity.
+@param step2: int. step number at which we end estimating velocity.
+@param project2D: bool. If true, the values are averaged over the x-axis to create a 2D array.
+@returns None.
+@post: creates a Numpy array sim.Theta which has shape equal to sim.num, and
+which stores spatial velocity fluctuation of the average velocity from step1 to step2.
+"""
+def computeGranularTemperature(sim, step1, step2, project2D = True):
+ #Create arrays to store grain velocities.
+ sim.s_vel_x = numpy.zeros(sim.num, dtype = float)
+ sim.s_vel_y = numpy.zeros(sim.num, dtype = float)
+ sim.s_vel_z = numpy.zeros(sim.num, dtype = float)
+ #Array for square of magnitude of velocity
+ sim.m_velsqr = numpy.zeros(sim.num, dtype = float)
+
+ #Read information at steps 1 and 2.
+ sim.readstep(step1)
+ t1 = sim.time_current
+ xyzsum1 = sim.xyzsum
+
+
+ sim.readstep(step2)
+ t2 = sim.time_current
+ xyzsum2 = sim.xyzsum
+
+ #Compute total displacement between t1 and t2.
+ xyzsum = xyzsum2 - xyzsum1
+
+ #Create array for counting grains in each cell, and a map indicating which cell each grain is in.
+ sim.count = numpy.zeros(sim.num, dtype = int)
+ sim.cellInd = sim.x//(sim.L[0]/sim.num[0])
+ sim.cellInd=sim.cellInd.astype(int)
+
+ for i in range(sim.np):
+ #Increase count of the grain for the respective cell location
+ sim.count[sim.cellInd[i,0],sim.cellInd[i,1],sim.cellInd[i,2]] += 1
+
+ #Compute velocities of grains. Will be divided by cell count at the end.
+ sim.s_vel_x[sim.cellInd[i,0],sim.cellInd[i,1],sim.cellInd[i,2]] += xyzsum[i,0]/(t2-t1)
+ sim.s_vel_y[sim.cellInd[i,0],sim.cellInd[i,1],sim.cellInd[i,2]] += xyzsum[i,1]/(t2-t1)
+ sim.s_vel_z[sim.cellInd[i,0],sim.cellInd[i,1],sim.cellInd[i,2]] += xyzsum[i,2]/(t2-t1)
+ 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)
+
+ #To avoid division by zero.
+ sim.count[sim.count==0] = 1
+
+ #Compute average velocities in each direction for each cell.
+ #If doing calculations averaged along x:
+ if project2D:
+ sim.s_vel_x = sim.s_vel_x.sum(axis = 0)/sim.count.sum(axis = 0)
+ sim.s_vel_y = sim.s_vel_y.sum(axis = 0)/sim.count.sum(axis = 0)
+ sim.s_vel_z = sim.s_vel_z.sum(axis = 0)/sim.count.sum(axis = 0)
+ sim.m_velsqr = sim.m_velsqr.sum(axis = 0)/sim.count(axis = 0)
+
+ else:
+ sim.s_vel_x = sim.s_vel_x/sim.count
+ sim.s_vel_y = sim.s_vel_y/sim.count
+ sim.s_vel_z = sim.s_vel_z/sim.count
+ sim.m_velsqr = sim.m_velsqr/sim.count
+
+ #Create sim.Theta.
+ sim.Theta = numpy.sqrt(sim.m_velsqr - sim.s_vel_x**2 - sim.s_vel_y**2 - sim.s_vel_z**2)