tdouble_gyre.jl - 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
       ---
       tdouble_gyre.jl (3099B)
       ---
            1 #!/usr/bin/env julia
            2 import Granular
            3 using Random
            4 
            5 sim = Granular.createSimulation(id="double_gyre")
            6 
            7 # Initialize ocean
            8 L = [100e3, 50e3, 1e3]
            9 Ly_constriction = 20e3
           10 #n = [750, 500, 2]  # high resolution
           11 n = [30, 15, 2]  # intermedite resolution
           12 #n = [8, 5, 2]  # coarse resolution
           13 sim.ocean = Granular.createRegularOceanGrid(n, L, name="double_gyre")
           14 
           15 epsilon = 0.25  # amplitude of periodic oscillations
           16 t = 0.
           17 a = epsilon*sin(2.0*pi*t)
           18 b = 1.0 - 2.0*epsilon*sin(2.0*pi*t)
           19 for i=1:size(sim.ocean.u, 1)
           20     for j=1:size(sim.ocean.u, 2)
           21 
           22         x = sim.ocean.xq[i, j]/(L[1]*0.5)  # x in [0;2]
           23         y = sim.ocean.yq[i, j]/L[2]       # y in [0;1]
           24 
           25         f = a*x^2.0 + b*x
           26         df_dx = 2.0*a*x + b
           27 
           28         sim.ocean.u[i, j, 1, 1] = -pi/10.0*sin(pi*f)*cos(pi*y) * 1e1
           29         sim.ocean.v[i, j, 1, 1] = pi/10.0*cos(pi*f)*sin(pi*y)*df_dx * 1e1
           30     end
           31 end
           32 
           33 # Initialize confining walls, which are ice floes that are fixed in space
           34 r = minimum(L[1:2]./n[1:2])/2.0
           35 h = 1.
           36 
           37 ## N-S wall segments
           38 for y in range(r, stop=L[2]-r, length=Int(round((L[2] - 2.0*r)/(r*2))))
           39     Granular.addGrainCylindrical!(sim, [r, y], r, h, fixed=true,
           40                                   verbose=false)
           41     Granular.addGrainCylindrical!(sim, [L[1]-r, y], r, h, fixed=true,
           42                                   verbose=false)
           43 end
           44 
           45 ## E-W wall segments
           46 for x in range(3.0*r, stop=L[1]-3.0*r,
           47                length=Int(round((L[1] - 6.0*r)/(r*2))))
           48     Granular.addGrainCylindrical!(sim, [x, r], r, h, fixed=true,
           49                                   verbose=false)
           50     Granular.addGrainCylindrical!(sim, [x, L[2]-r], r, h, fixed=true,
           51                                   verbose=false)
           52 end
           53 
           54 n_walls = length(sim.grains)
           55 @info "added $(n_walls) fixed ice floes as walls"
           56 
           57 
           58 
           59 # Initialize ice floes everywhere
           60 floe_padding = 0.5*r
           61 noise_amplitude = 0.8*floe_padding
           62 Random.seed!(1)
           63 for y in (4.0*r + noise_amplitude):(2.0*r + floe_padding):(L[2] - 4.0*r - 
           64                                                            noise_amplitude)
           65                                                          
           66     for x in (4.0*r + noise_amplitude):(2.0*r + floe_padding):(L[1] - 4.0*r - 
           67                                                                noise_amplitude)
           68         #if iy % 2 == 0
           69             #x += 1.5*r
           70         #end
           71 
           72         x_ = x + noise_amplitude*(0.5 - rand())
           73         y_ = y + noise_amplitude*(0.5 - rand())
           74 
           75         Granular.addGrainCylindrical!(sim, [x_, y_], r, h, verbose=false)
           76     end
           77 end
           78 n = length(sim.grains) - n_walls
           79 @info "added $n ice floes"
           80 
           81 # Remove old simulation files
           82 Granular.removeSimulationFiles(sim)
           83 
           84 k_n = 1e6  # N/m
           85 gamma_t = 1e7  # N/(m/s)
           86 mu_d = 0.7
           87 rotating = false
           88 for i=1:length(sim.grains)
           89     sim.grains[i].contact_stiffness_normal = k_n
           90     sim.grains[i].contact_stiffness_tangential = k_n
           91     sim.grains[i].contact_viscosity_tangential = gamma_t
           92     sim.grains[i].contact_dynamic_friction = mu_d
           93     sim.grains[i].rotating = rotating
           94 end
           95 
           96 # Set temporal parameters
           97 Granular.setTotalTime!(sim, 12.0*60.0*60.0)
           98 Granular.setOutputFileInterval!(sim, 60.0)
           99 Granular.setTimeStep!(sim)
          100 
          101 Granular.run!(sim)