tatmosphere.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
---
tatmosphere.jl (8067B)
---
1 using Test
2 using LinearAlgebra
3
4 export createEmptyAtmosphere
5 "Returns empty ocean type for initialization purposes."
6 function createEmptyAtmosphere()
7 return Atmosphere(false,
8
9 zeros(1),
10
11 zeros(1,1),
12 zeros(1,1),
13 zeros(1,1),
14 zeros(1,1),
15
16 zeros(1),
17
18 zeros(1,1,1,1),
19 zeros(1,1,1,1),
20
21 Array{Vector{Int}}(undef, 1, 1),
22 zeros(1,1),
23
24 1, 1, 1, 1,
25
26 false,
27
28 false, [0.,0.,0.], [1.,1.,1.], [1,1,1], [1.,1.,1.])
29 end
30
31 export interpolateAtmosphereState
32 """
33 Atmosphere data is containted in `Atmosphere` type at discrete times
34 (`Atmosphere.time`). This function performs linear interpolation between time
35 steps to get the approximate atmosphere state at any point in time. If the
36 `Atmosphere` data set only contains a single time step, values from that time
37 are returned.
38 """
39 function interpolateAtmosphereState(atmosphere::Atmosphere, t::Float64)
40 if length(atmosphere.time) == 1
41 return atmosphere.u, atmosphere.v
42 elseif t < atmosphere.time[1] || t > atmosphere.time[end]
43 error("selected time (t = $(t)) is outside the range of " *
44 "time steps in the atmosphere data")
45 end
46 end
47
48 export createRegularAtmosphereGrid
49 """
50 createRegularAtmosphereGrid(n, L[, origo, time, name,
51 bc_west, bc_south, bc_east, bc_north])
52
53 Initialize and return a regular, Cartesian `Atmosphere` grid with `n[1]` by
54 `n[2]` cells in the horizontal dimension, and `n[3]` vertical cells. The cell
55 corner and center coordinates will be set according to the grid spatial
56 dimensions `L[1]`, `L[2]`, and `L[3]`. The grid `u`, `v`, `h`, and `e` fields
57 will contain one 4-th dimension matrix per `time` step. Sea surface will be at
58 `z=0.` with the atmosphere spanning `z<0.`. Vertical indexing starts with `k=0`
59 at the sea surface, and increases downwards.
60
61 # Arguments
62 * `n::Vector{Int}`: number of cells along each dimension [-].
63 * `L::Vector{Float64}`: domain length along each dimension [m].
64 * `origo::Vector{Float64}`: domain offset in each dimension [m] (default =
65 `[0.0, 0.0]`).
66 * `time::Vector{Float64}`: vector of time stamps for the grid [s].
67 * `name::String`: grid name (default = `"unnamed"`).
68 * `bc_west::Integer`: grid boundary condition for the grains.
69 * `bc_south::Integer`: grid boundary condition for the grains.
70 * `bc_east::Integer`: grid boundary condition for the grains.
71 * `bc_north::Integer`: grid boundary condition for the grains.
72 """
73 function createRegularAtmosphereGrid(n::Vector{Int},
74 L::Vector{Float64};
75 origo::Vector{Float64} = zeros(2),
76 time::Array{Float64, 1} = zeros(1),
77 name::String = "unnamed",
78 bc_west::Integer = 1,
79 bc_south::Integer = 1,
80 bc_east::Integer = 1,
81 bc_north::Integer = 1)
82
83 xq = repeat(range(origo[1], stop=origo[1] + L[1], length=n[1] + 1),
84 outer=[1, n[2] + 1])
85 yq = repeat(range(origo[2], stop=origo[2] + L[2], length=n[2] + 1)',
86 outer=[n[1] + 1, 1])
87
88 dx = L./n
89 xh = repeat(range(origo[1] + .5*dx[1],
90 stop=origo[1] + L[1] - .5*dx[1],
91 length=n[1]),
92 outer=[1, n[2]])
93 yh = repeat(range(origo[2] + .5*dx[2],
94 stop=origo[1] + L[2] - .5*dx[2],
95 length=n[2])',
96 outer=[n[1], 1])
97
98 zl = -range(.5*dx[3], stop=L[3] - .5*dx[3], length=n[3])
99
100 u = zeros(n[1] + 1, n[2] + 1, n[3], length(time))
101 v = zeros(n[1] + 1, n[2] + 1, n[3], length(time))
102
103 return Atmosphere(name,
104 time,
105 xq, yq,
106 xh, yh,
107 zl,
108 u, v,
109 Array{Array{Int, 1}}(undef, size(xh, 1), size(xh, 2)),
110 zeros(size(xh)),
111 bc_west, bc_south, bc_east, bc_north,
112 false,
113 true, origo, L, n, dx)
114 end
115
116 export addAtmosphereDrag!
117 """
118 Add drag from linear and angular velocity difference between atmosphere and all
119 grains.
120 """
121 function addAtmosphereDrag!(simulation::Simulation)
122 if typeof(simulation.atmosphere.input_file) == Bool
123 error("no atmosphere data read")
124 end
125
126 u, v = interpolateAtmosphereState(simulation.atmosphere, simulation.time)
127 uv_interp = Vector{Float64}(undef, 2)
128 sw = Vector{Float64}(undef, 2)
129 se = Vector{Float64}(undef, 2)
130 ne = Vector{Float64}(undef, 2)
131 nw = Vector{Float64}(undef, 2)
132
133 for grain in simulation.grains
134
135 if !grain.enabled
136 continue
137 end
138
139 i, j = grain.atmosphere_grid_pos
140 k = 1
141
142 x_tilde, y_tilde = getNonDimensionalCellCoordinates(simulation.
143 atmosphere,
144 i, j,
145 grain.lin_pos)
146 x_tilde = clamp(x_tilde, 0., 1.)
147 y_tilde = clamp(y_tilde, 0., 1.)
148
149 bilinearInterpolation!(uv_interp, u, v, x_tilde, y_tilde, i, j, k, 1)
150 applyAtmosphereDragToGrain!(grain, uv_interp[1], uv_interp[2])
151 applyAtmosphereVorticityToGrain!(grain,
152 curl(simulation.atmosphere,
153 x_tilde, y_tilde,
154 i, j, k, 1, sw, se, ne, nw))
155 end
156 nothing
157 end
158
159 export applyAtmosphereDragToGrain!
160 """
161 Add Stokes-type drag from velocity difference between atmosphere and a single
162 grain.
163 """
164 function applyAtmosphereDragToGrain!(grain::GrainCylindrical,
165 u::Float64, v::Float64)
166 ρ_a = 1.2754 # atmosphere density
167
168 drag_force = ρ_a * π *
169 (2.0*grain.ocean_drag_coeff_vert*grain.areal_radius*.1*grain.thickness +
170 grain.atmosphere_drag_coeff_horiz*grain.areal_radius^2.0) *
171 ([u, v] - grain.lin_vel[1:2])*norm([u, v] - grain.lin_vel[1:2])
172
173 grain.force += vecTo3d(drag_force)
174 grain.atmosphere_stress = vecTo3d(drag_force/grain.horizontal_surface_area)
175 nothing
176 end
177
178 export applyAtmosphereVorticityToGrain!
179 """
180 Add Stokes-type torque from angular velocity difference between atmosphere and a
181 single grain. See Eq. 9.28 in "Introduction to Fluid Mechanics" by Nakayama
182 and Boucher, 1999.
183 """
184 function applyAtmosphereVorticityToGrain!(grain::GrainCylindrical,
185 atmosphere_curl::Float64)
186 ρ_a = 1.2754 # atmosphere density
187
188 grain.torque[3] +=
189 π * grain.areal_radius^4. * ρ_a *
190 (grain.areal_radius / 5. * grain.atmosphere_drag_coeff_horiz +
191 .1 * grain.thickness * grain.atmosphere_drag_coeff_vert) *
192 norm(.5 * atmosphere_curl - grain.ang_vel[3]) *
193 (.5 * atmosphere_curl - grain.ang_vel[3])
194 nothing
195 end
196
197 export compareAtmospheres
198 """
199 compareAtmospheres(atmosphere1::atmosphere, atmosphere2::atmosphere)
200
201 Compare values of two `atmosphere` objects using the `Base.Test` framework.
202 """
203 function compareAtmospheres(atmosphere1::Atmosphere, atmosphere2::Atmosphere)
204
205 @test atmosphere1.input_file == atmosphere2.input_file
206 @test atmosphere1.time ≈ atmosphere2.time
207
208 @test atmosphere1.xq ≈ atmosphere2.xq
209 @test atmosphere1.yq ≈ atmosphere2.yq
210
211 @test atmosphere1.xh ≈ atmosphere2.xh
212 @test atmosphere1.yh ≈ atmosphere2.yh
213
214 @test atmosphere1.zl ≈ atmosphere2.zl
215
216 @test atmosphere1.u ≈ atmosphere2.u
217 @test atmosphere1.v ≈ atmosphere2.v
218
219 if isassigned(atmosphere1.grain_list, 1)
220 @test atmosphere1.grain_list == atmosphere2.grain_list
221 end
222 nothing
223 end