tadded deleteAllParticles method, debugging fluid solver - 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 478b4e1cdaff74aa7a3269eb7dfc95341d10f368
(DIR) parent 42c1877b521f4be8d2dc40582d20c12cd7cb7b03
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Wed, 2 Apr 2014 15:29:47 +0200
added deleteAllParticles method, debugging fluid solver
Diffstat:
M python/fluidshear.py | 13 ++++++++-----
M python/sphere.py | 24 ++++++++++++++++++++++--
M src/integration.cuh | 1 -
M src/navierstokes.cuh | 16 ++++++++--------
4 files changed, 38 insertions(+), 16 deletions(-)
---
(DIR) diff --git a/python/fluidshear.py b/python/fluidshear.py
t@@ -7,14 +7,14 @@ sid = 'fluidshear'
## Initialization from loose packing to a gravitationally collapsed state
## without fluids
sim = sphere.sim(sid + '-init', np = 2400, fluid = False)
-#sim.cleanup()
+sim.cleanup()
sim.radius[:] = 0.05
sim.initRandomGridPos(gridnum = [12, 12, 9000])
sim.initTemporal(total = 5.0, file_dt = 0.05)
sim.g[2] = -9.81
sim.run(dry=True)
-#sim.run()
-#sim.writeVTKall()
+sim.run()
+sim.writeVTKall()
## Consolidation from a top wall without fluids
sim.readlast()
t@@ -24,10 +24,13 @@ sim.consolidate(normal_stress = 10.0e3)
#sim.initFluid(mu = 17.87e-4, p = 1.0e5, hydrostatic = True) # mu = water at 0 deg C
sim.initTemporal(total = 1.0, file_dt = 0.01)
sim.run(dry=True)
-#sim.run()
-#sim.writeVTKall()
+sim.run()
+sim.writeVTKall()
sim.visualize('walls')
+import sys
+sys.exit()
+
## Shear with fluids
sim.readlast()
sim.sid = sid + '-shear'
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -306,7 +306,7 @@ class sim:
# Smoothing parameter, should be in the range [0.0;1.0[.
# 0.0 = no smoothing.
- self.gamma = numpy.array(0.0)
+ self.gamma = numpy.array(0.5)
# Under-relaxation parameter, should be in the range ]0.0;1.0].
# 1.0 = no under-relaxation
t@@ -641,6 +641,26 @@ class sim:
self.ev = numpy.append(self.ev, ev)
self.p = numpy.append(self.p, p)
+ def deleteAllParticles(self):
+ '''
+ Deletes all particles in the simulation object.
+ '''
+ self.np[0] = 0
+ self.x = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.radius = numpy.ones(self.np, dtype=numpy.float64)
+ self.xysum = numpy.zeros((self.np, 2), dtype=numpy.float64)
+ self.vel = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.fixvel = numpy.zeros(self.np, dtype=numpy.float64)
+ self.force = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.angpos = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.angvel = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.torque = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.es_dot = numpy.zeros(self.np, dtype=numpy.float64)
+ self.es = numpy.zeros(self.np, dtype=numpy.float64)
+ self.ev_dot = numpy.zeros(self.np, dtype=numpy.float64)
+ self.ev = numpy.zeros(self.np, dtype=numpy.float64)
+ self.p = numpy.zeros(self.np, dtype=numpy.float64)
+
def readbin(self, targetbin, verbose = True, bonds = True, devsmod = True,
esysparticle = False):
'''
t@@ -2332,7 +2352,7 @@ class sim:
self.free_slip_bot = numpy.ones(1, dtype=numpy.int32)
self.free_slip_top = numpy.ones(1, dtype=numpy.int32)
- self.gamma = numpy.array(0.0)
+ self.gamma = numpy.array(0.5)
self.theta = numpy.array(1.0)
self.beta = numpy.array(0.0)
self.tolerance = numpy.array(1.0e-8)
(DIR) diff --git a/src/integration.cuh b/src/integration.cuh
t@@ -106,7 +106,6 @@ __global__ void integrate(Float4* dev_x_sorted, Float4* dev_vel_sorted, // Input
angacc = MAKE_FLOAT4(0.0, 0.0, 0.0, 0.0);
}
-
#ifdef EULER
// Forward Euler
// Truncation error O(dt^2) for positions and velocities
(DIR) diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
t@@ -998,7 +998,7 @@ __global__ void findPorositiesVelocitiesDiametersSpherical(
// Save porosity, porosity change, average velocity and average diameter
__syncthreads();
- //phi = 1.0; dphi = 0.0; // disable porosity effects
+ phi = 0.5; dphi = 0.0; // disable porosity effects
const unsigned int cellidx = idx(x,y,z);
dev_ns_phi[cellidx] = phi;
dev_ns_dphi[cellidx] = dphi;
t@@ -1700,11 +1700,11 @@ __global__ void findPredNSvelocities(
// Calculate the predicted velocity
Float3 v_p = v
+ pressure_term
- + 1.0/devC_params.rho_f*div_phi_tau*devC_dt/phi
+ + 0.0*1.0/devC_params.rho_f*div_phi_tau*devC_dt/phi
+ devC_dt*(f_g) // uncomment this line to disable gravity
- - devC_dt*(f_i)
- - v*dphi/phi
- - div_phi_vi_v*devC_dt/phi;
+ - 0.0*devC_dt*(f_i)
+ - 0.0*v*dphi/phi
+ - 0.0*div_phi_vi_v*devC_dt/phi;
// Report velocity components to stdout for debugging
/*const Float3 dv_pres = -ns.beta/devC_params.rho_f*grad_p*devC_dt/phi;
t@@ -1806,9 +1806,9 @@ __global__ void findNSforcing(
// Find forcing function coefficients
//f1 = 0.0;
f1 = div_v_p*devC_params.rho_f/devC_dt
- + dot(grad_phi, v_p)*devC_params.rho_f/(devC_dt*phi)
- + dphi*devC_params.rho_f/(devC_dt*devC_dt*phi);
- f2 = grad_phi/phi;
+ + 0.0*dot(grad_phi, v_p)*devC_params.rho_f/(devC_dt*phi)
+ + 0.0*dphi*devC_params.rho_f/(devC_dt*devC_dt*phi);
+ f2 = 0.0*grad_phi/phi;
// Report values terms in the forcing function for debugging
/*