Granular.jl

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

commit 0440db5e1befd66a7f4c9fb98061d8f4604d56f1
parent 784fc626abb1dde66a6aed241677ad8d5d8b53ac
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Mon,  8 May 2017 16:58:56 -0400

Merge branch 'master' of github.com:anders-dc/SeaIce.jl

Diffstat:
Aexamples/double_gyre.jl | 100+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Mexamples/nares_strait.jl | 51++++++++++++++++++++++++++++++++++++++-------------
Mtest/collision-2floes-oblique.jl | 28++++++++++++++--------------
3 files changed, 152 insertions(+), 27 deletions(-)

diff --git a/examples/double_gyre.jl b/examples/double_gyre.jl @@ -0,0 +1,100 @@ +#!/usr/bin/env julia +import SeaIce + +sim = SeaIce.createSimulation(id="double_gyre") + +# Initialize ocean +L = [100e3, 50e3, 1e3] +Ly_constriction = 20e3 +#n = [750, 500, 2] # high resolution +n = [30, 15, 2] # intermedite resolution +#n = [8, 5, 2] # coarse resolution +sim.ocean = SeaIce.createRegularOceanGrid(n, L, name="double_gyre") + +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]/(L[1]*.5) # x in [0;2] + y = sim.ocean.yq[i, j]/L[2] # 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) * 1e1 + sim.ocean.v[i, j, 1, 1] = pi/10.*cos(pi*f)*sin(pi*y)*df_dx * 1e1 + end +end + +# Initialize confining walls, which are ice floes that are fixed in space +r = minimum(L[1:2]/n[1:2])/2. +h = 1. + +## N-S wall segments +for y in linspace(r, L[2]-r, Int(round((L[2] - 2.*r)/(r*2)))) + SeaIce.addIceFloeCylindrical(sim, [r, y], r, h, fixed=true, + verbose=false) + SeaIce.addIceFloeCylindrical(sim, [L[1]-r, y], r, h, fixed=true, + verbose=false) +end + +## E-W wall segments +for x in linspace(3.*r, L[1]-3.*r, Int(round((L[1] - 6.*r)/(r*2)))) + SeaIce.addIceFloeCylindrical(sim, [x, r], r, h, fixed=true, + verbose=false) + SeaIce.addIceFloeCylindrical(sim, [x, L[2]-r], r, h, fixed=true, + verbose=false) +end + +n_walls = length(sim.ice_floes) +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) +for y in (4.*r + noise_amplitude):(2.*r + floe_padding):(L[2] - 4.*r - + noise_amplitude) + + for x in (4.*r + noise_amplitude):(2.*r + floe_padding):(L[1] - 4.*r - + noise_amplitude) + #if iy % 2 == 0 + #x += 1.5*r + #end + + x_ = x + noise_amplitude*(0.5 - Base.Random.rand()) + y_ = y + noise_amplitude*(0.5 - Base.Random.rand()) + + SeaIce.addIceFloeCylindrical(sim, [x_, y_], r, h, verbose=false) + end +end +n = length(sim.ice_floes) - n_walls +info("added $(n) ice floes") + +# Remove old simulation files +SeaIce.removeSimulationFiles(sim) + +k_n = 1e6 # N/m +gamma_t = 1e7 # N/(m/s) +mu_d = 0.7 +rotating = false +for i=1:length(sim.ice_floes) + sim.ice_floes[i].contact_stiffness_normal = k_n + sim.ice_floes[i].contact_stiffness_tangential = k_n + sim.ice_floes[i].contact_viscosity_tangential = gamma_t + sim.ice_floes[i].contact_dynamic_friction = mu_d + sim.ice_floes[i].rotating = rotating +end + +# Set temporal parameters +SeaIce.setTotalTime!(sim, 12.*60.*60.) +SeaIce.setOutputFileInterval!(sim, 60.) +SeaIce.setTimeStep!(sim) + +SeaIce.run!(sim, status_interval=1, + contact_tangential_rheology="Linear Viscous Frictional") diff --git a/examples/nares_strait.jl b/examples/nares_strait.jl @@ -5,12 +5,13 @@ sim = SeaIce.createSimulation(id="nares_strait") # Initialize ocean Lx = 50.e3 -Lx_constriction = Lx*.25 +Lx_constriction = 10e3 L = [Lx, Lx*1.5, 1e3] -Ly_constriction = L[2]*.33 +Ly_constriction = 20e3 #n = [100, 100, 2] # high resolution #n = [50, 50, 2] # intermedite resolution -n = [6, 6, 2] # coarse resolution +n = [25, 25, 2] +#n = [6, 6, 2] # coarse resolution 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.) @@ -108,16 +109,24 @@ info("added $(n) ice floes") # Remove old simulation files SeaIce.removeSimulationFiles(sim) -# Set temporal parameters -SeaIce.setTotalTime!(sim, 24.*60.*60.) -SeaIce.setOutputFileInterval!(sim, 60.) -SeaIce.setTimeStep!(sim) - -gamma_t = 1e4 # N/(m/s) +k_n = 1e6 # N/m +k_t = k_n +gamma_t = 1e7 # N/(m/s) +mu_d = 0.7 +rotating = true for i=1:length(sim.ice_floes) + sim.ice_floes[i].contact_stiffness_normal = k_n + sim.ice_floes[i].contact_stiffness_tangential = k_t sim.ice_floes[i].contact_viscosity_tangential = gamma_t + sim.ice_floes[i].contact_dynamic_friction = mu_d + sim.ice_floes[i].rotating = rotating end +# Set temporal parameters +SeaIce.setTotalTime!(sim, 12.*60.*60.) +SeaIce.setOutputFileInterval!(sim, 60.) +SeaIce.setTimeStep!(sim) + # Run simulation for 10 time steps, then add new icefloes from the top while sim.time < sim.time_total for it=1:10 @@ -132,13 +141,29 @@ while sim.time < sim.time_total # Enable for high mass flux SeaIce.addIceFloeCylindrical(sim, [x-r, y-r], r, h, verbose=false, - contact_viscosity_tangential=gamma_t) + contact_stiffness_normal=k_n, + contact_stiffness_normal=k_n, + contact_viscosity_tangential=gamma_t, + contact_dynamic_friction = mu_d, + rotating=rotating) SeaIce.addIceFloeCylindrical(sim, [x+r, y-r], r, h, verbose=false, - contact_viscosity_tangential=gamma_t) + contact_stiffness_normal=k_n, + contact_stiffness_normal=k_n, + contact_viscosity_tangential=gamma_t, + contact_dynamic_friction = mu_d, + rotating=rotating) SeaIce.addIceFloeCylindrical(sim, [x+r, y+r], r, h, verbose=false, - contact_viscosity_tangential=gamma_t) + contact_stiffness_normal=k_n, + contact_stiffness_normal=k_n, + contact_viscosity_tangential=gamma_t, + contact_dynamic_friction = mu_d, + rotating=rotating) SeaIce.addIceFloeCylindrical(sim, [x-r, y+r], r, h, verbose=false, - contact_viscosity_tangential=gamma_t) + contact_stiffness_normal=k_n, + contact_stiffness_normal=k_n, + contact_viscosity_tangential=gamma_t, + contact_dynamic_friction = mu_d, + rotating=rotating) # Enable for low mass flux #x += noise_amplitude*(0.5 - Base.Random.rand()) diff --git a/test/collision-2floes-oblique.jl b/test/collision-2floes-oblique.jl @@ -204,8 +204,8 @@ SeaIce.run!(sim, temporal_integration_method="Two-term Taylor", contact_tangential_rheology="Linear Viscous Frictional", verbose=verbose) -@test sim.ice_floes[1].ang_pos > 0. -@test sim.ice_floes[1].ang_vel > 0. +@test sim.ice_floes[1].ang_pos < 0. +@test sim.ice_floes[1].ang_vel < 0. @test sim.ice_floes[2].ang_pos ≈ 0. @test sim.ice_floes[2].ang_vel ≈ 0. E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) @@ -241,8 +241,8 @@ SeaIce.run!(sim, temporal_integration_method="Two-term Taylor", contact_tangential_rheology="Linear Viscous Frictional", verbose=verbose) -@test sim.ice_floes[1].ang_pos > 0. -@test sim.ice_floes[1].ang_vel > 0. +@test sim.ice_floes[1].ang_pos < 0. +@test sim.ice_floes[1].ang_vel < 0. @test sim.ice_floes[2].ang_pos ≈ 0. @test sim.ice_floes[2].ang_vel ≈ 0. E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) @@ -259,8 +259,8 @@ SeaIce.run!(sim, temporal_integration_method="Three-term Taylor", contact_tangential_rheology="Linear Viscous Frictional", verbose=verbose) -@test sim.ice_floes[1].ang_pos > 0. -@test sim.ice_floes[1].ang_vel > 0. +@test sim.ice_floes[1].ang_pos < 0. +@test sim.ice_floes[1].ang_vel < 0. @test sim.ice_floes[2].ang_pos ≈ 0. @test sim.ice_floes[2].ang_vel ≈ 0. E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) @@ -293,10 +293,10 @@ SeaIce.run!(sim, temporal_integration_method="Two-term Taylor", contact_tangential_rheology="Linear Viscous Frictional", verbose=verbose) -@test sim.ice_floes[1].ang_pos > 0. -@test sim.ice_floes[1].ang_vel > 0. -@test sim.ice_floes[2].ang_pos > 0. -@test sim.ice_floes[2].ang_vel > 0. +@test sim.ice_floes[1].ang_pos < 0. +@test sim.ice_floes[1].ang_vel < 0. +@test sim.ice_floes[2].ang_pos < 0. +@test sim.ice_floes[2].ang_vel < 0. E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim) @test_approx_eq_eps E_kin_lin_init+E_kin_rot_init E_kin_lin_final+E_kin_rot_final E_kin_lin_init*tol @@ -390,10 +390,10 @@ SeaIce.run!(sim, temporal_integration_method="Three-term Taylor", contact_tangential_rheology="Linear Viscous Frictional", verbose=verbose) -@test sim.ice_floes[1].ang_pos > 0. -@test sim.ice_floes[1].ang_vel > 0. -@test sim.ice_floes[2].ang_pos > 0. -@test sim.ice_floes[2].ang_vel > 0. +@test sim.ice_floes[1].ang_pos < 0. +@test sim.ice_floes[1].ang_vel < 0. +@test sim.ice_floes[2].ang_pos < 0. +@test sim.ice_floes[2].ang_vel < 0. E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim) @test_approx_eq_eps E_kin_lin_init+E_kin_rot_init E_kin_lin_final+E_kin_rot_final E_kin_lin_init*tol