timplemented correction for spatial advection of porosity - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
HTML git clone git://src.adamsgaard.dk/sphere
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
DIR commit f3f8a7ca3858ddfea3d16e2c20abb1b983478231
DIR parent 6613f7d7b2f2e3a9e149fadc5296e484f617496a
HTML Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Tue, 31 Mar 2015 13:22:32 +0200
implemented correction for spatial advection of porosity
Diffstat:
M python/halfshear-darcy-combined.py | 22 ++++++++++++----------
M python/halfshear-darcy-rate-depend… | 32 +++++++++++++++++++++----------
M python/halfshear-darcy-sigma0=8000… | 6 +++---
M src/darcy.cuh | 69 ++++++++++++++++++++++++-------
M src/device.cu | 5 ++++-
M src/sphere.h | 7 ++++---
6 files changed, 99 insertions(+), 42 deletions(-)
---
DIR diff --git a/python/halfshear-darcy-combined.py b/python/halfshear-darcy-combined.py
t@@ -18,7 +18,8 @@ sid = 'halfshear-darcy-sigma0=80000.0-k_c=3.5e-13-mu=1.04e-07-ss=10000.0-A=70000
outformat = 'pdf'
fluid = True
#threshold = 100.0 # [N]
-calculateforcechains = False
+calculateforcechains = True
+calculateforcechainhistory = False
legend_alpha=0.7
linewidth=0.5
t@@ -106,7 +107,7 @@ for i in numpy.arange(nsteps):
#'''
if calculateforcechains:
- if i > 0:
+ if i > 0 and calculateforcechainhistory:
loaded_contacts_prev = numpy.copy(loaded_contacts)
pairs_prev = numpy.copy(sim.pairs)
t@@ -118,17 +119,18 @@ for i in numpy.arange(nsteps):
sim.findCoordinationNumber()
coordinationnumber[i] = sim.findMeanCoordinationNumber()
- nfound = 0
- if i > 0:
- for a in loaded_contacts[0]:
- for b in loaded_contacts_prev[0]:
- if (sim.pairs[:,a] == pairs_prev[:,b]).all():
- nfound += 1;
+ if calculateforcechainhistory:
+ nfound = 0
+ if i > 0:
+ for a in loaded_contacts[0]:
+ for b in loaded_contacts_prev[0]:
+ if (sim.pairs[:,a] == pairs_prev[:,b]).all():
+ nfound += 1;
- nkept[i] = nfound
+ nkept[i] = nfound
+ print nfound
print coordinationnumber[i]
- print nfound
if calculateforcechains:
DIR diff --git a/python/halfshear-darcy-rate-dependence.py b/python/halfshear-darcy-rate-dependence.py
t@@ -22,10 +22,14 @@ fluid = True
#threshold = 100.0 # [N]
-def creep_rheology(friction, n):
+def creep_rheology1(friction, n):
''' return strain rate from friction (tau/N) value '''
return friction**n
+def creep_rheology2(friction, n, A):
+ ''' return strain rate from friction (tau/N) value '''
+ return A*friction**n
+
for sid in sids:
t@@ -79,12 +83,14 @@ for sid in sids:
#popt, pvoc = scipy.optimize.curve_fit(
#creep_rheology, tau/N, shearstrainrate)
popt, pvoc = scipy.optimize.curve_fit(
- creep_rheology, tau_nonzero[idxfit]/N_nonzero[idxfit],
+ creep_rheology1, tau_nonzero[idxfit]/N_nonzero[idxfit],
shearstrainrate_nonzero[idxfit])
print '# Critical state'
print popt
print pvoc
n = popt[0] # stress exponent
+ #A = popt[1] # stress exponent
+ A = 1.
friction = tau_nonzero/N_nonzero
x_min = numpy.floor(numpy.min(friction))
t@@ -93,27 +99,30 @@ for sid in sids:
numpy.linspace(numpy.min(tau_nonzero[idxfit]/N_nonzero[idxfit]),
numpy.max(tau_nonzero[idxfit]/N_nonzero[idxfit]),
100)
- strainrate_fit = friction_fit**n
+ strainrate_fit = A*friction_fit**n
### Consolidated state fit
- idxfit2 = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
- (shearstrainrate_nonzero < 0.1) &
- ((t_nonzero > 0.0) & (t_nonzero < 3.0)))
- #popt, pvoc = scipy.optimize.curve_fit(
- #creep_rheology, tau/N, shearstrainrate)
+ #idxfit2 = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
+ #(shearstrainrate_nonzero < 0.1) &
+ #((t_nonzero > 0.0) & (t_nonzero < 2.5)))
+ idxfit2 = numpy.nonzero((shearstrain_nonzero < 0.1) &
+ (tau_nonzero/N_nonzero < 0.38))
+ '''
popt2, pvoc2 = scipy.optimize.curve_fit(
- creep_rheology, tau_nonzero[idxfit2]/N_nonzero[idxfit2],
+ creep_rheology2, tau_nonzero[idxfit2]/N_nonzero[idxfit2],
shearstrainrate_nonzero[idxfit2])
print '# Consolidated state'
print popt2
print pvoc2
n2 = popt2[0] # stress exponent
+ A2 = popt2[1] # prefactor
friction_fit2 =\
numpy.linspace(numpy.min(tau_nonzero[idxfit2]/N_nonzero[idxfit2]),
numpy.max(tau_nonzero[idxfit2]/N_nonzero[idxfit2]),
100)
- strainrate_fit2 = friction_fit2**n2
+ strainrate_fit2 = A2*friction_fit2**n2
+ '''
### Plotting
t@@ -182,6 +191,9 @@ for sid in sids:
fig = plt.figure(figsize=(3.5,2.5))
ax1 = plt.subplot(111)
CS = ax1.scatter(friction, dilation[idx],
+ #tau_nonzero[idxfit2]/N_nonzero[idxfit2],
+ shearstrainrate_nonzero[idxfit2],
+ c=shearstrain_nonzero[idxfit2], linewidth=0.05,
c=shearstrain_nonzero, linewidth=0.05,
cmap=matplotlib.cm.get_cmap('afmhot'))
DIR diff --git a/python/halfshear-darcy-sigma0=80000.0-k_c=3.5e-13-mu=1.04e-07-ss=10000.0-A=70000.0-f=0.2-fc.txt b/python/halfshear-darcy-sigma0=80000.0-k_c=3.5e-13-mu=1.04e-07-ss=10000.0-A=70000.0-f=0.2-fc.txt
t@@ -1,3 +1,3 @@
mx1.adamsgaard.dk:70 /src/sphere/commit/f3f8a7ca3858ddfea3d16e2c20abb1b983478231.gph:156: line too long