tgrid.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
---
tgrid.jl (40906B)
---
1 import Random
2 using LinearAlgebra
3 using Random
4
5 """
6 bilinearInterpolation(field, x_tilde, y_tilde, i, j, k, it)
7
8 Use bilinear interpolation to interpolate a staggered grid to an arbitrary
9 position in a cell. Assumes south-west convention, i.e. (i,j) is located at the
10 south-west (-x, -y)-facing corner.
11
12 # Arguments
13 * `field::Array{Float64, 4}`: a scalar field to interpolate from
14 * `x_tilde::Float64`: x point position [0;1]
15 * `y_tilde::Float64`: y point position [0;1]
16 * `i::Int`: i-index of cell containing point
17 * `j::Int`: j-index of scalar field to interpolate from
18 * `it::Int`: time step from scalar field to interpolate from
19 """
20 @inline function bilinearInterpolation!(interp_val::Vector{Float64},
21 field_x::Array{Float64, 2},
22 field_y::Array{Float64, 2},
23 x_tilde::Float64,
24 y_tilde::Float64,
25 i::Int,
26 j::Int)
27
28 #if x_tilde < 0. || x_tilde > 1. || y_tilde < 0. || y_tilde > 1.
29 #error("relative coordinates outside bounds ($(x_tilde), $(y_tilde))")
30 #end
31
32 @views interp_val .=
33 (field_x[i+1, j+1]*x_tilde + field_x[i, j+1]*(1. - x_tilde))*y_tilde +
34 (field_x[i+1, j]*x_tilde + field_x[i, j]*(1. - x_tilde))*(1. - y_tilde),
35 (field_y[i+1, j+1]*x_tilde + field_y[i, j+1]*(1. - x_tilde))*y_tilde +
36 (field_y[i+1, j]*x_tilde + field_y[i, j]*(1. - x_tilde))*(1. - y_tilde)
37
38 nothing
39 end
40 @inbounds @inline function bilinearInterpolation!(interp_val::Vector{Float64},
41 field_x::Array{Float64, 4},
42 field_y::Array{Float64, 4},
43 x_tilde::Float64,
44 y_tilde::Float64,
45 i::Int,
46 j::Int,
47 k::Int,
48 it::Int)
49
50 #if x_tilde < 0. || x_tilde > 1. || y_tilde < 0. || y_tilde > 1.
51 #error("relative coordinates outside bounds ($(x_tilde), $(y_tilde))")
52 #end
53
54 @views interp_val .=
55 (field_x[i+1, j+1, k, it]*x_tilde +
56 field_x[i, j+1, k, it]*(1. - x_tilde))*y_tilde +
57 (field_x[i+1, j, k, it]*x_tilde +
58 field_x[i, j, k, it]*(1. - x_tilde))*(1. - y_tilde),
59 (field_y[i+1, j+1, k, it]*x_tilde +
60 field_y[i, j+1, k, it]*(1. - x_tilde))*y_tilde +
61 (field_y[i+1, j, k, it]*x_tilde +
62 field_y[i, j, k, it]*(1. - x_tilde))*(1. - y_tilde)
63
64 nothing
65 end
66
67 """
68 curl(grid, x_tilde, y_tilde, i, j, k, it)
69
70 Use bilinear interpolation to interpolate curl value for a staggered velocity
71 grid to an arbitrary position in a cell. Assumes south-west convention, i.e.
72 (i,j) is located at the south-west (-x, -y)-facing corner.
73
74 # Arguments
75 * `grid::Any`: grid for which to determine curl
76 * `x_tilde::Float64`: x point position [0;1]
77 * `y_tilde::Float64`: y point position [0;1]
78 * `i::Int`: i-index of cell containing point
79 * `j::Int`: j-index of scalar field to interpolate from
80 * `it::Int`: time step from scalar field to interpolate from
81 """
82 function curl(grid::Any,
83 x_tilde::Float64,
84 y_tilde::Float64,
85 i::Int,
86 j::Int,
87 k::Int,
88 it::Int,
89 sw::Vector{Float64} = Vector{Float64}(undef, 2),
90 se::Vector{Float64} = Vector{Float64}(undef, 2),
91 ne::Vector{Float64} = Vector{Float64}(undef, 2),
92 nw::Vector{Float64} = Vector{Float64}(undef, 2))
93
94 #sw, se, ne, nw = getCellCornerCoordinates(grid.xq, grid.yq, i, j)
95 sw[1] = grid.xq[ i, j]
96 sw[2] = grid.yq[ i, j]
97 se[1] = grid.xq[i+1, j]
98 se[2] = grid.yq[i+1, j]
99 ne[1] = grid.xq[i+1, j+1]
100 ne[2] = grid.yq[i+1, j+1]
101 nw[1] = grid.xq[ i, j+1]
102 nw[2] = grid.yq[ i, j+1]
103 sw_se = norm(sw - se)
104 se_ne = norm(se - ne)
105 nw_ne = norm(nw - ne)
106 sw_nw = norm(sw - nw)
107
108 @views @inbounds return (
109 ((grid.v[i+1, j , k,it] - grid.v[i , j , k,it])/sw_se*(1. - y_tilde) +
110 ((grid.v[i+1, j+1, k,it] - grid.v[i , j+1, k,it])/nw_ne)*y_tilde) -
111 ((grid.u[i , j+1, k,it] - grid.u[i , j , k,it])/sw_nw*(1. - x_tilde) +
112 ((grid.u[i+1, j+1, k,it] - grid.u[i+1, j , k,it])/se_ne)*x_tilde))
113 end
114
115 export sortGrainsInGrid!
116 """
117 Find grain positions in grid, based on their center positions.
118 """
119 function sortGrainsInGrid!(simulation::Simulation, grid::Any; verbose=true)
120
121 if simulation.time_iteration == 0
122 grid.grain_list =
123 Array{Array{Int, 1}}(undef, size(grid.xh, 1), size(grid.xh, 2))
124
125 for i=1:size(grid.xh, 1)
126 for j=1:size(grid.xh, 2)
127 @inbounds grid.grain_list[i, j] = Int[]
128 end
129 end
130 else
131 for i=1:size(grid.xh, 1)
132 for j=1:size(grid.xh, 2)
133 @inbounds empty!(grid.grain_list[i, j])
134 end
135 end
136 end
137
138 sw = Vector{Float64}(undef, 2)
139 se = Vector{Float64}(undef, 2)
140 ne = Vector{Float64}(undef, 2)
141 nw = Vector{Float64}(undef, 2)
142
143 for idx=1:length(simulation.grains)
144
145 @inbounds if !simulation.grains[idx].enabled
146 continue
147 end
148
149 # After first iteration, check if grain is in same cell before
150 # traversing entire grid
151 if typeof(grid) == Ocean
152 @inbounds i_old, j_old = simulation.grains[idx].ocean_grid_pos
153 elseif typeof(grid) == Atmosphere
154 @inbounds i_old, j_old =
155 simulation.grains[idx].atmosphere_grid_pos
156 else
157 error("grid type not understood.")
158 end
159
160 if simulation.time > 0. &&
161 !grid.regular_grid &&
162 i_old > 0 && j_old > 0 &&
163 isPointInCell(grid, i_old, j_old,
164 simulation.grains[idx].lin_pos[1:2], sw, se, ne, nw)
165 i = i_old
166 j = j_old
167
168 else
169
170 if grid.regular_grid
171 i, j = Int.(floor.((simulation.grains[idx].lin_pos[1:2]
172 - grid.origo)
173 ./ grid.dx[1:2])) + [1,1]
174 else
175
176 # Search for point in 8 neighboring cells
177 nx = size(grid.xh, 1)
178 ny = size(grid.xh, 2)
179 found = false
180 for i_rel=-1:1
181 for j_rel=-1:1
182 if i_rel == 0 && j_rel == 0
183 continue # cell previously searched
184 end
185 i_t = max(min(i_old + i_rel, nx), 1)
186 j_t = max(min(j_old + j_rel, ny), 1)
187
188 @inbounds if isPointInCell(grid, i_t, j_t,
189 simulation.grains[idx].
190 lin_pos[1:2],
191 sw, se, ne, nw)
192 i = i_t
193 j = j_t
194 found = true
195 break
196 end
197 end
198 if found
199 break
200 end
201 end
202
203 if !found
204 i, j = findCellContainingPoint(grid,
205 simulation.grains[idx].
206 lin_pos[1:2],
207 sw, se, ne, nw)
208 end
209 end
210
211 # remove grain if it is outside of the grid
212 if (!grid.regular_grid &&
213 (i < 1 || j < 1 ||
214 i > size(grid.xh, 1) || j > size(grid.xh, 2))) ||
215 (grid.regular_grid &&
216 (i < 1 || j < 1 ||
217 i > grid.n[1] || j > grid.n[2]))
218
219 if verbose
220 @info "Disabling grain $idx at pos (" *
221 "$(simulation.grains[idx].lin_pos))"
222 end
223 disableGrain!(simulation, idx)
224 continue
225 end
226
227 # add cell to grain
228 if typeof(grid) == Ocean
229 @inbounds simulation.grains[idx].ocean_grid_pos[1] = i
230 @inbounds simulation.grains[idx].ocean_grid_pos[2] = j
231 elseif typeof(grid) == Atmosphere
232 @inbounds simulation.grains[idx].atmosphere_grid_pos[1] = i
233 @inbounds simulation.grains[idx].atmosphere_grid_pos[2] = j
234 else
235 error("grid type not understood.")
236 end
237 end
238
239 # add grain to cell
240 @inbounds push!(grid.grain_list[i, j], idx)
241 end
242 nothing
243 end
244
245 export findCellContainingPoint
246 """
247 findCellContainingPoint(grid, point[, method])
248
249 Returns the `i`, `j` index of the grid cell containing the `point`.
250 The function uses either an area-based approach (`method = "Area"`), or a
251 conformal mapping approach (`method = "Conformal"`). The area-based approach is
252 more robust. This function returns the coordinates of the cell. If no match is
253 found the function returns `(0,0)`.
254
255 # Arguments
256 * `grid::Any`: grid object containing ocean or atmosphere data.
257 * `point::Vector{Float64}`: two-dimensional vector of point to check.
258 * `method::String`: approach to use for determining if point is inside cell or
259 not, can be "Conformal" (default) or "Area".
260 """
261 function findCellContainingPoint(grid::Any, point::Vector{Float64},
262 sw = Vector{Float64}(undef, 2),
263 se = Vector{Float64}(undef, 2),
264 ne = Vector{Float64}(undef, 2),
265 nw = Vector{Float64}(undef, 2);
266 method::String="Conformal")
267 for i=1:size(grid.xh, 1)
268 for j=1:size(grid.yh, 2)
269 if isPointInCell(grid, i, j, point,
270 sw, se, ne, nw,
271 method=method)
272 return i, j
273 end
274 end
275 end
276 return 0, 0
277 end
278
279 export getNonDimensionalCellCoordinates
280 """
281 Returns the non-dimensional conformal mapped coordinates for point `point` in
282 cell `i,j`, based off the coordinates in the grid.
283
284 This function is a wrapper for `getCellCornerCoordinates()` and
285 `conformalQuadrilateralCoordinates()`.
286 """
287 function getNonDimensionalCellCoordinates(grid::Any, i::Int, j::Int,
288 point::Vector{Float64})
289 if grid.regular_grid
290 return (point[1:2] - Float64.([i-1,j-1]).*grid.dx[1:2])./grid.dx[1:2]
291 else
292 sw, se, ne, nw = getCellCornerCoordinates(grid.xq, grid.yq, i, j)
293 return conformalQuadrilateralCoordinates(sw, se, ne, nw, point[1:2])
294 end
295 end
296
297 export isPointInCell
298 """
299 Check if a 2d point is contained inside a cell from the supplied grid.
300 The function uses either an area-based approach (`method = "Area"`), or a
301 conformal mapping approach (`method = "Conformal"`). The area-based approach is
302 more robust. This function returns `true` or `false`.
303 """
304 function isPointInCell(grid::Any, i::Int, j::Int,
305 point::Vector{Float64},
306 sw::Vector{Float64} = Vector{Float64}(undef, 2),
307 se::Vector{Float64} = Vector{Float64}(undef, 2),
308 ne::Vector{Float64} = Vector{Float64}(undef, 2),
309 nw::Vector{Float64} = Vector{Float64}(undef, 2);
310 method::String="Conformal")
311
312 if grid.regular_grid
313 if [i,j] == Int.(floor.((point[1:2] - grid.origo) ./ grid.dx[1:2])) + [1,1]
314 return true
315 else
316 return false
317 end
318 end
319
320 @views sw .= grid.xq[ i, j], grid.yq[ i, j]
321 @views se .= grid.xq[ i+1, j], grid.yq[ i+1, j]
322 @views ne .= grid.xq[ i+1, j+1], grid.yq[ i+1, j+1]
323 @views nw .= grid.xq[ i, j+1], grid.yq[ i, j+1]
324
325 if method == "Area"
326 if areaOfQuadrilateral(sw, se, ne, nw) ≈
327 areaOfTriangle(point, sw, se) +
328 areaOfTriangle(point, se, ne) +
329 areaOfTriangle(point, ne, nw) +
330 areaOfTriangle(point, nw, sw)
331 return true
332 else
333 return false
334 end
335
336 elseif method == "Conformal"
337 x_tilde, y_tilde = conformalQuadrilateralCoordinates(sw, se, ne, nw,
338 point)
339 if x_tilde >= 0. && x_tilde <= 1. && y_tilde >= 0. && y_tilde <= 1.
340 return true
341 else
342 return false
343 end
344 else
345 error("method not understood")
346 end
347 end
348
349 export isPointInGrid
350 """
351 Check if a 2d point is contained inside the grid. The function uses either an
352 area-based approach (`method = "Area"`), or a conformal mapping approach
353 (`method = "Conformal"`). The area-based approach is more robust. This
354 function returns `true` or `false`.
355 """
356 function isPointInGrid(grid::Any, point::Vector{Float64},
357 sw::Vector{Float64} = Vector{Float64}(undef, 2),
358 se::Vector{Float64} = Vector{Float64}(undef, 2),
359 ne::Vector{Float64} = Vector{Float64}(undef, 2),
360 nw::Vector{Float64} = Vector{Float64}(undef, 2);
361 method::String="Conformal")
362
363 #sw, se, ne, nw = getCellCornerCoordinates(grid.xq, grid.yq, i, j)
364 nx, ny = size(grid.xq)
365 @views sw .= grid.xq[ 1, 1], grid.yq[ 1, 1]
366 @views se .= grid.xq[ nx, 1], grid.yq[ nx, 1]
367 @views ne .= grid.xq[ nx, ny], grid.yq[ nx, ny]
368 @views nw .= grid.xq[ 1, ny], grid.yq[ 1, ny]
369
370 if method == "Area"
371 if areaOfQuadrilateral(sw, se, ne, nw) ≈
372 areaOfTriangle(point, sw, se) +
373 areaOfTriangle(point, se, ne) +
374 areaOfTriangle(point, ne, nw) +
375 areaOfTriangle(point, nw, sw)
376 return true
377 else
378 return false
379 end
380
381 elseif method == "Conformal"
382 x_tilde, y_tilde = conformalQuadrilateralCoordinates(sw, se, ne, nw,
383 point)
384 if x_tilde >= 0. && x_tilde <= 1. && y_tilde >= 0. && y_tilde <= 1.
385 return true
386 else
387 return false
388 end
389 else
390 error("method not understood")
391 end
392 end
393
394 export getGridCornerCoordinates
395 """
396 getGridCornerCoordinates(xq, yq)
397
398 Returns grid corner coordinates in the following order (south-west corner,
399 south-east corner, north-east corner, north-west corner).
400
401 # Arguments
402 * `xq::Array{Float64, 2}`: nominal longitude of q-points [degrees_E]
403 * `yq::Array{Float64, 2}`: nominal latitude of q-points [degrees_N]
404 """
405 @inline function getGridCornerCoordinates(xq::Array{Float64, 2},
406 yq::Array{Float64, 2})
407 nx, ny = size(xq)
408 @inbounds return Float64[xq[ 1, 1], yq[ 1, 1]],
409 Float64[xq[ nx, 1], yq[ nx, 1]],
410 Float64[xq[ nx, ny], yq[ nx, ny]],
411 Float64[xq[ 1, ny], yq[ 1, ny]]
412 end
413
414
415 export getCellCornerCoordinates
416 """
417 getCellCornerCoordinates(xq, yq, i, j)
418
419 Returns grid-cell corner coordinates in the following order (south-west corner,
420 south-east corner, north-east corner, north-west corner).
421
422 # Arguments
423 * `xq::Array{Float64, 2}`: nominal longitude of q-points [degrees_E]
424 * `yq::Array{Float64, 2}`: nominal latitude of q-points [degrees_N]
425 * `i::Int`: x-index of cell.
426 * `j::Int`: y-index of cell.
427 """
428 @inline function getCellCornerCoordinates(xq::Array{Float64, 2},
429 yq::Array{Float64, 2},
430 i::Int, j::Int)
431 @inbounds return Float64[xq[ i, j], yq[ i, j]],
432 Float64[xq[i+1, j], yq[i+1, j]],
433 Float64[xq[i+1, j+1], yq[i+1, j+1]],
434 Float64[xq[ i, j+1], yq[ i, j+1]]
435 end
436
437 export getCellCenterCoordinates
438 """
439 getCellCenterCoordinates(grid.xh, grid.yh, i, j)
440
441 Returns grid center coordinates (h-point).
442
443 # Arguments
444 * `xh::Array{Float64, 2}`: nominal longitude of h-points [degrees_E]
445 * `yh::Array{Float64, 2}`: nominal latitude of h-points [degrees_N]
446 * `i::Int`: x-index of cell.
447 * `j::Int`: y-index of cell.
448 """
449 function getCellCenterCoordinates(xh::Array{Float64, 2}, yh::Array{Float64, 2},
450 i::Int, j::Int)
451 return [xh[i, j], yh[i, j]]
452 end
453
454 export areaOfTriangle
455 "Returns the area of an triangle with corner coordinates `a`, `b`, and `c`."
456 function areaOfTriangle(a::Vector{Float64},
457 b::Vector{Float64},
458 c::Vector{Float64})
459 return abs(
460 (a[1]*(b[2] - c[2]) +
461 b[1]*(c[2] - a[2]) +
462 c[1]*(a[2] - b[2]))/2.
463 )
464 end
465
466 export areaOfQuadrilateral
467 """
468 Returns the area of a quadrilateral with corner coordinates `a`, `b`, `c`, and
469 `d`. Corners `a` and `c` should be opposite of each other, the same must be
470 true for `b` and `d`. This is true if the four corners are passed as arguments
471 in a "clockwise" or "counter-clockwise" manner.
472 """
473 function areaOfQuadrilateral(a::Vector{Float64},
474 b::Vector{Float64},
475 c::Vector{Float64},
476 d::Vector{Float64})
477 return areaOfTriangle(a, b, c) + areaOfTriangle(c, d, a)
478 end
479
480 export conformalQuadrilateralCoordinates
481 """
482 Returns the non-dimensional coordinates `[x_tilde, y_tilde]` of a point `p`
483 within a quadrilateral with corner coordinates `A`, `B`, `C`, and `D`.
484 Points must be ordered in counter-clockwise order, starting from south-west
485 corner.
486 """
487 function conformalQuadrilateralCoordinates(A::Vector{Float64},
488 B::Vector{Float64},
489 C::Vector{Float64},
490 D::Vector{Float64},
491 p::Vector{Float64})
492
493 if !(A[1] < B[1] && B[2] < C[2] && C[1] > D[1])
494 error("corner coordinates are not passed in the correct order")
495 end
496 alpha = B[1] - A[1]
497 delta = B[2] - A[2]
498 beta = D[1] - A[1]
499 epsilon = D[2] - A[2]
500 gamma = C[1] - A[1] - (alpha + beta)
501 kappa = C[2] - A[2] - (delta + epsilon)
502 a = kappa*beta - gamma*epsilon
503 dx = p[1] - A[1]
504 dy = p[2] - A[2]
505 b = (delta*beta - alpha*epsilon) - (kappa*dx - gamma*dy)
506 c = alpha*dy - delta*dx
507 if abs(a) > 0.
508 d = b^2. / 4. - a*c
509 if d >= 0.
510 yy1 = -(b / 2. + sqrt(d)) / a
511 yy2 = -(b / 2. - sqrt(d)) / a
512 if abs(yy1 - .5) < abs(yy2 - .5)
513 y_tilde = yy1
514 else
515 y_tilde = yy2
516 end
517 else
518 error("could not perform conformal mapping\n" *
519 "A = $(A), B = $(B), C = $(C), D = $(D), point = $(p),\n" *
520 "alpha = $(alpha), beta = $(beta), gamma = $(gamma), " *
521 "delta = $(delta), epsilon = $(epsilon), kappa = $(kappa)")
522 end
523 else
524 if !(b ≈ 0.)
525 y_tilde = -c/b
526 else
527 y_tilde = 0.
528 end
529 end
530 a = alpha + gamma*y_tilde
531 b = delta + kappa*y_tilde
532 if !(a ≈ 0.)
533 x_tilde = (dx - beta*y_tilde)/a
534 elseif !(b ≈ 0.)
535 x_tilde = (dy - epsilon*y_tilde)/b
536 else
537 error("could not determine non-dimensional position in quadrilateral " *
538 "(a = 0. and b = 0.)\n" *
539 "A = $(A), B = $(B), C = $(C), D = $(D), point = $(p),\n" *
540 "alpha = $(alpha), beta = $(beta), gamma = $(gamma), " *
541 "delta = $(delta), epsilon = $(epsilon), kappa = $(kappa)")
542 end
543 return Float64[x_tilde, y_tilde]
544 end
545
546 export findEmptyPositionInGridCell
547 """
548 findEmptyPositionInGridCell(simulation, grid, i, j, r[, n_iter, seed,
549 verbose])
550
551 Attempt locate an empty spot for an grain with radius `r` with center
552 coordinates in a specified grid cell (`i`, `j`) without overlapping any other
553 grains in that cell or the neighboring cells. This function will stop
554 attempting after `n_iter` iterations, each with randomly generated positions.
555
556 This function assumes that existing grains have been binned according to the
557 grid (e.g., using `sortGrainsInGrid()`).
558
559 If the function sucessfully finds a position it will be returned as a
560 two-component Vector{Float64}. If a position is not found, the function will
561 return `false`.
562
563 # Arguments
564 * `simulation::Simulation`: the simulation object to add grains to.
565 * `grid::Any`: the grid to use for position search.
566 * `i::Int`: the grid-cell index along x.
567 * `j::Int`: the grid-cell index along y.
568 * `r::Float64`: the desired grain radius to fit into the cell.
569 * `n_iter::Int = 30`: the number of attempts for finding an empty spot.
570 * `seed::Int = 1`: seed for the pseudo-random number generator.
571 * `verbose::Bool = false`: print diagnostic information.
572 """
573 function findEmptyPositionInGridCell(simulation::Simulation,
574 grid::Any,
575 i::Int,
576 j::Int,
577 r::Float64;
578 n_iter::Int = 30,
579 seed::Int = 1,
580 verbose::Bool = false)
581 overlap_found = false
582 spot_found = false
583 i_iter = 0
584 pos = [NaN, NaN]
585
586 nx, ny = size(grid.xh)
587
588 for i_iter=1:n_iter
589
590 overlap_found = false
591 Random.seed!(i*j*seed*i_iter)
592 # generate random candidate position
593 x_tilde = rand()
594 y_tilde = rand()
595 bilinearInterpolation!(pos, grid.xq, grid.yq, x_tilde, y_tilde, i, j)
596 if verbose
597 @info "trying position $pos in cell $i,$j"
598 end
599
600 # do not penetrate outside of grid boundaries
601 if i == 1 && pos[1] - r < grid.xq[1,1] ||
602 j == 1 && pos[2] - r < grid.yq[1,1] ||
603 i == nx && pos[1] + r > grid.xq[end,end] ||
604 j == ny && pos[2] + r > grid.yq[end,end]
605 overlap_found = true
606 end
607
608 # search for contacts in current and eight neighboring cells
609 if !overlap_found
610 for i_neighbor_corr=[0 -1 1]
611 for j_neighbor_corr=[0 -1 1]
612
613 # target cell index
614 it = i + i_neighbor_corr
615 jt = j + j_neighbor_corr
616
617 # do not search outside grid boundaries
618 if it < 1 || it > nx || jt < 1 || jt > ny
619 continue
620 end
621
622 # traverse list of grains in the target cell and check
623 # for overlaps
624 for grain_idx in grid.grain_list[it, jt]
625 overlap = norm(simulation.grains[grain_idx].
626 lin_pos[1:2] - pos) -
627 (simulation.grains[grain_idx].contact_radius + r)
628
629 if overlap < 0.
630 if verbose
631 @info "overlap with $grain_idx in cell $i,$j"
632 end
633 overlap_found = true
634 break
635 end
636 end
637 end
638 if overlap_found == true
639 break
640 end
641 end
642 end
643 if overlap_found == false
644 break
645 end
646 end
647 if overlap_found == false
648 spot_found = true
649 end
650
651 if spot_found
652 if verbose
653 @info "Found position $pos in cell $i,$j"
654 end
655 return pos
656 else
657 if verbose
658 @warn "could not insert an grain into " *
659 "$(typeof(grid)) grid cell ($i, $j)"
660 end
661 return false
662 end
663 end
664
665 """
666 Copy grain related information from ocean to atmosphere grid. This is useful
667 when the two grids are of identical geometry, meaning only only one sorting
668 phase is necessary.
669 """
670 function copyGridSortingInfo!(ocean::Ocean, atmosphere::Atmosphere,
671 grains::Array{GrainCylindrical, 1})
672
673 for grain in grains
674 grain.atmosphere_grid_pos = deepcopy(grain.ocean_grid_pos)
675 end
676 atmosphere.grain_list = deepcopy(ocean.grain_list)
677 nothing
678 end
679
680 export setGridBoundaryConditions!
681 """
682 setGridBoundaryConditions!(grid, grid_face, mode)
683
684 Set boundary conditions for the granular phase at the edges of `Ocean` or
685 `Atmosphere` grids. The target boundary can be selected through the `grid_face`
686 argument, or the same boundary condition can be applied to all grid boundaries
687 at once.
688
689 When the center coordinate of grains crosses an inactive boundary (`mode =
690 "inactive"`), the grain is disabled (`GrainCylindrical.enabled = false`). This
691 keeps the grain in memory, but stops it from moving or interacting with other
692 grains. *By default, all boundaries are inactive*.
693
694 If the center coordinate of a grain crosses a periodic boundary (`mode =
695 periodic`), the grain is repositioned to the opposite side of the model domain.
696 Grains can interact mechanically across the periodic boundary.
697
698 # Arguments
699 * `grid::Any`: `Ocean` or `Atmosphere` grid to apply the boundary condition to.
700 * `grid_face::String`: Grid face to apply the boundary condition to. Valid
701 values are any combination and sequence of `"west"` (-x), `"south"` (-y),
702 `"east"` (+x), `"north"` (+y), or simply any combination of `"-x"`, `"+x"`,
703 `"-y"`, and `"+y"`. The specifiers may be delimited in any way.
704 Also, and by default, all boundaries can be selected with `"all"` (-x, -y,
705 +x, +y), which overrides any other face selection.
706 * `mode::String`: Boundary behavior, accepted values are `"inactive"`,
707 `"periodic"`, and `"impermeable"`. You cannot specify more than one mode at
708 a time, so if several modes are desired as boundary conditions for the grid,
709 several calls to this function should be made.
710 * `verbose::Bool`: Confirm boundary conditions by reporting values to console.
711
712 # Examples
713 Set all boundaries for the ocean grid to be periodic:
714
715 setGridBoundaryConditions!(ocean, "periodic", "all")
716
717 Set the south-north boundaries to be inactive, but the west-east boundaries to
718 be periodic:
719
720 setGridBoundaryConditions!(ocean, "inactive", "south north")
721 setGridBoundaryConditions!(ocean, "periodic", "west east")
722
723 or specify the conditions from the coordinate system axes:
724
725 setGridBoundaryConditions!(ocean, "inactive", "-y +y")
726 setGridBoundaryConditions!(ocean, "periodic", "-x +x")
727
728 """
729 function setGridBoundaryConditions!(grid::Any,
730 mode::String,
731 grid_face::String = "all";
732 verbose::Bool=true)
733
734 something_changed = false
735
736 if length(mode) <= 1
737 error("The mode string is required ('$mode')")
738 end
739
740 if !(mode in grid_bc_strings)
741 error("Mode '$mode' not recognized as a valid boundary condition type")
742 end
743
744 if occursin("west", grid_face) || occursin("-x", grid_face)
745 grid.bc_west = grid_bc_flags[mode]
746 something_changed = true
747 end
748
749 if occursin("south", grid_face) || occursin("-y", grid_face)
750 grid.bc_south = grid_bc_flags[mode]
751 something_changed = true
752 end
753
754 if occursin("east", grid_face) || occursin("+x", grid_face)
755 grid.bc_east = grid_bc_flags[mode]
756 something_changed = true
757 end
758
759 if occursin("north", grid_face) || occursin("+y", grid_face)
760 grid.bc_north = grid_bc_flags[mode]
761 something_changed = true
762 end
763
764 if grid_face == "all"
765 grid.bc_west = grid_bc_flags[mode]
766 grid.bc_south = grid_bc_flags[mode]
767 grid.bc_east = grid_bc_flags[mode]
768 grid.bc_north = grid_bc_flags[mode]
769 something_changed = true
770 end
771
772 if !something_changed
773 error("grid_face string '$grid_face' not understood, " *
774 "must be east, west, north, south, -x, +x, -y, and/or +y.")
775 end
776
777 if verbose
778 reportGridBoundaryConditions(grid)
779 end
780 nothing
781 end
782
783 export reportGridBoundaryConditions
784 """
785 reportGridBoundaryConditions(grid)
786
787 Report the boundary conditions for the grid to the console.
788 """
789 function reportGridBoundaryConditions(grid::Any)
790 println("West (-x): " * grid_bc_strings[grid.bc_west] *
791 "\t($(grid.bc_west))")
792 println("East (+x): " * grid_bc_strings[grid.bc_east] *
793 "\t($(grid.bc_east))")
794 println("South (-y): " * grid_bc_strings[grid.bc_south] *
795 "\t($(grid.bc_south))")
796 println("North (+y): " * grid_bc_strings[grid.bc_north] *
797 "\t($(grid.bc_north))")
798 nothing
799 end
800
801 """
802 moveGrainsAcrossPeriodicBoundaries!(simulation::Simulation)
803
804 If the ocean or atmosphere grids are periodic, move grains that are placed
805 outside the domain correspondingly across the domain. This function is to be
806 called after temporal integration of the grain positions.
807 """
808 function moveGrainsAcrossPeriodicBoundaries!(sim::Simulation)
809
810 # return if grids are not enabled
811 if typeof(sim.ocean.input_file) == Bool &&
812 typeof(sim.atmosphere.input_file) == Bool
813 return nothing
814 end
815
816 # return immediately if no boundaries are periodic
817 if sim.ocean.bc_west != 2 &&
818 sim.ocean.bc_south != 2 &&
819 sim.ocean.bc_east != 2 &&
820 sim.ocean.bc_north != 2
821 return nothing
822 end
823
824 # throw error if ocean and atmosphere grid BCs are different and both are
825 # enabled
826 if (typeof(sim.ocean.input_file) != Bool &&
827 typeof(sim.atmosphere.input_file) != Bool) &&
828 (sim.ocean.bc_west != sim.atmosphere.bc_west &&
829 sim.ocean.bc_south != sim.atmosphere.bc_south &&
830 sim.ocean.bc_east != sim.atmosphere.bc_east &&
831 sim.ocean.bc_north != sim.atmosphere.bc_north)
832 error("Ocean and Atmosphere grid boundary conditions differ")
833 end
834
835 for grain in sim.grains
836
837 # -x -> +x
838 if sim.ocean.bc_west == 2 && grain.lin_pos[1] < sim.ocean.xq[1]
839 grain.lin_pos[1] += sim.ocean.xq[end] - sim.ocean.xq[1]
840 end
841
842 # -y -> +y
843 if sim.ocean.bc_south == 2 && grain.lin_pos[2] < sim.ocean.yq[1]
844 grain.lin_pos[2] += sim.ocean.yq[end] - sim.ocean.yq[1]
845 end
846
847 # +x -> -x
848 if sim.ocean.bc_east == 2 && grain.lin_pos[1] > sim.ocean.xq[end]
849 grain.lin_pos[1] -= sim.ocean.xq[end] - sim.ocean.xq[1]
850 end
851
852 # +y -> -y
853 if sim.ocean.bc_north == 2 && grain.lin_pos[2] > sim.ocean.yq[end]
854 grain.lin_pos[2] -= sim.ocean.yq[end] - sim.ocean.yq[1]
855 end
856 end
857 nothing
858 end
859
860 export reflectGrainsFromImpermeableBoundaries!
861 """
862 reflectGrainsFromImpermeableBoundaries!(simulation::Simulation)
863
864 If the ocean or atmosphere grids are impermeable, reflect grain trajectories by
865 reversing the velocity vectors normal to the boundary. This function is to be
866 called after temporal integration of the grain positions.
867 """
868 function reflectGrainsFromImpermeableBoundaries!(sim::Simulation)
869
870 # return if grids are not enabled
871 if typeof(sim.ocean.input_file) == Bool &&
872 typeof(sim.atmosphere.input_file) == Bool
873 return nothing
874 end
875
876 # return immediately if no boundaries are periodic
877 if sim.ocean.bc_west != 3 &&
878 sim.ocean.bc_south != 3 &&
879 sim.ocean.bc_east != 3 &&
880 sim.ocean.bc_north != 3
881 return nothing
882 end
883
884 # throw error if ocean and atmosphere grid BCs are different and both are
885 # enabled
886 if (typeof(sim.ocean.input_file) != Bool &&
887 typeof(sim.atmosphere.input_file) != Bool) &&
888 (sim.ocean.bc_west != sim.atmosphere.bc_west &&
889 sim.ocean.bc_south != sim.atmosphere.bc_south &&
890 sim.ocean.bc_east != sim.atmosphere.bc_east &&
891 sim.ocean.bc_north != sim.atmosphere.bc_north)
892 error("Ocean and Atmosphere grid boundary conditions differ")
893 end
894
895 for grain in sim.grains
896
897 # -x
898 if sim.ocean.bc_west == 3 &&
899 grain.lin_pos[1] - grain.contact_radius < sim.ocean.xq[1]
900
901 grain.lin_vel[1] *= -1.
902 end
903
904 # -y
905 if sim.ocean.bc_south == 3 &&
906 grain.lin_pos[2] - grain.contact_radius < sim.ocean.yq[1]
907
908 grain.lin_vel[2] *= -1.
909 end
910
911 # +x
912 if sim.ocean.bc_east == 3 &&
913 grain.lin_pos[1] + grain.contact_radius > sim.ocean.xq[end]
914
915 grain.lin_vel[1] *= -1.
916 end
917
918 # +y
919 if sim.ocean.bc_east == 3 &&
920 grain.lin_pos[2] + grain.contact_radius > sim.ocean.yq[end]
921
922 grain.lin_vel[2] *= -1.
923 end
924 end
925 nothing
926 end
927
928 export fitGridToGrains!
929 """
930 fitGridToGrains!(simulation, grid[, padding])
931
932 Fit the ocean or atmosphere grid for a simulation to the current grains and
933 their positions.
934
935 # Arguments
936 * `simulation::Simulation`: simulation object to manipulate.
937 * `grid::Any`: Ocean or Atmosphere grid to manipulate.
938 * `padding::Real`: optional padding around edges [m].
939 * `verbose::Bool`: show grid information when function completes.
940 """
941 function fitGridToGrains!(simulation::Simulation, grid::Any;
942 padding::Real=0., verbose::Bool=true)
943
944 if typeof(grid) != Ocean && typeof(grid) != Atmosphere
945 error("grid must be of Ocean or Atmosphere type")
946 end
947
948 min_x = Inf
949 min_y = Inf
950 max_x = -Inf
951 max_y = -Inf
952 max_radius = 0.
953
954 if length(simulation.grains) < 1
955 error("Grains need to be initialized before calling fitGridToGrains")
956 end
957
958 r = 0.
959 for grain in simulation.grains
960 r = grain.contact_radius
961
962 if grid.bc_west == grid_bc_flags["periodic"]
963 if grain.lin_pos[1] < min_x
964 min_x = grain.lin_pos[1] - r
965 end
966 else
967 if grain.lin_pos[1] - r < min_x
968 min_x = grain.lin_pos[1] - r
969 end
970 end
971
972 if grid.bc_east == grid_bc_flags["periodic"]
973 if grain.lin_pos[1] > max_x
974 max_x = grain.lin_pos[1] + grain.contact_radius
975 end
976 else
977 if grain.lin_pos[1] + r > max_x
978 max_x = grain.lin_pos[1] + grain.contact_radius
979 end
980 end
981
982 if grid.bc_south == grid_bc_flags["periodic"]
983 if grain.lin_pos[2] < min_y
984 min_y = grain.lin_pos[2] - grain.contact_radius
985 end
986 else
987 if grain.lin_pos[2] - r < min_y
988 min_y = grain.lin_pos[2] - grain.contact_radius
989 end
990 end
991
992 if grid.bc_north == grid_bc_flags["periodic"]
993 if grain.lin_pos[2] > max_y
994 max_y = grain.lin_pos[2] + grain.contact_radius
995 end
996 else
997 if grain.lin_pos[2] + r > max_y
998 max_y = grain.lin_pos[2] + grain.contact_radius
999 end
1000 end
1001
1002 if r > max_radius
1003 max_radius = r
1004 end
1005 end
1006 min_x -= padding
1007 min_y -= padding
1008 max_x += padding
1009 max_y += padding
1010
1011 L::Vector{Float64} = [max_x - min_x, max_y - min_y]
1012 dx::Float64 = 2. * max_radius
1013 n = convert(Vector{Int}, floor.(L./dx))
1014 if 0 in n || 1 in n
1015 println("L = $L")
1016 println("dx = $dx")
1017 println("n = $n")
1018 error("Grid is too small compared to grain size (n = $n). " *
1019 "Use all-to-all contact search instead.")
1020 end
1021
1022
1023 if typeof(grid) == Ocean
1024 simulation.ocean = createRegularOceanGrid(vcat(n, 1), vcat(L, 1.),
1025 origo=[min_x, min_y],
1026 time=[0.], name="fitted",
1027 bc_west = grid.bc_west,
1028 bc_south = grid.bc_south,
1029 bc_east = grid.bc_east,
1030 bc_north = grid.bc_north)
1031 elseif typeof(grid) == Atmosphere
1032 simulation.atmosphere = createRegularAtmosphereGrid(vcat(n, 1),
1033 vcat(L, 1.),
1034 origo=[min_x,
1035 min_y],
1036 time=[0.],
1037 name="fitted",
1038 bc_west = grid.bc_west,
1039 bc_south = grid.bc_south,
1040 bc_east = grid.bc_east,
1041 bc_north = grid.bc_north)
1042 end
1043
1044 if verbose
1045 @info "Created regular $(typeof(grid)) grid from " *
1046 "[$min_x, $min_y] to [$max_x, $max_y] " *
1047 "with a cell size of $dx ($n)."
1048 end
1049
1050 nothing
1051 end
1052
1053 function findPorosity!(sim::Simulation, grid::Any; verbose::Bool=true)
1054
1055 if !isassigned(grid.grain_list)
1056 @info "Sorting grains in grid"
1057 sortGrainsInGrid!(sim, grid, verbose=verbose)
1058 end
1059
1060 sw = Vector{Float64}(undef, 2)
1061 se = Vector{Float64}(undef, 2)
1062 ne = Vector{Float64}(undef, 2)
1063 nw = Vector{Float64}(undef, 2)
1064 cell_area = 0.0
1065 p = zeros(2)
1066 r = 0.0
1067 A = 0.0
1068
1069 for ix in 1:size(grid.xh, 1)
1070 for iy in 1:size(grid.xh, 2)
1071
1072 @views sw .= grid.xq[ ix, iy], grid.yq[ ix, iy]
1073 @views se .= grid.xq[ ix+1, iy], grid.yq[ ix+1, iy]
1074 @views ne .= grid.xq[ ix+1, iy+1], grid.yq[ ix+1, iy+1]
1075 @views nw .= grid.xq[ ix, iy+1], grid.yq[ ix, iy+1]
1076 cell_area = areaOfQuadrilateral(sw, se, ne, nw)
1077
1078 # Subtract grain area from cell area
1079 particle_area = 0.0
1080 for ix_ = -1:1
1081 for iy_ = -1:1
1082
1083 # Make sure cell check is within grid
1084 if ix + ix_ < 1 || ix + ix_ > size(grid.xh, 1) ||
1085 iy + iy_ < 1 || iy + iy_ > size(grid.xh, 2)
1086 continue
1087 end
1088
1089 # Traverse grain list
1090 for i in grid.grain_list[ix + ix_, iy + iy_]
1091
1092 # Grain geometry
1093 p = sim.grains[i].lin_pos
1094 r = sim.grains[i].areal_radius
1095 A = grainHorizontalSurfaceArea(sim.grains[i])
1096
1097 #= if sw[1] <= p[1] - r && =#
1098 #= sw[2] <= p[2] - r && =#
1099 #= ne[1] >= p[1] + r && =#
1100 #= ne[2] >= p[2] + r =#
1101 if sw[1] <= p[1] &&
1102 sw[2] <= p[2] &&
1103 ne[1] >= p[1] &&
1104 ne[2] >= p[2]
1105 # If particle is entirely contained within cell,
1106 # assuming a regular and orthogonal grid
1107 # TODO: Adjust coordinates with conformal mapping
1108 # for irregular grids.
1109 particle_area += A
1110
1111 #= elseif sw[1] >= p[1] + r || =#
1112 #= sw[2] >= p[2] + r || =#
1113 #= ne[1] <= p[1] - r || =#
1114 #= ne[2] <= p[2] - r =#
1115 #= # particle does not intersect with cell [ix,iy] =#
1116 #= continue =#
1117
1118
1119 #= else =#
1120 #= continue =#
1121 # (likely) intersection between grid and grain
1122
1123 # 1. There is an intersection if one of the cell
1124 # corners lies within the circle. This occurs if the
1125 # distance between the cell center and corner is
1126 # less than the radii.
1127
1128 # 2. There is an intersection if one of the cell
1129 # edges comes to a closer distance to the cell
1130 # center than the radius.
1131
1132 end
1133 end
1134 end
1135 end
1136
1137 grid.porosity[ix, iy] = (cell_area - particle_area)/cell_area
1138 end
1139 end
1140 end