tMerge branch 'master' of src.adamsgaard.dk:src/sphere - sphere - GPU-based 3D … | |
git clone git://src.adamsgaard.dk/sphere | |
Log | |
Files | |
Refs | |
LICENSE | |
--- | |
commit 4ad7f4d8bd40be3988fec86661de48699a1ca69b | |
parent a944c7e38af5c26553f8add9df01b56fe98c4c01 | |
Author: Anders Damsgaard <[email protected]> | |
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(-) | |
--- | |
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… | |
+ | |
+ # 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) | |
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() | |
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, dpd… | |
+ 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$ = ' + … | |
+ | |
+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, dpd… | |
+ print('### ' + sim.id()) | |
+ sim.readlast() | |
+ | |
+ N_time = min(100, sim.status()) | |
+ timesteps = np.linspace(sim.time_file_dt[0], sim.time_current[0] - sim… | |
+ 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() | |
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\… | |
+ "--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\… | |
+ #"--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35… | |
#"--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35… | |
#"--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35… | |
#"--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35… | |
ELSE() # Fermi | |
SET(CUDA_NVCC_FLAGS | |
- "--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20\… | |
+ "--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20\… | |
+ #"--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20… | |
#"--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20… | |
#"--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20… | |
ENDIF (GPU_GENERATION EQUAL 1) |