Granular.jl

Julia package for granular dynamics simulation
git clone git://src.adamsgaard.dk/Granular.jl
Log | Files | Refs | README | LICENSE

commit 4a79e477b83eb068b4024f6695bca0295a5a06c4
parent eeb31284d08583ecaa038ef8fb55d77692c48e71
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Sun, 30 Apr 2017 22:14:46 -0400

change example to include in namespace, add warning if point is outside cell during ocean drag

Diffstat:
Mexamples/nares_strait.jl | 43+++++++++++++++++++++++++------------------
Msrc/grid.jl | 2++
Msrc/io.jl | 6+++---
Msrc/ocean.jl | 7+++++++
Mtest/vtk.jl | 2+-
5 files changed, 38 insertions(+), 22 deletions(-)

diff --git a/examples/nares_strait.jl b/examples/nares_strait.jl @@ -1,15 +1,16 @@ #!/usr/bin/env julia -using SeaIce +import SeaIce -sim = createSimulation(id="nares_strait") +sim = SeaIce.createSimulation(id="nares_strait") # Initialize ocean Lx = 50.e3 Lx_constriction = Lx*.25 L = [Lx, Lx*1.5, 1e3] Ly_constriction = L[2]*.33 -n = [100, 100, 2] -sim.ocean = createRegularOceanGrid(n, L, name="poiseuille_flow") +#n = [100, 100, 2] +n = [50, 50, 2] +sim.ocean = SeaIce.createRegularOceanGrid(n, L, name="poiseuille_flow") sim.ocean.v[:, :, 1, 1] = 1e-8*((sim.ocean.xq - Lx/2.).^2 - Lx^2./4.) # Initialize confining walls, which are ice floes that are fixed in space @@ -21,14 +22,15 @@ h = 1. for y in linspace((L[2] - Ly_constriction)/2., Ly_constriction + (L[2] - Ly_constriction)/2., Int(floor(Ly_constriction/(r*2)))) - addIceFloeCylindrical(sim, [(Lx - Lx_constriction)/2., y], r, h, fixed=true, -verbose=false) + SeaIce.addIceFloeCylindrical(sim, [(Lx - Lx_constriction)/2., y], r, h, + fixed=true, verbose=false) end for y in linspace((L[2] - Ly_constriction)/2., Ly_constriction + (L[2] - Ly_constriction)/2., Int(floor(Ly_constriction/(r*2)))) - addIceFloeCylindrical(sim, [Lx_constriction + (L[1] - Lx_constriction)/2., - y], r, h, fixed=true, verbose=false) + SeaIce.addIceFloeCylindrical(sim, + [Lx_constriction + (L[1] - Lx_constriction)/2., + y], r, h, fixed=true, verbose=false) end dx = 2.*r*sin(atan((Lx - Lx_constriction)/(L[2] - Ly_constriction))) @@ -38,7 +40,8 @@ x = r:dx:((Lx - Lx_constriction)/2.) y = linspace(L[2] - r, (L[2] - Ly_constriction)/2. + Ly_constriction + r, length(x)) for i in 1:length(x) - addIceFloeCylindrical(sim, [x[i], y[i]], r, h, fixed=true, verbose=false) + SeaIce.addIceFloeCylindrical(sim, [x[i], y[i]], r, h, fixed=true, + verbose=false) end ## NE diagonal @@ -46,21 +49,24 @@ x = (L[1] - r):(-dx):((Lx - Lx_constriction)/2. + Lx_constriction) y = linspace(L[2] - r, (L[2] - Ly_constriction)/2. + Ly_constriction + r, length(x)) for i in 1:length(x) - addIceFloeCylindrical(sim, [x[i], y[i]], r, h, fixed=true, verbose=false) + SeaIce.addIceFloeCylindrical(sim, [x[i], y[i]], r, h, fixed=true, + verbose=false) end ## SW diagonal x = r:dx:((Lx - Lx_constriction)/2.) y = linspace(r, (L[2] - Ly_constriction)/2. - r, length(x)) for i in 1:length(x) - addIceFloeCylindrical(sim, [x[i], y[i]], r, h, fixed=true, verbose=false) + SeaIce.addIceFloeCylindrical(sim, [x[i], y[i]], r, h, fixed=true, + verbose=false) end ## SE diagonal x = (L[1] - r):(-dx):((Lx - Lx_constriction)/2. + Lx_constriction) y = linspace(r, (L[2] - Ly_constriction)/2. - r, length(x)) for i in 1:length(x) - addIceFloeCylindrical(sim, [x[i], y[i]], r, h, fixed=true, verbose=false) + SeaIce.addIceFloeCylindrical(sim, [x[i], y[i]], r, h, fixed=true, + verbose=false) end n_walls = length(sim.ice_floes) @@ -81,7 +87,7 @@ for y in (L[2] - r - noise_amplitude):(-(2.*r + floe_padding)):((L[2] - end x_ = x + noise_amplitude*(0.5 - Base.Random.rand()) - y_ = y +noise_amplitude*(0.5 - Base.Random.rand()) + y_ = y + noise_amplitude*(0.5 - Base.Random.rand()) if y_ < -dy/dx*x_ + L[2] + spacing_to_boundaries continue @@ -91,16 +97,17 @@ for y in (L[2] - r - noise_amplitude):(-(2.*r + floe_padding)):((L[2] - continue end - addIceFloeCylindrical(sim, [x_, y_], r, h, verbose=false) + SeaIce.addIceFloeCylindrical(sim, [x_, y_], r, h, verbose=false) end iy += 1 end n = length(sim.ice_floes) - n_walls info("added $(n) ice floes") -writeVTK(sim) +SeaIce.writeVTK(sim) # Run temporal loop -setTotalTime!(sim, 24.*60.*60.) -setTimeStep!(sim) -run!(sim) +SeaIce.setTotalTime!(sim, 24.*60.*60.) +SeaIce.setOutputFileInterval!(sim, 60.) +SeaIce.setTimeStep!(sim) +SeaIce.run!(sim, status_interval=1) diff --git a/src/grid.jl b/src/grid.jl @@ -1,4 +1,6 @@ """ + bilinearInterpolation(field, x_tilde, y_tilde, i, j, k, it) + Use bilinear interpolation to interpolate a staggered grid to an arbitrary position in a cell. Assumes south-west convention, i.e. (i,j) is located at the south-west (-x, -y)-facing corner. diff --git a/src/io.jl b/src/io.jl @@ -15,7 +15,7 @@ written to separate `.vtu` files. This can be disabled by setting the argument """ function writeVTK(simulation::Simulation; folder::String=".", - verbose::Bool=false, + verbose::Bool=true, ocean::Bool=true) simulation.file_number += 1 @@ -97,7 +97,7 @@ function writeIceFloeVTK(simulation::Simulation, outfiles = WriteVTK.vtk_save(vtkfile) if verbose - println("Output file: " * outfiles[1]) + info("Output file: " * outfiles[1]) else return nothing end @@ -151,7 +151,7 @@ function writeOceanVTK(ocean::Ocean, outfiles = WriteVTK.vtk_save(vtkfile) if verbose - println("Output file: " * outfiles[1]) + info("Output file: " * outfiles[1]) else return nothing end diff --git a/src/ocean.jl b/src/ocean.jl @@ -261,6 +261,13 @@ function addOceanDrag!(simulation::Simulation) x_tilde, y_tilde = getNonDimensionalCellCoordinates(simulation.ocean, i, j, ice_floe.lin_pos) + if x_tilde < 0. || x_tilde > 1. || y_tilde < 0. || y_tilde > 1. + warn(""" + relative coordinates outside bounds ($(x_tilde), $(y_tilde)), + pos = $(ice_floe.lin_pos) at i,j = $(i), $(j). + + """) + end u_local = bilinearInterpolation(u, x_tilde, y_tilde, i, j, k, 1) v_local = bilinearInterpolation(v, x_tilde, y_tilde, i, j, k, 1) diff --git a/test/vtk.jl b/test/vtk.jl @@ -9,7 +9,7 @@ sim = SeaIce.createSimulation(id="test") SeaIce.addIceFloeCylindrical(sim, [ 0., 0.], 10., 1., verbose=false) SeaIce.addIceFloeCylindrical(sim, [18., 0.], 10., 1., verbose=false) sim.ocean = SeaIce.createRegularOceanGrid([10, 20, 5], [10., 25., 2.]) -SeaIce.writeVTK(sim) +SeaIce.writeVTK(sim, verbose=false) if Base.is_linux() cmd = "sha256sum"