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