commit e07ff2fd09377cda10a4da88a74ea2e4b06f5e62
parent 55b6301591042d961a0592758304217d4afc3397
Author: Anders Damsgaard <andersd@riseup.net>
Date: Sun, 30 Apr 2017 13:04:58 -0400
change to south-west convention for corner indexing for consistency with memory layout
Diffstat:
4 files changed, 62 insertions(+), 49 deletions(-)
diff --git a/src/contact_search.jl b/src/contact_search.jl
@@ -1,17 +1,29 @@
## Contact mapping
export findContacts!
"""
+ findContacts!(simulation[, method])
+
Top-level function to perform an inter-ice floe contact search, based on ice
floe linear positions and contact radii.
The simplest contact search algorithm (`method="all to all"`) is the most
-computationally expensive (O(n^2)).
+computationally expensive (O(n^2)). The method "ocean grid" bins the ice floes
+into their corresponding cells on the ocean grid and searches for contacts only
+within the vicinity. When this method is applied, it is assumed that the
+`contact_radius` values of the ice floes are *smaller than half the cell size*.
+
+# Arguments
+* `simulation::Simulation`: the simulation object containing the ice floes.
+* `method::String`: the contact-search method to apply. Valid options are "all
+ to all" and "ocean grid".
"""
function findContacts!(simulation::Simulation,
method::String = "all to all")
if method == "all to all"
findContactsAllToAll!(simulation)
+ elseif method == "ocean grid"
+ findContactsOceanGrid!(simulation)
else
error("Unknown contact search method '$method'")
end
diff --git a/src/datatypes.jl b/src/datatypes.jl
@@ -94,11 +94,11 @@ end
#=
Type containing all relevant data from MOM6 NetCDF files. The ocean grid is a
-staggered of Arakawa-C type, with north-east convention centered on the
- h-points. During read, the velocities are interpolated to the cell corners
- (q-points).
+staggered of Arakawa-C type, with south-west convention centered on the
+h-points. During read, the velocities are interpolated to the cell corners
+(q-points).
- q(i-1, j) ------------------ q( i, j)
+ q( i,j+1) ------------------ q(i+1,j+1)
| |
| |
| |
@@ -108,10 +108,7 @@ staggered of Arakawa-C type, with north-east convention centered on the
| |
| |
| |
- q(i-1,j-1) ------------------ q( i,j-1)
-
-Source:
-https://mom6.readthedocs.io/en/latest/api/generated/pages/Horizontal_indexing.html
+ q( i, j) ------------------ q(i+1, j)
# Fields
* `input_file::String`: path to input NetCDF file
diff --git a/src/grid.jl b/src/grid.jl
@@ -1,7 +1,7 @@
"""
Use bilinear interpolation to interpolate a staggered grid to an arbitrary
-position in a cell. Assumes north-east convention, i.e. (i,j) is located at the
-north-east corner.
+position in a cell. Assumes south-west convention, i.e. (i,j) is located at the
+south-west (-x, -y)-facing corner.
# Arguments
* `field::Array{Float64, 4}`: a scalar field to interpolate from
@@ -22,10 +22,10 @@ function bilinearInterpolation(field::Array{Float64, 4},
error("relative coordinates outside bounds ($(x_tilde), $(y_tilde))")
end
- return (field[i, j, k, it]*x_tilde +
- field[i-1, j, k, it]*(1. - x_tilde))*y_tilde +
- (field[i, j-1, k, it]*x_tilde +
- field[i-1, j-1, k, it]*(1. - x_tilde))*(1. - y_tilde)
+ return (field[i+1, j+1, k, it]*x_tilde +
+ field[i, j+1, k, it]*(1. - x_tilde))*y_tilde +
+ (field[i+1, j, k, it]*x_tilde +
+ field[i, j, k, it]*(1. - x_tilde))*(1. - y_tilde)
end
export sortIceFloesInOceanGrid!
@@ -58,8 +58,8 @@ Returns the `i`, `j` index of the ocean grid cell containing the `point`.
function findCellContainingPoint(ocean::Ocean, point::Array{float, 1})
found = false
- for i=2:(size(ocean.h)[1] + 1)
- for j=2:(size(ocean.h)[2] + 1)
+ for i=1:size(ocean.h, 1)
+ for j=1:size(ocean.h, 2)
if isPointInCell(ocean, i, j, point)
return i, j
end
@@ -131,10 +131,14 @@ Returns ocean-grid corner coordinates in the following order (south-west corner,
south-east corner, north-east corner, north-west corner).
"""
function getCellCornerCoordinates(ocean::Ocean, i::Int, j::Int)
- sw = [ocean.xq[i-1, j-1], ocean.yq[i-1, j-1]]
- se = [ocean.xq[ i, j-1], ocean.yq[ i, j-1]]
- ne = [ocean.xq[ i, j], ocean.yq[ i, j]]
- nw = [ocean.xq[i-1, j], ocean.yq[i-1, j]]
+ #sw = [ocean.xq[i-1, j-1], ocean.yq[i-1, j-1]]
+ #se = [ocean.xq[ i, j-1], ocean.yq[ i, j-1]]
+ #ne = [ocean.xq[ i, j], ocean.yq[ i, j]]
+ #nw = [ocean.xq[i-1, j], ocean.yq[i-1, j]]
+ sw = [ocean.xq[ i, j], ocean.yq[ i, j]]
+ se = [ocean.xq[i+1, j], ocean.yq[i+1, j]]
+ ne = [ocean.xq[i+1, j+1], ocean.yq[i+1, j+1]]
+ nw = [ocean.xq[ i, j+1], ocean.yq[ i, j+1]]
return sw, se, ne, nw
end
diff --git a/test/grid.jl b/test/grid.jl
@@ -13,23 +13,23 @@ info("Testing area-determination methods")
@test SeaIce.areaOfQuadrilateral([1., 0.], [0., 1.], [0., 0.], [1., 1.]) ≈ 1.
info("Testing area-based cell content determination")
-@test SeaIce.isPointInCell(ocean, 2, 2, [6.5, 53.5]) == true
-@test SeaIce.getNonDimensionalCellCoordinates(ocean, 2, 2, [6.5, 53.5]) ≈
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.5, 53.5]) == true
+@test SeaIce.getNonDimensionalCellCoordinates(ocean, 1, 1, [6.5, 53.5]) ≈
[.5, .5]
-@test SeaIce.isPointInCell(ocean, 2, 2, [6.1, 53.5]) == true
-@test SeaIce.getNonDimensionalCellCoordinates(ocean, 2, 2, [6.1, 53.5]) ≈
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.5]) == true
+@test SeaIce.getNonDimensionalCellCoordinates(ocean, 1, 1, [6.1, 53.5]) ≈
[.1, .5]
-@test SeaIce.isPointInCell(ocean, 2, 2, [6.0, 53.5]) == true
-@test SeaIce.getNonDimensionalCellCoordinates(ocean, 2, 2, [6.0, 53.5]) ≈
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.0, 53.5]) == true
+@test SeaIce.getNonDimensionalCellCoordinates(ocean, 1, 1, [6.0, 53.5]) ≈
[.0, .5]
-@test SeaIce.isPointInCell(ocean, 2, 2, [6.1, 53.7]) == true
-@test SeaIce.isPointInCell(ocean, 2, 2, [6.1, 53.9]) == true
-@test SeaIce.isPointInCell(ocean, 2, 2, [6.1, 53.99999]) == true
-@test SeaIce.getNonDimensionalCellCoordinates(ocean, 2, 2, [6.1, 53.99999]) ≈
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.7]) == true
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.9]) == true
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.99999]) == true
+@test SeaIce.getNonDimensionalCellCoordinates(ocean, 1, 1, [6.1, 53.99999]) ≈
[.1, .99999]
-@test SeaIce.isPointInCell(ocean, 2, 2, [7.5, 53.5]) == false
-@test SeaIce.isPointInCell(ocean, 2, 2, [0.0, 53.5]) == false
-x_tilde, _ = SeaIce.getNonDimensionalCellCoordinates(ocean, 2, 2, [0., 53.5])
+@test SeaIce.isPointInCell(ocean, 1, 1, [7.5, 53.5]) == false
+@test SeaIce.isPointInCell(ocean, 1, 1, [0.0, 53.5]) == false
+x_tilde, _ = SeaIce.getNonDimensionalCellCoordinates(ocean, 1, 1, [0., 53.5])
@test x_tilde < 0.
info("Testing conformal mapping methods")
@@ -55,16 +55,16 @@ info("Testing conformal mapping methods")
[7.5,-1.5])
info("Checking cell content using conformal mapping methods")
-@test SeaIce.isPointInCell(ocean, 2, 2, [6.4, 53.4], method="Conformal") == true
-@test SeaIce.isPointInCell(ocean, 2, 2, [6.1, 53.5], method="Conformal") == true
-@test SeaIce.isPointInCell(ocean, 2, 2, [6.0, 53.5], method="Conformal") == true
-@test SeaIce.isPointInCell(ocean, 2, 2, [6.1, 53.7], method="Conformal") == true
-@test SeaIce.isPointInCell(ocean, 2, 2, [6.1, 53.9], method="Conformal") == true
-@test SeaIce.isPointInCell(ocean, 2, 2, [6.1, 53.99999],
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.4, 53.4], method="Conformal") == true
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.5], method="Conformal") == true
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.0, 53.5], method="Conformal") == true
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.7], method="Conformal") == true
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.9], method="Conformal") == true
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.99999],
method="Conformal") == true
-@test SeaIce.isPointInCell(ocean, 2, 2, [7.5, 53.5],
+@test SeaIce.isPointInCell(ocean, 1, 1, [7.5, 53.5],
method="Conformal") == false
-@test SeaIce.isPointInCell(ocean, 2, 2, [0.0, 53.5],
+@test SeaIce.isPointInCell(ocean, 1, 1, [0.0, 53.5],
method="Conformal") == false
info("Testing bilinear interpolation scheme on conformal mapping")
@@ -72,13 +72,13 @@ ocean.u[1, 1, 1, 1] = 1.0
ocean.u[2, 1, 1, 1] = 1.0
ocean.u[2, 2, 1, 1] = 0.0
ocean.u[1, 2, 1, 1] = 0.0
-@test SeaIce.bilinearInterpolation(ocean.u, .5, .5, 2, 2, 1, 1) ≈ .5
-@test SeaIce.bilinearInterpolation(ocean.u, 1., 1., 2, 2, 1, 1) ≈ .0
-@test SeaIce.bilinearInterpolation(ocean.u, 0., 0., 2, 2, 1, 1) ≈ 1.
-@test SeaIce.bilinearInterpolation(ocean.u, .25, .25, 2, 2, 1, 1) ≈ .75
-@test SeaIce.bilinearInterpolation(ocean.u, .75, .75, 2, 2, 1, 1) ≈ .25
+@test SeaIce.bilinearInterpolation(ocean.u, .5, .5, 1, 1, 1, 1) ≈ .5
+@test SeaIce.bilinearInterpolation(ocean.u, 1., 1., 1, 1, 1, 1) ≈ .0
+@test SeaIce.bilinearInterpolation(ocean.u, 0., 0., 1, 1, 1, 1) ≈ 1.
+@test SeaIce.bilinearInterpolation(ocean.u, .25, .25, 1, 1, 1, 1) ≈ .75
+@test SeaIce.bilinearInterpolation(ocean.u, .75, .75, 1, 1, 1, 1) ≈ .25
info("Testing cell binning")
-@test SeaIce.findCellContainingPoint(ocean, [6.2, 53.4]) == (2, 2)
-@test SeaIce.findCellContainingPoint(ocean, [7.2, 53.4]) == (3, 2)
+@test SeaIce.findCellContainingPoint(ocean, [6.2, 53.4]) == (1, 1)
+@test SeaIce.findCellContainingPoint(ocean, [7.2, 53.4]) == (2, 1)
@test_throws ErrorException SeaIce.findCellContainingPoint(ocean, [0.2, 53.4])
