atmosphere.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