tocean.jl - Granular.jl - Julia package for granular dynamics simulation
HTML git clone git://src.adamsgaard.dk/Granular.jl
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
tocean.jl (15284B)
---
1 import Pkg
2 using Test
3 using LinearAlgebra
4
5 hasNetCDF = false
6 if VERSION < v"0.7.0-alpha"
7 if typeof(Pkg.installed("NetCDF")) == VersionNumber
8 import NetCDF
9 hasNetCDF = true
10 end
11 else
12 import Pkg
13 if haskey(Pkg.installed(), "NetCDF")
14 import NetCDF
15 hasNetCDF = true
16 end
17 end
18 if !hasNetCDF
19 @warn "Package NetCDF not found. " *
20 "Ocean/atmosphere grid read not supported. " *
21 "If required, install NetCDF and its " *
22 "requirements with `Pkg.add(\"NetCDF\")`."
23 end
24
25 export createEmptyOcean
26 "Returns empty ocean type for initialization purposes."
27 function createEmptyOcean()
28 return Ocean(false,
29
30 zeros(1),
31
32 zeros(1,1),
33 zeros(1,1),
34 zeros(1,1),
35 zeros(1,1),
36
37 zeros(1),
38 zeros(1),
39
40 zeros(1,1,1,1),
41 zeros(1,1,1,1),
42 zeros(1,1,1,1),
43 zeros(1,1,1,1),
44 Array{Array{Int, 1}}(undef, 1, 1),
45 zeros(1,1),
46 1, 1, 1, 1,
47 false, [0.,0.,0.], [1.,1.,1.], [1,1,1], [1.,1.,1.])
48 end
49
50 export readOceanNetCDF
51 """
52 Read ocean NetCDF files generated by MOM6 from disk and return as `Ocean` data
53 structure.
54
55 # Arguments
56 * `velocity_file::String`: path to NetCDF file containing ocean velocities,
57 etc., (e.g. `prog__####_###.nc`).
58 * `grid_file::String`: path to NetCDF file containing ocean super-grid
59 information (typically `INPUT/ocean_hgrid.nc`).
60 * `regular_grid::Bool=false`: `true` if the grid is regular (all cells
61 equal and grid is Cartesian) or `false` (default).
62 """
63 function readOceanNetCDF(velocity_file::String, grid_file::String;
64 regular_grid::Bool=false)
65
66 if !hasNetCDF
67 @warn "Package NetCDF not found. " *
68 "Ocean/atmosphere grid read not supported. " *
69 "Please install NetCDF and its " *
70 "requirements with `Pkg.add(\"NetCDF\")`."
71 else
72
73 time, u, v, h, e, zl, zi = readOceanStateNetCDF(velocity_file)
74 xh, yh, xq, yq = readOceanGridNetCDF(grid_file)
75
76 if size(u[:,:,1,1]) != size(xq) || size(v[:,:,1,1]) != size(xq) ||
77 size(xq) != size(yq)
78 error("size mismatch between velocities and grid
79 (u: $(size(u[:,:,1,1])), v: $(size(v[:,:,1,1])),
80 xq: $(size(xq)), yq: $(size(yq)))")
81 end
82
83 ocean = Ocean([grid_file, velocity_file],
84
85 time,
86
87 xq,
88 yq,
89 xh,
90 yh,
91
92 zl,
93 zi,
94
95 u,
96 v,
97 h,
98 e,
99 Array{Array{Int, 1}}(undef, size(xh, 1), size(xh, 2)),
100 zeros(size(xh)),
101 1, 1, 1, 1,
102
103 false, [0.,0.,0.], [1.,1.,1.], [1,1,1], [1.,1.,1.]
104 )
105 return ocean
106 end
107 end
108
109 export readOceanStateNetCDF
110 """
111 Read NetCDF file with ocean state generated by MOM6 (e.g. `prog__####_###.nc`
112 or `########.ocean_month.nc`) from disk and return time stamps, velocity fields,
113 layer thicknesses, interface heights, and vertical coordinates.
114
115 # Returns
116 * `time::Vector{Float64}`: Time [s]
117 * `u::Array{Float64, 2}`: Cell corner zonal velocity [m/s],
118 dimensions correspond to placement in `[xq, yq, zl, time]`
119 * `v::Array{Float64, 2}`: Cell corner meridional velocity [m/s],
120 dimensions correspond to placement in `[xq, yq, zl, time]`
121 * `h::Array{Float64, 2}`: layer thickness [m], dimensions correspond to
122 placement in `[xh, yh, zl, time]`
123 * `e::Array{Float64, 2}`: interface height relative to mean sea level [m],
124 dimensions correspond to placement in `[xh, yh, zi, time]`
125 * `zl::Vector{Float64}`: layer target potential density [kg m^-3]
126 * `zi::Vector{Float64}`: interface target potential density [kg m^-3]
127 """
128 function readOceanStateNetCDF(filename::String)
129
130 if !hasNetCDF
131 @warn "Package NetCDF not found. " *
132 "Ocean/atmosphere grid read not supported. " *
133 "Please install NetCDF and its " *
134 "requirements with `Pkg.add(\"NetCDF\")`."
135 else
136
137 if !isfile(filename)
138 error("$(filename) could not be opened")
139 end
140
141 u_staggered = convert(Array{Float64, 4}, NetCDF.ncread(filename, "u"))
142 v_staggered = convert(Array{Float64, 4}, NetCDF.ncread(filename, "v"))
143 u, v = interpolateOceanVelocitiesToCorners(u_staggered, v_staggered)
144
145 time = convert(Vector{Float64},
146 NetCDF.ncread(filename, "time") .* 24. * 60. * 60.)
147 h = convert(Array{Float64, 4}, NetCDF.ncread(filename, "h"))
148 e = convert(Array{Float64, 4}, NetCDF.ncread(filename, "e"))
149
150 zl = convert(Vector{Float64}, NetCDF.ncread(filename, "zl"))
151 zi = convert(Vector{Float64}, NetCDF.ncread(filename, "zi"))
152
153 return time, u, v, h, e, zl, zi
154 end
155 end
156
157 export readOceanGridNetCDF
158 """
159 Read NetCDF file with ocean *supergrid* information generated by MOM6 (e.g.
160 `ocean_hrid.nc`) from disk and return as `Ocean` data structure. This file is
161 located in the simulation `INPUT/` subdirectory.
162
163 # Returns
164 * `xh::Array{Float64, 2}`: Longitude for cell centers [deg]
165 * `yh::Array{Float64, 2}`: Latitude for cell centers [deg]
166 * `xq::Array{Float64, 2}`: Longitude for cell corners [deg]
167 * `yq::Array{Float64, 2}`: Latitude for cell corners [deg]
168 """
169 function readOceanGridNetCDF(filename::String)
170
171 if !hasNetCDF
172 @warn "Package NetCDF not found. " *
173 "Ocean/atmosphere grid read not supported. " *
174 "Please install NetCDF and its " *
175 "requirements with `Pkg.add(\"NetCDF\")`."
176 else
177
178 if !isfile(filename)
179 error("$(filename) could not be opened")
180 end
181 x = convert(Array{Float64, 2}, NetCDF.ncread(filename, "x"))
182 y = convert(Array{Float64, 2}, NetCDF.ncread(filename, "y"))
183
184 xh = x[2:2:end, 2:2:end]
185 yh = y[2:2:end, 2:2:end]
186
187 xq = x[1:2:end, 1:2:end]
188 yq = y[1:2:end, 1:2:end]
189
190 return xh, yh, xq, yq
191 end
192 end
193
194 export interpolateOceanVelocitiesToCorners
195 """
196 Convert gridded data from Arakawa-C type (decomposed velocities at faces) to
197 Arakawa-B type (velocities at corners) through interpolation.
198 """
199 function interpolateOceanVelocitiesToCorners(u_in::Array{Float64, 4},
200 v_in::Array{Float64, 4})
201
202 if size(u_in) != size(v_in)
203 error("size of u_in ($(size(u_in))) must match v_in ($(size(v_in)))")
204 end
205
206 nx, ny, nz, nt = size(u_in)
207 u = zeros(nx+1, ny+1, nz, nt)
208 v = zeros(nx+1, ny+1, nz, nt)
209 for i=1:nx
210 for j=1:ny
211 if j < ny - 1
212 u[i, j, :, :] = (u_in[i, j, :, :] + u_in[i, j+1, :, :])/2.
213 else
214 u[i, j, :, :] = u_in[i, j, :, :]
215 end
216 if i < nx - 1
217 v[i, j, :, :] = (v_in[i, j, :, :] + v_in[i+1, j, :, :])/2.
218 else
219 v[i, j, :, :] = v_in[i, j, :, :]
220 end
221 end
222 end
223 return u, v
224 end
225
226 export interpolateOceanState
227 """
228 Ocean data is containted in `Ocean` type at discrete times (`Ocean.time`). This
229 function performs linear interpolation between time steps to get the approximate
230 ocean state at any point in time. If the `Ocean` data set only contains a
231 single time step, values from that time are returned.
232 """
233 function interpolateOceanState(ocean::Ocean, t::Float64)
234 if length(ocean.time) == 1
235 return ocean.u, ocean.v, ocean.h, ocean.e
236 elseif t < ocean.time[1] || t > ocean.time[end]
237 error("selected time (t = $(t)) is outside the range of time steps " *
238 "in the ocean data")
239 end
240
241 i = 1
242 rel_time = 0.
243 while i < length(ocean.time)
244 if ocean.time[i+1] < t
245 i += 1
246 continue
247 end
248
249 dt = ocean.time[i+1] - ocean.time[i]
250 rel_time = (t - ocean.time[i])/dt
251 if rel_time < 0. || rel_time > 1.
252 error("time bounds error")
253 end
254 break
255 end
256
257 return ocean.u[:,:,:,i]*(1. - rel_time) + ocean.u[:,:,:,i+1]*rel_time,
258 ocean.v[:,:,:,i]*(1. - rel_time) + ocean.v[:,:,:,i+1]*rel_time,
259 ocean.h[:,:,:,i]*(1. - rel_time) + ocean.h[:,:,:,i+1]*rel_time,
260 ocean.e[:,:,:,i]*(1. - rel_time) + ocean.e[:,:,:,i+1]*rel_time
261 end
262
263 export createRegularOceanGrid
264 """
265
266 createRegularOceanGrid(n, L[, origo, time, name,
267 bc_west, bc_south, bc_east, bc_north])
268
269 Initialize and return a regular, Cartesian `Ocean` grid with `n[1]` by `n[2]`
270 cells in the horizontal dimension, and `n[3]` vertical cells. The cell corner
271 and center coordinates will be set according to the grid spatial dimensions
272 `L[1]`, `L[2]`, and `L[3]`. The grid `u`, `v`, `h`, and `e` fields will contain
273 one 4-th dimension matrix per `time` step. Sea surface will be at `z=0.` with
274 the ocean spanning `z<0.`. Vertical indexing starts with `k=0` at the sea
275 surface, and increases downwards.
276
277 # Arguments
278 * `n::Vector{Int}`: number of cells along each dimension [-].
279 * `L::Vector{Float64}`: domain length along each dimension [m].
280 * `origo::Vector{Float64}`: domain offset in each dimension [m] (default =
281 `[0.0, 0.0]`).
282 * `time::Vector{Float64}`: vector of time stamps for the grid [s].
283 * `name::String`: grid name (default = `"unnamed"`).
284 * `bc_west::Integer`: grid boundary condition for the grains.
285 * `bc_south::Integer`: grid boundary condition for the grains.
286 * `bc_east::Integer`: grid boundary condition for the grains.
287 * `bc_north::Integer`: grid boundary condition for the grains.
288 """
289 function createRegularOceanGrid(n::Vector{Int},
290 L::Vector{Float64};
291 origo::Vector{Float64} = zeros(2),
292 time::Vector{Float64} = zeros(1),
293 name::String = "unnamed",
294 bc_west::Integer = 1,
295 bc_south::Integer = 1,
296 bc_east::Integer = 1,
297 bc_north::Integer = 1)
298
299 xq = repeat(range(origo[1], stop=origo[1] + L[1],
300 length=n[1] + 1),
301 outer=[1, n[2] + 1])
302 yq = repeat(range(origo[2], stop=origo[2] + L[2],
303 length=n[2] + 1)',
304 outer=[n[1] + 1, 1])
305
306 dx = L./n
307 xh = repeat(range(origo[1] + .5*dx[1],
308 stop=origo[1] + L[1] - .5*dx[1],
309 length=n[1]),
310 outer=[1, n[2]])
311 yh = repeat(range(origo[2] + .5*dx[2],
312 stop=origo[2] + L[2] - .5*dx[2],
313 length=n[2])',
314 outer=[n[1], 1])
315
316 zl = -range(.5*dx[3], stop=L[3] - .5*dx[3], length=n[3])
317 zi = -range(0., stop=L[3], length=n[3] + 1)
318
319 u = zeros(n[1] + 1, n[2] + 1, n[3], length(time))
320 v = zeros(n[1] + 1, n[2] + 1, n[3], length(time))
321 h = zeros(n[1] + 1, n[2] + 1, n[3], length(time))
322 e = zeros(n[1] + 1, n[2] + 1, n[3], length(time))
323
324 return Ocean(name,
325 time,
326 xq, yq,
327 xh, yh,
328 zl, zi,
329 u, v, h, e,
330 Array{Array{Int, 1}}(undef, size(xh, 1), size(xh, 2)),
331 zeros(size(xh)),
332 bc_west, bc_south, bc_east, bc_north,
333 true, origo, L, n, dx)
334 end
335
336 export addOceanDrag!
337 """
338 Add drag from linear and angular velocity difference between ocean and all ice
339 floes.
340 """
341 function addOceanDrag!(simulation::Simulation)
342 if typeof(simulation.ocean.input_file) == Bool
343 error("no ocean data read")
344 end
345
346 u, v, h, e = interpolateOceanState(simulation.ocean, simulation.time)
347 uv_interp = Vector{Float64}(undef, 2)
348 sw = Vector{Float64}(undef, 2)
349 se = Vector{Float64}(undef, 2)
350 ne = Vector{Float64}(undef, 2)
351 nw = Vector{Float64}(undef, 2)
352
353 for grain in simulation.grains
354
355 if !grain.enabled
356 continue
357 end
358
359 i, j = grain.ocean_grid_pos
360 k = 1
361
362 x_tilde, y_tilde = getNonDimensionalCellCoordinates(simulation.ocean,
363 i, j,
364 grain.lin_pos)
365 x_tilde = clamp(x_tilde, 0., 1.)
366 y_tilde = clamp(y_tilde, 0., 1.)
367
368 bilinearInterpolation!(uv_interp, u, v, x_tilde, y_tilde, i, j, k, 1)
369 applyOceanDragToGrain!(grain, uv_interp[1], uv_interp[2])
370 applyOceanVorticityToGrain!(grain,
371 curl(simulation.ocean, x_tilde, y_tilde,
372 i, j, k, 1, sw, se, ne, nw))
373 end
374 nothing
375 end
376
377 export applyOceanDragToGrain!
378 """
379 Add Stokes-type drag from velocity difference between ocean and a single ice
380 floe.
381 """
382 function applyOceanDragToGrain!(grain::GrainCylindrical,
383 u::Float64, v::Float64)
384 freeboard = .1*grain.thickness # height above water
385 ρ_o = 1000. # ocean density
386 draft = grain.thickness - freeboard # height of submerged thickness
387
388 drag_force = ρ_o * π *
389 (2.0*grain.ocean_drag_coeff_vert*grain.areal_radius*draft +
390 grain.ocean_drag_coeff_horiz*grain.areal_radius^2.0) *
391 ([u, v] - grain.lin_vel[1:2])*norm([u, v] - grain.lin_vel[1:2])
392
393 grain.force += vecTo3d(drag_force)
394 grain.ocean_stress = vecTo3d(drag_force/grain.horizontal_surface_area)
395 nothing
396 end
397
398 export applyOceanVorticityToGrain!
399 """
400 Add Stokes-type torque from angular velocity difference between ocean and a
401 single grain. See Eq. 9.28 in "Introduction to Fluid Mechanics" by Nakayama
402 and Boucher, 1999.
403 """
404 function applyOceanVorticityToGrain!(grain::GrainCylindrical,
405 ocean_curl::Float64)
406 freeboard = .1*grain.thickness # height above water
407 ρ_o = 1000. # ocean density
408 draft = grain.thickness - freeboard # height of submerged thickness
409
410 grain.torque[3] +=
411 π * grain.areal_radius^4. * ρ_o *
412 (grain.areal_radius/5. * grain.ocean_drag_coeff_horiz +
413 draft * grain.ocean_drag_coeff_vert) *
414 norm(.5 * ocean_curl - grain.ang_vel[3]) *
415 (.5 * ocean_curl - grain.ang_vel[3])
416 nothing
417 end
418
419 export compareOceans
420 """
421 compareOceans(ocean1::Ocean, ocean2::Ocean)
422
423 Compare values of two `Ocean` objects using the `Base.Test` framework.
424 """
425 function compareOceans(ocean1::Ocean, ocean2::Ocean)
426
427 @test ocean1.input_file == ocean2.input_file
428 @test ocean1.time ≈ ocean2.time
429
430 @test ocean1.xq ≈ ocean2.xq
431 @test ocean1.yq ≈ ocean2.yq
432
433 @test ocean1.xh ≈ ocean2.xh
434 @test ocean1.yh ≈ ocean2.yh
435
436 @test ocean1.zl ≈ ocean2.zl
437 @test ocean1.zi ≈ ocean2.zi
438
439 @test ocean1.u ≈ ocean2.u
440 @test ocean1.v ≈ ocean2.v
441 @test ocean1.h ≈ ocean2.h
442 @test ocean1.e ≈ ocean2.e
443
444 if isassigned(ocean1.grain_list, 1)
445 @test ocean1.grain_list == ocean2.grain_list
446 end
447 nothing
448 end