tadd example where packing is initialized from a PNG image. Use shorthand for random functions - Granular.jl - Julia package for granular dynamics simulation
(HTM) git clone git://src.adamsgaard.dk/Granular.jl
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
(DIR) commit e5fe4e54a81b423af217a9eaa8ddd442a6f3ba4a
(DIR) parent 82f95ef1f309fbbb0bbdfc40fa99eaf1219e7bf2
(HTM) Author: Anders Damsgaard <andersd@riseup.net>
Date: Fri, 22 Sep 2017 09:55:05 -0400
add example where packing is initialized from a PNG image. Use shorthand for random functions
Diffstat:
M examples/double_gyre.jl | 6 +++---
A examples/image.jl | 137 +++++++++++++++++++++++++++++++
M examples/logo.jl | 8 ++++----
M examples/strait.jl | 20 ++++++++++----------
M src/grid.jl | 6 +++---
M src/util.jl | 2 +-
M test/util.jl | 2 +-
7 files changed, 159 insertions(+), 22 deletions(-)
---
(DIR) diff --git a/examples/double_gyre.jl b/examples/double_gyre.jl
t@@ -57,7 +57,7 @@ info("added $(n_walls) fixed ice floes as walls")
# Initialize ice floes everywhere
floe_padding = .5*r
noise_amplitude = .8*floe_padding
-Base.Random.srand(1)
+srand(1)
for y in (4.*r + noise_amplitude):(2.*r + floe_padding):(L[2] - 4.*r -
noise_amplitude)
t@@ -67,8 +67,8 @@ for y in (4.*r + noise_amplitude):(2.*r + floe_padding):(L[2] - 4.*r -
#x += 1.5*r
#end
- x_ = x + noise_amplitude*(0.5 - Base.Random.rand())
- y_ = y + noise_amplitude*(0.5 - Base.Random.rand())
+ x_ = x + noise_amplitude*(0.5 - rand())
+ y_ = y + noise_amplitude*(0.5 - rand())
SeaIce.addIceFloeCylindrical!(sim, [x_, y_], r, h, verbose=false)
end
(DIR) diff --git a/examples/image.jl b/examples/image.jl
t@@ -0,0 +1,137 @@
+#!/usr/bin/env julia
+
+import SeaIce
+import FileIO
+import Colors
+
+const verbose = true
+
+const img_file = "Cundall2008.png"
+
+img = FileIO.load(img_file)
+
+# resize the image if it is too large, preceed with lopass to avoid antialias
+max_pixels = 100^2
+if size(img, 1)*size(img, 2) > max_pixels
+ cp(img_file, "backup-" * img_file, remove_destination=true)
+ run(`convert $(img_file) -resize "$(max_pixels)@>" $(img_file)`)
+ img = FileIO.load(img_file)
+end
+
+const img_bw = Colors.Gray.(img)
+
+#const forcing = "gyres"
+#const forcing = "down"
+#const forcing = "convergent"
+const forcing = "sandpile"
+
+const dx = 1.
+const dy = dx
+
+const nx = size(img_bw, 2) + 1
+const ny = size(img_bw, 1) + 1
+
+Lx = nx*dx
+if forcing == "sandpile"
+ Lx *= 1.5
+end
+const Ly = ny*dy
+
+const youngs_modulus = 2e7
+const tensile_strength = 0e3
+const h = .5
+
+sim = SeaIce.createSimulation(id="image")
+
+info("nx = $nx, ny = $ny")
+
+for iy=1:size(img_bw, 1)
+ for ix=1:size(img_bw, 2)
+
+ x = ix*dx - dx
+ if forcing == "sandpile"
+ x += Lx/6.
+ end
+ y = Ly - (iy*dy - dy)
+ r = .5*dx*((1. - Float64(img_bw[iy, ix])))
+
+ if r > .1*dx
+ SeaIce.addIceFloeCylindrical!(sim, [x + dx, y - dy], r, h,
+ tensile_strength=tensile_strength,
+ youngs_modulus=youngs_modulus,
+ verbose=verbose)
+ end
+ end
+end
+
+# set ocean forcing
+sim.ocean = SeaIce.createRegularOceanGrid([nx, ny, 1], [Lx, Ly, 1.],
+ name="image_ocean")
+
+if forcing == "gyres"
+ epsilon = 0.25 # amplitude of periodic oscillations
+ t = 0.
+ a = epsilon*sin(2.*pi*t)
+ b = 1. - 2.*epsilon*sin(2.*pi*t)
+ for i=1:size(sim.ocean.u, 1)
+ for j=1:size(sim.ocean.u, 2)
+
+ x = sim.ocean.xq[i, j]/(Lx*.5) # x in [0;2]
+ y = sim.ocean.yq[i, j]/Ly # y in [0;1]
+
+ f = a*x^2. + b*x
+ df_dx = 2.*a*x + b
+
+ sim.ocean.u[i, j, 1, 1] = -pi/10.*sin(pi*f)*cos(pi*y) * 4e1
+ sim.ocean.v[i, j, 1, 1] = pi/10.*cos(pi*f)*sin(pi*y)*df_dx * 4e1
+ end
+ end
+
+elseif forcing == "down" || forcing == "sandpile"
+ srand(1)
+ sim.ocean.u[:, :, 1, 1] = (rand(nx+1, ny+1) - .5)*.1
+ sim.ocean.v[:, :, 1, 1] = -Ly/5.
+
+elseif forcing == "convergent"
+ srand(1)
+ sim.ocean.u[:, :, 1, 1] = (rand(nx+1, ny+1) - .5)*.1
+ for j=1:size(sim.ocean.u, 2)
+ sim.ocean.v[:, j, 1, 1] = -(j/ny - .5)*10.
+ end
+
+else
+ error("Forcing not understood")
+end
+
+# Initialize confining walls, which are ice floes that are fixed in space
+r = dx/4.
+
+## N-S wall segments
+for y in linspace(r, Ly-r, Int(round((Ly - 2.*r)/(r*2))))
+ SeaIce.addIceFloeCylindrical!(sim, [r, y], r, h, fixed=true,
+ youngs_modulus=youngs_modulus,
+ verbose=false)
+ SeaIce.addIceFloeCylindrical!(sim, [Lx-r, y], r, h, fixed=true,
+ youngs_modulus=youngs_modulus,
+ verbose=false)
+end
+
+## E-W wall segments
+for x in linspace(3.*r, Lx-3.*r, Int(round((Lx - 6.*r)/(r*2))))
+ SeaIce.addIceFloeCylindrical!(sim, [x, r], r, h, fixed=true,
+ youngs_modulus=youngs_modulus,
+ verbose=false)
+ SeaIce.addIceFloeCylindrical!(sim, [x, Ly-r], r, h, fixed=true,
+ youngs_modulus=youngs_modulus,
+ verbose=false)
+end
+
+
+# Finalize setup and start simulation
+SeaIce.setTimeStep!(sim, verbose=verbose)
+
+SeaIce.setTotalTime!(sim, 5.)
+SeaIce.setOutputFileInterval!(sim, .1)
+
+SeaIce.removeSimulationFiles(sim)
+SeaIce.run!(sim, verbose=verbose)
(DIR) diff --git a/examples/logo.jl b/examples/logo.jl
t@@ -118,13 +118,13 @@ if forcing == "gyres"
end
elseif forcing == "down"
- Base.Random.srand(1)
- sim.ocean.u[:, :, 1, 1] = (Base.Random.rand(nx+1, ny+1) - .5)*.1
+ srand(1)
+ sim.ocean.u[:, :, 1, 1] = (rand(nx+1, ny+1) - .5)*.1
sim.ocean.v[:, :, 1, 1] = -5.
elseif forcing == "convergent"
- Base.Random.srand(1)
- sim.ocean.u[:, :, 1, 1] = (Base.Random.rand(nx+1, ny+1) - .5)*.1
+ srand(1)
+ sim.ocean.u[:, :, 1, 1] = (rand(nx+1, ny+1) - .5)*.1
for j=1:size(sim.ocean.u, 2)
sim.ocean.v[:, j, 1, 1] = -(j/ny - .5)*10.
end
(DIR) diff --git a/examples/strait.jl b/examples/strait.jl
t@@ -82,7 +82,7 @@ dy = sqrt((2.*r_walls)^2. - dx^2.)
spacing_to_boundaries = 4.*r
floe_padding = .5*r
noise_amplitude = floe_padding
-Base.Random.srand(1)
+srand(1)
for y in (L[2] - r - noise_amplitude):(-(2.*r + floe_padding)):((L[2] -
Ly_constriction)/2. + Ly_constriction)
for x in (r + noise_amplitude):(2.*r + floe_padding):(Lx - r -
t@@ -91,8 +91,8 @@ for y in (L[2] - r - noise_amplitude):(-(2.*r + floe_padding)):((L[2] -
x += 1.5*r
end
- x_ = x + noise_amplitude*(0.5 - Base.Random.rand())
- y_ = y + noise_amplitude*(0.5 - Base.Random.rand())
+ x_ = x + noise_amplitude*(0.5 - rand())
+ y_ = y + noise_amplitude*(0.5 - rand())
if y_ < -dy/dx*x_ + L[2] + spacing_to_boundaries
continue
t@@ -102,7 +102,7 @@ for y in (L[2] - r - noise_amplitude):(-(2.*r + floe_padding)):((L[2] -
continue
end
- r_rand = r_min + Base.Random.rand()*(r - r_min)
+ r_rand = r_min + rand()*(r - r_min)
SeaIce.addIceFloeCylindrical!(sim, [x_, y_], r_rand, h, verbose=false)
end
iy += 1
t@@ -144,7 +144,7 @@ while sim.time < sim.time_total
size(sim.ocean.xh, 2))
# Enable for high mass flux
- r_rand = r_min + Base.Random.rand()*(r - r_min)
+ r_rand = r_min + rand()*(r - r_min)
SeaIce.addIceFloeCylindrical!(sim, [x-r, y-r], r_rand, h,
verbose=false,
contact_stiffness_normal=k_n,
t@@ -152,7 +152,7 @@ while sim.time < sim.time_total
contact_viscosity_tangential=gamma_t,
contact_dynamic_friction = mu_d,
rotating=rotating)
- r_rand = r_min + Base.Random.rand()*(r - r_min)
+ r_rand = r_min + rand()*(r - r_min)
SeaIce.addIceFloeCylindrical!(sim, [x+r, y-r], r_rand, h,
verbose=false,
contact_stiffness_normal=k_n,
t@@ -160,7 +160,7 @@ while sim.time < sim.time_total
contact_viscosity_tangential=gamma_t,
contact_dynamic_friction = mu_d,
rotating=rotating)
- r_rand = r_min + Base.Random.rand()*(r - r_min)
+ r_rand = r_min + rand()*(r - r_min)
SeaIce.addIceFloeCylindrical!(sim, [x+r, y+r], r_rand, h,
verbose=false,
contact_stiffness_normal=k_n,
t@@ -168,7 +168,7 @@ while sim.time < sim.time_total
contact_viscosity_tangential=gamma_t,
contact_dynamic_friction = mu_d,
rotating=rotating)
- r_rand = r_min + Base.Random.rand()*(r - r_min)
+ r_rand = r_min + rand()*(r - r_min)
SeaIce.addIceFloeCylindrical!(sim, [x-r, y+r], r_rand, h,
verbose=false,
contact_stiffness_normal=k_n,
t@@ -178,8 +178,8 @@ while sim.time < sim.time_total
rotating=rotating)
# Enable for low mass flux
- #x += noise_amplitude*(0.5 - Base.Random.rand())
- #y += noise_amplitude*(0.5 - Base.Random.rand())
+ #x += noise_amplitude*(0.5 - rand())
+ #y += noise_amplitude*(0.5 - rand())
#SeaIce.addIceFloeCylindrical!(sim, [x, y], r, h, verbose=false)
end
end
(DIR) diff --git a/src/grid.jl b/src/grid.jl
t@@ -484,10 +484,10 @@ function findEmptyPositionInGridCell(simulation::Simulation,
for i_iter=1:n_iter
overlap_found = false
- Base.Random.srand(i*j*seed*i_iter)
+ srand(i*j*seed*i_iter)
# generate random candidate position
- x_tilde = Base.Random.rand()
- y_tilde = Base.Random.rand()
+ x_tilde = rand()
+ y_tilde = rand()
bilinearInterpolation!(pos, grid.xq, grid.yq, x_tilde, y_tilde, i, j)
if verbose
info("trying poisition $pos in cell $i,$j")
(DIR) diff --git a/src/util.jl b/src/util.jl
t@@ -18,7 +18,7 @@ Returns one or more random numbers from a power-law probability distribution.
max_val::Number = 1.)
val = ((max_val^(distribution_power + 1.) -
- min_val^(distribution_power + 1.))*Base.Random.rand(dims) +
+ min_val^(distribution_power + 1.))*rand(dims) +
min_val^(distribution_power + 1.)).^(1./(distribution_power + 1.))
if dims == 1
(DIR) diff --git a/test/util.jl b/test/util.jl
t@@ -15,7 +15,7 @@ info("Testing power-law RNG")
@test 5 == length(SeaIce.randpower(5))
@test (5,) == size(SeaIce.randpower(5))
-Base.Random.srand(1)
+srand(1)
for i=1:10^5
@test 0. <= SeaIce.randpower() <= 1.
@test 0. <= SeaIce.randpower(1, 1., 0., 1.) <= 1.