double_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)