tadd c_grad_p to forcing function - 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 5fb2209cc878db151655af8f0681a723063a5957
(DIR) parent 00f5e0875755ad2b080b0bb70d4090450f336c49
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Thu, 9 Oct 2014 15:35:11 +0200
add c_grad_p to forcing function
Diffstat:
M python/shear-results-pressures.py | 14 +++++++++-----
M python/shear-results-strain.py | 10 +++++-----
M python/shear-results-strength.py | 6 +++---
M src/device.cu | 1 +
M src/navierstokes.cuh | 7 ++++---
5 files changed, 22 insertions(+), 16 deletions(-)
---
(DIR) diff --git a/python/shear-results-pressures.py b/python/shear-results-pressures.py
t@@ -66,9 +66,6 @@ for c in numpy.arange(len(c_grad_p)):
#fig = plt.figure(figsize=(8,12))
fig = plt.figure(figsize=(8,15))
-min_p = numpy.min(dev_pres[0])/1000.0
-#max_p = numpy.min(dev_pres)
-max_p = numpy.abs(min_p)
#cmap = matplotlib.colors.ListedColormap(['b', 'w', 'r'])
#bounds = [min_p, 0, max_p]
t@@ -79,12 +76,19 @@ for c in numpy.arange(len(c_grad_p)):
ax.append(plt.subplot(len(c_grad_p), 1, c+1))
+ max_p_dev = numpy.max((numpy.abs(numpy.min(dev_pres[c])),
+ numpy.max(dev_pres[c])))
+ #max_p = numpy.min(dev_pres)
+ min_p = -max_p_dev/1000.0
+ max_p = max_p_dev/1000.0
+
#im1 = ax[c].pcolormesh(shear_strain[c], zpos_c[c], dev_pres[c]/1000.0,
#vmin=min_p, vmax=max_p, rasterized=True)
im1 = ax[c].pcolormesh(shear_strain[c], zpos_c[c], dev_pres[c]/1000.0,
rasterized=True)
- ax[c].set_xlim([0, shear_strain[c,-1]])
- ax[c].set_ylim([zpos_c[0,0], sim.w_x[0]])
+ if c == 0:
+ ax[c].set_xlim([0, numpy.max(shear_strain[c])])
+ ax[c].set_ylim([zpos_c[0,0], sim.w_x[0]])
ax[c].set_xlabel('Shear strain $\\gamma$ [-]')
ax[c].set_ylabel('Vertical position $z$ [m]')
(DIR) diff --git a/python/shear-results-strain.py b/python/shear-results-strain.py
t@@ -13,10 +13,10 @@ from permeabilitycalculator import *
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
-#cvals = ['dry', 1.0, 0.1]
sigma0 = 20000.0
-cvals = ['dry', 1.0]
-step = 1000
+cvals = ['dry', 1.0, 0.1]
+#cvals = ['dry', 1.0]
+step = 800
nsteps_avg = 1 # no. of steps to average over
sim = sphere.sim('halfshear-sigma0=' + str(sigma0) + '-shear')
t@@ -101,9 +101,9 @@ for s in numpy.arange(len(cvals)):
legend = 'wet, c = ' + str(cvals[s])
#ax[0].plot(xdisp[s], zpos_p[s], ',', color = '#888888')
- ax[0].plot(xdisp[s], zpos_p[s], ',', color=color[s], alpha=0.5)
+ #ax[0].plot(xdisp[s], zpos_p[s], ',', color=color[s], alpha=0.5)
ax[0].plot(xdisp_mean[s], zpos_c[s], linetype[s], color=color[s],
- label=legend, linewidth=2)
+ label=legend, linewidth=1)
ax[0].set_ylabel('Vertical position $z$ [m]')
ax[0].set_xlabel('$\\boldsymbol{x}^x_\\text{p}$ [m]')
(DIR) diff --git a/python/shear-results-strength.py b/python/shear-results-strength.py
t@@ -22,9 +22,9 @@ def print_strengths(sid, fluid=False, c=0.0):
tau_ultimate = numpy.average(friction[-500:-1])
if fluid:
- print('%.2f \t %.3f \t %.3f' % (c, tau_peak, tau_ultimate))
+ print('%.2f \t %.2f \t %.2f' % (c, tau_peak, tau_ultimate))
else:
- print('dry \t %.3f \t %.3f' % (tau_peak, tau_ultimate))
+ print('dry \t %.2f \t %.2f' % (tau_peak, tau_ultimate))
return friction
t@@ -32,7 +32,7 @@ def print_strengths(sid, fluid=False, c=0.0):
# print header
-print('$c$ [-] \t Peak \\tau/\\sigma\' [-] \t Ultimate \\tau/\\sigma\' [-]')
+print('$c$ [-] \t \\tau/\\sigma\' (peak) [-] \t \\tau/\\sigma\' (ultimate) [-]')
f = print_strengths(baseid + '-shear', fluid=False)
f = print_strengths(baseid + '-c=1.0-shear', fluid=True, c=1.0)
f = print_strengths(baseid + '-c=0.1-shear', fluid=True, c=0.1)
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -1458,6 +1458,7 @@ __host__ void DEM::startTime()
dev_ns_v_p_z,
nijac,
ns.ndem,
+ ns.c_grad_p,
dev_ns_f1,
dev_ns_f2,
dev_ns_f);
(DIR) diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
t@@ -2493,6 +2493,7 @@ __global__ void findNSforcing(
const Float* __restrict__ dev_ns_v_p_z, // in
const unsigned int nijac, // in
const unsigned int ndem, // in
+ const Float c_grad_p, // in
Float* __restrict__ dev_ns_f1, // out
Float3* __restrict__ dev_ns_f2, // out
Float* __restrict__ dev_ns_f) // out
t@@ -2544,9 +2545,9 @@ __global__ void findNSforcing(
// Find forcing function terms
#ifdef SET_1
- const Float t1 = phi*devC_params.rho_f*div_v_p/dt;
- const Float t2 = devC_params.rho_f*dot(v_p, grad_phi)/dt;
- const Float t4 = dphi*devC_params.rho_f/(dt*dt);
+ const Float t1 = phi*devC_params.rho_f*div_v_p/(c_grad_p*dt);
+ const Float t2 = devC_params.rho_f*dot(v_p, grad_phi)/(c_grad_p*dt);
+ const Float t4 = dphi*devC_params.rho_f/(c_grad_p*dt*dt);
#endif
#ifdef SET_2