tMerge branch 'master' of src.adamsgaard.dk:src/sphere - 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 4ad7f4d8bd40be3988fec86661de48699a1ca69b
DIR parent a944c7e38af5c26553f8add9df01b56fe98c4c01
HTML Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Tue, 10 Dec 2019 13:42:45 +0100
Merge branch 'master' of src.adamsgaard.dk:src/sphere
Diffstat:
A python/hansen-zoet.py | 152 +++++++++++++++++++++++++++++++
A python/supraglacial-master.py | 159 +++++++++++++++++++++++++++++++
A python/supraglacial-plots.py | 150 +++++++++++++++++++++++++++++++
M src/CMakeLists.txt | 6 ++++--
4 files changed, 465 insertions(+), 2 deletions(-)
---
DIR diff --git a/python/hansen-zoet.py b/python/hansen-zoet.py
t@@ -0,0 +1,152 @@
+#!/usr/bin/env python
+
+# Import sphere functionality
+import sphere
+
+### EXPERIMENT SETUP ###
+initialization = False
+consolidation = False
+shearing = True
+rendering = False
+plots = True
+
+# Number of particles
+np = 1e4
+
+# Common simulation id
+sim_id = "hz"
+
+# Deviatoric stress [Pa]
+Nlist = [51e3, 101e3, 202e3, 303e3, 404e3]
+
+### INITIALIZATION ###
+
+# New class
+init = sphere.sim(np = np, nd = 3, nw = 0, sid = sim_id + "-init")
+
+# Save radii
+init.generateRadii(mean = 800e-5)
+
+# Use default params
+init.defaultParams(gamma_n = 100.0, mu_s = 0.6, mu_d = 0.6)
+
+# Add gravity
+init.g[2] = -9.81
+
+# Periodic x and y boundaries
+init.periodicBoundariesXY()
+
+# Initialize positions in random grid (also sets world size)
+hcells = np**(1.0/3.0)
+init.initRandomGridPos(gridnum = [hcells, hcells, 1e9])
+
+# Set duration of simulation
+init.initTemporal(total = 5.0)
+
+if (initialization == True):
+
+ # Run sphere
+ init.run(dry = True)
+ init.run()
+
+ if (plots == True):
+ # Make a graph of energies
+ init.visualize('energy')
+
+ init.writeVTKall()
+
+ if (rendering == True):
+ # Render images with raytracer
+ init.render(method = "angvel", max_val = 0.3, verbose = False)
+
+
+
+# For each normal stress, consolidate and subsequently shear the material
+for N in Nlist:
+
+ ### CONSOLIDATION ###
+
+ # New class
+ cons = sphere.sim(np = init.np, nw = 1, sid = sim_id +
+ "-cons-N{}".format(N))
+
+ # Read last output file of initialization step
+ lastf = sphere.status(sim_id + "-init")
+ cons.readbin("../output/" + sim_id + "-init.output{:0=5}.bin".format(lastf), verbose=False)
+
+ # Periodic x and y boundaries
+ cons.periodicBoundariesXY()
+
+ # Setup consolidation experiment
+ cons.consolidate(normal_stress = N)
+ cons.adaptiveGrid()
+
+ # Set duration of simulation
+ cons.initTemporal(total = 1.5)
+
+ """
+ cons.w_m[0] *= 0.001
+ cons.mu_s[0] = 0.0
+ cons.mu_d[0] = 0.0
+ cons.gamma_wn[0] = 1e4
+ cons.gamma_wt[0] = 1e4
+ cons.contactmodel[0] = 1
+ """
+
+ if (consolidation == True):
+
+ # Run sphere
+ cons.run(dry = True) # show values, don't run
+ cons.run() # run
+
+ if (plots == True):
+ # Make a graph of energies
+ cons.visualize('energy')
+ cons.visualize('walls')
+
+ cons.writeVTKall()
+
+ if (rendering == True):
+ # Render images with raytracer
+ cons.render(method = "pres", max_val = 2.0*N, verbose = False)
+
+
+ ### SHEARING ###
+
+ # New class
+ shear = sphere.sim(np = cons.np, nw = cons.nw, sid = sim_id +
+ "-shear-N{}".format(N))
+
+ # Read last output file of initialization step
+ lastf = sphere.status(sim_id + "-cons-N{}".format(N))
+ shear.readbin("../output/" + sim_id +
+ "-cons-N{}.output{:0=5}.bin".format(N, lastf),
+ verbose = False)
+
+ # Periodic x and y boundaries
+ shear.periodicBoundariesXY()
+
+ # Setup shear experiment
+ shear.shear(shear_strain_rate = 0.10)
+ shear.adaptiveGrid()
+ #shear.initFluid(mu=17.87e-4, p=0.0, hydrostatic=True, cfd_solver=1)
+
+ # Set duration of simulation
+ shear.initTemporal(total = 10.0)
+
+ if (shearing == True):
+
+ # Run sphere
+ shear.run(dry = True)
+ shear.run()
+
+ if (plots == True):
+ # Make a graph of energies
+ shear.visualize('energy')
+ shear.visualize('shear')
+
+ shear.writeVTKall()
+
+ if (rendering == True):
+ # Render images with raytracer
+ shear.render(method = "pres", max_val = 2.0*N, verbose = False)
DIR diff --git a/python/supraglacial-master.py b/python/supraglacial-master.py
t@@ -0,0 +1,159 @@
+#!/usr/bin/env python
+
+# sphere grain/fluid simulation: https://src.adamsgaard.dk/sphere
+import sphere
+import numpy
+
+### EXPERIMENT SETUP ###
+initialization = False
+creeping = True
+plots = True
+
+# Common simulation id
+sim_id = "supraglacial"
+
+# Fluid-pressure gradient [Pa/m]
+dpdz = 0.0
+
+# Grain density
+rho_g = 3600.0
+
+# Fluid density
+rho_f = 1000.0
+
+# Gravitational acceleration
+g = 9.8
+
+# Slope
+slope_angle = 20.0
+
+# Number of particles
+np = 1e4
+
+device = 0 # automatically choose best GPU
+
+
+### INITIALIZATION ###
+
+# New class
+init = sphere.sim(np = np, nd = 3, nw = 0, sid = sim_id + "-init")
+
+# Uniform diameters from 0.3 cm to 0.7 cm
+init.generateRadii(psd = 'uni', mean = 0.0025, variance = 0.001)
+
+# Use default params
+init.defaultParams(gamma_n = 100.0, mu_s = 0.6, mu_d = 0.6)
+init.setYoungsModulus(1e8)
+
+# Disable wall viscosities
+init.gamma_wn[0] = 0.0
+init.gamma_wt[0] = 0.0
+
+# Add gravity
+init.g[2] = -g
+
+# Periodic x and y boundaries
+init.periodicBoundariesXY()
+
+# Initialize positions in random grid (also sets world size)
+hcells = np**(1.0/3.0)
+init.initRandomGridPos(gridnum = [hcells-2, hcells-2, 1e9])
+
+# Set duration of simulation
+init.initTemporal(total = 5.0)
+
+if (initialization == True):
+
+ # Run sphere
+ init.run(dry = True)
+ init.run(device=device)
+
+ if (plots == True):
+ # Make a graph of energies
+ init.visualize('energy')
+
+ init.writeVTKall()
+
+
+### CREEP ###
+
+# New class
+creep = sphere.sim(np = init.np,
+ sid = sim_id + "-slope{}-dpdz{}".format(slope_angle, dpdz))
+
+# Read last output file of initialization step
+creep.readbin("../output/" + sim_id + "-init.output{:0=5}.bin"
+ .format(init.status()))
+
+# Tilt gravity
+creep.g[2] = -g*numpy.cos(numpy.deg2rad(slope_angle))
+creep.g[0] = g*numpy.sin(numpy.deg2rad(slope_angle))
+
+# Disable particle contact viscosities
+creep.gamma_n[0] = 0.0
+creep.gamma_t[0] = 0.0
+
+# zero all velocities and accelerations
+creep.zeroKinematics()
+
+# Periodic x and y boundaries
+creep.periodicBoundariesXY()
+
+# Fit grid to grains
+creep.adjustUpperWall(z_adjust=1.2)
+creep.nw = 0 # no dynamic wall on top
+
+# Fix bottom grains
+z_min = numpy.min(creep.x[:,2] - creep.radius)
+z_max = numpy.max(creep.x[:,2] + creep.radius)
+d_max_below = numpy.max(creep.radius[numpy.nonzero(creep.x[:,2] <
+ (z_max-z_min)*0.3)])*2.0
+I = numpy.nonzero(creep.x[:,2] < (z_min + d_max_below))
+creep.fixvel[I] = 1
+
+# set fluid and solver properties
+creep.initFluid(mu=8.9e-4, p=0.0, rho=rho_f, cfd_solver=1) # water at 25 C
+creep.setMaxIterations(2e5)
+creep.setPermeabilityGrainSize()
+creep.setFluidCompressibility(4.6e-10) # water at 25 C
+
+# set fluid BCs
+# creep.setFluidTopNoFlow()
+# creep.setFluidBottomNoFlow()
+# creep.setFluidXFixedPressure()
+# creep.adaptiveGrid()
+creep.setFluidTopFixedPressure()
+creep.setFluidBottomFixedPressure()
+creep.setFluidXPeriodic()
+creep.setFluidYPeriodic()
+
+# set fluid pressures at the boundaries and internally
+dz = creep.L[2]/creep.num[2]
+for iz in range(creep.num[2]):
+ z = iz + 0.5*dz
+ creep.p_f[:,:,iz] = numpy.abs(creep.L[2]*dpdz) + z*dpdz
+
+# Remove fixvel constraint from uppermost grains
+#creep.fixvel[numpy.nonzero(creep.x[:,2] > 0.5*creep.L[2])] = 0
+
+# Produce regular coloring pattern
+creep.checkerboardColors(creep.num[0], creep.num[1], creep.num[2])
+creep.color[numpy.nonzero(creep.fixvel == 1)] == -1
+
+# Adapt grid size during progressive deformation
+#creep.adaptiveGrid()
+
+# Set duration of simulation
+creep.initTemporal(total=5.0, file_dt=0.01)
+
+if (creeping == True):
+
+ # Run sphere
+ creep.run(dry = True)
+ creep.run(device=device)
+
+ if (plots == True):
+ # Make a graph of energies
+ creep.visualize('energy')
+
+ creep.writeVTKall()
DIR diff --git a/python/supraglacial-plots.py b/python/supraglacial-plots.py
t@@ -0,0 +1,150 @@
+#!/usr/bin/env python
+
+import sphere
+import numpy as np
+import matplotlib
+matplotlib.rcParams.update({'font.size': 12})
+import matplotlib.pyplot as plt
+import subprocess
+
+dpdz_values = [0.0, -50.0, -100.0]
+slope_angle_values = [5.0, 10.0, 15.0, 20.0]
+
+scatter_color = '#666666'
+scatter_alpha = 0.6
+dp_dz_vals = np.empty(len(dpdz_values)*len(slope_angle_values))
+slope_angles = np.empty_like(dp_dz_vals)
+fluxes = np.empty_like(dp_dz_vals)
+
+'''
+j = 0
+for dpdz in dpdz_values:
+ outfiles = ''
+ for slope_angle in slope_angle_values:
+
+ sim = sphere.sim("supraglacial-slope{}-dpdz{}".format(slope_angle, dpdz),
+ fluid=True)
+ print('### ' + sim.id())
+ print('Last output file: ' + str(sim.status()))
+ sim.readTime(4.99)
+
+ fig = plt.figure()
+
+ title = 'slope = ' + str(slope_angle) + '$^\circ$, ' + \
+ '$dp/dz$ = -' + str(dpdz) + ' Pa/m'
+
+ z = np.zeros(20)
+ v_x_space_avg = np.empty_like(z)
+ xsum_space_avg = np.empty_like(z)
+ dz = np.max(sim.x[:,2])/len(z)
+ for i in range(len(z)):
+ z[i] = i*dz + 0.5*dz
+ I = np.nonzero((sim.x[:,2] >= i*dz) & (sim.x[:,2] < (i+1)*dz))
+ v_x_space_avg[i] = np.mean(sim.vel[I,0])
+ xsum_space_avg[i] = np.mean(sim.xyzsum[I,0])
+
+ plt.plot(sim.vel[:,0], sim.x[:,2], '.',
+ color=scatter_color, alpha=scatter_alpha)
+ plt.plot(v_x_space_avg, z, '+-k')
+ plt.title(title)
+ plt.xlabel('Horizontal particle velocity $v_x$ [m/s]')
+ plt.ylabel('Vertical position $z$ [m]')
+ plt.savefig(sim.id() + '-vel.pdf')
+ plt.savefig(sim.id() + '-vel.png')
+ plt.clf()
+
+ plt.plot(sim.xyzsum[:,0]/sim.time_current, sim.x[:,2], '.',
+ color=scatter_color, alpha=scatter_alpha)
+ plt.plot(xsum_space_avg/sim.time_current, z, '+-k')
+ plt.title(title)
+ plt.xlabel('Average horizontal particle velocity $\\bar{v}_x$ [m/s]')
+ plt.ylabel('Vertical position $z$ [m]')
+ plt.savefig(sim.id() + '-avg_vel.pdf')
+ plt.savefig(sim.id() + '-avg_vel.png')
+ plt.clf()
+
+ plt.close()
+ outfiles += sim.id() + '-avg_vel.png '
+
+ dp_dz_vals[j] = dpdz
+ slope_angles[j] = slope_angle
+ fluxes[j] = np.trapz(xsum_space_avg/sim.time_current, dx=dz)
+ j += 1
+
+ subprocess.call('montage ' + outfiles +
+ '-geometry +0+0 ' +
+ 'supraglacial-avg_vel-dpdz_' + str(dpdz) + '.png ',
+ shell=True)
+
+for dpdz in dpdz_values:
+ I = np.nonzero(dp_dz_vals == dpdz)
+ plt.semilogy(slope_angles[I[0]], fluxes[I[0]], '+-', label='$dp/dz$ = ' + str(dpdz) + ' Pa/m')
+
+plt.legend()
+plt.xlabel('Slope [$^\circ$]')
+plt.ylabel('Specific sediment flux [m$^2$/s]')
+plt.savefig('supraglacial_flux.png')
+plt.savefig('supraglacial_flux.pdf')
+'''
+
+# time series
+for dpdz in dpdz_values:
+ for slope_angle in slope_angle_values:
+ sim = sphere.sim("supraglacial-slope{}-dpdz{}".format(slope_angle, dpdz), fluid=True)
+ print('### ' + sim.id())
+ sim.readlast()
+
+ N_time = min(100, sim.status())
+ timesteps = np.linspace(sim.time_file_dt[0], sim.time_current[0] - sim.time_file_dt[0], N_time)
+ porosity = np.empty(N_time)
+ velocity = np.empty(N_time)
+ displacement = np.empty(N_time)
+ flux = np.empty(N_time)
+ z = np.zeros(20)
+ v_x_space_avg = np.empty_like(z)
+ xsum_space_avg = np.empty_like(z)
+
+ for it in np.arange(N_time):
+ sim.readTime(timesteps[it])
+
+ dz = np.max(sim.x[:,2])/len(z)
+ for i in range(len(z)):
+ z[i] = i*dz + 0.5*dz
+ I = np.nonzero((sim.x[:,2] >= i*dz) & (sim.x[:,2] < (i+1)*dz))
+ xsum_space_avg[i] = np.mean(sim.xyzsum[I,0])
+
+ porosity[it] = np.mean(sim.phi)
+ velocity[it] = np.mean(sim.vel[:,0])
+ displacement[it] = np.mean(sim.xyzsum[:,0])
+ flux[it] = np.trapz(xsum_space_avg/sim.time_current, dx=dz)
+
+ fig = plt.figure(figsize=(6,8))
+ fig.suptitle('slope = ' + str(slope_angle) + '$^\circ$, ' + \
+ '$dp/dz$ = -' + str(dpdz) + ' Pa/m',
+ horizontalalignment='left')
+
+ ax1 = plt.subplot(3,1,1)
+ plt.plot(timesteps, porosity, '-')
+ plt.ylabel('Porosity [-]')
+ plt.setp(ax1.get_xticklabels(), visible=False)
+
+ #ax2 = plt.subplot(3,1,2)
+ #plt.semilogy(timesteps, velocity, '-')
+ #plt.ylabel('Avg. velocity [m/s]')
+ #plt.setp(ax2.get_xticklabels(), visible=False)
+
+ ax3 = plt.subplot(3,1,2)
+ plt.semilogy(timesteps, displacement, '-')
+ plt.ylabel('Avg. displacement [m]')
+ plt.setp(ax3.get_xticklabels(), visible=False)
+
+ ax1 = plt.subplot(3,1,3)
+ plt.semilogy(timesteps, flux, '-')
+ plt.ylabel('Cumulative flux [m$^2$/s]')
+ plt.xlabel('Time [s]')
+
+ plt.tight_layout()
+ plt.savefig(sim.id() + '-timeseries.png')
+ plt.savefig(sim.id() + '-timeseries.pdf')
+ plt.clf()
+ plt.close()
DIR diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
t@@ -19,13 +19,15 @@ INCLUDE(FindCUDA)
# NOTE: Multiple arguments must be semi-colon selimited
IF (GPU_GENERATION EQUAL 1) # Kepler
SET(CUDA_NVCC_FLAGS
- "--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\";--fmad=false -ccbin gcc-8.3")
+ "--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\";--fmad=false -ccbin gcc")
+ #"--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\";--fmad=false -ccbin gcc-8.3")
#"--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\";--fmad=false -ccbin gcc-4.6")
#"--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\";--fmad=false -ccbin gcc;-Xcompiler -fPIC")
#"--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\";--fmad=false;-ccbin clang-3.8")
ELSE() # Fermi
SET(CUDA_NVCC_FLAGS
- "--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20\";--fmad=false -ccbin gcc-8.3")
+ "--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20\";--fmad=false -ccbin gcc")
+ #"--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20\";--fmad=false -ccbin gcc-8.3")
#"--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20\";--fmad=false -ccbin gcc;-Xcompiler -fPIC")
#"--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20\";--fmad=false;-ccbin clang-3.8")
ENDIF (GPU_GENERATION EQUAL 1)